Optimal Strategies for Dengue Prevention and Control during Daily Commuting between Two Residential Areas

: In this paper, we report an application for the mathematical theory of dynamic optimization for design of optimal strategies that account for daily commuting of human residents, aiming to reduce vector-borne infections (dengue) among human populations. Our analysis is based on a two-patch dengue transmission model amended with control variables that represent personal protection measures aimed at reduction of the number of contacts between mosquitoes and human hosts (e.g., the use of repellents, mosquito nets, or insecticide-treated clothing). As a result, we have proposed and numerically solved an optimal control problem to minimize the costs associated with the application of control measures, while also minimizing the total number of dengue-infected people in both residential areas. Our principal goal was to identify an optimal strategy for personal protection that renders the maximal number of averted human infections per unit of invested cost, and this goal has been accomplished on the grounds of cost-e ﬀ ectiveness analysis.


Introduction
Population mobility is one of the factors that have historically influenced the spread of epidemics. An infection that affects individuals in some geographically isolated area (patch or zone) can reach other locations due to people traveling. This is true for many infectious diseases, including those that are not transmitted by direct people-to-people contacts, but rather by means of an infectious agent (transmitting vector), such as dengue, Chikungunya, Zika, etc.
In the present study, we focus on the transmission of dengue infections that are spread by the female mosquito, mainly of the species Aedes aegypti, while accounting for population mobility. This mosquito species is closely associated with human habitation due to its blood-feeding habits and the presence of breeding sites widely available around households (desert coolers, flower vases, potted plants, water tanks, cisterns, and other stored water). Adult mosquitoes usually spend all their life in a radius not exceeding 100 meters around the breeding site they have emerged from, provided there is food and other resources for their subsequent reproduction.
Aedes aegypti females usually bite people during daylight, and they rest indoors after darkness. Therefore, the use of bed nets does little to protect people from their bites. Aedes aegypti females need to ingest human blood to mature their eggs and may acquire the dengue virus (DENV) during the blood-feeding on a DENV-infected person. When taking a subsequent blood meal on another (uninfected) person, an infected female mosquito injects her saliva to prevent the host's blood from clotting and to ease feeding. This injection of saliva infects the host with the dengue virus. Infected mosquitoes continue to transmit dengue with each blood meal for the rest of their lives, while infected human hosts usually remain infectious for about 5-12 days, and acquire life-long immunity for homologous DENV strains and temporal cross-immunity for heterologous DENV strains. For more information regarding the mosquito biology and ecology, as well as epidemiology of dengue, the reader is invited to revise the book [1].
In the absence of an effective vaccine [2], dengue morbidity among human hosts can only be reduced by appropriate vector control measures. Among these measures, we can detect two main approaches: • Suppression of the mosquito population using chemical substances (e.g., insecticides or larvicides), biological control agents (natural predators of mosquito larvae and pupae), environmental management (e.g., lethal ovitraps or elimination of mosquito breeding sites in and around households), and the release of genetically modified mosquitoes (sterile males or insects carrying a dominant lethal gene).

•
Reduction of effective contact between female mosquitoes and humans (or mosquito bites) through the use of repellents, mosquito nets, insecticide-treated clothing, and other measures targeting personal protection.
In this work, we address the second approach and illustrate it using an example of our home city (Cali, Colombia) together with its suburbs, since this area is considered hyper-endemic regarding dengue morbidity among human residents [3].
To formalize our study, we present a two-patch dengue transmission model, where human hosts residing in both patches may commute between them while mosquitoes do not relocate from their zones. This metapopulation model describes the transmission of dengue virus between vectors and human hosts residing in the city (Patch 1) and its suburbs (Patch 2), while accounting for population mobility or people commuting between two zones under the two most probable scenarios. The first scenario is focused on one-way people flow from suburbs to the city, and the second one accounts for population movements in both directions; that is, between the city and its suburbs.
Further, we introduce the personal protection measures that can be assumed by (a share of) human populations residing in both patches in order to avoid mosquito bites, and thus to reduce the risk of infection with dengue virus. Such measures are modeled by two exogenous dynamical variables widely known as control functions, and these variables are patch-specific.
Using the framework of optimal control theory [4], we propose possible strategies for implementation of the personal protection measures by the residents of both patches and analyze them from the standpoint of cost-effectiveness. As a result, we identify an optimal strategy, which is capable of avoiding a greater number of human infections in both patches per unit of invested costs. This paper is organized as follows. In Section 2, we construct the two-patch model, parting from the simplest dengue transmission dynamics (that combines the "Susceptible-Infected-Susceptible" dynamics for human hosts with "Susceptible-Infected" dynamics for vectors and is usually referred to as "SIS-SI" transmission model) at each patch and introducing population mobility between patches. Subsequently, we calculate the basic reproduction number (R 2 0 ) associated with our model. In Section 3, control variables are introduced into the two-patch model, and the underlying optimal control problem is formulated. Section 4.1 contains numerical solutions of the optimal control problem under different scenarios of population mobility and application of personal protection measures by the residents of one or both patches. Finally, in Section 4.2 we perform the cost-effectiveness analysis in order to identify an optimal strategy for personal protection that allows avoidance of the maximal number of human infections in both patches per unit of invested costs. Section 5 provides some final remarks and conclusions.

1.
Human and vector populations remain essentially invariant in time.

2.
Populations of human hosts and vectors are homogeneous in terms of susceptibility, attraction, and exposure. 3.
Virus incubation periods within both humans and mosquitoes are ignored.

4.
Once infected, mosquitoes do not recover and die being infectious.

