Abstract

We propose a metapopulation model for malaria with two control variables, treatment and prevention, distributed between 𝑛 different patches (localities). Malaria spreads between these localities through human travel. We used the theory of optimal control and applied a mathematical model for three connected patches. From previous studies with the same data, two patches were identified as reservoirs of malaria infection, namely, the patches that sustain malaria epidemic in the other patches. We argue that to reduce the number of infections and semi-immunes (i.e., asymptomatic carriers of parasites) in overall population, two considerations are needed, (a) For the reservoir patches, we need to apply both treatment and prevention to reduce the number of infections and to reduce the number of semi-immunes; neither the treatment nor prevention were specified at the beginning of the control application, except prevention that seems to be effective at the end. (b) For unreservoir patches, we should apply the treatment to reduce the number of infections, and the same strategy should be applied to semi-immune as in (a).

1. Introduction

Malaria is a mosquito-borne infection caused by protozoa of the genus plasmodium. Parasites are transmitted indirectly from humans to humans by the bite of infected female mosquitoes of the genus Anopheles. Malaria is a public health problem for tropical countries, which has negative impacts on development. The fight against mosquitoes passes through the draining of marshes or conversion to running water and elimination of stagnant water especially around houses. These measures are difficult to apply where health facilities are inadequate [1].

Mathematical models coupled to microeconomic concepts can be applied to malaria control using the theory of optimal control [2–5]. The latter has already been used to discuss strategies to reduce or eradicate other diseases such as chronic myeloid leukemia [6], AIDS [7, 8], tuberculosis [9, 10], smoking [11], West Nile virus [12], and Chikungunya disease [13].

Human movements play a key role in the spatiotemporal spread of malaria [14, 15]. Models have already been proposed to study malaria spread, but neither control variables such as treatment and prevention nor human mobility were considered [16, 17]. In [18], a model with control variables (synergy prevention and treatment) was developed for malaria spread. Because human population was assumed to be motionless, their model is only applicable within a small geographical region. A model that takes into account human mobility was developed in [19] to analyze the impact of human migration within 𝑛 geographical patches (localities) on the malaria spread.

Our model is not an extension of the model developed in [18] but rather the one developed in [19].

First, in [19], the authors have considered two infectious classes in the human population: infectious and semi-immune individuals (i.e., asymptomatic carriers). This consideration is very important because an experimental evidence showed that 60–90% of humans in endemic area are semi-immune [1, 16, 17]. Introduction of semi-immunes in their model presents the difficulty to introduce the control variable, mainly the treatment variable because treated individuals may become either susceptible or semi-immune depending on the type of the used drugs. In this paper, we will introduce a parameter which has a role of regulation denoted by πœƒ to derive a biological meaningful model (see Figure 1).


Figure 1 
A conceptual mathematical model for malaria transmission involving human hosts and vector mosquitoes in each patch 𝑖 , 𝑖 = 1 , … , 𝑛 . The dotted arrow shows the direction of the transmission from humans to mosquitoes (through infectious humans to susceptible mosquitoes) or from mosquitoes to humans (through infectious mosquitoes to susceptible humans); 𝑒 2 𝑖 and 𝑣 2 𝑖 represent the prevention and treatment control over time, respectively. Ξ¦ 2 𝑖 = ( 1 βˆ’ 𝑒 2 𝑖 )  Ξ¦ 2 𝑖 and Ξ¦ 2 𝑖 βˆ’ 1 = ( 1 βˆ’ 𝑒 2 𝑖 )  Ξ¦ 2 𝑖 βˆ’ 1 represent the force of infection from mosquitoes to humans and from humans to mosquitoes, respectively, where  Ξ¦ 2 𝑖 and  Ξ¦ 2 𝑖 βˆ’ 1 are defined in (2.2). The other parameters are described in Table 1.

Second, the model developed in [19] is of the metapopulation type that considers the explicit movement of humans between many patches. Because the mathematical analysis of their model has provided a methodology to identify the spatial reservoirs of malaria infection (i.e., the patches that sustain malaria epidemic in the other patches) based on the theory of the type reproductive number, our main objective in this paper is to extend their results by introducing control variables (treatment and prevention) within each patch. With these innovations, the simulations identify the best strategy of control and answer the following question: what control should be used when the patch is (or is not) a reservoir? This question is not trivial because the infectious individuals can migrate in all the other patches.

The paper is structured as follows: in Section 2, we summarize the main points of the metapopulation model from [19] by introducing prevention and treatment controls. Furthermore, we show that it is mathematically well posed. Section 3 includes the formulation of the objective function with the discount rate, and properties of optimal control existence follow its characterization. In Section 4, we present the results of simulations and discussion for three connected patches by migration according to the type of reservoir of infections. The last section includes the conclusion and perspectives.

2. Mathematical Modelling

2.1. Model Description

In this section, a metapopulation model with control variables (prevention and treatment) is developed. In the sequel, we use even and odd index to represent the human and mosquito variables, respectively, as in [1]. Within each small patch, the human hosts are split into three subclasses: susceptible 𝑆2𝑖, infectious 𝐼2𝑖, and semi-immune 𝑅2𝑖.𝑁2𝑖(𝑑)=𝑆2𝑖(𝑑)+𝐼2𝑖(𝑑)+𝑅2𝑖(𝑑) denotes the total size of the human population in the patch 𝑖 at time 𝑑. The mosquito population is split into two subclasses: susceptible 𝑆2π‘–βˆ’1 and infectious 𝐼2π‘–βˆ’1 in the patch 𝑖,𝑖=1,…,𝑛. The total size of the mosquito population is denoted by 𝑁2π‘–βˆ’1(𝑑)=𝑆2π‘–βˆ’1(𝑑)+𝐼2π‘–βˆ’1(𝑑) at time 𝑑. π‘β„Žβˆ‘(𝑑)=𝑛𝑖=1𝑁2𝑖(𝑑) and π‘π‘£βˆ‘(𝑑)=𝑛𝑖=1𝑁2π‘–βˆ’1(𝑑) denote the total size of the human and mosquito population for the complete system, respectively, at any time 𝑑 (See Figure 1). The model with control reads as follows: for all 𝑖=1,…,𝑛,𝑑𝑆2𝑖𝑑𝑑=Ξ›2𝑖+𝛽2𝑖𝑅2𝑖+𝜌2𝑖𝐼2π‘–βˆ’ξ€·πœ‡2𝑖+Ξ¦2𝑖𝑆2𝑖+πœƒπœ†2𝑖𝑣2𝑖𝐼2𝑖+𝑛𝑗=1ξ€·π‘šπ‘†π‘–π‘—π‘†2π‘—βˆ’π‘šπ‘†π‘—π‘–π‘†2𝑖,(2.1a)𝑑𝐼2𝑖𝑑𝑑=Ξ¦2𝑖𝑆2π‘–βˆ’πœ–2𝑖𝐼2π‘–βˆ’πœ†2𝑖𝑣2𝑖𝐼2𝑖+𝑛𝑗=1ξ€·π‘šπΌπ‘–π‘—πΌ2π‘—βˆ’π‘šπΌπ‘—π‘–πΌ2𝑖,(2.1b)𝑑𝑅2𝑖𝑑𝑑=𝛼2𝑖𝐼2π‘–βˆ’π›Ώ2𝑖𝑅2𝑖+(1βˆ’πœƒ)πœ†2𝑖𝑣2𝑖𝐼2𝑖+𝑛𝑗=1ξ€·π‘šπ‘…π‘–π‘—π‘…2π‘—βˆ’π‘šπ‘…π‘—π‘–π‘…2𝑖,(2.1c)𝑑𝑆2π‘–βˆ’1𝑑𝑑=Ξ›2π‘–βˆ’1βˆ’πœ‡2π‘–βˆ’1𝑆2π‘–βˆ’1βˆ’Ξ¦2π‘–βˆ’1𝑆2π‘–βˆ’1,(2.1d)𝑑𝐼2π‘–βˆ’1𝑑𝑑=Ξ¦2π‘–βˆ’1𝑆2π‘–βˆ’1βˆ’πœ‡2π‘–βˆ’1𝐼2π‘–βˆ’1,(2.1e)where πœ–2𝑖=𝛼2𝑖+𝛾2𝑖+𝜌2𝑖+πœ‡2𝑖 and 𝛿2𝑖=𝛽2𝑖+πœ‡2𝑖. Initial conditions are assumed to satisfy 𝑆2𝑖(0)>0,𝑆2π‘–βˆ’1(0)>0,𝐼2𝑖(0)β‰₯0,𝑅2𝑖(0)β‰₯0, and 𝐼2π‘–βˆ’1(0)β‰₯0 for 𝑖=1,…,𝑛.

