Research Article  Open Access
Aristide G. Lambura, Gasper G. Mwanga, Livingstone Luboobi, Dmitry Kuznetsov, "Mathematical Model for Optimal Control of SoilTransmitted Helminth Infection", Computational and Mathematical Methods in Medicine, vol. 2020, Article ID 6721919, 15 pages, 2020. https://doi.org/10.1155/2020/6721919
Mathematical Model for Optimal Control of SoilTransmitted Helminth Infection
Abstract
In this paper, we study the dynamics of soiltransmitted helminth infection. We formulate and analyse a deterministic compartmental model using nonlinear differential equations. The basic reproduction number is obtained and both diseasefree and endemic equilibrium points are shown to be asymptotically stable under given threshold conditions. The model may exhibit backward bifurcation for some parameter values, and the sensitivity indices of the basic reproduction number with respect to the parameters are determined. We extend the model to include control measures for eradication of the infection from the community. Pontryagian’s maximum principle is used to formulate the optimal control problem using three control strategies, namely, health education through provision of educational materials, educational messages to improve the awareness of the susceptible population, and treatment by mass drug administration that target the entire population(preschool and schoolaged children) and sanitation through provision of clean water and personal hygiene. Numerical simulations were done using MATLAB and graphical results are displayed. The cost effectiveness of the control measures were done using incremental costeffective ratio, and results reveal that the combination of health education and sanitation is the best strategy to combat the helminth infection. Therefore, in order to completely eradicate soiltransmitted helminths, we advise investment efforts on health education and sanitation controls.
1. Introduction
The worldwide burden of helminth (worm) infection cannot be underestimated. It is estimated that over 1.5 billion people worldwide are infected with helminths [1]. Soiltransmitted helminths are among the NTDs which mostly affect the poorest and the resourceconstrained populations of the world. Soiltransmitted helminths are caused by interstinal parasitic nematode hookworm, Ascaries lumbriocoides and Trichuris trichiura species through ingesting eggs in unwashed undercooked vegetables and unpeeled fruits in contaminated water sources or ingestion of eggs by children who play in contaminated environment. More than 267 million preschoolaged children and more than 568 million schoolaged children (SAC) live in endemic areas where these parasites are highly transmitted. This establishes chronic infections over time which hinders and impairs the developmental and cognitive growth of both preschool and schoolaged children which in turn affect their academic performance. It also has social and economic deficits [2]. The possibility of elimination of the helminth infection is by reducing the number of worms/parasites in a host, but this cannot be achieved without application of other means of reducing the density of parasites from the contaminated environment. Hence proper treatment of the helminth infection and effective control strategy that target this group and infected adults is needed. Currently, there are control programmes to combat the disease which include regular treatment by mass drug administration (MDA), water sanitation, and health education. To reduce the morbidity of the soiltransmitted helminth, antihelminthic drugs like albendazole and benzimidazole have been used for decades but still there has not been an optimal application of the existing control measures with these treatment alternatives in endemic areas [3, 4].
Over years, mathematical models have helped to improve our understanding of population dynamics and provide tools for the assessment and evaluation of a number of control programmes in place. The mathematical model for the helminth infection can be traced back to the papers by [4–6]. Deterministic mathematical models inform policymakers to realize the effort required to increase control coverage for the achievement of less than one percent prevalence and intensity of soiltransmitted helminths by 2020 [7]. These models have been used intensively but none have applied optimal control for the soiltransmitted helminth infection.
However, optimization of the existing interventions tools for the helminth infection has been the priority research areas to examine their impact and sustainability [8]. The study on the implementation of these control measures and how to deliver them optimally is of great importance. Thus, in this study, we consider optimal control of the helminth mathematical model with the preventive measures (health education) to sensitize the susceptible population and treatment by mass drug administration and sanitation. The preventive measures include provision of education materials and educational messages and sanitation include proper provision of clean water, personal hygiene, and installation of public toilets while MDA is administered to the entire population (preSAC and SAC), but for the healthy individuals (susceptible and recovered), this treatment procedure becomes placebo.
2. Materials and Methods
2.1. Model Formulation
For modeling, the total human population is divided into four subpopulations , , , and . represents the number of individuals who are not infected but can be infected by helminth parasites. represents the number of individuals exposed to helminth infection but not infectious. represents the number of individuals infested with parasitic worms. represents the number of individuals recovered from the helminth infection. represents the parasite population in the environment. Individuals in the susceptible class contract infection after a sufficient contact with the contaminated environment. All the newborns enter a susceptible class at the rate and are moved to the exposed class at the rate . Individuals in the exposed class are being reinfected as they interact with contaminated environment but remain latent for a period of , where is the progression rate from the exposed to infective class. Again, we assume that an infected individual contaminate the environment at rate but do not acquire new infection; hence, individuals in the infected class are moved to the recovered class as the result of natural immunity at the rate . Also, we assume that the recovered individuals are aware of infection so that individuals in this class may not interact with contaminated environment; furthermore, they have temporary immunity against the infection. Hence, individuals in the recovered class return to the susceptible class at a rate due to loss of immunity and awareness.
The infected individuals contaminate the environment at the rate as they defecate outside the latrines to increase the concentration of parasitic worms in the soil. The parasitic worms in the soil infect human population at the rate where is the intake rate of eggs in contaminated food or larvae that has penetrated the skin to lead the transmission of infection and is the number (density) of the helminths at which the infection rate is half the maximum rate. The natural mortality in each human class is denoted by , and the infectious individuals have added the helminthinduced death rate denoted by . Further, we assume that the parasitic worms die naturaly at the rate . The aforementioned assumptions and explanation give the dynamics of the helminth model as illustrated in Figure 1. Model parameters are fully described in Table 1. The model has five nonlinear ordinary differential equations given by system (2). with the condition

