Mathematical Model for Optimal Control of Soil-Transmitted Helminth Infection

In this paper, we study the dynamics of soil-transmitted helminth infection. We formulate and analyse a deterministic compartmental model using nonlinear differential equations. The basic reproduction number is obtained and both disease-free 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 school-aged 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 cost-effective 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 soil-transmitted helminths, we advise investment efforts on health education and sanitation controls.


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 resource-constrained populations of the world. Soil-transmitted 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 preschool-aged children and more than 568 million school-aged 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 school-aged 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 soil-transmitted 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][5][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 soil-transmitted helminths by 2020 [7]. These models have been used intensively but none have applied optimal control for the soil-transmitted 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 (pre-SAC and SAC), but for the healthy individuals (susceptible and recovered), this treatment procedure becomes placebo.

Model Formulation.
For modeling, the total human population is divided into four subpopulations SðtÞ, EðtÞ, IðtÞ, and RðtÞ. SðtÞ represents the number of individuals who are not infected but can be infected by helminth parasites. EðtÞ represents the number of individuals exposed to helminth infection but not infectious. IðtÞ represents the number of individuals infested with parasitic worms. RðtÞ represents the number of individuals recovered from the helminth infection. MðtÞ 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 b 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 1/ρ, 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 q. 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 K 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 helminth-induced death rate denoted by d. Further, we assume that the parasitic worms die naturaly at the rate μ M . 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 S + E + I + R = N: From equation (2), we have Hence the dynamics of the helminth model can be investigated by normalizing model (2).
Let s = S/N, e = E/N, i = I/N, r = R/N, and m = M/K. Then, equation (7) can be written in terms of fractional variables as Computational and Mathematical Methods in Medicine Using equation (8) and differentiating the fractional variables with respect to t 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 Ω ∈ R 5 + is positively invariant with respect to system (9) and nonnegative solutions exist for all future time.
In the first case, consider the first equation of system (9) which is a contradiction implying that sðtÞ > 0.
In the second case, we have if this is true both sðt 1 Þ ≥ 0 and mðt 1 Þ ≥ 0. It follows that eð t 1 Þ ≥ 0 which is a contradiction implying that eðtÞ > 0.
In the third case, we have which is a contradiction implying that iðtÞ > 0.
In the fourth case, we have which is a contradiction implying that rðtÞ > 0.
In the fifth case, we have which is a contradiction implying that mðtÞ > 0.
Hence, in all cases, sðtÞ, eðtÞ, iðtÞ, rðtÞ, and m are all positive for all t > 0. Thus, the components of the solutions to system (9) remain positive for all t > 0.
Thus, the helminth model is well posed epidemiologically and mathematically. Thus, it is sufficient to consider the dynamics of (9) in Ω [10].

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.

The Disease-Free Equilibrium and Basic Reproduction
Number. The disease-free equilibrium of model (9) is obtained by equating all the derivatives to zero and then setting e = 0, i = 0, r = 0, and m = 0. Thus, we obtain the DFE as The basic reproduction number denoted by R 0 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,   The next-generation approach is used to obtain R 0 following Diekmann et al. [11]. Using the principles of a next-generation matrix, let be the appearance of new infections and transfer on individuals between compartments, respectively. Applying the linearization technique, the Jacobian matrices of F and V at DFE E 0 are The basic reproduction number R 0 is the spectral radius of the next-generation matrix FV −1 . Thus, From (23), the fraction β/ðb + d + q + εÞ represents the average number of susceptible individuals being infected during the infectious period, and the fraction ρ/ðb + ρÞ represents the proportion of individuals that survives the latent period, while ε/μ M represents the fraction of parasites diminished from the environment.

Stability of DFE.
We state Theorem 2 for the stability of the DFE.
Proof. To investigate the local stability of the DFE at E 0 , we obtain the Jacobian matrix of equation (9), i.e., Clearly, λ 1 = −b and λ 2 = −ðb + γÞ are the first and second eigenvalues of the Jacobian matrix J E 0 which are strictly Immunity waning rate for recovered individuals γ 1 7 day −1 [27] Clearance rate for parasitic worms μ M 13 37500 day −1 [27] Shading rate of parasites to the environment ε 0.09 day −1 Assumed Progression rate from exposed to infected class ρ 1 10 day −1 [27] 4 Computational and Mathematical Methods in Medicine negative. The remaining three can be obtained by considering the submatrix J 1 E 0 : From equation (25), the characteristic equation is given by where Applying the Routh-Hurwitz criteria [12], the conditions in (27) Thus, the DFE, E 0 is locally asymptotically stable if R 0 < 1.

Global Stability of the Disease-Free Equilibrium.
For the proof of global stability, we consider the theorem developed by Castillo-Chavez et al. [13] and later applied in [14,15].
To apply the theorem, system (9) is written in the form where X = ðs, rÞ ∈ ℝ 2 + denotes the number of uninfected components (individuals) and Z = ðe, i, mÞ ∈ ℝ 3 + denotes infected components (individuals) including the latent and infectious individuals.
The conditions (H1) and (H2) in system (32) must be satisfied to guarantee local asymptotic stability: where A = D Z GðX * , 0Þ is the Metzeler-matrix (the offdiagonal elements of A 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 E 0 = ðX * , 0Þ is globally asymptotically stable of system (30) provided that R 0 < 1 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, sðtÞ = 1 + ðsð0Þ − 1Þe −bt . As t ⟶ ∞, the solution sðtÞ ⟶ 1 implying global convergence of solution of (35) in Ω. Thus, condition H1 is satisfied. Let Then, From above,ĜðX, ZÞ ≥ 0 if and only if d = ε = 0. Thus, E 0 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 E * 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 λ * = βm/1 + m at steady states and let z * = b − ðd + εÞi * for any choice of 5 Computational and Mathematical Methods in Medicine i * at the endemic equilibrium; then, the endemic equilibrium can be expressed in terms of the force of infection as Using equation five of (38) and the force of infection λ * = βm * /ð1 + m * Þ with the model parameters, then (38) satisfies the following polynomial: where In equation (44), λ * = 0 corresponds to the disease free equilibrium and the polynomial Pðλ * Þ = 0 corresponds to the existence of endemic equilibrium which co-exist with disease free equilibrium point when R 0 < 1.

Bifurcation Analysis.
To investigate the type of bifurcation model (9) exhibits, we use the center manifold theory of Castillo-Chavez and Song [16]. To apply the center manifold theory we make the change of variables for the normalized helminth model (9). Let x 1 = s, x 2 = e, x 3 = i, x 4 = r, and If we choose β = β * as the bifurcation parameter and solving for R 0 = 1 we have: where Now, the matrix W at DFE E 0 of the equations (46) is Computational and Mathematical Methods in Medicine At the value R 0 = 1, the Jacobian W 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 a and b where Now, evaluating the partial derivatives of system (46) at ðE 0 , β * Þ, we obtain Furthermore, Therefore, where The coefficient b is positive. By Castillo-Chavez and Song [16], coefficient a decides the local dynamics of E 0 . Therefore, if a 0 < 0, then system (9) will exhibit forward bifurcation, and if a 0 > 0, it undergoes backward bifurcation.

Sensitivity Analysis of the Model Parameters.
Numerical sensitivity analysis is done by computing sensitivity indices of basic reproduction number R 0 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.
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 μ M , b, d, and q have negative indices and when each of these are decreased, while keeping the other parameters constantly decreases the value of R 0 implying that they decrease 7 Computational and Mathematical Methods in Medicine 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.

Extension of the Model to Optimal Control.
In this section, we extend the heminth model (2) by incorporating three time-dependent 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 u 1 ðtÞ 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 preschool-aged and school-aged children [18]. Thus. the force of infection is modified as λ = ðð1 − κ 1 u 1 ÞβMÞ/ðK + MÞ where κ 1 is the effectiveness measure of health education control, where κ 1 ∈ ð0, 1Þ such that if κ 1 = 0, then health education is not effective and κ 1 = 1 corresponds to completely effective health education, 0 < κ 1 < 1 implies that health education is effective to some degree (2) The use of mass drug administration (deworming) κ 2 u 2 ðtÞ, which is administered regularly without individual diagnosis is applied to the entire population [19]. Application of this control moves the indi-viduals in exposed class to the susceptible class and those in the infected class to the recovered class (3) Sanitation control ðu 3 ðtÞÞ 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 μ M + κ 3 u 1 , where κ 3 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 Sð0Þ = S 0 , Eð0Þ = E 0 , Ið0Þ = I 0 , Rð0Þ = R 0 , and Mð0Þ = M 0 . In this study, we assume that the control functions u 1 ðtÞ, u 2 ðtÞ, and u 3 ðtÞ are Lebesque integrable functions with 0 ≤ u i ≤ 1 for i = 1, 2, 3 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), A 1 is the relative weight of the infected individuals that accounts for the social cost while A 2 u 2 is the individual cost for MDA of the entire population(preschool-and school-aged children) because it the most vulnerable population for helminth infection. Furthermore, w i is the relative costs weight associated with background costs for each control measure u i such as advisement costs, production of posters, ordering and transportation of drugs and water treatment and t f is the final fixed time. In this work, the cost is not linear as  Computational and Mathematical Methods in Medicine 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 ð1/2Þw i u 2 i (see [20]). The problem becomes of determining the optimal triplet (u * 1 ,u * 2 ,u * 3 ) such that where U = fðu 1 , u 2 , u 3 Þ | each u i is measurable with 0 ≤ u i ≤ 1 for 0 ≤ t ≤ t f g is the admissible set.

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][22][23]. For detailed proof, (see [14], Theorem 4.1, p 68-69).