In the above model, Ξ¦2𝑖 and Ξ¦2π‘–βˆ’1 denote the force of infection from mosquitoes to humans and from humans to mosquitoes, respectively. Therefore, infection only involves vectors or hosts present in the patch; there is no between-patch infection. These forces of infection are modeled to take into account the prevention as in [18] as follows: Ξ¦2𝑖=(1βˆ’π‘’2𝑖)Φ2𝑖 and Ξ¦2π‘–βˆ’1=(1βˆ’π‘’2𝑖)Φ2π‘–βˆ’1 where Φ2𝑖=Μƒπ‘Ž2π‘–βˆ’1Μƒπ‘Ž2𝑖𝑁2π‘–βˆ’1Μƒπ‘Ž2π‘–βˆ’1𝑁2π‘–βˆ’1+Μƒπ‘Ž2𝑖𝑁2π‘–πœŽ2π‘–βˆ’1,2𝑖𝐼2π‘–βˆ’1𝑁2π‘–βˆ’1,Φ2π‘–βˆ’1=Μƒπ‘Ž2π‘–βˆ’1Μƒπ‘Ž2𝑖𝑁2π‘–Μƒπ‘Ž2π‘–βˆ’1𝑁2π‘–βˆ’1+Μƒπ‘Ž2𝑖𝑁2π‘–ξ‚΅πœŽ2𝑖,2π‘–βˆ’1𝐼2𝑖𝑁2𝑖+𝜎2𝑖,2π‘–βˆ’1𝑅2𝑖𝑁2𝑖(2.2) are defined in [17] for one patch, Μƒπ‘Ž2𝑖,Μƒπ‘Ž2π‘–βˆ’1,𝜎2π‘–βˆ’1,2𝑖,𝜎2𝑖,2π‘–βˆ’1βˆˆβ„+,𝜎2𝑖,2π‘–βˆ’1βˆˆβ„βˆ—+. The infection force Ξ¦2𝑖 and Ξ¦2π‘–βˆ’1 depends on the individuals within patch 𝑖 and not in another patch 𝑗≠𝑖: infection that only involves those individuals (vectors or hosts) present in the patch is no between-patch infection.

𝑒2𝑖 is the prevention effort for humans which reduces the infection rate with a failure probability 1βˆ’π‘’2𝑖 if prevention controls are introduced. The control function 𝑒2𝑖 represents time-dependent efforts of prevention on human and practiced on a time interval [0,𝑇]. Prevention could come from surveillance, treating vector-breading ground, and reducing vector-host contacts. Note that when 𝑒2𝑖=0, then Ξ¦2𝑖 and Ξ¦2π‘–βˆ’1 correspond to those used in [19].

πœ†2𝑖𝑣2𝑖 is the per capita recovery rate of humans. 0β‰€πœ†2𝑖≀1 is the proportion of effective treatment of humans. The control function 𝑣2𝑖 represents the measure of the rate at which infected humans are cured by drugs or vaccination on a time interval [0,𝑇].

The difference between the effects of drugs is related to the fact that they act at different stages of the parasite-cell mutation in the human body. There are drugs that act against preerythrocytic stages, against the asexual blood stages, and against antigens of sexual stages that prevent fertilization in the stomach of the mosquito [20]. We thought that each drug has its own effect on the mode of healing. Therefore, there are drugs that favor the individuals who immunize quickly, while there are others that favor the total healing without being immunized. We then introduced the model of parameter πœƒ, which regulates the process. πœƒ is the probability that the treated infectious humans pass the sensitive compartment, and 1βˆ’πœƒ is the probability to pass the semi-immune compartment. When using treatments that immunize the majority of patients, πœƒ tends to 1, and these patients will go into the compartment of the semi-immune. Otherwise, πœƒ tends to 0, and the patients will move into the susceptible compartment.

We also provide an insight into the major assumptions made in the original model in [19] as follows: disease-induced death rate of semi-immune was assumed to be negligible because this host acquires some immunity. Human mobility from one patch to another was considered, while immigration of mosquitoes was neglected because they can explore only a few kilometers during their lives. During the travel, humans do not change status. π‘šπœ‹π‘–π‘—,πœ‹=𝑆,𝐼,𝑅 denote the constant rate of travel of humans from an area 𝑗 to an area 𝑖 for all 𝑖≠𝑗 with π‘€πœ‹=[π‘šπœ‹π‘–π‘—], and πœ‹=𝑆,𝐼,𝑅 denote the travel rate matrices. The matrices 𝑀𝑆 were assumed to be irreducible and π‘šπœ‹π‘–π‘–=0,πœ‹=𝑆,𝐼,𝑅;𝑖=1,…,𝑛.

Table 1 summarizes the parameters and their biological description that will be used in the metapopulation model.

Table 1 
Parameters for the model described in any patch 𝑖,𝑖=1,…,𝑛.

By adding up (2.1a)–(2.1c) and (2.1d)-(2.1e), we obtain expressions for the total human and mosquito populations, respectively, in patch 𝑖=1,…,𝑛: 𝑑𝑁2𝑖𝑑𝑑=Ξ›2π‘–βˆ’πœ‡2𝑖𝑁2π‘–βˆ’π›Ύ2𝑖𝐼2𝑖+ξ“πœ‹=𝑆,𝐼,𝑅𝑛𝑗=1π‘šπœ‹π‘–π‘—πœ‹2π‘—βˆ’π‘›ξ“π‘—=1π‘šπœ‹π‘—π‘–πœ‹2𝑖ξƒͺ,𝑑𝑁2π‘–βˆ’1𝑑𝑑=Ξ›2π‘–βˆ’1βˆ’πœ‡2π‘–βˆ’1𝑁2π‘–βˆ’1.(2.3) Let Ξ©=(ℝ+⧡{0})2𝑛×ℝ+3𝑛, and denote the points in Ξ© by (𝑆,𝐼)𝑇, where 𝑆=(𝑆2,𝑆1,…,𝑆2𝑛,𝑆2π‘›βˆ’1) and 𝐼=(𝐼2,𝑅2,𝐼1,…,𝐼2𝑛,𝑅2𝑛,𝐼2π‘›βˆ’1). Then we rewrite the system (2.1a)–(2.1e) in compact form 𝑑𝑆𝑑𝑑=Ξ¨1(𝑆,𝐼),𝑑𝐼𝑑𝑑=Ξ¨2(𝑆,𝐼).(2.4) For any initial condition (𝑆(0),𝐼(0)) in Ξ©, system (2.1a)–(2.1e) has a unique globally defined solution (𝑆(𝑑),𝐼(𝑑)) which remains in Ξ©. Moreover, the total human population, π‘β„Ž(𝑑), and mosquitoes, 𝑁𝑣(𝑑), are bounded for all 𝑑β‰₯0. This latter result was proved in [19].

