# Modeling the Effects of Spatial Heterogeneity and Seasonality on Guinea Worm Disease Transmission

Authors: Anthony A. E. Losio; Steady Mushayabasa
Journal: Journal of Applied Mathematics
 (2018)

---

## Abstract

Guinea worm disease is one of the neglected tropical diseases that is on the verge of elimination. Currently the disease is endemic in four countries, namely, Ethiopia, Mali, Chad, and South Sudan. Prior studies have demonstrated that climate factors and limited access to safe drinking water have a significant impact on transmission and control of Guinea worm disease. In this paper, we present a new mathematical model to understand the transmission dynamics of Guinea worm disease in South Sudan. The model incorporates seasonal variations, educational campaigns, and spatial heterogeneity. Both qualitative and quantitative analysis of the model have been carried out. Utilizing Guinea worm disease surveillance data of South Sudan (2007-2013) we estimate the model parameters. Meanwhile, we perform an optimal control study to evaluate the implications of vector control on long-term Guinea worm infection dynamics. Our results demonstrate that vector control could play a significant role on Guinea worm disease eradication.


---

## Body

## 1. Introduction

Guinea worm disease, also known as Dracunculiasis, is one of the neglected tropical diseases that is on the verge of elimination. The disease is transmitted to human via drinking contaminated water. Since the inception of the Guinea worm eradication program in the 1980s, the number of reported cases reduced from 3.5 million in 20 countries in 1986 to only 22 cases in 2015 from only 4 countries, namely, South Sudan, Chad, Mali, and Ethiopia [1]. Precisely, 99% of these cases were from South Sudan [2]. The disease has neither a vaccine to prevent nor medication to treat it. Despite being rarely fatal, individuals infected with the disease become nonfunctional for about 8.5 weeks [3].Effective control of the disease can be achieved through provision of safe drinking water, behavioral changes in patients, and communities and vector control. Although access to safe drinking water alone is known to be extremely integral on Guinea worm disease eradication [1], its provision remains a significant challenge in Guinea worm disease endemic countries. Recently the world health organization reported that only 60% of the residence of one village in Ethiopia had access to safe drinking water [4]. Limited access to safe drinking water was also noted to be a strong factor on the spread of Guinea worm disease in South Sudan [1].Such heterogeneity in individuals’ degree of susceptibility to infection has a strong impact on the transmission and control of Guinea worm disease. So far this aspect of heterogeneity has not been taken into account, leading to inadequate understanding of the influence of the spatial factors in the transmission and spread of Guinea worm disease. Another limitation in Guinea worm disease modeling is that the effects of seasonal variation has been inadequately addressed. In the Sahelian zone, transmission of Guinea worm disease has been observed to occur in the rainy season (May to August) while in the humid and forest zone, the peak occurs in the dry season (September to January) [5]. Such seasonal variations need to be incorporated in models that aim to inform policy makers effective and efficient ways to attain Guinea worm disease eradication.Mathematical models have become essential tools in mapping the spread and control of infectious diseases. By setting up a suitable epidemic model, we can improve our qualitative and quantitative understanding on the effects of preventative and control measures to aid the elimination of the disease. So far, few mathematical models have been proposed to understand the spread and control of Guinea worm disease (see, for example, [6–9]). Smith et al. [6] proposed a deterministic model for Guinea worm disease transmission and explored the role of disease preventative measures, namely, education, filtration, and chlorination. Utilizing the Latin Hypercube Sampling technique Smith and coworkers established that education is more effective toward the eradication of the disease. This work and the other aforementioned studies have undeniably produced many significant insights and improved our existing knowledge on Guinea worm disease dynamics.In this article, we propose a new nonautonomous model to explore the transmission and control of Guinea worm disease in a heterogeneous population of South Sudan. While our model takes a clue from the one by Smith et al. [6], in this article we have gone a step further by incorporating spatial heterogeneity and seasonal variations into a single framework for a comprehensive modeling of Guinea worm disease dynamics. To that end, the proposed model subdivides the susceptible human population into two categories, namely, the high risk and low risk. In addition, all epidemiological stages of the disease that are strongly influenced by seasonal variations have been modeled as periodic functions.We organize the remainder of the paper as follows. In the next section we present the methods and results of the study. In particular, we will formulate the model, compute the basic reproduction number, analyze the stability of the steady states, fit the model with Guinea worm disease data of South Sudan (2007-2013), and perform an optimal control study to explore the effects of vector control on Guinea worm disease eradication. Finally, we conclude the paper with some discussion in Section3.


## 2. Methods and Results

### 2.1. Model Framework

The model proposed in [6] subdivides the host population into categories of susceptible S(t), exposed E(t), infected individuals displaying clinical signs of the disease I(t), and the concentration of the bacteria in the environment W(t). We now extend this model by incorporating seasonal variations and variation in susceptibility to infection:(i)
Seasonal variations: to account for seasonal variations on Guinea worm disease dynamics we modeled the contact rate, incubation rate, pathogen shedding rate, and pathogen decay rate by the following periodic functions, respectively:(1)βt=β01+sin⁡2πtω-1,γt=γ01+sin⁡2πtω-1,αt=α01+sin⁡2πtω-1,δt=δ01+sin⁡2πtω-1,whereβ0, γ0,α0, and δ0 denote the contact rate, incubation rate, pathogen shedding rate, and pathogen decay rate, respectively, without seasonal forcing and ω>0 is the period.(ii)
Variation in susceptibility: in order to account for heterogeneity in our study, we subdivided the susceptible population into high and low risk individuals. High risk refers to susceptible individuals who have limited access to safe drinking water. As in [13, 14] we assume that low risk susceptible individuals have low chances of acquiring the infection.Keeping the above facts in mind, the dynamics of Guinea worm disease in this study are governed by the following nonautonomous system:(2)S˙1t=Aπ1-βtWtS1t-μ+a1S1t+a2S2t+pκσIt,S˙2t=Aπ2-βtWt1-θS2t-μ+a2S2t+a1S1t+pκ1-σIt,E˙t=βtWtS1t+1-θS2t-γt+μEt,I˙t=γtEt-μ+κIt,R˙t=1-pκIt-μRt,W˙t=αtIt-δtWt.All model variables and parameters are considered positive. Further the variablesS1(t), S2(t), E(t), I(t), and R(t) represent the susceptible high risk, susceptible low risk, exposed, infectious, and recovered individuals at time t, respectively, such that the total human population is given by N(t)=S1(t)+S2(t)+E(t)+I(t)+R(t). Meanwhile W(t) represents the population of the pathogen in the environment. Model parameters defined in (1) retain the same definitions.Model parameterA is the recruitment rate, π1 and π2 represent the fraction of new recruits into classes S1 and S2, respectively, μ is the natural mortality rate, a1 is the transfer from susceptible high risk to susceptible low risk, a2 is the transfer from susceptible low risk to susceptible high risk, θ is a modification factor that accounts for the impact of access of safe drinking water on disease transmission, and κ-1 is the average infectious period.In addition,p is the proportion of recovered individuals who become susceptible to infection again based on their behavior and the complementary part (1-p) represent recovered individuals who have extremely minimal chances of reinfection. People with limited access safe drinking water can minimize chances of infection by using a cloth filter or a pipe filter, to remove the copepods. Thus we assume that the pain and experience of recovered individuals have an impact on one’s behavior which in turn lead to one being reinfected or not. We further assume that a fraction σ of recovered individuals who move become susceptible to infection join S1(t) and the complementary (1-σ) join S2(t).Figure1 shows a flowchart depicting the dynamics of Guinea worm disease in this study.Figure 1
Diagram of the model structure.


### 2.2. The Reproduction Number

