Discrete Dynamics in Nature and Society

Volume 2013 (2013), Article ID 681689, 9 pages

http://dx.doi.org/10.1155/2013/681689

Research Article

Disease Spread in Coupled Populations: Minimizing Response Strategies Costs in Discrete Time Models

Geisel Alpízar¹ and Luis F. Gordillo²

¹Departamento de Matemáticas, Instituto Tecnológico de Costa Rica, Cartago 30101, Costa Rica

²Department of Mathematics and Statistics, Utah State University, Logan, UT 84322, USA

Received 6 September 2012; Accepted 19 April 2013

Academic Editor: Francisco Solís Lozano

Copyright © 2013 Geisel Alpízar and Luis F. Gordillo. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Social distancing, vaccination, and medical treatments have been extensively studied and widely used to control the spread of infectious diseases. However, it is still a difficult task for health administrators to determine the optimal combination of these strategies when confronting disease outbreaks with limited resources, especially in the case of interconnected populations, where the flow of individuals is usually restricted with the hope of avoiding further contamination. We consider two coupled populations and examine them independently under two variants of well-known discrete time disease models. In both examples we compute approximations for the control levels necessary to minimize costs and quickly contain outbreaks. The main technique used is simulated annealing, a stochastic search optimization tool that, in contrast with traditional analytical methods, allows easy implementation to any number of patches with different kinds of couplings and internal dynamics.

1. Introduction

In contrast with the tremendous benefits that modern transportation systems have brought to our society, their potential as contributors to the global and fast spread of infectious agents has become a major concern [1, 2]. This has been evidenced, for instance, with the recent SARS and swine flu outbreaks in 2003 and 2009, respectively. In those occasions, careful monitoring and efficient control of individuals’ flow among cities were implemented, although it has been determined that introducing travel restrictions does not have a drastic impact on a pandemic development, [3, 4]. However, this kind of measures may introduce a delay in the spread of the disease [5], providing additional time to implement nonpharmaceutical interventions. The important role that transportation has on the spread of infectious diseases has been extensively studied in particular cases. Remarkable efforts are currently made to understand the spread of pandemic and seasonal influenza in spatial structured populations, [6–8]. The detailed study made in [9] shows how social distancing policies implanted in Mexico in 2009, combined with control in the transportation of individuals, lead to the effective interruption of the first two waves of the influenza’s emergence.

Health administrators face the nontrivial task of controlling disease spread through optimal responses during the first stages of an outbreak, which include the implementation of social distancing, travel restrictions, medical treatments, and vaccination. Coupled populations scenarios may be difficult to assess if confronted with scarce resources. One of our goals is to show that simulated annealing, a well-known computational optimization technique, may be useful in the search of appropriate responses for some diseases modeled in discrete time. We use variants of two extensively studied models: SIR and SIS in coupled populations, which are introduced in Section 2. The computational results obtained by application of simulated annealing are presented in Section 3. This technique has the advantage of combining easy implementation and simple theoretical principles [10].

Although theoretical models that describe disease dynamics in metapopulations have been studied in detail, the problem of optimal deployment of control combining epidemiological and economic factors has barely been touched. For continuous time, the SIS formulation with two interconnected patches is successfully addressed in [11], and the optimal allocation of resources for SIRS models in metapopulations has been recently studied in [12]. For discrete time epidemic models, optimal control techniques have been employed only for one patch populations [13, 14]. The usual approach in this cases is to use Pontryagin’s maximum principle adapted to discrete systems [14, 15], which establishes necessary conditions for the existence of an optimal condition that can be recovered using a forward-backward algorithm. This theoretical formulation might turn out to be cumbersome if a large number of interconnected patches combined with a high number of model parameters in each patch are involved. In contrast, the alternative of approximating optimal solutions with numerical simulations only requires a very general framework, easily adaptable to any situation.

2. Coupled Populations: Two Examples

For the sake of simplicity, we consider the case of only two coupled populations, for example and . The disease dynamics in each patch are described with discrete time equations, capturing the processes of infection and migration at each unit in time (day). We assume that the system does not allow the introduction of individuals from outside these patches and that, within each population, individuals are homogeneously mixed.

