The role of residence times in two-patch dengue transmission dynamics and optimal strategies

Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 Q1 16 17 Q2 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 Q3 61 62 63 64 65 66

Contents lists available at ScienceDirect

Journal of Theoretical Biology journal homepage: www.elsevier.com/locate/yjtbi

The role of residence times in two-patch dengue transmission dynamics and optimal strategies Sunmi Lee a, Carlos Castillo-Chavez b a

Department of Applied Mathematics, Kyung Hee University, Yongin-si, Gyeonggi-do 446-701, Republic of Korea Mathematical, Computational, and Modeling Sciences Center, School of Human Evolution and Social Change, Arizona State University, Tempe, AZ 85287-1804, USA

b

H I G H L I G H T S








The two-patch dengue transmission dynamics with the human movement. The human movement is modeled by a residence-time. Optimal control theory is used to identify and evaluate patch-specific control. Targeting intervention on the areas where individuals spend “most” time or where transmissibility is higher. Reducing traffic is likely to take a host–vector system into the world of manageable outbreaks.

art ic l e i nf o

a b s t r a c t

Article history: Received 20 October 2014 Received in revised form 8 February 2015 Accepted 3 March 2015

The reemergence and geographical dispersal of vector-borne diseases challenge global health experts around the world and in particular, dengue poses increasing difficulties in the Americas, due in part to explosive urban and semi-urban growth, increases of within and between region mobility, the absence of a vaccine, and the limited resources available for public health services. In this work, a simple deterministic two-patch model is introduced to assess the impact of dengue transmission dynamics in heterogeneous environments. The two-patch system models the movement (e.g. urban versus rural areas residence times) of individuals between and within patches/environments using residence-time matrices with entries that budget within and between host patch relative residence times, under the assumption that only the human budgets their residence time across regions. Three single-outbreak scenarios are considered: (i) resident hosts in Patch i visit patch j, where i a j but not the other way around, a scenario referred to as unidirectional two-patch system motion; (ii) symmetric bi-directional system motion; and (iii) asymmetric bi-directional patch system motion. Optimal control theory is used to identify and evaluate patch-specific control measures aimed at reducing dengue prevalence in humans and vectors at a minimal cost. Optimal policies are computed under different residence-matrix configurations mentioned above as well as transmissibility scenarios characterized by the magnitude of the basic reproduction number. Due to the complexity of the model, our results are mostly based on numerical simulations. Optimal patch-specific polices can ameliorate the impact of epidemic outbreaks substantially when the basic reproduction number is moderate. The final patchspecific epidemic size variation increases as the residence time matrix moves away from the symmetric case (asymmetry). As expected, the patch where individuals spend most of their time or in the patch where transmissibility is higher tend to support larger patch-specific final epidemic sizes. Hence, focusing on intervention that target areas where individuals spend “most” time or where transmissibility is higher turn out to be optimal. Therefore, reducing traffic is likely to take a host–vector system into the world of manageable outbreaks. & 2015 Elsevier Ltd. All rights reserved.

Keywords: Dengue dynamics Two-patch model Residence-time matrix Optimal patch-specific strategies

1. Introduction Dengue fever, a reemerging vector-borne disease, affects between 50 and 100 million people around the world each year (Halstead et al., 2004). Dengue is transmitted by the vector Aedes aegypti, which carries four different virus serotypes of the genus Flavivirus. Dengue fever is mainly a mild acute febrile disease but it may progress to

severe stages including dengue hemorrhagic fever and/or dengue shock syndrome (World Health Organization, 2013). Increases in severity are typically attributed to secondary infections generated by the delayed interactions of different serotypes, or generic strain-specific variations in virulence, or variations in host susceptibility (Gubler and Kuno, 1997; Homes and Twiddy, 2003). The role of distinct serotypes and/or prior infections on the severity of clinical outcomes

http://dx.doi.org/10.1016/j.jtbi.2015.03.005 0022-5193/& 2015 Elsevier Ltd. All rights reserved.

Please cite this article as: Lee, S., Castillo-Chavez, C., The role of residence times in two-patch dengue transmission dynamics and optimal strategies. J. Theor. Biol. (2015), http://dx.doi.org/10.1016/j.jtbi.2015.03.005i

67 68 69 70 71 72 73 74 75 76 77

2

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

S. Lee, C. Castillo-Chavez / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎

has been established (Endy et al., 2004). Since no effective vaccine is available for dengue fever yet (but see Guy, 2011; Guirakhoo, 2000), disease control and prevention tend to focus on vector control activities that incorporate community participation (Halstead and Deen, 2002; Siqueira et al., 2005) while researchers work on exploring the use value of introducing genetically modified vectors (Sparagano and de Luna, 2008). The sustainability of multiple dengue control measures requires a prolonged and constant efforts by health care and government officials and members of the community, not an easy task. Despite heightened public awareness and instituted intervention efforts, over the past 20 years, outbreaks continue to occur in highly urbanized areas in Central and South America and various other regions in the world (Guzman and Kouri, 2003). Nations with highly effective public health establishments still must face tremendous challenges in dealing with recurrent dengue outbreaks (Burattini et al., 2008; Goh, 1998; Ooi et al., 2006). A. aegypti eradication programs started in the 1950s were stopped in the 1970s in the Americas and as a result, Central and South America have experienced dengue reemergence (Center for Disease Control; Hayes et al., 1996; Morrison et al., 2004; Phillips, 1992; Reiskind et al., 2001). Complex factors contribute greatly to dengue recurrence and its geographical expansion, factors that include population growth, socioeconomic changes, urbanization and transportation (Gubler, 1998, 1997; Harrington et al., 2005; Hayes et al., 1996). Global, regional and local traffic dynamics have contributed to dramatic increases in human mobility with several countries moving from the non-endemic (no virus continuously present) or hypoendemic (one virus present) category to the hyperendemic (multiple virus serotypes co-circulating) class, transitions that have raised global public health efforts to control vector born diseases (Gubler and Kuno, 1997). Sir Ronald Ross pioneered the use of vector born disease transmission models with a multitude of applications and extensions found in the literature (Anderson and May, 1991; Brauer and Castillo-Chavez, 2012; CastilloChavez and Lee, 2013; Favier et al., 2006; MacDonald, 1957; Nishiura, 2006; Ross, 1910). The challenges posed by the incorporation of geographic heterogeneity in the study the spatial dynamics of dengue fever in Peru, Mexico and Brazil have begun to be addressed (see e.g. Chowell et al., 2008, 2007; Favier et al., 2005 and references there in) and, not surprisingly, the role of human movement plays a significant role on disease reemergence and persistence (Adams and Kapan, 2009; Arino et al., 2007; Brauer and Castillo-Chavez, 2012; Cosnera et al., 2009; Martens and Hall, 2000). There are two standard approaches to study the spatial dynamics of infectious disease such as partial differential equations and meta-population models (CastilloChavez et al., 2003; Herrera-Valdez et al., 2011; Murray, 2003; Rodriguez and Tores-Sorando, 2001; Sattenspiel, 2009). There have been various studies to understand spatial and temporal dengue dynamics especially in South America. The recurrent 1990s dengue outbreaks in Peru have been investigated (Hayes et al., 1996; Reiskind et al., 2001). They investigated the dengue outbreaks occurred in Iquitos, the largest city of Peruvian Amazon and the surrounding areas. These results indicated that most infections of residents of the rural villages were acquired while visiting the city of Iquitos. Recently, the 2000s dengue outbreaks in Peru were examined by using two-patch models where jungle areas were always endemic observing how human movement caused epidemics in coast areas (Carlos, 2009). Likewise, the recent dengue outbreaks in Colombia have been explained by human travels (Carrasquilla et al., 2003; Restrepo et al., 2014). They pointed out that many dengue cases were obtained by visiting touristic cities and therefore, control measures should be focused on touristic cities to reduce the spread of dengue cases more effectively. The review on the worldwide dengue transmission has been investigated (Murray et al., 2013). The Global Strategy for Dengue Prevention and Control, 2012–2020 was declared by WHO (World Health Organization, 2013). They have suggested that all possible countermeasures including chemical and non-chemical

controls as well as an effective integrated vector control strategy should target areas of high human–vector contacts in order to reduce the dengue transmission and the resulting disease burden. In this paper, we introduce a two-patch meta-population model to study patterns of dengue spread scenarios between possibly highly distinct environments (e.g. rural versus urban). The patches are coupled by a matrix whose entries represent the proportion of time that residents (under the assumption that only humans can “move”) spend (or have budgeted) “visiting” another patch (Agusto, 2014; Arino et al., 2012; Cosnera et al., 2009; Sattenspiel, 2009). This twopatch model of dengue fever is an extension of the work in Carlos (2009), where a deterministic two-patch model was used to model the impact of uniderectional motions (from Patch 1 to Patch 2). The model extension of Carlos (2009) involves bi-directional motions between patches. Specifically, three single-outbreak scenarios are considered in this paper: (i) resident hosts in Patch i visit Patch j, where i a j but not the other way around, a scenario referred to as unidirectional two-patch system motion; (ii) symmetric bi-directional; and (iii) asymmetric bi-directional patch system motion. We explore the coupling effect on the overall transient dynamics such as the peak size, the peak timing and the final epidemic sizes. Moreover, we vary the level of transmissibility (measured by the local R0i for i¼1, 2) to investigate their impact on the global basic reproductive number, R0. Many applications of optimal control theory have studied the impact of control measures in epidemiological and biological models (Lenhart and Workman, 2007; Lee et al., 2010, 2011, 2012). Rowthorn et al. (2008) have addressed the dilemma that epidemiologists and health administrators must overcome when deciding how to deploy limited resources in a two-patch susceptible-infected-susceptible model for sexually transmitted diseases. Optimal control theory has been used to identify control strategies for vector-borne diseases (Agusto and Lenhart, 2013; Agusto et al., 2012; Blayneh et al., 2009, 2010). Here, we formulate an optimal control problem in a host–vector population living in two-patches and under the assumption that humans budget their time differentially across and between patches. The goal is to identify strategies that minimize the total proportion of infected humans and vectors under cost-effective preventive intervention measures. We explore how patch-specific preventive measures affect dengue spread in both patches by assessing their impact on the transient dynamics including the peak size and the peak timing as well as on the final epidemic sizes under various coupling and transmissibility scenarios. Due to the complexity of the model, our results are mostly based on numerical simulations.