Now we explore the dynamical behavior of Guinea worm disease. Since the variableR(t) described by differential equation R˙(t)=(1-p)κI(t)-μR(t) does not appear in all the other five equations of system (2) it is sufficient to study the dynamics of Guinea worm disease from system:(3)S˙1t=Aπ1-βtWtS1t-μ+a1S1t+a2S2t+pκσIt,S˙2t=Aπ2-βtWt1-θS2t-μ+a2S2t+a1S1t+pκ1-σIt,E˙t=βtWtS1t+1-θS2t-γt+μEt,I˙t=γtEt-μ+κIt,W˙1t=αtIt-δtWt.System (3) has a disease-free equilibrium given by E0=(S10,S20,0,0,0), with(4)S10=Aa2+μπ1μμ+a1+a2,S20=Aa1+μπ2μμ+a1+a2.The severity of a disease can be measured by the basic reproduction number,R0, which is defined as the average number of secondary infections caused by a single infected individual in a completely susceptible population. Utilizing the next generation matrix approach [15] and adopting the notations therein, the matrices for new infections (denoted by F(t)) and the transfer terms (denoted by V(t)) at disease-free equilibrium are given by(5)Ft=00βtAa2+μπ1+1-θa1+μπ2μμ+a1+a20000αt0,Vt=μ+γt00-γtμ+κ000δt.Thus, the basic reproduction number of the time-averaged autonomous systems is(6)R0=ρF0V-10=β0γ0α0Aa2+μπ1+1-θa1+μπ2δ0μ+γ0κ+μμμ+a1+a20.5.In order to analyze the threshold dynamics of epidemiological models in periodic environments, Wang and Zhao [16] extended the framework in [15] by introducing the next infection operator(7)Lϕt=∫0∞Yt,t-sFt-sϕt-sdswhereY(t,s), t≥s, is the evolution operator of the linear ω-periodic system dy/dt=-V(t)y and ϕ(t), the initial distribution of infectious humans, is ω-periodic and always positive. The effective reproduction number for a periodic model is then determined by calculating the spectral radius of the next infection operator,(8)R0=ρL.Through direct calculation, the evolution operatorY(t,s), for model (3), is given by(9)Yt,s=Y~1t,s00Y~2t,se-μ+κt-s000Y~3t,s,where(10)Y~1t,s=e-μ+M1t-s,Y~2t,s=γ01M1-κ-11+M1-κ2M1-κsin⁡2πtω+cos⁡2πtω×e-μ+κt-s-e-μ+M1t-s,Y~3t,s=exp⁡-δ0t-s-ω2πcos⁡2πtω-cos⁡2πsωwith(11)M1=exp-γ01-ω2πcos2πωt-s.We can numerically evaluate the next infection operator by(12)Lϕt=∫0∞Yt,t-sFt-sϕt-sds=∫0ωGt-sϕt-sdswhere(13)Gt,s≈∑k=0MYt,t-s-kωFt-s-kω(14)Gt,s≈∑k=0M00βt-s-kωAa2+μπ1+1-θa1+μπ2μμ+a1+a2Y~1t,s00βt-s-kωAa2+μπ1+1-θa1+μπ2μμ+a1+a2Y~2t,s0αt-s-kωY~3t,s0.


### 2.3. Threshold Dynamics