The technique used below to find optimal control parameter values is very general and independent of the type of discrete time model describing the disease dynamics within each patch. Thus the method could be, in principle, adjusted to numberless possible modeling situations. In this paper, we consider separately for analysis two instances of simple but plausible schemes largely employed to describe dynamics of infectious diseases. For both models we assume that the epidemic evolves fast enough to ignore demographic effects, and consequently our attention is focused on the disease spread between populations only for the first few days after an outbreak appears in one of the patches. The first model, referred to here as Model A, is a mechanistic extension of a discrete approximation to the basic Susceptible-Infected-Removed (SIR), which is valid for short periods of time and examined here due to its simplicity. The motivation for this model comes from the ideas originally explored in [16]. The second model, Model B, is a Susceptible-Infected-Susceptible (SIS) type, derived originally in [17], which might be useful indescribing infectious diseases where almost immediate host reinfection is possible.

In addition to each model we define functionals that represent the economic impact caused by the disease, which include the costs of direct intervention and those generated by reducing susceptibility of individuals. Direct intervention practices include the reduction of contact rates through social distancing, decreasing individual infectiousness with appropriate treatment, or isolation. Vaccination and preventive care, which in contrast reduce the susceptibility of a population [18, 19], are not going to be taken into account here.

2.1. Model A: Dispersal in Two Patches with SIR

In this case, the disease dynamics in each patch are described with a basic discrete numerical approximation to the SIR model, valid for a short period of time. Susceptible and infected individuals are allowed to leave their original (home) group at a certain rate and return to it, taking an active part in the infectious process while visiting the second (temporary) group. Under these assumptions [16] presents a description for the mixing patterns among patches in continuous time, where the original mechanistic model (which includes demography) generates results that are successfully linked to the classical phenomenological formulations.

Following the ideas in [16], we formulate the discrete time model: let be the class to which an individual may belong (: susceptible or : infected), and let represent the number of individuals that belong to patch being currently located in patch at time . Sometimes, the flow of individuals from one patch to another may not be symmetrical and it would be desirable to keep track of the amount of individuals from one population that are in a different location. Thus, let be the rate at which individuals of class   return from group to group and the rate at which individuals leave group towards group . Similar definitions are made for and . Let be the contact rate for both groups and the natural removal rate of infectives. It is important to notice that variations in patch population may cause changes in the effective average contact rates. However, we assume that each patch has a large number of individuals and that travelers excursions to other patches are made only by a small amount of individuals in each population. Thus, we neglect fluctuations on .

The following equations capture the process of infection and subsequent movement of individuals between patches, in that order, at each time step for patch :

Similar equations for patch are obtained interchanging and . Although at each step in time we assume that the infection process runs first and then the dispersal, the order does not matter.

Figure 1 shows a reemergence of the total number of infected individuals from both patches, caused by a decreased flow of infected individuals from patch to patch , which is assumed initially to be disease-free. Relaxing the control on the infected individuals flow towards group causes an apparent merge of the two peaks, as it would be expected.

(a)

(b)

(c)

Figure 1: The impact of the variation in (here = ) on the total number of infected individuals along 30 days for different initial values of infected individuals in patch : (a) , (b) , and (c) . Varying individuals’ flow intensity from one patch to another causes a notably increased appearance of infectives in the SIR model, due to the delayed entrance of infecteds. For this example we assume that patch is initially to be disease-free and that we are interested in controlling the flow of infected individuals towards it. Accordingly, we take . The other initial values are , , and . The parameter values are chosen from a particular influenza episode modeled with an SIR scheme [20], , . The values for the fractions of individuals moving between patches are chosen in a range estimated for patches with high interconnectivity [9]: , , and .

2.1.1. Cost Functional for Model A

For each time step, the flow of infected individuals is described by the associated parameters and . Thus, we define the control vectors and to keep track in time of the movement allowed between patches, corresponding to the fraction of individuals that originally belong to group and are returning and the fraction of individuals that belong to group and are leaving towards group , respectively. We define and similarly. For this model, it is assumed that we have only control on the flow of individuals, and no other intervention strategy is adopted.

Let be the prevalence in the group at time . Also, define the cost functional associated with the epidemic impact and control of dispersal for the patch y by