The Hamiltonian and Optimality Systems.
To obtain the Hamiltonian H, we apply "Pontryagin's maximum principle" [24] to derive the necessary conditions for the optimal pair. Therefore, the Hamiltonian is expressed as with λ i , i = 1, 2, 3, 4, 5 denotes the adjoint variables to be determined by taking the negative derivative of the Hamiltnian function.
Theorem 5. Given the optimal set u 1 , u 2 , and u 3 that minimizes J over U, there exist an adjoint variables λ 1 , λ 2 , λ 3 , λ 4 , and λ 5 such that with transversality conditions, λ i ðt f Þ = 0, i = 1, 2, 3, 4, 5: The characterization of the control set (u * 1 , u * 2 , and u * 3 ) using the technique of control bounds gives where Proof. By using ( [24]), the adjoint system is the Euler-Lagrange 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 ∂H/∂u i = 0 at u * 1 for i = 1, 2, 3, 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 u 1 ðtÞ, u 2 ðtÞ, and u 3 ðtÞ, state system (55) is solved numerically forward in time using the fourth-order Runge-Kutta method with initial conditions Sð0Þ = 1000, Eð0Þ = 150, Ið0Þ = 50, Rð0Þ = 5, and Mð0Þ = 200. Then, co state (adjoint) system (63) is solved numerically backward in time using the fourth-order Runge-Kutta 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: A 1 = 800, A 2 = 0:15, w 1 = 15, w 2 = 30, and w 3 = 10 with the efficacy rates κ 1 = 0:7, κ 2 = 0:8, and κ 3 = 0:6.
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 J 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 helminth-infected 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.
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 school-aged 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.
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 J 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.
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 J. 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 helminth-infected 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 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.
3.9. Cost-Effectiveness 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 ICER = Difference in costs between strategies Difference in the totalnumber of infection averted : 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.
Comparing strategy 3.8.1 and 3.8.4, strategy 3.8.1 is more costly and less effective than strategy 3.8.4. Hence, we omit strategy 3.8.1 and recompute ICER values as indicated in Table 3.
From Table 3, we observe that strategy 3.8.2 has higher cost compared to strategy 3.8.4, Hence, we omit 3.8.2 and recompute ICER as indicated in Table 4.     Table 4, we conclude that strategy 3.8.3 (control with health education and sanitation) is more cost effective than strategy 3.8.4 (control with health education, MDA and sanitation), implying that it is the cheapest of all the control strategies.