5.
Disease-induced death in humans or in vectors is not considered. 6.
Superinfection does not occur in either humans or mosquitoes; only susceptible or fully recovered individuals may get infected. 7.
Gradual acquisition of immunity in human hosts is ignored; they become susceptible immediately after recovery.
The first assumption is rather typical for the majority of epidemiological models, since it allows for their mathematical tractability. While it is quite natural to suppose that human populations N hi , i = 1, 2 remain "essentially invariant" in time (that is, the numbers of births and deaths among human residents is about the same per unit of time), this is not the case for mosquito populations N vi , i = 1, 2. In reality, mosquito populations exhibit daily variabilities induced by environmental and climatic factors, since their reproduction and longevity essentially depend upon current temperature and humidity.
For the sake of our modeling framework, nonetheless, we assume that both mosquito populations do not surpass their maximal sizes N vi , i = 1, 2, and thus focus on the "worst" scenario by supposing that mosquito populations in both patches are always close to their saturated sizes. It is worthwhile to note that the third and seventh assumptions are also related to the "worst" scenario: the transmissibility of the virus is not delayed either by its incubation within humans and vectors, or by the gradual acquisition of (cross)immunity in human hosts. These are the major limitations of our modeling approach, which we acknowledge to the reader.
It is worth pointing out that the last assumption holds for dengue-endemic areas where different serotypes of dengue virus (DENV1-DENV4) circulate simultaneously. This is exactly the case we are interested in, since it fits the conditions of the city of Cali, Colombia (to be further regarded as Patch 1), and its suburbs (to be regarded as Patch 2) in our two-patch model.
According to the Ross-Macdonald modeling framework, human and vector individuals at each patch are considered to be either susceptible or infected. Since the virus incubation periods are ignored, both human and mosquito individuals are considered infectious (i.e., capable of transmitting the virus) from the moment that they become carriers of the virus. Table 1 provides descriptions and notations for eight interacting sub-population groups of our model. For both patches, we assume that human residents in each zone have similar characteristics regarding their susceptibility, exposure, and recovery from the disease. Additionally, we suppose that both zones have the same climatic conditions, where mosquitoes have the same entomological characteristics. Table 2 presents notations and descriptions of the key parameters of the model, as well as their respective units. Generally speaking, the mosquito biting rate α expresses the number of bites taken by one female mosquito on human hosts in average per unit time. Parameters p h and p v describe the probabilities for a human host and vector to become infected after effective vector-to-human and human-to-vector contacts, respectively.
Human recovery rate γ indicates that a virus-carrying human host remains infectious during 1/γ days from the contact event until full recovery, and then becomes susceptible to acquiring infections caused by other DENV serotypes.
Mosquito mortality rate µ v is an inverse of its average lifespan, which corresponds to 1/µ v days. Here we ignore demographic changes (birth and death rates) in human populations and suppose that mosquito recruitment rate is equal to µ v in order to meet the first assumption of the model.
Parameters 0 ≤ q ij ≤ 1, i,j = 1,2, are elements of the residence time matrix, which is defined as in accordance with a previous study [7]. Each q ij expresses the fraction of time a person residing in Patch i spends, in average, in Patch j. In addition, the relationship 2 j=1 q ij = 1 holds for i = 1, 2.
Thus, people commuting between two patches can be modeled using the elements q ij of matrix (1) and applying the approach developed in a previous study [8] to a particular case of a SIS-SI metapopulation model, where human hosts and vectors are co-residents in both patches. It is essential to recall that mosquitoes do not relocate from their zones and only human hosts may commute between two zones. The dynamics of dengue transmission is schematically illustrated in Figure 1, and the differential equations describing our SIS-SI two-patch model are: Processes 2019, 7, x FOR PEER REVIEW 5 of 23 Since the populations of vectors and human hosts remain essentially invariant in both patches, the above model can be significantly simplified by using the following relationships: In other words, system (2) can be reduced to four equations with four state variables corresponding to sub-populations of infectious human hosts and vectors in both patches (i.e., second, fourth, sixth, and eighth equations from (2) with and replaced according to relationships in (3)). Since the populations of vectors and human hosts remain essentially invariant in both patches, the above model can be significantly simplified by using the following relationships: 6 of 24 In other words, system (2) can be reduced to four equations with four state variables corresponding to sub-populations of infectious human hosts and vectors in both patches (i.e., second, fourth, sixth, and eighth equations from (2) with s hi and s vi replaced according to relationships in (3)).
Moreover, taking the fractions of the four remaining sub-populations and performing the corresponding change of variables, the normalized version of model (2) can be obtained: Equation (5) shows that population mobility (expressed by means of the residence times q ij ) alters the fractions of infected human hosts and vectors residing in both patches. To see such alterations in a more comprehensive way, let us introduce the following matrix Its elements, P ij , satisfy the condition 2 i=1 P ij = 1 for j = 1, 2 and express the proportion of residents from Patch i effectively present in Patch j.
The matrix P defined by (6) is the matrix of residence times in terms of effective populations [8,9] and its components are generated by the elements of the residence time matrix (1).
In model (5), the parameter β v = αp v represents the rate of human-to-vector contact in both zones. On the other hand, for the vector-to-human contact rate, we must take into account that the number of people effectively present at each patch may vary, as we are considering the mobility of people between the two patches. Therefore, the vectorial density (i.e., an average number of vectors per one human host) must fit the context of the model in both patches. Instead of considering an average vectorial density at each patch (i.e., an average number of female mosquitoes per one human resident of the patch), we should focus on the effective vectorial density, which corresponds to the total number of female mosquitos N vi present in Patch i divided by the effective size of the human population in Patch i, that is, by q 1i N h1 + q 2i N h2 .
Thus, we can introduce the quantities N v1 q 11 N h1 + q 21 N h2 and N v2 q 12 N h1 + q 22 N h2 , which stand for effective vectorial densities in both patches. Using these quantities, the vector-to-human contact rates β h1 and β h2 become patch-specific, and can be written as follows: Finally, using the forms of components P ij and β hi defined in (6) and (8), respectively, the normalized version of the two-patch dengue transmission model (5) can be rewritten even more succinctly: It is worthwhile noting that all solutions X(t) = (I h1 (t), I v1 (t), I h2 (t), I v2 (t)) generated by any initial condition X(0) = X 0 ∈ Ω = [0, 1] 4 remain in Ω = [0, 1] 4 ⊂ R 4 + for all t ≥ 0. This conclusion is straightforward and follows from (4).

