Mathematical Modeling and Analysis of Transmission Dynamics and Control of Schistosomiasis

Department of Mathematics, Pan African University Institute for Basic Sciences, Technology and Innovation, Box 62000-00200, Nairobi, Kenya Department of Mathematics, School of Arts and Sciences, University of The Gambia, P.O. Box 3530, Serekunda, Gambia Department of Mathematical Sciences, Federal University of Technology Akure, P.M.B. 704, Akure, Ondo State, Nigeria School of Mathematics, College of Biological and Physical Science, University of Nairobi, Box 30197-00100, Nairobi, Kenya


Introduction
Schistosomiasis is a disease that is common in the tropics, and it is caused by a group of parasitic worms known as schistosomes. These collections of parasitic worms are, in some instances, called blood flukes. Several species of schistosomes cause illnesses in human beings. The three major ones are Schistosoma haematobium, Schistosoma mansoni, and Schistosoma japonicum. These species of schistosomes can be found in parts of Africa, Asia, and South America, [1]. Schistosoma mansoni and Schistosoma japonicum mainly cause diseases in the liver and bowels but Schistosoma haematobium mostly affects the urinary and genital areas [2].
The life cycles of schistosomes are very complex and involve fresh water and (mostly) two hosts: a definite human host (also cattle in the case of Schistosoma japonicum) and an intermediate snail host. Schistosomes reproduce both sexually and asexually [2]. The sexual reproduction takes place in the definite hosts, and the asexual reproduction takes place in an intermediate snail host. The life cycle begins in the human portal system. Schistosome eggs are released to freshwater by infected humans through their feces or urine. These eggs hatch into larvae once the eggs enter freshwater. These larvae are known as miracidia. The miracidia, within a day, find and infect the intermediate host, the snails. After a period of development in the snail, the miracidia further mature and transform into cercariae through asexual reproduction and are released into the water in the presence of light. The cercariae (with a life expectancy of up to 2 days) will further swim and penetrate into the skin of a definite host (humans and also cattle in the case of Schistosoma Japonicum). The cercariae find their way to the liver through the lungs and develop into matured adult worms. Thereafter, they move into the blood vessels and the organs around it. These worms can live in these organs for many years. The adult schistosomes can lay up to 300 eggs within a definite host daily, and these eggs can penetrate into the bladder or intestines. Some of these eggs get released during urination or defecation and may find their way back into fresh waters, and their life cycles will begin afresh. However, other eggs will be trapped in some of these organs, and their activities cause damage to the organs. It takes between 4 and 6 weeks for an exposed person to become infectious [2].
The disease is contracted by humans through contact with infested freshwater which is inhabited by infected snails; there is no direct human to human transmission. People with schistosomiasis can suffer from rash or itchy skin, cough, fever, blood in urine or stool, muscle ache, stunted growth, etc. Later in their lives, they can experience malfunctioning of some major organs like the liver, intestines, lungs, and bladder due to damage caused by the eggs trapped in their body. These eggs can be found in the spinal cord and even the brain in some rare cases. It may result in paralysis, inflammation of the spinal cord, or seizures. Children who are reinfected with this disease may experience malnutrition, anemia, and even difficulty of learning [1].
There is still no commercially available vaccine for schistosomiasis. However, praziquantel is the WHO-approved medication for the treatment of the disease. It works by killing the adult worms, hence preventing the production of new eggs in the human host, but does not prevent reinfection as treated individuals can easily get reinfected [3]. Improved sanitation and snail control are also mechanisms used to put the disease under control. Snail control involves the use of molluscicides or biological control methods like use of competitor snails (see [4,5]); however, their effectiveness varies and the disease continues to spread.
Mathematical models can be used as guiding tools to effectively study the spread of diseases, and they have over the years helped both policymakers and public health officials in their decision-making processes with respect to key intervention programs. Its role in the study and understanding of the transmission dynamics as well as the effectiveness of the various control strategies of many infectious diseases has been considered to be very effective. Several mathematical models were designed over the years to study the transmission dynamics and control of the spread of schistosomiasis [4][5][6][7][8][9][10][11].
A review of a number of existing mathematical models of schistosomiasis, their merits, and demerits was carried out by [12]. [13] suggested ways by which controls could be incorporated into mathematical models in order to eradicate or  [27] β h Rate of transmission of humans from susceptible to exposed 0.09753 L cer −1 d −1 [28] β s Rate of transmission of snails from susceptible to exposed 0.615 L mir −1 d −1 [7,29] ξ h Preventive factor due to WASH 0- ε Limitation of the growth velocity 0.2 [7] reduce the prevalence of schistosomiasis in an endemic area. [14] presented some of the measures adopted for the control schistosomiasis in the People's Republic of China. [15] integrated a mathematical model for schistosomiasis to obtain a set of nonlinear integral equations. The Contraction Mapping Principle was then used to establish the existence of a unique solution of the schistosomiasis model. The transmission dynamics of schistosomiasis among humans, cattle, and snails was presented in a mathematical model by [16]. Sensitivity analysis of the model parameters on the basic reproduction number showed that cattle-snail transmission of S. japonicum played an important role in the endemicity of the disease in Hubei Province. [17] performed stability analysis on a schistosomiasis model, and the result indicated that the prevalence of schistosomiasis was affected by the latent period of infection and that effective public health education campaign on drug treatment of the disease could help to curtail it.

