Mathematical Modeling and Analysis of Transmission Dynamics and Control of Schistosomiasis
This paper presents a mathematical model that describes the transmission dynamics of schistosomiasis for humans, snails, and the free living miracidia and cercariae. The model incorporates the treated compartment and a preventive factor due to water sanitation and hygiene (WASH) for the human subpopulation. A qualitative analysis was performed to examine the invariant regions, positivity of solutions, and disease equilibrium points together with their stabilities. The basic reproduction number, , is computed and used as a threshold value to determine the existence and stability of the equilibrium points. It is established that, under a specific condition, the disease-free equilibrium exists and there is a unique endemic equilibrium when . It is shown that the disease-free equilibrium point is both locally and globally asymptotically stable provided , and the unique endemic equilibrium point is locally asymptotically stable whenever using the concept of the Center Manifold Theory. A numerical simulation carried out showed that at , the model exhibits a forward bifurcation which, thus, validates the analytic results. Numerical analyses of the control strategies were performed and discussed. Further, a sensitivity analysis of was carried out to determine the contribution of the main parameters towards the die out of the disease. Finally, the effects that these parameters have on the infected humans were numerically examined, and the results indicated that combined application of treatment and WASH will be effective in eradicating schistosomiasis.
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, . Schistosoma mansoni and Schistosoma japonicum mainly cause diseases in the liver and bowels but Schistosoma haematobium mostly affects the urinary and genital areas .
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 . 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 .
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 .
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 . 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–11].
A review of a number of existing mathematical models of schistosomiasis, their merits, and demerits was carried out by .  suggested ways by which controls could be incorporated into mathematical models in order to eradicate or reduce the prevalence of schistosomiasis in an endemic area.  presented some of the measures adopted for the control schistosomiasis in the People’s Republic of China.  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 . 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.  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.
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 . 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.
2. 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 , exposed humans , infected humans , and treated humans ; the snail subpopulation into three compartments: susceptible snails , exposed snails , and infected snails ; and compartments for the free living miracidia and free living cercariae 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(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
where and .
The model is represented by the following system of ODEs:
With initial conditions
3. Model Analysis
3.1. 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, , is given by
In the absence of mortality due to infection, we have
Similarly, the total snail subpopulation is given by
In the absence of mortality due to infection, we get
Thus, , , , and are bounded, and the regionis positively invariant, and accordingly, the system equation (1) is epidemiologically meaningful and analytically well posed in the region .
Theorem 1. Suppose thatis the initial condition for the model system equation (1), then the set of solutionsis nonnegative .
Proof. Given that the set of initial conditions is nonnegative, consider the model system equation (1). Then,Separating the variables and integrating, we haveSimilarly,
3.2. Equilibrium Points
At the equilibrium point, all states are considered to have a constant solution. So, it is obtained by setting It is given by
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 which is denoted as . Accordingly,and a combination of the above equation (15) and the second equation of the model system 1, the following relation was obtained for : where
3.2.1. Disease-Free Equilibrium Point
At the disease-free equilibrium, there are no infections or recovery, and thus, no miracidia or cercariae is produced. Accordingly, at this point, in the system equation (16) must be zero. Hence, the point is obtained aswhere .
3.2.2. Endemic Equilibrium Point
At the endemic equilibrium point, all disease states in equation (14) are considered to be positive, and consequently, must be greater than zero for all the other states to be positive. Therefore, the unique positive endemic equilibrium point exists only when in the system equation (16) is positive, which is also in this case the same as . And thus, it is the pointwherewhenever .
3.3. Stability Analysis
3.3.1. Basic Reproduction Number,
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 . It is denoted by . We considered the next generation matrix approach according to Chavez et al.  (for details, see Appendix A). The next generation matrix is denoted by and is obtained as
Accordingly, the basic reproduction number, , which is the dominant eigenvalue of , is obtained as
For simplicity, since denotes the possible number of uninfected a particular infected individual (in this case either human or snail) can infect, we can rewrite as where represents the number of snails an infected human is capable of infecting through the production of miracidia and denotes the number of humans an infected snail is capable of infecting through the production of cercariae.
3.3.2. Stability of Disease-Free Equilibrium Point
Theorem 2. The disease-free equilibrium point, , is locally asymptotically stable provided and unstable otherwise.
Proof. To verify the validity of the above proposition, see Theorem 2 in Van den Driessche and Watmough .
The method according to Chavez et al.  is used to determine the global stability of the disease-free equilibrium point within . Hence, we denote the system equation (1) bywhere comprises of the uninfected subpopulations and comprises of the infected and infectious subpopulations. is guaranteed to be globally asymptotically stable (GAS) if and the following two conditions (C1) and (C2) hold:
(C1): For is GAS
where and is the invariant region of the model.
Theorem 3. The disease-free equilibrium point, , is GAS for the system equation (1) provided that and conditions (C1) and (C2) both hold.
Proof. From the system of equation (1), we haveFirst, suppose that the time, , we need to show that .
But we know that for ,Also,Thus, all points with respect to this condition converge atNow, from the system equation (1), we know thatwhereand , whereClearly, , sinceMoreover, is an M-matrix, since the nondiagonal entries are nonnegative . Therefore, the DFE, denoted , is globally asymptotically stable provided that .
3.3.3. Stability of the Endemic Equilibrium Point
We consider the concept of Center Manifold Theory as described in Theorem 4.1 by Chavez and Song . The theorem helps to determine the course of the bifurcation at a critical value of the parameter . The theorem is reproduced in Appendix B.
Theorem 4. The endemic equilibrium point, based on the Theorem 4.1 by  is locally asymptotically stable for (but near 1).
Proof. Now, consider the model system equation (1). Using Theorem 4.1 by Chavez and Song , let , , , , , , , , and . Further, denote and . Then, the model system equation (1) can be expressed in the form , whereThe role of the parameter in Theorem 4.1 by Chavez and Song  is played by with the critical value obtained from [7, 24]. Hence, , regarded as the bifurcation point, can be obtained by evaluating it at . SoLet the Jacobian of the transformed system equation (33) at when is denoted by . Then, it can be shown that 0 is a simple eigenvalue of the matrix , whereand if the vector is a right eigenvector corresponding to the eigenvalue, 0 of , then the components of can be obtained aswhere is defined in equation (34).
Similarly, let be a left eigenvector corresponding to the zero eigenvalue of . Then, the components of areNow, to get the values of and as in Theorem 4.1 by Chavez and Song , we have to get the second partial derivatives in the system equation (33). Clearly, we do not need the partial derivatives of , and , since . Thus, the second (nonzero) partial derivatives from the linearized system equation (33) for areand for isTherefore,So .So
Clearly, and , and with the help of condition (iv) in Theorem 4.1 by Chavez and Song , it is confirmed that the system follows a forward bifurcation at (see Figure 2), and hence, the endemic equilibrium is locally asymptotically stable for , but sufficiently closed to 1.
Generally, also based on the principles of the next generation matrix according to Van den Driessche and Watmough , the disease-free equilibrium point is stable when and unstable when . To validate these analytic results, a numerical simulation was carried out. Figure 2 shows a graph of against the steady state of , denoted by . The figure shows that at , the model system equation (1) exhibits a forward bifurcation. Forward bifurcations have been observed in many other epidemic models [7, 23, 25]. It shows an exchange of stability between the disease-free equilibrium point and the disease endemic equilibrium point. To this effect, it can be concluded that the model system equation (1) has a disease-free equilibrium point that is globally asymptotically stable when , and when , it posses an 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 , and when , a locally asymptotically stable unique endemic equilibrium point exists.
(a) Evolution of the human subpopulation with time
(b) Evolution of the snail subpopulation with time
(c) Evolution of the free living miracidia with time
(d) Evolution of the free living cercariae with time
(a) Effects of varying WASH on susceptible humans
(b) Effects of varying WASH on the infected humans
(c) Effects of varying WASH on miracidia
(d) Effects of varying WASH on cercariae
(a) Effects of varying treatments on the infected humans
(b) Effects of varying treatments on miracidia
(c) Effects of treatments on cercariae
(a) Effects of varying on infected humans
(b) Effects of varying on infected humans
(c) Effects of varying on infected humans
(d) Effects of varying on infected humans
(e) Effects of varying on infected humans
(f) Effects of varying on infected humans
(g) Effects of varying on infected humans
4. Numerical Results
Table 1 gives the parameter values used in the numerical simulation of the model equation (1). Many of the parameters are obtained from existing literature and are rightly cited while some are taken from real life data. Other parameters are estimated based on convenience using what is generally known about the dynamics of the studied population and the disease. In particular, the model is designed to study the dynamics of schistosomiasis within the Gambian population. In 2020, The Gambia has a population of , with a birth rate of per population and life expectancy of up to years. Thus, and . 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 . Hence, is estimated to lie between and . The incubation period for snails ranges from 4 weeks to 5 weeks . Thus, lies between and . The parasite can live within a definite human host for a period ranging from years to years . So, is varied between values corresponding to and years. Consequently, lie between and .
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).
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.
4.1. 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.
4.2. 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.
5. Sensitivity Analysis
Sensitivity indices permit us to estimate the perturbations in a variable with respect to changes in a parameter. The 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 , which depends differentiably on a parameter , is defined byTable 2 contains all model parameters that are partially differentiable with respect to , their values, and sensitivity indices.
5.1. Interpretation of the Sensitivity Indices on
The sensitivity indices for the model parameters to the basic reproduction number, , 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 , , , , , , , and have positive indices and, therefore, are each directly proportional to the value of . 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 , , , , , , and have each a negative sensitivity index, which means that each of them varies inversely as the value of . Therefore, increasing the value of any of these parameters with negative index, while holding all other parameters constant, will reduce the value of and, hence, contribute to the elimination of the disease and vice versa. For example, increasing the recruitment rate of either humans or snails by will lead to a increase on , and increasing the rate of treatment by will result to a reduction of by .
5.2. Effects of Parameters on Infected Humans
Figure 6 shows the contributions of the parameters , , , , , , and on the infected human subpopulation.
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.
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, , as the threshold and given an , 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 was established and using the Center Manifold Theory, it is confirmed that the endemic equilibrium is locally asymptotically stable when 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 , and hence, the disease-free equilibrium is globally asymptotically stable with , and the endemic equilibrium is locally asymptotically stable when .
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 , we examined the contribution of the model parameters on the basic reproduction number, . The results demonstrate that the production rates of the free living miracidia and the free living cercariae are very effective contributors to 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.
A. Computation of
The next generation matrix approach according to Van Den Driessche and Watmough  is used to compute . Let be the infected or infectious class of the system of differential equation (1), with . Consider the following:
Let be the number of new infections coming into the population. Then,
Also, let be the number of infectives leaving the system. Then,
The Jacobian of at is denoted by with
Similarly, the Jacobian of at is denoted by withwith
Therefore, the next generation matrix is obtained aswhere , , and Accordingly, the spectral radius of the next generation matrix [20, 22] which is given by and defined as the basic reproductive number, , is obtained as
B. Theorem on Center Manifold Theory
Here, we reproduce Theorem 4.1 of Chavez and Song .
Theorem 6 (see Chavez and Song ). Consider the following general system of ordinary differential equations with a parameter .where 0 is an equilibrium of the system, that is, for all . Assume the following:
A1: is the linearization of the system 24 around the equilibrium 0 with evaluated at 0. Zero is a simple eigenvalue of , and other eigenvalues have negative real parts;
A2: the matrix has a right eigenvector w and a left eigenvector v corresponding to the zero eigenvalue.
Let be the component of andThe local dynamics of systems 24 around 0 are completely determined by the signs of a and b:(i). When with , 0 is locally asymptotically stable and there exists a positive unstable equilibrium; when , 0 is unstable and there exists a negative and locally asymptotically stable equilibrium;(ii). When with , 0 is unstable; when , 0 is locally asymptotically stable, and there exists a positive unstable equilibrium;(iii). When with , 0 is unstable, and there exists a locally asymptotically stable negative equilibrium; when , 0 is stable, and a positive unstable equilibrium appears;(iv). When