A Two-Patch Mathematical Model for Temperature- Dependent Dengue Transmission Dynamics

Dengue fever has been a threat to public health not only in tropical regions but non-tropical regions due to recent climate change. Motivated by a recent dengue outbreak in Japan, we develop a two-patch model for dengue transmission associated with temperature-dependent parameters. The two patches represent a park area where mosquitoes prevail and a residential area where people live. Based on climate change scenarios, we investigate the dengue transmission dynamics between the patches. We employ an optimal control method to implement proper control measures in the two-patch model. We find that blockage between two patches for a short-term period is effective in a certain degree for the disease control, but to obtain a significant control effect of the disease, a long-term blockage should be implemented. Moreover, the control strategies such as vector control and transmission control are very effective, if they are implemented right before the summer outbreak. We also investigate the cost-effectiveness of control strategies such as vaccination, vector control and virus transmission control. We find that vector control and virus transmission control are more cost-effective than vaccination in case of Korea.


Introduction
Dengue fever is a vector-borne disease spread by Aedes type mosquitoes such as Aedes aegypti and Aedes albopictus. Since Aedes mosquitoes were generally found in tropical regions, dengue fever has been known as a tropical disease [1]. However, recent dengue outbreaks are expanding beyond the tropic regions by climate change due to global warming [2]. It has been reported that dengue transmission is affected by the climate environment [3][4][5] and in particular, the temperature strongly affects the dengue dynamics [6,7].
Recently, 160 cases of confirmed autochthonous dengue fever were reported in Tokyo, Japan, and most of the confirmed cases have been exposed to mosquito bites at Yoyogi Park in the city [8,9]. In case of Korea, a neighboring country of Japan, although there is no autochthonous dengue case yet, the dengue fever has been predicted to be one of the most probable major threats to public health in the near future [10], and it has been shown that the frequency of the imported dengue cases in Korea and Japan has a similar pattern [11]. Moreover, the number of the imported dengue cases have been increasing recently in Korea [12]. In particular, Seoul, the most populated capital city of Korea, has several big parks where mosquitoes prevail like Tokyo, and the city would be at a risk from dengue transmission in the future [13]. 2 of 26 In this paper, we develop a mathematical model associated with temperature-dependent parameters for describing dengue transmission between two patches which represent a park area where the dengue vector inhabits and an urban area where humans reside. Based on the Representative Concentration Pathway (RCP) climate change scenarios, we investigate the effect of control strategies for the dengue transmission in the two-patch model using the optimal control method and cost-effectiveness analysis.

Two-Patch Dengue Transmission Model
In this section, we develop a two-patch dengue transmission model by applying differential equation approach. It is assumed that patch 1 is a park area where mosquitoes prevail, and patch 2 is a residential area where people live. The focus area for the model is Seoul Forest Park (patch 1) and the residential area (patch 2) around the park in Seoul, Korea. A schematic diagram of the full two-patch model is shown in Figure 1. The model considers the states of mosquito larvae (susceptible (S ei ) and infectious (I ei ) by vertical infection), female adult mosquitoes (susceptible (S vi ), exposed (E vi ) and infectious (I vi )) and humans (susceptible (S hi ), exposed (E hi ), infectious (I hi ) and recovered (R hi )), for patch i = 1, 2. We denote the total larvae population, female adult mosquito population and human population by N ei , N vi and N hi for patch i = 1, 2. That is, N ei = S ei + I ei , N vi = S vi + E vi + I vi and N hi = S hi + E hi + I hi + R hi . To describe the transmission dynamics in patch 2, we use the dengue model in [14]. In our two-patch model, we assume that humans can move between patches, but mosquitoes cannot. We write the governing equations of the model as follows: Host (1)

Patch 2
Vectoṙ In the governing Equation (1), the parameters relevant to larvae and mosquitoes are described as follows: ω is the maturation rate of pre-adult mosquitoes, and µ v and µ e are the mortality rate of adult mosquitoes and larvae, respectively. ν and 1/ε denote the rate of vertical infection from infected mosquitoes to eggs and the extrinsic incubation period, respectively, and δ i is the number of new recruits in the larva stage for patch i = 1, 2. The parameters β vh = bb h and β hv = bb v are the transmissible rates from mosquito to human and from human to mosquito, respectively, where b is the daily biting rate of a mosquito and b v and b h are the probability of infection from human to mosquito per bite and the probability of infection from mosquito to human per bite, respectively [14]. The parameters µ hb and µ hd represent the human birth rate and death rate, respectively, and the two rates are assumed to be equal. The parameters 1/α and 1/γ are the latent period and infectious period for humans, respectively. The inflow rate of infection due to international travelers is defined by η [14]. p ij refers to the human movement rate from patch i to j, where ∑ 2 j=1 p ij = 1 and 0 ≤ p ij ≤ 1 for i = 1, 2. Since there are about 7,500,000 visitors to Seoul Forest Park each year [15], approximately 20,550 people visit the park daily. Hence, assuming N h1 (0) = 20, 000, N h2 (0) = 480, 000, i.e., the total human population of both patches is 500, 000, we compute p 21 = 20, 550/500, 000 = 0.0411. Moreover, we assume that p 11 = 0.001, which represents that a small number of people such as park keepers and homeless people stay in the park, and p 12 = 1 − p 11 = 0.999.
The parameters in the system (1) are described with their values in Table 1.