Journal of Applied Mathematics
Schistosomiasis is declared by the WHO as one of the most deadly of the neglected tropical diseases. Despite numerous intervention programs to curb the spread of the disease, it still remains endemic in many tropical and subtropical countries particularly in the poverty-stricken African countries [1]. Mathematical models have been used to model the control and spread of many diseases [18,19]. Hence, there is the need to develop a mathematical model that will help to explain the spread and control of schistosomiasis. To the very best of our knowledge, the proposed model is the first to include the treated human class in a schistosomiasis model and has incorporated the exposed period for both human and snail subpopulations alongside a parameter for water, sanitation, and hygiene (WASH). In order to capture the complex nature of the disease dynamics and proffer effective strategies for its control, compartments that have to do with the incubation periods for the hosts that are infected but not releasing pathogen into the environment alongside the treated humans should be incorporated. Consequently, a deterministic model that divides the populations into susceptible humans and snails, infected humans and snails, exposed humans and snails, and treated humans, with the dynamics of free living miracidia and cercariae together with their interactions is proposed.

Model Formulation
Here, a deterministic model that describes the transmission dynamics of schistosomiasis is proposed. We considered   the fact that once a pathogen is introduced within a population, it divides the population into fractions called compartments or components. The human subpopulation is divided into four compartments: susceptible humans ðS h Þ, exposed humans ðE h Þ, infected humans ðI h Þ, and treated humans ð T h Þ; the snail subpopulation into three compartments: susceptible snails ðS s Þ, exposed snails ðE s Þ, and infected snails ð I s Þ; and compartments for the free living miracidia ðN m Þ and free living cercariae ðN c Þ based on the life cycle of the parasite.
The susceptible humans comprise those individuals capable of contracting the disease but have yet to do so. The exposed humans include those individuals who contract cercariae, undergoing and incubation period and have yet to begin producing (miracidia) eggs. The infected humans are those individuals that are symptomatic and have (miracidia) eggs appearing in their feces. The treated humans include individuals who have received treatment. Praziquantel works only by killing the adult worms; thus, only the infected humans are moved to the treated class, and the treated   humans are assumed not to be infectious, i.e., they do not produce eggs for miracidia. The susceptible snails are those snails that have the possibility of getting exposed to miracidia but have yet to do so. The exposed snails are those snails that get infected by miracidia but have yet to shed cercariae. The infected snails are those snails that contract miracidia, survived the exposed period, and are shedding cercariae. Moreover, a parameter for prevention through water sanitation and hygiene (WASH) has been incorporated. Table 1) gives a description of various parameters used in the model, and Figure 1) is the schematic diagram for the proposed model.
The model is derived based on these assumptions: (1) There is no mother-to-child transmission of the disease among humans, and treated humans have no immunity   Journal of Applied Mathematics (2) The rates of depletion of the parasites by humans or snails have no effect on the population dynamics of the parasites (3) Contact with free living cercariae is the only mean through which susceptible humans become exposed to the disease (4) There are no introduction of infectious (humans or snails) into the studied population (5) Treated humans do not produce eggs for the production of miracidia (6) Infected snails do not get treated and, hence, do not recover The model is represented by the following system of ODEs:

