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

,


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][7][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 wellknown 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.

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.

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.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   = 1, 500,   =   = 0, and   = 1, 000.The parameter values are chosen from a particular influenza episode modeled with an SIR scheme [20],  = 10 −3 /day,  = 0.44/day.The values for the fractions of individuals moving between patches are chosen in a range estimated for patches with high interconnectivity [9]:    =    = 0.03,    =    = 0.02, and    =    = 0.02.

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   = (   (0), . . .,    ( − 1)) and   = (   (0), . . .,    ( − 1)) 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  1 = 1 and consider the costs of the spread,  2 and  3 , 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   +   .

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    =    +    > 0 the total population in patch  at time .Let 0 ≤ 1 −   ≤ 1 represent the fraction of individuals that recover naturally from the disease at each time step and 0 ≤  ≤ 1 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].) 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 (    −    ) 2 : 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:  1 is the associated cost produced by a new infection,  2 is the cost of the implementation of social distancing,  3 is the relative cost associated with the application of treatment, and  4 is the cost of applying restrictions on infected individuals movement.We set  1 = 1, and let  2 ,  3 ,  4 be the relative costs in respect to  1 [14].

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).

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,    ,    ≥ 0.001, 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  2 =  3 = 0.2 were used for the simulations.

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 (1−  ) = (1−  ) = 1/7 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 [2.4,3.2], which agrees with known estimates [26].The corresponding transmission rate for the worst case  = 1.94 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:   ,   ∈ [0, 0.2] and   ,   ∈ [0.95, 1].Remember that 1 −   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,    ∈ [0.001,  For the simulations we choose the values  2 = 0.04,  4 = 0.0004, and  3 = 0.004, 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   1 = 100 and   1 = 0 and also   1 = 1, 500 and   1 = 1, 500.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   1 = 100, 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 .

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.

Figure 1 :
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)   0 = 1, (b)   0 = 10, and (c)   0 = 100.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 ] and    ∈ [0.001,    ], where    =    = 0.03.

Figure 2 :
Figure2: (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.

Figure 3 :
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.
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   = (1 −   ) be the transmission parameter within group  after intervention strategies have been applied to reduce transmission, which are measured by the parameter .Let 1 −  be the fraction of infected individuals that recover naturally from the disease.Suppose that a fraction 1 −  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.