Calculation of Basic Reproductive Number R 0
In epidemiological modeling, the basic reproduction number (R 0 ) usually represents the average number of new secondary infections produced by one infected individual introduced to a fully susceptible population [10].
In the context of the two-patch model (9), the basic reproduction number has a more global sense. It expresses the average number of secondary infections produced by one infectious individual (mosquito or human host residing either in Patch 1 or Patch 2) when such an individual is introduced into a totally susceptible community comprising human and vector populations of both patches, where only people are allowed to commute between two patches.
To calculate the basic reproduction number for our two-patch dengue transmission model (8), we have followed the standard procedure described in a previous study [10]. First, we note that the state vector X = (I h1 , I v1 , I h2 , I v2 ) of the dynamical system (9) contains four infectious classes of vectors and human hosts, while X 0 = (0, 0, 0, 0) denotes the disease-free equilibrium of the system (9).
Let us rewrite the right-hand side of the model (9) as where F (X) ≥ 0 represents the rate of the disease transmission (i.e., rate of appearance of new infections), while V(X) ≥ 0 stands for the rate of the disease transition. The Jacobian matrices F = D X F (X) and V = D X V(X) of F (X) and V(X), evaluated in the disease-free equilibrium X 0 = (0, 0, 0, 0), have the following forms: According to a previous study [9], the basic reproduction number is defined by the largest eigenvalue (or spectral radius) of the next-generation matrix FV −1 , which is Omitting some tedious calculations, we arrive to the following form of the basic reproductive number for our model (9): where It should be emphasized that under the next-generation approach [10], the spectral radius of FV −1 defines the number of secondary infections generated per stage [11], and that dengue, as well as other vector-borne diseases, involves two stages of virus transmission from one human host to another (that is, human-to-vector and vector-to-human transmission stages). Therefore, the basic reproductive number obtained by this approach in the form (10) provides the average number of secondary infections for each transmission stage, without specifying the initial source of infection (from vector to human host or vice versa).
However, we are interested in the average number of secondary human infections produced by one infectious human host in an entirely susceptible host community, including human populations residing in both patches. Therefore, we should use the square of the "per stage" reproductive number defined by (10).
Therefore, the global basic reproductive number corresponding to our two-patch dengue transmission model (9) is given by the square of R 0 , that is, Remark 1. It is easy to see that in absence of people commuting between two patches we have where R 2 01 and R 2 02 denote patch-specific (or local) basic reproductive numbers corresponding to each particular patch. Namely, when q 11 = q 22 = 1 and q 12 = q 21 = 0, we have Consequently, Processes 2019, 7, 197 9 of 24 which stand for the basic reproductive numbers corresponding to the traditional Ross-Macdonald model for a single patch [5,6]). The latter also stays in line with the results of previous studies [7,9], where a two-patch dengue transmission model with additional compartments of "Exposed-Recovered" human hosts and "Exposed" vectors (known as SEIR(S)-SEI type and containing 14 state variables) was proposed and analyzed. Figure 2 displays the plot of the global reproductive number R 2 0 as a function of residence times q 12 and q 21 when Patch 2 (suburbs) has higher average vectorial density than Patch 1 (the city); that is, while all other parameters of the model remain the same for both patches. Consequently, which stand for the basic reproductive numbers corresponding to the traditional Ross-Macdonald model for a single patch [5,6]). The latter also stays in line with the results of previous studies [7,9], where a two-patch dengue transmission model with additional compartments of "Exposed-Recovered" human hosts and "Exposed" vectors (known as SEIR(S)-SEI type and containing 14 state variables) was proposed and analyzed. Figure 2 displays the plot of the global reproductive number ℛ 0 2 as a function of residence times 12 and 21 when Patch 2 (suburbs) has higher average vectorial density than Patch 1 (the city); that is, while all other parameters of the model remain the same for both patches. As shown in Figure 2, the global basic reproductive number ℛ 0 2 exhibits a slight rise as 12 increases, while 21 stays close to zero. In other words, more disease cases can be expected if the residents of suburbs (Patch 2) decide to spend a greater fraction of their time in the city (Patch 1) when the inverse flow (from Patch 1 to Patch 2) is very low or absent. Here we refer to the average As shown in Figure 2, the global basic reproductive number R 2 0 exhibits a slight rise as q 12 increases, while q 21 stays close to zero. In other words, more disease cases can be expected if the residents of suburbs (Patch 2) decide to spend a greater fraction of their time in the city (Patch 1) when the inverse flow (from Patch 1 to Patch 2) is very low or absent. Here we refer to the average number of expected secondary dengue infections produced by one infected human host in the total human population corresponding to both patches.
Conversely, the global basic reproductive number R 2 0 declines with respect to q 12 and regardless of q 21 . Therefore, a lesser number of disease cases can be expected if the city residents (Patch 1) decide to spend a greater fraction of their time in the suburbs (Patch 2). Thus, Figure 2 plainly illustrates that population mobility may either increase or reduce the number of new infections. This conclusion fully agrees with other results obtained by analysis of more sophisticated dengue transmission models [7,9].