Invariant Region and Positivity of Solutions.
Since a population model always requires positive solutions, it is relevant that all solutions be nonnegative. Consider the model as described by the model system equation (1). The total human subpopulation at any time, t, is given by In the absence of mortality due to infection, we have Thus, Similarly, the total snail subpopulation is given by In the absence of mortality due to infection, we get So, N s ≤ α s /μ s . From equation (18) of the state system equation (1), we have Thus, N h , N s , N m , and N c are bounded, and the region

Theorem 1. Suppose that
is the initial condition for the model system equation (1), then the set of solutions is nonnegative ∀t ≥ 0.
Proof. Given that the set of initial conditions is nonnegative, consider the model system equation (1). Then, Separating the variables and integrating, we have Similarly, Due to the complexity of the model system 1, all other state variables at the steady state are expressed in terms of the steady state of I h which is denoted as I * h . Accordingly, 10 Journal of Applied Mathematics and a combination of the above equation (15) and the second equation of the model system 1, the following relation was obtained for I * h : where I * h = 0.

Endemic Equilibrium Point.
At the endemic equilibrium point, all disease states in equation (14) are considered to be positive, and consequently, I * h must be greater than zero for all the other states to be positive. Therefore, the unique positive endemic equilibrium point exists only when I * h in the system equation (16) is positive, which is also in this case the same as R 0 > 1. And thus, it is the point where whenever R 0 > 0.

Stability Analysis
3.3.1. Basic Reproduction Number, R 0 . The basic reproduction number is defined as the average number of secondary infections caused by the emergence of an infectious individual into a completely susceptible population [20]. It is denoted by R 0 . We considered the next generation matrix approach according to Chavez et al. [21] (for details, see Appendix A). The next generation matrix is denoted by G and is obtained as 11 Journal of Applied Mathematics Accordingly, the basic reproduction number, R 0 , which is the dominant eigenvalue of G, is obtained as For simplicity, since R 0 denotes the possible number of uninfected a particular infected individual (in this case either human or snail) can infect, we can rewrite R 0 as represents the number of snails an infected human is capable of infecting through the production of miracidia and R SH 0 = ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi α s β s σ s λ 2 /c 0 μ s μ c ðδ s + μ s Þðσ s + δ s + μ s Þ p denotes the number of humans an infected snail is capable of infecting through the production of cercariae.

Stability of Disease-Free Equilibrium Point
Theorem 2. The disease-free equilibrium point, E 0 , is locally asymptotically stable provided R 0 < 1 and unstable otherwise.
Proof. To verify the validity of the above proposition, see Theorem 2 in Van den Driessche and Watmough [22].
Proof. From the system of equation (1), we have First, suppose that the time, t ⟶ ∞, we need to show that X ⟶ E 0 .
But we know that for FðX, 0Þ, Also, Thus, all points with respect to this condition converge at Now, from the system equation (1), we know that where and GðX, YÞ = AY − G * ðX, YÞ, where Journal of Applied Mathematics Clearly, G * ðX, YÞ ≥ 0, ∀ðX, YÞ ∈ Ω, since Moreover, A is an M-matrix, since the nondiagonal entries are nonnegative [21]. Therefore, the DFE, denoted E 0 , is globally asymptotically stable provided that R 0 < 1.