For simplicity, we set and consider the costs of the spread, and , relative to those generated directly by each infected individual. A similar expression for the cost functional can be written for patch , , and the functional to minimize is .

2.2. Model B: Dispersal in Two Patches with SIS

The second model considered is a slight variant of the SIS model for two coupled patches originally presented and analyzed in [17]. The underlying motivation to consider this example is given by infectious diseases that are likely to behave as in a pandemic influenza scenario: an individual may be reinfected with the same virus strain if an infectious contact occurs before her/his primary antibody response matures, a stage that usually should be reached after three to four weeks since the first infection. Reinfections are more likely to occur during pandemic influenza because the virus circulates more extensively than in nonpandemic events [21].

For patch , let and represent the population of susceptible and infected individuals at time . Recall that we do not take into accoutn the effects of demographic changes and denote by the total population in patch at time . Let represent the fraction of individuals that recover naturally from the disease at each time step and a constant that weights the role of prevalence on disease transmission at time for both groups. First, we write the equations for the dynamics of the disease within a unit of time:

Similar equations are obtained for patch replacing by . The dynamics of more general equations than these have been studied for a single population in [22]. Now, we introduce the diffusion of individuals between patches: assume that fractions and of susceptible and infected individuals from each population are exchanged at each step time. We assign superscripts to distinguish the diffusing fractions associated to each patch:

Notice that although the disease may disappear locally in one location, global persistence may be maintained throughout the dispersion of infected individuals [17, 23].

2.2.1. Cost Functional for Model B

The impact of direct intervention measures, like social distancing, produces changes in the value of . More precisely, let be the transmission parameter within group after intervention strategies have been applied to reduce transmission, which are measured by the parameter . Let be the fraction of infected individuals that recover naturally from the disease. Suppose that a fraction of the infected individuals that did not recover by natural means is effectively treated and returns to the susceptible class. The motivation behind this is the treatments that only inhibit virus replication but immune response is still needed to control the infection.

After the introduction of these controls, the equations for patch read as

Similar equations result for patch . As in Model A, our interest is to include the possible restriction of infected individuals movements between patches. The impact of these measures is directly captured by the parameters and .

Let and be the vectors that have in each entry the level of direct intervention at each time step. Then, a cost functional can be written as

where hats indicate the prevalence of the disease in the patch considered at time . The proportion of infected individuals who move from one patch to another is low when the restrictions of traveling are high, and consequently a high cost should be expected. For this reason the term in the functional corresponding to the application of restrictions on the dispersion of infected individuals is defined by : we require that (1) is at most equal to when no restrictions are applied and (2) the costs should increase when decreases from to 0.

In the objective functional, the constants represent the relative costs of the control measures: is the associated cost produced by a new infection, is the cost of the implementation of social distancing, is the relative cost associated with the application of treatment, and is the cost of applying restrictions on infected individuals movement. We set , and let , , be the relative costs in respect to [14].

3. Numerical Results

In this section we use simulated annealing to assess the challenge of finding appropriate control efforts levels that minimize the cost functionals for each model previously described. Simulated annealing is a stochastic search technique derived from the well-known Metropolis algorithm see the (appendix for a brief review or [10, 24, 25] for more detailed and complete accounts). This computational tool becomes very helpful for finding approximations to function extrema when the number of parameters in the system is large and consequently the model turns out to be analytically intractable.

Simulations for Model A are made only considering movement restrictions on individuals, and the results are compared with the worst case scenario (no movement restriction). In practice, it is very likely that the introduction of direct intervention measures during the first days of an epidemic development would be delayed because of practical or economic constrains. This observation motivates excluding this type of measures for this case.

For Model B, simulations combine movement restrictions with direct intervention, motivated by the practices and results obtained for pandemic influenza, and compare the results to the worst case scenario (no control of any sort applied).

3.1. Numerical Results for Model A

The simulations are obtained considering the particular case of one patch starting at the disease-free state, say population , and then applying control on the dispersal of infected individuals towards this patch from , which initially had a positive number of infectives. We assign an arbitrary small positive lower bound to the spread of infectives, , , considering that complete restriction of infected individuals flow is frequently an unavoidable fact in real scenarios. Other parameter values are taken the same as in Figure 1. The simulation results are summarized in Figure 2: the worst case scenario (no movement restriction) is compared to the disease evolution obtained by using the optimal control efforts computed by the algorithm. It is observed that restriction in the movement from patch to is imposed only several days after the start of the epidemic in patch . In contrast, the restriction of movement from to is complete since the beginning; see Figure 2(c). The values were used for the simulations.

