Printer Friendly

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

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 S(t), exposed E(t), infected individuals displaying clinical signs of the disease 7(f), and the concentration of the bacteria in the environment W(f). 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:

[mathematical expression not reproducible] (1)

where [[beta].sub.0], [[gamma].sub.0], [[alpha].sub.0], and [[delta].sub.0] denote the contact rate, incubation rate, pathogen shedding rate, and pathogen decay rate, respectively, without seasonal forcing and w > 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:

[mathematical expression not reproducible] (2)

All model variables and parameters are considered positive. Further the variables [S.sub.1](t), [S.sub.2](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) = [S.sub.1] (t) + [S.sub.2](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 parameter A is the recruitment rate, [[pi].sub.1] and [[pi].sub.2] represent the fraction of new recruits into classes [S.sub.1] and [S.sub.2], respectively, [mu] is the natural mortality rate, [a.sub.1] is the transfer from susceptible high risk to susceptible low risk, [a.sub.2] is the transfer from susceptible low risk to susceptible high risk, [theta] is a modification factor that accounts for the impact of access of safe drinking water on disease transmission, and [k.sup.-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 (l - 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 a of recovered individuals who move become susceptible to infection join [S.sub.1] (t) and the complementary (1 - [sigma]) join [S.sub.2](t).

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 R(t) described by differential equation [??](t) = (1- p)kl(t) - [mu]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:

[mathematical expression not reproducible] (3)

System (3) has a disease-free equilibrium given by [E.sup.0] = (S01, S [degrees],0,0,0), with

[mathematical expression not reproducible] (4)

The severity of a disease can be measured by the basic reproduction number, [R.sub.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 F(t)) and the transfer terms (denoted by V(t)) at disease-free equilibrium are given by

[mathematical expression not reproducible] (5)

Thus, the basic reproduction number of the time-averaged autonomous systems is

[mathematical expression not reproducible] (6)

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

[mathematical expression not reproducible] (7)

where Y(t,s), t [greater than or equal to] s, is the evolution operator of the linear w-periodic system dy/dt = -V(t)y and [phi](t), the initial distribution of infectious humans, is w-periodic and always positive. The effective reproduction number for a periodic model is then determined by calculating the spectral radius ofthe next infection operator,

[R.sub.0] = p(L). (8)

Through direct calculation, the evolution operator Y(t, s), for model (3), is given by

[mathematical expression not reproducible] (9)


[mathematical expression not reproducible] (10)


[mathematical expression not reproducible] (11)

We can numerically evaluate the next infection operator by

[mathematical expression not reproducible] (12)


[mathematical expression not reproducible] (13)

[mathematical expression not reproducible] (14)

2.3. Threshold Dynamics. In this section we will use the basic reproduction number [R.sup.0] defined in (8) to establish the threshold results, presented in Theorem 1, for model (3). To that end, we first note that R+ is positively invariant for the following cooperative system:

[mathematical expression not reproducible] (15)

and that ([S.sup.0.sub.1], [S.sup.0.sub.2]) is the unique equilibrium solution which is globally attractive in [R.sup.2.sub.+].

Theorem 1.

(i) If [R.sup.0] < 1, then the disease-free equilibrium [E.sub.0] of system (3) is globally asymptotically stable.

(ii) If [R.sup.0] > 1, then system (3) admits at least one positive w-periodic solution, and solutions of system (3) are uniformly persistent.

Proof. If ([S.sub.1],(t), [S.sub.2](t),E(t),I(t),W(t)) is a nonnegative solution of (3), then we have

[mathematical expression not reproducible] (16)

Note that any nonnegative solution ([S.sub.1] (t), [S.sub.2](t)) of system (3) approaches ([S.sup.0.sub.1], [S.sup.0.sub.2]) as t [right arrow] [infinity]. It then follows from the standard comparison theorem (see, e.g., [17, Theorem A.4]) that, for any [epsilon] > 0, there is a T > 0 such that

[S.sub.i] (t) < [S.sup.0.sub.i] + [member of], i = -1,2, for t > T. (17)

Thus, for t > T, we have

[mathematical expression not reproducible] (18)


[mathematical expression not reproducible] (19)

By [16, Thorem 2.2], we have [mathematical expression not reproducible] is the spectral radius of [[phi].sub.F-V] ([omega]) and [[phi].sub.F-V] ([omega]) is the monodromy matrix of the linear w-periodic system dy/dt = (F - V)y. Then we can set e sufficiently small such that [mathematical expression not reproducible]. As a consequence, the trivial solution (0,0, 0) of the following linear w-periodic system

[mathematical expression not reproducible] (20)

is globally asymptotically stable. Again by the comparison theorem, we know that E(t) [right arrow] 0, I(t) [right arrow] 0 and W(t) [right arrow] 0 as t [right arrow] [infinity]. Finally, from the first and second equations of system (3) it follows that [S.sub.i](t) [right arrow] [S.sup.0.sub.i] as t [right arrow] [infinity] for i - 1,2. This proves the result in part (i) of Theorem 1.

Next we shall focus on the case [R.sup.0] > 1. We define [mathematical expression not reproducible]. It can easily be verified that both X and [X.sub.0] are positively invariant. Let P : [R.sup.5.sub.+] [right arrow] [R.sup.5.sub.+] be the Poincare map associated with system (3); that is, P([x.sub.0]) = u(w, [x.sub.0]) for all [x.sub.0] [member of] [R.sup.5.sub.+], where u(t, [x.sub.0]) is the unique solution of (3) with u(0, [x.sub.0]) = [x.sub.0]. Set

[mathematical expression not reproducible] (21)

We first demonstrate that

[M.sub.[partial derivative]] = M. (22)

It is evident that M [subset or equal to] [M.sub.[partial derivative]]. For any ([S.sub.1](0), [S.sub.2](0), E(0), 1(0), W(0)) [member of] [partial derivative][X.sub.0] \M, if E(0) = I(0) = 0 and W(0) = 0, then [mathematical expression not reproducible]. By the positive invariance of [X.sub.0], we know that [mathematical expression not reproducible], and thus (22) holds.

Now consider the fixed point [mathematical expression not reproducible] of the Poincare map P. Define [mathematical expression not reproducible]. We show that

[W.sup.S] ([M.sub.0]) [intersection] [X.sub.0] = [theta]. (23)

Based on the continuity of solutions with respect to the initial conditions, for any [member of] >0, there exists [??] > 0 small enough such that for all ([S.sub.1](0), [S.sub.2](0),E(0),I(0),W(0)) [member of] [X.sub.0] with

[mathematical expression not reproducible], (24)

we have

[mathematical expression not reproducible] (25)

To obtain (23), we claim that

[mathematical expression not reproducible] (26)

We prove this claim by contradiction; that is, we suppose [mathematical expression not reproducible]. Without loss of generality, we assume that [mathematical expression not reproducible]. Thus,

[mathematical expression not reproducible] (27)

Moreover, for any t [greater than or equal to] 0, we write t = [t.sub.0] + nw with [t.sub.0] [member of] [0, w) and n = [t/w], the greatest integer less than or equal to t/w. Then we obtain

[mathematical expression not reproducible] (28)

for any t [greater than or equal to] 0. Let [mathematical expression not reproducible]. It follows that

[mathematical expression not reproducible] (29)

[mathematical expression not reproducible] (30)

Hence we obtain

[mathematical expression not reproducible] (31)

where F - V is given by (5) and

[mathematical expression not reproducible] (32)

Again based on [16, Theorem 2.2], [R.sup.0] > 1 if and only if [rho]([[PHI].sub.F-V]([omega])) > 1. Thus, for e small enough, we have [mathematical expression not reproducible] 1 which immediately yields the contradiction as

[mathematical expression not reproducible] (33)

Let [P.sub.1] : [R.sup.2.sub.+] [right arrow] [R.sup.2.sub.+] be the Poincaree map associated with (15). Then ([S.sup.0.sub.1], [S.sup.0.sub.2]) is globally attractive in [R.sup.2.sub.+] \ {0} for [P.sub.1]. It follows that [M.sub.0] is isolated invariant set in X, and notice that [W.sup.S] ([M.sub.0]) [intersection] [X.sub.0] = [theta]. Hence, every orbit in [M.sub.[partial derivative]] converges to [M.sub.0] and [M.sub.0] is acyclic in [M.sub.[partial derivative]]. By [18, Theorem 1.3.1], for a stronger repelling property of [partial derivative][X.sub.0], we conclude that P is uniformly persistent with respect to ([X.sub.0], [partial derivative][X.sub.0]), which implies the uniform persistence of the solutions of system (3) with respect to ([X.sub.0], [partial derivative] [X.sub.0]) [18, Theorem 3.1.1]. Consequently, based on [18, Theorem 1.3.6], the Poincare: map P has a fixed point [mathematical expression not reproducible] is a positive w-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 [[pi].sub.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 [[pi].sub.1] = 0.9 and it follows that [[pi].sub.2] = 0.1.

(ii) The natural mortality rate [mu] = 1/life expectancy: the united nations children's fund report of 2013 [11] highlighted that life expectancy in South Sudan is 55 years; thus [mu] = 0.0015 [month.sup.-1].

(iii) The mean infectious period [k.sup.-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 k = 8.3 x [10.sup.-2] [month.sup.-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 [S.sub.1] to [S.sub.2]; thus we set [a.sub.1] = 9.0 x [10.sup.-6] [month.sup.-1]. In contrast, due to civil unrest which characterize South Sudan population we assume that the rate of transfer from low risk is higher ([a.sub.1] < [a.sub.2]); thus we set [a.sub.2] = 0.009 [month.sup.-1].

Table 1 presents the description of model parameters and their baseline values.

Now, we need to estimate the periodic functions [beta](t), [alpha] (t), [gamma](t), and [delta](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

g(t) = [gamma](t) E(t). (34)

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

[mathematical expression not reproducible], (35)

where [c.sub.j] and [d.sub.j] (j = 1,2,...,8) represent coefficients to be determined and L is the fundamental frequency. We used MATLAB software to determine the coefficients [c.sub.j] and [d.sub.j] (j = 1,2,..., 8) and we obtained

[mathematical expression not reproducible] (36)

with L = 4[pi]/7 and to be precise and exact [omega] = 7/2. The comparison of the data with curve of g(t) is shown in Figure 2.

Next, we proceed to use (36) to estimate the periodic functions [beta](t) (disease transmission rate), [gamma](t) (incubation rate), [alpha](t) (pathogen shedding rate), and [delta](t) (pathogen decay rate). First, we assume the initial population levels as follows: [S.sub.1](0) = S2(0) = 10,000, W(0) = 5000. Further, we set E(0) = 1(0) = 192 based on the number of cases reported in January 2007. Let

[mathematical expression not reproducible] (37)

where [mathematical expression not reproducible] are related to the magnitude of seasonal fluctuations and [[[beta].sub.1] (t), [[gamma].sub.1](t), [[alpha].sub.1] (t), and [[delta].sub.1] (t) are periodic functions to be determined. After simulations and comparisons we obtained

[mathematical expression not reproducible] (38)

with [mathematical expression not reproducible]. Using direct calculation and (37) one can easily derive the expressions for [beta](1), [gamma](l), [alpha](l), and [delta](l).

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 m(1) that represent the effect of vector control through the application of chemicals to water. The extended model takes the form

[mathematical expression not reproducible] (39).

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

J [m (1)] = [[integra].sup.T.sub.0] [E (1) + I(t) + Bu (1)] dt (40)

subject 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

[mathematical expression not reproducible] (41)

Given an optimal control u(t), there exist adjoint functions, [[lambda].sub.1], [[lambda].sub.2], [[lambda].sub.3], [[lambda].sub.4], [[lambda].sub.5] corresponding to the states [S.sub.1], [S.sub.2], B, /, and W such that

[mathematical expression not reproducible] (42)

where [[lambda].sub.1](T) = 0, [[lambda].sub.2] (T) = 0, [[lambda].sub.3] (T) = 0, [[lambda].sub.4] (T) = 0, and [[lambda].sub.5] (T) = 0 are the transversality conditions.

The Hamiltonian is minimized with respect to the control variable u*. 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

[partial derivative]H/[partial derivative]u = B - [[lambda].sub.5]W (43)

is zero on nontrivial interval of time. Note that the optimal control would be at its upper bound or its lower bound according to

[partial derivative]H/[partial derivative]u < 0 or > 0. (44)

To investigate the singular case, let us suppose [partial derivative]H/[partial derivative]m = 0 on some nontrivial interval. In this case we calculate

d/dt ([partial derivative]H/[partial derivative]u) = 0 (45)

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

[d.sup.2]/[dt.sup.2] ([partial derivative]H/[partial derivative]u) = 0. (46)

We simplify the time derivative of [partial derivative]H/[partial derivative]u,

[mathematical expression not reproducible] (47)

We calculate both sums separately and add them together. The first sum can be written as

[mathematical expression not reproducible] (48)

The second sum can be written as

[mathematical expression not reproducible] (49)

Thus combining, we have

[mathematical expression not reproducible] (50)

We see that the control does not explicitly show in this expression, so next we calculate the second derivative with respect to time.

[mathematical expression not reproducible] (51)

Using systems (3) and (41), we simplify (47) as follows:

[mathematical expression not reproducible] (52)

The above equation can be written in the form

[mathematical expression not reproducible] (53)

and we can solve for the singular control as

[mathematical expression not reproducible] (54)


[mathematical expression not reproducible] (55)


[mathematical expression not reproducible] (56)


[mathematical expression not reproducible] (57)

To check the generalized Legendre-Clebsch condition for the singular control to be optimal, we require (d/du)([d.sup.2]/[dt.sup.2]) ([partial derivative]H/[partial derivative]u) = [[phi].sub.1] (t) to be negative [26, 27]. To summarize, our control characterization is as follows: on a nontrivial interval,

[mathematical expression not reproducible] (58)

Hence, our control is optimal at t provided [[phi].sub.1] (f) < 0 and [mathematical expression not reproducible].

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 t [greater than or equal to] 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 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 5 shows the optimal control profile as a function of time, with B = 2 x [10.sup.-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 x [10.sup.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.

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

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.


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.


[1] B. H. Beyene, A. Bekele, A. Shifara, and Y. Ebsite, "Elimination of Guinea worm in Ethopia; current status of the diseases eradictaion strategies and challenges to the end game. Ethiop Md J," in current status of the disease's eradictaion strategies and challenges to the end game. Ethiop Md J, vol. 55, supplement 1, pp. 15-31, 55,15-31, 2017.

[2] World Health Organization, World Health Assembly, Resolution WHA 57.9. Elimination of Dracunculiasis: resolution of the 57th World Health Assembly, Geneva, Switzerland, 2016.

[3] E. Ruiz-Tiben and D. R. Hopkins, "Dracunculiasis (Guinea Worm Disease) Eradication," Advances in Parasitology, vol. 61, pp. 275-309, 2006.

[4] O. O. Alew and N. Peters, Ethiopian Dracunculiasis Eradication Program (EDEP) GOG WOREDA Meeting, National Review, Gambella, Ethopia, 2015.

[5] K. A. Mitra and R. A. Mawson, "Neglected tropical diseases: epidemiology and global burden," Trop. Med. Infect. Dis, vol. 2, p. 36, 2017.

[6] J. R. Smith, P Cloutier, J. Harrison, and A. Desforges, A Mathematical Model for the Eradication of Guinea Worm Disease. Understanding the Dynamics of Emerging and Re-Emerging Infectious Diseases Using Mathematical Models, 2012.

[7] R. Netshikweta and W. Garira, "A multiscale model for the world's first parasitic disease targeted for eradication: guinea worm disease," Comput Math Methods Med, vol. 2017, 2017.

[8] I. A. Adetunde, "The Epidemiology of guinea worm infection in Tamale District, in the Northern Region of Ghana," Journal of Modern Mathematics and Statistics, pp. 50-54, 2008.

[9] L. Kathryh, "Guinea Worm Disease (Dracunculiasis): Opening a mathematical can of worms," Tech. Rep., Bryn Mawr College, 2012.

[10] World Health Organization, Public Health Risk Assessment and Interventions - Conflict and Humanitarian Crisis, Republic of South Sudan, 2014.

[11] United Nations Children's Fund (UNICEF), State of the World'S Children, 2013.

[12] G. Biswas, D. P. Sankara, J. Agua-Agum, and A. Maiga, "Dracunculiasis (guinea worm disease): eradication without a drug or a vaccine.," Philosophical Transactions of the Royal Society B: Biological Sciences, vol. 368, no. 1623, p. 20120146, 2013.

[13] O. C. Collins and K. S. Govinder, "Incorporating heterogeneity into the transmission dynamics of a waterborne disease model," Journal of Theoretical Biology, vol. 356, pp. 133-143, 2014.

[14] S. L. Robertson, M. C. Eisenberg, and J. H. Tien, "Heterogeneity in multiple transmission pathways: Modelling the spread of cholera and other waterborne disease in networks with a common water source," Journal of Biological Dynamics, vol. 7, no. 1, pp. 254-275, 2013.

[15] P van den Driessche and J. Watmough, "Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission," Mathematical Biosciences, vol. 180, pp. 29-48, 2002.

[16] W. Wang and X. Zhao, "Threshold dynamics for compartmental epidemic models in periodic environments," Journal of Dynamics and Differential Equations, vol. 20, no. 3, pp. 699-717, 2008.

[17] H. L. Smith and P Waltman, The Theory of the Chemostat, Cambridge University Press, 1995.

[18] X. Q. Zhao, Dynamical Systems in Population Biology, Springer, New York, NY, USA, 2003.

[19] L. Liu, X.-Q. Zhao, and Y. Zhou, "A tuberculosis model with seasonality," Bulletin of Mathematical Biology, vol. 72, no. 4, pp. 931-952, 2010.

[20] S. Bowong and J. Kurths, "Modeling and analysis of the transmission dynamics of tuberculosis without and with seasonality," Nonlinear Dynamics, vol. 67, no. 3, pp. 2027-2051, 2012.

[21] W. Fleming and R. Rishel, Deterministic and Stochastic Optimal Control, Springer New York, New York, NY, 1975.

[22] C. J. Silva and D. F. Torres, "A TB-HIV/AIDS coinfection model and optimal control treatment," Discrete and Continuous Dynamical Systems - Series A,vol. 35, no. 9,pp. 4639-4663,2015.

[23] H. Schattler, U. Ledzewicz, and H. Maurer, "Sufficient conditions for strong local optimality in optimal control problems with l2 -type objectives and control constraints, discrete contin," Discrete and Continuous Dynamical Systems - Series B, vol. 19, no. 8, pp. 2657-2679, 2014.

[24] C. Silva, H. Maurer, and D. Torres, "Optimal control of a tuberculosis model with state and control delays," Mathematical Biosciences and Engineering, vol. 14, no. 1, pp. 321-337, 2017

[25] L. S. Pontryagin, V. G. Boltyanskii, R. V. Gamkrelidze, and E. F. Mishchenko, The Mathematical Theory of Optimal Processes, Wiley, New York, NY, USA, 1962.

[26] A. J. Krener, "The high order maximal principle and its application to singular extremals," SIAM Journal on Control and Optimization, vol. 15, no. 2, pp. 256-293,1977

[27] H. R. Joshi, S. Lenhart, S. Hota, and F. Agusto, "Optimal control of an SIR model with changing behavior through an education campaign," Electronic Journal of Differential Equations, No. 50, 14 pages, 2015.

[28] S. M. Lenhart and J. T. Workman, Optimal Control Applied to Biological Models, CRC Press, 2007

[29] D. Molyneux and D. P Sankara, "Guinea worm eradication: Progress and challenges--should we beware of the dog?" PLOS Neglected Tropical Diseases, vol. 11, no. 4, Article ID e0005495, 2017.

Anthony A. E. Losio (1,2) and Steady Mushayabasa (01)

(1) University of Zimbabwe, Department of Mathematics, P.O. Box MP167, Harare, Zimbabwe

(2) University of Juba, Department of Mathematics, P.O. Box 82, Juba, Central Equatoria, South Sudan

Correspondence should be addressed to Steady Mushayabasa;

Received 12 February 2018; Revised 18 May 2018; Accepted 11 June 2018; Published 5 July 2018

Academic Editor: Keshlan S. Govinder

Caption: Figure 1: Diagram of the model structure.

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

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

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

Caption: Figure 5: The control profile.

Caption: Figure 6: The control profile with high costs.
TABLE 1: Description of model parameters for system (3).

Symbol                        Definition                    Value

p                       Proportion of recovered
                     who become susceptible to the           0.3

[sigma]                 Proportion of recovered
                        who are at high risk to              0.9

[a.sub.1]           Rate of transfer from high risk
                       to low risk susceptible        9.0 x [10.sup.-6]

[a.sub.1]           Rate of transfer from low risk
                       to high risk susceptible             0.009

[[pi].sub.1]          Proportion of new births in
                        high risk susceptible                0.9

[[pi].sub.1]          Proportion of new births in
                    low risk susceptible population          0.1

[mu]                  Natural mortality rate for           0.0015
[[kappa].sup.-l]       Average infectious period      8.3 x [10.sup.-2]

[theta]                   Modification factor               0.98

A                          Human birth rate                  100

[[beta].sub.0]       Averaged environment-to-host
                             transmission             1.2 x [10.sup.-5]

[[gamma].sub.0]        Averaged incubation rate       1.0 x [10.sup.-1]

[[alpha].sub.0]     Averaged parasite shedding rate    2 x [10.sup.-1]

[[delta].sub.0]           Parasite death rate               2.16

Symbol                         Unit              Source


                                --               Assumed


                                --                [10]

                          [month.sup.-1]         Assumed

                          [month.sup.-1]         Assumed

                                --                [10]

                                --                [10]

[mu]                      [month.sup.-1]          [11]

[[kappa].sup.-l]              month              Assumed

[theta]                         --               Assumed

A                   Individuals [month.sup.-1]   Assumed

                      Larvae [month.sup.-1]      Fitting

[[gamma].sub.0]               month               [12]

[[alpha].sub.0]       Larvae [month.sup.-1]      Fitting

[[delta].sub.0]           [month.sup.-1]           [6]
COPYRIGHT 2018 Hindawi Limited
No portion of this article can be reproduced without the express written permission from the copyright holder.
Copyright 2018 Gale, Cengage Learning. All rights reserved.

Article Details
Printer friendly Cite/link Email Feedback
Title Annotation:Research Article
Author:Losio, Anthony A.E.; Mushayabasa, Steady
Publication:Journal of Applied Mathematics
Article Type:Report
Date:Jan 1, 2018
Previous Article:Analytical Synthesis of Regulators for Nonlinear Systems with a Terminal State Method on Examples of Motion Control of a Wheeled Robot and a Vessel.
Next Article:Solution of Quadratic Programming with Interval Variables Using a Two-Level Programming Approach.

Terms of use | Privacy policy | Copyright © 2022 Farlex, Inc. | Feedback | For webmasters |