Stability of the Endemic Equilibrium Point.
We consider the concept of Center Manifold Theory as described in Theorem 4.1 by Chavez and Song [23]. The theorem helps to determine the course of the bifurcation at a critical value of the parameter [24]. The theorem is reproduced in Appendix B.  [23] is locally asymptotically stable for R 0 > 1 (but near 1).
Proof. Now, consider the model system equation (1). Using The role of the parameter ϕ in Theorem 4.1 by Chavez and Song [23] is played by β h with the critical value obtained from R 0 = 1 [7,24]. Hence, β h = β * h , regarded as the bifurcation point, can be obtained by evaluating it at R 0 = 1. So 13 Journal of Applied Mathematics Let the Jacobian of the transformed system equation (33) at E 0 when β h = β * h is denoted by A. Then, it can be shown that 0 is a simple eigenvalue of the matrix A, where and if the vector w = ðw 1 , w 2 , w 3 , w 4 , w 5 , w 6 , w 7 , w 8 , w 9 Þ T is a right eigenvector corresponding to the eigenvalue, 0 of A, then the components of w can be obtained as endemic equilibrium point that is locally asymptotically stable. This is actually a validation of the results from the analytic processes which also confirms that the disease-free equilibrium point is both locally and globally asymptotically stable when R 0 < 1, and when R 0 > 1, a locally asymptotically stable unique endemic equilibrium point exists. ≈ 0:00004379. The human exposed period for schistosomiasis, i.e., the time period of humans from having contact with free living cercariae to the period eggs begin to appear in the human feces ranges from 5 to 6 weeks [2]. Hence, σ h is estimated to lie between 1/6 × 7 ≈ 0:0238 and 1/5 × 7 ≈ 0:0286. The incubation period for snails ranges from 4 weeks to 5 weeks [2]. Thus, σ s lies between 1/5 × 7 ≈ 0:0286 and 1/4 × 7 ≈ 0:0357. The parasite can live within a definite human host for a period ranging from 3 years to 10 years [26]. So, δ h is varied between values corresponding to 3 and 10 years. Consequently, δ h lie between 1/10 × 365 ≈ 0:000274 and 1/3 × 365 ≈ 0:000913. Using ode45 in MATLAB, a numerical simulation of the model equation (1) was carried out with the help of the parameter values presented in Table 1, and the results illustrated graphically show the behavior of the various subpopulations: the human subpopulation (Figure 3(a)), the snail subpopulation (Figure 3(b)), and the subpopulations of the free living miracidia (Figure 3(c)) and cercariae (Figure 3(d)), against time (days).

Numerical Results
It can be observed from Figure 3(a) that the susceptible humans have risen within the first few days but gradually fall afterwards as they continuously come in contact with free living cercariae and later became stable. We experienced a steady and then a quick rise in the exposed and infected human compartments due to the introduction of the disease, both of which later had a slow rise as the disease continue to spread. Equally, the treated human compartment shows a steady rise within the first few days and, afterwards, gains stability. The susceptible snail and the exposed snail compartments, from Figure 3(b), had a little rise within the first few days but become stable shortly and then fall gradually afterwards due to more of the susceptible snails getting in contact with free living miracidia and the exposed snails shedding cercariae. However, there is a high rise in the number of infected snails whose population became stable later on. The infected snails rise above both the susceptible and exposed snails causing the number of susceptible and exposed snails to fall gradually as a result of the fact that infected snails do not recover. Hence, it demonstrates the dominance of the cercariae shedding snails over all other snail compartments. The free living miracidia, from Figure 3(c), also rises within the first few days but reaches a peak and slowly changed course but continues to rise slowly as more infected individuals continue to release eggs for the production of miracidia. The dominance of the cercariae shedding snails explains why more and more cercariae are produced as seen in Figure 3(d). The cercariae subpopulation rises very quickly and gain stability after a number of days, and this can be interpreted to be the result of the dominance of the infected snails over all other snail compartments.

Effects of WASH.
Here, we examined the effectiveness of the application of WASH on the dynamics of the human subpopulation and the dynamics of the free living miracidia and cercariae. It can be observed from the subplots in Figure 4 that WASH as a preventive mechanism in controlling the spread of schistosomiasis is effective if nearly the entire human subpopulation practices it. The results indicate that if more and more individuals within the studied population are encouraged to practice WASH, it will mitigate the number of people contracting the disease, thereby lowering the number of infected individuals and, consequently, reducing the production rate of miracidia. However, it may effectively slow down the rate of production of cercariae, but it will not stop the number of cercariae produced from reaching the point of endemicity as infected snails will continue to produce free living cercariae. Thus, practicing WASH alongside water chlorination (which is considered to be effective in killing cercariae) should be considered.