2.2. Formulation of the Objective Functional

In this section, we formulate the optimal control problem with the following functional objective (cost): 𝐽𝑒2,𝑣2,…,𝑒2𝑛,𝑣2𝑛=𝑛𝑖=1ξ‚Έξ€œπ‘‡0π‘’βˆ’π‘Ÿπ‘‘ξ‚΅πΌ2𝑖+𝑅2𝑖+𝐴2𝑖2𝑒2𝑖2+𝐡2𝑖2𝑣2𝑖2ξ‚Άπ‘‘π‘‘βˆ’Ξ₯𝑖𝑆2𝑖.(𝑇)(2.5)βˆ‘π‘›π‘–=1𝐼2𝑖 and βˆ‘π‘›π‘–=1𝑅2𝑖 are the number of infected and semi-immune of 𝑛 patches, respectively. The term (𝐴2𝑖/2)(𝑒2𝑖)2+(𝐡2𝑖/2)(𝑣2𝑖)2 is the cost of prevention and treatment where 𝐴2𝑖,𝐡2𝑖>0 are the weight factor in the cost of control. Ξ₯𝑖(𝑆2𝑖(𝑇)) is the fitness of the susceptibles at the end of the process as a result of the treatment and prevention efforts for the patch 𝑖=1,…,𝑛. We also take the same form of the Ξ₯𝑖(𝑆2𝑖(𝑇))=π‘Šπ‘†2𝑖𝑆2𝑖(𝑇), π‘Šπ‘†2𝑖β‰₯0 as in [18]. π‘Ÿ is the discount rate. The discount rate is included to allow for long-term changes, thus giving greater emphasis to control in the short rather than the long term [21]. In the above formulation, one can note that the time 𝑑=0 is the time when treatment and prevention are initiated, and the time 𝑑=𝑇 is the time when treatment and prevention are stopped.

Additionally to the above assumptions, we assume that finance for treatment and prevention is not transferable through time, so that money which is not spent immediately cannot be saved for the future purchase of treatment and prevention.

Basically, we assume that the costs are proportional to the square of the corresponding control function due to some mathematics properties (positivity, convexity…).

2.3. Existence of an Optimal Control

The basic framework of this section is to characterize the optimal control and to prove the existence and uniqueness of the optimal control. We begin to simplify the writing by noting (𝑒2,…,𝑒2𝑛,𝑣2,…,𝑣2𝑛)=(𝑒,𝑣) and (π‘’βˆ—2,…,π‘’βˆ—2𝑛,π‘£βˆ—2,…,π‘£βˆ—2𝑛)=(π‘’βˆ—,π‘£βˆ—). Because the model is linear with respect to the control variables and bounded by a linear system with respect to the state variables, the conditions for the existence of an optimal control are satisfied. While applying the Fleming and Rishel theorem [22], the existence of the 2𝑛-upplet optimal control can be obtained in our case.

Given 𝒰=(𝑒,𝑣),𝑒,𝑣,measurable0β‰€π‘Ž2𝑖≀𝑒2𝑖≀𝑏2𝑖≀1,0≀𝑐2𝑖≀𝑣2𝑖≀𝑑2𝑖≀1,(2.6) therefore, one can state the following theorem.

Theorem 2.1. Given the objective functional 𝐽(𝑒,𝑣) defined by 𝐽(𝑒,𝑣)∢=𝑛𝑖=1ξ‚Έξ€œπ‘‡0π‘’βˆ’π‘Ÿπ‘‘ξ‚΅πΌ2𝑖+𝑅2𝑖+𝐴2𝑖2𝑒2𝑖2+𝐡2𝑖2𝑣2𝑖2ξ‚Άπ‘‘π‘‘βˆ’Ξ₯𝑖𝑆2𝑖,(𝑇)(2.7) for all π‘‘βˆˆ[0;𝑇] subject to the equations of system (2.1a)–(2.1e) with 𝑆2𝑖(0)>0,𝑆2π‘–βˆ’1(0)>0,𝑅2𝑖(0)β‰₯0,𝐼2𝑖(0)β‰₯0, and 𝐼2π‘–βˆ’1(0)β‰₯0 for 𝑖=1,…,𝑛, then there exists 2𝑛-upplet optimal control (π‘’βˆ—,π‘£βˆ—) such that π½ξ€·π‘’βˆ—,π‘£βˆ—ξ€Έ=min(𝑒,𝑣)βˆˆπ’°π½(𝑒,𝑣),(2.8) when the following conditions are satisfied:(i)the class of all initial conditions with the 2𝑛-upplet optimal control in the admissible control set and corresponding state variables is nonempty,(ii)the admissible control set 𝒰 is convex and closed,(iii)the right-hand side of the state system is bounded by a linear function in the state and control,(iv)the integrand of the objective functional is convex on 𝒰 and is bounded below by βˆ‘π‘›π‘–=1(β„Ž2𝑖(|𝑒2𝑖|2+|𝑣2𝑖|2)𝜚2𝑖/2βˆ’π‘˜2𝑖), where β„Ž2𝑖,π‘˜2𝑖>0, and 𝜚2𝑖>1,(v)the function βˆ‘π‘›π‘–=1Ξ₯𝑖(𝑆2𝑖(𝑇)) is continuous with respect to the variable 𝑆2𝑖.

Proof. (i)It is obtained by definition.(ii)By definition, the admissible control set 𝒰 is convex and closed.(iii)The right-hand side of the state system (2.1a)–(2.1e) is bounded by a linear function in the state (refer to Theorem  1 of [19]). Our state system is bilinear in the control variable.(iv)To show that the integrand of the objective functional is convex on 𝒰, we must prove that 𝐹2𝑖𝑑,𝐼2𝑖,𝑅2𝑖,2𝑗=1πœ‚2𝑗𝑋2𝑗ξƒͺ≀2𝑗=1πœ‚2𝑗𝐹2𝑖𝑑,𝐼2𝑖,𝑅2𝑖,𝑋2𝑗,(2.9) where βˆ‘2𝑗=1πœ‚2𝑗=1,𝑋2𝑗=(𝑒2𝑗,𝑣2𝑗) and 𝐽(𝑒,𝑣)+𝑛𝑖=1𝑆2𝑖(𝑇)=𝑛𝑖=1ξ€œπ‘‡0π‘’βˆ’π‘Ÿπ‘‘ξ‚΅πΌ2𝑖+𝑅2𝑖+𝐴2𝑖2𝑒2𝑖2+𝐡2𝑖2𝑣2𝑖2ξ‚Ά=𝑑𝑑𝑛𝑖=1ξ€œπ‘‡0𝐹2𝑖𝑑,𝐼2𝑖,𝑅2𝑖,𝑋2𝑖,(2.10) where 𝐹2𝑖𝑑,𝐼2𝑖,𝑅2𝑖,𝑋2𝑖=π‘’βˆ’π‘Ÿπ‘‘ξ‚΅πΌ2𝑖+𝑅2𝑖+𝐴𝑖2𝑒2𝑖2+𝐡2𝑖2𝑣2𝑖2𝐹,(2.11)2𝑖𝑑,𝐼2𝑖,𝑅2𝑖,2𝑗=1πœ‚2𝑗𝑋2𝑗ξƒͺ=𝐼2𝑖+𝑅2i+𝐴2𝑖22𝑗=1πœ‚2𝑗𝑒2𝑗ξƒͺ2+𝐡2𝑖22𝑗=1πœ‚2𝑗𝑣2𝑗ξƒͺ2≀𝐼2𝑖+𝑅2𝑖+𝐴2𝑖22𝑗=1πœ‚2𝑗𝑒2𝑗2+𝐡2𝑖22𝑗=1πœ‚2𝑗𝑣2𝑗2≀2𝑗=1πœ‚2𝑗𝐼2𝑖+𝑅2𝑖+𝐴2𝑖2𝑒2𝑗2+𝐡2𝑖2𝑣2𝑗2ξ‚Ά=2𝑗=1πœ‚2𝑗𝐹2𝑖𝑑,𝐼2𝑖,𝑅2𝑖,𝑋2𝑗.(2.12) Since the sum of convex functions in the domain convex is convex, then there exists β„Ž2𝑖,π‘˜2𝑖,𝜚2𝑖>1 satisfying π‘’βˆ’π‘Ÿπ‘‘ξ‚΅πΌ2𝑖+𝑅2𝑖+𝐴2𝑖2𝑒2𝑖2+𝐡2𝑖2𝑣2𝑖2ξ‚Άβ‰₯ξ‚΅β„Ž2𝑖||𝑒2𝑖||2+||𝑣2𝑖||2ξ‚πœš2𝑖/2βˆ’π‘˜2𝑖,(2.13) because the state variable is bounded. So summing member to member, one obtains the result.(v)The function Ξ₯𝑖(𝑆2𝑖(𝑇)) is continuous so that βˆ‘π‘›π‘–=1Ξ₯𝑖(𝑆2𝑖(𝑇)) is also continuous.