2. Two-patch dengue transmission model One patch dengue model that captures the nonlinear dynamics of humans and vectors has been studied (Carlos, 2009). Since the time scale of interest in this study is short (less than a year) and as a consequence, the demographic dynamics of humans (birth and death rates of the human) are assumed to be negligible. Furthermore, human and vector populations are assumed to be constant, which implies that deaths due to dengue are also assumed to be negligible. Hence, it is understood that all epidemiological variables have been normalized, that is, we deal with population proportions. The variables of the vectors are Sv for the susceptible class, Ev for the incubating class, and Iv for the infected class with Nv  Sv þ Ev þ I v ¼ 1. The variables for the human classes are Sh for the susceptible class, Eh for the incubating class, Ih for the infected class, and Rh for the recovered class with Nh  Sh þEh þ I h þ Rh ¼ 1. Dengue is assumed to be transmitted through two types of interactions: Susceptible vectors may acquire infections through contacts with infected human individuals at the per-infective and per-capita rate βv. Dengue is also transmitted when susceptible individuals are infected via contacts with infected vectors at the per-infective and per-capita rate βh. The incidence rates at which

Please cite this article as: Lee, S., Castillo-Chavez, C., The role of residence times in two-patch dengue transmission dynamics and optimal strategies. J. Theor. Biol. (2015), http://dx.doi.org/10.1016/j.jtbi.2015.03.005i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

S. Lee, C. Castillo-Chavez / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

vectors and humans get infected are therefore βv Sv I h and β h Sh I v ; respectively. The mosquito (or vector) population birth and death rates μ are not ignored but assumed to be equal, in order to support the assumption that the vector population size does not change significantly. Analytical and numerical studies for one patch dengue dynamics can be found in Carlos (2009). A two-patch dengue transmission model is introduced to explore the role of mobility and environmental heterogeneity on the transmission dynamics and control of dengue. Each patch's equations follow the assumptions specified for a single patch model with the local classes being denoted by Svi, Evi, Ivi, Shi, Ehi, Ihi and Rhi, for Patch i for i ¼ 1; 2. The patches are coupled via the resident budgeting time matrix P ¼ ðpij Þ22 for i; j ¼ 1; 2, where each pij is a constant in ½0; 1 P and 2j ¼ 1 pij ¼ 1 for i ¼ 1; 2. More precisely, the coupling parameter pij represents the proportion of time that a person residing in Patch i budgets and spends in Patch j. It is assumed that vectors do not relocate, that is, only humans can transit across patches. Therefore, the P expression, 2j ¼ 1 pij ¼ 1 holds for only humans (Shi and Ehi). The two patch dynamics are captured by the following patch-specific system of nonlinear ordinary differential equations: S_ vi ¼ μ  βvi ð1  ui ÞSvi

2 X

2 X

pji I hj  μEvi  κ Evi

j¼1

I_vi ¼ κ Evi  μI vi S_ hi ¼  Shi

2 X

βhj ð1  u j Þpij Ivj

j¼1

E_ hi ¼ Shi

2 X

method to compute the basic reproduction number, R0 (Brauer and Castillo-Chavez, 2012; Diekmann and Heesterbeek). The global basic reproductive number for the two patch system is given as sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffi κ ð2Þ ðϕ1 þ ϕ2 Þ R0 ¼ 2μγ ðκ þ μÞ where

ϕ1 ¼ p211 βh1 βv1 þ p221 βh2 βv1 þp212 βh1 βv2 þ p222 βv2 βh2 ϕ2 ¼ ðp211 βh1 βv1 þ p221 βh2 βv1 Þ2 þðp212 βh1 βv2 þ p222 βv2 βh2 Þ2 2 2 þ2β v1 β v2 ðp211 p212 β h1 þ p221 p222 βh2 Þ þ2β h1 βh2 β v1 β v2 ð4p11 p12 p21 p22  p211 p222  p212 p221 Þ Uncoupling the system, that is, taking p11 ¼ 1, p22 ¼ 1, gives qffiffiffiffiffiffiffiffiffiffiffiffiffiffi βv1 βh1 κ what we refer to as the local reproductive numbers: R01 ¼ μγ ðκ þ μÞ qffiffiffiffiffiffiffiffiffiffiffiffiffiffi β β κ for Patch 1 and R02 ¼ μγv2ðκ þh2μÞ for Patch 2 (details are shown in Appendix). 2.1. Two-patch transmission dynamics under different residencetime configurations

pji I hj  μSvi

j¼1

E_ vi ¼ βvi ð1  ui ÞSvi

3

βhj ð1  u j Þpij I vj kEhi

j¼1

I_hi ¼ kEhi  γ I hi R_ hi ¼ γ I hi :

ð1Þ

for i¼1, 2. This two-patch model has been previously studied for the case of a unidirectional case (from Patch 1 to Patch 2) (Carlos, 2009). This paper revisits the results and extends the model to more general scenarios including bi-directional motions. A susceptible vector in Patch i can be infected by an infected person from Patch i as well as by an infected person visiting Patch i from Patch j. Susceptible vectors in Patch i get infected through contacts at the per capita rate βvi per infected (host) individual. Susceptible individuals in Patch i get infected through contacts at the per capita rate βhi per infected vector. Hence, the rate at which vectors and humans get infected in Patch 1 is given by βv1 ðp11 I h1 þ p21 I h2 ÞSv1 and ðβh1 p11 Iv1 þ βh2 p12 I v2 ÞSh1 ; respectively. Similarly, the incidence rates in Patch 2 are β v2 ðp12 I h1 þ p22 I h2 ÞSv2 and ðβh1 p21 I v1 þ βh2 p22 I v2 ÞSh2 . An optimal control problem of interest is formulated through the incorporation of the patch-specific control functions ui(t) as ð1 ui ðtÞÞ in the transmission rates for each patch (the detailed descriptions are given in the next section). The basic reproductive number, R0, is the average number of secondary infectious cases when one infectious individual is introduced in a whole susceptible population. For the two-patch model, it is necessary to use the next generation operator

We explore the role of the residence-time matrix on the transmission dynamics, in particular, on the human incidences. We consider three distinct coupling scenarios. First, one-way coupling when people visit only from Patch 1 in Patch 2 (p22 ¼1). Second, symmetric coupling where people follow symmetrical patterns in each patch ðp11 ¼ p22 Þ. Thirdly, asymmetric coupling when visitors from each patch spend their visiting time in one patch more than the other. Therefore, we have p11 o p22 and as j p11  p22 j gets larger, the coupling becomes more asymmetric. We assume that p11 ; p22 A ½0:5; 1 with p12 ¼ 1  p11 and p21 ¼ 1  p22 , which implies that people spend more time (or not less) in their home patch than other patch. Furthermore, we consider different coupling intensities, (i) “uncoupled”(P1), (ii) weakly coupled (P2) and (iii) strongly coupled (P3). When the patches are “weakly” coupled, that is, a small value for p12 ðp11 ¼ 0:999; p12 ¼ 0:001Þ. While with “strong” coupling, visitors from Patch 1 spend quite an amount of time in Patch 2 ðp11 ¼ 0:6; p12 ¼ 0:4 c 0Þ. For all simulations (unless we mention), we use the initial conditions, Sh1 ð0Þ ¼ 0:999, I h1 ð0Þ ¼ 0:001, Sh2 ¼ Sv1 ð0Þ ¼ Sv2 ð0Þ ¼ 1 and the rest of initial conditions are zero. Other parameter values are found in Table 2. The epidemic starts in Patch 1 while Patch 2 gets infected as a result of coupling. Again, since only human can move between patches, the conditions mentioned above for the residence-time matrix hold for only humans in Shi and Ehi: P1 ¼

p11

p12

p21

p22

!

 ¼

1

0

0

1

 ð3Þ

Fig. 1 shows the infected and cumulative proportions of human incidence under the unidirectional case varying the coupling intensity: uncoupled (solid), weak coupling (dashed p11 ¼ 0:99) and strong coupling (dotted p11 ¼ 0:9). This illustrates how the coupling intensity of patches affects the patch-specific transient dynamics such as the peak time, the peak size (top) and final epidemic sizes (bottom). Here, the peak size is the maximum incidence and the

Table 1 Parameter values in the residence-time matrix. Coupling intensity

P2 (weak coupling)

P3 (strong coupling)

One way coupling Symmetric coupling Asymmetric coupling

p11 ¼ 0.999, p12 ¼ 0.001, p21 ¼ 0, p22 ¼ 1 p11 ¼ 0.99, p12 ¼ 0.01, p21 ¼ 0.01, p22 ¼ 0.99 p11 ¼ 0.9, p12 ¼ 0.1, p21 ¼ 0.001, p22 ¼0.999

p11 ¼ 0.9, p12 ¼0.1, p21 ¼0, p22 ¼1 p11 ¼ 0.6, p12 ¼0.4, p21 ¼0.4, p22 ¼0.6 p11 ¼ 0.6, p12 ¼0.4, p21 ¼0.001, p22 ¼ 0.999

Please cite this article as: Lee, S., Castillo-Chavez, C., The role of residence times in two-patch dengue transmission dynamics and optimal strategies. J. Theor. Biol. (2015), http://dx.doi.org/10.1016/j.jtbi.2015.03.005i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

S. Lee, C. Castillo-Chavez / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎

4

Table 2 Parameter definitions and baseline values used in numerical simulations. Parameter Description

βhi k κ γ μ tf I hi ð0Þ b Wi

Infected host proportion

One−way coupling

Patch 1

−3

15

Value

