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, [69]). 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



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 [2224].

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