Effects of Treatment.
Here, we look into how increase in the number of treated infected humans can affect the number of infected humans and also the production rates of the pathogens. It can be observed from Figure 5(a) that increasing the rate of treatment of humans will be efficient in reducing the number of infected humans. Further, treatment can help in reducing the production rate of miracidia as seen in Figure 5(b). However, it does not seem to have an impact in reducing the production rate of cercariae, and thus, by not combining it with preventive strategies like WASH, individuals will continue to get infected or reinfected (upon treatment) (see Figure 5(c)). This is an indication that the effectiveness of treatment is a short-term measure. As a result, the susceptible and treated individuals should be educated on WASH procedures to effectively curtail the spread of the disease.

Sensitivity Analysis
Sensitivity indices permit us to estimate the perturbations in a variable with respect to changes in a parameter. The 16 Journal of Applied Mathematics normalized forward sensitivity index of a variable in relation to a parameter is the proportion of the perturbations in the variable with respect to changes in the given parameter. Suppose that a variable is partially differentiable with respect to the parameter, then the sensitivity index can be defined using partial derivatives.
Definition 5. The normalized forward sensitivity index of R 0 , which depends differentiably on a parameter ε, is defined by Table 2 contains all model parameters that are partially differentiable with respect to R 0 , their values, and sensitivity indices.

5.1.
Interpretation of the Sensitivity Indices on R 0 . The sensitivity indices for the model parameters to the basic reproduction number, R 0 , are presented in Table 2. The most sensitive parameters include the recruitment rates of humans and snails and the production rate of free living miracidia and cercariae as produced by the infected humans and infected snails, respectively. In general, the results indicate that the parameters α h , α s , β h , β s , σ h , σ s , λ 1 , and λ 2 have positive indices and, therefore, are each directly proportional to the value of R 0 . Thus, increasing the value of any of these parameters will lead to the disease remaining endemic within the studied population and vice versa. However, the parameters γ h , δ h , δ s , μ h , μ s , μ m , and μ c have each a negative sensitivity index, which means that each of them varies inversely as the value of R 0 . Therefore, increasing the value of any of these parameters with negative index, while holding all other parameters constant, will reduce the value of R 0 and, hence, contribute to the elimination of the disease and vice versa. For example, increasing the recruitment rate of either humans or snails by 10% will lead to a 5% increase on R 0 , and increasing the rate of treatment by 10% will result to a reduction of R 0 by 4:948%. Figure 6 shows the contributions of the parameters α h , α s , β h , β s , σ h , λ 1 , and λ 2 on the infected human subpopulation.

Effects of Parameters on Infected Humans.
Figures 6(a) and 6(b) indicate that reducing the recruitment rates of either humans or (most importantly) snails will significantly reduce the number of humans getting infected. The effects of reducing the forces of infection as seen in Figures 6(c) and 6(d) on the infected humans are an indication that vaccination or provision of quality water and improved sanitation system are highly essential in containing the spread of the disease. Figure 6(e) shows that reducing the rate exposed individuals become infected will significantly reduce the number of individuals that become infected eventually. This can be done through provision of early diagnostic materials that will be able to detect the disease in an exposed individual before the production of eggs and also invention of drugs that can treat the disease at the very onset of infection. Figures 6(f) and 6(g) show that reducing the production rate of either miracidia or cercariae will drastically reduce the number of infected individuals. Miracidia production can be reduced by treatment of infected individuals (see Figure 5(b)) or practice of WASH (see Figure 4(c)) or possibly through the introduction of vaccines for humans. Although miracidia can be reduced by treatment of infected individuals, the high possibility of them getting reinfected will make it very difficult to control it unless the individuals are vaccinated or sensitized on the key preventive mechanisms. Cercariae production, on the other hand, can be reduced through chlorination of fresh waters or snail control using chemical molluscicides.