Discussions and Conclusions
In this paper, we formulated and analysed deterministic model for the helminth infection (soil-transmitted helminth).
In Section 3, we analysed the model by obtaining the feasible region and disease-free and endemic equilibrium points with their stability. Moreover, bifurcation analysis were carried out and the model exhibited backward bifurcation for some parameter values, and sensitivity analysis of R 0 with respect to parameters β, ρ, ε, μ M , b, d, q was done and results revealed that β was the most positive sensitive parameter while μ M was the most negative sensitivity parameter. Furthermore, the helminth model developed in Section 2 was extended to include health education, treatment by mass drug administration and sanitation as control measures. We solved the optimal control problem using Pontryagin's maximum principle by first obtaining the Hamiltonian, the adjoint variables, and the characterization of the controls. Applying single-control measure (health education or MDA or sanitation) alone is not sufficient to control the infection; thus, combination of control measures were used to optimize the objective functional J by considering combination of health education with mass drug administration, combination of mass drug administration with sanitation, combination of health education with sanitation, and all the three controls. Finally, we investigated the cost effectiveness of the control strategies to ascertain the least and most expensive strategies. Based on our results, we conclude that the combination of health education and sanitation is the best strategy.

Data Availability
The data are available inside the manuscript and were obtained from related published articles. There is no restriction on the availability of the data.

Conflicts of Interest
The authors declare no conflict of interest.