(a)

(b)

(c)

Figure 2: (a) The evolution of the worst case scenario, that is, when no control is applied, is shown in bullets () for susceptibles and infected classes in each patch. Using the parameter values found by the algorithm provides the scenario shown in (▲). (b) The total amount of infected individuals in both populations. Again, the worst case scenario is in () and the minimized results in (▲). The appearance of a second bump is due to the delay of introducing infective individuals in the initially disease-free population. (c) These are the strategy levels that minimize the total cost.

3.2. Numerical Results for Model B

For the simulations in this example we use the parameter and baseline values given in [14] for pandemic influenza. The model studied in [14] describes the disease dynamics for a long period time, without taking into account reinfection during the first weeks, and includes additional compartments in the population (asymptomatic, treated, recovered, and dead). We emphasize that our interest is on the first weeks after the pandemic started, when individuals’ reinfection is still plausible. We consider the values for the fraction of infected individuals recovered by natural means. In absence of controls, the basic reproductive number range considered in [14] for high transmission is , which agrees with known estimates [26]. The corresponding transmission rate for the worst case is used here. Assuming large populations allows for regarding possible variations in the transmission rate as negligible. The control bounds are interpreted as the maximum and minimum daily rates, also chosen here as presented in [14] for the social distancing and treatment: and . Remember that is the fraction of the infected individuals that recover by the application of treatment. The bounds for the spread of infectives are the same used in Model A, and , where . For the simulations we choose the values , , and , suggested by [14], and assume that the costs associated with restrictions on the dispersion of infected individuals are higher than those associated with social distancing and treatment.

We fix the initial conditions and and also and . With these values we observe in Figure 3(a) how optimal levels of control efforts drastically reduce infected individuals in the population compared to the worst case scenario. The impact of the optimal control efforts on the infected population size is presented in Figure 3(b): implementing control measures delays the appearance of a reduced maximum number of total infecteds. Figure 3(c) shows the optimal control values. In population , where initially , it is required to apply maximum treatment, social distancing and restriction in movement towards the population . Population should apply maximum social distancing, but treatment should gradually grow and reach its maximum around the tenth day. No restrictions apply for infecteds’ flow starting at .

(a)

(b)

(c)

Figure 3: (a) The population of susceptible and infected individuals for the worst case () and the results of applying the minimizing algorithm (▲). (b) Total of infected individuals in both populations presents the worst case () and the minimized results (▲). The delay in the introduction of infective individuals produces an apparent reemergence. (c) Strategy that minimizes the combined costs of the disease impact and the control of the disease spread.

4. Discussion

The increased extensivity and intensity of global interconnections jointly with the extreme dynamic character of current transportation systems turn out to be important to understand the role of mobility in the spread of infectious diseases, as illustrated with the SARS outbreak in 2003 [1]. Currently, many countries have developed or adapted response plans for pandemics of infectious diseases, which generally include the monitoring and control of travelers flow from abroad and between its own cities. Generally, in the early stages of a pandemic development, local health departments implement domestic travel-related measures to slow disease spread [27], like determining travelers’ condition through fast laboratory tests followed by movement restriction. Fast disease spread may require quick computation of good approximations to the resulting dynamics from the application of combined strategies under conditions of scarce resources. These approximations may turn into a valuable tool to assess the risk of disease spread and gain understanding on the consequences of policies impact within the contemporary contexts of globalization.

In this note we show how stochastic search techniques may provide valuable insight into the quest of finding optimal strategies for metapopulations epidemic control with discrete time models, when a quick response is needed and analytical tools are not viable. We provide approximate solutions computed for two independent recurrent systems where restriction on individuals’ flow between populations is regulated and combined in one of the examples with direct intervention measures.