2.4. Characterization of the 2𝑛-Upplet Optimal Control

Since there exists 2𝑛-upplet optimal control for minimizing the functional, (2.7), subject to system (2.1a)–(2.1e), we derive the necessary conditions on the optimal control. We discuss the theorem that relates to the characterization of the optimal control. In order to derive the necessary conditions for this optimal control, we use Pontryagin’s maximum principle [23]. The Lagrangian, sometimes called the Hamiltonian, augmented with penalty terms for control constraints is defined as 𝐿=𝑛𝑖=1𝐼2𝑖+𝑅2𝑖+𝐴2𝑖2𝑒2𝑖2+𝐡2𝑖2𝑣2𝑖2ξ‚Ά+𝑛𝑖=1πœ†π‘†2𝑖Λ2𝑖+𝛽2𝑖𝑅2𝑖+𝜌2𝑖𝐼2π‘–βˆ’ξ€·πœ‡2𝑖+Ξ¦2𝑖𝑆2𝑖+πœƒπœ†2𝑖𝑣2𝑖𝐼2𝑖+𝑛𝑗=1ξ€·π‘šπ‘†π‘–π‘—π‘†2π‘—βˆ’π‘šπ‘†π‘—π‘–π‘†2𝑖ξƒͺ+𝑛𝑖=1πœ†πΌ2𝑖Φ2𝑖𝑆2π‘–βˆ’πœ–2𝑖𝐼2π‘–βˆ’πœ†2𝑖𝑣2𝑖𝐼2𝑖+𝑛𝑗=1ξ€·π‘šπΌπ‘–π‘—πΌ2π‘—βˆ’π‘šπΌπ‘—π‘–πΌ2𝑖ξƒͺ+𝑛𝑖=1πœ†π‘…2𝑖𝛼2𝑖𝐼2π‘–βˆ’π›Ώ2𝑖𝑅2𝑖+(1βˆ’πœƒ)πœ†2𝑖𝑣2𝑖𝐼2𝑖+𝑛𝑗=1ξ€·π‘šπ‘…π‘–π‘—π‘…2π‘—βˆ’π‘šπ‘…π‘—π‘–π‘…2𝑖ξƒͺ+𝑛𝑖=1πœ†π‘†2π‘–βˆ’1ξ€·Ξ›2π‘–βˆ’1βˆ’πœ‡2π‘–βˆ’1𝑆2π‘–βˆ’1βˆ’Ξ¦2π‘–βˆ’1𝑆2π‘–βˆ’1ξ€Έ+𝑛𝑖=1πœ†πΌ2π‘–βˆ’1ξ€·Ξ¦2π‘–βˆ’1𝑆2π‘–βˆ’1βˆ’πœ‡2π‘–βˆ’1𝐼2π‘–βˆ’1ξ€Έβˆ’π‘›ξ“π‘–=1πœ”2𝑖𝑒2π‘–βˆ’π‘Ž2π‘–ξ€Έβˆ’π‘›ξ“π‘–=1πœ›2𝑖𝑏2π‘–βˆ’π‘’2π‘–ξ€Έβˆ’π‘›ξ“π‘–=1𝜁2𝑖𝑣2π‘–βˆ’π‘2π‘–ξ€Έβˆ’π‘›ξ“π‘–=1πœ‰2𝑖𝑑2π‘–βˆ’π‘£2𝑖,(2.14) where πœ†πœ‹, πœ‹=𝑆2𝑖, 𝐼2𝑖, 𝑅2𝑖, 𝑆2π‘–βˆ’1, 𝐼2π‘–βˆ’1 is the costate variable to the state variable (𝑆,𝐼), respectively, for patch 𝑖=1,…,𝑛. We can interpret πœ†πœ‹(𝑑) as the marginal value or shadow price of the last unit of 𝑆2𝑖, 𝑆2π‘–βˆ’1, 𝐼2𝑖, 𝐼2π‘–βˆ’1, and 𝑅2𝑖 was evaluated at time 𝑑 [3]. For example, πœ†π‘†2𝑖 is the increase in welfare if the number of susceptible is exogenously increased at time 𝑑. πœ†πœ‹ can be negative. The parameters πœ”2𝑖, πœ›2𝑖, 𝜁2𝑖, πœ‰2𝑖 with 𝑖=1,…,𝑛 are the penalty multipliers satisfying these conditions:πœ”2𝑖β‰₯0,𝑒2π‘–βˆ’π‘Ž2𝑖β‰₯0,πœ”2𝑖𝑒2π‘–βˆ’π‘Ž2π‘–ξ€Έπœ›=0,𝑖=1,…,𝑛,2𝑖β‰₯0,𝑏2π‘–βˆ’π‘’2𝑖β‰₯0,πœ›2𝑖𝑏2π‘–βˆ’π‘’2π‘–ξ€Έπœ=0,𝑖=1,…,𝑛,2𝑖β‰₯0,𝑣2π‘–βˆ’π‘2𝑖β‰₯0,𝜁2𝑖𝑣2π‘–βˆ’π‘2π‘–ξ€Έπœ‰=0,𝑖=1,…,𝑛,2𝑖β‰₯0,π‘‘π‘–π‘—βˆ’π‘£π‘–π‘—β‰₯0,πœ‰2𝑖𝑑2π‘–βˆ’π‘£2𝑖=0,𝑖=1,…,𝑛.(2.15) The supplementary condition at the first and second line of the system (2.15) realized at optimal control π‘’βˆ—2𝑖 and the last two lines of this system is realized at the optimal control π‘£βˆ—2𝑖.