Probability person from Patch i is visiting Patch j 0–1 Per capita transmission rate from host to vector in Patch i (days  1) 0.15– 0.6 Per capita transmission rate from vector to host in Patch i (days  1) 0.15– 0.6 Host rate of progression from latent to infectious (days  1) 1/5.5 Vector rate of progression from latent to infectious (days  1) 1/5.5 Host recovery rate (days  1) for infectious class (days  1) 1/5.5 Vector per capita vector birth/death rate (days  1) 1/10.5 The simulated duration (days) 300 The initial values of infected host (i¼1,2) 0.001 The upper bound of control (efficient reduction rates, days  1) 0.35 Weight constants on controls i¼ 3,4 0.1–2

pij βvi

Cumulative infected proprotion

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

x 10

P1 P2

10

P3

5 0

15 10 5

0

100

200

300

0

1

1

0.5

0.5

0

Patch 2

−3

x 10

0

100

200

300

Time (days)

0

0

100

200

300

0

100

200

300

Time (days)

Fig. 1. Incidence and cumulative infected for human are displayed under the unidirectional motion using the local basic reproduction number, R0i ¼ 1:8, i ¼ 1; 2. The impact of the coupling intensity in the residence-time matrix on the transient dynamics is investigated: uncoupled (P1), weakly coupled (P2) and strongly coupled (P3) given in Table 1.

peak time is the time where the maximum incidence occurs. The peak time in Patch 2 is heavily dependent on the coupling intensity. Weak coupling (dashed p11 ¼ 0:99) causes the largest delay of the epidemic peak of Patch 2 and also the peak size in Patch 1 is larger. However, with strong coupling, the peak of Patch 2 occurs even earlier than Patch 1's peak (dotted curves) and the peak size in Patch 2 is larger. As expected, the final epidemic sizes in both patches are same since the local basic reproductive numbers R01 ¼ R02 are same. Fig. 2 shows the results of two-way coupling: symmetric (top) and asymmetric (bottom) under different coupling intensities. The peak size in two patches are the same for the symmetric coupling case (top) regardless of coupling intensities. Like the unidirectional case, as the intensity increases, the peak time in Patch 2 occurs faster. On the other hand, the results of asymmetric coupling (bottom) indicate that the peaks of Patch 2 occur earlier (or not later) than Patch 1 (see dashed and dotted curves) and Patch 2 has the larger peak size regardless of coupling intensities due to the fact that Patch 2 is the dense-patch. Under all coupling cases, as the coupling becomes weaker, that is, people from Patch 1 spent most of their time in Patch 1, the peak to peak distance (peaks in Patch 2) increases greatly. When the residence-time matrix becomes more asymmetric, the epidemics in both patches occur closer and eventually Patch 2 peaks

earlier than Patch 1 despite the epidemic starting in Patch 1 and also the dense-patch has the larger peak size. 2.2. The global basic reproduction number, R0 and final epidemic size We present the impact of the local R0i (i¼ 1, 2) on the global R0 under different coupling cases. Fig. 3 displays the global R0 given in (2) as functions of R01 and R02. For the uncoupled case, the result shows that R0 ¼max fR01 ; R02 g (top left). For one way coupling R02 slightly dominates R01 since R0 ¼ R02 when R02 4R01 , but R0 o R01 when R01 4 R02 (top right). For both cases, results indicate that R0 r maxfR01 ; R02 g. However, for the two-way coupling cases, the value of R0 is slightly reduced for symmetric coupling (bottom left) and slightly increased for asymmetric coupling (bottom right). It is worth to notice that the global R0 gets larger as coupling becomes more asymmetric. Fig. 4 illustrates the impact of the residence-time matrix component on the global basic reproductive number R0, which is displayed as a function of the p11 under three different values of p22. Now, the local basic reproductive number in each patch is fixed as 1.8 in the left and as 2.0 in the left. For both cases, the global basic reproductive number R0 gets the largest value at p11 ¼ 0:5 for p22 ¼ 1 while R0 gets the largest at p11 ¼ 1 for p22 ¼ 0:5 and p22 ¼ 0:75. Also, in both

Please cite this article as: Lee, S., Castillo-Chavez, C., The role of residence times in two-patch dengue transmission dynamics and optimal strategies. J. Theor. Biol. (2015), http://dx.doi.org/10.1016/j.jtbi.2015.03.005i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

S. Lee, C. Castillo-Chavez / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎

Symmetric coupling Patch 1

Cumulative infected proprotion

Infected host proportion

−3

15

x 10

P1 P2

10

P3

5 0

15 10 5

0

100

200

300

0

1

1

0.5

0.5

0

Patch 2

−3

x 10

0

100

200

300

0

0

100

200

300

0

100

200

300

Time (days)

Time (days)

Asymmetric coupling Patch 1

Infected host proportion

−3

Cumulative infected proprotion

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

5

15

x 10

Patch 2

−3

P1

15

x 10

P

2

10

P

10

3

5

0

5

0

100

200

300

0

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

100

200

300

Time (days)

0

0

100

200

300

0

100

200

300

Time (days)

Fig. 2. Incidence and cumulative infected for human are displayed under bi-directional motions using the local basic reproduction number, R0i ¼ 1:8, i ¼ 1; 2: (i) symmetric coupling (top panel) and (ii) asymmetric coupling (bottom panel). For each case, three different coupling intensities are considered: uncoupled (P1), weakly coupled (P2) and strongly coupled (P3) given in Table 1.

graphs, R0 has the minimum value when p11 ¼ p22 . This reconfirms that the global R0 gets larger as coupling becomes more asymmetric (or j p11  p22 j becomes larger) as seen in Fig. 3. Next, we evaluate how coupling affects the patch-specific final epidemic size. Fig. 5 highlights the final epidemic size as a function of local basic reproduction numbers R0i (i¼ 1, 2) under three distinct coupling scenarios as described above. For all cases, obviously, when R0i o 1 for i¼ 1, 2, there is no outbreak in both patches. For one-way coupling, since coupling is only from Patch 1 to Patch 2, the final epidemic size in Patch 2 is larger than Patch 1 (blue solid curves). This become more clear for asymmetric coupling, the red triangle curves show the largest final epidemic sizes in both patches while the symmetric coupling gives the least final epidemic size among different coupling scenarios. Again, this is consistent with the results for the global R0 explained in Figs. 3 and 4.

3. Two-patch model with prevention interventions An optimal control problem is formulated in the two-patch model by incorporating the patch-specific control functions ð1  ui ðtÞÞ into the incidence rates which humans and vectors get infected for Patch i (i¼ 1, 2) in the state system (1). The preventive control efforts may involve the application of pesticide (sprays) or mosquito repellents, netting, window screens or the reduction of the impact of vector breeding grounds or the result of education campaigns that increases personal protection. We assume that these preventive interventions do not reduce the total vector population significantly. Also, we assume that the effect of these interventions implicitly translates in reductions of contacts between vectors and humans per unit of time. Specifically, incidence rates including controls are modified as P P βvi ð1  ui ðtÞÞSvi 2j ¼ 1 pji I hj and Shi 2j ¼ 1 βhj ð1  uj ðtÞÞpij Ivj for i¼1,2.

Please cite this article as: Lee, S., Castillo-Chavez, C., The role of residence times in two-patch dengue transmission dynamics and optimal strategies. J. Theor. Biol. (2015), http://dx.doi.org/10.1016/j.jtbi.2015.03.005i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

S. Lee, C. Castillo-Chavez / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎

6

One way

Uncoupled 3.5 3.5 3 3 2.5

R0

R0

2.5

2

2 1.5

1.5

1 3 2.5

1 3

3 2.5

2

2

2

1.5 R02

1.5

1 1

R02

R01

1

Symmetric

3

2.5

2 R01

1.5

1

Asymmetric 3

3.5

3

2.8 3

2.6

2.5

2.4

R0

2.5 R0

2

2.2 2

1.5

2 1.8

1.5

1.6 1 3

2.5

2

1.5

1

R02

2

1.5

1

2.5

1 3

3

1.4 2 R02

R01

1.5

1 1

2.5

2

3

1.2

R01

Fig. 3. The global basic reproduction number, R0 is plotted as a function of the local R0i for i ¼ 1; 2 under four different residence-time configurations. The right on the bottom shows the largest values of R0 for the asymmetric residence-time matrix.

2.15

2.35

2.1

2.3

2.05

2.25 The global R0

The global R0

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

2

1.95 1.9 1.85

2.1 2.05

1.8 1.75 0.5

2.2 2.15

2 0.6

0.7

0.8

0.9

p11

1

1.95 0.5

0.6

0.7

0.8

0.9

1

p11

Fig. 4. The global basic reproductive number R0 is displayed as a function of the p11 under three different values of p22 (the component of the residence-time matrix). The local basic reproduction number R0i is fixed as R0i ¼ 1:8 on the left and R0i ¼ 2:0 on the right for i¼ 1,2.

If ui ðtÞ ¼ 0, there is no control and as ui ðtÞ-1 (ideal efforts or perfect preventions), we see that the transmission rate goes to zero. The goal of the optimal problem is to minimize the infected proportions of humans and vectors in both patches at a minimal cost of implementation over a finite time horizon. The objective functional to be minimized is given by Z tf Jðu1 ðtÞ; u2 ðtÞÞ ¼ W 1 ðI h1 ðtÞ þ I v1 ðtÞÞ þ W 2 ðI h2 ðtÞ þ I v2 ðtÞÞ 0

þ 12 W 3 u21 ðtÞ þ 12 W 4 u22 ðtÞ dt where W1, W2 are weight constants on the infected human and vector for Patch 1 and Patch 2, respectively. W3 and W4 are the weight constants for the prevention effort or the relative cost of the

implementation of the preventive control for Patch 1 and Patch 2, respectively. We seek an optimal pair (Un, Xn) such that JðU n Þ ¼ minfJðUÞj U A Ωg;

ð4Þ