Formulation of the Optimal Control Problem
To thwart the spread of the virus, intervention strategies could focus on reducing the number of contacts (or mosquito bites) between female mosquitoes and human hosts, for example, through the use of repellents, mosquito nets, and insecticide-treated clothing. Many scholars provide evidence that topical repellents constitute an important tool for prevention of infections caused by vector-borne pathogens, and may offer nearly 100% protection when applied as a spray, lotion, or cream directly on exposed skin [12]. A sole application of topic repellent may provide either a short-term protection (an hour or even less), in the case of plant-derived non-allergenic oils, or a prolonged complete protection (up to 12 h), in the case of commercially available chemical substances, such as DEET (or diethyltoluamide), that are rejected by some people due to their allergenic potential [13,14]. The efficacy of insect repellents depends not only on their active ingredients, but also on the application frequency.
There are some interesting studies that provide a solid background on the variability of the type of insect repellent used, factors influencing their effectiveness, possible attitudinal responses by individuals to their use, and a presentation of evidence regarding the correlation and some inconsistencies between the efficacy of mosquito bite reduction and mosquito-borne disease prevention [12][13][14][15].
To formalize the modeling of repellent application by the residents of both patches, we are going to introduce two control strategies that affect the rate of effective contacts of mosquitoes with susceptible and infected human individuals.
To incorporate these interventions into model (9), we define two exogenous variables that are time dependent and are independent of the other components in the model (9).
Let us denote u 1 (t) and u 2 (t), t ∈ [0, T] the control variables that act upon the rate of effective contacts in Patches 1 and 2, respectively. Here T > 0 is the finite time of control action and u 1 , u 2 : [0, T] → [0, 1] are two piecewise continuous real functions that represent the proportion of human residents from Patch 1 and Patch 2, respectively, that should apply repellent.
These control variables have the same goal that consists of reducing the number of bites taken by female mosquitoes on human residents in both zones.
Incorporating these two control measures; our model (9) can be written as: where η expresses the effectiveness of the repellent, which, in our case, is considered the same in both patches. We now formulate the problem of optimal control with the goal of minimizing the number of infected human hosts residing in both zones, while also minimizing the cost of applying the two controls. This goal can be defined by the following objective functional: where nonnegative constants A 1 , A 2 , A 3 , and A 4 represent the weight coefficients and express the priorities of each particular objective. They can also be considered as (relative) societal costs associated with each summand [6,16]. In this study, we assume that the marginal cost associated with reducing the number of infected human hosts (first two summands in (13)) is constant, whereas the marginal costs of control measures (last two summands in (13)) are not constant, that is, they depend on the underlying control actions. The resulting optimal control problem is to find the optimal strategies u * 1 (t) and u * 2 (t) that minimize the objective functional (13) subject to dynamical system (12) with assigned initial conditions 0 < I h1 (0) < 1, 0 < I h2 (0) < 1, 0 < I v1 (0) < 1, 0 < I v2 (0) < 1 (14) and under the following constraints imposed on both control variables In the above expression, u i , i = 1, 2, represent the maximum proportions of human hosts residing in Patch 1 (the city) and Patch 2 (the suburbs), respectively, that should use repellent for protection from the mosquito bites. From relationship (15), we define the set U = [0, u 1 ] × [0, u 2 ] of all admissible controls, which are piecewise continuous functions taking values in U for each t ∈ [0, T].

Existence of Optimal Controls
The problem of optimal control defined by (12)-(15) only makes sense when the disease persists. Therefore, we assume that the basic reproduction number corresponding to model (9) without controls is greater than one (that is, R 2 0 > 1).
Proof. First, we prove that all trajectories of the system (12) are bounded for any admissible pair of controls (u(t) 1 , u 2 (t)) ∈ U and for all t ∈ [0, T]. A mere glance at the equations of the system (12) reveals that their right-hand sides are Lipschitz continuous with respect to state variables X = (I h1 , I v1 , I h2 , I v2 ). Therefore, for any initial condition X(0) = X 0 ∈ Ω = [0, 1] 4 , there exists a unique solution X(t) = (I h1 (t), I v1 (t), I h2 (t), I v2 (t)) ∈ Ω corresponding to any admissible pair of controls (u 1 (t), u 2 (t)) ∈ U that remains in Ω for all t ≥ 0. Let G i (I h1 , I v1 , I h2 , I v2 ; u 1 , u 2 ), i = 1, 2, 3, 4 denote the right-hand side of the system (12). It is easy to see that the vector field G = (G 1 , G 2 , G 3 , G 4 ) satisfies the following conditions: where (X 1 , X 2 , X 3 , X 4 ) = (I h1 , I v1 , I h2 , I v2 ). From the above relationships, we can conclude that the controlled epidemiological dynamical model (12) is cooperative according to the definition given in a previous study [17]. Therefore, its trajectories X(t; u 1 , u 2 ) corresponding to any admissible pair of controls (u 1 (t), u 2 (t)) ∈ U are bounded from below by the so-called super-solution X(t) = X(t; 0, 0), and bounded from above by the so-called sub-solution X(t) = X(t; u 1 , u 2 ) [17].
The proof of existence of the optimal controls is based on the standard existence result given in a previous study [18] (see Theorem 4.1 and Corollary 4.1 on p. 68). In this context, let us emphasize that all the hypotheses established in this result are satisfied for the optimal control (12)-(15) problem, namely:

2.
The integrand of the objective functional (13) is convex with respect to the state and control variables.

3.
The right-hand side of the dynamical system (12) is linear with respect to both controls u 1 (t) and u 2 (t).

4.
By definition, the set of admissible controls U is compact.

Characterization of Optimal Controls
The optimal control problem defined by (12)-(15) can be formally solved via direct application of Pontryagin's maximum principle [4]. To proceed, we define the Hamiltonian function that depends on four state variables X = (I h1 , I h2 , I v1 , I v2 ), two control variables u 1 , u 2 , and four adjoint variables Λ = (λ 1 , λ 2 , λ 3 , λ 4 ) or co-states. The adjoint vector depends on time t ∈ [0, T] and satisfies the so-called adjoint system of differential equations where the endpoint condition at t = T is usually referred to as the transversality condition. The components λ i (t), i = 1, 2, 3, 4 of Λ(t) are the so-called "shadow prices" associated with the respective state variables. Generally speaking, they express the change in the objective value J(u 1 , u 2 ) calculated for the optimal solutions u * 1 (t), u * 2 (t) when the dynamic constraints represented by the controlled system (12) are relaxed by one unit [4].

Remark 2.
Necessary conditions for optimality (19) of u * i (t), i = 1, 2 can be written as In the above relationships, the positive terms A 3 u * 1 and A 4 u * 2 express the marginal costs of two control actions u * 1 (t) and u * 2 (t), respectively, whereas the negative terms of the form −η[. . .] refer to the marginal benefits of u * 1 (t) and u * 2 (t).
Thus, from the economics standpoint, Pontryagin's maximum principle states that the marginal cost of the optimal strategy u * i (t), i = 1, 2 must be equal to the marginl benefit rendered by this strategy (that is, ∂H ∂u i = 0). It also follows from (20) that if the marginal cost of u * i (t), i = 1, 2 is higher than its marginal benefit (that is, ∂H ∂u i > 0), then it is optimal not to use this strategy at all, and we have u * i (t) = 0, i = 1, 2. Alternatively, if the marginal cost of u * i (t), i = 1, 2 is lower than its marginal benefit (that is, ∂H ∂u i < 0), then it is optimal to use all available resources, and we have u * i (t) = u i , i = 1, 2. It is worthwhile to point out that transversality conditions λ i (T) = 0, i = 1, 2, 3, 4 from the adjoint system (17) guarantee that u * 1 (T) = u * 2 (T) = 0, meaning that both control actions must be suspended by the end of the observation period [0, T]. This property is attributed to the continuity of adjoint variables λ i (t) = 0, i = 1, 2, 3, 4 when t → T and to the closed forms of optimal control characterizations (18). Additionally, the optimal controls given by (18) are minimizers, because the Hessian matrix D 2 u ij H of the Hamiltonian function H taken with respect to control variables and evaluated in (u * 1 (t), u * 2 (t)) is constant and positive definite: Using the closed forms of optimal controls (18), the original optimal control problem (12)-(15) can be reduced to a two-point boundary value problem (known as "optimality system"), which is composed of eight differential equations with eight endpoint conditions, namely: • Four direct equations (12) plus four adjoint equations (17), where (u 1 , u 2 ) are replaced by their characterizations (18); • Four initial conditions (14) specified at t = 0, plus four transversality conditions from (17) specified at t = T.
The solution for the optimality system is a rather challenging task, due to its non-linearity and high dimension. Therefore, it can only be solved numerically, and the following section provides the underlying results.

Numerical Simulations
To solve the optimality system resulting from the optimal control problem (12)-(15), we have employed the software package GPOPS-II (Next-Generation Optimal Control Software, Version 2.4, RP Optimization Research LLC, Gainesville, FL, USA), which implements the numerical technique based on the direct orthogonal collocation [19].
Numerical values of parameters for the optimal control problem formulated in the previous section are given in Table 3. The greater part of these values are borrowed from published papers [5,6,9] dealing with dengue epidemic studies conducted in the city of Cali and its suburbs. Some parameter values have been assumed or estimated.
Let us provide some explanations concerning the entries of Table 3, which are assumed or estimated. There is no published data regarding estimations of the average vectorial density in Patch 2 (corresponding to suburbs), and we have assumed it is about two times higher than that in Patch 1 (the city) for the following reasons:

1.
Smaller towns and settlements in Colombia are known to have problems with sanitation and intermittent water supply. Therefore, water storage tanks kept by suburban residents contribute to mosquito proliferation.