Theorem 2.2. Given 2𝑛-upplet optimal controls (π‘’βˆ—,π‘£βˆ—) and solutions (𝑆,𝐼) of the corresponding state system (2.1a)–(2.1e), there exists adjoint variables πœ†πœ‹, with πœ‹=𝑆2𝑖, 𝑆2π‘–βˆ’1, 𝑅2𝑖, 𝐼2𝑖, 𝐼2π‘–βˆ’1 where 𝑖=1,…,𝑛 satisfying the following canonical equations: π‘‘πœ†π‘†2𝑖𝑑𝑑=π‘Ÿπœ†π‘†2π‘–βˆ’πœ•πΏπœ•π‘†2𝑖=π‘Ÿπœ†π‘†2𝑖+πœ†π‘†2π‘–ξƒ©πœ‡2𝑖+Ξ¦2𝑖1βˆ’Μƒπ‘Ž2𝑖𝑆2π‘–Μƒπ‘Ž2π‘–βˆ’1𝑁2π‘–βˆ’1+Μƒπ‘Ž2𝑖𝑁2𝑖+𝑛𝑗=1π‘šπ‘†π‘—π‘–ξƒͺβˆ’πœ†πΌ2𝑖Φ2𝑖1βˆ’Μƒπ‘Ž2𝑖𝑆2π‘–Μƒπ‘Ž2π‘–βˆ’1𝑁2π‘–βˆ’1+Μƒπ‘Ž2𝑖𝑁2𝑖+ξ€·πœ†πΌ2π‘–βˆ’1βˆ’πœ†π‘†2π‘–βˆ’1ξ€ΈΜƒπ‘Ž2𝑖𝑆2π‘–βˆ’1Ξ¦2π‘–βˆ’1Μƒπ‘Ž2π‘–βˆ’1𝑁2π‘–βˆ’1+Μƒπ‘Ž2𝑖𝑁2𝑖,π‘‘πœ†πΌ2𝑖𝑑𝑑=π‘Ÿπœ†πΌ2π‘–βˆ’πœ•πΏπœ•πΌ2𝑖=π‘Ÿπœ†πΌ2π‘–βˆ’1βˆ’πœ†π‘†2π‘–ξ‚΅πœŒ2𝑖+Μƒπ‘Ž2𝑖𝑆2𝑖Φ2π‘–Μƒπ‘Ž2π‘–βˆ’1𝑁2π‘–βˆ’1+Μƒπ‘Ž2𝑖𝑁2𝑖+πœƒπœ†2𝑖𝑣2𝑖+πœ†πΌ2π‘–ξƒ©Μƒπ‘Ž2𝑖𝑆2𝑖Φ2π‘–Μƒπ‘Ž2π‘–βˆ’1𝑁2π‘–βˆ’1+Μƒπ‘Ž2𝑖𝑁2𝑖+πœ–2𝑖+πœ†2𝑖𝑣2𝑖+𝑛𝑗=1π‘šπΌπ‘—π‘–ξƒͺβˆ’πœ†π‘…2𝑖𝛼2𝑖+(1βˆ’πœƒ)𝛼2𝑖𝑣2𝑖+ξ€·πœ†π‘†2π‘–βˆ’1βˆ’πœ†πΌ2π‘–βˆ’1ξ€ΈΜƒπ‘Ž2π‘–ξ€·Μƒπ‘Ž2π‘–βˆ’1ξ€·1βˆ’π‘’2π‘–ξ€ΈπœŽ2𝑖,2π‘–βˆ’1βˆ’Ξ¦2π‘–βˆ’1𝑆2π‘–βˆ’1Μƒπ‘Ž2π‘–βˆ’1𝑁2π‘–βˆ’1+Μƒπ‘Ž2𝑖𝑁2𝑖,π‘‘πœ†π‘…2𝑖𝑑𝑑=π‘Ÿπœ†π‘…2π‘–βˆ’πœ•πΏπœ•π‘…2𝑖=π‘Ÿπœ†π‘…2π‘–βˆ’1βˆ’πœ†π‘†2𝑖𝛽2𝑖+Μƒπ‘Ž2𝑖𝑆2𝑖Φ2π‘–Μƒπ‘Ž2π‘–βˆ’1𝑁2π‘–βˆ’1+Μƒπ‘Ž2𝑖𝑁2𝑖+πœ†πΌ2π‘–Μƒπ‘Ž2𝑖𝑆2𝑖Φ2π‘–Μƒπ‘Ž2π‘–βˆ’1𝑁2π‘–βˆ’1+Μƒπ‘Ž2𝑖𝑁2𝑖+πœ†π‘…2𝑖𝛿2𝑖+𝑛𝑗=1π‘šπ‘…π‘—π‘–ξƒͺ+ξ€·πœ†π‘†2π‘–βˆ’1βˆ’πœ†πΌ2π‘–βˆ’1ξ€ΈΜƒπ‘Ž2π‘–ξ€·Μƒπ‘Ž2π‘–βˆ’1ξ€·1βˆ’π‘’2π‘–ξ€ΈξπœŽ2𝑖,2π‘–βˆ’1βˆ’Ξ¦2π‘–βˆ’1𝑆2π‘–βˆ’1Μƒπ‘Ž2π‘–βˆ’1𝑁2π‘–βˆ’1+Μƒπ‘Ž2𝑖𝑁2𝑖,π‘‘πœ†π‘†2π‘–βˆ’1𝑑𝑑=π‘Ÿπœ†π‘†2π‘–βˆ’1βˆ’πœ•πΏπœ•π‘†2π‘–βˆ’1=π‘Ÿπœ†π‘†2π‘–βˆ’1+ξ€·πœ†πΌ2π‘–βˆ’πœ†π‘†2π‘–ξ€ΈΜƒπ‘Ž2π‘–βˆ’1𝑆2𝑖Φ2π‘–Μƒπ‘Ž2π‘–βˆ’1𝑁2π‘–βˆ’1+Μƒπ‘Ž2𝑖𝑁2𝑖+πœ†π‘†2π‘–βˆ’1πœ‡2π‘–βˆ’1+ξ€·πœ†πΌ2π‘–βˆ’1βˆ’πœ†π‘†2π‘–βˆ’1ξ€ΈΞ¦2π‘–βˆ’1ξ‚΅Μƒπ‘Ž2π‘–βˆ’1𝑆2π‘–βˆ’1Μƒπ‘Ž2π‘–βˆ’1𝑁2π‘–βˆ’1+Μƒπ‘Ž2𝑖𝑁2𝑖,βˆ’1π‘‘πœ†πΌ2π‘–βˆ’1𝑑𝑑=π‘Ÿπœ†πΌ2π‘–βˆ’1βˆ’πœ•πΏπœ•πΌ2π‘–βˆ’1=π‘Ÿπœ†πΌ2π‘–βˆ’1+ξ€·πœ†π‘†2π‘–βˆ’πœ†πΌ2π‘–ξ€ΈΜƒπ‘Ž2π‘–βˆ’1ξ€·Μƒπ‘Ž2𝑖1βˆ’π‘’2π‘–ξ€ΈπœŽ2π‘–βˆ’1,2π‘–βˆ’Ξ¦2𝑖𝑆2π‘–Μƒπ‘Ž2π‘–βˆ’1𝑁2π‘–βˆ’1+Μƒπ‘Ž2𝑖𝑁2𝑖+πœ†πΌ2π‘–βˆ’1πœ‡2π‘–βˆ’1βˆ’ξ€·πœ†π‘†2π‘–βˆ’1βˆ’πœ†πΌ2π‘–βˆ’1ξ€ΈΜƒπ‘Ž2π‘–βˆ’1𝑆2π‘–βˆ’1Ξ¦2π‘–βˆ’1Μƒπ‘Ž2π‘–βˆ’1𝑁2π‘–βˆ’1+Μƒπ‘Ž2𝑖𝑁2𝑖,(2.16) with the transversality conditions (terminal conditions): πœ†π‘†2𝑖(𝑇)=πœ•Ξ₯π‘–πœ•π‘†2𝑖||||𝑑=𝑇,πœ†πœ‹(𝑇)=0,𝑖=1,…,𝑛forπœ‹=𝐼2𝑖,𝐼2π‘–βˆ’1,𝑅2𝑖,𝑆2π‘–βˆ’1.(2.17) Furthermore, the following characterization of optimal control holds:π‘’βˆ—2π‘–ξƒ©π‘Ž=max2𝑖𝑏,min2𝑖,Μƒπ‘Ž2π‘–βˆ’1Μƒπ‘Ž2𝑖𝐴2π‘–ξƒ¬πœ†π‘†πΌ2π‘–πœŽ2π‘–βˆ’1,2𝑖𝐼2π‘–βˆ’1𝑆2π‘–Μƒπ‘Ž2π‘–βˆ’1𝑁2π‘–βˆ’1+Μƒπ‘Ž2𝑖𝑁2𝑖+πœ†π‘†πΌ2π‘–βˆ’1ξ€·πœŽ2𝑖,2π‘–βˆ’1𝐼2𝑖+𝜎2𝑖,2π‘–βˆ’1𝑅2𝑖𝑆2π‘–βˆ’1Μƒπ‘Ž2π‘–βˆ’1𝑁2π‘–βˆ’1+Μƒπ‘Ž2𝑖𝑁2𝑖,𝑣ξƒͺξƒͺβˆ—2𝑖𝑐=max2𝑖𝑑,min2𝑖,βˆ’πœ†2𝑖𝐼2π‘–ξ€·πœƒπœ†π‘†2π‘–βˆ’πœ†πΌ2𝑖+(1βˆ’πœƒ)πœ†π‘…2𝑖𝐡2𝑖,ξƒͺξƒͺ(2.18) where πœ†π‘†π‘—βˆ’πœ†πΌπ‘—=βˆ’πœ†π‘†πΌπ‘—.