The huge complexity involved in the spread of diseases among populations is reflected in the analytically untractable number of equations that would result in the modeling process. Optimizing functionals over these systems may therefore turn into a very complicated task. As initially remarked, optimal control techniques have been used only to explore one patch populations [13, 14], evidencing analytical difficulties in more complicated models. Our approach to overcome this problem consists in suggesting the use of stochastic search techniques. Stochastic search allows not only for considering different population sizes for each patch as well as different ways of coupling but aslo any number of populations (if provided with enough computational power), each with its own dynamics. This advantage is not found in the current analytical methods.

The numerical results for the SIR extension presented here suggest that the optimal travel restrictions create a sudden increase in the total number of infected individuals: the dynamics of the disease in the second patch are delayed for several days. The delay in the transmission of the disease to the second patch can be used to prepare and implement more efficient nonpharmaceutical strategies to reduce the disease impact, like developing public awareness, instituting social distancing, or organizing vaccination centers [5]. For the SIS example, the total number of infected individuals stabilizes at a lower level than the worst case scenario, after the appearance of a little bump generated by the delayed introduction of infected individuals into the second population.

Appendix

The Simulated Annealing Algorithm

For the sake of completeness we include a brief description of how simulated annealing works. For more detailed accounts the reader is encouraged to have look at the cited references. Essentially, the algorithm runs a search in a space of states for an element that globally minimize (or maximize) the value of a function. In our case, the space of states is given by a multidimensional lattice where each point is a vector with the feasible controls at each time as entries. For instance, in Model A we look for the -dimensional vector:

that minimizes the cost functional defined by ; see (2). The grid is evidently chosen in such a way that has practical significance and computational feasibility.

Once the cost functional and the finite set of states have been defined, the following algorithm can be used to find the minimum. First choose an initial state and define an inhomogeneous Markov chain that evolves as follows: if is state , that is, , choose randomly a neighbor state that can be reached in one transition step. If , then set . Otherwise, if , then set with probability

and with probability . It can be shown [10] that the probability of choosing an element from state space that minimizes approaches 1 when the temperature parameter tends to 0. The choice of an appropriate sequence , called cooling schedule, is important: if the cooling is too fast, the chain may get stuck in a local minimum. So there is a balance that has to be reached, where you want a reasonable time of computation to observe convergence but slow enough to avoid convergence to an element that is not a global minimizer. For the analysis of cooling functions see [25].

References

S. H. Ali and R. Keil, “Contagious cities,” Geography Compass, vol. 1, no. 5, pp. 1207–1226, 2007.

V. Colizza, A. Barrat, M. Barthélemy, and A. Vespignani, “The role of the airline transportation network in the prediction and predictability of global epidemics,” Proceedings of the National Academy of Sciences of the United States of America, vol. 103, no. 7, pp. 2015–2020, 2006.

B. S. Cooper, R. J. Pitman, W. J. Edmunds, and N. J. Gay, “Delaying the international spread of pandemic influenza,” PLoS Medicine, vol. 3, no. 6, articl e212, 2006.

T. D. Hollingsworth, N. M. Ferguson, and R. M. Anderson, “Will travel restrictions control the international spread of pandemic influenza?” Nature Medicine, vol. 12, no. 5, pp. 497–499, 2006.

J. M. Epstein, D. M. Goedecke, F. Yu, R. J. Morris, D. K. Wagener, and G. V. Bobashev, “Controlling pandemic flu: the value of international air travel restrictions,” PLoS ONE, vol. 2, no. 5, article e401, 2007.

T. D. Hollingsworth, N. M. Ferguson, and R. M. Anderson, “Frequent travelers and rate of spread of epidemics,” Emerging Infectious Diseases, vol. 13, no. 9, pp. 1288–1294, 2007.

J. M. Hyman and T. LaForce, “Modeling the spread of inuenza among cities,” in Bioterrorism, Mathematical Modeling Applications in Homeland Security, H. T. Banks and C. Castillo-Chavez, Eds., pp. 215–240, SIAM Frontiers in Applied Mathematics, 2003.

L. A. Rvachev and I. M. Longini Jr., “A mathematical model for the global spread of influenza,” Mathematical Biosciences, vol. 75, pp. 3–22, 1985.

M. A. Herrera-Valdez, M. Cruz-Aponte, and C. Castillo-Chávez, “Multiple outbreaks for the same pandemic: local transportation and social distancing explain the different “waves” of A-H1N1pdm cases observed in México during 2009,” Mathematical Biosciences and Engineering, vol. 8, no. 1, pp. 21–48, 2011.