Parameter Estimation
The temperature-sensitive parameters for the dengue mosquitoes have been studied in previous researches [14,16,[21][22][23][24][25]. Using the previous results, we describe the parameters sensitive to the temperature as the following temperature-dependent functions.
(8) The number of new recruits in the larvae stage δ for patch i = 1, 2 is computed as [14] δ i = µ v N vi + µ e N ei . Figure 2 illustrates the the graph of the temperature-dependent parameters. For the temperature data, we utilize RCP scenarios, which provide four representative scenarios such as the low level scenario (RCP 2.6), the two medium level scenarios (RCP 4.5/6.0) and the high level scenario (RCP 8.5) [27]. Since patch 1, Seoul Forest Park, is located in Seongdong-gu, Seoul, the RCP temperature data for Seongdong-gu is used in the simulation. One can see the tendency of the temperature rise between 2030 and 2100 according to the RCP scenarios (refer to Appendix A).

The Seasonal Reproduction Number
The basic reproduction number is important in epidemiology since it measures the expected number of infectious cases directly caused by one infectious case in a susceptible population. It is known that if R 0 < 1, the system has a locally asymptotically stable disease-free equilibrium, but if R 0 > 1, it has an unstable disease-free equilibrium [28]. However, when some model parameters are time-dependent as in Table 1, one has to use the seasonal reproduction number R s instead of the basic reproduction number [29]. Theorem 1. The seasonal reproduction number R s corresponding to a single patch model with only patch 2 without the inflow rate η is computed as Proof. The proof of the theorem can be found in Appendix B.
Theorem 2. The seasonal reproduction number R s in the two-patch model (1) is the spectral radius ρ of the matrix G, i.e., Proof. The proof of the theorem can be found in Appendix B. Figure 3 shows the seasonal reproduction number R s for three years from 1 January 2030 for each RCP scenario. It is observed that the value of R s is much higher than 1 in the summer season for all RCP scenarios. This implies that it is very likely that the dengue outbreak will occur during the summer.

Dengue Transmission Dynamics Based on Rcp Scenarios
In this section, we perform numerical simulations for the two patch model based on RCP scenarios. For each simulation, we assume the initial condition as N h1 (0) = 20, 000, N h2 (0) = 480, 000, N v1 (0) = 0.5 × N h (0), N v2 (0) = N h (0), where N h (0) = N h1 (0) + N h2 (0). Moreover, we assume that there are no infected mosquitoes and humans initially, that is, I ei (0) = E vi (0) = I vi (0) = 0 and E hi (0) = I hi (0) = 0 for i = 1, 2, so that the first infection in the model is initiated by the inflow of infected international travelers. Figure 4 shows the time evolution of human incidence and cumulative human incidence from January 1 in 2030 for 10 years in the model without control for each RCP scenario. One can see that there will be more dengue incidences for RCP 2.6 and 8.5 than RCP 4.5 and 6.0, which implies that the two extreme level RCP scenarios might provide more favorable temperature environment for dengue virus transmission than the two medium level RCP scenarios.

The Effects of Human and Vector Controls
In the case of the dengue outbreak in Tokyo, Yoyogi park in the city was known as an infection hub and the closure of the park turned out to be very effective control strategies [8]. In accordance with the case in Tokyo, it is worth investigating the effects of the control strategies including the park closure as well as vector and human controls for our model (1). We assume that the park closure begins when the cumulative incidence is over 10. In order to see the control effect of the park closure, we set p 12 = p 21 = 0, and all human population stays in the city area (patch 2). Figure 5 shows the effect of park closure for duration 3, 5, 10, 30, 60 days. It is observed that the park closure for a short-term period such as 3 and 5 days would be effective in a certain degree and the closure for a long-term period such as 30 and 60 days would make a significant control effect for all RCP scenarios. Human Incidence Concerning the control for mosquitoes and humans, it was estimated that the current level of control for mosquitoes in Korea is about 2% increase of mosquito death rate and 2% decrease of transmissible rate between mosquitoes and humans, respectively, and these control measures were implemented between May and October in each year [29]. Now we compare the control effects of the vector death rate, dengue transmissible rate and park closure. In Figure 6, it is assumed that the controls of the vector death rate and dengue transmissible rate are implemented as a 2% increase and a 2% decrease of the rates, respectively, and the park closure is made for only 30 days in the year 2030 at the early stage of the dengue outbreak. Figure 6 shows that the vector control is more effective than the transmission control, and the combination of the park closure and vector and human control is most effective. Now we compare the control effects of the vector death rate, dengue transmissible rate and park closure. In Figure 6, it is assumed that the controls of the vector death rate and dengue transmissible rate are implemented as a 2% increase and a 2% decrease of the rates, respectively, and the park closure is made for only 30 days in the year 2030 at the early stage of the dengue outbreak. Figure 6 shows that the vector control is more effective than the transmission control, and the combination of the park closure and vector and human control is most effective.

Optimal Control
In this section, we implement effective control measures by formulating an optimal control problem for the two-patch dengue model. By incorporating the control functions (1 − u 1 ) and (1 + u 2 ) into the transmissible rate between the vector and human and the mortality rate of the vector in each patch, respectively, in the model Equation (1), we obtain the controlled two-patch system (2) as follows:

Optimal Control
In this section, we implement effective control measures by formulating an optimal control problem for the two-patch dengue model. By incorporating the control functions (1 − u 1 ) and (1 + u 2 ) into the transmissible rate between the vector and human and the mortality rate of the vector in each patch, respectively, in the model Equation (1), we obtain the controlled two-patch system (2) as follows: Now we set up an optimal control problem for the two-patch model in order to minimize the proportions of infected vectors and humans in both patches for a finite time interval at a minimal cost of implementation. We first define the classical objective functional [30,31] where W 1 and W 2 denote the weight constants on the infected humans and vectors and the total vectors, respectively, and W 3 and W 4 denote the weight constants that are the relative costs of the implementation of the preventive controls for decreasing the transmissible rate between vector and human and increasing the mortality rate of the vector, respectively. Then, we find an optimal solution (U * , X * ) that satisfies where It is known that the standard results of optimal control theory guarantees the existence of optimal controls, and the necessary conditions of optimal solutions can be derived from Pontryagin maximum principle [31,32]. The Pontryagin maximum principle converts the system (2) into the problem of minimizing the Hamiltonian H given by Using the Hamiltonian H and the Pontryagin maximum principle, we obtain the theorem. Theorem 3. There exist optimal controls U * (t) and state solutions X * (t) which minimize J(U) over Ω in (3). In order for the above statement to be true, it is necessary that there exist continuous functions λ j (t) such thaṫ with the transversality conditions λ j (t f ) = 0 for j = 1, ..., 16 and the optimality conditions Proof. The proof of the theorem can be found in Appendix B.
We assume the control duration as 5 years throughout the simulations, and the upper bound for u i , i = 1, 2 is 0.1, since the control resources are limited. In Figures 7-9, we simulate the effects on incidence and optimal control functions from different control strategies when only the transmissible rates β hv , β vh are controlled, only the mortality rate µ v is controlled and both the transmissible rates β hv , β vh and the mortality rate µ v are controlled, respectively. Considering the ratio of the number of infected vectors and humans to the total number of vectors, we use the weight constants W 1 = 1, W 2 = 0.0001, W 3 = 1000 and W 4 = 2000 for Figures 7-9. Figure 7 suggests that when the transmission control is considered, it is effective to focus on the control during the summer. Moreover, Figure 8 shows that the peaks of u 2 occasionally occurred, and Figure 9 implies that when all β hv , β vh and µ v are controlled, the control period of µ v is longer than β hv and β vh . These results imply that it is more effective if the controls are concentrated right before the summer outbreak, since seasonal patterns are observed in all cases. For more cases with different weight constants, refer to Appendix C.

Vaccination Model and Cost-Effectiveness of Control Strategies
Recently the vaccine for dengue such as Dengvaxia has been used to prevent dengue transmission in the dengue-endemic countries including Mexico, Philippines, Indonesia and Brazil [33]. A recent research investigated the cost-effectiveness of dengue vaccination in Mexico [20], which concluded that a proper dengue vaccination program would be very cost-effective and also highly reduce the dengue cases and casualties. Thus, it is worth investigating the cost-effectiveness of control strategies including vaccination in our two-patch model. We assume that the vaccination will begin at the year 2030 and any susceptible individual, either seropositive or seronegative, can be vaccinated.
In order to consider the effect of the vaccination, we first construct a two-patch dengue transmission model with vaccination by modifying the model (1). The model with vaccination includes the new compartments such as S V hi , E V hi , I V hi and R V hi which denote the vaccinated susceptible, exposed, infected and recovered human class in patch i = 1, 2, respectively. The schematic diagram for the model is shown in Figure 10. The governing equation for the model is written as follows: where the parameter ψ denotes the rate at which vaccine wanes off and the vaccination rate φ followed by antibody formation is computed by , where a is the proportion of the vaccinated humans and b is the vaccination period. Here, we assume a = 0.3 and b = 120 days between June and September. The parameters relevant to vaccination are described in Table 2. We evaluate the cost-effectiveness of the control measures such as vector control, dengue virus transmission control, and vaccination by using the incremental cost-effectiveness ratio (ICER) in terms of dollars per quality-adjusted life year (QALY) gained for a range of control costs. The cost-effectiveness of a control strategy is related to the costs per disability-adjusted life years (DALY) and gross domestic product (GDP) per capita; (i) a control strategy is very cost-effective if the costs per DALY are less than GDP per capita, (ii) cost-effective if the costs per DALY are between GDP per capita and 3× GDP per capita and (iii) not cost-effective if the costs per DALY are greater than 3× GDP per capita [20,36]. The QALY function Q is computed as follows [20,37]: where D is the disability weight for dengue fever(DF), dengue hemorrhagic fever(DHF) and death, which are denoted by D DF , D DHF and D Death , respectively, C is the age-weighting correction constant, h is the parameter from the age-weighting function, a is the average age, L is the duration of the disability or the years of life lost due to premature death expressed in years such as L DF , L DHF or L Death and r is the social discount rate [37]. The total QALYs lost (TQL) is computed as follows [20]: where the rates of new DF, DHF and death cases for the vector and transmission control are and the rates for the vaccination case are The total cost is obtained by the sum of the cost for direct control along the control strategies and the cost associated with dengue infection under the observation period T f . The direct controls we consider are vector control, transmission control and vaccination where the costs of each direct control are C µ v , C β and C V , respectively. The costs associated with dengue infection DF and DHF are C DF and C DHF , respectively, where each cost is estimated from [20]. The total cost (TC) for each control strategy is computed as follows: (i) Total costs for vector control: Total costs for vaccination: The parameters relevant to cost-effectiveness are described in Table 3. In order to evaluate the cost-effectiveness of controls, we use the incremental cost-effectiveness ratio (ICER) which is calculated by dividing the difference in total costs (incremental cost) by the difference in the QALY with and without control cases. Figure 11 shows ICER per DALY averted for the vaccination, vector control and transmission control cases. Here we assume, as in Section 3.2, that the vector control and transmission control are implemented as a 2% increase of mosquito death rate and a 2% decrease of transmissible rate between May and October, respectively. According to recent researches, it was estimated that the cost of vaccination per capita is hundreds of USD [20,43], while the cost per capita of the vector control and the transmission control ranges from tens of cents to a few dollars in USD [44,45]. Thus, the results of Figure 11 imply that vector control and transmission control are relatively more cost-effective than vaccination, and vaccination is not a suitable control strategy in Korea in terms of cost-effectiveness. Transmission control (C =0.98) Figure 11. Incremental cost-effectiveness ratio (ICER) in log 10 scale for the vaccinate, vector control and transmission control cases.

Discussion and Conclusions
In this paper, we developed the two-patch dengue transmission model associated with temperature-dependent parameters. The focus area for the model was Seoul Forest Park (patch 1) and the residential area (patch 2) around the park in Seoul, the most populated city in Korea. In the model, we represented the parameters sensitive to temperature as the temperature-dependent nonlinear functions using the previous literatures. Using the temperature data under RCP climate change scenarios, we investigated the dengue transmission dynamics within and between patches. The simulation results for the model showed that if a dengue infection is initiated by the inflow of infected international travelers into the focus area, there will be thousands of infected humans within 10 years in the case of no controls for the dengue disease.
We derived the formulas for the seasonal reproduction number R s for the single-patch (patch 2) and the two-patch model by using the next generation matrix. The simulation results ( Figure 3) for R s showed that the value of R s is much bigger than 1 in the summer season for all RCP scenarios. This implies that it is very likely that the dengue outbreak will occur during the summer in the near future, if there is no proper control strategy. To reduce the potential of the dengue outbreak, proper control strategies should be implemented.
We studied optimal control strategies by using an optimal control framework under different scenarios. We found that the control strategies are effective if they are implemented right before the summer outbreak. Concerning the park closure, we found that the closure for a short-term period such as 3 and 5 days would be effective in a certain degree, but the closure for a long-term period such as 30 and 60 days would make a substantial control effect.
By incorporating the vaccination policy into the two-patch model, we constructed the two-patch dengue transmission model (5) with vaccination. Concerning the vaccination, currently Dengvaxia is vaccinated for seropositive cases and also other vaccines are now in clinical development [46,47]. Since our model simulation begins at the year 2030, there is a possibility of the development of other vaccines which can be used for any susceptible cases, either seropositive or seronegative. In light of this aspect, in the vaccination model (5), we assume that any susceptible individual, seropositive or seronegative, can be vaccinated. We investigated the cost-effectiveness of the control policies such as vaccination, vector control, and transmission control, using the incremental cost-effectiveness ratio (ICER) in terms of dollars per quality-adjusted life year (QALY). We found that the transmissible rate control and the vector control are cost-effective while the vaccination is less cost-effective. This result is not compatible with the result in Mexico [20] because there are not a sufficient number of infective humans in Korea, compared to the case of Mexico.
Since there have been no autochthonous dengue cases in Korea yet, the parameters relevant to the cost of the dengue vaccination were adapted from the previous studies in dengue-endemic regions [20,39,40]. Moreover, to the best of our knowledge, there have been no studies about the relation between the cost and effectiveness of the dengue control in Korea. Although these factors may result in a limitation in the accurate estimation of the cost-effectiveness for the control strategies, the simulation results ( Figure 11) clearly show the difference in cost-effectiveness between different strategies; when the control resources are limited, it is more effective to implement vector control and transmission control rather than vaccination.
In this research, we used the regional data such as temperature and human movement rate for Seoul, Korea. Thus, most of the results presented in this paper may not be applied directly to the area with different environment for mosquitoes and humans, but we expect that the modeling approach presented in this work will be applied to other cases, especially when the temperature-dependent transmission dynamics between endemic and non-endemic regions are investigated.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A. Temperature Data under RCP Scenarios
The Intergovernmental Panel on Climate Change has developed RCP scenarios in 2014, and four representative scenarios are the lowest-level scenario (RCP 2.6), the two medium level scenarios (RCP 4.5/6.0) and the high-level scenario (RCP 8.5) [27]. In this paper, we use the daily climate data estimated by the Korea Meteorological Administration under the four RCP scenarios. Figure A1 illustrates the 5-year averages of daily temperature and the ranges from the mean temperature in the summer (June to August) to the mean temperature in the winter (December to February) for five years in Seongdong-gu, Seoul from year 2030 to 2099.
Here F (x) denotes all of the new infections and V (x) denotes the net transition rates of the corresponding compartment. F and V are 5 × 5 matrices at x 0 given by Hence, the next generation matrix G is computed as Since the seasonal reproduction number R s for the single-patch model is the dominant eigenvalue of the matrix G, R s is obtained as Proof of Theorem 2. The system (1) has the disease-free state x 0 = (S ei , 0, S vi , 0, 0, S hi , 0, 0, 0) with η = 0. Let x = (I ei , E vi , I vi , E hi , I hi ) T for i = 1, 2. If F (x) and V (x) denote the functions for all of the new infections and the net transition rates of the corresponding compartment, respectively, then one can obtain Thus, F and V are 10 × 10 matrices at x 0 given by Hence, one can obtain the next generation matrix G as = αS v2 β hv (p 12 (α + γ + p 21 + p 12 ) + αγ − p 21 µ hd ) N h2 (α + p 12 + p 21 )(γ + p 12 + (1 − g)p 21 )(α + µ hd )(γ + µ hd ) G 4,9 = (p 12 − µ hd )S v2 β hv N h2 (γ + p 12 + (1 − g)p 21 )(γ + µ hd ) G 4,10 = (γ + p 12 )S v2 β hv N h2 (γ + p 12 + (1 − g)p 21 )(γ + µ hd ) Finally, the seasonal reproduction number R s for the two-patch model is computed as the spectral radius ρ of the next generation matrix G, i.e., R s = ρ(G).
Solving for u i (t), one can obtain By using the standard argument for bounds a ≤ u i ≤ b for i = 1, 2, we have the optimality conditions.

Appendix C. Optimal Control Result with Different Weight Constants
We consider various weight constants for the case that all of β hv , β hv and µ v are controlled. Figure A2a-c show the plots of u 1 and u 2 for W 4 = 1000, 3000, 5000, when W 1 = 1, W 2 = 0.0001, W 3 = 1000 are fixed. One can see the similar results as in Section 3.3.