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.
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 Section 3.
2. Methods and Results
2.1. 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 , , and denote the contact rate, incubation rate, pathogen shedding rate, and pathogen decay rate, respectively, without seasonal forcing and 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 , , , , 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 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, and represent the fraction of new recruits into classes and , respectively, is the natural mortality rate, is the transfer from susceptible high risk to susceptible low risk, 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 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 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 and the complementary join
Figure 1 shows a flowchart depicting the dynamics of Guinea worm disease in this study.
2.2. The Reproduction Number
Now we explore the dynamical behavior of Guinea worm disease. Since the variable described by differential equation 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 , with
The severity of a disease can be measured by the basic reproduction number, , 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 disease-free equilibrium are given by
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 bywhere
2.3. Threshold Dynamics
In this section we will use the basic reproduction number defined in (8) to establish the threshold results, presented in Theorem 1, for model (3). To that end, we first note that is positively invariant for the following cooperative system:and that is the unique equilibrium solution which is globally attractive in .
Theorem 1. (i)If , then the disease-free equilibrium of system (3) is globally asymptotically stable.(ii)If , then system (3) admits at least one positive -periodic solution, and solutions of system (3) are uniformly persistent.
Proof. If is a nonnegative solution of (3), then we haveNote that any nonnegative solution of system (3) approaches as . It then follows from the standard comparison theorem (see, e.g., [17, Theorem A.4]) that, for any , there is a such thatThus, for , we haveDefine By [16, Thorem 2.2], we have , where is the spectral radius of and is the monodromy matrix of the linear -periodic system . Then we can set sufficiently small such that . As a consequence, the trivial solution of the following linear -periodic systemis globally asymptotically stable. Again by the comparison theorem, we know that and as . Finally, from the first and second equations of system (3) it follows that as for . This proves the result in part (i) of Theorem 1.
Next we shall focus on the case . We define . It can easily be verified that both and are positively invariant. Let be the Poincaré map associated with system (3); that is, for all , where is the unique solution of (3) with . Set We first demonstrate thatIt is evident that . For any , if and , then . If and , then and If and then and It follows that for . By the positive invariance of , we know that for ; hence , and thus (22) holds.
Now consider the fixed point of the Poincaré map . Define . We show thatBased on the continuity of solutions with respect to the initial conditions, for any , there exists small enough such that for all with we haveTo obtain (23), we claim thatWe prove this claim by contradiction; that is, we suppose for some . Without loss of generality, we assume that . Thus,Moreover, for any , we write with and , the greatest integer less than or equal to . Then we obtainfor any . Let . It follows thatHence we obtainwhere is given by (5) andAgain based on [16, Theorem 2.2], if and only if . Thus, for small enough, we have which immediately yields the contradiction as Let be the Poincareé map associated with (15). Then is globally attractive in for . It follows that is isolated invariant set in , and notice that . Hence, every orbit in converges to and is acyclic in . By [18, Theorem 1.3.1], for a stronger repelling property of , we conclude that is uniformly persistent with respect to , which implies the uniform persistence of the solutions of system (3) with respect to [18, Theorem 3.1.1]. Consequently, based on [18, Theorem 1.3.6], the Poincaré map has a fixed point , with . Thus, and , 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 : according to the world health organization report of 2014, more than 90% of the people in South Sudan lives on less than US1 a day [10]. Based on this assertion we assume that and it follows that (ii)The natural mortality rate : the united nations children’s fund report of 2013 [11] highlighted that life expectancy in South Sudan is 55 years; thus (iii)The mean infectious period 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 .(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 to ; thus we set In contrast, due to civil unrest which characterize South Sudan population we assume that the rate of transfer from low risk is higher ; thus we set .
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
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 and represent coefficients to be determined and is the fundamental frequency. We used MATLAB software to determine the coefficients and and we obtained
with and to be precise and exact The comparison of the data with curve of is shown in Figure 2.
Next, we proceed to use (36) to estimate the periodic functions (disease transmission rate), (incubation rate), (pathogen shedding rate), and (pathogen decay rate). First, we assume the initial population levels as follows: , Further, we set based on the number of cases reported in January 2007. Let
where , , , and are related to the magnitude of seasonal fluctuations and , , , and are periodic functions to be determined. After simulations and comparisons we obtained
with , , , and Using direct calculation and (37) one can easily derive the expressions for , , , and
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 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 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–24].
By using Pontryagin’s maximum principle [21, 25] we derive necessary conditions for our optimal control and corresponding states. The Hamiltonian is
Given an optimal control , there exist adjoint functions, corresponding to the states , and such that
where , , , , and are the transversality conditions.
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 on some nontrivial interval. In this case we calculate
and 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
We simplify the time derivative of ,
We calculate both sums separately and add them together. The first sum can be written as
The second sum can be written as
Thus combining, we have
We see that the control does not explicitly show in this expression, so next we calculate the second derivative with respect to time.
Using systems (3) and (41), we simplify (47) as follows:
The above equation can be written in the form
and we can solve for the singular control as
if
with
and
To check the generalized Legendre-Clebsch condition for the singular control to be optimal, we require to be negative [26, 27]. To summarize, our control characterization is as follows: on a nontrivial interval,
Hence, our control is optimal at provided and .
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 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 months. More precisely, the number of exposed and infectious individuals will be reduced to a level close to 0 when and , 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.
(a)
(b)
Figure 5 shows the optimal control profile as a function of time, with . As we can observe, the control profile starts from the maximum () initially and stays there for a period of 60 months and then it switches to its minimum () where it remains till the final time. Figure 6 shows the optimal control profile as a function of time, with high costs (). 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.
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 number and demonstrated that when the model is globally asymptotically stable. We also demonstrated that for 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.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors have no conflicts of interest.
Acknowledgments
The authors are grateful to the handling editor Dr. Keshlan S. Govinder for valuable comments and suggestions that led to an improvement of this paper.