Proof. The adjoint equations and transversality conditions are standard results from Pontryagin’s maximum principle. Also, solutions to the adjoint system exist and are bounded. To determine the interior optimum of our Lagrangian, we take the partial derivatives of Lagrangian 𝐿 with respect to π‘’βˆ—2𝑖 and π‘£βˆ—2𝑖 and set it equal to zero: πœ•πΏπœ•π‘’2𝑖=βˆ’Μƒπ‘Ž2π‘–βˆ’1Μƒπ‘Ž2π‘–πœ†π‘†πΌ2π‘–πœŽ2π‘–βˆ’1,2𝑖𝐼2π‘–βˆ’1𝑆2π‘–Μƒπ‘Ž2π‘–βˆ’1𝑁2π‘–βˆ’1+Μƒπ‘Ž2𝑖𝑁2π‘–βˆ’Μƒπ‘Ž2π‘–βˆ’1Μƒπ‘Ž2𝑖𝑆2π‘–βˆ’1πœ†π‘†πΌ2π‘–βˆ’1ξ€·πœŽ2𝑖,2π‘–βˆ’1𝐼2𝑖+𝜎2𝑖,2π‘–βˆ’1𝑅2π‘–ξ€ΈΜƒπ‘Ž2π‘–βˆ’1𝑁2π‘–βˆ’1+Μƒπ‘Ž2𝑖𝑁2𝑖+πœ›2π‘–βˆ’πœ”2𝑖+𝐴2𝑖𝑒2𝑖=0,πœ•πΏπœ•π‘£2𝑖=πœ†2𝑖𝐼2π‘–ξ€·πœƒπœ†π‘†2π‘–βˆ’πœ†πΌ2𝑖+(1βˆ’πœƒ)πœ†π‘…2𝑖+πœ‰2π‘–βˆ’πœ2𝑖+𝐡2𝑖𝑣2𝑖=0.(2.19) Solving for optimal control, we have π‘’βˆ—2𝑖=1𝐴2π‘–ξƒ¬Μƒπ‘Ž2π‘–βˆ’1Μƒπ‘Ž2π‘–πœ†π‘†πΌ2π‘–πœŽ2π‘–βˆ’1,2𝑖𝐼2π‘–βˆ’1𝑆2π‘–Μƒπ‘Ž2π‘–βˆ’1𝑁2π‘–βˆ’1+Μƒπ‘Ž2𝑖𝑁2𝑖+Μƒπ‘Ž2π‘–βˆ’1Μƒπ‘Ž2𝑖𝑆2π‘–βˆ’1πœ†π‘†πΌ2π‘–βˆ’1ξ€·πœŽ2𝑖,2π‘–βˆ’1𝐼2𝑖+𝜎2𝑖,2π‘–βˆ’1𝑅2π‘–ξ€ΈΜƒπ‘Ž2π‘–βˆ’1𝑁2π‘–βˆ’1+Μƒπ‘Ž2𝑖𝑁2π‘–βˆ’πœ›2𝑖+πœ”2𝑖,π‘£βˆ—2𝑖=1𝐡2π‘–ξ€Ίβˆ’πœ†2𝑖𝐼2π‘–ξ€·πœƒπœ†π‘†2π‘–βˆ’πœ†πΌ2𝑖+(1βˆ’πœƒ)πœ†π‘…2π‘–ξ€Έβˆ’πœ‰2𝑖+𝜁2𝑖.(2.20) To determine an explicit expression for the optimal control without the penalty multipliers πœ”2𝑖,πœ›2𝑖,𝜁2𝑖,πœ‰2𝑖, a standard optimality technique is used. We consider the following cases to discuss the control: case of the prevention or case of the treatment.(i) Case of prevention:(1)on the set ξ€½π‘‘π‘Ž2𝑖<π‘’βˆ—2𝑖<𝑏2𝑖,𝑖=1,…,𝑛,(2.21) we have πœ”2𝑖=πœ›2𝑖=0. Hence, the optimal control is π‘’βˆ—2𝑖=1𝐴2π‘–ξƒ¬Μƒπ‘Ž2π‘–βˆ’1Μƒπ‘Ž2π‘–πœ†π‘†πΌ2π‘–πœŽ2π‘–βˆ’1,2𝑖𝐼2π‘–βˆ’1𝑆2π‘–Μƒπ‘Ž2π‘–βˆ’1𝑁2π‘–βˆ’1+Μƒπ‘Ž2𝑖𝑁2𝑖+Μƒπ‘Ž2π‘–βˆ’1Μƒπ‘Ž2𝑖𝑆2π‘–βˆ’1πœ†π‘†πΌ2π‘–βˆ’1ξ€·πœŽ2𝑖,2π‘–βˆ’1𝐼2𝑖+𝜎2𝑖,2π‘–βˆ’1𝑅2π‘–ξ€ΈΜƒπ‘Ž2π‘–βˆ’1𝑁2π‘–βˆ’1+Μƒπ‘Ž2𝑖𝑁2𝑖,(2.22)(2)on the set ξ€½π‘‘π‘’βˆ—2𝑖=𝑏2𝑖,𝑖=1,…,𝑛,(2.23) we have πœ”2𝑖(𝑑)=0. Hence, 𝑏2𝑖=π‘’βˆ—2𝑖=1𝐴2π‘–ξƒ¬Μƒπ‘Ž2π‘–βˆ’1Μƒπ‘Ž2π‘–πœ†π‘†πΌ2π‘–πœŽ2π‘–βˆ’1,2𝑖𝐼2π‘–βˆ’1𝑆2π‘–Μƒπ‘Ž2π‘–βˆ’1𝑁2π‘–βˆ’1+Μƒπ‘Ž2𝑖𝑁2𝑖+Μƒπ‘Ž2π‘–βˆ’1Μƒπ‘Ž2𝑖𝑆2π‘–βˆ’1πœ†π‘†πΌ2π‘–βˆ’1ξ€·πœŽ2𝑖,2π‘–βˆ’1𝐼2𝑖+𝜎2𝑖,2π‘–βˆ’1𝑅2π‘–ξ€ΈΜƒπ‘Ž2π‘–βˆ’1𝑁2π‘–βˆ’1+Μƒπ‘Ž2𝑖𝑁2π‘–βˆ’πœ›2𝑖.(2.24) This implies that 1𝐴2π‘–ξƒ¬Μƒπ‘Ž2π‘–βˆ’1Μƒπ‘Ž2π‘–πœ†π‘†πΌ2π‘–πœŽ2π‘–βˆ’1,2𝑖𝐼2π‘–βˆ’1𝑆2π‘–Μƒπ‘Ž2π‘–βˆ’1𝑁2π‘–βˆ’1+Μƒπ‘Ž2𝑖𝑁2𝑖+Μƒπ‘Ž2π‘–βˆ’1Μƒπ‘Ž2𝑖𝑆2π‘–βˆ’1πœ†π‘†πΌ2π‘–βˆ’1ξ€·πœŽ2𝑖,2π‘–βˆ’1𝐼2𝑖+𝜎2𝑖,2π‘–βˆ’1𝑅2π‘–ξ€ΈΜƒπ‘Ž2π‘–βˆ’1𝑁2π‘–βˆ’1+Μƒπ‘Ž2𝑖𝑁2𝑖β‰₯𝑏2𝑖,(2.25) since πœ›2𝑖(𝑑)β‰₯0,(3)on the set ξ€½π‘‘π‘’βˆ—2𝑖=π‘Ž2𝑖,𝑖=1,…,𝑛,(2.26) we have πœ›2𝑖(𝑑)=0. Hence, π‘’βˆ—2𝑖=1𝐴2π‘–ξƒ¬Μƒπ‘Ž2π‘–βˆ’1Μƒπ‘Ž2π‘–πœ†π‘†πΌ2π‘–πœŽ2π‘–βˆ’1,2𝑖𝐼2π‘–βˆ’1𝑆2π‘–Μƒπ‘Ž2π‘–βˆ’1𝑁2π‘–βˆ’1+Μƒπ‘Ž2𝑖𝑁2𝑖+Μƒπ‘Ž2π‘–βˆ’1Μƒπ‘Ž2𝑖𝑆2π‘–βˆ’1πœ†π‘†πΌ2π‘–βˆ’1ξ€·πœŽ2𝑖,2π‘–βˆ’1𝐼2𝑖+𝜎2𝑖,2π‘–βˆ’1𝑅2π‘–ξ€ΈΜƒπ‘Ž2π‘–βˆ’1𝑁2π‘–βˆ’1+Μƒπ‘Ž2𝑖𝑁2𝑖+πœ”2𝑖.(2.27) This implies that 1𝐴2π‘–ξƒ¬Μƒπ‘Ž2π‘–βˆ’1Μƒπ‘Ž2π‘–πœ†π‘†πΌ2π‘–πœŽ2π‘–βˆ’1,2𝑖𝐼2π‘–βˆ’1𝑆2π‘–Μƒπ‘Ž2π‘–βˆ’1𝑁2π‘–βˆ’1+Μƒπ‘Ž2𝑖𝑁2𝑖+Μƒπ‘Ž2π‘–βˆ’1Μƒπ‘Ž2𝑖𝑆2π‘–βˆ’1πœ†π‘†πΌ2π‘–βˆ’1ξ€·πœŽ2𝑖,2π‘–βˆ’1𝐼2𝑖+𝜎2𝑖,2π‘–βˆ’1𝑅2π‘–ξ€ΈΜƒπ‘Ž2π‘–βˆ’1𝑁2π‘–βˆ’1+Μƒπ‘Ž2𝑖𝑁2π‘–ξƒ­β‰€π‘Ž2𝑖.(2.28) Combining these cases, the optimal control π‘’βˆ—2𝑖 for 𝑖=1,…,𝑛 is characterized as π‘’βˆ—2π‘–ξƒ©π‘Ž=max2𝑖𝑏,min2𝑖,Μƒπ‘Ž2π‘–βˆ’1Μƒπ‘Ž2𝑖𝐴2π‘–ξƒ¬πœ†π‘†πΌ2π‘–πœŽ2π‘–βˆ’1,2𝑖𝐼2π‘–βˆ’1𝑆2π‘–Μƒπ‘Ž2π‘–βˆ’1𝑁2π‘–βˆ’1+Μƒπ‘Ž2𝑖𝑁2𝑖+πœ†π‘†πΌ2π‘–βˆ’1ξ€·πœŽ2𝑖,2π‘–βˆ’1𝐼2𝑖+𝜎2𝑖,2π‘–βˆ’1𝑅2𝑖𝑆2π‘–βˆ’1Μƒπ‘Ž2π‘–βˆ’1𝑁2π‘–βˆ’1+Μƒπ‘Ž2𝑖𝑁2𝑖.ξƒ­ξƒͺξƒͺ(2.29)(ii) Case of treatment: using similar arguments as in the case of prevention, we also obtain the second optimal π‘£βˆ—2𝑖 with 𝑖=1,…,𝑛 control function is characterized by π‘£βˆ—2𝑖𝑐=max2𝑖𝑑,min2𝑖,βˆ’πœ†2𝑖𝐼2π‘–ξ€·πœƒπœ†π‘†2π‘–βˆ’πœ†πΌ2𝑖+(1βˆ’πœƒ)πœ†π‘…2𝑖𝐡2𝑖ξƒͺξƒͺ.(2.30)