From equation (2), we have
Hence the dynamics of the helminth model can be investigated by normalizing model (2).
Let , , , , and .
Then, equation (7) can be written in terms of fractional variables as
Using equation (8) and differentiating the fractional variables with respect to leads to the following normalized model:
The solutions of system (9) enter the region that is positively invariant given by
The following theorem guarantees the well posedness of system (9) such that the solutions with nonnegative initial conditions will remain nonnegative for all future time.
Theorem 1. The region is positively invariant with respect to system (9) and nonnegative solutions exist for all future time.
Proof. Given the nonnegative initial conditions, it can be shown that the solutions to system (9), are nonnegative for otherwise we assume a contradiction that there exists a first time such that , and , , , , for ; or there exists such that , and , , , , for ; or there exists such that , and , , , , for ; or there exists such that , and , , , , for ; or there exists such that , and , , , , for as in [9].
In the first case, consider the first equation of system (9) which is a contradiction implying that .
In the second case, we have if this is true both and . It follows that which is a contradiction implying that .
In the third case, we have which is a contradiction implying that .
In the fourth case, we have which is a contradiction implying that .
In the fifth case, we have which is a contradiction implying that .
Hence, in all cases, , , , , and are all positive for all . Thus, the components of the solutions to system (9) remain positive for all .
Thus, the helminth model is well posed epidemiologically and mathematically. Thus, it is sufficient to consider the dynamics of (9) in [10].
3. Model Analysis
The qualitative analysis of the normalized model (9) is considered in order to get an understanding of the dynamics of the helminth infection.
3.1. The DiseaseFree Equilibrium and Basic Reproduction Number
The diseasefree equilibrium of model (9) is obtained by equating all the derivatives to zero and then setting , , , and . Thus, we obtain the DFE as
The basic reproduction number denoted by is the quantity that quantifies the ability of the disease to invade the community, and it is the quantity that governs the spread of the disease. Furthermore, it helps in deciding on control strategies to be adopted for disease eradication. Therefore, in our case, determines the role of parasitic worms/larvae present in the soil to produce secondary infections to the individual hosts.
The nextgeneration approach is used to obtain following Diekmann et al. [11]. Using the principles of a nextgeneration matrix, let be the appearance of new infections and transfer on individuals between compartments, respectively. Applying the linearization technique, the Jacobian matrices of and at DFE are
The basic reproduction number is the spectral radius of the nextgeneration matrix . Thus,
From (23), the fraction represents the average number of susceptible individuals being infected during the infectious period, and the fraction represents the proportion of individuals that survives the latent period, while represents the fraction of parasites diminished from the environment.
3.2. Stability of DFE
We state Theorem 2 for the stability of the DFE.
Theorem 2. The diseasefree equilibrium of model (9) is LAS if and unstable if .
Proof. To investigate the local stability of the DFE at , we obtain the Jacobian matrix of equation (9), i.e.,
Clearly, and are the first and second eigenvalues of the Jacobian matrix which are strictly negative. The remaining three can be obtained by considering the submatrix :
From equation (25), the characteristic equation is given by where
Applying the RouthHurwitz criteria [12], the conditions in (27) , , and . Here, , , and when , and then . Hence, according to the RouthHurwitz criteria, the submatrix has negative real parts whenever .
Thus, the DFE, is locally asymptotically stable if .
3.3. Global Stability of the DiseaseFree Equilibrium
For the proof of global stability, we consider the theorem developed by CastilloChavez et al. [13] and later applied in [14, 15].
To apply the theorem, system (9) is written in the form where denotes the number of uninfected components (individuals) and denotes infected components (individuals) including the latent and infectious individuals.
The diseasefree equilibrium is now denoted by .
The conditions (H1) and (H2) in system (32) must be satisfied to guarantee local asymptotic stability: where is the Metzelermatrix (the offdiagonal elements of are nonnegative) and is the region where the model makes biological sense.
If system (30) satisfies the two conditions in (32), then the following theorem holds.
Theorem 3. The disease fixed point is globally asymptotically stable of system (30) provided that and that conditions H1 and H2 are satisfied.
Proof. From system (9) we have, Then, system (33) can be reduced to
Now, condition H1 is satisfied as it can be observed from the solution, namely, . As , the solution implying global convergence of solution of (35) in . Thus, condition H1 is satisfied.
Let
Then,
From above, if and only if . Thus, may not be globally asymptotically stable for some parameter values. This suggests the possibility of backward bifurcation of Section 3.5.
3.4. Existence of Endemic Equilibrium
The endemic equilibrium of model (9) is defined as the steady state solution which occurs when disease persists in the community and is obtained by equating the derivatives of system (9) to zero. If we denote the force of infection as at steady states and let for any choice of at the endemic equilibrium; then, the endemic equilibrium can be expressed in terms of the force of infection as where
Using equation five of (38) and the force of infection with the model parameters, then (38) satisfies the following polynomial: where
In equation (44), corresponds to the disease free equilibrium and the polynomial corresponds to the existence of endemic equilibrium which coexist with disease free equilibrium point when .
3.5. Bifurcation Analysis
To investigate the type of bifurcation model (9) exhibits, we use the center manifold theory of CastilloChavez and Song [16]. To apply the center manifold theory we make the change of variables for the normalized helminth model (9). Let , , , and . Using the vector notation , then (9) can be written in the form with . If we choose as the bifurcation parameter and solving for we have: where
Now, the matrix at DFE of the equations (46) is
At the value , the Jacobian has a simple zero eigenvalue whose corresponding left and right eigenvectors are expressed as
Using the result in [16], we need to compute the coefficients and where
Now, evaluating the partial derivatives of system (46) at , we obtain
Furthermore,
Therefore, where
The coefficient is positive. By CastilloChavez and Song [16], coefficient decides the local dynamics of . Therefore, if , then system (9) will exhibit forward bifurcation, and if , it undergoes backward bifurcation.
3.6. Sensitivity Analysis of the Model Parameters
Numerical sensitivity analysis is done by computing sensitivity indices of basic reproduction number which measures the model robustness to parameter values [17]. The forward sensitivity index of a variable to a parameter is a ratio of the relative change in the variable to the relative change in the parameter.
Definition 4. The normalized forward sensitivity index of a variable that depends differentiable on a parameter is defined as .
From equation (23), we compute the sensitivity indices of each parameter involved in using the equation . For example, the sensitivity index of parameter value with respect to is given by and other indices are , , , , , and which were obtained and evaluated using parameter values in Table 1 as illustrated in Figure 2.
To determine the most sensitive parameters for the dynamics of the model, sensitivity analysis was carried out to determine the model robustness to parameter values. Clearly , , and have positive indices with most sensitive parameter be and the least positively sensitive parameter . This imply that the endemicity of the disease increases when the parameters , , and are increased while keeping other parameters constant. On the other hand, the parameters , , , and have negative indices and when each of these are decreased, while keeping the other parameters constantly decreases the value of implying that they decrease the endemicity of the disease. Thus, based on these results, we suggest the following interventions to be made for the eradication of the helminth infection from the endemic communities. The first intervention is the mass cleanliness to reduce the concentration of the parasites from the contaminated environment. The other one is the use of anthelmintic drugs including MDA for both healthy and infected individuals which will in turn reduce the shedding rate of the parasite to the environment.
3.7. Extension of the Model to Optimal Control
In this section, we extend the heminth model (2) by incorporating three timedependent control measures. The best intervention strategies will be identified from combination of these control measures so as to minimize the infection from the community. The optimal control model for helminth infection is formulated by considering the following controls: (1)Health education control that is aimed at sensitizing the susceptible population. This is done through provision of education materials (TV messages, radio, posters, and leaflets), educational messages which can be delivered by health practitioners or teachers in schools, education that improve hygiene awareness and behavior of people who defecate outside the latrines, and purchase of shoes to both preschoolaged and schoolaged children [18]. Thus. the force of infection is modified as where is the effectiveness measure of health education control, where such that if , then health education is not effective and corresponds to completely effective health education, implies that health education is effective to some degree(2)The use of mass drug administration (deworming) , which is administered regularly without individual diagnosis is applied to the entire population [19]. Application of this control moves the individuals in exposed class to the susceptible class and those in the infected class to the recovered class(3)Sanitation control which measures the efforts to prevent susceptible from contracting the disease. This include drinking clean water, personal hygiene, avoiding eating raw vegetables, cleaning of fruits, and intensive cleanliness of the environment which helps to break the helminth transmission cycle (see [19]). Thus, the clearance of the parasite from the environment will be accelerated by , where is the effectiveness of the control to protect susceptible population from contracting helminth infection
After incorporating the controls and using the same parameters and variables as in (2), we obtain the modified model with controls for helminth infection given by with , , , , and . In this study, we assume that the control functions , , and are Lebesque integrable functions with for to mean that when the control is zero, there is no any control implemented and when the control is one there is maximum implementation of the control. The objective is to obtain the optimal levels of the control and state variables that optimize the objective functional given by
Here, we want to minimize the number of infected individual from the population while keeping the cost of controls low. In equation (60), is the relative weight of the infected individuals that accounts for the social cost while is the individual cost for MDA of the entire population(preschool and schoolaged children) because it the most vulnerable population for helminth infection. Furthermore, is the relative costs weight associated with background costs for each control measure such as advisement costs, production of posters, ordering and transportation of drugs and water treatment and is the final fixed time. In this work, the cost is not linear as may result to bang bang which means that the optimal control take values at only upper and lower control set; thus, we choose to have quadratic cost on the control of the form (see [20]). The problem becomes of determining the optimal triplet ( ) such that where is the admissible set.
3.7.1. Existence of an Optimal Control
The solution of the helminth model (9) is proved to be bounded so we use this result to prove the existence of the optimal control. However, existence of optimal control is shown using the approach of Fleming and Richel [14] and later applied in [15, 21–23]. For detailed proof, (see [14], Theorem 4.1, p 6869).
3.7.2. The Hamiltonian and Optimality Systems
To obtain the Hamiltonian , we apply “Pontryagin’s maximum principle” [24] to derive the necessary conditions for the optimal pair. Therefore, the Hamiltonian is expressed as with denotes the adjoint variables to be determined by taking the negative derivative of the Hamiltnian function.
Theorem 5. Given the optimal set , , and that minimizes over , there exist an adjoint variables such that with transversality conditions, The characterization of the control set (, , and ) using the technique of control bounds gives where
Proof. By using ([24]), the adjoint system is the EulerLagrange equations obtained taking the negative of the partial derivative of the Hamiltonian in (62) with respect to the state variables. Thus, the adjoint system can be written as
Using the same principle, we get the controls by solving the equation obtained by taking the partial derivative of the Hamiltonian with respect to each control. For example,
Now, setting at for , we obtain the controls set:
The controls can then be written using standard control arguments with bounds on the controls as
In compact notation, we writ, where
The optimality system consists of the state system (55) coupled with adjoint system (63) together with the characterization of the optimal control with the initial and transversality conditions.
The state equations in (55) and the costate (adjoint) equation in system (63) are bounded and satisfies the Lipchitz condition, thus uniqueness of the optimality systems follows using the approaches in [25, 26].
3.8. Numerical Simulation
We perform numerical simulation of the optimal control model using the parameter values given in Table 1. Various optimal control strategies are applied to control the burden of helminth infection in the community. Starting with an initial guess for the optimal controls , , and , state system (55) is solved numerically forward in time using the fourthorder RungeKutta method with initial conditions , , , , and . Then, co state (adjoint) system (63) is solved numerically backward in time using the fourthorder RungeKutta method with the state variables and initial guess of the controls obtained earlier. All controls are updated using convex combination of the previous controls and its characterization then the process is repeated until convergence is achieved. The following cost weights were used in the simulation: , , , , and with the efficacy rates , , and .
3.8.1. Control with Health Education and MDA
In this strategy, we used the combination of health education and MDA interventions to optimize the objective functional while sanitation control is set to zero. Figure 3(a) shows that with this strategy there is a significant effect in reducing the number of exposed individuals compared with the case when there is no control but does not bring the exposed individuals to zero at the final control period. The helminthinfected individuals in Figure 3(b) decreases to almost zero when the optimal combination of health education and MDA are in place compared to when there is no control. In Figure 3(c), we observed that the parasite population in the environment is reduced but does not reach zero at the final control period. This suggests that using combination of control measures, the infection can be lowered in the community but not completely removed in the specified control period. The control profile in Figure 3(d) suggests that control combination with health education and MDA should be maximized and applied throughout the intervention period.
(a)
(b)
(c)
(d)
3.8.2. Control with MDA and Sanitation
Next, we consider the combination of MDA and sanitation interventions while the health education is set to zero. Figures 4(a) and 4(b) show a similar trend as in Figures 3(a) and 3(b), respectively, but with this strategy more effective than the combination of health education and MDA. This is due to the fact that with MDA, anthelminthic drugs are regularly administered to the entire population (preschool and schoolaged children) while using sanitation clears the density of parasites from the environment. The control profile in Figure 4 suggests the control with MDA should be seized after 20 days while sanitation control should be kept at a maximum for approximately 20 days and then decreased to zero. Therefore, these results show that combination of these two controls measures can as well manage to control helminth infection in the community.
(a)
(b)
(c)
(d)
3.8.3. Control with Health Education and Sanitation
We next consider the control with combination of health education and sanitation in the absence of MDA. The health education and sanitation were used to optimize the objective functional while MDA control was set to zero. Figure 5(a) shows that the number of exposed individuals were reduced and were totally controlled after 50 days. Figure 5(b) shows that the number of infected individuals were reduced and totally controlled at the final period. Figure 5(c) shows that the parasite population is reduced to the minimum level. The control profile in Figure 5(d) suggests that the control with health education should be maximized for 70 days while sanitation should be maximized for 65 days.
(a)
(b)
(c)
(d)
3.8.4. Control with Health Education, MDA, and Sanitation
Finally, we consider the case where all controls were in place. In this strategy the controls’ health education, MDA, and sanitation were used to optimize the objective functional . Figure 6(a) shows that the number of exposed individuals is totally controlled at the final control period. Figure 6(b) shows that the number of helminthinfected individuals is also totally controlled at the final control period. Figure 6(c) shows that the parasite population is more reduced to zero at the final control time. Figure 6(d) suggests that the health education control should be maximized for 15 days, the MDA control be seized after 15 days, while the sanitation should be maximized for 18 days and reduced to zero due to the fact that the infection has been cleared in the community.
(a)
(b)
(c)
(d)
3.9. CostEffectiveness Analysis
Cost effectiveness analysis is performed to quantify the cost effectiveness of different combinations of control strategies using Incremental Cost Effective Ratio (ICER). In this approach, when comparing two competing intervention strategies, one intervention is compared with the next less effective alternative [28, 30, 31]. The ICER formula is given by
To obtain the total number of infections averted, we take the sum of the difference between the total number of exposed and infected individuals without and with control. Also to obtain the total cost for each strategy, we used the cost function as shown in Table 2.
 
, , 