where Ω ¼ fðu1 ðtÞ; u2 ðtÞÞ A ðL ð0; t f ÞÞ J a r ui ðtÞ rb; t A ½0; t f ; i ¼ 1; 2g subject to the state system (1) with X ¼ ðSv1 ; Ev1 ; I v1 ; Sh1 ; Eh1 ; I h1 ; Rh1 ; Sv2 ; Ev2 ; I v2 ; Sh2 ; Eh2 ; I h2 ; Rh2 Þ and U ¼ ðu1 ; u2 Þ. The existence of optimal controls is guaranteed from standard results on optimal control theory (Fleming and Rishel, 1975). Pontryagin et al.'s (1962) Maximum Principle is used to establish necessary conditions that must be satisfied by an optimal solution. Derivations of the necessary conditions are shown in Appendix. The standard scheme (a two point boundary method Lenhart and Workman, 2007) is employed to find numerical solutions to (4). First, 1

2

Please cite this article as: Lee, S., Castillo-Chavez, C., The role of residence times in two-patch dengue transmission dynamics and optimal strategies. J. Theor. Biol. (2015), http://dx.doi.org/10.1016/j.jtbi.2015.03.005i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

S. Lee, C. Castillo-Chavez / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎

level of control effort or the effectiveness of control measures. It is not easy to estimate the accurate control upper bound since it strongly depends on various economic, environmental and geographical factors. In this optimal formulation, we assume that the best possible reduction of transmission by implementing preventive control is 35% (i. e. b¼0.35 is used) (Gubler and Kuno, 1997). The numerical sensitivity analysis for these two parameters is carried out in the last subsection.

the state system (1) is solved forward in time with initial conditions and an initial guess for the control. Second, the adjoint system with transversality conditions is solved backward in time. Third, the optimality condition is updated and finally, the three steps above are iterated until convergence is achieved. Parameter values used are collected in Table 2 with the initial conditions, Sh1 ð0Þ ¼ 0:999, I h1 ð0Þ ¼ 0:001, Sh2 ð0Þ ¼ 1, Svi ð0Þ ¼ 1 for i¼1,2 and the rest of them are zero. Since our main interest is a single outbreak, the final time is taken as 300 days. There are two critical parameter values which affect optimal solutions greatly (control weight constants and control upper bounds). For weight constants of the infected humans and infected vectors, W 1 ¼ W 2 ¼ 1 are used throughout our simulations in the objective functional (3). Next, the control weight constants W 3 ¼ W 4 represent the relative cost of implementation of preventive controls. As we increase (decrease) them, the relative cost of control increases (decreases) therefore, the magnitude of the optimal control functions decreases (increases). We use W 3 ¼ W 4 ¼ 0:1 throughout our simulations unless we mention. Another important control parameter value is the control upper bound, b, which represents the maximum

3.1. Implications of prevention strategies under three distinct coupling We present the numerical simulations associated with implementing optimal prevention control functions as well as their effect on human incidence under different coupling parameters. Figs. 6, 7 and 8 compare incidence and cumulative curves in the presence of controls (dashed curves) with the results in the absence of controls (solid curves). We consider the first case where the epidemic first starts in Patch 1 and moves to Patch 2 (one-way coupling). Fig. 6 shows the proportion of incidence and cumulative infected for human and optimal control functions. Optimal preventive efforts should be implemented with the maximum level in Patch 1 and with the moderate level in Patch 2. There is almost no outbreak under control in Patch 2 while Patch 1 starts a slight outbreak. For one-way coupling, as the coupling becomes stronger, the level of controls is increased in Patch 2. The results of symmetric and asymmetric coupling are displayed in Figs. 7 and 8, respectively. Fig. 7 displays the results for symmetric coupling. It is obvious that optimal prevention efforts in both patches require the maximum level in Patch 1 and Patch 2 and optimal interventions can almost stop epidemics. Fig. 8 shows the results for asymmetric coupling. Again, due to asymmetric coupling, optimal prevention efforts should be implemented with the maximum level in both patches during the entire time window. There is almost no outbreak under controls in Patch 1 and Patch 2. Regardless of coupling scenarios, under the implementation of patch-specific optimal controls can “almost” stop an outbreak and lead to dramatic reductions in the human incidence (both patches) for the moderate value of R0 (this case R0i ¼ 1:8 for i ¼ 1, 2). As pointed in the previous section, more asymmetric coupling results in the larger final epidemic size, therefore, it requires more intensive controls in both patches. We evaluate how patch-specific optimal interventions affect the final epidemic size. In Fig. 9, the final epidemic sizes in both patches are plotted under three different coupling cases compared with the ones in the absence of controls (solid curves). For all cases, there is a critical value of Rc (the controlled reproduction number Brauer and Castillo-

Final epidemic size

Patch 1

1

0.5

0

1

1.2

1.4

1.6

1.8

2

R

01

1 Patch 2

One way Symmetric Asymmetric

0.5

0

1

1.2

1.4

1.6

1.8

2

R02 Fig. 5. The impact of the residence-time matrix on the final epidemic size is illustrated in the absence of controls. The final epidemic size is displayed as a function of R01 ¼ R02 under three different residence-time configurations. The total final epidemic size becomes larger as the residence-time matrix becomes more asymmetric (the red triangle curve shows the largest final epidemic size in both patches). (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

Infected human proportion Patch 1

0.015

No control

Cumulative infected proportion 1

Control effort 0.4

With control

0.01

0.5

u1(t)

0.2

0.005 0

0

100

200

300

0.015 Patch 2

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 Q6 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

7

0

0

200

300

1 No control With control

0.5 0.005 0

100

200

Time (days)

300

0

0

0

100

200

300

0.4

0.01

0

100

0

100

200

Time (days)

0.2

u (t) 2

300

0

0

100

200

300

Time (days)

Fig. 6. Incidence, cumulative infected for human and optimal patch-specific controls are displayed under the unidirectional motion using the local basic reproduction number, R0i ¼ 1:8, i ¼ 1; 2 and the residence-time matrix P2 (weakly coupled) in Table 1.

Please cite this article as: Lee, S., Castillo-Chavez, C., The role of residence times in two-patch dengue transmission dynamics and optimal strategies. J. Theor. Biol. (2015), http://dx.doi.org/10.1016/j.jtbi.2015.03.005i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

S. Lee, C. Castillo-Chavez / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎

8

Infected human proportion Patch 1

0.015

Cumulative infected proportion 1

No control

0.5

u (t) 1

0.2

0.005 0

100

200

0

300

0.015 Patch 2

Control effort 0.4

With control

0.01

0

0

100

200

300

0

1

0.4

0.5

0.2

0

100

200

300

u (t)

0.01 No control With control

0.005 0

0

100

200

0

300

0

100

Time (days)

200

300

0

2

0

100

Time (days)

200

300

Time (days)

Fig. 7. Incidence, cumulative infected for human and optimal patch-specific controls are displayed under the symmetric bi-directional motion using the local basic reproduction number, R0i ¼ 1:8, i ¼ 1; 2 and the residence-time matrix P2 (weakly coupled) in Table 1.

Infected human proportion

Patch 1

0.015

Cumulative infected proportion 1

No control

Control effort 0.4

With control

0.01

0.5

0.2

u (t) 1

0.005 0

0

100

200

0.015 Patch 2

300 No control

0

0

100

200

300

0

1

0.4

0.5

0.2

0

100

200

With control

0.01

300

u (t) 2

0.005 0

0

100 200 Time (days)

300

0

0

100 200 Time (days)

300

0

0

100 200 Time (days)

300

Fig. 8. Incidence, cumulative infected for human and optimal patch-specific controls are displayed under the asymmetric bi-directional motion using the local basic reproduction number, R0i ¼ 1:8, i ¼ 1; 2 and the residence-time matrix P2 (weakly coupled) in Table 1.

3.2. Implications of prevention strategies under different transmissibility We also explore the case where two patches have a different level of transmissibility. The results for the one-way coupling case are

Final epidemic size 1 0.8 Patch 1

Chavez, 2012) where an epidemic occurs in both patches. As we can see, there are huge reductions in the final epidemic size of both patches under the implementation of control measures. For small R0i o ðRc Þ, the implementation of optimal control significantly prevents the epidemic from occurring regardless of coupling in both patches. Above an R0i 4 Rc , the final epidemic size starts to increase, eventually equaling the final epidemic size with no control when R0i ð ¼ 3Þ. For the one-way coupling, the value of Rc is  1:6 for both patches (see dashed curves). For the symmetric coupling case, Rc  1:6 for Patch 1 and Rc  1:7 for Patch 2 (see dotted curves). For the asymmetric coupling case, since there are higher proportions of human moving in between two patches, the value of Rc is getting smaller Rc  1:5 for both patches (see triangle curves). It is clear to observe significant reductions in the total final epidemic sizes when the optimal prevention controls are implemented regardless of coupling. Again, as coupling becomes more asymmetric, the final epidemic size gets larger, hence, the value of Rc gets smaller.

0.6

No control (uncoupled) One way Symmetric Asymmetric

0.4 0.2 0 1

1.2

1.4

1.6

1.8

2

2.2

2.4

2.6

2.8

3

2.2

2.4

2.6

2.8

3

R

01

1 Patch 2

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

0.5

0 1

1.2

1.4

1.6

1.8

2 R02

Fig. 9. The impact of the residence-time matrix on the final epidemic size is illustrated in the presence of controls. The final epidemic size is displayed as a function of R01 ¼ R02 under four different cases (Patch 1 on top and Patch 2 on bottom). Optimal patch-specific policies can bring the final epidemic size down substantially (The curves of triangle, dashed, starred curves in the presence of controls are compared with the solid black curve in the absence of controls).

Please cite this article as: Lee, S., Castillo-Chavez, C., The role of residence times in two-patch dengue transmission dynamics and optimal strategies. J. Theor. Biol. (2015), http://dx.doi.org/10.1016/j.jtbi.2015.03.005i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

S. Lee, C. Castillo-Chavez / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎

unidirectional case). For the bi-directional cases, it is noted that the impact of the optimal prevention controls becomes less efficient as the transmissibility becomes more asymmetric in two patches, in particular, when the residence matrix becomes more asymmetric. Hence, the results indicate that where the outbreak starts and how fast a disease spreads is critical to determine optimal strategies to prevent a large outbreak. This implies that reducing traffics between patches can lead to more manageable levels of dengue outbreaks.