3. Numerical Results and Discussion

3.1. Parameters

We fix the probability for treatment of infectious humans for patches 2 and 3 at πœƒ=0.5. Also we take πœ†2𝑗=0.5, the weights of prevention and treatment 𝐴2𝑗=𝐡2𝑗=50, and the bounds of all control π‘Ž2𝑗=𝑐2𝑗=0,𝑏2𝑗=𝑑2𝑗=1 with 𝑗=2,3. We fix the coefficient of fitness π‘Šπ‘†2𝑖=1.

We take a very small discount rate π‘Ÿ=0.0001 because the daily discounting of the cost decreases very slowly. The other parameters of the model were obtained from [1] as well as the following value of the migration matrix: data for migration matrix for the semi-immune, 𝑀𝑅, 𝑀𝑅=⎑⎒⎒⎒⎒⎣00.7Γ—10βˆ’10.8Γ—10βˆ’10.1Γ—10βˆ’100.1Γ—10βˆ’40.2Γ—10βˆ’10.1Γ—10βˆ’40⎀βŽ₯βŽ₯βŽ₯βŽ₯⎦,(3.1)𝑀𝑆, data for migration matrix for susceptible, and 𝑀𝐼, data for migration matrix for the infectious 𝑀𝑆=𝑀𝐼=⎑⎒⎒⎒⎒⎣00.7Γ—10βˆ’30.8Γ—10βˆ’30.1Γ—10βˆ’300.1Γ—10βˆ’60.2Γ—10βˆ’30.1Γ—10βˆ’60⎀βŽ₯βŽ₯βŽ₯βŽ₯⎦.(3.2)

3.2. Implementation

To solve our problem of optimal control, we used the program MATLAB dynamic optimisation code (DYNOPT), which is a set of MATLAB functions for the determination of optimal control trajectory by describing the process, the cost to be minimized, subject to equality and inequality constraints, and using orthogonal collocation on the finite elements method [24]. For more information about this algorithm, we can see the user’s guide in [24].

We implemented the model with the initial condition: 𝑆2(0)=15000;  𝑆4(0)=50;𝑆6(0)=1000;  𝐼2(0)=1000;  𝐼4(0)=50;  𝐼6(0)=100;  𝑅2(0)=100;  𝑅4(0)=250;  𝑅6(0)=10;  𝑆1(0)=5000;  𝑆3(0)=8000;  𝑆5(0)=5000;  𝐼1(0)=50;  𝐼3(0)=6000;  𝐼5(0)=4000.

