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

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 𝑥: 𝑆^𝑥𝑥_𝑘+1=1−𝜌^𝑆_𝑥𝑦𝑆_𝑘^𝑥𝑥−𝛽𝑆_𝑘^𝑥𝑥𝐼_𝑘^𝑥𝑥+𝐼_𝑘^𝑦𝑥+𝜏^𝑆_𝑦𝑥𝑆_𝑘^𝑥𝑦−𝛽𝑆_𝑘^𝑥𝑦𝐼_𝑘^𝑥𝑦+𝐼_𝑘^𝑦𝑦,𝑆^𝑥𝑦_𝑘+1=1−𝜏^𝑆_𝑦𝑥𝑆_𝑘^𝑥𝑦−𝛽𝑆_𝑘^𝑥𝑦𝐼_𝑘^𝑥𝑦+𝐼_𝑘^𝑦𝑦+𝜌^𝑆_𝑥𝑦𝑆_𝑘^𝑥𝑥−𝛽𝑆_𝑘^𝑥𝑥𝐼_𝑘^𝑥𝑥+𝐼_𝑘^𝑦𝑥,𝐼^𝑥𝑥_𝑘+1=1−𝜌^𝐼_𝑥𝑦𝐼_𝑘^𝑥𝑥+𝛽𝑆_𝑘^𝑥𝑥𝐼_𝑘^𝑥𝑥+𝐼_𝑘^𝑦𝑥−𝑑𝐼_𝑘^𝑥𝑥+𝜏^𝐼_𝑦𝑥𝐼_𝑘^𝑥𝑦+𝛽𝑆_𝑘^𝑥𝑦𝐼_𝑘^𝑦𝑦+𝐼_𝑘^𝑥𝑦−𝑑𝐼_𝑘^𝑥𝑦,𝐼^𝑥𝑦_𝑘+1=1−𝜏^𝐼_𝑦𝑥𝐼_𝑘^𝑥𝑦+𝛽𝑆_𝑘^𝑥𝑦𝐼_𝑘^𝑥𝑦+𝐼_𝑘^𝑦𝑦−𝑑𝐼_𝑘^𝑥𝑦+𝜌^𝐼_𝑥𝑦𝐼_𝑘^𝑥𝑥+𝛽𝑆_𝑘^𝑥𝑥𝐼_𝑘^𝑥𝑥+𝐼_𝑘^𝑦𝑥−𝑑𝐼_𝑘^𝑥𝑥. (1) 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) 𝐼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 𝐼0𝑦𝑦=𝐼0𝑥𝑦=𝐼0𝑦𝑥=0. 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.

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 𝜏_𝑥𝑦=(𝜏^𝐼_𝑥𝑦(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 𝐽_𝑦𝐼^𝑦̂,𝜏_𝑥𝑦,𝜌_𝑥𝑦=^𝑇−1_𝑘=0𝐵₁𝐼̂^𝑦_𝑘²+𝐵₂𝜏^𝐼_𝑥𝑦(𝑘)−𝜏^𝑆_𝑥𝑦²+𝐵₃𝜌^𝐼_𝑥𝑦(𝑘)−𝜌^𝑆_𝑥𝑦²+𝐼̂^𝑦_𝑇². (2) For simplicity, we set 𝐵₁=1 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 𝑇^𝑥_𝑘=𝑆^𝑥_𝑘+𝐼^𝑥_𝑘>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: 𝑆^𝑥_𝑘=𝑆^𝑥_𝑘exp−𝛽𝐼^𝑥_𝑘𝑇^𝑥_𝑘+1−𝜎_𝑥𝐼^𝑥_𝑘,𝐼^𝑥_𝑘̃=1−exp−𝛽𝐼^𝑥_𝑘𝑇^𝑥_𝑘𝑆^𝑥_𝑘+𝜎_𝑥𝐼^𝑥_𝑘. (3) 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: 𝑆^𝑥_𝑘+1=1−𝐷^𝑥_𝑆𝑆^𝑥_𝑘+𝐷^𝑦_𝑆𝑆^𝑦_𝑘,𝐼^𝑥_𝑘+1=1−𝐷^𝑥_𝐼𝐼^𝑥_𝑘̃+𝐷^𝑦_𝐼𝐼^𝑦_𝑘̃,𝑆^𝑦_𝑘+1=𝐷^𝑥_𝑆𝑆^𝑥_𝑘+1−𝐷^𝑦_𝑆𝑆^𝑦_𝑘,𝐼^𝑦_𝑘+1=𝐷^𝑥_𝐼𝐼^𝑥_𝑘̃+1−𝐷^𝑦_𝐼𝐼^𝑦_𝑘̃. (4) 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 𝛽_𝑥=(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.

After the introduction of these controls, the equations for patch 𝑥 read as 𝑆^𝑥_𝑘=𝑆^𝑥_𝑘exp−𝛽_𝑥𝐼^𝑥_𝑘𝑇^𝑥_𝑘+1−𝜏𝜎_𝑥𝐼^𝑥_𝑘,𝐼^𝑥_𝑘̃=1−exp−𝛽_𝑥𝐼^𝑥_𝑘𝑇^𝑥_𝑘𝑆^𝑥_𝑘+𝜏𝜎_𝑥𝐼^𝑥_𝑘. (5) 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 𝐟_𝑥=(𝑓^𝑥₁,…,𝑓^𝑥_𝑇−1) and 𝐟_𝑦=(𝑓^𝑦₁,…,𝑓^𝑦_𝑇−1) 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𝐽𝐼,̂𝐟_𝑥,𝐟_𝑦,𝜏=_{𝑤=𝑥,𝑦}^𝑇−1_𝑘=0𝐵₁𝐼̂^𝑤_𝑘²+𝐵₂𝑓^𝑤_𝑘²+𝐵₃1−𝜏_𝑘²+𝐵₄𝐷^𝑤_𝐼(𝑘)−𝐷^𝑤_𝑆²+𝐼̂^𝑤_𝑇², (6) 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 𝐵₁=1, 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, 𝜌^𝐼_𝑥𝑦, 𝜏^𝐼_𝑥𝑦≥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 𝐵₂=𝐵₃=0.2 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 (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,𝐷^𝑦_𝑆] and 𝐷^𝑥_𝐼∈[0.001,𝐷^𝑦_𝑆], where 𝐷^𝑥_𝑆=𝐷^𝑦_𝑆=0.03. For the simulations we choose the values 𝐵₂=0.04, 𝐵₄=0.0004, and 𝐵₃=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 𝐼^𝑥₁=100 and 𝐼^𝑦₁=0 and also 𝑆^𝑥₁=1,500 and 𝑆^𝑦₁=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 𝐼^𝑥₁=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 𝑦.

(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 4𝑇-dimensional vector: 𝜏^𝐼_𝑥𝑦(0),…,𝜏^𝐼_𝑥𝑦(𝑇−1),𝜌^𝐼_𝑥𝑦(0),…,𝜌^𝐼_𝑥𝑦(),𝜏𝑇−1_𝐼𝑦𝑥(0),…,𝜏^𝐼_𝑦𝑥(𝑇−1),𝜌^𝐼_𝑦𝑥(0),…,𝜌^𝐼_𝑦𝑥(𝑇−1) (A.1) 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 𝑋(𝑘+1)=𝑗. Otherwise, if 𝐽(𝑖)>𝐽(𝑖), then set 𝑋(𝑘+1)=𝑗 with probability −𝐽𝑗(𝑖)𝑞=exp−𝐽𝑇(𝑘),where𝑇(𝑘)>0isatemperatureparameter, (A.2) and 𝑋(𝑘+1)=𝑖 with probability 1−𝑞. 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].