In this section we will use the basic reproduction numberR0 defined in (8) to establish the threshold results, presented in Theorem 1, for model (3). To that end, we first note that R+2 is positively invariant for the following cooperative system:(15)S˙1t=Aπ1-μ+a1S1t+a2S2t,S˙2t=Aπ2-μ+a2S2t+a1S1t,and that (S10,S20) is the unique equilibrium solution which is globally attractive in R+2.Theorem 1.
(i)
IfR0<1, then the disease-free equilibrium E0 of system (3) is globally asymptotically stable.(ii)
IfR0>1, then system (3) admits at least one positive ω-periodic solution, and solutions of system (3) are uniformly persistent.Proof.
If(S1(t),S2(t),E(t),I(t),W(t)) is a nonnegative solution of (3), then we have(16)S˙1t≤Aπ1-μ+a1S1t+a2S2t,S˙2t≤Aπ2-μ+a2S2t+a1S1t.
Note that any nonnegative solution(S1(t),S2(t)) of system (3) approaches (S10,S20) as t→∞. It then follows from the standard comparison theorem (see, e.g., [17, Theorem A.4]) that, for any ϵ>0, there is a T>0 such that(17)Sit<Si0+ϵ,i=1,2,fort>T.Thus, for t>T, we have(18)E˙t≤βtWtS10+ϵ+1-θS20+ϵ-γt+μEt,I˙t≤γtEt-μ+κIt,W˙t≤αtIt-δtWt.
Define(19)Fϵt=00βtS10+ϵ+1-θS20+ϵ0000αt0.
By [16, Thorem 2.2], we have R0<1⇔ρ(ϕF-V(ω))<1, where ρ(ϕF-V(ω)) is the spectral radius of ϕF-V(ω) and ϕF-V(ω) is the monodromy matrix of the linear ω-periodic system dy/dt=(F-V)y. Then we can set ϵ sufficiently small such that ρ(ϕFϵ-V(ω))<1. As a consequence, the trivial solution (0,0,0) of the following linear ω-periodic system(20)E˙t=βtWtS10+ϵ+1-θS20+ϵ-γt+μEt,I˙t=γtEt-μ+κIt,W˙t=αtIt-δtWt
is globally asymptotically stable. Again by the comparison theorem, we know thatE(t)→0,I(t)→0 and W(t)→0 as t→∞. Finally, from the first and second equations of system (3) it follows that Si(t)→Si0 as t→∞ for i=1,2. This proves the result in part (i) of Theorem 1.
Next we shall focus on the caseR0>1. We define X=R+5,X0=R+2×Int(R+3),∂X0=X∖X0. It can easily be verified that both X and X0 are positively invariant. Let P:R+5→R+5 be the Poincaré map associated with system (3); that is, P(x0)=u(ω,x0) for all x0∈R+5, where u(t,x0) is the unique solution of (3) with u(0,x0)=x0. Set (21)M∂=S10,S20,E0,I0,W0∈∂X0:PmS10,S20,E0,I0,W0∈∂X0,∀m≥0,M=S1,S2,0,0,0:S1≥0,S2≥0.
We first demonstrate that(22)M∂=M.
It is evident thatM⊆M∂. For any (S1(0),S2(0),E(0),I(0),W(0))∈∂X0∖M, if E(0)=I(0)=0 and W(0)=0, then E˙(0)=β(t)W(t)[(S10+ϵ)+(1-θ)(S20+ϵ)]=β(0)W(0)[S10(0)+(1-θ)S2(0)]>0. If E(0)=W(0)=0 and I(0)>0, then E˙(0)>0 and W˙(0)=α(0)I(0)>0. If I(0)=W(0)=0 and E(0)>0 then E˙(0)>0 and I˙(0)=γ(0)E(0)>0. It follows that (S1(t),S2(t),E(t),I(t),W(t))∉∂X0 for 0<t≪1. By the positive invariance of X0, we know that Pm(S1(0),S2(0),E(0),I(0),W(0))∉∂X0 for m≥1; hence (S1(0),S2(0),E(0),I(0),W(0))∉M∂, and thus (22) holds.
Now consider the fixed pointM0=(S10,S20,0,0,0) of the Poincaré map P. Define WS(M0)={x0:Pm(x0)→M0,m→∞}. We show that(23)WSM0∩X0=∅.
Based on the continuity of solutions with respect to the initial conditions, for anyϵ>0, there exists ϑ>0 small enough such that for all (S1(0),S2(0),E(0),I(0),W(0))∈X0 with (24)S10,S20,E0,I0,W0-M0≤ϑ,we have(25)u(t,S10,S20,E0,I0,W0-ut,M0<ϵ,∀t∈0,ω.
To obtain (23), we claim that(26)limsupm→∞PmS10,S20,E0,I0,W0-M0≥ϑ,∀S10,S20,E0,I0,W0∈X0.
We prove this claim by contradiction; that is, we supposelimsupm→∞Pm(S1(0),S2(0),E(0),I(0),W(0))-M0<ϑ for some (S1(0),S2(0),E(0),I(0),W(0))∈X0. Without loss of generality, we assume that Pm(S1(0),S2(0),E(0),I(0),W(0))-M0<ϑ,∀m≥0. Thus,(27)u(t,PmS10,S20,E0,I0,W0-ut,M0<ϵ,∀t∈0,ω,m≥0.
Moreover, for anyt≥0, we write t=t0+nω with t0∈[0,ω) and n=[t/ω], the greatest integer less than or equal to t/ω. Then we obtain(28)u(t,S10,S20,E0,I0,W0-ut,M0=u(t0,PmS10,S20,E0,I0,W0-ut0,M0<ϵ
for anyt≥0. Let (S1(t),S2(t),E(t),I(t),W(t))=u(t,(S1(0),S2(0),E(0),I(0),W(0)). It follows that(29)-ϵ<S1t-S10<ϵ,-ϵ<S2t-S20<ϵ,0<Et<ϵ,0<It<ϵ,0<Wt<ϵ.(30)E˙t=βtWtS1t+1-θS2t-γt+μEt≥βtWtS10-ϵ+1-θS20-ϵ-γt+μEt=-γt+μEt-βtWtϵ2-θ+βtWtS10+1-θS20.
Hence we obtain(31)ddtEIW≥F-V-ϵKEIW,
whereF-V is given by (5) and(32)ϵ·K=ϵ·002-θβt000000.
Again based on [16, Theorem 2.2],  R0>1 if and only if ρ(ΦF-V(ω))>1. Thus, for ϵ small enough, we have ρ(ΦF-V-ϵ·K(ω))>1 which immediately yields the contradiction as (33)limt→∞Et=∞,limt→∞It=∞,limt→∞Wt=∞.
LetP1:R+2→R+2 be the Poincareé map associated with (15). Then (S10,S20) is globally attractive in R+2∖{0} for P1. It follows that M0 is isolated invariant set in X, and notice that WS(M0)∩X0=∅. Hence, every orbit in M∂ converges to M0 and M0 is acyclic in M∂. By [18, Theorem 1.3.1], for a stronger repelling property of ∂X0, we conclude that P is uniformly persistent with respect to (X0,∂X0), which implies the uniform persistence of the solutions of system (3) with respect to (X0,∂X0) [18, Theorem 3.1.1]. Consequently, based on [18, Theorem 1.3.6], the Poincaré map P has a fixed point (S~1(0),S~2(0),E~(0),I~(0),W¯(0))∈X0, with S~1(0)≠0,S~2(0)≠0. Thus, S~10,S~20,E~0,I~0,W¯0∈Int(R+5) and (S~1(t),S~2(t),E~(t),I~(t),W¯(t))=u(t,(S~1(0),S~2(0),E~(0), I~(0),W¯(0))) is a positive ω-periodic solution of the system.


### 2.4. Estimation of Parameter Values

In this section, we aim to use the Guinea worm disease surveillance data of South Sudan (2007-2013) to estimate periodic functions defined in model (3), while the baseline values for constant model parameters will be drawn from literature:(i)
The fraction of individuals recruited into susceptible risk populationπ1: according to the world health organization report of 2014, more than 90% of the people in South Sudan lives on less than US$1 a day [10]. Based on this assertion we assume that π1=0.9 and it follows that π2=0.1.(ii)
The natural mortality rateμ=1/lifeexpectancy: the united nations children’s fund report of 2013 [11] highlighted that life expectancy in South Sudan is 55 years; thus μ=0.0015month-1.(iii)
The mean infectious periodκ-1: as discussed in the introduction section, Guinea worm disease infectious individuals shed larvae into the environment during a short time and more often this occurs when they physically submerge their wound in open water sources. Based on findings in prior studies [6], we set κ=8.3×10-2month-1.(iv)
The rate of transfer from high risk susceptible population to low risk susceptible population and vice versa: since more than half of the population live below the international poverty line [11], we assume an extremely low transfer rate of individuals from class S1 to S2; thus we set a1=9.0×10-6month-1. In contrast, due to civil unrest which characterize South Sudan population we assume that the rate of transfer from low risk is higher (a1<a2); thus we set a2=0.009month-1.Table1 presents the description of model parameters and their baseline values.Table 1
Description of model parameters for system (3).
Symbol Definition Value Unit Source p Proportion of recovered individuals who become susceptible to the disease 0.3 - Assumed σ Proportion of recovered individuals who are at high risk to reinfection 0.9 - [10] a 1 Rate of transfer from high risk to low risk susceptible population 9.0 × 1 0 - 6 m o n t h - 1 Assumed a 2 Rate of transfer from low risk to high risk susceptible population 0.009 month-1 Assumed π 1 Proportion of new births in high risk susceptible population 0.9 - [10] π 2 Proportion of new births in low risk susceptible population 0.1 - [10] μ Natural mortality rate for humans 0.0015 month-1 [11] κ-1 Average infectious period 8.3 × 1 0 - 2 month Assumed θ Modification factor 0.98 - Assumed A Human birth rate 100 Individualsmonth-1 Assumed β0 Averaged environment–to–host transmission 1.2×10-5 L a r v a e     m o n t h - 1 Fitting γ0 Averaged incubation rate 1.0 × 1 0 - 1 m o n t h [12] α0 Averaged parasite shedding rate 2×10-1 Larvaemonth-1individual-1 Fitting δ0 Parasite death rate 2.16 month-1 [6]Now, we need to estimate the periodic functionsβ(t), α(t), γ(t), and δ(t). This will be done via the Fourier series analysis method as in [19, 20]. Observe that the monthly numbers of new Guinea worm cases in system (3) correspond to the term(34)gt=γtEt.Further, since variables and the periodic parameters in model system (3) are continuous functions of t, we will make use of trigonometric functions to fit g(t). Define(35)gt=c0+∑j=18cjcos⁡jLt+djsin⁡jLt,wherecj and dj(j=1,2,…,8) represent coefficients to be determined and L is the fundamental frequency. We used MATLAB software to determine the coefficients cj and dj(j=1,2,…,8) and we obtained(36)gt≈171.20-44.82cos⁡1.78t+8.38cos⁡3.56t-29.07cos⁡5.34t-3.08cos⁡7.12t-16.63cos⁡8.90t-51.20cos⁡10.68t+172.20cos⁡12.46+21.91cos⁡14.24t-153.3sin⁡1.78t+73.84sin⁡3.56t-59.83sin⁡5.34t+68.95sin⁡7.12t-62.73sin⁡8.90t+105.60sin⁡10.68t-2.17sin⁡12.46t-78.16sin⁡14.24t,withL=4π/7 and to be precise and exact ω=7/2. The comparison of the data with curve of g(t) is shown in Figure 2.Figure 2
The monthly numbers of new Guinea worm cases and its fitted curve.Next, we proceed to use (36) to estimate the periodic functions β(t) (disease transmission rate), γ(t) (incubation rate), α(t) (pathogen shedding rate), and δ(t) (pathogen decay rate). First, we assume the initial population levels as follows: S1(0)=S2(0)=10,000, W(0)=5000. Further, we set E(0)=I(0)=192 based on the number of cases reported in January 2007. Let(37)βt=β0~β1t,γt=γ0~γ1t,αt=α0~α1t,δt=δ0~δ1t,whereβ0~, γ0~, α0~, and δ0~ are related to the magnitude of seasonal fluctuations and β1(t), γ1(t), α1(t), and δ1(t) are periodic functions to be determined. After simulations and comparisons we obtained(38)β1t=10-3205-53.78cos⁡1.78t+10.06cos⁡3.56t-34.88cos⁡5.34t-3.70cos⁡7.12t-19.20cos⁡8.90t-61.44cos⁡10.68t+206.64cos⁡12.46+26.29cos⁡14.24t-183.96sin⁡1.78t+88.608sin⁡3.56t-71.80sin⁡5.34t+82.74sin⁡7.12t-75.28sin⁡8.90t+126.72sin⁡10.68t-2.60sin⁡12.46t-93.79sin⁡14.24t,α1t=10-3342.4-89.64cos⁡1.78t+16.76cos⁡3.56t-58.14cos⁡5.34t-6.16cos⁡7.12t-33.26cos⁡8.90t-102.4cos⁡10.68t+344.4cos⁡12.46+43.82cos⁡14.24t-306.6sin⁡1.78t+147.68sin⁡3.56t-119.66sin⁡5.34t+137.9sin⁡7.12t-125.46sin⁡8.90t+211.2sin⁡10.68t-4.34sin⁡12.46t-156.32sin⁡14.24t,γ1t=10-2256.8-67.23cos⁡1.78t+12.57cos⁡3.56t-43.61cos⁡5.34t-4.62cos⁡7.12t-24.95cos⁡8.90t-76.8cos⁡10.68t+258.3cos⁡12.46+32.87cos⁡14.24t-229.95sin⁡1.78t+110.76sin⁡3.56t-89.75sin⁡5.34t+103.43sin⁡7.12t-94.10sin⁡8.90t+158.4sin⁡10.68t-3.26sin⁡12.46t-117.24sin⁡14.24t,δ1t=10-4[369.79-96.81cos⁡1.78t+18.10cos⁡3.56t-62.79cos⁡5.34t-6.65cos⁡7.12t-35.92cos⁡8.90t-110.59cos⁡10.68t+371.95cos⁡12.46+47.33cos⁡14.24t-331.13sin⁡1.78t+159.49sin⁡3.56t-129.23sin⁡5.34t+148.93sin⁡7.12t-135.50sin⁡8.90t+228.10sin⁡10.68t-4.69sin⁡12.46t-168,83sin⁡14.24t.withβ0~=0.001, γ0~=0.15, α0~=0.006, and δ0~=0.00012. Using direct calculation and (37) one can easily derive the expressions for β(t), γ(t), α(t), and δ(t).


### 2.5. Optimal Control

Guinea worm disease has no treatment; however, the disease can be prevented by applying a chemical called temephos (which is an organophosphate), to unsafe drinking water sources in order to kill the copepods in water [7]. This approach reduces the parasite population and in turn it minimizes the likelihood of disease transmission. Now, we extend our initial model to incorporate a control function u(t) that represent the effect of vector control through the application of chemicals to water. The extended model takes the form(39)S˙1t=Aπ1-βtWtS1t-μ+a1S1t+a2S2t+pκσIt,S˙2t=Aπ2-βtWt1-θS2t-μ+a2S2t+a1S1t+pκ1-σIt,E˙t=βtWtS1t+1-θS2t-γt+μEt,I˙t=γtEt-μ+κIt,W˙t=αtIt-ut+δtWt.A successful control strategy is one that reduces the numbers of exposed and infected individuals over a finite time horizon[0,T] at minimal cost. Our objective functional is therefore formulated as (40)Jut=∫0TEt+It+Butdtsubject to the constraints of the ordinary differential equations in system (39) and the coefficient B (positive) represents a weight in the cost of the control. This objective functional and the differential equations are linear in the control with bounded states, and one can demonstrate by standard results that an optimal control and corresponding optimal states exist [21]. Although majority of applications of optimal control theory to infectious disease control consider controls in a quadratic form due to a number of mathematical advantages (e.g., if the control set is a compact and convex polyhedron it follows that the Hamiltonian attains its minimum over the control set at a unique point), here we considered a linear control since it is regarded as a more realistic function to use in a biological framework compared to the quadratic [22–24].By using Pontryagin’s maximum principle [21, 25] we derive necessary conditions for our optimal control and corresponding states. The Hamiltonian is(41)H=Et+It+But+λ1Aπ1-βtWtS1t-μ+a1S1t+a2S2t+pκσIt+λ2Aπ2-βtWt1-θS2t-μ+a2S2t+a1S1t+pκ1-σIt+λ3βtWtS1t+1-θS2t-μ+γEt+λ4γEt-μ+κIt+λ5αIt-ut+δWt.Given an optimal controlu(t), there exist adjoint functions, λ1,λ2,λ3,λ4,λ5 corresponding to the states S1,S2,E,I, and W such that(42)λ˙1=λ1βtW+μ+a1-λ2a1-λ3βtW,λ˙2=λ2βtW1-θ+μ+a2-λ1a2-λ3βtW1-θ,λ˙3=-1+λ3μ+γ-λ4γt,λ˙4=-1-λ1pκσ-λ2pκ1-σ+λ4μ+κ-λ5α,λ˙5=λ1βtS1+λ2βtS21-θ-λ3βtS1+1-θS2+λ5ut+δ,whereλ1(T)=0, λ2(T)=0, λ3(T)=0, λ4(T)=0, and λ5(T)=0 are the transversality conditions.The Hamiltonian is minimized with respect to the control variableu∗. Since the Hamiltonian is linear in the control, we must consider if the optimal control is bang-bang (at it is lower or upper bound), singular, or a combination. The singular case could occur if the slope or the switching function(43)∂H∂u=B-λ5Wis zero on nontrivial interval of time. Note that the optimal control would be at its upper bound or its lower bound according to(44)∂H∂u<0or>0.To investigate the singular case, let us suppose∂H/∂u=0 on some nontrivial interval. In this case we calculate(45)ddt∂H∂u=0and then we will show that control is not present in that equation. To solve for the value of the singular control, we will further calculate(46)d2dt2∂H∂u=0.We simplify the time derivative of∂H/∂u,(47)0=ddt∂H∂u=ddtB-λ5W=-λ5W′-λ5′W.We calculate both sums separately and add them together. The first sum can be written as(48)-λ5W′=-λ5αI-ut+δW=-αλ5I+ut+δλ5W.The second sum can be written as(49)-λ5′W=λ1-βtS1+λ2-βtS21-θ+λ3βtS1+1-θS2+λ5-ut-δW=-βtλ1S1W-βt1-θλ2S2W+βtλ3S1W+βt1-θλ3S2W-ut+δλ5W.Thus combining, we have(50)0=ddt∂H∂u=-αλ5I-βλ1S1W-β1-θλ2S2W+βλ3S1W+β1-θλ3S2W=-αλ5I+βλ3-λ1S1W+β1-θλ3-λ2S2W.We see that the control does not explicitly show in this expression, so next we calculate the second derivative with respect to time.(51)0=d2dt2∂H∂u=-αλ5′I+-αλ5I′+βλ3′-λ1′S1W+βλ3-λ1S1W′+S1′W+β1-θλ3′-λ2′S2W+β1-θλ3-λ2S2W′+S2′W.Using systems (3) and (41), we simplify (47) as follows:(52)0=d2dt2∂H∂u=α2βtλ3-λ1S1+λ3-λ21-θS2+μ+κλ5-δλ5I-αγλ5E+βtAλ3-λ1π1+λ3-λ21-θπ2+pκλ3-λ1σ+λ3-λ21-θ1-σIW-αλ5I+βλ3-λ1S1W+βλ3-λ21-θS2Wut+βt-1+γλ3-λ4-δλ3-λ1+a1λ2-λ3+λ3-λ21-θS1W+βt1-θ-1+γλ3-λ4-δλ3-λ2+a2λ1-λ31-θ+λ3-λ1S2W.The above equation can be written in the form(53)d2dt2∂H∂u=Φ1tut+Φ2t=0.and we can solve for the singular control as(54)usingulart=-Φ2tΦ1tif(55)Φ1t≠0,a≤-Φ2tΦ1t≤bwith(56)Φ1t=-αλ5I+βtλ3-λ1S1W+βtλ3-λ21-θS2W=-αBWI-βtλ3-λ1S1W-βtλ3-λ21-θS2Wand(57)Φ2t=α2βtλ3-λ1S1+λ3-λ21-θS2+μ+κλ5-δλ5I-αγBWE+βtAπ1λ3-λ1+λ3-λ21-θπ2W+βtpκσλ3-λ1+λ3-λ21-θ1-σIW+βt-1+γλ3-λ4-δλ3-λ1+a1λ2-λ3+a1λ3-λ21-θS1W+βt1-θ-1+γλ3-λ4-δλ3-λ2+a2λ1-λ31-θ+λ3-λ1S2W.To check the generalized Legendre-Clebsch condition for the singular control to be optimal, we required/dud2/dt2(∂H/∂u)=Φ1(t) to be negative [26, 27]. To summarize, our control characterization is as follows: on a nontrivial interval,(58)if∂H∂u<0attthenu∗t=bif∂H∂u>0attthenu∗t=aif∂H∂u<0attthenusingulart=-Φ2Φ1Hence, our control is optimal att provided Φ1(t)<0 and a≤-Φ2(t)/Φ1(t)≤b.Using parameter values from Table1 we investigate the effects of vector control on Guinea worm eradication in South Sudan. We solved the optimality system numerically using forward-backward sweep method [28]. Starting with an initial guess for the control, the state system is solved forward in time. Using those new state values, the adjoint system is solved backward in time. The control is updated using a convex combination of the old control values and the new control values from the characterization. The iterative method is repeated until convergence.Figure3 shows the concentration of the pathogen in the environment as a function of time, in both the presence and absence of the control. It is clear that the presence of optimal control can lead to eradication of the pathogen from the environment. More precisely, we note that, in the presence of vector control, the concentration of the pathogen within the environment will be reduced to a level close to 0 when t≥10 months. We also note that in the absence of the control the bacteria population has an oscillatory pattern and the amplitude of the oscillations appears to be constant.Figure 3
The concentration of the pathogen in the environment with and without the optimal control.Figure4 illustrates the numbers of exposed and infectious individuals over time with and without optimal control. The results clearly depict that the optimal control strategy significantly reduces the exposed and infectious populations (compared to the case without control), to a level close to 0 when t>30 months. More precisely, the number of exposed and infectious individuals will be reduced to a level close to 0 when t>30 and t>60, respectively. Meanwhile, we observe that the numbers of the pathogen in the environment and the infected population will not converge to zero at the same period. This is due to the long incubation period associated with Guinea worm disease.Figure 4
The numbers of exposed and infectious human population with and without the optimal control.
(a) (b)Figure5 shows the optimal control profile as a function of time, with B=2×10-5. As we can observe, the control profile starts from the maximum (u=0.8) initially and stays there for a period of 60 months and then it switches to its minimum (u=0.2) where it remains till the final time. Figure 6 shows the optimal control profile as a function of time, with high costs (B=2×105). It is evident that when the cost of the control is high, the control will have to be implemented at its maximum for a short period of time compared to a scenario when the costs are low.Figure 5
The control profile.Figure 6
The control profile with high costs.


## 2.1. Model Framework

The model proposed in [6] subdivides the host population into categories of susceptible S(t), exposed E(t), infected individuals displaying clinical signs of the disease I(t), and the concentration of the bacteria in the environment W(t). We now extend this model by incorporating seasonal variations and variation in susceptibility to infection:(i)
Seasonal variations: to account for seasonal variations on Guinea worm disease dynamics we modeled the contact rate, incubation rate, pathogen shedding rate, and pathogen decay rate by the following periodic functions, respectively:(1)βt=β01+sin⁡2πtω-1,γt=γ01+sin⁡2πtω-1,αt=α01+sin⁡2πtω-1,δt=δ01+sin⁡2πtω-1,whereβ0, γ0,α0, and δ0 denote the contact rate, incubation rate, pathogen shedding rate, and pathogen decay rate, respectively, without seasonal forcing and ω>0 is the period.(ii)
Variation in susceptibility: in order to account for heterogeneity in our study, we subdivided the susceptible population into high and low risk individuals. High risk refers to susceptible individuals who have limited access to safe drinking water. As in [13, 14] we assume that low risk susceptible individuals have low chances of acquiring the infection.Keeping the above facts in mind, the dynamics of Guinea worm disease in this study are governed by the following nonautonomous system:(2)S˙1t=Aπ1-βtWtS1t-μ+a1S1t+a2S2t+pκσIt,S˙2t=Aπ2-βtWt1-θS2t-μ+a2S2t+a1S1t+pκ1-σIt,E˙t=βtWtS1t+1-θS2t-γt+μEt,I˙t=γtEt-μ+κIt,R˙t=1-pκIt-μRt,W˙t=αtIt-δtWt.All model variables and parameters are considered positive. Further the variablesS1(t), S2(t), E(t), I(t), and R(t) represent the susceptible high risk, susceptible low risk, exposed, infectious, and recovered individuals at time t, respectively, such that the total human population is given by N(t)=S1(t)+S2(t)+E(t)+I(t)+R(t). Meanwhile W(t) represents the population of the pathogen in the environment. Model parameters defined in (1) retain the same definitions.Model parameterA is the recruitment rate, π1 and π2 represent the fraction of new recruits into classes S1 and S2, respectively, μ is the natural mortality rate, a1 is the transfer from susceptible high risk to susceptible low risk, a2 is the transfer from susceptible low risk to susceptible high risk, θ is a modification factor that accounts for the impact of access of safe drinking water on disease transmission, and κ-1 is the average infectious period.In addition,p is the proportion of recovered individuals who become susceptible to infection again based on their behavior and the complementary part (1-p) represent recovered individuals who have extremely minimal chances of reinfection. People with limited access safe drinking water can minimize chances of infection by using a cloth filter or a pipe filter, to remove the copepods. Thus we assume that the pain and experience of recovered individuals have an impact on one’s behavior which in turn lead to one being reinfected or not. We further assume that a fraction σ of recovered individuals who move become susceptible to infection join S1(t) and the complementary (1-σ) join S2(t).Figure1 shows a flowchart depicting the dynamics of Guinea worm disease in this study.Figure 1
Diagram of the model structure.


## 2.2. The Reproduction Number

Now we explore the dynamical behavior of Guinea worm disease. Since the variableR(t) described by differential equation R˙(t)=(1-p)κI(t)-μR(t) does not appear in all the other five equations of system (2) it is sufficient to study the dynamics of Guinea worm disease from system:(3)S˙1t=Aπ1-βtWtS1t-μ+a1S1t+a2S2t+pκσIt,S˙2t=Aπ2-βtWt1-θS2t-μ+a2S2t+a1S1t+pκ1-σIt,E˙t=βtWtS1t+1-θS2t-γt+μEt,I˙t=γtEt-μ+κIt,W˙1t=αtIt-δtWt.System (3) has a disease-free equilibrium given by E0=(S10,S20,0,0,0), with(4)S10=Aa2+μπ1μμ+a1+a2,S20=Aa1+μπ2μμ+a1+a2.The severity of a disease can be measured by the basic reproduction number,R0, which is defined as the average number of secondary infections caused by a single infected individual in a completely susceptible population. Utilizing the next generation matrix approach [15] and adopting the notations therein, the matrices for new infections (denoted by F(t)) and the transfer terms (denoted by V(t)) at disease-free equilibrium are given by(5)Ft=00βtAa2+μπ1+1-θa1+μπ2μμ+a1+a20000αt0,Vt=μ+γt00-γtμ+κ000δt.Thus, the basic reproduction number of the time-averaged autonomous systems is(6)R0=ρF0V-10=β0γ0α0Aa2+μπ1+1-θa1+μπ2δ0μ+γ0κ+μμμ+a1+a20.5.In order to analyze the threshold dynamics of epidemiological models in periodic environments, Wang and Zhao [16] extended the framework in [15] by introducing the next infection operator(7)Lϕt=∫0∞Yt,t-sFt-sϕt-sdswhereY(t,s), t≥s, is the evolution operator of the linear ω-periodic system dy/dt=-V(t)y and ϕ(t), the initial distribution of infectious humans, is ω-periodic and always positive. The effective reproduction number for a periodic model is then determined by calculating the spectral radius of the next infection operator,(8)R0=ρL.Through direct calculation, the evolution operatorY(t,s), for model (3), is given by(9)Yt,s=Y~1t,s00Y~2t,se-μ+κt-s000Y~3t,s,where(10)Y~1t,s=e-μ+M1t-s,Y~2t,s=γ01M1-κ-11+M1-κ2M1-κsin⁡2πtω+cos⁡2πtω×e-μ+κt-s-e-μ+M1t-s,Y~3t,s=exp⁡-δ0t-s-ω2πcos⁡2πtω-cos⁡2πsωwith(11)M1=exp-γ01-ω2πcos2πωt-s.We can numerically evaluate the next infection operator by(12)Lϕt=∫0∞Yt,t-sFt-sϕt-sds=∫0ωGt-sϕt-sdswhere(13)Gt,s≈∑k=0MYt,t-s-kωFt-s-kω(14)Gt,s≈∑k=0M00βt-s-kωAa2+μπ1+1-θa1+μπ2μμ+a1+a2Y~1t,s00βt-s-kωAa2+μπ1+1-θa1+μπ2μμ+a1+a2Y~2t,s0αt-s-kωY~3t,s0.


## 2.3. Threshold Dynamics

In this section we will use the basic reproduction numberR0 defined in (8) to establish the threshold results, presented in Theorem 1, for model (3). To that end, we first note that R+2 is positively invariant for the following cooperative system:(15)S˙1t=Aπ1-μ+a1S1t+a2S2t,S˙2t=Aπ2-μ+a2S2t+a1S1t,and that (S10,S20) is the unique equilibrium solution which is globally attractive in R+2.Theorem 1.
(i)
IfR0<1, then the disease-free equilibrium E0 of system (3) is globally asymptotically stable.(ii)
IfR0>1, then system (3) admits at least one positive ω-periodic solution, and solutions of system (3) are uniformly persistent.Proof.
If(S1(t),S2(t),E(t),I(t),W(t)) is a nonnegative solution of (3), then we have(16)S˙1t≤Aπ1-μ+a1S1t+a2S2t,S˙2t≤Aπ2-μ+a2S2t+a1S1t.
Note that any nonnegative solution(S1(t),S2(t)) of system (3) approaches (S10,S20) as t→∞. It then follows from the standard comparison theorem (see, e.g., [17, Theorem A.4]) that, for any ϵ>0, there is a T>0 such that(17)Sit<Si0+ϵ,i=1,2,fort>T.Thus, for t>T, we have(18)E˙t≤βtWtS10+ϵ+1-θS20+ϵ-γt+μEt,I˙t≤γtEt-μ+κIt,W˙t≤αtIt-δtWt.
Define(19)Fϵt=00βtS10+ϵ+1-θS20+ϵ0000αt0.
By [16, Thorem 2.2], we have R0<1⇔ρ(ϕF-V(ω))<1, where ρ(ϕF-V(ω)) is the spectral radius of ϕF-V(ω) and ϕF-V(ω) is the monodromy matrix of the linear ω-periodic system dy/dt=(F-V)y. Then we can set ϵ sufficiently small such that ρ(ϕFϵ-V(ω))<1. As a consequence, the trivial solution (0,0,0) of the following linear ω-periodic system(20)E˙t=βtWtS10+ϵ+1-θS20+ϵ-γt+μEt,I˙t=γtEt-μ+κIt,W˙t=αtIt-δtWt
is globally asymptotically stable. Again by the comparison theorem, we know thatE(t)→0,I(t)→0 and W(t)→0 as t→∞. Finally, from the first and second equations of system (3) it follows that Si(t)→Si0 as t→∞ for i=1,2. This proves the result in part (i) of Theorem 1.
Next we shall focus on the caseR0>1. We define X=R+5,X0=R+2×Int(R+3),∂X0=X∖X0. It can easily be verified that both X and X0 are positively invariant. Let P:R+5→R+5 be the Poincaré map associated with system (3); that is, P(x0)=u(ω,x0) for all x0∈R+5, where u(t,x0) is the unique solution of (3) with u(0,x0)=x0. Set (21)M∂=S10,S20,E0,I0,W0∈∂X0:PmS10,S20,E0,I0,W0∈∂X0,∀m≥0,M=S1,S2,0,0,0:S1≥0,S2≥0.
We first demonstrate that(22)M∂=M.
It is evident thatM⊆M∂. For any (S1(0),S2(0),E(0),I(0),W(0))∈∂X0∖M, if E(0)=I(0)=0 and W(0)=0, then E˙(0)=β(t)W(t)[(S10+ϵ)+(1-θ)(S20+ϵ)]=β(0)W(0)[S10(0)+(1-θ)S2(0)]>0. If E(0)=W(0)=0 and I(0)>0, then E˙(0)>0 and W˙(0)=α(0)I(0)>0. If I(0)=W(0)=0 and E(0)>0 then E˙(0)>0 and I˙(0)=γ(0)E(0)>0. It follows that (S1(t),S2(t),E(t),I(t),W(t))∉∂X0 for 0<t≪1. By the positive invariance of X0, we know that Pm(S1(0),S2(0),E(0),I(0),W(0))∉∂X0 for m≥1; hence (S1(0),S2(0),E(0),I(0),W(0))∉M∂, and thus (22) holds.
Now consider the fixed pointM0=(S10,S20,0,0,0) of the Poincaré map P. Define WS(M0)={x0:Pm(x0)→M0,m→∞}. We show that(23)WSM0∩X0=∅.
Based on the continuity of solutions with respect to the initial conditions, for anyϵ>0, there exists ϑ>0 small enough such that for all (S1(0),S2(0),E(0),I(0),W(0))∈X0 with (24)S10,S20,E0,I0,W0-M0≤ϑ,we have(25)u(t,S10,S20,E0,I0,W0-ut,M0<ϵ,∀t∈0,ω.
To obtain (23), we claim that(26)limsupm→∞PmS10,S20,E0,I0,W0-M0≥ϑ,∀S10,S20,E0,I0,W0∈X0.
We prove this claim by contradiction; that is, we supposelimsupm→∞Pm(S1(0),S2(0),E(0),I(0),W(0))-M0<ϑ for some (S1(0),S2(0),E(0),I(0),W(0))∈X0. Without loss of generality, we assume that Pm(S1(0),S2(0),E(0),I(0),W(0))-M0<ϑ,∀m≥0. Thus,(27)u(t,PmS10,S20,E0,I0,W0-ut,M0<ϵ,∀t∈0,ω,m≥0.
Moreover, for anyt≥0, we write t=t0+nω with t0∈[0,ω) and n=[t/ω], the greatest integer less than or equal to t/ω. Then we obtain(28)u(t,S10,S20,E0,I0,W0-ut,M0=u(t0,PmS10,S20,E0,I0,W0-ut0,M0<ϵ
for anyt≥0. Let (S1(t),S2(t),E(t),I(t),W(t))=u(t,(S1(0),S2(0),E(0),I(0),W(0)). It follows that(29)-ϵ<S1t-S10<ϵ,-ϵ<S2t-S20<ϵ,0<Et<ϵ,0<It<ϵ,0<Wt<ϵ.(30)E˙t=βtWtS1t+1-θS2t-γt+μEt≥βtWtS10-ϵ+1-θS20-ϵ-γt+μEt=-γt+μEt-βtWtϵ2-θ+βtWtS10+1-θS20.
Hence we obtain(31)ddtEIW≥F-V-ϵKEIW,
whereF-V is given by (5) and(32)ϵ·K=ϵ·002-θβt000000.
Again based on [16, Theorem 2.2],  R0>1 if and only if ρ(ΦF-V(ω))>1. Thus, for ϵ small enough, we have ρ(ΦF-V-ϵ·K(ω))>1 which immediately yields the contradiction as (33)limt→∞Et=∞,limt→∞It=∞,limt→∞Wt=∞.
LetP1:R+2→R+2 be the Poincareé map associated with (15). Then (S10,S20) is globally attractive in R+2∖{0} for P1. It follows that M0 is isolated invariant set in X, and notice that WS(M0)∩X0=∅. Hence, every orbit in M∂ converges to M0 and M0 is acyclic in M∂. By [18, Theorem 1.3.1], for a stronger repelling property of ∂X0, we conclude that P is uniformly persistent with respect to (X0,∂X0), which implies the uniform persistence of the solutions of system (3) with respect to (X0,∂X0) [18, Theorem 3.1.1]. Consequently, based on [18, Theorem 1.3.6], the Poincaré map P has a fixed point (S~1(0),S~2(0),E~(0),I~(0),W¯(0))∈X0, with S~1(0)≠0,S~2(0)≠0. Thus, S~10,S~20,E~0,I~0,W¯0∈Int(R+5) and (S~1(t),S~2(t),E~(t),I~(t),W¯(t))=u(t,(S~1(0),S~2(0),E~(0), I~(0),W¯(0))) is a positive ω-periodic solution of the system.


## 2.4. Estimation of Parameter Values

In this section, we aim to use the Guinea worm disease surveillance data of South Sudan (2007-2013) to estimate periodic functions defined in model (3), while the baseline values for constant model parameters will be drawn from literature:(i)
The fraction of individuals recruited into susceptible risk populationπ1: according to the world health organization report of 2014, more than 90% of the people in South Sudan lives on less than US$1 a day [10]. Based on this assertion we assume that π1=0.9 and it follows that π2=0.1.(ii)
The natural mortality rateμ=1/lifeexpectancy: the united nations children’s fund report of 2013 [11] highlighted that life expectancy in South Sudan is 55 years; thus μ=0.0015month-1.(iii)
The mean infectious periodκ-1: as discussed in the introduction section, Guinea worm disease infectious individuals shed larvae into the environment during a short time and more often this occurs when they physically submerge their wound in open water sources. Based on findings in prior studies [6], we set κ=8.3×10-2month-1.(iv)
The rate of transfer from high risk susceptible population to low risk susceptible population and vice versa: since more than half of the population live below the international poverty line [11], we assume an extremely low transfer rate of individuals from class S1 to S2; thus we set a1=9.0×10-6month-1. In contrast, due to civil unrest which characterize South Sudan population we assume that the rate of transfer from low risk is higher (a1<a2); thus we set a2=0.009month-1.Table1 presents the description of model parameters and their baseline values.Table 1
Description of model parameters for system (3).
Symbol Definition Value Unit Source p Proportion of recovered individuals who become susceptible to the disease 0.3 - Assumed σ Proportion of recovered individuals who are at high risk to reinfection 0.9 - [10] a 1 Rate of transfer from high risk to low risk susceptible population 9.0 × 1 0 - 6 m o n t h - 1 Assumed a 2 Rate of transfer from low risk to high risk susceptible population 0.009 month-1 Assumed π 1 Proportion of new births in high risk susceptible population 0.9 - [10] π 2 Proportion of new births in low risk susceptible population 0.1 - [10] μ Natural mortality rate for humans 0.0015 month-1 [11] κ-1 Average infectious period 8.3 × 1 0 - 2 month Assumed θ Modification factor 0.98 - Assumed A Human birth rate 100 Individualsmonth-1 Assumed β0 Averaged environment–to–host transmission 1.2×10-5 L a r v a e     m o n t h - 1 Fitting γ0 Averaged incubation rate 1.0 × 1 0 - 1 m o n t h [12] α0 Averaged parasite shedding rate 2×10-1 Larvaemonth-1individual-1 Fitting δ0 Parasite death rate 2.16 month-1 [6]Now, we need to estimate the periodic functionsβ(t), α(t), γ(t), and δ(t). This will be done via the Fourier series analysis method as in [19, 20]. Observe that the monthly numbers of new Guinea worm cases in system (3) correspond to the term(34)gt=γtEt.Further, since variables and the periodic parameters in model system (3) are continuous functions of t, we will make use of trigonometric functions to fit g(t). Define(35)gt=c0+∑j=18cjcos⁡jLt+djsin⁡jLt,wherecj and dj(j=1,2,…,8) represent coefficients to be determined and L is the fundamental frequency. We used MATLAB software to determine the coefficients cj and dj(j=1,2,…,8) and we obtained(36)gt≈171.20-44.82cos⁡1.78t+8.38cos⁡3.56t-29.07cos⁡5.34t-3.08cos⁡7.12t-16.63cos⁡8.90t-51.20cos⁡10.68t+172.20cos⁡12.46+21.91cos⁡14.24t-153.3sin⁡1.78t+73.84sin⁡3.56t-59.83sin⁡5.34t+68.95sin⁡7.12t-62.73sin⁡8.90t+105.60sin⁡10.68t-2.17sin⁡12.46t-78.16sin⁡14.24t,withL=4π/7 and to be precise and exact ω=7/2. The comparison of the data with curve of g(t) is shown in Figure 2.Figure 2
The monthly numbers of new Guinea worm cases and its fitted curve.Next, we proceed to use (36) to estimate the periodic functions β(t) (disease transmission rate), γ(t) (incubation rate), α(t) (pathogen shedding rate), and δ(t) (pathogen decay rate). First, we assume the initial population levels as follows: S1(0)=S2(0)=10,000, W(0)=5000. Further, we set E(0)=I(0)=192 based on the number of cases reported in January 2007. Let(37)βt=β0~β1t,γt=γ0~γ1t,αt=α0~α1t,δt=δ0~δ1t,whereβ0~, γ0~, α0~, and δ0~ are related to the magnitude of seasonal fluctuations and β1(t), γ1(t), α1(t), and δ1(t) are periodic functions to be determined. After simulations and comparisons we obtained(38)β1t=10-3205-53.78cos⁡1.78t+10.06cos⁡3.56t-34.88cos⁡5.34t-3.70cos⁡7.12t-19.20cos⁡8.90t-61.44cos⁡10.68t+206.64cos⁡12.46+26.29cos⁡14.24t-183.96sin⁡1.78t+88.608sin⁡3.56t-71.80sin⁡5.34t+82.74sin⁡7.12t-75.28sin⁡8.90t+126.72sin⁡10.68t-2.60sin⁡12.46t-93.79sin⁡14.24t,α1t=10-3342.4-89.64cos⁡1.78t+16.76cos⁡3.56t-58.14cos⁡5.34t-6.16cos⁡7.12t-33.26cos⁡8.90t-102.4cos⁡10.68t+344.4cos⁡12.46+43.82cos⁡14.24t-306.6sin⁡1.78t+147.68sin⁡3.56t-119.66sin⁡5.34t+137.9sin⁡7.12t-125.46sin⁡8.90t+211.2sin⁡10.68t-4.34sin⁡12.46t-156.32sin⁡14.24t,γ1t=10-2256.8-67.23cos⁡1.78t+12.57cos⁡3.56t-43.61cos⁡5.34t-4.62cos⁡7.12t-24.95cos⁡8.90t-76.8cos⁡10.68t+258.3cos⁡12.46+32.87cos⁡14.24t-229.95sin⁡1.78t+110.76sin⁡3.56t-89.75sin⁡5.34t+103.43sin⁡7.12t-94.10sin⁡8.90t+158.4sin⁡10.68t-3.26sin⁡12.46t-117.24sin⁡14.24t,δ1t=10-4[369.79-96.81cos⁡1.78t+18.10cos⁡3.56t-62.79cos⁡5.34t-6.65cos⁡7.12t-35.92cos⁡8.90t-110.59cos⁡10.68t+371.95cos⁡12.46+47.33cos⁡14.24t-331.13sin⁡1.78t+159.49sin⁡3.56t-129.23sin⁡5.34t+148.93sin⁡7.12t-135.50sin⁡8.90t+228.10sin⁡10.68t-4.69sin⁡12.46t-168,83sin⁡14.24t.withβ0~=0.001, γ0~=0.15, α0~=0.006, and δ0~=0.00012. Using direct calculation and (37) one can easily derive the expressions for β(t), γ(t), α(t), and δ(t).


## 2.5. Optimal Control

Guinea worm disease has no treatment; however, the disease can be prevented by applying a chemical called temephos (which is an organophosphate), to unsafe drinking water sources in order to kill the copepods in water [7]. This approach reduces the parasite population and in turn it minimizes the likelihood of disease transmission. Now, we extend our initial model to incorporate a control function u(t) that represent the effect of vector control through the application of chemicals to water. The extended model takes the form(39)S˙1t=Aπ1-βtWtS1t-μ+a1S1t+a2S2t+pκσIt,S˙2t=Aπ2-βtWt1-θS2t-μ+a2S2t+a1S1t+pκ1-σIt,E˙t=βtWtS1t+1-θS2t-γt+μEt,I˙t=γtEt-μ+κIt,W˙t=αtIt-ut+δtWt.A successful control strategy is one that reduces the numbers of exposed and infected individuals over a finite time horizon[0,T] at minimal cost. Our objective functional is therefore formulated as (40)Jut=∫0TEt+It+Butdtsubject to the constraints of the ordinary differential equations in system (39) and the coefficient B (positive) represents a weight in the cost of the control. This objective functional and the differential equations are linear in the control with bounded states, and one can demonstrate by standard results that an optimal control and corresponding optimal states exist [21]. Although majority of applications of optimal control theory to infectious disease control consider controls in a quadratic form due to a number of mathematical advantages (e.g., if the control set is a compact and convex polyhedron it follows that the Hamiltonian attains its minimum over the control set at a unique point), here we considered a linear control since it is regarded as a more realistic function to use in a biological framework compared to the quadratic [22–24].By using Pontryagin’s maximum principle [21, 25] we derive necessary conditions for our optimal control and corresponding states. The Hamiltonian is(41)H=Et+It+But+λ1Aπ1-βtWtS1t-μ+a1S1t+a2S2t+pκσIt+λ2Aπ2-βtWt1-θS2t-μ+a2S2t+a1S1t+pκ1-σIt+λ3βtWtS1t+1-θS2t-μ+γEt+λ4γEt-μ+κIt+λ5αIt-ut+δWt.Given an optimal controlu(t), there exist adjoint functions, λ1,λ2,λ3,λ4,λ5 corresponding to the states S1,S2,E,I, and W such that(42)λ˙1=λ1βtW+μ+a1-λ2a1-λ3βtW,λ˙2=λ2βtW1-θ+μ+a2-λ1a2-λ3βtW1-θ,λ˙3=-1+λ3μ+γ-λ4γt,λ˙4=-1-λ1pκσ-λ2pκ1-σ+λ4μ+κ-λ5α,λ˙5=λ1βtS1+λ2βtS21-θ-λ3βtS1+1-θS2+λ5ut+δ,whereλ1(T)=0, λ2(T)=0, λ3(T)=0, λ4(T)=0, and λ5(T)=0 are the transversality conditions.The Hamiltonian is minimized with respect to the control variableu∗. Since the Hamiltonian is linear in the control, we must consider if the optimal control is bang-bang (at it is lower or upper bound), singular, or a combination. The singular case could occur if the slope or the switching function(43)∂H∂u=B-λ5Wis zero on nontrivial interval of time. Note that the optimal control would be at its upper bound or its lower bound according to(44)∂H∂u<0or>0.To investigate the singular case, let us suppose∂H/∂u=0 on some nontrivial interval. In this case we calculate(45)ddt∂H∂u=0and then we will show that control is not present in that equation. To solve for the value of the singular control, we will further calculate(46)d2dt2∂H∂u=0.We simplify the time derivative of∂H/∂u,(47)0=ddt∂H∂u=ddtB-λ5W=-λ5W′-λ5′W.We calculate both sums separately and add them together. The first sum can be written as(48)-λ5W′=-λ5αI-ut+δW=-αλ5I+ut+δλ5W.The second sum can be written as(49)-λ5′W=λ1-βtS1+λ2-βtS21-θ+λ3βtS1+1-θS2+λ5-ut-δW=-βtλ1S1W-βt1-θλ2S2W+βtλ3S1W+βt1-θλ3S2W-ut+δλ5W.Thus combining, we have(50)0=ddt∂H∂u=-αλ5I-βλ1S1W-β1-θλ2S2W+βλ3S1W+β1-θλ3S2W=-αλ5I+βλ3-λ1S1W+β1-θλ3-λ2S2W.We see that the control does not explicitly show in this expression, so next we calculate the second derivative with respect to time.(51)0=d2dt2∂H∂u=-αλ5′I+-αλ5I′+βλ3′-λ1′S1W+βλ3-λ1S1W′+S1′W+β1-θλ3′-λ2′S2W+β1-θλ3-λ2S2W′+S2′W.Using systems (3) and (41), we simplify (47) as follows:(52)0=d2dt2∂H∂u=α2βtλ3-λ1S1+λ3-λ21-θS2+μ+κλ5-δλ5I-αγλ5E+βtAλ3-λ1π1+λ3-λ21-θπ2+pκλ3-λ1σ+λ3-λ21-θ1-σIW-αλ5I+βλ3-λ1S1W+βλ3-λ21-θS2Wut+βt-1+γλ3-λ4-δλ3-λ1+a1λ2-λ3+λ3-λ21-θS1W+βt1-θ-1+γλ3-λ4-δλ3-λ2+a2λ1-λ31-θ+λ3-λ1S2W.The above equation can be written in the form(53)d2dt2∂H∂u=Φ1tut+Φ2t=0.and we can solve for the singular control as(54)usingulart=-Φ2tΦ1tif(55)Φ1t≠0,a≤-Φ2tΦ1t≤bwith(56)Φ1t=-αλ5I+βtλ3-λ1S1W+βtλ3-λ21-θS2W=-αBWI-βtλ3-λ1S1W-βtλ3-λ21-θS2Wand(57)Φ2t=α2βtλ3-λ1S1+λ3-λ21-θS2+μ+κλ5-δλ5I-αγBWE+βtAπ1λ3-λ1+λ3-λ21-θπ2W+βtpκσλ3-λ1+λ3-λ21-θ1-σIW+βt-1+γλ3-λ4-δλ3-λ1+a1λ2-λ3+a1λ3-λ21-θS1W+βt1-θ-1+γλ3-λ4-δλ3-λ2+a2λ1-λ31-θ+λ3-λ1S2W.To check the generalized Legendre-Clebsch condition for the singular control to be optimal, we required/dud2/dt2(∂H/∂u)=Φ1(t) to be negative [26, 27]. To summarize, our control characterization is as follows: on a nontrivial interval,(58)if∂H∂u<0attthenu∗t=bif∂H∂u>0attthenu∗t=aif∂H∂u<0attthenusingulart=-Φ2Φ1Hence, our control is optimal att provided Φ1(t)<0 and a≤-Φ2(t)/Φ1(t)≤b.Using parameter values from Table1 we investigate the effects of vector control on Guinea worm eradication in South Sudan. We solved the optimality system numerically using forward-backward sweep method [28]. Starting with an initial guess for the control, the state system is solved forward in time. Using those new state values, the adjoint system is solved backward in time. The control is updated using a convex combination of the old control values and the new control values from the characterization. The iterative method is repeated until convergence.Figure3 shows the concentration of the pathogen in the environment as a function of time, in both the presence and absence of the control. It is clear that the presence of optimal control can lead to eradication of the pathogen from the environment. More precisely, we note that, in the presence of vector control, the concentration of the pathogen within the environment will be reduced to a level close to 0 when t≥10 months. We also note that in the absence of the control the bacteria population has an oscillatory pattern and the amplitude of the oscillations appears to be constant.Figure 3
The concentration of the pathogen in the environment with and without the optimal control.Figure4 illustrates the numbers of exposed and infectious individuals over time with and without optimal control. The results clearly depict that the optimal control strategy significantly reduces the exposed and infectious populations (compared to the case without control), to a level close to 0 when t>30 months. More precisely, the number of exposed and infectious individuals will be reduced to a level close to 0 when t>30 and t>60, respectively. Meanwhile, we observe that the numbers of the pathogen in the environment and the infected population will not converge to zero at the same period. This is due to the long incubation period associated with Guinea worm disease.Figure 4
The numbers of exposed and infectious human population with and without the optimal control.
(a) (b)Figure5 shows the optimal control profile as a function of time, with B=2×10-5. As we can observe, the control profile starts from the maximum (u=0.8) initially and stays there for a period of 60 months and then it switches to its minimum (u=0.2) where it remains till the final time. Figure 6 shows the optimal control profile as a function of time, with high costs (B=2×105). It is evident that when the cost of the control is high, the control will have to be implemented at its maximum for a short period of time compared to a scenario when the costs are low.Figure 5
The control profile.Figure 6
The control profile with high costs.


## 3. Discussion and Conclusions

Guinea worm disease is parasitic waterborne disease that is on the verge of eradication. The numbers of individuals afflicted by this disease has declined from 3.5 million cases in 1986 to only 25 cases in 2016, thanks to the Global Guinea worm eradication program [29].In this paper, we have proposed, analyzed, and simulated a new mathematical model for Guinea worm disease that incorporates indirect disease transmission, low and high risk susceptible population, and seasonality. Our model is motivated by the fact that observed incidence of Guinea worm disease exhibits strong seasonal fluctuations in the Sahelian zone and the humid and forest zone [5], with morbidity burden of the disease concentrated in a few months each year. In addition, we are also motivated by the fact that in countries where the disease is still a menace (South Sudan, Chad, Mali, and Ethiopia) access to safe water remains a formidable challenge [1].To account for seasonal fluctuations we modeled the transmission rate, incubation rate, pathogen shedding rate, and pathogen decay rate as periodic functions. To include in heterogeneity on disease transmission we subdivided the susceptible population into two compartments, the low and high risk. We computed the reproduction numberR0 and demonstrated that when R0<1 the model is globally asymptotically stable. We also demonstrated that for R0>1 the disease persists and there exists at least one positive periodic solution. This implies that the disease can be eradicated from the community whenever the basic reproduction number is less than unity; otherwise the disease persists.Using the Fourier series analysis method we fitted our proposed model to the reported data on symptomatic cases of Guinea worm disease in South Sudan in order to estimate the periodic function for the disease transmission rate, incubation rate, pathogen shedding rate, and pathogen decay rate. Meanwhile, we performed an optimal control study to assess the implications of vector control on Guinea worm disease eradication. Our optimal control aims to minimize the numbers of latently infected and infections individuals at minimal costs. Our results demonstrate that optimal control has the potential to ensure successful and timely elimination of the disease in South Sudan. Further, we also observed that with low costs the control will be carried out at maximum for a longer period of time compared to when the costs are high.This study has some limitations. In our mathematical model we did not incorporate the life cycle of the parasite, a factor which can significantly enhance our understanding of Guinea worm disease dynamics. In future we plan to extend this work to incorporate this aspect.



---

*Source: 5084687-2018-07-05.xml*