The weight assigned to the controls is much higher than the weight assigned to the state variables because the two functions are not expressed in the same scale. The controls are expressed in terms of cost, while infections and semi-immune are expressed in term of number of individuals. We chose the time at 𝑇=300 days for our simulation.

3.3. Results and Discussions
3.3.1. Basic Reproductive Number and Reservoir of Infection

The basic reproduction number generally denoted by β„›0 is the expected number of secondary cases produced by a typical infective individual introduced into a completely susceptible population, in the absence of any control measure [25, 26]. Using the data on Table 2 which were compiled in [1] without control variable, β„›0 was equal to 3.864. Therefore, there is a persistence of the disease in the whole population (patches 1, 2, and 3). In [1], it was shown that only patches 2 and 3 constitute a reservoir of infection. Indeed, a subgroup of patches is said to be a reservoir when only targeting a control on the reservoir is sufficient to eliminate the malaria in the whole population (all the three patches). As such the patch 1 cannot sustain an epidemic by itself.

Table 2 
Value compiled in [1]: patches 2 and 3 correspond to rural areas, while patch 1 corresponds to urban area.
3.3.2. Evolution over Time of the Optimal Control in the Three Patches

We seek the optimal solution by minimizing the number of infectious hosts and semi-immune, in all patches by considering four cases: the first case where we seek the optimal solution when we consider simultaneously prevention and treatment in the two patches (see Figure 2(a)), the second case where we seek the optimal solution only with prevention without treatment in two patches (see Figure 2(b)), the third case where we seek the optimal solution only with the treatment without prevention in two patches (see Figure 2(c)), and finally the fourth case where no strategy of prevention and treatment is applied. Figure 2 shows a strong preventive action early in the process of elimination of the disease and a high processing action at the end of the process. Between these two strategies, prevention and treatment are preferred to reduce the number of infections and semi-immunes in all patches. Interestingly, these results show that the dynamic of controls depends on the bounds that we choose for the controls.

(a)
(a)
(b)
(b)
(c)
(c)
Figure 2 
Results of the simulations achieved using data from Table 2 showing the evolution over time of the optimal control in the two reservoirs of infection. (a) Optimal control for prevention and treatment; (b) Optimal control for prevention control with no treatment control; (c) optimal control for treatment control with no prevention control.
3.3.3. Dynamics of Human Infection in the Three Patches

Figures 3 and 4 show that the increase of susceptible hosts involves a decrease of infectious hosts. Figure 4(a) shows that no action should be taken during the half time interval in a patch which is our urban area. The second half time should be considered the treatment of infections from two other patches. This is because it takes time (𝑇/2) for production of sick people from the reservoir area of infection, and after this time (𝑇/2), we realize that all the patches contain enough sick people. However, we must now apply a treatment in the area that is not a reservoir of infection initially. This treatment will be done by setting up at the entrances to the urban area by the distribution of drugs to fight malaria infection before accessing it. These measures of treatment become necessary to prevent the urban area, constitutes a reservoir of infection.

(a)
(a)
(b)
(b)
(c)
(c)
Figure 3 
Results of the simulations achieved using data from Table 2 showing the evolution over time of susceptible hosts for the three patches. We show the four cases: black line optimal solution solved with prevention and treatment in the 2 patches; blue line optimal solution solved with prevention control in the 2 patches; red line optimal solution solved with treatment control in the 2 patches; green line solution without control. (a) Evolution of the susceptible hosts for patch 1; (b) evolution of the susceptible hosts for patch 2; (c) evolution of the susceptible hosts for patch 3.
(a)
(a)
(b)
(b)
(c)
(c)
Figure 4 
Results of the simulations achieved using data from Table 2 showing the evolution over time of infectious host for the three patches. We show the four cases: black line optimal solution solved with prevention and treatment in the 2 patches; blue line optimal solution solved with prevention control in the 2 patches; red line optimal solution solved with treatment control in the 2 patches; green line solution without control. (a) Evolution of the infectious hosts for patch 1; (b) evolution of the infectious hosts for patch 2; (c) evolution of the infectious hosts for patch 3.

Figure 4(b) shows that the treatment is effective for the infectious 50(𝑇/6) first days because patch 2 is a reservoir of infection. Between days 50(𝑇/6) and 250(5𝑇/6), we must apply simultaneously prevention and treatment, and after 250(5𝑇/6) days, only prevention can be applied in this patch.

Figure 4(c) shows that during the first 200(2𝑇/3) days we must apply simultaneously the prevention and the treatment, and after this time, only treatment should be applied to reduce the number of infections in patch 3.

3.3.4. Dynamics of Semi-Immune in the Three Patches

Figure 5(a) shows that no strategies must be applied during the 𝑇/4 first days for the semi-immune in patch 1. Between 𝑇/4 and 3𝑇/4 days, we must apply simultaneously prevention and treatment, and during the remaining period, only prevention should be applied.

(a)
(a)
(b)
(b)
(c)
(c)
Figure 5 
Results of the simulations achieved using data from Table 2 showing the evolution over time of semi-immune host for the three patches. We show the four cases: black line optimal solution solved with prevention and treatment in the 2 patches; blue line optimal solution solved with prevention control in the 2 patches; red line optimal solution solved with treatment control in the 2 patches; green line solution without control. (a) Evolution of the semi-immune hosts for patch 1; (b) evolution of the semi-immune hosts for patch 2; (c) evolution of the semi-immune hosts for patch 3.

Figures 5(b) and 5(c) show that the same strategies of control should be considered during the same period in patch 1 to reduce, respectively, the number of semi-immunes in patches 2 and 3.

To summarize, we used a recent technique of identification of the spatial infection of connected patches, to design a strategy based on the status of infection of the reservoirs of infection. We show that it is better to treat people only in areas that do not constitute a reservoir of infection and use simultaneously the prevention and the treatment to reduce the number of infections in all patches constituting a reservoir of infection. While reducing the number of semi-immunes, no differences in control strategies is made based on the type of infection reservoir. Whatever the level of infection of the reservoir of infection, the strategy to reduce the number of semi-immunes remains the same: no strategy is adopted in the early stage of malaria control, then both treatment and prevention are implemented, and in the last period, only prevention is implemented.

4. Conclusion

A mathematical model has been developed for malaria using the theory of optimal control. The formulation of the optimal control includes 𝑛 control variables for prevention and 𝑛 variables for treatment. The mathematical analysis proved the existence of an optimal control for 𝑛 connected patches under suitable conditions using the Fleming and Rishel theorem. Furthermore, using Pontryagin’s maximum principle, a characterization of the optimal control was given. Numerical simulations were also performed showing the evolution over time of the optimal control as well as the different health status of humans and mosquitoes within each patch. These results underline the usefulness of a synergy control rather than only the prevention or the treatment. The results of our simulation show that we must choose a strategy based on the infectious status of the reservoir of infection. We show that it is better to treat people only in areas that do not constitute a reservoir of infection and use simultaneously the prevention and the treatment to reduce the number of infections in all patches constituting a reservoir of infection. While reducing the number of semi-immunes, no differences in control strategies is made based on the type of infection reservoir. Whatever the level of infection of the reservoir of infection, the strategy to reduce the number of semi-immunes remains the same: no strategy is adopted in the early stage of malaria control, then both treatment and prevention are implemented, and in the last period, only prevention is implemented.

To state our main perspectives, we will include a budget constraint in our optimal problem. Before characterizing the optimal prevention and treatment, two cases may arise under budget constraints: (i) when the budget allocation for the prevention and treatment is sufficient; (ii) when the budget is insufficient. Moreover, it would be interesting to apply a sensitivity analysis for some key parameters of the model.