Conclusion
A mathematical model that explains the transmission dynamics of the disease involving humans, snails, and the dynamics of the free living miracidia and the free living cercariae has been presented. Attention is paid to controls within the human subpopulation by incorporating the treated human compartment and preventive mechanism through water sanitation and hygiene (WASH). The equilibrium points were obtained and analyzed for both local and global stabilities. Considering the basic reproduction number, R 0 , as the threshold and given an R 0 < 1, it is established that the disease-free equilibrium point is both locally and globally asymptotically stable. The existence of a unique endemic equilibrium point when R 0 > 1 was established and using the Center Manifold Theory, it is confirmed that the endemic equilibrium is locally asymptotically stable when R 0 > 1 but close to 1. To validate these analytic results, a numerical simulation was carried out. The results showed that the model system follows a forward bifurcation at R 0 = 1, and hence, the disease-free equilibrium is globally asymptotically stable with R 0 < 1, and the endemic equilibrium is locally asymptotically stable when R 0 > 1.
The effectiveness of control by WASH and treatment of infected individuals was analyzed numerically. The results indicated that a combination of both controls is required for the proper eradication of the disease. However, WASH is a more effective control mechanism especially if it is exercised by the entire population.
Using the normalized sensitivity index of R 0 , we examined the contribution of the model parameters on the basic reproduction number, R 0 . The results demonstrate that the production rates of the free living miracidia and the free living cercariae are very effective contributors to R 0 and hence very essential in the spread and control of the disease. Thus, control mechanisms that can significantly mitigate the production of either pathogen will most definitely reduce the disease endemicity. It was shown that this could be achieved through the application of the combined control measures, i.e., treatment of infected individuals and WASH. Hence, infected individuals should be treated, and at the same time, communities should be educated regarding best practices on schistosomiasis prevention. Adequate provision of clean water and proper sanitation facilities will curtail the spread of the disease. Also, let v be the number of infectives leaving the system.
The Jacobian of f at E 0 is denoted by F with Similarly, the Jacobian of v at E 0 is denoted by V with ðA:6Þ Therefore, the next generation matrix is obtained as 18 Journal of Applied Mathematics where M = γ h + δ h + μ h , N = ρ + ð1 − ρÞδ h + μ h , and G = F V −1 : Accordingly, the spectral radius of the next generation matrix [20,22] which is given by ρðGÞ and defined as the basic reproductive number, R 0 , is obtained as

B. Theorem on Center Manifold Theory
Here, we reproduce Theorem 4.1 of Chavez and Song [23].
Theorem 6 (see Chavez and Song [23]). Consider the following general system of ordinary differential equations with a parameter ϕ.
where 0 is an equilibrium of the system, that is, f ð0, ϕÞ = 0 for all ϕ. Assume the following: A1: A = D x f ð0, 0Þ = ðð∂f i /∂x j Þð0, 0ÞÞ is the linearization of the system 24 around the equilibrium 0 with ϕ evaluated at 0. Zero is a simple eigenvalue of A, and other eigenvalues have negative real parts; A2: the matrix A has a right eigenvector w and a left eigenvector v corresponding to the zero eigenvalue.
Let The local dynamics of systems 24 around 0 are completely determined by the signs of a and b: (i) a > 0, b > 0. When ϕ < 0 with jϕj ≪ 1, 0 is locally asymptotically stable and there exists a positive unstable equilibrium; when 0 < ϕ ≪ 1, 0 is unstable and there exists a negative and locally asymptotically stable equilibrium; (ii) a < 0, b < 0. When ϕ < 0 with jϕj ≪ 1, 0 is unstable; when 0 < ϕ ≪ 1, 0 is locally asymptotically stable, and there exists a positive unstable equilibrium; (iii) a > 0, b < 0. When ϕ < 0 with jϕj ≪ 1, 0 is unstable, and there exists a locally asymptotically stable negative equilibrium; when 0 < ϕ ≪ 1, 0 is stable, and a positive unstable equilibrium appears; (iv) a < 0, b > 0. When ϕ changes from negative to positive, 0 changes its stability from stable to unstable. Correspondingly, a negative unstable equilibrium becomes positive and locally asymptotically stable.

Data Availability
The data used in this study are included within the article.

Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.