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

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 onGuinea worm disease eradication.


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

Model
Framework.The model proposed in [6] subdivides the host population into categories of susceptible (), exposed (), infected individuals displaying clinical signs of the disease (), and the concentration of the bacteria in the environment ().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: 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: ( All model variables and parameters are considered positive.Further the variables  1 (),  2 (), (), (), and () represent the susceptible high risk, susceptible low risk, exposed, infectious, and recovered individuals at time , respectively, such that the total human population is given by () =  1 () +  2 () + () + () + ().Meanwhile () represents the population of the pathogen in the environment.Model parameters defined in (1) retain the same definitions.
Model parameter  is the recruitment rate,  1 and  2 represent the fraction of new recruits into classes  1 and  2 , respectively,  is the natural mortality rate,  1 is the transfer from susceptible high risk to susceptible low risk,  2 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,  is the proportion of recovered individuals who become susceptible to infection again based on their behavior and the complementary part (1 − ) 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  1 () and the complementary (1 − ) join  2 ().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: ( System (3) has a disease-free equilibrium given by E 0 = ( 0 1 ,  0 2 , 0, 0, 0), with The severity of a disease can be measured by the basic reproduction number, R 0 , 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 ()) and the transfer terms (denoted by ()) at diseasefree equilibrium are given by ) . (5) Thus, the basic reproduction number of the time-averaged autonomous systems is 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 where (, ),  ≥ , is the evolution operator of the linear -periodic system / = −() and (), 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, Through direct calculation, the evolution operator (, ), for model (3), is given by where with We can numerically evaluate the next infection operator by where 2.3.Threshold Dynamics.In this section we will use the basic reproduction number R 0 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: and that ( 0 1 ,  0 2 ) is the unique equilibrium solution which is globally attractive in R 2 + .
(i) If  0 < 1, then the disease-free equilibrium E 0 of system ( 3) is globally asymptotically stable.
(ii) If  0 > 1, then system (3) admits at least one positive -periodic solution, and solutions of system ( 3) are uniformly persistent.

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/life expectancy: the united nations children's fund report of 2013 [11] highlighted that life expectancy in South Sudan is 55 years; thus  = 0.0015 month −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 −2 month −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  1 to  2 ; thus we set  1 = 9.0 × 10 −6 month −1 .In contrast, due to civil unrest which characterize South Sudan population we assume that the rate of transfer from low risk is higher ( 1 <  2 ); thus we set  2 = 0.009 month −1 .
Table 1 presents the description of model parameters and their baseline values.
Now, we need to estimate the periodic functions (), (), (), and ().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) Further, since variables and the periodic parameters in model system (3) are continuous functions of , we will make use of trigonometric functions to fit ().Define where with  = 4/7 and to be precise and exact  = 7/2.The comparison of the data with curve of () is shown in Figure 2.

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 () that represent the effect of vector control through the application of chemicals to water.The extended model takes the form A successful control strategy is one that reduces the numbers of exposed and infected individuals over a finite time horizon [0, ] at minimal cost.Our objective functional is therefore formulated as subject to the constraints of the ordinary differential equations in system (39) and the coefficient  (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][23][24].By using Pontryagin's maximum principle [21,25] we derive necessary conditions for our optimal control and corresponding states.
The Hamiltonian is minimized with respect to the control variable  * .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 is zero on nontrivial interval of time.Note that the optimal control would be at its upper bound or its lower bound according to To investigate the singular case, let us suppose / = 0 on some nontrivial interval.In this case we calculate To check the generalized Legendre-Clebsch condition for the singular control to be optimal, we require (/)( 2 / 2 )(/) = Φ 1 () to be negative [26,27].To summarize, our control characterization is as follows: on a nontrivial interval, Hence, our control is optimal at  provided Φ 1 () < 0 and  ≤ −Φ 2 ()/Φ 1 () ≤ .
Using parameter values from Table 1 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.
Figure 3 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  ≥ 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 4 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  > 30 months.More precisely, the number of exposed and infectious individuals will be reduced to a level close to 0 when  > 30 and  > 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 5 shows the optimal control profile as a function of time, with  = 2×10 −5 .As we can observe, the control profile starts from the maximum ( = 0.8) initially and stays there for a period of 60 months and then it switches to its minimum ( = 0.2) where it remains till the final time.Figure 6 shows the optimal control profile as a function of time, with high costs ( = 2 × 10 5 ).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.

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 number R 0 and demonstrated that when R 0 < 1 the model is globally asymptotically stable.We also demonstrated that for R 0 > 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.

Figure 1 :
Figure 1: Diagram of the model structure.

Figure 1
Figure1shows a flowchart depicting the dynamics of Guinea worm disease in this study.

Figure 2 :
Figure 2: The monthly numbers of new Guinea worm cases and its fitted curve.

Figure 3 :
Figure 3: The concentration of the pathogen in the environment with and without the optimal control.

Figure 4 :
Figure 4: The numbers of exposed and infectious human population with and without the optimal control.

Figure 6 :
Figure 6: The control profile with high costs.

Table 1 :
Description of model parameters for system (3).