shown here under two scenarios: First, Patch 1 has a higher transmissibility than Patch 2 and secondly, Patch 2 is higher than Patch 1. Figs. 10 and 11 compare incidence and cumulative curves in the presence of controls (dashed curves) with the results in the absence of controls (solid curves). In Fig. 10, using R01 ¼ 2 for Patch 1and R02 ¼ 1:2 for Patch 2, we observe that despite implementing the maximum level of control throughout the whole time window in Patch 1, there is still a large outbreak in Patch 1 (see dashed curves in Patch 1: the peak time is delayed and the peak size is reduced greatly) while Patch 2 is under control. Fig. 11 displays the results using R01 ¼ 1:2 for Patch 1 and R02 ¼ 2:0 for Patch 2. Clearly, (almost) the maximum level of control is required in Patch 2 while the decreasing control is required in Patch 1. In this case, the control prevents an outbreak in both patches even though Patch 2 has higher transmissibility. When an epidemic starts at Patch 1 with a higher transmission rate, even the maximum level of control cannot stop the outbreak. However, the patch originates an epidemic with a lower transmission rate, it is still manageable to stop the outbreak (for this

Infected human proportion 0.015

3.3. Sensitivity analysis As we mentioned earlier, there are two critical parameter values which affect optimal solutions greatly (control weight constants and control upper bounds). We present the impact of control weight constants and control upper bounds on patch-specific final epidemic sizes under the asymmetric case. The left two panels in Fig. 12 show the final epidemic size as a function of R0i under three different values of weight constants ðW 3 ¼ W 4 ¼ 0:1; 1; 2Þ. The weight constant repr-

Cumulative infected proportion 1

Control effort 0.4

Patch 1

No control With control

0.01

0.5

u (t)

0.2

1

0.005 0

0

100

200

300

0

0

100

200

300

0

0

100

200

300

−3

x 10

Patch 2

1

0.06

0.4

0.04 0.5

u2(t)

0.2 0.02

0

0

100

200

300

0

0

Time (days)

100

200

300

0

0

Time (days)

100

200

300

Time (days)

Fig. 10. Incidence, cumulative infected for human and optimal patch-specific controls are displayed under the unidirectional motion using the local basic reproduction number, R01 ¼ 2:0 4R02 ¼ 1:2 and p11 ¼ 0:999, p22 ¼ 1.

Infected human proportion

Cumulative infected proportion

Control effort

−3

Patch 1

2

x 10

0.2

0.4

0.15

0.3

1

0.1

0.2

0.5

0.05

0.1

1.5

0 0

u (t)

No control With control

100

200

300

0.015 Patch 2

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

9

0

0

100

200

300

1

0

100

200

300

0.4 0.3

0.01 0.5

0.2

0.005 0 0

0

1

u2(t)

0.1 100

200

Time (days)

300

0

0

100

200

Time (days)

300

0

0

100

200

300

Time (days)

Fig. 11. Incidence, cumulative infected for human and optimal patch-specific controls are displayed under the unidirectional motion using the local basic reproduction number, R01 ¼ 1:2 o R02 ¼ 2:0 and p11 ¼ 0:999, p22 ¼ 1.

Please cite this article as: Lee, S., Castillo-Chavez, C., The role of residence times in two-patch dengue transmission dynamics and optimal strategies. J. Theor. Biol. (2015), http://dx.doi.org/10.1016/j.jtbi.2015.03.005i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

S. Lee, C. Castillo-Chavez / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎

10

esents the relative cost of the control implementation and a smaller (larger) value means a lower (higher) relative cost. If the relative cost is low, the control level increases, which would keep the epidemic from growing for a moderate value of the basic reproductive number (see dashed curves in the left panel using W 3 ¼ W 4 ¼ 0:1). Higher weight constants mean that the relative cost of control is higher, thus the control level decreases, which leads to larger outbreaks in both patches. For a higher cost, the final epidemic size increases greatly compared with the result for a lower cost (see triangle curves using W 3 ¼ W 4 ¼ 2). The lower the relative cost, the epidemic can be controlled more significantly in both patches. Next, we evaluate the impact of different control upper bounds on the patch-specific final epidemic sizes under the asymmetric case. The constant b represents the maximum level of efficiency of the control implementation and three different control upper bounds are used (b ¼ 0.2, 0.35, 0.5). The right panels in Fig. 12 show the final epidemic size as a function of R0i. The smaller (larger) control upper bound implies the less (more) efficient controls. If the control upper bound is smaller, the constant control (with the maximum level) is required in both patches, however, even though the maximum level of control is implemented in both patches throughout the whole time period, there are still large outbreaks (see dashed curves using b¼0.2). Larger upper bounds (in this case, when b 4 0:35) implies that the effectiveness of control is higher, thus it is able to keep the epidemic from growing in both patches (see dotted or triangle curves using b¼ 0.35 or 0.5). Therefore, the final epidemic size can be reduced dramatically in both patches for a moderate value of R0 (up to around 2) when the effectiveness of control measures is large enough (b¼0.5).

4. Discussion We have studied the dynamics of dengue transmission in a twopatch dengue model in which the patches are connected via a residence-time matrix. The model proposed here could represent two locations (urban versus rural areas) that have a constant and well-defined visiting relationship. We assume that the epidemic starts in Patch 1 and then spreads to Patch 2 only by human visitors. Three coupling scenarios are considered: one-way coupling, symmetric coupling and asymmetric coupling. First, the local reproductive numbers are computed when the patches are decoupled. Then, the impact of these coupling on the global R0 for the patch system is illustrated by varying each patch's uncoupled local R0i for i ¼ 1, 2. As

No control W =W =0.1

0.5

3

4

W3=W4=1 W3=W4=2 1.5

2

2.5

0.5

No control b=0.2 b=0.35 b=0.5

0

0 1

Final epidemic size

1

Patch 1

Patch 1

coupling becomes asymmetric, the value of R0 increases, therefore, the value of R0 for asymmetric coupling is the largest while symmetric coupling leads to the minimum global R0. This is also consistent with the final epidemic size (more asymmetry leads to larger final epidemic sizes). For all coupling scenarios, when the probability (that one infected individual who resides in Patch i is visiting Patch j) is higher, the epidemic peaks in both patches occur closer and eventually, Patch 2 peaked earlier than Patch 1. In particular, asymmetric coupling results in the peak of Patch 2 occurring earlier than Patch 1, which in turn, causes a larger final epidemic size in Patch 2 compared to Patch 1. This can be critical because one has to be careful when determining where an epidemic starts under this uncertainty. We have developed an optimal control framework to access how patch-specific optimal control strategies affect the spread of dengue in a coupled two-patch model. First, we identify the optimal strategy to prevent or mitigate an epidemic under three distinct coupling scenarios and compare with the results in the absence of controls with the same transmission rate in both patches. Results show that the transient details including time and magnitude of the epidemic's peak as well as control measures are different for each coupling case. A more symmetric coupling matrix results in the smaller final epidemic size (or smaller global R0). Therefore, the optimal strategies under symmetric coupling require the minimum level of controls among three coupling cases. In one-way coupling, stronger coupling requires more intensive controls focused on Patch 2. For two-way coupling with more asymmetric coupling, optimal prevention efforts should be implemented with the maximum level in both patches. With a diffe; rent transmission rate, an epidemic in Patch 1 (with a higher transmission rate than Patch 2) cannot be prevented despite the maximum level of control employed. However, the patch which originates an epidemic with a lower transmission rate can have success with a manageable elimination of the outbreak. These results highlight that the start of an outbreak and rate of disease spread between patches are critical to determine the optimal strategies necessary to prevent a large outbreak. The impact of coupling on the final epidemic size as a function of R0 is illustrated. Overall, there is a certain threshold (Rc, the controlled reproductive number) where the epidemic cannot take off under the implementation of optimal control. In this threshold region ½1; Rc , the optimal control is successful to prevent the outbreak, however, above Rc, the epidemic starts to grow and no control can stop the epidemic. This threshold depends on the residence-time matrix. As

Final epidemic size

1

1

3

1.5

2

2.5

3

2.5

3

R01

R01 1 Patch 2

1

Patch 2

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

0.5

0.5

0

0 1

1.5

2 R

02

2.5

3

1

1.5

2 R02

Fig. 12. The impact of the control parameter values on the final epidemic size is illustrated. The final epidemic size is displayed as a function of the local basic reproductive number under different values of the control weight (left panel) and the control upper bound (right panel).

Please cite this article as: Lee, S., Castillo-Chavez, C., The role of residence times in two-patch dengue transmission dynamics and optimal strategies. J. Theor. Biol. (2015), http://dx.doi.org/10.1016/j.jtbi.2015.03.005i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

S. Lee, C. Castillo-Chavez / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

coupling becomes more asymmetric, the window of R0 which can be controlled completely by optimal preventions gets smaller. In general (the overall profiles of controls regardless of coupling) for a smaller value of R0 o 1:5, the constant control is required for a longer period of time. This is because small R0 produces smaller, but long lasting epidemics that could last longer that the time control measures are implemented. For a larger value of R0 generates a larger and earlier epidemic peak therefore, the control strategy takes the maximum level (intensive level) in a shorter (earlier) time of period. We take into account the importance of efficiency in implementation of the control by means of the relative cost (weight constant) and the effectiveness (upper bound). If the implementation of control is not too costly (weight constants), the control can stop the outbreak for a moderate value of R0 ¼ 1.8. Also, sensitivity results for control upper bounds show that the maximum level or efficiency is required (in the range of 0.35) to prevent the outbreak for a moderate value of R0 ¼1.8. In this formulation, if the control parameters lie in the reasonable range (W 3 ; W 4 r 1 and b Z 0:35), then the epidemic can be prevented significantly in both patches. These results apply directly to “real world” issues like a government's ability to control and epidemic strongly depends on its resources and efficiency. If there are enough resources, control measures can be applied to both patches, and if resources are more scarce, it is better to directly target the location of interest. Extensive simulations show that optimal patch-specific interventions can limit the severity of outbreaks when the a priori R0 is brought below a certain threshold (the controlled reproduction number, Rc). For the higher values of R0, the use of optimal-cost effective intervention turns out to be insufficient to prevent large outbreaks. Note that this threshold (Rc) is dependent on the residence-time matrix. As the resident time matrix becomes more asymmetric (or time spent from each part becomes more asymmetric), Rc gets smaller which in turn leads to large outbreaks. Furthermore, a high-dense patch resulting from elevated time spent by individuals or increased transmissibility may have a larger final epidemic size compared to the less-dense patch. Thus, intervention implementation that focuses on the patch acquiring “most” time spent or higher transmissibility may have the strongest public health impact required to minimize dengue disease spread. This is consistent with the Global Strategy for Dengue Prevention and Control, which targets the areas of high human-vector contacts (Murray et al., 2013; World Health Organization, 2013). This study allows us to appreciate the challenges in taking action to prevent or decrease the severity of an epidemic of dengue fever. Even when an optimal control strategy is agreed upon, how efficient the implementation of the strategy affects the outcome greatly. These models show the challenges that public health officials have when deciding on strategies to prevent an epidemic. Models can get close to capturing the dynamics of the spread of dengue, but at the same time using simpler models, in our case a two-patch model, allow us to isolate and analyze the main factors that can affect transmission such as coupling parameters, basic reproductive numbers, and final epidemic sizes. Our results can suggest promising directions for future research on dengue transmissions in a multi-patch model in an optimal control framework (Asano et al., 2008; Potts and Kimbrell, 2009). There still much work remain to study further effective control measures between patches with the aid of multiple interdisciplinary approaches. While much more attention has been paid on mathematical models which assume a relatively homogeneous population even though considerable heterogeneity such as space, socio-economics and environments is critical to understand dengue transmission (Gudelj and White, 2004). Such heterogeneities can significantly alter parameter estimates and model predictions therefore, preventive interventions relative to random mixing assumed in simpler models (Rodriguez and Tores-Sorando, 2001; Woolhouse et al., 1997).

11

Acknowledgments This publication was made possible by Grant no. 1R01GM10047101 from the National Institute of General Medical Sciences (NIGMS) at the National Institutes of Health. Its contents are solely the responsibility of the authors and do not necessarily represent the official views of NIGMS.

Appendix A The basic reproductive number is calculated by using the methodology (the next generation matrix approach) outlined in van den Driessche and Watmough (2002). Let x ¼ ðEv1 ; I v1 ; Eh1 ; I h1 ; Ev2 ; I v2 ; Eh2 ; I h2 ÞT . Now, we let F j represent the rate of appearance of new infections in the component j's class, thus F ðxÞ represents all the new infections. The net transition rates out of the corresponding compartment are represented by VðxÞ: 0 1 0 1 κ Ev1 þ μEv1 βv1 ðp11 I h1 þ p21 Ih2 ÞSv1 B μI  κ E C B C B v1 0 v1 C B C B C B C B γ Eh1 C B β h1 ðp11 I v1 þ p12 I v2 ÞSh1 C B C B C B C B C B δI h1  γ Eh1 C 0 B C C F ðxÞ ¼ B C; VðxÞ ¼ B B C B β v2 ðp12 I h1 þ p22 I h2 ÞSv2 C B κ Ev2 þ μEv2 C B C B C B C B μI v2  κ Ev2 C 0 B C B C B C B γ Eh2 C @ β h2 ðp21 I v1 þ p22 I v2 ÞSh2 A @ A δIh2  γ Eh2 0 where we define 8  8 matrices F ¼

h

∂F n ∂xj ðx Þ

i

and V ¼

h

∂V n ∂xj ðx Þ

i

eval-

uated at the disease-free equilibrium, xn , which composes of Svi ¼

Shi ¼ 1 for i¼1, 2 and the rest of them zero. The spectral radius ρ of matrix FV  1 yields the basic reproductive number, where 2 3 0 0 0 0 0 0 0 ðκ þ μÞ  1 6 7 κ 6 μ1 0 0 0 0 0 0 7 6 7 6 ðκ þ μÞμ 7 6 7 1 6 7 γ 0 0 0 0 0 0 0 6 7 6 7 1 1 6 0 0 δ δ 0 0 0 0 7 1 6 7 ¼6 V 1 0 0 0 7 0 0 0 0 ðκ þ μÞ 6 7 6 7 κ 6 7 0 7 0 0 0 0 μ1 0 6 6 7 ðκ þ μÞμ 6 7 1 6 7 γ 0 0 0 0 0 0 0 4 5 1 1 0 0 0 0 0 0 δ δ and 2 0 6 6 6 0 6 6 6 κβh1 p11 6 6 μð κ þ μÞ 6 6 0 1 FV ¼6 6 6 0 6 6 6 6 0 6 6 κβh2 p21 6 4 μð κ þ μÞ 0

0

βv1 p11 γ

βv1 p11 γ

0

0

0

0

0

βh1 p11 μ 0

0

0 0

κβh1 p12 μð κ þ μÞ

0

0

βv1 p21 γ

0

0

βh1 p12 μ 0 0

βv2 p22 γ

βv2 p22 γ

0

0

0

βv2 p12 γ

0

0

0

0

0

κβh2 p22 μð κ þ μÞ

βh2 p22 μ

0

0

0

0

0

0

0

βv2 p12 γ

0

0

0

0

0

βh2 p21 μ

0

3

βv1 p21 γ 7 7

0

0

0

0

0

7 7 7 7 7 7 7 7 7: 7 7 7 7 7 7 7 7 7 5

Thus, the global basic reproductive number for the two patch system is sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffi κ ðϕ1 þ ϕ2 Þ R0 ¼ ð5Þ 2μγ ðκ þ μÞ where

ϕ1 ¼ p211 βh1 βv1 þ p221 βh2 βv1 þ p212 βh1 βv2 þ p222 βv2 βh2 ϕ2 ¼ ðp211 βh1 βv1 þ p221 βh2 βv1 Þ2 þ ðp212 βh1 βv2 þp222 βv2 βh2 Þ2 þ2β v1 β v2 ðp211 p212 β h1 þp221 p222 βh2 Þ 2

2

Please cite this article as: Lee, S., Castillo-Chavez, C., The role of residence times in two-patch dengue transmission dynamics and optimal strategies. J. Theor. Biol. (2015), http://dx.doi.org/10.1016/j.jtbi.2015.03.005i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

S. Lee, C. Castillo-Chavez / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎

12

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

with the transversality conditions λi ðt f Þ ¼ 0 for i ¼ 1; …; 12 and the optimality conditions

þ 2β h1 βh2 β v1 β v2 ð4p11 p12 p21 p22  p211 p222  p212 p221 Þ Next, the optimal control problem for the two-patch model is formulated to minimize the proportions of infected vectors and humans in both patches for a finite time interval at a minimal cost of implementation. We define our objective functional as follows: Z

tf

Jðu1 ðtÞ; u2 ðtÞÞ ¼ 0

8 9 9 8 2 < < X λ8  λ7 λ5  λ4 λ11  λ10 = = pj2 I hj Sv2 þ βh2 p12 I v2 Sh1 þ βh2 p22 I v2 Sh2 ;b un2 ¼ min max a; βv2 : : W4 W4 W4 ; ; j¼1

W 1 ðI h1 ðtÞþ I v1 ðtÞÞ þ W 2 ðI h2 ðtÞþ I v2 ðtÞÞ þ 12 W 3 u21 ðtÞ þ 12 W 4 u22 ðtÞ dt

Then, we seek an optimal pair (Un, Xn) such that JðU n Þ ¼ minfJðUÞj U A Ωg;

ð6Þ

where Ω ¼ fðui ðtÞ A ðL1 ð0; t f ÞÞ2 J a r ui ðtÞ r b; t A ½0; t f ; i ¼ 1; 2g subject to the state equations in (1) with X ¼ ðSv1 ; Ev1 ; I v1 ; Sh1 ; Eh1 ; I h1 ; Rh1 ; Sv2 ; Ev2 ; I v2 ; Sh2 ; Eh2 ; I h2 ; Rh2 Þ and U ¼ ðu1 ; u2 Þ. The existence of optimal controls is guaranteed by standard results of optimal control theory (Fleming and Rishel, 1975). The necessary conditions of optimal solutions are derived from Pontryagin et al.'s (1962) Maximum Principle. This principle converts the system (1), (9), (10) into the problem of minimizing the Hamiltonian H given by H ¼ W 1 ðI h1 ðtÞ þ I v1 ðtÞÞ þ W 2 ðI h2 ðtÞ þ I v2 ðtÞÞ þ 12 W 3 u21 ðtÞ þ 12 W 4 u22 ðtÞ þ λ1 ½μ  βv1 ð1  u1 ðtÞÞðp11 I h1 þ p21 I h2 ÞSv1  μSv1 

þ λ2 ½βv1 ð1  u1 ðtÞÞðp11 I h1 þ p21 I h2 ÞSv1  μEv1  κ Ev1 

Proof. The existence of optimal controls follows from Corollary 4.1 of Fleming and Rishel (1975) since the integrand of J is a convex function of U(t) and the state system satisfies the Lipschitz property with respect to the state variables. The following can be derived from the Pontryagin et al.'s (1962) Maximum Principle: dλ1 ðtÞ ∂H dλ2 ðtÞ ∂H dλ3 ðtÞ ∂H ¼ ¼ ¼ ; ; ; dt ∂Sv1 dt ∂Ev1 dt ∂I v1 dλ4 ðtÞ ∂H dλ5 ðtÞ ∂H dλ6 ðtÞ ∂H ¼ ¼ ¼ ; ; ; dt ∂Sh1 dt ∂Eh1 dt ∂I h1 dλ7 ðtÞ ∂H dλ8 ðtÞ ∂H dλ9 ðtÞ ∂H ¼ ¼ ¼ ; ; ; dt ∂Sv2 dt ∂Ev2 dt ∂I v2 dλ10 ðtÞ ∂H dλ11 ðtÞ ∂H dλ12 ðtÞ ∂H ¼ ¼ ¼ ; ; ; dt ∂Sh2 dt ∂Eh2 dt ∂I h2

þ λ6 ½γ Eh1  δI h1 

with λi ðt f Þ ¼ 0 for i ¼ 1; …; 12 and evaluated at the optimal controls and corresponding states, which results in the adjoint system (8). The Hamiltonian H is minimized with respect to the controls, so we differentiate H with respect to ui on the set Ω, respectively, giving the following optimality conditions:

þ λ8 ½βv2 ð1  u2 ðtÞÞðp12 I h1 þ p22 I h2 ÞSv2  μEv2  κ Ev2 

0¼

þ λ3 ½κ Ev1  μI v1    þ λ4 ½  Sh1 βh1 ð1  u1 ðtÞÞp11 I v1 þ β h2 ð1  u2 ðtÞÞp12 I v2  þ λ5 ½Sh1 ðβh1 ð1  u1 ðtÞÞp11 I v1 þ β h2 ð1 u2 ðtÞÞp12 I v2 Þ  γ Eh1  þ λ7 ½μ  βv2 ð1  u2 ðtÞÞðp12 I h1 þ p22 I h2 ÞSv2  μSv2  þ λ9 ½κ Ev2  μI v2 

þ λ10 ½  Sh2 ðβh1 ð1  u1 ðtÞÞp21 I v1 þ β h2 ð1  u2 ðtÞÞp22 I v2 Þ þ λ12 ½γ Eh2  δI h2 :

2 X ∂H λ4  λ5 ¼ W 3 u1 þ βv1 pj1 I hj Sv1 ðλ1  λ2 Þ þ β h1 p11 I v1 Sh1 ∂u1 W3 j¼1

þ β h1 p21 I v1 Sh2

þ λ11 ½Sh2 ðβ h1 ð1  u1 ðtÞÞp21 I v1 þ βh2 ð1  u2 ðtÞÞp22 I v2 Þ  γ Eh2  ð7Þ

From this Hamiltonian and Pontryagin et al.'s (1962) Maximum Principle, we obtain the following theorem: Theorem 1. There exist optimal controls U n ðtÞ and corresponding state solutions X n ðtÞ that minimize J(U) over Ω. In order for the above statement to be true, it is necessary that there exist continuous functions λi ðtÞ such that 0 1

8 9 9 8 2 < < X λ2  λ1 λ5  λ4 λ11  λ10 = = un1 ¼ min max a; βv1 pj1 I hj Sv1 þ βh1 p11 I v1 Sh1 þ βh1 p21 I v1 Sh2 ;b : : W3 W3 W3 ; ; j¼1

λ ¼ λ1 ½βv1 ð1  u1 ðtÞÞðp11 Ih1 þp21 Ih2 Þ þ μ  λ2 βv1 ð1  u1 ðtÞÞðp11 I h1 þ p21 I h2 Þ 0 λ2 ¼ λ2 ðμ þ κ Þ  λ3 κ λ03 ¼  W 1 þ λ3 μ þ λ4 βh1 ð1  u1 ðtÞÞp11 Sh1  λ5 βh1 ð1  u1 ðtÞÞp11 Sh1 þ λ10 β h1 ð1  u1 ðtÞÞp21 Sh2  λ11 βh1 ð1  u1 ðtÞÞp21 Sh2 λ04 ¼ ðλ4  λ5 Þðβh1 ð1  u1 ðtÞÞp11 I v1 þ βh2 ð1  u2 ðtÞÞp12 Iv2 Þ λ05 ¼ λ5 γ  λ6 γ λ06 ¼  W 1 þ ðλ1  λ2 Þβv1 ð1  u1 ðtÞÞp11 Sv1 þ λ6 δ þ ðλ7  λ8 Þβv2 ð1  u2 ðtÞÞp12 Sv2 λ07 ¼ λ7 ½βv2 ð1  u2 ðtÞÞðp12 Ih1 þp22 Ih2 Þ þ μ  λ8 βv2 ð1  u2 ðtÞÞðp12 I h1 þ p22 I h2 Þ 0 λ8 ¼ λ8 ðμ þ κ Þ  λ9 κ λ09 ¼  W 2 þ λ9 μ þ λ4 βh2 ð1  u2 ðtÞÞp21 Sh1  λ5 βh2 ð1  u2 ðtÞÞp21 Sh1 þ λ10 β h2 ð1  u2 ðtÞÞp22 Sh2  λ11 βh2 ð1  u2 ðtÞÞp22 Sh2 λ010 ¼ ðλ10  λ11 Þðβh1 ð1  u1 ðtÞÞp21 I v1 þ βh2 ð1  u2 ðtÞÞp22 Iv2 Þ λ011 ¼ λ11 γ  λ12 γ λ012 ¼  W 2 þ ðλ1  λ2 Þβv1 ð1  u1 ðtÞÞp21 Sv1 þ λ12 δ þðλ7  λ8 Þβ v2 ð1 u2 ðtÞÞp22 Sv2 ð8Þ

λ10  λ11

W3 2 X ∂H λ4  λ5 0¼ ¼ W 4 u2 þ βv2 pj2 I hj Sv2 ðλ7  λ8 Þ þ β h2 p12 I v2 Sh1 ∂u2 W4 j¼1 þ βh2 p22 I v2 Sh2

λ10  λ11

ð9Þ

W4

Solving for uni ðtÞ we obtain un1 ¼ β v1

2 X

pj1 I hj Sv1

j¼1

þ β h1 p21 I v1 Sh2 un2 ¼ β v2

2 X

þ β h2 p22 I v2 Sh2

W3

þ β h1 p11 I v1 Sh1

λ5  λ4 W3

λ11  λ10 W3

pj2 I hj Sv2

j¼1

λ2  λ1

λ8  λ7 W4

λ11  λ10 W4

þ β h2 p12 I v2 Sh1

λ5  λ4 W4 ð10Þ

By using the standard argument for bounds a r ui ðtÞ r b for i¼ 1,2, we have the optimality conditions. □ References Adams, B., Kapan, D., 2009. Man bites mosquito: understanding the contribution of human movement to vector-borne disease dynamics. Plos One 4 (8), e6763. Agusto, F.B., 2014. Malaria drug resistance: the impact of human movement and spatial heterogeneity. Bull. Math. Biol. 76, 1607–1641. Agusto, F.B., Lenhart, S., 2013. Optimal control of the spread of malaria superinfectivity. J. Biol. Syst. Spec. Issue Infect. Dis. Model. 21 (04). Agusto, Folashade B., Marcus, Nizar, Okosun, Kazeem O., 2012. Application of optimal control to the epidemiology of malaria disease. Electron. J. Differ. Equ. 81, 1–22. Anderson, R.M., May, R.M., 1991. Infectious Diseases of Humans: Dynamics and Control. Oxford University Press, Oxford.

Please cite this article as: Lee, S., Castillo-Chavez, C., The role of residence times in two-patch dengue transmission dynamics and optimal strategies. J. Theor. Biol. (2015), http://dx.doi.org/10.1016/j.jtbi.2015.03.005i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

S. Lee, C. Castillo-Chavez / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎

1 2 3 4 5 6 7 8 9 10 11 12 13 14 Q4 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 Q5 57 58

Arino, J., Ducrot, A., Zongo., P., 2012. A meta-population model for malaria with transmission blocking partial immunity in hosts. J. Math. Biol. 64 (3), 423–448. Arino, J., Jordan, R., van den Driessche, P., 2007. Quarantine in a multi-species epidemic model with spatial dynamics. Math. Biosci. 206, 46–60. Asano, E., Gross, L.J., Lenhart, S., Real, L.A., 2008. Optimal control of vaccine distribution in a rabies metapopulation model. Math. Biosci. Eng. 5, 219–238. Blayneh, K., Cao, Y., Kwon, H.D., 2009. Optimal control of vector-borne diseases: treatment and prevention. Discret. Contin. Dyn. Syst. Ser. B 11, 587–611. Blayneh, K., Gumel, A.B., Lenhart, S., Clayton, T., 2010. Backward bifurcation and optimal control in transmission dynamics. Bull. Math. Biol. 72, 1006–1028. Brauer, F., Castillo-Chavez, C., 2012. Mathematical Models in Population Biology and Epidemiology. Springer Verlag, New York. Burattini, M.N., Chen, M., Chow, A., Coutinho, F.A.B., Goh, K.T., Lopez, L.F., Ma, S., Massad, E., 2008. Modelling the control strategies against dengue in Singapore. Epidemiol. Infect. 136 (3), 309–319. Carlos, T.A., 2009. Deterministic and Stochastic Metapopulation Models for Dengue Fever (thesis), ASU. Center for Disease Control, Dengue. 〈http://www.cdc.gov/NCIDOD/DVBID/DENGUE〉. Chowell, G., Torre, C.A., Munayco-Escate, C., Suarez-Ognio, L., Lopez-Cruz, R., Hyman, J.M., Castillo-Chavez, C., 2008. Spatial and temporal dynamics of dengue fever in Peru: 1994–2006. Epidemiol. Infect. 8, 1–11. Chowell, G., Diaz-Duenas, P., Miller, J.C., 2007. Estimation of the reproduction number of dengue fever from spatial epidemic data. Math. Biosci. 208, 571–589. Cosnera, C., Beierb, J.C., Cantrella, R.S., Impoinvilc, D., Kapitanskia, L., Pottsd, M.D., Troyoe, A., Ruan, S., 2009. The effects of human movement on the persistence of vector-borne diseases. J. Theor. Biol. 258, 550–560. Carrasquilla, G., Quintero, J., Surez, R., Gonzlez, C., 2003. Dengue in Colombia, National Institute of Health (Dengue Colombia, Este documento fue elaborado en junio del 2003, con el apoyo del Centro Internacional de Investigaciones para el Desarrollo (IDRC, por sus siglas en ingls) del Canad y del Instituto Salud y Trabajo del Per. Castillo-Chavez, C., Song, B., Zhang, J., 2003. An Epidemic Model with Virtual Mass Transportation: The Case of Smallpox in a Large City In: SIAM's Frontiers in Applied Mathematics Series. Castillo-Chavez, C., Lee, S., 2013. Epidemiology modeling. In: Engquist, B., et al. (Eds.), Encyclopedia of Applied and Computational Mathematics. SpringerVerlag, Berlin and Heidelberg. Diekmann, O., Heesterbeek, J.A.P., Mathematical Epidemiology of Infectious Diseases: Model Building, Analysis and Interpretation. John Wiley & Sons, Ltd., New York. Endy, T.P., Nisalak, A., Chunsuttitwat, S., Vaughn, D.W., Green, S., Ennis, F.A., Rothman, A.L., Libraty, D.H., 2004. Relationship of preexisting dengue virus (DV) neutralizing antibody levels to viremia and severity of disease in a prospective cohort study of DV infection in Thailand. J. Infect. Dis. 189, 990–1000. Favier, C., Schmit, D., Muller-Graf, C.D.M., Cazelles, B., Degallier, N., Mondet, B., Dubois, M.A., 2005. Influence of spatial heterogeneity on an emerging infectious disease: the case of dengue epidemics. Proc. R. Soc. B 272, 1171–1177. Favier, C., Degallier, N., Rosa-Freitas, M.G., Boulanger, J.P., Costa Lima, J.R., LuitgardsMoura, J.F., Menkes, C.E., Mondet, B., Oliveira, C., Weimann, E.T.S., Tsouris, P., 2006. Early determination of the reproductive number for vector-borne diseases: the case of dengue in Brazil. Trop. Med. Health 11, 332–340. Fleming, W.H., Rishel, R.W., 1975. Deterministic and Stochastic Optimal Control. Springer Verlag, New York. Goh, K.T., 1998. Dengue—a re-emerging infectious disease in Singapore. Dengue Singap., 33–49. Gubler, D.J., Kuno, G., 1997. Dengue and Dengue Hemorrhagic Fever. CABI Publishing, Massachusetts. Gubler, D.J., 1998. Resurgent vector-borne diseases as a global health problem. Emerg. Infect. Dis. 4 (3), 442–450. Gubler, D.J., 1997. Human behaviour and cultural context in disease control. Trop. Med. Int. Health. 2 (11), A1–2. Gudelj, I., White, K.A.J., 2004. Spatial heterogeneity, social structure and disease dynamics of animal populations. Theor. Popul. Biol. 66, 139–149. Guy, B., et al., 2011. From research to phase III: preclinical, industrial and clinical development of the Sanofi Pasteur tetravalent dengue vaccineVaccine 29, 7229–7241. Guirakhoo, F., et al., 2000. From recombinant chimeric yellow fever-dengue type 2 virus is immunogenic and protective in nonhuman primates. J. Virol., 5477–5485. Guzman, M.G., Kouri, G., 2003. Dengue and dengue hemorrhagic fever in the Americas: lessons and challenges. J. Clin. Virol. 27, 1–13. Halstead, S.B., 2004. Dengue fever/dengue hemorrhagic fever. In: Cohen, J., Powderly, W.G. (Eds.), Infectious Diseases. Mosby, pp. 1681–1689. Halstead, S.B., Deen, J., 2002. The future of dengue vaccines. Lancet 360, 1243–1245.

13

Harrington, L.C., Scott, T.W., Lerdthusnee, K., Coleman, R.C., Costero, A., Clark, G.G., Jones, J.J., Kitthawee, S., Kittayapong, P., Sithiprasasna, R., Edman, J.D., 2005. Dispersal of the dengue vector Aedes aegypti within and between rural communities. Am. J. Trop. Med. Hyg. 72 (2), 209–220. Hayes, C.G., Phillips, I.A., Callanhan, J.D., Griebenow, W.F., Hyams, K.C., Wu, S.J., Waits, D.M., 1996. The epidemiology of dengue virus infection among urban, jungle, and rural populations in the Amazon region of Peru. Am. J. Trop. Med. Hyg. 55 (4), 459–463. Herrera-Valdez, M.A., Cruz-Aponte, M., Castillo-Chavez, C., 2011. Multiple outbreaks for the same pandemic: Local transportation and social distancing explain the different “waves” of A-H1N1pdm cases observed in Mexico during 2009. Math. Biosci. Eng. 8 (1), 21–48. Homes, E.C., Twiddy, S.S., 2003. The origin, emergence and evolutionary genetics of dengue virus. Infect. Genet. Evol. 1, 19–28. Lenhart, S., Workman, J., 2007. Optimal Control Applied to Biological Models. Chapman and Hall/CRC Press, New York. Lee, S., Chowell, G., Castillo-Chavez, C., 2010. Optimal control of influenza pandemics: the role of antiviral treatment and isolation. J. Theor. Biol. 265, 136–150. Lee, S., Morales, R., Castillo-Chavez, C., 2011. A note on the use of influenza vaccination strategies when supply is limited. Math. Biosci. Eng. Spec. Issue H1N1 2009 Influenza Models 8.1, 172–182. Lee, S., Golinski, M., Chowell, G., 2012. Modeling optimal age-specific vaccination strategies against pandemic influenza. Bull. Math. Biol. 74, 958–980. MacDonald, G., 1957. The Epidemiology of Control of Malaria. Oxford University Press, London. Martens, P., Hall, L., 2000. Malaria on the move: human population movement and malaria transmission. Emer. Infect. Dis. 6, 103–109. Morrison, A., Gray, K., Getis, A., Astete, H., Sihuincha, M., Focks, D., Watts, D., Stancil, J.D., Olson, J.G., Blair, P., Scott, T.W., 2004. Temporal and geographic patterns of Aedes aegypti (Diptera: Culicidae) production in Iquitos, Peru. J. Med. Entomol. 41, 1123–1142. Murray, J.D., 2003. Mathematical Biology. II: Spatial Models and Biomedical Applications. Springer. Murray, N.A., Quam, M.B., Wilder-Smith, A., 2013. Epidemiology of dengue: past, present and future prospects. Clin. Epidemiol. 5, 299–309. Nishiura, H., 2006. Mathematical and statistical analyses of the spread of dengue. Dengue Bull., 30–51. Ooi, E.E., Goh, K.T., Gublert, D.J., 2006. Dengue prevention and 35 years of vector control in Singapore. Emerg. Infect. Dis. 12, 887–893. Phillips, I., et al., 1992. First documented outbreak of Dengue in the Peruvian Amazon Region. Bull. PAHO 26 (3), 38–42. Pontryagin, L.S., Boltyanskii, V.G., Gamkrelidze, R.V., Mishchenko, E.F., 1962. The Mathematical Theory of Optimal Processes. Wiley, New Jersey. Potts, M.D., Kimbrell, T., 2009. In: Cosner, Chris, Cantrell, Stephen, Ruan, Shigui (Eds.), Spatial Ecology. Chapman and Hall, CRC, New York, pp. 273–291. Reiskind, M.H., Baisley, K.J., Calampa, C., Sharp, T.W., Watts, D.M., Wilson, M.L., 2001. Epidemiological and ecological characteristics of past dengue virus Infection in Santa Clara, PeruTrop. Med. Int. Health 6, 212–218. Restrepo, A.C., Baker, P., Clements, A.C.A., 2014. National spatial and temporal patterns of notified dengue cases, Colombia 2007–2010. Trop. Med. Int. Health 7, 863–871. Ross, R., 1910. The Prevention of Malaria. John Murray, London. Rodriguez, D.J., Tores-Sorando, L., 2001. Models of infectious diseases in spatially heterogeneous environments. Bull. Math. Biol. 63, 547–571. Rowthorn, R.E., Laxminarayan, R., Gilligan, C.A., 2008. Optimal control of epidemics in metapopulations. J. R. Soc. Interface 6, 1135–1144. Sparagano, O.A., de Luna, C.J., 2008. From population structure to geneticallyengineered vectors: new ways to control vector-borne diseases? Infect. Genet. Evol. 8 (4), 520–525. Sattenspiel, L., 2009. The Geographic Spread of Infectious Disease. Princeton University Press, New Jersey. Siqueira, J.B., Martelli, C.M.T., Coelho, G.E., da Rocha Simplicio, A.C., Hatch, D.L., 2005. Dengue and dengue hemorrhagic fever, Brazil, 1981–2002. Emerg. Infect. Dis. 11, 48–53. van den Driessche, P., Watmough, J., 2002. Reproduction numbers and subthreshold endemic equilibria for compartmental models of disease transmission. Math. Biosci. 180, 29–48. World Health Organization. Dengue Control Strategies [webpage on the Internet, cited March 2, 2013]. Woolhouse, M.E., Dye, C., Etard, J.F., Smith, T., Charlwood, J.D., Garnett, G.P., Hagan, P., Hii, J.L., Ndhlovu, P.D., Quinnell, R.J., Watts, C.H., Chandiwana, S.K., Anderson, R.M., 1997. Heterogeneities in the transmission of infectious agents: implications for the design of control programs. Proc. Natl. Acad. Sci. USA 94, 338–342.

Please cite this article as: Lee, S., Castillo-Chavez, C., The role of residence times in two-patch dengue transmission dynamics and optimal strategies. J. Theor. Biol. (2015), http://dx.doi.org/10.1016/j.jtbi.2015.03.005i

59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115