2.
In contrast to major Colombian cities located in dengue endemic areas (such as Cali-Patch 1), vector control measures in suburban areas are irregular or absent [20]. Since we do not possess information regarding the total sizes of vector populations present in both patches, these entries have been estimated as N v1 = 1.59691 × N h1 and N v2 = 3 × N h2 , respectively. The upper values u 1 and u 2 for control variables (representing the maximum proportion of residents in both patches who use repellent for personal protection) was set to unity in order to fit the assumption of homogeneous mixing of all human hosts residing in both patches. On the other hand, we set η = 0.7 despite knowing that most common repellents can reduce the number of mosquito bites by up to 95% when applied with adequate frequency [21]. In this manner, we can admit that not all human hosts may apply repellents, and that the number of effective contacts (mosquito ⇔ human) can be reduced by 70% at most.
Knowing that dengue epidemics in Cali and its suburbs usually lasts for 1-3 months [5,6,9], the observation period T > 0 was set to 30 days.
To define the values of weights A i , i = 1, 2, 3, 4 in the objective functional (13), we used the arguments provided in [6]. The values of A 1 and A 2 are associated with the total daily cost of one infected individual (accounting for treatment and temporary disability leave) residing at either patch. Since the societal cost of one dengue case in Colombia is $600 [6], and considering that it takes 10 days to recover (1/γ), it was supposed that A 1 = A 2 = 60.
The values of A 3 and A 4 represent the estimated expenditure for educational campaigns aiming to motivate the human population (residing in both patches) to take personal protection measures. From previous studies [6,16], the unit cost estimated for these highly efficient campaigns is approximately 50 times less than the total medical and social unit cost of one infected person when only one patch is considered. However, when dealing with two patches, we should account for different human population sizes of these patches (N h1 > N h2 ), and make the underlying adjustments. Thus, we assume As a result, we can assess the total cost related to implementation of the control measures u i (t), t ∈ [0, T], i = 1, 2 for both patches in accordance with a previous study [16], making use of the marginal instantaneous costs: To characterize the current level of the disease in our two-patch system, we assign the initial conditions (expressed in proportions) to all four state variables of the dynamical system (12): The first condition is realistic and comes from a previous study [9], while the second one is chosen to be compatible with current statistics on yearly dengue incidence rates in Colombia [20]. According to this source, major Colombian cities report 2-4 times more dengue cases per thousand inhabitants than their suburban areas. Since there is no reliable data regarding the fractions of infectious vectors in both patches, we suppose the fractions of infectious human hosts residing in underlying patches to be proportional (the last two initial conditions).
Finally, we should assign numerical values to all parameters dealing with population mobility (i.e., define the elements of time residence matrix Q given by (1)). Before considering possible scenarios of population mobility between two patches, it is useful to revisit Figure 2, which displays the global basic reproductive number R 2 0 (see Formula (11)) as a function of residence times q 12 and q 21 , while keeping in mind that q 11 = 1 − q 12 and q 22 = 1 − q 21 . This figure is plotted using the values of parameters corresponding to the city of Cali (Patch 1) and its suburbs (Patch 2), which are given in Table 3. This figure also helps visualize that people commuting between two patches may enhance or reduce the value of R 2 0 , and thus have an impact on the disease propagation. However, we seek to keep our model as close as possible to realistic situations regarding people commuting between the city (Patch 1) and its suburbs (Patch 2). We are also interested in exploring the effects of the personal protection measures on the transient dengue morbidities in both patches and to estimate the efficiency of these control measures.
Thus, we assume that Cali residents stay mostly at their home patch (q 12 = 0) or rarely commute to another patch (q 12 > 0 is rather small), while the opposite is true for the suburban residents (q 21 > 0). This assumption leads to two basic types of population mobility between the two patches:

1.
One way: Suburban residents (Patch 2) spend, on average, 40% of their time in the city (Patch 1), while the city residents do not commute to the suburbs: The One way option is rather realistic since a significant share of suburban residents commute to the city for work, study, shopping, etc., on a regular basis. However, some companies are starting to move their headquarters and plants to the suburban areas; therefore, the Both ways option may also become realistic in the near future.
Before proceeding to the numerical solution to the optimal control problem (12)-(15), let us briefly revise the outcome of our two-patch dengue transmission model without control intervention, i.e., with u 1 (t) = u 2 (t) = 0, t ∈ [0, T], and under both options for commuting (One way and Both ways). For that purpose, it is helpful to introduce two auxiliary variables, C hi (t), i = 1, 2 (linked to each patch), known as cumulative incidences [9]. These variables express the cumulative proportion of all human infections occurring in each patch during the observation period. In mathematical terms, these additional variables are defined by the following differential equations with corresponding initial conditions Figure 3 presents the cumulative incidences C hi (t), i = 1, 2 for both commuting options (One way and Both ways) obtained by numerical integration of six differential equations, (12) and (22), with = (1 − 1 )(1 − ℎ1 )( 11 ℎ1 1 + 12 ℎ2 2 ), ℎ1 (0) = ℎ1 (0), ℎ2 = (1 − 2 )(1 − ℎ2 )( 21 ℎ1 1 + 22 ℎ2 2 ), ℎ2 (0) = ℎ2 (0). Figure 3 presents the cumulative incidences ℎ ( ), = 1,2 for both commuting options (One way and Both ways) obtained by numerical integration of six differential equations, (12) and (22)  It is worth noting that One way option generates less infections in Patch 1 and more infections in Patch 2 than Both ways option (see Figure 3). This outcome can be explained using the concept of effective vectorial density (see Equation (7)). In effect, under the One way option, Patch 1 receives visitors, and this reduces its average vectorial density (i.e., number of vectors per one resident of It is worth noting that One way option generates less infections in Patch 1 and more infections in Patch 2 than Both ways option (see Figure 3). This outcome can be explained using the concept of effective vectorial density (see Equation (7)). In effect, under the One way option, Patch 1 receives visitors, and this reduces its average vectorial density (i.e., number of vectors per one resident of Patch 1) to a lesser value of effective vectorial density (i.e., number of vectors per one human host effectively present in Patch 1): Therefore, residents of Patch 1 receive less mosquito bites under the One way option, and this is reflected in a reduced number of infections occurring in Patch 1. In case of Both ways option, Patch 2 exhibits a reduction in the number of infections, since it holds that The latter decreases the number of infectious bites received by the residents of Patch 2 under Both ways option. To corroborate the above deduction, we can also calculate the absolute cumulative numbers of infections C hi (T) = N hi × C hi (T), i = 1, 2 acquired during T = 30 days by human hosts in both patches: These estimations will help us in evaluation of preventive control policies, which are further obtained by the numerical solution of the optimal control problem (12)-(15), amended with auxiliary variables C hi (t), i = 1, 2 introduced in (22).
Subsequently, we can define six scenarios resulting from combinations of three strategies (S1, S2, S3) described above with two options of commuting (One way and Both ways). These scenarios can be summarized as follows: Finally, the problem of optimal control (12)-(15), (22) is solved numerically (using GPOPS-II software package) under six different scenarios (23) and its solutions u * i (t), C * hi (t), i = 1, 2 are displayed in Figure 4 (Scenarios 1 and 4) and Figure 5 (Scenarios 2 and 5). It should be noted that we omit here the plots for Scenarios 3 and 6 because they bear almost no difference to plots presented in Figures 4 and 5 for each patch, and this difference is visually undetectable.    On the other hand, this small difference can be viewed by calculating the absolute cumulative numbers of infections C * hi (T) = N hi × C * hi (T), i = 1, 2 acquired by the human hosts in both patches under the personal protection measures u * 1 (t), u * 2 (t) corresponding to each scenario. Additionally, we can estimate the number of infections that can be averted in each patch whenever the residents of one or both patches use personal protection measures in accordance with the optimal control strategy given by u * i (t), i = 1, 2. Let P i u * 1 (t), u * 2 (t) , i = 1, 2 express the absolute number of human infections prevented by the optimal control policy u * 1 (t), u * 2 (t) , t ∈ [0, T] in Patch i, i = 1, 2 for each scenario given in (23). Then we have where C hi (T) corresponds to the total number of infections acquired by the residents of Patch i, i = 1, 2 during the observation period [0, T] without control (i.e., supposing u 1 (t) = u 2 (t) = 0, t ∈ [0, T]), while C * hi (T) stands for the total number of human infections in the same patch under the optimal control policy u * 1 (t), u * 2 (t) , t ∈ [0, T]. Table 4 presents the summary of results. Scenarios 1 and 4 deal with application of the personal protection measure only by the city residents (strategy S1 in (23)), and the bottom row of Figure 4 provides the optimal profiles of u * 1 (t), t ∈ [0, T] under two mobility options: One way (left chart) and Both ways (right chart). In this case, we have u * 2 (t) = 0, t ∈ [0, T], meaning that suburban residents do not use personal protection. The structure of u * 1 (t), t ∈ [0, T] in the lower row of Figure 4 indicates that until approximately day 25, all residents of Patch 1 should protect themselves with repellent; then, the fraction of human hosts who use the measures of personal protection gradually decrease towards zero during the last 5 days of the observation period.
The effect of optimal control on the reduction of dengue morbidity is illustrated on the two upper charts of Figure 4 (for Patch 1) and the two upper charts of Figure 5 (for Patch 2) under two mobility options: One way and Both ways (left and right charts, respectively). For Patch 1 (the city), the impact of control action is clearly visible (compare black and red dashed curves in Figure 4). However, the application of protective measures only by the residents of Patch 1 has very little effect on the disease suppression among the residents of suburbs (Patch 2, see Figure 5). Yet, Scenario 4 (mobility option: Both ways) yields more prevented human infections among suburban residents than Scenario 1 (mobility option: One way). The estimated numbers of averted human infections given in Table 4 confirms this conclusion.
Scenarios 2 and 5 deal with application of the personal protection measure only by the suburban residents (strategy S2 in (23)), and the bottom row of Figure 5 provides the optimal profiles of u * 2 (t), t ∈ [0, T] under two mobility options: One way (left chart) and Both ways (right chart). In this case, we have u * 1 (t) = 0, t ∈ [0, T], meaning that city residents do not use personal protection. The structure of u * 2 (t), t ∈ [0, T] in the lower row of Figure 5 indicates that until approximately day 28, all residents of Patch 2 (both commuters and non-commuters) should protect themselves with repellent. Then, the fraction of human hosts who use the measures of personal protection abruptly decreases towards zero during the last 2 days of the observation period.
The effect of optimal control on the reduction in dengue morbidity is illustrated in the two upper charts of Figure 4 (for Patch 1) and the two upper charts of Figure 5 (for Patch 2) under two mobility options: One way and Both ways (left and right charts, respectively). For Patch 2 (suburban areas), the impact of control action is clearly visible (see the difference between black dashed and blue solid curves in Figure 5).
However, the application of protective measures only by the residents of Patch 2 has very little effect on the disease suppression among the city residents (Patch 1, see Figure 4), even though Scenario 5 (mobility option: Both ways) yields more prevented human infections among city residents than Scenario 2 (mobility option: One way). The latter is reflected in the estimated numbers of averted human infections given in Table 4. Scenarios 3 and 6 deal with application of the personal protection measure for the residents of both patches (strategy S3 in (23)), and the optimal control profiles u * 1 (t), u * 2 (t) , t ∈ [0, T] corresponding to these scenarios have almost the same structures as u * 1 (t) and u * 2 (t) plotted in the lower charts of Figures 4 and 5 under two mobility options: One way (left charts) and Both ways (right charts). Naturally, and quite expectedly, the effect of u * 1 (t), u * 2 (t) on the disease suppression in both patches is very explicit and can be visualized in the upper chart of Figures 4 and 5 for two mobility options (see the difference between black dashed and green dotted curve in Figures 4 and 5). Additionally, the estimations given in Table 4 indicate that u * 1 (t), u * 2 (t) are capable of preventing between 50% and 60% of human infections in both patches.
Thus, from the standpoint of potential benefits, the control intervention policies u * 1 (t), u * 2 (t) corresponding to Scenarios 3 and 6 are definitely the best for each underlying option of population mobility.
However, in order to decide which control strategy (S1, S2, S3) is the most efficient (or cost-effective), we have to take into account not only the benefits rendered by each strategy (expressed in terms of averted human infections in both patches) but also the underlying costs.

Cost-Effectiveness Analysis
The cost-effectiveness analysis is an economic assessment tool that aims to compare the costs and the effects of two or more control intervention policies in order to determine which particular policy renders higher benefits per unit cost. In the healthcare management, the principal measure of cost-effectiveness is the indicator known as Average Cost-Effectiveness Ratio (or ACER) that expresses the cost per one averted disease case, and can be formally defined [22] in the following way: The total costs of the optimal control policies u * 1 (t), u * 2 (t) , t ∈ [0, T] corresponding to each scenario described in (23) can be estimated using the underlying formula for their marginal costs (21). The benefit of each control strategy is then obtained by summing up the estimated number of human infections averted by applying this strategy in both patches, that is, Benefit u * 1 (t), u * 2 (t) = P 1 u * 1 (t), u * 2 (t) + P 2 u * 1 (t), u * 2 (t) where P 1 and P 2 are calculated according to (24). Table 5 presents the summary of results for three types of control strategies (S1, S2, S3) and two options for population mobility (One way and Both ways). A mere glance at the last column of Table 5 reveals that Scenarios 2 and 5 possess the lowest ACER for each mobility option. This suggests that the control intervention policies corresponding to these two scenarios render higher benefits per unit cost invested. Therefore, strategy S2 is the most efficient (or cost-effective) under both mobility options (One way and Both ways).
In other words, it is more "cost effective" that all suburban residents (Patch 2), commuters, and non-commuters, susceptible and infected (possibly asymptomatic), protect themselves with repellents.
This conclusion, obtained from the mathematical standpoint, can also be attributed to the following "common-sense" factors:

1.
The human population size of Patch 2 is smaller than that of Patch 1. Therefore, the cost associated with educational campaigns promoting personal protection measures among residents of Patch 2 is lower than the cost of similar campaigns targeting residents of Patch 1.

2.
Both average and effective vectorial density is higher in Patch 2 than in Patch 1. Therefore, residents of Patch 2 are at greater risk of receiving mosquito bites, and reducing their contact with vectors (bites) by applying repellent should have a greater effect on residents of Patch 2 than on residents of Patch 1 (where both average and effective vectorial density is lower).

Conclusions
In this paper, we have addressed the role of personal protection measures on dengue transmission, while considering people commuting between two zones, both located in a dengue-endemic area (city of Cali, Colombia, and its suburbs). We have tried to model two realistic situations regarding daily commuting and to analyze strategies of personal protection aiming to reduce the number of contacts between human hosts and vector transmitters of the disease (mosquito bites).
We have considered the modeling framework of an SIS-SI for the dynamics together with the so-called Lagrangian approach [8]. Under this approach, the dispersal of human hosts between two patches caused by daily commuting was modeled using the fractions of time that human individuals spend at each patch, which are also known as residence times [7][8][9].
Our two-patch model for dengue transmission preserves the key properties of more sophisticated models where multiple patches are considered [8], or where the dengue transmission modeling is more detailed (dengue dynamics of SEIRS-SEI type, [7,9]). Namely, it shows that daily mobility affects the disease transmission in both patches quite differently, and such differences can be explained using the concepts of effective human population sizes [8] and effective vectorial densities [9].
Using the two-patch dengue transmission model, we have developed an optimal control framework in order to design patch-specific strategies for personal protection of human hosts from mosquito bites that are capable of mitigating the disease transmission in both patches.
Optimal policies for application of personal protection measures have been designed under six different scenarios resulting from combinations of two options for population mobility (One way and Both ways) with three strategies (protection used only by the city residents, only by suburban residents, and by both human populations).
All designed policies have been evaluated from the standpoint of cost-effectiveness, where the potential benefits of each policy was assessed via the number of prevented human infections. As a result of the cost-effectiveness analysis, it was concluded that application of personal protection measures by all suburban residents (both susceptible and infected, commuters, and non-commuters) renders higher benefits (expressed via number of avoided human infection in both patches) per unit of cost invested in promotion of personal protection measures.
This sends a clear message to local healthcare authorities that may help them in appropriate preparation and scheduling of the educational campaigns, which seek to motivate local populations in daily use of personal protection from mosquito bites. If there are no sufficient funds to hold a large-scale campaign both in the city and its suburbs, then it is better to strongly promote the use of personal protection measures among the residents of suburban areas, including free distribution of repellent or mosquito nets to the suburban population.
The role of personal protection via the use of repellent may become even more visible if combined with other vector control measures, such as elimination of mosquito breeding sites inside and around households. However, implementation of successful programs based on personal protection would require cooperation and interaction between the residential communities and healthcare authorities, which is not easy to reach. A starting point could be to organize educational campaigns and to inform the local residents that combinations of personal protection measures are likely to be more effective than single methods.