O. Häggström, Finite Markov Chains and Algorithmic Applications, vol. 52, Cambridge University Press, Cambridge, UK, 2002.

R. E. Rowthorn, R. Laxminarayan, and C. A. Gilligan, “Optimal control of epidemics in metapopulations,” Journal of the Royal Society Interface, vol. 6, no. 41, pp. 1135–1144, 2009.

M. L. Ndeffo Mbah and C. A. Gilligan, “Resource allocation for epidemic control in metapopulations,” PLoS ONE, vol. 6, no. 9, Article ID e24577, 2011.

W. Ding and S. Lenhart, “Introduction to optimal control for discrete time models with an application to disease modeling,” DIMACS Series in Discrete Mathematics and Theoretical Computer Science, vol. 75, pp. 109–119, 2010.

P. A. González-Parra, S. Lee, L. Velázquez, and C. Castillo-Chávez, “A note on the use of optimal control on a discrete time model of influenza dynamics,” Mathematical Biosciences and Engineering, vol. 8, no. 1, pp. 183–197, 2011.

S. Lenhart and J. T. Workman, Optimal Control Applied to Biological Models, CRC Mathematical and Computational Biology Series, Chapman & Hall, Boca Raton, Fla, USA, 2007.

M. J. Keeling and P. Rohani, “Estimating spatial coupling in epidemiological systems: a mechanistic approach,” Ecology Letters, vol. 5, no. 1, pp. 20–29, 2002.

C. Castillo-Chávez and A.-A. Yakubu, “Intraspecific competition, dispersal and disease dynamics in discrete time patchy environments,” in Mathematical Approaches for Emerging and reemerging Infectious Diseases: Models, Methods and Theory, C. Castillo-Chávez, P. van den Driessche, D. Kirschner, and A.-A. Yakubu, Eds., pp. 165–181, Springer, New York, NY, USA, 1999.

N. M. Ferguson, D. A. T. Cummings, S. Cauchemez et al., “Strategies for containing an emerging influenza pandemic in Southeast Asia,” Nature, vol. 437, no. 7056, pp. 209–214, 2005.

M. Lipsitch, T. Cohen, B. Cooper et al., “Transmission dynamics and control of severe acute respiratory syndrome,” Science, vol. 300, no. 5627, pp. 1966–1970, 2003.

J. D. Murray, Mathematical Biology. I. An introduction, vol. 17 of Interdisciplinary Applied Mathematics, Springer, New York, NY, USA, 3rd edition, 2002.

C. M. Perez, M. Ferres, and J. A. Labarca, “Pandemic (H1N1) 2009 reinfection, Chile,” Emerging Infectious Diseases, vol. 16, no. 1, pp. 156–157, 2010.

C. Castillo-Chávez and A.-A. Yakubu, “Discrete-time SIS models with simple and complex population dynamics,” in Mathematical Approaches for Emerging and Reemerging Infectious Diseases: Models, Methods and Theory, C. Castillo-Chávez, P. van den Driessche, D. Kirschner, and A.-A. Yakubu, Eds., pp. 153–163, Springer, New York, NY, USA, 1999.

M. J. Keeling and C. A. Gilligan, “Metapopulation dynamics of bubonic plague,” Nature, vol. 407, no. 6806, pp. 903–906, 2000.

S. Kirkpatrick, C. D. Gelatt, Jr., and M. P. Vecchi, “Optimization by simulated annealing,” Science, vol. 220, no. 4598, pp. 671–680, 1983.

J. C. Spall, Introduction to Stochastic Search and Optimization: Estimation, Simulation, and Control, Wiley-Interscience Series in Discrete Mathematics, Wiley-Interscience, Hoboken, NJ, USA, 2003.

G. Chowell, H. Nishiura, and L. M. A. Bettencourt, “Comparative estimation of the reproduction number for pandemic influenza from daily case notification data,” Journal of the Royal Society Interface, vol. 4, no. 12, pp. 154–166, 2007.

United States Department of Health and Human Services, “HHS Pandemic Inuenza Plan,” http://www.flu.gov/planning-preparedness/federal/hhspandemicinfluenzaplan.pdf