Electronic Journal of Differential Equations, Vol. 2020 (2020), No. 70, pp. 1–23. ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu OPTIMAL CONTROL ANALYSIS APPLIED TO A TWO-PATCH MODEL FOR GUINEA WORM DISEASE STEADY MUSHAYABASA, ANTHONY A. E. LOSIO, CHAIRAT MODNAK, JIN WANG Communicated by Suzanne M. Lenhart Abstract. We applied optimal control theory to a mathematical model for guinea worm disease, to determine the effectiveness of optimal education cam- paigns on long-term dynamics of the disease. Our model is concerned with two different host populations, represented by two patches, sharing a common water source. We computed the basic reproduction number of the model and demonstrated that whenever the reproduction number is less than unity the disease dies out in the community. Also we established that when the basic re- production number is greater than unity the disease persists. Utilizing optimal control theory, we explored the potential of time dependent education to elim- inate the disease within 120 months. The model showed that time dependent education can be successful to minimize disease prevalence in the two patches, however, its success strongly depends on the total cost of implementation as well as its maximum strength. 1. Introduction In the previous three decades the global community witnessed a major decline on the number of reported Guinea worm disease (GWD) cases, from 3.5 million in 20 countries in 1986 to only 22 cases in 2015 from only 4 countries, namely South Sudan, Mali, Chad and Ethiopia [3]. Prior studies have attributed this remarkable decline to the Guinea worm eradication program launched in the 1980s [3]. Despite this remarkable achievement, there are a couple of challenges that may thwart global eradication of the disease. Among many other challenges, continued use and sharing of open water sources by individuals from differential communities in GWD endemic countries remains a stumbling block to GWD eradication [3, 9, 24]. Individuals from different communities are known to have different behavior and knowledge about the disease [23]. Such heterogeneities have strong impacts on the dynamics of GWD, hence there is need to evaluate and quantify the possibility of GWD eradication under these circumstances. GWD is a waterborne disease which infects humans when they drink unfiltered water containing copepods (small crustaceans) which are infected with larvae of Dracunculus medinensis. Once ingested, the copepods are destroyed by gastric 2010 Mathematics Subject Classification. 92D30, 49M05. Key words and phrases. Mathematical model; Guinea worm disease; optimal control theory. c©2020 Texas State University. Submitted January 2, 2018. Published July 2, 2020. 1 2 S. MUSHAYABASA, A. A. E. LOSIO, C. MODNAK, J. WANG EJDE-2020/70 juice in the human stomach, releasing the infective larvae [3, 7]. The released larvae will then migrate to small intestines where they penetrate through the duodenal wall and become adult worms in connective tissues. After maturation into adults and copulation, the male worm subsequently dies while the female worm grow into a full-size adult (length: 70 to 120 cm) [4]. After a year long incubation period, during which time the human host shows almost no symptoms, the female worm will migrate to other parts of the body, mostly on the distal lower extremity and induces a blister on the skin which ruptures leading to the emergence of the worm [6]. The itching and pain of the induced blister prompts the use of water as therapy, whereby infected individuals in rural communities often immerse their affected body part into water sources to get relief from the pain. Subsequently, the female worm releases thousands of its immature larvae (L1) [12]. These larvae are then ingested by copepods also called cyclops. Cyclops serve as intermediate host for Dracunculus medinensis [3]. Although GWD is rarely fatal, it may cause permanent disability and may result in loss of family income and school absenteeism [3, 13]. Since GWD has no vaccine or medicine the eradication programme has been based on the following preventive measures: provision of improved water sources, use of water filtration using different types of cloth or filament filters, health edu- cation to inform populations of how the infection is acquired and can be prevented, control of the intermediate host copepod using the larvicide, Abate, (temephos), containment of cases before they have an opportunity to contaminate water sources, active surveillance in endemic or previously endemic [24, 26]. In this paper, we aim to utilize a mathematical model to provide a comprehensive assessment on the impact of time dependent educational campaigns on controlling the spread of GWD in heterogeneous environment characterised by communities sharing a common open water source. Mathematical models play an important role in understanding and providing solutions to phenomena which are difficult to measure in the field [19, 22]. Only a few mathematical models have been con- tributed by researchers on the study of GWD (see, for example [1, 17, 18, 26, 31]). For example Adetunde, [1] investigated the current pattern of GWD in the northern region of Ghana. He analyzed the data from the region and proposed a time series model for the purpose of prediction. From his analysis it was observed that the number of infection cases reduces with time, implying that if the current trend con- tinues then there is a likelihood that GWD will be completely eradicated. Recently, Link [17] proposed a compartment framework to study the dynamics of GWD. The author determined an algebraic solution to disease-free equilibrium and performed numerical stability analysis of the solution as well. Using the next generation ma- trix, Link determined the reproduction number R0 and found that the disease-free equilibrium is stable provided that the host visitation rate to the contaminated water sources is reduced. Smith? et al. [31] also developed a mathematical model for GWD. Impulsive differential equations were used to evaluate the effectiveness of chlorination. Latin Hypercube Sampling was used to determine model parame- ters that were highly sensitive to the basic reproduction number. From their work Smith? and co-workers established that education is the most effective intervention method, but a combination of education, chlorination and filtration will likely be required to achieve the final steps in the long journey to eradication. Although these studies produced many useful results and improved the exist- ing knowledge on GWD dynamics, however, none of these studies evaluated the EJDE-2020/70 TWO-PATCH GUINEA WORM DISEASE MODEL 3 impact of time dependent educational campaigns on long term GWD dynamics in an environment where different communities share a common water source. This paper aims to fill this gap. The literature on epidemic models of the concurrent spread of infectious diseases and information is quite substantial (see, for example [5, 14, 25, 29, 30]). In [14], Joshi et al. used a deterministic compartmental model to demonstrate the effects of time dependent educational campaigns on lowering susceptibility classes and its impact on long term disease dynamics. Also, Samanta et al. [29], utilized a mathematical model to illustrate the effects of awareness programs by media on epidemic outbreaks. These studies, indeed, produced useful results and we will utilize some of the results to develop and analyze our framework. The remainder of this paper is structured as follows. In Section 2, we present the methods and main results of the study. Precisely, we will analyze the dynamical behavior of the model, and perform an optimal control study for the implementation of educational campaigns. Analytical results in this section will be supported by numerical results constructed using MATLAB software. Section 3 provides a brief discussion to conclude the paper. 2. Mathematical model We developed a mathematical model for GWD that comprises of two popula- tions (patches) sharing a common water source. Although there are many possible network configurations here we considered a network of non-mixing patches with a single common water source, allowing patches to differ in the level of GWD awareness which in-turn affects their degree of risk to infection as well as their contribution to the pathogen population in the environment. This simple frame- work is motivated by living conditions in rural GWD endemic countries such as Ethiopia [3]. In rural Ethiopia, the provision of potable water supply remains a big challenge, for instance in one recent published study it was noted that only 42 out of 70 villages had access to safe drinking water [2]. Due to the unavailability of potable water supply, it was also noted that villagers in Ethiopia and other GWD endemic countries highly depend on water from open sources which are often shared by several communities [9]. According to Mari et al [23] different communities may be endowed with differential infection risk due to their geographical and/or socioe- conomic factors. Heterogeneity in health education is one of the socio-economic factors that leads to differential in infection risk among individuals living within or in different communities. To account for heterogeneity in this study we let the host population at time t in patch i by Ni(t), i = 1, 2, such that N(t) = N1(t) +N2(t). Further, we subdivided the total population of individuals in patch i into categories of: susceptible Si(t), ex- posed (latently infected) Ei(t) and infectious Ii(t), individuals. Meanwhile, let the dynamics of infected copepods in the environment be represented by compartment W (t). Note that, here the intrinsic dynamics of the parasite have been neglected, as our aim is to provide a comprehensive, yet as simple as possible, account of the impact of heterogeneity in awareness on GWD dynamics. Therefore the proposed 4 S. MUSHAYABASA, A. A. E. LOSIO, C. MODNAK, J. WANG EJDE-2020/70 model is governed by the following system of ordinary differential equations: Ṡi(t) = Λi − βi(1− ui)SiW − µiSi, Ėi(t) = βi(1− ui)SiW − (αi + µi)Ei, İi(t) = αiEi − (κi + µi)Ii, Ẇ (t) = (1− u1)γ1I1 + (1− u2)γ2I2 − εW. (2.1) All model variables and parameters are considered to be positive. Model parameter Λi denotes the constant recruitment rate of the host through birth and movement of susceptible individuals, µi is the natural mortality rate of the host; i.e., 1 µi is the average life span, βi account for indirect disease transmission rate, ui represents the effects of health education on disease transmission, αi is the incubation rate and it follows that 1 αi is the incubation period-for GWD and it ranges between 10–14 months [3], κi is the recovery rate; i.e., 1 κi is the average infectious period. Meanwhile, clinically infected individuals of patch i contribute to the parasite pop- ulation in the environment at rate γi when they immerse their affected body part into water sources to get relief from the pain. Infected copepods die naturally at rate ε. Table 1 presents a summary of the model parameters and their baseline values. Table 1. Description of model parameters and their baseline values Var. Definition Value Units Source β1 Copepods ingestion rate for patch 1 1.2 × 10−5 month−1 [18] β2 Copepods ingestion rate for patch 2 3β1 month−1 µi Host natural death rate 0.0015 month−1 [18] ε Parasite death rate 2.16 month−1 [18] αi Incubation rate 0.9996 month−1 [18] κi Recovery rate 0.0083 month−1 [18] Λi Host birth rate 100 Individ. month−1 [18] γ1 Parasite shedding rate for patch 1 0.2 month−1 [18] γ2 Parasite shedding rate for patch 2 0.2 month−1 ui Efficacy of educational campaigns in Path i varied Dimensionless Since model (2.1) represents human and pathogen population, all parameters in the model are non-negative and one can easily verify that all the solutions of the system are non-negative, given non-negative initial values. The model (2.1) will be analyzed in a biologically feasible region Ω defined as{ (N,W ) ∈ R7 + : N ≤ (Λ1 + Λ2) [min(µ1, µ2)] , W ≤ (Λ1 + Λ2)((1− u1)γ1 + (1− u2)γ2) [εmin(µ1, µ2)] } . Thus, the compact set Ω is positively invariant and attracting for model (2.1). Table 1 presents the model parameters and their baseline values. Now, we compute the basic reproduction number of the model (2.1). The basic reproduction number R0, is arguably the most important quantity in infectious disease epidemiology. It can be utilized to account for the severity of the disease in the community whether epidemic or pandemic. To compute the basic reproduction EJDE-2020/70 TWO-PATCH GUINEA WORM DISEASE MODEL 5 number of system (2.1) we need to determine the disease-free equilibrium (DFE) of the system. Thus, one can easily verify that system (2.1) has the following DFE E0 = ( S0 1 , S 0 2 , E 0 1 , E 0 2 , I 0 1 , I 0 2 ,W 0 ) = (Λ1 µ1 , Λ2 µ2 , 0, 0, 0, 0, 0 ) . Using the next-generation matrix notations in [32], the non-negative matrix F that denotes the generation of new infections and the non-singular matrix V that denotes the disease transfer among compartments, are respectively given by (at DFE) F =  0 0 0 0 β1(1− u1)S0 1 0 0 0 0 β2(1− u2)S0 2 0 0 0 0 0 0 0 0 0 0  , V =  m1 0 0 0 0 0 m2 0 0 0 −α1 0 m3 0 0 0 −α2 0 m4 0 0 0 −m5 −m6 ε  , where m1 = (α1 + µ1), m2 = (α2 + µ2), m3 = (κ1 + µ1), m4 = (κ2 + µ2), m5 = (1− u1)γ1, m6 = (1− u2)γ2. It follows that the spectral radius of system (2.1) is determined from ρ(FV −1) and is given by R0 = β1α1(1− u1)2γ1Λ1 εµ1m1m3 + β2α2(1− u2)2γ2Λ2 εµ2m2m4 . (2.2) The quantity R0 corresponds to the average number of secondary infections through environment-to-host transmission caused by one infectious individual in its infec- tious lifetime. Next, we investigate the existence of the endemic equilibrium point of model system (2.1). Theorem 2.1. Suppose R0 > 1. Then model system (2.1) admits a unique endemic equilibrium point. Proof. Let the endemic equilibrium of (2.1) be denoted by E∗(S∗i , E∗i , I∗i ,W ∗), i = 1, 2. This steady state is determined by solving the system of equations Λi − βi(1− ui)S∗iW ∗ − µiS∗i = 0, βi(1− ui)S∗iW ∗ − (αi + µi)E ∗ i = 0, αiE ∗ i − (κi + µi)I ∗ i = 0, 2∑ i=1 (1− ui)γiI∗i − εW ∗ = 0. (2.3) For the third and fourth equations of (2.3) we have Ei = (κi + µi) αi I∗i , W = 2∑ i=1 (1− ui)γi ε I∗i . (2.4) Substituting (2.4) into the second equation of (2.3) leads to 2∑ i=1 (βi(1− ui)γiS∗i ε − (κi + µi)(αi + µi) αi ) I∗i = 0. 6 S. MUSHAYABASA, A. A. E. LOSIO, C. MODNAK, J. WANG EJDE-2020/70 We can write 2∑ i=1 ( βi(1− ui)2γiαiS ∗ i ε(κi + µi)(αi + µi) − 1 ) (κi + µi)(αi + µi) αi I∗i = 0. (2.5) Equation (2.5) implies that I∗i = 0, or S∗i = ε(κi + µi)(αi + µi) βi(1− ui)2γiαi = Λi µiRi . From the first equation of (2.3) we have 2∑ i=1 ( Λi − βi(1− ui)2 Λi µiRi γi ε I∗i − Λi Ri ) = 0, and this gives I∗i = εµiRi βi(1− ui)2γi ( 1− 1 Ri ) . From (2.4) we have E∗i = εµi(κi + µi)Ri βi(1− ui)2αiγi ( 1− 1 Ri ) , W ∗ = 2∑ i=1 µiRi βi(1− ui)2 ( 1− 1 Ri ) . As we can note, E∗ exists whenever, S∗i > 0, E∗i > 0, I∗i > 0 and W ∗ > 0 and this only feasible if R0 > 1. Therefore we conclude that the endemic equilibrium E∗ exists and is unique whenever R0 > 1. � Now, we are aware that model (2.1) has two unique equilibrium points and the existence of a unique endemic equilibrium implies that the disease will persist in the community if R0 > 1. Next, we carry out sensitivity analysis of R0 to infer the role of awareness campaigns on the average number of secondary cases that will be generated in the community. From the expression of the basic reproduction number (2.2) shows that changes in behavior of individuals from either patch due to awareness can lead to extinction or persistence of the disease, since lim (u1→1,u2→1) R0 = 0, and the reverse case leads to lim (u1→0,u2→0) R0 = R∗0 (2.6) with R∗0 = β1α1γ1Λ1 εµ1m1m3 + β2α2γ2Λ2 εµ2m2m4 , note that R0 < R∗0, for ui 6= 0. In what follows, we utilize numerical simulations to demonstrate the effects of different awareness levels on the basic reproduction number R0. As presented on Table 1, let patch 1 represent the less risk population while patch 2 represents the high risk population, such that disease transmission in the risk patch is three times that of low risk patch, that is, β2 = 3β1. In Figure 1, a contour plot of the basic reproduction number R0 as a function of u1 and u2 is presented. Baseline values of other model parameters are based on Table 1. We can observe that there are many combinations of u1 and u2 that can lead to disease eradication, for instance if u1 = 0.4, then the disease can be EJDE-2020/70 TWO-PATCH GUINEA WORM DISEASE MODEL 7 Figure 1. Contour plot illustrating the effects of different levels of awareness per patch on the basic reproduction number. Baseline values for parameters used are in Table 1 eliminated from the community if u2 ≥ 0.3, however, if u2 < 0.3 then the infection will persist. Overall, we can conclude that if ui ≥ 0.4 in both patches then the infection will not persist in either of the patches. In Figure 2 shows the effects of varying disease transmission rate and educational campaigns on the magnitude of the basic reproduction number. By varying ui (with u1 = u2) and β1, we observed that an increase in transmission rate in both patches will require an increase in educational campaigns for the disease to be eradicated, for instance when β1 = 4 × 10−4 then educational campaigns equivalent to 0.7 or above will be needed to reduce the basic reproduction number to a value below unity. Figure 2. Effects of varying disease transmission rate and edu- cational campaigns on the basic reproduction number. Parameter values used are in Table 1 Next, we determine the global stability of the disease-free equilibrium and the endemic equilibrium. The global stability of these equilibrium points will be inves- tigated with the aid of Lyapunov functionals. 8 S. MUSHAYABASA, A. A. E. LOSIO, C. MODNAK, J. WANG EJDE-2020/70 Theorem 2.2. If R0 ≤ 1, then the DFE is globally asymptotically stable in Ω. Proof. Consider the Lyapunov function L(t) = 2∑ i=1 [aiEi(t) + biIi(t)] + cW (t), (2.7) where a1 = α1γ1m2m4, a2 = α2γ2m1m3, b1 = γ1m1m2m4, b2 = γ2m1m2m3, c = m1m2m3m4. Differentiating L(t) along the solutions of the system (2.1) gives L̇(t) = 2∑ i=1 ( aiĖi(t) + biİi(t) ) + cẆ (t). Since Si ≤ Ni ≤ Λi/µi for i = 1, 2, it follows that L̇(t) = εm1m2m3m4 ( 2∑ i=1 βi(1− ui)2αiγiSi ε(µi + αi)(µi + κi) − 1 ) W (t) ≤ εm1m2m3m4 ( 2∑ i=1 βi(1− ui)2αiγiΛi εµi(µi + αi)(µi + κi) − 1 ) W (t) = εm1m2m3m4 (R0 − 1)W (t) ≤ 0 provided that R0 ≤ 1. When R0 < 1, L̇ = 0 yields W = 0. Then it can be easily observed from the system (2.1) that as t → ∞, Ei → 0, Ii → 0, and Si → S0 i , for i = 1, 2. Hence, the only invariant set when L̇ = 0 is the singleton E0 = (S0 1 , S 0 2 , 0, 0, 0, 0, 0). It follows from Lasalle’s Invariance Principle [20] that every solution of the system (2.1), with initial conditions in Ω, approaches E0 as t→∞. Thus the equilibrium point E0 is globally asymptotically stable if R0 < 1. When R0 = 1, L̇ = 0 implies either W = 0, or 1 = 2∑ i=1 βi(1− ui)2αiγiSi ε(µi + αi)(µi + κi) ≤ 2∑ i=1 βi(1− ui)2αiγiΛi εµi(µi + αi)(µi + κi) = R0 = 1. The latter case yields Si = S0 i and, consequently, Ei = Ii = W = 0 for i = 1, 2. Hence, in either case, the only invariant set for L̇ = 0 is the singleton E0 = (S0 1 , S 0 2 , 0, 0, 0, 0, 0). Therefore the equilibrium point E0 is globally asymptotically stable if R0 = 1. � In contrast, if R0 > 1, then by continuity, L̇(t) > 0 in a neighbourhood of E0 in Ω̂. Solutions in Ω̂ sufficiently close to E0 move away from E0 implying that the DFE is unstable. Applying a uniform persistence result from [11] and an argument as in the proof of [16, Prop. 3.3], it can be shown that, if R0 > 1, instability of E0 implies uniform persistence of the system (2.1). Next, we study the global stability of the positive endemic equilibrium E∗. The problem is generally difficult for systems with dimensions higher than two. Follow- ing the approach used in [8], we introduce the following assumptions for our model (2.1): EJDE-2020/70 TWO-PATCH GUINEA WORM DISEASE MODEL 9 (A1) There exist a family of functions Φi(Si) : (0, Λi µi ]→ R+, i = 1, 2, Such that for all Si, Ei, W > 0, [Si − S∗i ] [Φi(Si)− Φi(S ∗ i )] ≥ 0,[ fi(Si,W )Φi(S ∗ i ) fi(S∗i ,W ∗)Φi(Si) − 1 ][ 1− fi(S ∗ i ,W ∗)Φi(Si)Ei fi(Si,W )Φi(S∗i )E∗i ] ≤ 0. (A2) For all Ei, Ii > 0, 1 ≤ i ≤ 2,[Ei E∗i − 1 ][ 1− E∗i Ii EiI∗i ] ≤ 0. (2.8) (A3) For all Ii, W > 0, 1 ≤ i ≤ 2,[ gi(Ii) gi(I∗i ) − 1 ][ 1− gi(I ∗ i )W gi(Ii)W ∗ ] ≤ 0. These assumptions are motivated by the proof of [8, Theorem 6.1]. The functions Φi, for example, can be taken as Φi(Si) = Si for i = 1, 2. We now prove that if R0 > 1 then the endemic equilibrium point is globally asymptotically stable. Theorem 2.3. If R0 > 1 then the endemic equilibrium point of system E∗ exists and is globally asymptotically stable. Proof. Consider the Lyapunov function V(t) = 2∑ i=1 {∫ Si S∗ i Φi(ξ)− Φi(S ∗ i ) Φi(ξ) dξ + [ Ei − E∗i − E∗i ln Ei E∗i ] + [ Ii − I∗i − I∗i ln Ii I∗i ]} + [ W −W ∗ −W ∗ ln W W ∗ ] . At endemic equilibrium we have the following identities Λi = fi(S ∗ i ,W ∗) + µiS ∗ i , (αi + µi)E ∗ i = fi(S ∗ i ,W ∗), (κi + µi)I ∗ i = αiE ∗ i , εW ∗ = 2∑ i=1 g(I∗i ). (2.9) Differentiating V along the solutions of (2.1) and using the equilibrium equations (2.9), we obtain dV dt = 2∑ i=1 {[ 1− Φi(S ∗ i ) Φi(Si) ]dSi dt + [ 1− E∗i Ei ]dEi dt + [ 1− I∗i Ii ]dIi dt } + [ 1− W ∗ W ]dW dt = 2∑ i=1 {[ 1− Φi(S ∗ i ) Φi(Si) ] [fi(S ∗ i ,W ∗) + µiS ∗ i − fi(Si,W )− µiSi] + [ fi(Si,W )− fi(Si,W ) E∗i Ei − fi(S∗i ,W ∗) Ei E∗i + fi(S ∗ i ,W ∗) ] + [ αiEi − αiEi I∗i Ii − αiE∗i Ii I∗i + αiE ∗ i ] 10 S. MUSHAYABASA, A. A. E. LOSIO, C. MODNAK, J. WANG EJDE-2020/70 + [ gi(Ii)− gi(Ii) W ∗ W − gi(I∗i ) W W ∗ + gi(I ∗ i ) ]} . After some algebraic manipulations, we have dV dt = 2∑ i=1 { −µi Φi(Si) [Si − S∗i ][Φi(Si)− Φi(S ∗ i )] + fi(S ∗ i ,W ∗) [ fi(Si,W )Φi(S ∗ i ) fi(S∗i ,W ∗)Φi(Si) − 1 ][ 1− fi(S ∗ i ,W ∗)Φi(Si)Ei fi(Si,W )Φi(S∗i )E∗i ] + fi(S ∗ i ,W ∗) [ 3− Φi(S ∗ i ) Φi(Si) − fi(Si,W )E∗i fi(S∗i ,W ∗)Ei − fi(S ∗ i ,W ∗)Φi(Si)Ei fi(Si,W )Φi(S∗i )E∗i ] + αiE ∗ i [Ei E∗i − 1 ][ 1− E∗i Ii EiI∗i ] + αiE ∗ i [ 2− EiI ∗ i E∗i Ii − E∗i Ii EiI∗i ] + gi(I ∗ i ) [ gi(Ii) gi(I∗i ) − 1 ][ 1− gi(I ∗ i )W gi(Ii)W ∗ ] + gi(I ∗ i ) [ 2− gi(Ii)W ∗ gi(I∗i )W − gi(I ∗ i )W gi(Ii)W ∗ ]} . It follows from assumptions (A1)–(A3) that dV dt ≤ 2∑ i=1 { fi(S ∗ i ,W ∗) [ 3− Φi(S ∗ i ) Φi(Si) − fi(Si,W )E∗i fi(S∗i ,W ∗)Ei − fi(S ∗ i ,W ∗)Φi(Si)Ei fi(Si,W )Φi(S∗i )E∗i ] + αiE ∗ i [ 2− EiI ∗ i E∗i Ii − E∗i Ii EiI∗i ] + gi(I ∗ i ) [ 2− gi(Ii)W ∗ gi(I∗i )W − gi(I ∗ i )W gi(Ii)W ∗ ]} ≤ 0. The terms inside the square brackets are non-positive due to the fact that the geometric mean is less than or equal to the arithmetic mean: Φi(S ∗ i ) Φi(Si) + fi(Si,W )E∗i fi(S∗i ,W ∗)Ei + fi(S ∗ i ,W ∗)Φi(Si)Ei fi(Si,W )Φi(S∗i )E∗i ≥ 3 3 √ Φi(S∗i ) Φi(Si) . fi(Si,W )E∗i fi(S∗i ,W ∗)Ei . fi(S∗i ,W ∗)Φi(Si)Ei fi(Si,W )Φi(S∗i )E∗i = 3, EiI ∗ i E∗i Ii + E∗i Ii EiI∗i ≥ 2 √ EiI∗i E∗i Ii · E ∗ i Ii EiI∗i = 2, gi(Ii)W ∗ gi(I∗i )W + gi(I ∗ i )W gi(Ii)W ∗ ≥ 2 √ gi(Ii)W ∗ gi(I∗i )W . gi(I∗i )W gi(Ii)W ∗ = 2. Meanwhile, it can be easily observed that the largest invariant set where V̇ = 0 is the singleton {E∗}. Hence, by LaSalle’s invariance principle, E∗ is globally asymptotically stable. This completes the proof of Theorem. � In Figure 3 phase portraits were constructed to illustrate the convergence of solutions to the disease-free and endemic equilibrium for R0 < 1 and R0 > 1, respectively. In Figure 3 (a) one can observe that all the three solution orbits converge to the disease-free equilibrium over time whenever R0 < 1 as guaranteed by Theorem 2.2. However, whenever R0 > 1 we can also note that the solution orbits converge to the endemic point, as illustrated in Figure 3 (b). EJDE-2020/70 TWO-PATCH GUINEA WORM DISEASE MODEL 11 0 2000 4000 6000 8000 10000 S 1 (t) 0 50 100 150 200 250 300 350 E 1(t ) E 1 (0)=150 E 1 (0)=250 E 1 (0)=350 0 100 200 300 400 500 E 1 (t) 0 50 100 150 200 250 300 350 I 1(t ) E 1 (0)=300 E 1 (0)=400 E 1 (0)=500 (a) (b) Figure 3. Phase portraits illustrating the converge of solutions to the (a) disease-free equilibrium when R0 < 1, (b) endemic equilib- rium when R0 > 1. To construct the phase portrait in (a) we set u1 = 0.55 and u2 = 0.45 to obtain R0 = 0.983 and for (b) we set u1 = 0.3 and u2 = 0.1 to obtain R0 = 2.56. In addition, the initial population levels were set to S(0)i = 1000, E2(0) = 0, I1 = I2 = 0 and W (0) = 500. In (a) we can observe that all the three solution orbits converge to the disease-free equilibrium demonstrating that whenever R0 < 1 the disease dies out in the community as guar- anteed by Theorem 2.2. However, if R0 > 1 we can see that the disease persists since the solution orbits converge to the endemic point. 3. Formulation and analysis of the optimal control problem In this section, we apply optimal control theory to identify optimal educational campaign strategies for GWD management. Thus the constant educational cam- paigns parameter ui considered earlier in model (2.1), is now assumed to be time dependent, i.e. ui(t) i = 1, 2. Utilizing the same variables and parameters names as before, our model with time dependent controls is given by Ṡi(t) = Λi − βi(1− ui(t))Si(t)W (t)− µiSi(t), Ėi(t) = βi(1− ui(t))Si(t)W (t)− (αi + µi)Ei(t), İi(t) = αiEi(t)− (κi + µi)Ii(t), Ẇ (t) = γ(1− u1(t))I1 + γ(1− u2(t))I2 − εW. (3.1) Our goal is to minimize the numbers of clinically infected individuals over a finite time horizon [0, T ] at minimal costs in each patch. Mathematically, we formulate an objective functional J(u(t)) as follows: J ( u1(t), u2(t) ) = ∫ T 0 ( A1I1(t) + B1 2 u2 1(t) +A2I2(t) + B2 2 u2 2(t) ) dt , (3.2) where Ai and Bi, i = 1, 2, are balancing coefficients (positive) transferring the integral into monetary quantity over a finite period of time. Note that the control efforts are assumed to be nonlinear, due to a number of advantages associated with a nonlinear function on the control. One of the advantages is that a nonlinear 12 S. MUSHAYABASA, A. A. E. LOSIO, C. MODNAK, J. WANG EJDE-2020/70 control allows the Hamiltonian to attain its minimum over the control set at a unique point. Moreover, a quadratic structure in the control has mathematical advantages. We seek an optimal control pair (u∗1, u ∗ 2) ∈ U such that J(u∗1, u ∗ 2) = inf (u1,u2)∈U J(u1, u2) (3.3) for the admissible set U = {(u1, u2) ∈ (L∞(0, T ))2 : 0 ≤ ui ≤ ai; ai ∈ R+, i = 1, 2}. The following theorem states the existence of the solution of the system (3.1) as well as their non-negativity and boundedness. Theorem 3.1. Given u = (u1, u2) ∈ U , there exists a non-negative bounded so- lution (Si, Ei, Ii,W ) i = 1, 2 to the state system (3.1) on the finite interval [0, T ] with given initial conditions. Proof. Since the controls and the state variables are uniformly bounded and non- negative on the finite interval [0, T ], there exists a minimizing sequence (un1 , u n 2 ) such that lim n→∞ J(un1 , u n 2 ) = inf (u1,u2)∈U J(u1, u2). We denote the corresponding sequence of the state variables by (Si, Ei, Ii,Wi) i = 1, 2. Further, since all state and control variables are bounded, then deriva- tives of the state variables are also bounded and it follows that all state variables are Lipschitz continuous with the same Lipschitz constant. Thus the sequence (Si, Ei, Ii,Wi) i = 1, 2 is uniformly equicontinuous in [0, T ]. Therefore by the Arzela-Ascoli Theorem [21, 27], the state sequence has a subsequence that con- verges uniformly to (S∗i , E ∗ i , I ∗ i ,W ∗ i ) i = 1, 2 in 0, T ]. Hence the control sequence unn = (un1 , u n 2 ) has a subsequence that converges weakly in L2(0, T ). Let (u∗1, u ∗ 2) ∈ U be such that uni ⇀ u∗i weakly in L2(0, T ) for i = 1, 2. utilizing the lower semi- continuity norms in weak L2, we obtain ‖u∗i ‖2L2 ≤ lim n→∞ ‖uni ‖2L2 , fori = 1, 2. (3.4) Hence, J(u∗1, u ∗ 2) ≤ lim n→∞ ∫ T 0 [ A1I n 1 (t) + B1 2 un1 (t) +A2I n 2 (t) + B2 2 un2 (t) ] dt = lim n→∞ J(u1, u2). Thus, we conclude that there exists a pair of control (u∗1, u ∗ 2) that minimizes the objective functional J(u1, u2). � By utilizing the result from Lukes [21] one can easily verify the existence and uniqueness of solutions for the state system (3.1) with a given control pair. Since there exists an optimal control from Theorem 3.1, now we can apply Pontryagin’s Maximum Principle [28]. Thus, the optimal control system (3.1) is converted into an equivalent problem of minimizing the Hamiltonian H(t) = A1I1(t) +A2I2(t) + B1 2 u2 1(t) + B2 2 u2 2(t) + λS1 { Λ1 − [1− u1(t)]β1S1W − µ1S1 } + λS2 { Λ2 − [1− u2(t)]β2S2W − µ2S2 } EJDE-2020/70 TWO-PATCH GUINEA WORM DISEASE MODEL 13 + λE1 { [1− u1(t)]β1S1W − (α1 + µ1)E1 } + λE2 { (1− u2(t))β2S2W − (α2 + µ2)E2 } + λI1 { α1E1 − (κ1 + µ1)I1 } + λI2 { α2E2 − (κ2 + µ2)I2 } + λW { (1− u1(t))γ1I1 + (1− u2(t))γ2I2 − εW } , where λSi (t), λEi (t), λIi(t) (i = 1, 2) and λW (t) denote the adjoint functions as- sociated with the states Si, Ei, Ii and W , i = 1, 2 respectively. Note that, in H, each adjoint function multiplies the right-hand side of the differential equation of its corresponding state function. The first term in H comes from the integrand of the objective functional. Given an optimal control pair u = (u1, u2) ∈ U and corresponding states (Si, Ei, Ii,W ), i = 1, 2 there exist adjoint functions [15] satisfying dλSi (t) dt = −∂H ∂Si , dλEi (t) dt = − ∂H ∂Ei , dλIi(t) dt = −∂H ∂Ii , dλW (t) dt = − ∂H ∂W . These equalities give λ̇S1 (t) = µ1λS1 (t) + β1[1− u1(t)][λS1 (t)− λE1 (t)]W (t), λ̇S2 (t) = µ2λS2 (t) + β2[1− u2(t)][λS2 (t)− λE2 (t)]W (t), λ̇E1(t) = µ1λE1(t) + α1[λE1(t)− λI1(t)], λ̇E2 (t) = µ2λE2 (t) + α2[λE2 (t)− λI2(t)], λ̇I1(t) = −A1 + (κ1 + µ1)λI1(t)− (1− u1)γ1λW (t), λ̇I2(t) = −A2 + (κ2 + µ2)λI2(t)− (1− u2)γ2λW (t), λ̇W (t) = ελW (t) + β1[1− u1(t)][λS1(t)− λE1(t)]S1(t) + β2[1− u2(t)][λS2(t)− λE2(t)]S2(t), (3.5) with transversality conditions λj(T ) = 0 for j = Si(t), Ei(t), Ii(t) (i = 1, 2) and W (t). Furthermore, the optimal controls are characterized by the optimality con- dition: u1(t) = min ( max ( 0, β1S1W (λE1 − λS1 ) + γ1I1λW B1 ) , a1 ) , u2 = min ( max ( 0, β2S2W (λE2 − λS2) + γ2I2λW B2 ) , a2 ) . (3.6) Finally, for sufficiently large T , one can show the uniqueness of an optimal control pair by following the approach in [10]. Our optimal control problem thus cou- ples the state system (3.1), the adjoint system (3.5), and the optimality condition (3.6). These equations are solved numerically using the forward-backward sweep- ing method [15], based on parameters listed in Table 1 and the following assumed initial population levels: Si(0) = 1000, Ei(0) = 50, Ii(0) = 100, W (0) = 500, i = 1, 2. In particular, we will find the simulation results under two different cases: Case A: Homogeneous setting, that is, disease characteristics are the same in the two patches. Case B: Heterogeneous setting, that is, different disease characteristics 14 S. MUSHAYABASA, A. A. E. LOSIO, C. MODNAK, J. WANG EJDE-2020/70 between the two patches. For each case, we find the total number of new infections in the two patches by the formula Tn = ∫ T 0 ( β1(1− u1)S1 + β2(1− u2)S2 ) Wdt. (3.7) Recall that the total cost associated with infected humans and the controls J , is given by equation (3.2). Case A: Homogeneous setting. To account for homogeneous setting, we assume equal disease transmission rate and educational campaigns, that is, we set β1 = β2 = 3.6 × 10−5, and 0 ≤ ui(t) ≤ 0.6. We set the initial guess of the controls as follows u1 = u2 = 0.3 and these values also correspond to educational campaign levels in the absence of time dependent intervention strategies. Furthermore, we set the weight constant to A1 = A2 = 1, B1 = B2 = 10, and the numerical illustrations are presented in Figure 4. From Figure 4, we observe that when disease characteristics (initial population levels and model parameters) are the same in the two patches then the dynamics of the disease will be similar in both patches. However, in the absence of optimal control we can see that the total number of exposed and infectious individuals will be higher compared to when there is optimal control. In particular, with optimal control implemented, the infected host population for both patches reduces remarkably to levels close to zero. Precisely, over a period of 120 months, the total number of new infections that will be generated from the two patches in the absence of optimal control is 1.2037×104, whereas, with optimal control approximately 408 infections and the total cost is J = 990. These results suggests that the presence of optimal control may lead to a reduction of approximately 1.1629×104 infections over a period of 120 months. Figure 4 (d) shows the optimal control profile for ui(t), i = 1, 2. As we can observe the control profile for ui starts from its maximum ui = 0.6 and remains there till the final time horizon, where it will drop to the origin. Thus, for effective GWD management time dependent educational campaigns will have to be maintained at maximum intensity for the entire horizon suggesting that perhaps, educational campaigns may need to be combined with other less cost intervention strategies (such as filtering water using pieces of cloth) in order to minimize cumulative costs in the long run. Table 2. Total number of newly infected humans over 120 months and the total cost J with respect to control bounds, under the homogeneous scenario. The balancing constants were set to Ai = 1 and Bi = 10, for i = 1, 2. Upper bounds of ai Total cost J Number of new infections Tn 0.4 2.0256× 104 1.8389× 104 0.6 990.0536 408 0.8 915.1205 52 From Table 2, we can observe that as the upper bound of ui approaches 1, the total number of newly infected individuals and total cost decreases significantly. EJDE-2020/70 TWO-PATCH GUINEA WORM DISEASE MODEL 15 (a) (b) (c) (d) Figure 4. Simulation results for case A: homogeneous setting. The solid curves in (a)-(c) represent the population levels of the host and vector in the presence of optimal control while the dotted lines depict the population levels in the absence of optimal control. The time varying educational rate is plotted in (d), 0 ≤ ui(t) ≤ 0.6. Note that since equality is assumed for both patches the diseases dynamics are similar, thus E1(t) = E2(t), I1(t) = I2(t), u1(t) = u2(t). Case B: Heterogeneous setting. Although case A (Homogeneous setting), is highly impractical it serves to demonstrate the likely outcome under homogeneous conditions and allows comparison to be made with the real world scenario – the heterogeneous setting. Here, we will investigate the impact of heterogeneous dis- eases characteristics between the two patches on long term disease dynamics. In particular we will investigate heterogeneity on (i) disease transmission rates, with the assumption that β2 = 6β1 (ii) bounds of the controls a2 < a1 (iii) both disease transmission rates and upper bounds of the controls. For all the simulations in case B we set the initial guess for the control as u1 = 0.4,and u2 = 0.3 and these values also correspond to education campaign levels in the absence of optimal control. Figure 5 depicts the time evolution of the infected host and pathogen in patches 1 and 2, for the scenario with and without the optimal control, under heterogeneous disease transmission rates (β2 = 6β1, with β1 = 3.6 × 10−4) and similar bounds of the controls (0 ≤ ui(t) ≤ 0.8). As we have already observed, the presence of 16 S. MUSHAYABASA, A. A. E. LOSIO, C. MODNAK, J. WANG EJDE-2020/70 (a) (b) (c) (d) (e) (f) Figure 5. Solution of the model (asymptomatic individuals) with and without optimal control under heterogeneous setting, with 0 ≤ ui(t) ≤ 0.8. The control profile for u2(t) has been omitted since it has exactly the same pattern as u1(t) optimal control will be associated with extremely low new infections as compared to when they is no optimal control. Precisely, in the absence of optimal control approximately a total of 3.535 × 103 new infections will generated compared to approximately 98, over a period of 120 months. The total cost of implementing the strategy will be J = 1.258× 103. EJDE-2020/70 TWO-PATCH GUINEA WORM DISEASE MODEL 17 Table 3. Total number of newly infected humans over 120 months and the total cost J with respect to control bounds, under the heterogeneous disease transmission between the two patches. The balancing constants were set to Ai = 1 and Bi = 10, 0 ≤ ui(t) ≤ 0.8 for i = 1, 2, with β1 = 3.6× 10−4. Transmission rate β2 Total cost J Number of new infections Tn 3β1 1.010× 103 62 6β1 1.258× 103 98 9β1 2.193× 103 256 12β1 5.735× 103 784 Table 3 depicts the effects of heterogeneous disease transmission rates between the two patches. As we can observe, an increase in disease transmission rate for the risk patch (patch 2) will lead to an increase on the total number of new infections and total cost over a period of 120 months. For instance, an increase on disease transmission rate from 6β1 to 9β1 will increase the total number of new infections by more than 100%. Furthermore, we can also observe that increasing disease transmission rate in patch 2 from 9β1 to 12β1 will also be associated with more than 100% increase on the total number of new infections and total costs over the same period. Numerical illustrations in Figure 6 illustrates the impact of different upper bound to controls u1 and u2 (that is, 0 ≤ u1 ≤ 0.8 and 0 ≤ u2 ≤ 0.5). As once observed earlier, the total number of new infections will be high in the absence of optimal control compared to when there is optimal control. In particular, approximately a total of 1.086 × 104 infections will be generated when there is no time dependent intervention strategies compared to approximately 260 in the presence of time de- pendent controls, over the same period of 120 months. The total cost associated with strategy implementation is J = 927.86. Figure 6 (d) shows the control pro- files for controls u1 and u2, and we can observe that all the control profiles starts at their maximum till the final time horizon where they will drop sharply to the origin. This suggests that for effective disease management educational campaigns will have to be maintained at their maximum intensities in both patches over 120 months. Simulations in Figure 7 demonstrates the effects of heterogeneous on: disease transmission rates, upper bound of the controls and initial population levels between the two patches, on long term disease dynamics. In particular, we have assumed higher disease transmission rate, lower upper bound to the control and higher initial population levels for patch 2-assumed to be the high risk community compared to patch 1. Under these circumstances we observed that approximately 8.51 × 103 new infections will be generated in the absence of optimal control compared to approximately 1.673×103 infections in the presence of optimal control over a period of 120 months. In addition, the total cost of implementing the time dependent strategies is J = 3.254× 103. Table 4 presents a summary of the variations in the weight parameters and their corresponding effects, under a heterogeneous scenario. Overall, we can observe that as these various weight constants increase the total cost J and the total number of new infections increases significantly. In particular, an increase in Bi leads 18 S. MUSHAYABASA, A. A. E. LOSIO, C. MODNAK, J. WANG EJDE-2020/70 (a) (b) (c) (d) (e) (f) Figure 6. Simulation results for case B, illustrating the dynamics of exposed individuals, infectious individuals and the pathogen in the environment, with and without optimal control, for 0 ≤ u1 ≤ 0.8, 0 ≤ u2 ≤ 0.5, β1 = β2 = 3.6 × 10−4. The solid curves in (a)-(e) represent the population levels of the host and pathogen in the presence of optimal control while the dotted lines depict the population levels without optimal control. to remarkable changes on the total cost and the total number of new infections compared to Ai, for i = 1, 2. EJDE-2020/70 TWO-PATCH GUINEA WORM DISEASE MODEL 19 (a) (b) (c) (d) (e) (f) Figure 7. Solution of the model (asymptomatic individuals) with and without optimal control under heterogeneous setting, with 0 ≤ u1(t) ≤ 0.8 and 0 ≤ u2(t) ≤ 0.5. Numerical results in Figure 8 suggest that when the cost of implementing the strategies is low Bi and the efforts are high Ai, then the disease can be effectively manage with control profiles at low intensity. Simulation results in Figure 9 depicts the impact of varying the weights. In particular, we varied the weights Ai, while Bi, i = 1, 2, is fixed to Bi = 10−2. As we can observe, increasing Ai may lead to varying optimal controls. Precisely, we can note that for value of A1 greater that 10, the control profile for u1 may not 20 S. MUSHAYABASA, A. A. E. LOSIO, C. MODNAK, J. WANG EJDE-2020/70 Table 4. Summary of variations in the weight parameters and their corresponding effects, we set β2 = 1.5β1, with β1 = 3.6 × 10−4, 0 ≤ u1 ≤ 0.8, 0 ≤ u2 ≤ 0.5. The initial population levels were assumed as follows, S1(0) = 1000, E1(0) = 50, I1(0) = 100, S2(0) = 1200, E1(0) = 250, I1(0) = 200 and W = 500. A1 A2 B1 B2 Total cost J New infections Tn 1 1 10 10 3.254× 103 1.6731× 103 2 2 10 10 5.9743× 103 1.6731× 103 1 1 102 102 8.0461× 103 1.7086× 103 2 2 102 102 1.0772× 104 1.6874× 103 2 2 103 103 5.8477× 104 1.8044× 103 5 5 103 103 6.6785× 104 1.7393× 103 need to be maintained at its maximum (0.8). However, in contrast we can note that the control profile for u2 may not need to be maintained at its maximum (0.5) whenever A2 is greater than 104. In a nutshell, one can conclude that whenever the control efforts are high and the cost of implementation is low, then the disease can be effectively controlled with control at varying and low intensity. 4. Conclusions In this article, a differential equation-based two-patch model for guinea worm disease, a water-borne infection, is proposed and analysed. The proposed model is is concerned with different host populations at the two patches, but with a common water source. Results obtained from the study demonstrate that, while each patch could contribute differently to the infection risk and the transmission dynamics, the overall reproduction number is determined by the sum of the two individual reproduction numbers. Further, analytical and numerical results of the study show that this overall reproduction number, R0, can provide a sharp threshold for the disease dynamics: when R0 ≤ 1, the disease-free equilibrium is globally asymptotically stable, indicating that the infection would die out; when R0 > 1, there exists a unique endemic equilibrium which is globally asymptotically stable, and the disease persists. Subsequently, we have performed an optimal control study on this two-patch model to effectively design control strategies for guinea worm disease. The main goal here, is to strengthen the efforts of educational campaigns so as to reduce the number of infectious human populations with minimal costs. From the analysis of the optimal control system observed that when the two patches have similar charac- teristics in disease transmission and control efforts, then the disease patterns will be same for the two patches. The reverse is true when the two patches have different disease characteristics. We also observed that the presence of optimal control will always be associated with lower total number of new infections compared to when there is no optimal control. However, the success of optimal control hinges on the total cost of implementing the strategies as well as its maximum strength. Further- more, we observed that an increase in weight constants may lead to a remarkable increase on the total number of new infections. Our study is not exhaustive and in future will consider the effects of combin- ing educational campaigns with other control measures for guinea worm disease, EJDE-2020/70 TWO-PATCH GUINEA WORM DISEASE MODEL 21 (a) (b) (c) (d) Figure 8. Simulation results for model system (2.1) with and without optimal control. We set 0 ≤ u1 ≤ 0.8, 0 ≤ u2 ≤ 0.5, Ai = 100 while Bi is fixed to 10−2, and the rest of the parameters are as in Table 1. The solid curves in (a)-(b) represent the population levels of the host in the presence of optimal control while the dotted lines depict the population levels without optimal control. Figure (c) and (d) shows the control profiles. such as water filter, larvicides, and surgery, in our optimal control study. These control methods, whether individually applied or combined together, might yield very different results with regard to the optimal control. Nevertheless, our findings illustrate that a strategic design of prevention and intervention methods for guinea worm disease that takes into account the unique setting of each different place may achieve the best outcome. Acknowledgments. The authors would like to thank the anonymous referee and the handling editor Professor S. Lenhart for their suggestions that have improved this article. References [1] Adetunde, I. A.; The Epidemiology of guinea worm infection in Tamale District, in the Northern Region of Ghana. Journalof Modern Mathematics and Statistics, (2008), 50-54. [2] Alew, O. O., Peters, N.; Ethiopian Dracunculiasis Eradication Program (EDEP) GOG WOREDA. Meeting, National Review; 14-15 December 2015; Gambella, Ethiopia: 2015. 22 S. MUSHAYABASA, A. A. E. LOSIO, C. MODNAK, J. WANG EJDE-2020/70 (a) (b) Figure 9. Simulation results for model system (2.1) illustrating the effects of varying the weights. We set 0 ≤ u1 ≤ 0.8, 0 ≤ u2 ≤ 0.5, and varied Ai while Bi is fixed to 10−2, and the rest of the parameters are as in Table 1. [3] Beyene, H. B.; Bekele, A; Shifara, A; Ebstie, Y. A.; Desalegn, Z; Kebede, Z; et al. Elimination of Guinea worm disease in Ethiopia; current status of the Disease’s, eradication strategies and challenges to the end game. Ethiop Med J, 55(Suppl 1), (2017), 15-31. [4] Biswas, G.; Sankara, D. P.; Agua-Agum, J.; Maiga, A.; Dracunculiasis (guinea worm disease): Eradication without a drug or a vaccine. Philos. Trans. R. Soc. B, 368, (2013), 20120146. [5] Bhunu, C. P.; Mushayabasa, S.; Kojouharov, H.; Tchuenche, J. M.; Mathematical Analysis of an HIV/AIDS Model: Impact of Educational Programs and Abstinence in Sub-Saharan Africa, J. Math. Model. Algorithms, 10 (2011), 31-55. [6] Cheesbrough, M.; Curtis, S.; District laboratory practice in tropical countries. second. Vol. 53, Journal of Chemical Information and Modeling. Cambridge: Cambrige University press, (2013), 1689-1699. [7] Cyclops an intermediate host for guinea worm. http: //www.who.int/water health/resources/vector324 to 336.pdf Accessed May 25, 2013. [8] Eisenberg, M. C.; Shuai. Z.; Tien, J. H.; van den Driessche, P.; A cholera model in a patchy environment with water and human movement, Math. Biosci. 246 (2013), 105-112. [9] https://water.org/our-impact/ethiopia/. Accessed 31 December 2018. [10] Fister, K. R.; Lenhart, S.; McNally, J. S.; Optimizing chemotherapy in an HIV model, Electron. J. Differ. Equ. textbf1998 no. 32 (1998), 1-12. [11] Freedman, H. I.; Ruan, S.; Tang, M.; Uniform persistence and flows near a closed positively invariant set, J. Dyn. Differ. Equat. 6 (1994), 583-600. [12] Hopkins, D. R.; Ruiz-Tiben, E.; Ruebush, T. K.; Diallo, N.; Agle, A.; Withers, C. J.; Dra- cunculiasis eradication: delayed, not denied. Am. J. Trop. Med. Hyg, 62 (2000), 163- 168 [13] Ilegbodu, V. A.; Kale, O. O.; Wise, R. A.; Christensen, B. L.; Steele, J. H.; Chambers, L. A.; Impact of guinea worm disease on children in Nigeria. The American Journal of Tropical Medicine and Hygiene, 35(5) (1986), 962-964. [14] Joshi, H. R.; Lenhart, S.; Hota, S.; Agusto, F.; Optimal control of an SIR model with changing behavior through an education campaign, Electron. J. Diff. Equ, 2015 no. 50 (2015), 1-4. [15] Lenhart, S.; Workman, J. T.; Optimal Control Applied to Biological Models, Chapman and Hall, 2007. [16] Li, M. Y.; Graef, J. R.; Wang, L.; Karsai. J.; Global dynamics of a SEIR model with varying total population size, Math. Biosci. 160 (1999), 191-213. [17] Link, K.; Guinea worm disease (Dracunculiasis): opening a mathematical can of worms. M.Sc. Thesis, Pennsylvania, USA, 2012. [18] Losio, A. E.; Mushayabasa, S.; Modeling the effects of spatial heterogeneity and seasonality on guinea worm disease transmission, Journal of Applied Mathematics, 2018 (2018), Article ID 5084687, 12 pages. EJDE-2020/70 TWO-PATCH GUINEA WORM DISEASE MODEL 23 [19] Lolika, O. P.; Modnak, C.; Mushayabasa, S.; On the dynamics of brucellosis infection in bison population with vertical transmission and culling, Math. Biosci, 305 (2018), 42-54. [20] LaSalle, J. S.; The stability of Dynamical Systems. CBMS-NSF Regional Conference Series in Applied Mathematics, SIAM: Philadelphia, 25 (1976), [21] Lukes, D. L.; Differential Equations: Classical to Controlled, Mathematics in Science and Engineering, Academic Press, New York, NY, USA, 1982. [22] Lutambi, A.; Penny, M.; Smith, T.; Chitnis, N.; Mathematical modelling of mosquito disper- sal in a heterogeneous environment, Math. Biosci, 214 (2013), 198-216. [23] Mari, L; Ciddio, M; Casagrandi, R; Perez-Saez, J; Bertuzzo, E; Rinaldo, A; et al.; Hetero- geneity in schistosomiasis transmission dynamics, emphJ. Theor. Biol, 432, (2017), 87-99. [24] Molyneux, D.; Sankara, D. P.; Guinea worm eradication: Progress and challenges-should we beware of the dog? PLoS Negl. Trop. Dis, 11, (2017), e0005495. [25] Mushayabasa, S.; Bhunu, C. P.; Smith?, R. J.; Assessing the impact of educational campaigns on controlling HCV among women in prison settings, Commun. Nonlinear. Sci. Numer. Simul, 17 (2011), 1714-1724. [26] Netshikweta, R; Garira, W.; A Multiscale Model for the World’s First Parasitic Disease Tar- geted for Eradication: Guinea Worm Disease, Comput. Math. Methods. Med., 2017 (2017), 1473287. [27] Pantha, B.; Day, J.; Lenhart, S.; Optimal control applied in an Anthrax epizootic model, J. Biol. Syst., 24 (2016), 495-517. [28] Pontryagin, L. S.; Boltyanskii, V. G.; Gamkrelidze, R. V.; Mishchenko, E. F.; The mathe- matical theory of optimal processes. Wiley, New Jersey, 1962. [29] Samanta, S.; Rana, S.; Sharma, A; Misra, A. K.; Chattopadhyay, J.; Effect of awareness programs by media on the epidemic outbreaks: A mathematical model, Appl. Math. Comp, 219 (2013), 6965-6977. [30] Sharma, A.; Misra, A. K.; Modeling the impact of awareness created by media campaigns on vaccination coverage in a variable population, J. Biol. Syst, 22 (2014), 249-270. [31] Smith, R. J.; Cloutier, P.; Harrison, J.; Desforges, A.; A Mathematical Model for the erad- ication of guinea worm disease. Understanding the dynamics of emerging and re-emerging infectious diseases using mathematical models, 2012: 133-156. [32] van den Driessche, P.; Watmough, J.; Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Math. Biosci., 180 (2002), 29- 48. Steady Mushayabasa (corresponding author) Department of Mathematics, University of Zimbabwe, P.O. Box MP 167, Harare, Zim- babwe Email address: steadymushaya@gmail.com Anthony A. E. Losio Department of Mathematics, University of Zimbabwe, P.O. Box MP 167, Harare, Zim- babwe Email address: abulelosio@yahoo.com Chairat Modnak Department of Mathematics, Faculty of Science, Naresuan University, Phitsanulok 65000, Thailand Email address: cmodn001@odu.edu Jin Wang Department of Mathematics, University of Tennessee at Chattanooga, 615 McCallie Ave., Chattanooga, TN 37403, USA Email address: Jin-Wang02@utc.edu 1. Introduction 2. Mathematical model 3. Formulation and analysis of the optimal control problem Case A: Homogeneous setting Case B: Heterogeneous setting 4. Conclusions Acknowledgments References