On the mathematical modeling of schistosomiasis transmission dynamics with heterogeneous intermediate host

Schistosomiasis is a neglected disease affecting almost every region of the world, with its endemicity mainly experience in sub-Saharan Africa. It remains difficult to eradicate due to heterogeneity associated with its transmission mode. A mathematical model of Schistosomiasis integrating heterogeneous host transmission pathways is thus formulated and analyzed to investigate the impact of the disease in the human population. Mathematical analyses are presented, including establishing the existence and uniqueness of solutions, computation of the model equilibria, and the basic reproduction number (R₀). Stability analyses of the model equilibrium states show that disease-free and endemic equilibrium points are locally and globally asymptotically stable whenever R₀ < 1 and R₀>1, respectively. Additionally, bifurcation analysis is carried out to establish the existence of a forward bifurcation around R₀ = 1. Using Latin-hypercube sampling, global sensitivity analysis was performed to examine and investigate the most significant model parameters in R₀ which drives the infection. The sensitivity analysis result indicates that the snail's natural death rate, cercariae, and miracidia decay rates are the most influential parameters. Furthermore, numerical simulations of the model were done to show time series plots, phase portraits, and 3-D representations of the model and also to visualize the impact of the most sensitive parameters on the disease dynamics. Our numerical findings suggest that reducing the snail population will directly reduce Schistosomiasis transmission within the human population and thus lead to its eradication.

1. Introduction

Schistosomiasis which is also known as Bilharzia, Bilharziasis, and snail fever, is one of the neglected tropical diseases (NTD) common in some tropical and subtropical countries, in areas where there is no access to safe drinking water, adequate hygiene, and appropriate sanitation (sewage disposal). It is a water-based disease that is transmitted through contact with water infested with cercariae and it affects about 240 million people worldwide and more than 700 million people live in an area where the disease is endemic [1]. Schistosomiasis infection is also known as a parasite disease that is caused by schistosomiasis organisms (blood fluke). The two main forms of schistosomiasis infection are urogenital and intestinal which are caused by five schistosomiasis species responsible for human or animal infections namely, Schistosoma mansoni, Schistosoma intercalatum, Schistosoma haematobium, Schistosoma mekonai, and Schistosoma japonicum [1]. There are two stages of schistosomiasis infection namely; acute and chronic stages. Acute schistosomiasis infection referred to as Katayama Syndrome is a clinical manifestation of infection caused by trematode worms [2].

When non-immune persons are exposed to schistosomiasis infection, there is an incubation period symptoms of infection include abdominal pain, blood in stool or urine, wasting, anemia, dysuria [2]. Chronic schistosomiasis infection is caused by schistosome, which has severe long-term implications [3]. The symptoms of the chronic stage are growth retardation, renal failure, ascites, grand mal epilepsy, decreased fitness, anemia, portal hypertension, transverse myelitis, bladder carcinoma, and infertility [2]. Normally, snails harbor cercariae and shed the infective stage (cercariae) of the parasite (schistosome) in the freshwater [4]. Hence, there is an increase in the prevalence of Schistosomiasis as a result of freshwater habitats such as dams and irrigation, which enhances snail habitation and thus, increases transmission. Despite the progressive achievements in schistosomiasis control over the years, Schistosomiasis is still one of the most devastating parasitic diseases ravaging the entire world. Guiro et al. [5] proposed a mathematical model for schistosomiasis control, considering two population hosts, humans, and snails with delays. Their results suggest that an effective education campaign and reasonable coverage level of drug treatment will help in the control of Schistosomiasis. A mathematical model for the human-cattle-snail transmission of Schistosomiasis in Hubei Province, China was studied by Chen et al. [6]. Their results show that to control or eliminate Schistosomiasis in Hubei Province, China. Environmental factors need to be considered to lower (narrow) cattle-snail transmission. Since chemotherapy alone cannot end the spread of the causative parasite (schistosome), additional interventions and control strategies must be combined to lower re-infection, reduce the prevalence, and move toward the elimination of Schistosomiasis [7]. Many researchers have considered the use of molluscicides and different kinds of environmental modifications to control snails to reduce or eliminate schistosomiasis prevalence [7–9]. Interested readers can also see the following references [10–20] for other published works that have used mathematical models to gain insight into the transmission dynamics of Schistosomiasis and other infectious diseases in the human population and various suggestion for their control/elimination.

Given the aforementioned above, our study presents a mathematical modeling work that examines qualitative mathematical analysis for the prevention and control of the spread of Schistosomiasis through the control of snails. This is in addition to chemotherapy, molluscicides, and other environmental factors already studied in the previous research works and also the inclusion of the acute and chronic infected populations (compartments). The present study also investigates if mammals such as cattle contribute to the spread of Schistosomiasis infection within the population. The model described in this work considers three population hosts that is: humans, cattle, and snails.

The paper is organized as follows: Section 2 gives the Model formulation, Section 3 the Mathematical model analysis is presented, Section 4 presents the Sensitivity analysis of parameters in the model, Section 5 presents the Results of the numerical simulations and finally, Section 6 gives a brief Discussion and Conclusion.

2. Model formulation

The life cycle of the schistosome parasite is formulated mathematically. The populations considered are the human population, other mammal population, and snail population with the Cercariae C(t) and Miracida M(t). Cercariae C(t) is a larva worm that the infected snails shed into the aquatic environment, while the Miracida M(t) are eggs that the infected mammals (human and other mammals) shed in the stream when they come to the river for human activities such as fishing, swimming, fetching water, and/or drinking of water by other mammals. In this study, other mammals include cattle, mice, and dogs that the schistosome parasite can infect. The total human population, N_h(t), at time, t, is divided into the sub-populations; susceptible, S_h(t), acute infected, I_A(t), and chronic infected, I_ch(t) humans. The other mammal population, N_a(t), is subdivided into susceptible mammals, S_a and infected mammals, I_a, while the snail population, N_s(t) is made up of susceptible snails, S_s(t) and infected snails I_s(t). Human recruitment and natural death rates are assumed to be Λ_h and μ_h, respectively. The susceptible humans, S_h(t), become infected through contact with fresh water contaminated with Cercariae from infected snails. The force of infection is given by

$\begin{array}{ll}\frac{{\beta }_h{\alpha }_{h}C}{C_0+{\epsilon }_{1}C},& (1)\end{array}$

where β_h is the human transmission rate, α_h is the adult parasite production rate for humans, C₀, is the half-saturation constant of Cercariae within the environment and ε₁ is the limitation of the growth velocity of the pathogens for Cercariae. Susceptible humans become acutely infected upon infection and enter into class, I_A(t). In the acutely infected sub-population, a proportion of them may recover with partial immunity depending on the mass drug administration program within the community and return to the susceptible class at a rate, ψ. On the other hand, some human may develop a chronic infection at a rate, k and move to class, I_ch. Chronically infected humans may succumb to infection at a rate, δ_h. Shedding of infection within the environment by acute infected and chronic infected humans is assumed to occur at rate N_Eτ₁γ_h and N_Eτ₁γ_hσ, respectively. Here, N_E is the number of eggs shed by mammals (both human and other mammals), τ₁ is the probability of eggs developing into Miracidium, γ_h is the infected human shedding rate and σ is a parameter that influences the shedding rate of the chronically infected humans. For the other mammals, the recruitment rate is assumed to be Λ_a while the natural death rate is μ_a. The susceptible mammals are infected through contact with fresh water contaminated with Cercariae from infected snails and the force of infection is given by

$\begin{array}{ll}\frac{{\beta }_a{\alpha }_{a}C}{C_0+{\epsilon }_{1}C},& (2)\end{array}$

where β_a is the other mammal transmission rate, α_a, is the adult parasite production probability for other mammals. The other mammals shed the eggs that hatch into the water at a rate, N_Eτ₁γ_a where γ_a is the other mammals shedding rate. Death due to infection among other infected mammals is assumed to be at a rate, δ_a. Furthermore, susceptible snails are recruited at a rate, Λ_s and death due to natural causes of snails is at a rate, μ_s. Susceptible snails become infected by contact with miracidia from the viral shedding of infected humans and other mammals and the force of infection is given by

$\begin{array}{ll}\frac{{\beta }_s{\alpha }_{s}M}{M_0+{\epsilon }_{2}M},& (3)\end{array}$

where β_s, is the snail transmission rate, α_s is the Miracidial penetration probability for snails, ε₂, is the limitation of the growth velocity for the Miracidium and M₀, is the half-saturation constant for Miracida within the environment. Death due to infection of snails is assumed to be at a rate, δ_s. Infected snails shed larva worm in the environment at a rate, τ₂γ_s where τ₂ is the density of Cercariae for snails and γ_s is the snail shedding rate. Finally, the decay rates for the Cercariae and Miracida are μ_c, μ_M, respectively. The study assumes that there is no immigration of infectious humans, other mammals, and snails within their populations. Also, susceptible humans, other mammals and snails are recruited at a constant rate. The dynamics of infection are presented in the model flow diagram in Figure 1. This leads to the following system of non-linear differential equations:

$\begin{array}{ll}\{\begin{array}{ll}\frac{d{S}_h}{dt}={\Lambda }_h-\frac{{\beta }_h{\alpha }_h{S}_{h}C}{C_0+{\epsilon }_{1}C}-{\mu }_h{S}_h+\psi I_A,\hfill & \hfill \\ \hfill & \hfill \\ \frac{d{I}_A}{dt}=\frac{{\beta }_h{\alpha }_h{S}_{h}C}{C_0+{\epsilon }_{1}C}-(k+{\mu }_h+\psi )I_A,\hfill & \hfill \\ \hfill & \hfill \\ \frac{d{I}_{ch}}{dt}=k{I}_A-({\mu }_h+{\delta }_h)I_{ch},\hfill & \hfill \\ \hfill & \hfill \\ \frac{d{S}_a}{dt}={\Lambda }_a-\frac{{\beta }_c{\alpha }_{a}C{S}_a}{C_0+{\epsilon }_{1}C}-{\mu }_a{S}_a,\hfill & \hfill \\ \hfill & \hfill \\ \frac{d{I}_a}{dt}=\frac{{\beta }_a{\alpha }_{a}C{S}_a}{C_0+{\epsilon }_{1}C}-({\mu }_a+{\delta }_a)I_a,\hfill & \hfill \\ \hfill & \hfill \\ \frac{dM}{dt}=N_E{\tau }_1({\gamma }_h(I_A+\sigma I_{ch})+{\gamma }_a{I}_a)-{\mu }_{m}M,\hfill & \hfill \\ \hfill & \hfill \\ \frac{d{S}_s}{dt}={\Lambda }_s-\frac{{\beta }_s{\alpha }_{s}M{S}_s}{M_0+{\epsilon }_{2}M}-{\mu }_s{S}_s,\hfill & \hfill \\ \hfill & \hfill \\ \frac{d{I}_s}{dt}=\frac{{\beta }_s{\alpha }_{s}M{S}_s}{M_0+{\epsilon }_{2}M}-({\mu }_s+{\delta }_s)I_s,\hfill & \hfill \\ \hfill & \hfill \\ \frac{dC}{dt}={\tau }_2{\gamma }_s{I}_s-{\mu }_{c}C\hfill & \hfill \end{array}& (4)\end{array}$

with initial conditions

$\begin{array}{l}{S}_h(0)>0,I_A(0)\ge 0,I_h(0)\ge 0,S_a(0)>0,I_a(0)\ge 0,\\ S_s(0)>0,I_s(0)\ge 0,M(0)\ge 0,C(0)\ge 0.\end{array}$

FIGURE 1

Figure 1. Flow diagram for the dynamics of Schistosomiasis infection.

All the parameters are assumed to be non-negative over the modeling time frame.

3. Mathematical methods

3.1. Invariant region

A biologically feasible population model such as model system (4) should always have positive solutions. To show that all solutions of model system (4) are non-negative, we let N_h, N_a and N_s be the total human population, other mammals population and snail populations, respectively, where, N_h(t) = S_h(t) + I_A(t) + I_ch(t), N_a(t) = S_a(t) + I_a(t), and N_s(t) = S_s(t) + I_s(t), respectively with the initial conditions N_h(0), N_a(0), N_s(0), M(0), and C(0).

Adding the human, other mammals and snail populations give the respective total populations as$N_h{}^{\prime }(t)={\Lambda }_h-{\delta }_h{I}_{ch}-N_h{\mu }_h,N_a{}^{\prime }(t)={\Lambda }_a-{\delta }_a{I}_a-N_a{\mu }_a$ and $N_s{}^{\prime }(t)={\Lambda }_s-{\delta }_s{I}_s-N_s{\mu }_s$. In the absence of infection which implies that the disease-related death rates for chronically infected human, δ_h, other infected mammals, δ_a and infected snail, δ_s are negligible, we assume that the total populations, N_h, N_a and N_s are asymptotically stable and yield

$\begin{array}{l}{N}_h{}^{\prime }(t)\le {\Lambda }_h-N_h{\mu }_h,\\ N_a{}^{\prime }(t)\le {\Lambda }_a-N_a{\mu }_a,\\ N_s{}^{\prime }(t)\le {\Lambda }_s-N_s{\mu }_s.& (5)\end{array}$

Solving (5) using the inequality theorem by Birkhoff and Rota [21] and integrating with the initial conditions N_h(0), N_a(0), N_s(0), C(0), and M(0) gives

$\begin{array}{l}{N}_h\le \frac{{\Lambda }_h}{{\mu }_h}-\left(\frac{{\Lambda }_h-{\mu }_h{N}_{h0}}{{\mu }_h}\right)\mathrm{ exp}(-{\mu }_{h}t),\\ N_a\le \frac{{\Lambda }_a}{{\mu }_a}-\left(\frac{{\Lambda }_a-{\mu }_a{N}_{a0}}{{\mu }_a}\right)\mathrm{ exp}(-{\mu }_{a}t),\\ N_s\le \frac{{\Lambda }_s}{{\mu }_s}-\left(\frac{{\Lambda }_s-{\mu }_s{N}_{s0}}{{\mu }_s}\right)\mathrm{ exp}(-{\mu }_{s}t).\end{array}$

As t → ∞, the total populations for humans, other mammals and snails approach $N_h\le \frac{{\Lambda }_h}{{\mu }_h}$, $N_a\le \frac{{\Lambda }_a}{{\mu }_a}$ and $N_s\le \frac{{\Lambda }_s}{{\mu }_s}$, respectively. Since I_A, I_ch ≤ N_h(t), I_a ≤ N_a(t), and I_s ≤ N_s(t) it implies that the Cercariae (C) and Miracidia (M) populations can be written as

$\begin{array}{l}{C}^{\prime }(t)\le \frac{{\gamma }_s{\tau }_2{\Lambda }_s}{{\mu }_s}-{\mu }_{c}C,\\ M^{\prime }(t)\le N_E{\tau }_1\left[\frac{{\gamma }_h(1+\sigma ){\Lambda }_h}{{\mu }_h}+\frac{{\gamma }_a{\Lambda }_a}{{\mu }_a}\right]-{\mu }_{m}M.& (6)\end{array}$

Applying the inequality theorem by Birkhoff and Rota [21] with the initial conditions C(0) = C₀₀ and M(0) = M₀₀ yields

$C\left(t\right)\le \frac{{\Lambda }_s{\gamma }_s{\tau }_2}{{\mu }_s{\mu }_c}-\left[\frac{{\gamma }_s{\tau }_2{\Lambda }_s-{\mu }_s{\mu }_c{C}_{00}}{{\mu }_s{\mu }_c}\right]\exp (-{\mu }_{c}t),$

and

$\begin{array}{l}M(t)\le N_E{\tau }_1[\frac{{\Lambda }_h{\gamma }_h(1+\sigma )}{{\mu }_h{\mu }_m}+\frac{{\Lambda }_a{\gamma }_a}{{\mu }_a{\mu }_m}]\\ -\left[\frac{N_E{\tau }_1({\Lambda }_h{\gamma }_h{\mu }_a(1+\sigma )+{\Lambda }_a{\gamma }_a{\mu }_h)-{\mu }_h{\mu }_a{\mu }_m{M}_{00}}{{\mu }_h{\mu }_a{\mu }_m}\right]\\ \exp (-{\mu }_{m}t).\end{array}$

As t → ∞, the Cercariae and Miracidia populations tends to $C(t)\le \frac{{\Lambda }_s{\gamma }_s{\tau }_2}{{\mu }_s{\mu }_c}andM(t)\le N_E{\tau }_1\left[\frac{{\Lambda }_h{\gamma }_h(1+\sigma )}{{\mu }_h{\mu }_m}+\frac{{\Lambda }_a{\gamma }_a}{{\mu }_a{\mu }_m}\right]$ respectively.

This shows that all the feasible solutions N_h(t), N_a(t), N_s(t), C(t) and M(t) are bounded and the region,

$\begin{array}{l}{Z}=\{(S_h,I_A,I_{ch},S_a,I_a,M,S_s,I_s,C)\in {ℝ}_{+}^9\mid N_h\le \frac{{\Lambda }_h}{{\mu }_h},\\ N_a\le \frac{{\Lambda }_a}{{\mu }_a}, N_s\le \frac{{\Lambda }_s}{{\mu }_s}, M\le M_{max},C\le C_{max}\}\end{array}$

is positively invariant. Thus, we conclude that for all t>0, model system (4) is biologically feasible and well-posed in ${Z}$.

3.2. Positivity of the solution

We show that the solutions of the model system (4) are positive for all time by proving the following theorem:

Theorem 3.1. With positive initial conditions, the solution set of model system (4), [S_h(t),I_A(t),I_ch(t), S_a(t), I_a(t), M(t), S_s(t), I_s(t), C(t)], is positive for all time, t > 0.

Proof. From the first equation of the model Equation (4), we have

$\begin{array}{l}{S}_h{}^{\prime }={\Lambda }_h-\frac{{\beta }_h{\alpha }_{h}C{S}_h}{C_0+{\epsilon }_{c}C}-{\mu }_h{S}_h+\psi I_A\ge -\left(\frac{{\beta }_h{\alpha }_{h}C}{C_0+{\epsilon }_{c}C}+{\mu }_h\right)S_h.(7)\end{array}$

Integrating with initial condition S_h(0), we obtain

$\begin{array}{l}{S}_h(t)\ge S_h(0)exp\{{\displaystyle {\int }_0^t}-(y_h+{\mu }_h)dt\}>0(8)\end{array}$

with $y_h=\frac{{\beta }_h{\alpha }_{h}C}{C_0+{\epsilon }_{1}C}.$

In a similar way with $y_a=\frac{{\beta }_a{\alpha }_{a}C}{C_0+{\epsilon }_{1}C}$ and $y_s=\frac{{\beta }_s{\alpha }_{s}M}{M_o+{\epsilon }_{2}M},$

$\begin{array}{l}{I}_A(t)\ge I_A(0)exp\{-(k+{\mu }_h+\psi )t\}>0,\\ I_{ch}\ge I_{ch}(0)exp\{-({\mu }_h+{\delta }_h)t\}>0,\\ S_a(t)\ge S_a(0)exp\{{\displaystyle {\int }_0^t}-(y_a+{\mu }_a)dt\}>0,\\ I_a(t)\ge I_a(0)exp\{-({\mu }_a+{\delta }_a)t\}>0,\\ S_s(t)\ge S_s(0)exp\{{\displaystyle {\int }_0^t}-(y_s+{\mu }_s)dt\}>0,\\ I_s(t)\ge I_s(0)exp\{-({\mu }_s+{\delta }_s)t\}>0,\\ C(t)\ge C(0)exp\{-{\mu }_{c}t\}>0,\end{array}$

$\begin{array}{l}M(t)\ge M(0)exp\{-{\mu }_{m}t\}>0.(9)\end{array}$

Therefore, we conclude that the solution set [S_h(t), I_A(t), I_ch(t), S_a(t), I_a(t), M(t), S_s(t), I_s(t), C(t)] of model Equation (4) is positive for all t > 0. This shows that the model is well-posed and biologically meaningful since the sub-population cannot be negative.

3.3. Existence of disease-free equilibrium

We examine the stability of the disease-free equilibrium state of the schistosomiasis model system (4) in terms of basic reproduction number. The disease-free equilibrium, E₀, is the equilibrium state in the absence of schistosomiasis disease in the environment. At equilibrium state,

$\begin{array}{l}\frac{d{S}_h}{dt}=\frac{d{I}_A}{dt}=\frac{d{I}_{ch}}{dt}=\frac{d{S}_a}{dt}=\frac{d{I}_a}{dt}=\frac{d{S}_s}{dt}=\frac{I_s}{dt}=\frac{dM}{dt}\\ =\frac{dC}{dt}=0.\end{array}$

This is solved simultaneously for the schistosomiasis model Equation (4) to yield the disease-free equilibrium,

$\begin{array}{l}{E}_0=(S_h^0,I_A^0,I_{ch}^0,S_a^0,I_a^0,S_s^0,I_s^0,M^0,C^0)\\ =\left(\frac{{\Lambda }_h}{{\mu }_h},0,0,\frac{{\Lambda }_a}{{\mu }_a},0,\frac{{\Lambda }_s}{{\mu }_s},0,0,0\right).(10)\end{array}$

3.3.1. The basic reproduction number

The basic reproduction number of the schistosomiasis model system (4) denoted by R₀, is a threshold parameter that is interpreted biologically as the average number of infected individuals created by one infected individual introduced in a susceptible population during the course of its infectious period [22, 23]. We compute R₀ using the next-generation matrix approach by Van den Driessche and Watmough [22]. Applying the approach, let ${F}$ and ${V}$ denote the Jacobian matrices for the rates of appearance of new infections and the transfer of infections in and out of the infected compartments, respectively at DFE, E₀, so that

${F}=\left(\begin{array}{cccccc}0& 0& 0& 0& 0& \frac{{\beta }_h{\alpha }_h{S}_h^0}{C_0}\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& \frac{{\beta }_a{\alpha }_a{S}_a^0}{C_0}\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& \frac{{\beta }_s{\alpha }_s{S}_s^0}{M_0}& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right) \text{and}$

$\begin{array}{ll}{V}=\left(\begin{array}{cccccc}f& 0& 0& 0& 0& 0\\ -k& g& 0& 0& 0& 0\\ 0& 0& h& 0& 0& 0\\ -A& -A\sigma & -B& {\mu }_m& 0& 0\\ 0& 0& 0& 0& q& 0\\ 0& 0& 0& 0& -D& {\mu }_c\end{array}\right),& (11)\end{array}$

where

$\begin{array}{l}f=k+{\mu }_h+\psi ,h={\mu }_a+{\delta }_a,g={\mu }_h+{\delta }_h,q={\mu }_s+{\delta }_s,\\ A=N_E{\gamma }_h{\tau }_1,B=N_E{\gamma }_a{\tau }_1,D={\gamma }_s{\tau }_2.(12)\end{array}$

The basic reproduction number, R₀, for the model system (4) is defined as the spectral radius of the operator, ${F}{{V}}^{-1}$, where ${{V}}^{-1}$ is the inverse matrix of ${V}$. This gives

$\begin{array}{l}{R}_0=\sqrt{\frac{{\beta }_s{\alpha }_s{\Lambda }_{s}D}{q{\mu }_c{\mu }_s{M}_0}(\frac{{\beta }_h{\alpha }_h{\Lambda }_{h}A(k\sigma +g)}{fg{\mu }_h{C}_0{\mu }_m}+\frac{{\beta }_a{\alpha }_a{\Lambda }_{a}B}{h{\mu }_a{C}_0{\mu }_m})}\\ =\sqrt{R_{sch}+R_{sA}+R_{sm}}(13)\end{array}$

for

$\begin{array}{l}{R}_{sch}=R_s\times R_{ch},R_{sA}=R_s\times R_A,R_{sm}=R_s\times R_m,\\ R_s=\frac{{\beta }_s{\alpha }_{s}D{\Lambda }_s}{{\mu }_{c}q{M}_0{\mu }_s},\\ R_m=\frac{{\beta }_a{\alpha }_{a}B{\Lambda }_a}{h{\mu }_m{C}_0{\mu }_a},R_{ch}=\frac{{\beta }_h{\alpha }_{h}Ak\sigma {\Lambda }_h}{fg{\mu }_m{C}_0{\mu }_h},\\ R_A=\frac{{\beta }_h{\alpha }_{h}A{\Lambda }_h}{f{\mu }_m{C}_0{\mu }_h}.& (14)\end{array}$

Here, the reproduction numbers, R_A, R_ch, R_m, and R_s, are the reproduction numbers for the acutely infected humans, chronic infected humans, other infected mammals and infected snails contribute, respectively. This shows that the basic reproduction number, R₀ of the schistosomiasis model system (4) comprises of four parts, that is, R_A, R_ch, R_m, R_s. Importantly, whenever R₀ < 1, it implies that the disease will fizzle out in the population and whenever R₀>1, it means that the disease will persist in the population.

3.4. Local stability of the disease-free equilibrium

We prove the local stability of the disease-free equilibrium, E₀, for the Schistosomiasis model system (4) using the Jacobian method.

Theorem 3.2. The disease-free equilibrium state, E₀, of the schistosomiasis model system (4) is locally asymptotically stable when R₀ < 1 and unstable for R₀ > 1.

Proof. We compute the Jacobian matrix of the model system (4) at the equilibrium state, E₀, as follows.

$\begin{array}{ll}J(E_0)=\left(\begin{array}{ccccccccc}-{\mu }_h& \psi & 0& 0& 0& 0& 0& 0& -\frac{{\beta }_h{\alpha }_h{S}_h^0}{C_0}\\ & & & & & & & & \\ 0& -f& 0& 0& 0& 0& 0& 0& \frac{{\beta }_h{\alpha }_h{S}_h^0}{C_0}\\ & & & & & & & & \\ 0& k& -g& 0& 0& 0& 0& 0& 0\\ & & & & & & & & \\ 0& 0& 0& -{\mu }_a& 0& 0& 0& 0& -\frac{{\beta }_a{\alpha }_a{S}_a^0}{C_0}\\ & & & & & & & & \\ 0& 0& 0& 0& -h& 0& 0& 0& \frac{{\beta }_a{\alpha }_a{S}_a^0}{C_0}\\ & & & & & & & & \\ 0& A& \sigma A& 0& B& -{\mu }_m& 0& 0& 0\\ & & & & & & & & \\ 0& 0& 0& 0& 0& -\frac{{\beta }_s{\alpha }_s{S}_s^0}{M_0}& -{\mu }_s& 0& 0\\ & & & & & & & & \\ 0& 0& 0& 0& 0& \frac{{\beta }_s{\alpha }_s{S}_s^0}{M_0}& 0& -q& 0\\ 0& 0& 0& 0& 0& 0& 0& D& -{\mu }_c\end{array}\right),& (15)\end{array}$

where f, g, h, q, A, B, D are defined in (12). The eigenvalues of (15) are λ₁ = −μ_h, λ₂ = −μ_a, λ₃ = −μ_s and the roots of the characteristics equation

$\begin{array}{l}{\lambda }^6+P{\lambda }^5+Q{\lambda }^4+R{\lambda }^3+U{\lambda }^2+V\lambda +W=0(16)\end{array}$

where

$\begin{array}{l}\{\begin{array}{l}P={\mu }_c+{\mu }_m+q+h+g+f,\hfill \\ Q={\mu }_c({\mu }_m+f+g+h+q)+{\mu }_m(f+g+h+q)\hfill \\ +f(g+h+q)+g(h+q)+hq,\hfill \\ R={\mu }_m{\mu }_c(f+g+h+q)+({\mu }_{c}f+{\mu }_{m}f)(g+h+q)\hfill \\ +({\mu }_{m}g+{\mu }_{c}g+fg)(h+q)+hq({\mu }_c+f+g),\hfill \\ U={\mu }_m{\mu }_{c}qf(1-R_{sA})+{\mu }_m{\mu }_{c}qh(1-R_{sm})\hfill \\ +{\mu }_m{\mu }_{c}f(g+h)\hfill \\ +({\mu }_m{\mu }_{c}g+{\mu }_{c}fg+{\mu }_{m}fg)(h+q)\hfill \\ +({\mu }_{c}hq+{\mu }_{m}hq)(f+g)+fghq,\hfill \\ V={\mu }_m{\mu }_{c}gq(1-R_{sch}-R_{sA})\hfill \\ +{\mu }_m{\mu }_{c}fhq(1-R_{sA}-R_{sm})\hfill \\ +{\mu }_m{\mu }_{c}qhg(1-R_{sm})+fgh({\mu }_m{\mu }_c\hfill \\ +{\mu }_{c}q+{\mu }_{m}q),\hfill \\ W=fghq{\mu }_m{\mu }_c{M}_0{C}_0{\mu }_s{\mu }_a{\mu }_h(1-R_0^2).\hfill \end{array}(17)\end{array}$

With the definitions of reproduction numbers in (14), we have by Routh-Hurwitz criteria in Kim et al. [24] and the conditions in Heffernan et al. [25] that the polynomial in (16) has negative real roots provided that P > 0, Q > 0, R > 0, U > 0, V > 0, W > 0, PQ > R, QU > W, UV > RW and this happens when R₀ < 1. With negative real eigenvalues, it implies that the disease-free equilibrium state, E₀, for Schistosomiasis model system (4) is locally asymptotically stable if R₀ < 1 otherwise it is unstable.

3.5. Existence of endemic equilibrium

For the endemic equilibrium, let

$\begin{array}{l}{\lambda }_c^{*}=\frac{C^{*}}{C_0+{\epsilon }_1{C}^{*}},{\lambda }_m^{*}=\frac{M^{*}}{M_0+{\epsilon }_2{M}^{*}},(18)\end{array}$

in model system (4). We have at equilibrium state that

$\begin{array}{l}{S}_h^{*}=\frac{{\Lambda }_h-({\mu }_h+k)I_A^{*}}{{\mu }_h},I_A^{*}=\frac{{\lambda }_c^{*}{\beta }_h{\alpha }_h{\Lambda }_h}{{\beta }_h{\alpha }_h({\mu }_h+k){\lambda }_c^{*}+f{\mu }_h},\\ I_{ch}^{*}=\frac{k}{g}{I}_A^{*},S_a^{*}=\frac{{\Lambda }_a-h{I}_a^{*}}{{\mu }_a},\\ I_a^{*}=\frac{{\lambda }_c^{*}{\beta }_a{\alpha }_a{\Lambda }_a}{{\beta }_a{\alpha }_a{\lambda }_c^{*}h+h{\mu }_a},S_s^{*}=\frac{{\Lambda }_s-q{I}_s^{*}}{{\mu }_s},\\ I_s^{*}=\frac{{\beta }_s{\alpha }_s{\Lambda }_s{\lambda }_m^{*}}{{\beta }_s{\alpha }_{s}q{\lambda }_m^{*}+{\mu }_{s}q},C^{*}=\frac{D}{{\mu }_c}{I}_s^{*},\end{array}$

$\begin{array}{l}{M}^{*}=\frac{A(I_A^{*}+\sigma I_{ch}^{*})+B{I}_a*}{{\mu }_m}.(19)\end{array}$

Substituting $S_h^{*},I_A^{*},I_{ch}^{*},S_a^{*},I_a^{*},S_s^{*},I_s^{*},M^{*}$, and C^* into ${\lambda }_c^{*}=\frac{C^{*}}{C_0+{\epsilon }_1{C}^{*}}$ and ${\lambda }_m^{*}=\frac{M^{*}}{M_0+{\epsilon }_2{M}^{*}}$ and solve simultaneously, we obtain

${\lambda }_m^{*}=\frac{{\mu }_c{\mu }_{s}q{C}_0{\lambda }_c^{*}}{{\beta }_s{\alpha }_s(D{\Lambda }_s-{\mu }_{c}q{C}_0{\lambda }_c^{*}-D{\Lambda }_s{\epsilon }_1{\lambda }_c^{*})}$

with ${\lambda }_c^{*}$ as solution of the polynomial

$\begin{array}{l}{\lambda }_c^{*}(X{\lambda }_c^{*2}+Y{\lambda }_c^{*}+Z)=0,(20)\end{array}$

where

$\begin{array}{l}X=Ah{\Lambda }_h{\beta }_a{\beta }_h{\alpha }_h{\alpha }_a(k\sigma +g)[{\mu }_c{\mu }_{s}q{C}_0{\epsilon }_1+{\beta }_s{\alpha }_s{\mu }_{c}q{C}_0\\ +{\beta }_s{\alpha }_{s}D{\Lambda }_s{\epsilon }_1]+{\beta }_a{\beta }_h{\alpha }_h{\alpha }_{a}g({\mu }_h+k)\\ \times [{\mu }_c{\mu }_{s}q{C}_0({\mu }_{m}h{M}_0+{\epsilon }_{2}B{\Lambda }_a)\\ +{\beta }_s{\alpha }_{s}B{\Lambda }_a(D{\Lambda }_s{\epsilon }_1+{\mu }_{c}q{C}_0)],\\ Y=(g+k\sigma ){\beta }_h{\alpha }_{h}Ah{\Lambda }_h[{\Lambda }_s{\beta }_s{\alpha }_{s}D({\epsilon }_1{\mu }_a-{\beta }_a{\alpha }_a)\\ +{\mu }_{c}q{C}_0{\mu }_a({\beta }_s{\alpha }_s+{\mu }_s{\epsilon }_2)]+({\mu }_h+k){\beta }_h{\alpha }_{h}g\\ \times ({\mu }_c{\mu }_{s}q{C}_0{\mu }_{m}h{M}_0{\mu }_a-{\beta }_s{\alpha }_{s}BD{\Lambda }_a{\Lambda }_s{\beta }_a{\alpha }_a)\\ +{\mu }_c{\mu }_{s}qfg{C}_0{\beta }_a{\alpha }_a{\mu }_h({\mu }_{m}h{M}_0+{\epsilon }_{2}B{\Lambda }_a)\\ +{\beta }_s{\alpha }_{s}Bfg{\Lambda }_a{\Lambda }_s{\beta }_a{\beta }_h{\alpha }_h{\alpha }_a\\ +{\mu }_c{\mu }_{s}qfg{C}_0{\beta }_a{\alpha }_a{\mu }_h({\mu }_{m}h{M}_0+{\epsilon }_{2}B{\Lambda }_a)\\ +{\beta }_s{\alpha }_{s}Bfg{\Lambda }_a{\beta }_a{\alpha }_a{\mu }_h\\ \times (D{\Lambda }_s{\epsilon }_1+{\mu }_{c}q{C}_0),\end{array}$

$\begin{array}{l}Z={\mu }_c{\mu }_{s}q{C}_0{\mu }_{m}fgh{M}_0{\mu }_a{\mu }_h(1-R_0^2).(21)\end{array}$

Solving (20) for ${\lambda }_c^{*}$, it can be seen that one of the solutions is ${\lambda }_c^{*}=0$ which represents the disease-free equilibrium state, while the other solution of the polynomial (20) is the endemic equilibrium state of the model system (4) that is given by

${\lambda }_c^{*}=\frac{-Y+\sqrt{Y^2-4XZ}}{2X},$

and it is positive whenever R₀ > 1. We further use Descarte's rule of signs to discuss the existence of possible positive roots of Equation (18) especially when R₀ < 1. The results are summarized in Table 1. Hence, the following theorem.

TABLE 1

Table 1. No. of possible positive real roots for (20).

Theorem 3.3. System (4) has the following endemic equilibrium:

• If R₀ > 1, then the system has one positive unique endemic equilibrium state.

• If R₀ < 1, then the system has zero or two endemic equilibria.

In case (4), there is possibility of backward bifurcation of the model system (4) as it may contain endemic equilibrium with DFE when R₀ < 1. This is further determine by carrying out a bifurcation analysis in the following subsection.

3.6. Bifurcation analysis

The stability of endemic equilibrium state is established using Center Manifold Theory by Chavez et al. [26]. The theory shows the direction of bifurcation, whether it is a forward or backward bifurcation. A forward bifurcation indicates that the disease-free and endemic equilibrium states are locally asymptotically stable if R₀ < 1 and R₀ > 1, respectively, which implies the existence of global stability of the equilibrium states. A backward bifurcation means there is a co-existence of disease-free and endemic equilibrium when R₀ < 1. Therefore, we apply the theory in Chavez et al. [26] by re-writing the state variables as follows

$\begin{array}{l}{S}_h=x_1,I_A=x_2,I_{ch}=x_3,M=x_4,S_c=x_5,I_c=x_6,\\ S_s=x_7,I_s=x_8,C=x_9\end{array}$

such that the model system (4) can be rewritten as

$\frac{dX}{dt}=F(X)$

where $F={(f_1,f_2,f_3,f_4,f_5,f_6,f_7,f_8,f_9)}^T$ is the right hand side of model system (4).

Let the bifurcation parameter be ${\beta }_s={\beta }_s^{*}$ at R₀ = 1 so that the Jacobian matrix, J(E₀), of (15) will have a simple zero eigenvalue and negative real eigenvalues. The left and right eigenvectors associated with a simple zero eigenvalue of the Jacobian matrix (15) are w = (w₁, w₂, w₃, w₄, w₅, w₆, w₇, w₈, w₉) and v = (v₁, v₂, v₃, v₄, v₅, v₆, v₇, v₈, v₉), respectively, where

$\begin{array}{l}{w}_1=-\frac{({\mu }_h+k)}{{\mu }_h}{w}_2,w_3=\frac{k}{g}{w}_2,\\ w_4=-\frac{{\beta }_a{\alpha }_a{S}_a^{0}f}{{\mu }_a{\beta }_h{\alpha }_h{S}_h^0}{w}_2,w_5=\frac{f{\beta }_a{\alpha }_a{S}_a^0}{h{\beta }_h{\alpha }_h{S}_h^0}{w}_2,\\ w_6=-\frac{qf{C}_0{\mu }_c}{{\beta }_h{\alpha }_h{\mu }_{s}D{S}_h^0}{w}_2,w_7=\frac{f{C}_0{\mu }_c}{{\beta }_h{\alpha }_h{S}_h^{0}D}{w}_2,\\ w_8=\frac{q{M}_{0}f{C}_0{\mu }_c}{{\beta }_h{\beta }_s{\alpha }_h{\alpha }_{s}D{S}_h^0{S}_s^0}{w}_2,w_9=\frac{f{C}_0}{{\beta }_h{\alpha }_h{S}_h^0}{w}_2,\\ v_1=v_4=v_6=0,v_3=\frac{\sigma f}{(k\sigma +g)}{v}_2,\\ v_5=\frac{Bfg}{hA(k\sigma +g)}{v}_2,v_7=\frac{D}{q}{v}_9,v_8=\frac{fg}{A(k\sigma +g)}{v}_2,\end{array}$

$\begin{array}{l}{v}_9=\frac{{\beta }_h{\alpha }_h{S}_h^{0}hA(k\sigma +g)+Bfg{\beta }_a{\alpha }_a{S}_a^0}{C_{0}hA(k\sigma +g)}{v}_2.(22)\end{array}$

We state and prove the following theorem

Theorem 3.4. Model system (4) undergoes a forward bifurcation at R₀ = 1.

Proof. To determine the direction of the bifurcation, Theorem 4.1 in Castillo-Chavez and Song [27] is used to compute the bifurcation coefficients, m and n, of model system (4). The non-zero associated partial derivatives of F at the disease-free equilibrium, E₀, are

$\begin{array}{l}\frac{{\partial }^2{f}_2(E_0)}{\partial x_1\partial x_9}=\frac{{\beta }_h{\alpha }_h}{C_0},\frac{{\partial }^2{f}_2(E_0)}{\partial x_9^2}=-\frac{2{\beta }_h{\alpha }_h{\epsilon }_1{S}_h^0}{C_0^2},,\\ \frac{{\partial }^2{f}_5(E_0)}{\partial x_4\partial x_9}=\frac{{\beta }_a{\alpha }_a}{C_0},\frac{{\partial }^2{f}_5(E_0)}{\partial x_9^2}=-\frac{2{\beta }_a{\alpha }_a{\epsilon }_1{S}_a^0}{C_0^2},\\ \frac{{\partial }^2{f}_7(E_0)}{\partial x_6\partial x_8}=\frac{{\beta }_s{\alpha }_s}{M_0},\frac{{\partial }^2{f}_7(E_0)}{\partial x_8^2}=-\frac{2{\beta }_s{\alpha }_s{\epsilon }_2{S}_s^0}{M_0^2},\end{array}$

$\begin{array}{l}\frac{{\partial }^2{f}_7}{\partial {\beta }_s\partial x_8}=\frac{{\alpha }_s{S}_s^0}{M_0}.\end{array}$

Now applying the concept of Castillo-Chavez and Song [27], the value of m gives

$\begin{array}{l}m=v_2[w_1{w}_9\frac{{\partial }^2{f}_2(E_0)}{\partial x_1\partial x_9}+w_9^2\frac{{\partial }^2{f}_2(E_0)}{\partial x_9^2}]\\ +v_5[w_4{w}_9\frac{{\partial }^2{f}_5(E_0)}{\partial x_4\partial x_9}+w_9^2\frac{{\partial }^2{f}_5(E_0)}{\partial x_9^2}]\\ +v_7[w_6{w}_8\frac{{\partial }^2{f}_7(E_0)}{\partial x_6\partial x_8}+w_8^2\frac{{\partial }^2{f}_7(E_0)}{\partial x_8^2}].(23)\end{array}$

Upon substituting the eigenvectors and partial derivatives, we obtain

$\begin{array}{l}m=-v_2{w}_2^{2}G(24)\end{array}$

where

$\begin{array}{l}G=\frac{({\mu }_h+k)f}{{\mu }_h}+\frac{2{\epsilon }_1{f}^2}{{\beta }_h{\alpha }_h{S}_h^0}\\ +\frac{Bfg}{hA(k\sigma +g)}[\frac{{\beta }_a^2{\alpha }_a^2{f}^2{S}_a^{02}}{{\mu }_h{\beta }_h^2{\alpha }_h^2{S}_h^{02}}+\frac{2{\beta }_a{\alpha }_a{\epsilon }_1{f}^2{S}_a^{02}}{{\beta }_h^2{\alpha }_h^2{S}_h^{02}}]\\ +\left[\frac{{\beta }_h{\alpha }_h{S}_h^{0}hA(k\sigma +g)+Bfg{\beta }_a{\alpha }_a{S}_c^0}{hA(k\sigma +g)}\right]\\ [\frac{q{C}_0{f}^2{\mu }_c^2{M}_0}{{\mu }_s{\beta }_h^2{\alpha }_h^{2}f{S}_h^{02}}+\frac{2{\epsilon }_{2}f{\mu }_c}{M_0}].(25)\end{array}$

Also for n, we have

$\begin{array}{l}n=v_7{w}_8\frac{{\partial }^2{f}_7(E_0)}{\partial {\beta }_s\partial x_8}\\ =\frac{f{\mu }_c}{{\beta }_h{\beta }_s{\alpha }_h{\alpha }_s{S}_s^0}\left[\frac{{\beta }_h{\alpha }_{h}hA(k\sigma +g)S_h^0+Bfg{\beta }_a{\alpha }_a{S}_a^0}{hA(k\sigma +g)}\right]v_2{w}_2>0.(26)\end{array}$

A forward bifurcation exists since m < 0 in (24) and n > 0 in (26). Therefore, a forward bifurcation exists for the model (4) at R₀ = 1. This completes the proof.

With the existence of forward bifurcation, a unique endemic equilibrium, E₁, exists and is locally asymptotically stable whenever R₀ > 1 and unstable when R₀ < 1. Hence, there is possibility of global stability of the equilibrium states. A forward bifurcation plot is shown in Figure 2.

FIGURE 2

Figure 2. Bifurcation plot for the dynamics of Schistosomiasis.

3.7. Global stability of the model

In this subsection, the global behavior of the disease-free and endemic equilibria of the system (4) is established by constructing a suitable Lyapunov functions and applying LaSalle's invariance principle.

3.7.1. Global stability of disease-free equilibrium, E₀

Theorem 3.5. The disease-free equilibrium, E₀, of the schistosomiasis model system (4) is globally asymptotically stable when R₀ < 1 otherwise unstable.

Proof. Constructing a Lyapunov function using a matrix-theoretic method by [28] gives

$\begin{array}{l}L=\frac{A(k\sigma +g)}{{\mu }_{m}fg}{I}_A+\frac{A\sigma }{g{\mu }_m}{I}_{ch}+\frac{B}{h{\mu }_m}{I}_a+\frac{D}{q{\mu }_c}{I}_s+\frac{1}{{\mu }_m}M+\frac{1}{{\mu }_c}C(27)\end{array}$

where the letters, A, B, and D are defined in (12).

Taking the derivative of Equation (27) along the trajectory yields

$\begin{array}{l}{L}^{\prime }=\frac{A(k\sigma +g)}{{\mu }_{m}fg}(\frac{{\beta }_h{\alpha }_{h}C{S}_h}{C_0+{\epsilon }_{1}C}-f{I}_A)+\frac{A\sigma }{g{\mu }_m}(k{I}_A-g{I}_{ch})\\ +\frac{B}{h{\mu }_m}(\frac{{\beta }_a{\alpha }_{a}C{S}_a}{C_0+{\epsilon }_{1}C}-h{I}_a)\\ +\frac{D}{q{\mu }_c}(\frac{{\beta }_s{\alpha }_s{S}_{s}M}{M_0+{\epsilon }_{2}M}-q{I}_s)\\ +\frac{1}{{\mu }_m}(A{I}_A+A\sigma I_{ch}+B{I}_a-{\mu }_{m}M)+\frac{1}{{\mu }_c}(D{I}_s-{\mu }_{c}C).\left(28\right)\end{array}$

Expanding Equation (28) and simplifying in terms of reproduction numbers (R_A, R_ch, R_m, R_s) gives

$\begin{array}{l}{L}^{\prime }=(R_A+R_{ch}+R_m-1)C+(R_s-1)M\\ -\frac{{\beta }_h{\alpha }_{h}A(k\sigma +g)}{{\mu }_{m}fg}\left(\frac{(S_h^0-S_h)C_0+{\epsilon }_1{S}_h^{0}C}{C_0(C_0+{\epsilon }_{1}C)}\right)\\ -\frac{B{\beta }_a{\alpha }_a}{h{\mu }_m}\left(\frac{(S_a^0-S_a)C_0+{\epsilon }_1{S}_a^{0}C}{C_0(C_0+{\epsilon }_{1}C)}\right)C\\ -\frac{D{\beta }_s{\alpha }_s}{q{\mu }_c}\left(\frac{(S_s^0-S_s)M_0+{\epsilon }_2{S}_s^{0}M}{M_0(M_0+{\epsilon }_{2}M)}\right)M.(29)\end{array}$

With $S_h^0\ge S_h$, $S_a^0\ge S_a$, $S_s^0\ge S_s$, we have

$L^{\prime }\le (R_A+R_{ch}+R_m-1)C+(R_s-1)M.$

Since R_A+R_ch+R_m ≤ 1, R_s ≤ 1 from the definition of the reproduction number in (13) and (14), this yields L′ ≤ 0 if R₀ < 1. Hence by LaSalle's invariance principle, the disease-free equilibrium is globally asymptotically stable if R₀ ≤ 1. This completes the proof.

3.7.2. Global stability of endemic equilibrium, E₁

To prove the stability of the endemic equilibrium (E₁) we state and prove the following theorem

Theorem 3.6. The endemic equilibrium, E₁, of the schistosomiasis model system (4) is globally asymptotically stable whenever R₀ > 1 and ψ = 0.

Proof. The global stability of the endemic equilibrium is proved by constructing the Volterra-Lyapunov function

$\begin{array}{l}{H}=\frac{I_A^{*}+\sigma I_{ch}^{*}}{\sigma I_{ch}^{*}}\left[S_h-S_h^{*}-S_h^{*}ln\frac{S_h}{S_h^{*}}\right]\\ +\frac{I_A^{*}+\sigma I_{ch}^{*}}{\sigma I_{ch}^{*}}\left[I_A-I_A^{*}-I_A^{*}ln\frac{I_A}{I_A^{*}}\right]+\frac{{\beta }_h{\alpha }_h{S}_h^{*}{\lambda }_c^{*}}{k{I}_A^{*}}\\ \times \left[I_{ch}-I_{ch}^{*}-I_{ch}^{*}ln\frac{I_{ch}}{I_{ch}^{*}}\right]\\ +\frac{B{I}_a^{*}{\beta }_h{\alpha }_h{S}_h^{*}{\lambda }_m^{*}}{{\beta }_a{\alpha }_a{S}_a^{*}{\lambda }_c^{*}A\sigma I_{ch}^{*}}\left[S_a-S_a^{*}-S_a^{*}ln\frac{S_a}{S_a^{*}}\right]\\ +\frac{B{I}_a^{*}{\beta }_h{\alpha }_h{S}_h^{*}{\lambda }_m^{*}}{{\beta }_a{\alpha }_a{S}_a^{*}{\lambda }_c^{*}A\sigma I_{ch}^{*}}\\ \times \left[I_a-I_a^{*}-I_a^{*}ln\frac{I_a}{I_a^{*}}\right]+\left[S_s-S_s^{*}-S_s^{*}ln\frac{S_s}{S_s^{*}}\right]\\ +\left[I_s-I_s^{*}-I_s^{*}ln\frac{I_s}{I_s^{*}}\right]+\frac{{\beta }_h{\alpha }_h{S}_h^{*}{\lambda }_c^{*}}{A\sigma I_{ch}^{*}}\\ \times \left[M-M^{*}-M^{*}ln\frac{M}{M^{*}}\right]\end{array}$

$\begin{array}{l}+\frac{{\beta }_s{\alpha }_s{\lambda }_m^{*}{S}_s^{*}}{D{I}_s^{*}}\left[C-C^{*}-C^{*}ln\frac{C}{C^{*}}\right],(30)\end{array}$

where ${\lambda }_c^{*}$, ${\lambda }_m^{*}$, A, B, and D are defined in (12) and (18).

Taking the time derivative of ${H}$ along the solutions of the schistosomiasis model Equation (4) yields

$\begin{array}{l}{{H}}^{\prime }=\frac{I_A^{*}+\sigma I_{ch}^{*}}{\sigma I_{ch}^{*}}\left[1-\frac{S_h^{*}}{S_h}\right]({\Lambda }_h-{\beta }_h{\alpha }_h{S}_h{\lambda }_c-{\mu }_h{S}_h+\psi I_A)\\ +\frac{I_A^{*}+\sigma I_{ch}^{*}}{\sigma I_{ch}^{*}}\left[1-\frac{I_A^{*}}{I_A}\right]({\beta }_h{\alpha }_h{S}_h{\lambda }_c-f{I}_A)\\ +\frac{{\beta }_h{\alpha }_h{S}_h^{*}{\lambda }_c^{*}}{k{I}_A^{*}}\left[1-\frac{I_{ch}^{*}}{I_{ch}}\right](k{I}_A-g{I}_{ch})\\ +\frac{B{I}_a^{*}{\beta }_h{\alpha }_h{S}_h^{*}{\lambda }_m^{*}}{{\beta }_a{\alpha }_a{S}_a^{*}{\lambda }_c^{*}A\sigma I_{ch}^{*}}\left[1-\frac{S_a^{*}}{S_a}\right]({\Lambda }_a-{\beta }_a{\alpha }_a{\lambda }_c{S}_a-{\mu }_a{S}_a)\\ +\frac{B{I}_a^{*}{\beta }_h{\alpha }_h{S}_h^{*}{\lambda }_m^{*}}{{\beta }_a{\alpha }_a{S}_a^{*}{\lambda }_c^{*}A\sigma I_{ch}^{*}}\left[1-\frac{I_a^{*}}{I_a}\right]({\beta }_a{\alpha }_a{\lambda }_c{S}_a-h{I}_a)\\ +\left[1-\frac{S_s^{*}}{S_s}\right]({\Lambda }_s-{\beta }_s{\alpha }_s{\lambda }_m{S}_s-{\mu }_s{S}_s)+\left[1-\frac{I_s^{*}}{I_s}\right]\\ \times ({\beta }_s{\alpha }_s{\lambda }_m{S}_s-q{I}_s)\\ +\frac{{\beta }_h{\alpha }_h{S}_h^{*}{\lambda }_c^{*}}{A\sigma I_{ch}^{*}}\left[1-\frac{M^{*}}{M}\right](A{I}_A+A\sigma I_{ch}\\ +B{I}_a-{\mu }_{m}M)+\frac{{\beta }_s{\alpha }_s{\lambda }_m^{*}{S}_s^{*}}{D{I}_s^{*}}.\end{array}$

$\begin{array}{l}\times \left[1-\frac{C^{*}}{C}\right](D{I}_s-{\mu }_{c}C).(31)\end{array}$

At endemic equilibrium state, we have

$\begin{array}{l}{\Lambda }_h={\beta }_h{\alpha }_h{S}_h^{*}{\lambda }_c^{*}+{\mu }_h{S}_h^{*}-\psi I_A^{*},f=\frac{{\beta }_h{\alpha }_h{S}_h^{*}{\lambda }_c^{*}}{I_A^{*}},\\ g=\frac{k{I}_A^{*}}{I_{ch}^{*}},{\Lambda }_c={\beta }_a{\alpha }_a{\lambda }_c^{*}{S}_a^{*}+{\mu }_a{S}_a^{*},\\ h=\frac{{\beta }_a{\alpha }_a{\lambda }_c^{*}{S}_a^{*}}{I_a^{*}},{\Lambda }_s={\beta }_s{\alpha }_s{\lambda }_m^{*}{S}_s^{*}+{\mu }_s{S}_s^{*},\\ q=\frac{{\beta }_s{\alpha }_s{\lambda }_m^{*}{S}_s^{*}}{I_s^{*}},{\mu }_m=\frac{A{I}_A^{*}}{M^{*}}+\frac{A\sigma I_{ch}^{*}}{M^{*}}+\frac{B{I}_a^{*}}{M^{*}},\end{array}$

$\begin{array}{l}{\mu }_c=\frac{D{I}_s^{*}}{C^{*}}.(32)\end{array}$

Substituting Equation (32) into (31) and expanding and simplifying gives

$\begin{array}{l}{H}^{\prime }=\frac{{\mu }_h{S}_h^{*}(I_A^{*}+\sigma I_{ch}^{*})}{\sigma I_{ch}^{*}}\left[2-\frac{S_h^{*}}{S_h}-\frac{S_h}{S_h^{*}}\right]\\ +\frac{B{I}_a^{*}{\beta }_h{\alpha }_h{\lambda }_m^{*}{S}_h^{*}{\mu }_a}{{\beta }_a{\alpha }_{a}A\sigma {\lambda }_c^{*}{I}_{ch}^{*}}\left[2-\frac{S_a^{*}}{S_a}-\frac{S_a}{S_a^{*}}\right]\\ +{\mu }_s{S}_s^{*}\left[2-\frac{S^{*}}{S_s}-\frac{S_s}{S_s^{*}}\right]\\ +\psi I_A^{*}\left(\frac{I_A^{*}+\sigma I_{ch}^{*}}{\sigma I_{ch}^{*}}\right)\left[\frac{S_h^{*}}{S_h}+\frac{I_A}{I_A^{*}}-1-\frac{I_A}{I_A^{*}}\frac{S_h^{*}}{S_h}\right]\\ +\frac{{\beta }_h{\alpha }_h{S}_h^{*}{\lambda }_c^{*}{I}_A^{*}}{\sigma I_{ch}^{*}}[3+\frac{{\lambda }_c}{{\lambda }_c^{*}}-\frac{S_h^{*}}{S_h}-\frac{S_h}{S_h^{*}}\frac{{\lambda }_c}{{\lambda }_c^{*}}\frac{I_A^{*}}{I_A}-\frac{M}{M^{*}}\\ -\frac{I_A}{I_A^{*}}\frac{M^{*}}{M}]+{\beta }_h{\alpha }_h{S}_h^{*}{\lambda }_c^{*}[4+\frac{{\lambda }_c}{{\lambda }_c^{*}}-\frac{S_h^{*}}{S_h}-\frac{S_h}{S_h^{*}}\frac{{\lambda }_c}{{\lambda }_c^{*}}\frac{I_A^{*}}{I_A}\\ -\frac{I_A}{I_A^{*}}\frac{I_{ch}^{*}}{I_{ch}}-\frac{M}{M^{*}}-\frac{I_{ch}}{I_{ch}}\frac{M^{*}}{M}]+\frac{B{I}_a^{*}{\beta }_h{\alpha }_h{\lambda }_m^{*}{S}_h^{*}}{A\sigma I_{ch}^{*}}\\ \times \left[3+\frac{{\lambda }_c}{{\lambda }_c^{*}}-\frac{S_a^{*}}{S_a}-\frac{S_a}{S_a^{*}}\frac{{\lambda }_c}{{\lambda }_c^{*}}\frac{I_a^{*}}{I_a}-\frac{M}{M^{*}}-\frac{I_a}{I_a^{*}}\frac{M^{*}}{M}\right]\\ +{\beta }_s{\alpha }_s{\lambda }_m^{*}{S}_s^{*}[3+\frac{{\lambda }_m}{{\lambda }_m^{*}}-\frac{S_s^{*}}{S_s}-\frac{S_s}{S_s^{*}}\frac{{\lambda }_m}{{\lambda }_m^{*}}\frac{I_s^{*}}{I_s}-\frac{C}{C^{*}}\\ -\frac{C^{*}}{C}\frac{I_s}{I_s^{*}}].(33)\end{array}$

Using the assumption that ψ = 0 and $\frac{{\lambda }_c}{{\lambda }_c^{*}}$, $\frac{{\lambda }_m}{{\lambda }_m^{*}}\le 1$, we have by Arithmetic-Geometric Theorem that

$\begin{array}{l}2-\frac{S_h^{*}}{S_h}-\frac{S_h}{S_h^{*}}\le 0,2-\frac{S_a^{*}}{S_a}-\frac{S_a}{S_a^{*}}\le 0,2-\frac{S_s^{*}}{S_s}-\frac{S_s}{S_s^{*}}\le 0,\\ 3+\frac{{\lambda }_c}{{\lambda }_c^{*}}-\frac{S_h^{*}}{S_h}-\frac{S_h}{S_h^{*}}\frac{{\lambda }_c}{{\lambda }_c^{*}}\frac{I_A^{*}}{I_A}-\frac{M}{M^{*}}-\frac{I_A}{I_A^{*}}\frac{M^{*}}{M}\le 4-\frac{S_h^{*}}{S_h}\\ -\frac{S_h}{S_h^{*}}\frac{I_A^{*}}{I_A}-\frac{M}{M^{*}}-\frac{I_A}{I_A^{*}}\frac{M^{*}}{M}\le 0,\\ 4+\frac{{\lambda }_c}{{\lambda }_c^{*}}-\frac{S_h^{*}}{S_h}-\frac{S_h}{S_h^{*}}\frac{{\lambda }_c}{{\lambda }_c^{*}}\frac{I_A^{*}}{I_A}-\frac{I_A}{I_A^{*}}\frac{I_{ch}^{*}}{I_{ch}}-\frac{M}{M^{*}}-\frac{I_{ch}}{I_{ch}^{*}}\frac{M^{*}}{M}\le 5\\ -\frac{S_h^{*}}{S_h}-\frac{S_h}{S_h^{*}}\frac{I_A^{*}}{I_A}-\frac{I_A}{I_A^{*}}\frac{I_{ch}^{*}}{I_{ch}}-\frac{M}{M^{*}}-\frac{I_{ch}}{I_{ch}}\frac{M^{*}}{M}\le 0,\\ 3+\frac{{\lambda }_c}{{\lambda }_c^{*}}-\frac{S_a^{*}}{S_a}-\frac{S_a}{S_a^{*}}\frac{{\lambda }_c}{{\lambda }_c^{*}}\frac{I_a^{*}}{I_a}-\frac{M}{M^{*}}-\frac{I_a}{I_a^{*}}\frac{M^{*}}{M}\le 4-\frac{S_a^{*}}{S_a}\\ -\frac{S_a}{S_a^{*}}\frac{I_a^{*}}{I_a}-\frac{M}{M^{*}}-\frac{I_a}{I_a^{*}}\frac{M^{*}}{M}\le 0,\\ 3+\frac{{\lambda }_m}{{\lambda }_m^{*}}-\frac{S_s^{*}}{S_s}-\frac{S_s}{S_s^{*}}\frac{{\lambda }_m}{{\lambda }_m^{*}}\frac{I_s^{*}}{I_s}-\frac{C}{C^{*}}-\frac{C^{*}}{C}\frac{I_s}{I_s^{*}}\le 4-\frac{S_s^{*}}{S_s}\end{array}$

$\begin{array}{l}-\frac{S_s}{S_s^{*}}\frac{I_s^{*}}{I_s}-\frac{C}{C^{*}}-\frac{C^{*}}{C}\frac{I_s}{I_s^{*}}\le 0.(34)\end{array}$

This implies that L′ ≤ 0 in the region, Z. Applying LaSalle's invariance principle, the endemic equilibrium state of the schistosomiasis model system (4) is globally asymptotically stable in Z if R₀ > 1 and ψ = 0. This means that with different initial conditions at ψ = 0, the global stability of the schistosomiasis model system (4) will always converge to the disease-free and endemic equilibrium states whenever R₀ < 1 and R₀ > 1, respectively. The numerical simulations are carried out in the next section.

4. Sensitivity analysis

4.1. Parameter values

To perform numerical simulations, parameters for the model are mainly sourced from literature, as shown in Table 2. Parameters that are represented by probabilities are assumed to be within the range of [0 − 1]. Any parameters that are outside the prescribed range are properly sampled from a uniform distribution within ranges defined in Table 2.

TABLE 2

Table 2. Parameter values used for the numerical simulation of model system (4).

4.2. Sensitivity analysis results

Sensitivity analysis examines the relative change of the mathematical model output (variable) when the model parameter(s) are altered or changed. It helps in determining and identifying the model parameters that require special attention. To carry out the sensitivity analysis of the basic reproduction number R₀, we adopt the Latin Hypercube Sampling (LHS) scheme and Partial Rank Correlation Coefficients (PRCCs) technique used by Blower and Dowlatabadi [39] to compute and identify the biological implication of each model parameter on R₀. See also the following published articles [40–45] with similar analyses for different epidemic models where the LHS techniques have been applied. We performed 1,000 simulations per run and examined and evaluated the PRCC of the model parameters concerning R₀. The PRCC values indicate the degree of monotonicity between the parameters of the model and R₀. Thus, juxtaposing the PRCC values gives us a clear picture on the contribution of each parameter of the model on R₀. Table 3 and the Tornado plot of Figure 3 give the model parameters with their PRCCs. The parameters with the potential of increasing(decreasing) the value of R₀, are those with positive(negative) PRCCs. Thereby, increasing parameters, β_s, β_h, β_a, γ_h, γ_a, γ_s, τ₁, τ₂, increases the value of R₀ which in turn increase the spread of schistosomiasis disease within the population (see Table 2 and Figure 3). On the other hand, parameters μ_c, μ_m, μ_s with negative PRCCs influence the reduction of the burden of schistosomiasis disease within the population. When these parameters are increased, we observe a reduction in the value of R₀, hence the endemicity of the Schistosomiasis within the population is also reduced. The natural death rate for the snail population is the most sensitive parameter with respect to R₀ with a PRCC value of 0.776046728. The implication of this is that increasing snail death rate reduces the spread of Schistosomiasis within the population. To examine the significance of the model parameters, we compute the p-value of PRCC of each model parameter using the Fisher Transformation method as shown in Table 3. Table 3 reveals that the parameters, β_s, β_h, γ_h, γ_s, μ_c, μ_m, μ_s, τ₁, τ₂ have significant p-values while the parameters, β_a, γ_a, have insignificant p-values. This implies that the shedding and transmission rates of the other mammals, γ_a and β_a, are non-monotonically related to R₀, though they can still produce changes in the transmission dynamics of schistosomiasis infection.

TABLE 3

Table 3. PRCC values and significance.

FIGURE 3

Figure 3. Tornado plot for the dynamics of Schistosomiasis infection on the R₀ showing the model parameters with their PRCC values.

The pairwise comparison of the significant parameters of Table 3, whose p < 0.05 is given in Tables 5, 6. The pairwise comparison is carried out to observe the model parameters with notable influence on increasing the burden of the schistosomiasis infection in the population, as well as to establish whether there exists any difference between the processes describing the compared parameters. We presented the outcome of the pairwise PRCC comparison between the unadjusted p-values (Table 4) and the false discovery rate (FDR) adjusted p-values in Tables 5, 6. According to Table 6, whenever the p-values of the compared pair of significant parameters are < 0.05, it is interpreted to be significant (TRUE) otherwise it is insignificant (FALSE). It is also important to note from Table 6 that the pairs with significantly different (TRUE) are more sensitive parameters.

TABLE 4

Table 4. Pairwise PRCC comparisons (Unadjusted P-values).

TABLE 5

Table 5. Pairwise PRCC comparisons (FDR adjusted P-values).

TABLE 6

Table 6. Parameters different after FDR adjustment.

The effect of the most sensitive parameters on R₀ is presented in Figures 3, 4 as scatter plots to support results in Table 2. These plots show how parameters, β_s, γ_s, γ_h, τ₁, μ_c, μ_m, μ_s, τ₂ monotonically increase or decrease R₀. The most noticeable result is the scatter plot for snail death μ_s which is the most sensitive parameter. The sensitivity result reveals that the spread of schistosomiasis infection will be controlled and eradicated in the population if more snails die. If the Miracidia and Cercariae (μ_m, μ_c) are cleared within the population, we will observe a reduction in the contact rates of human to Cercariae (β_h), snail to Miracidia (β_s), probability of egg developing into a Miracidium (τ₁) and the density of Cercariae in the environment (τ₂).

FIGURE 4

Figure 4. Monte Carlo simulations for the eight parameters that have PRCC values that are significant as evidenced by Table 3. Parameter values were taken from Table 2. Each run consists of 1,000 simulations of randomly drawn parameters.

5. Numerical simulation results

5.1. Time series trajectories

To support the analytical results presented in this paper, we present numerical simulations in the form of time series plots and phase plots on the dynamics of Schistosomiasis. The simulations presented are produced using MATLAB software. Figures 5A,B, present the time series plots for the acute infected class I_A(t), chronic infected class, I_ch(t), infected animals I_a(t), and infected snails, I_a(t), while varying the snail death rate, μ_s. The snail death rate, μ_s is considered for simulation as it is the most sensitive parameter. It is observed that when the death rate of the snail population is increasing, the peak of I_A(t), I_ch(t), and I_a(t), shifts to the right in a reduced manner. On the other hand, Figures 6A,B present the dynamics of Miracidia, M(t) and Cercariae, C(t). We observe that an increase in the snail death rates leads to a decrease in the population of the Miracidia. Also, the populations of I_s(t) and C(t) reduce drastically as the snail population's death rate increases. This indicates that reducing the snail population reduces the time it will take to infect human and mammal populations, including Cercariae that will be cleared in the population. This eventually reduces the number of Schistosomiasis cases within the population.

FIGURE 5

Figure 5. (A) Dynamics of Schistosomiasis infection among the infected human classes, I_a(t) and I_ch(t) for varying snail death rate, μ_s. (B) Dynamics of Schistosomiasis infection among the infected classes, I_a(t) and I_s(t) for varying snail death rate, μ_s. All other parameters for this simulation are given in Table 2.

FIGURE 6

Figure 6. Dynamics of Schistosomiasis infection among: (A) the Miracidia, (M(t)) and (B) the Cercariae classes, (C(t)), for varying snail death rate, μ_s. All other parameters for this simulation are given in Table 2.

5.2. Phase portraits

Figures 7, 8 present the phase portraits of the infection sub-populations when R₀ > 1 and R₀ < 1, respectively. When R₀ < 1, it can be seen that solutions tend to disease-free equilibrium over time. This means that Schistosomiasis infection can be contained within the population over time provided that R₀ < 1. On the other hand when R₀ > 1, the solutions lead to an endemic equilibrium as seen by a decrease in the susceptible population and an increase in the infected population. It is also observed that as Miracidia initially increase, Cercariae increases as well and reaches a peak. Thereafter, Miracidia decreases gradually as Cercariae continues to increase. This means that Schistosomiasis infection will persist within the population whenever R₀ > 1.

FIGURE 7

Figure 7. Phase portraits for the dynamics of Schistosomiasis infection when R₀ < 1. All other parameters in Table 2.

FIGURE 8

Figure 8. Phase portraits for the dynamics of Schistosomiasis infection when R₀ > 1. All other parameters in Table 2.

5.3. Effects of parameters on R₀

Using contour plots, Figures 9A,B depict that R₀ attains a lower value when the probability of eggs developing into Miracidium is minimal and the snail death rate is high as shown by the contour plot and the 3-D representation. Figure 9C shows R₀ vs. varying the natural decay rates of Miracidium, (μ_m) and Cercariae (μ_c). It can be seen that if the decay rates are effectively increased, the reproduction number, R₀ can be reduced to below one, implying that the disease-free equilibrium can be achieved. The implication of Figures 9A–C is that there is a possibility of approaching a disease-free equilibrium if these parameters are controlled accordingly. However, this may be difficult to attain due to the lack of proper control measures that can simultaneously increase the decay rates. Generally, it can be seen that the control of snails is very important in the reduction of the spread of Schistosomiasis within the population. Therefore, measures need to be put in place so that all breeding areas for snails are identified and dealt with appropriately without harming other aquatic organisms.

FIGURE 9

Figure 9. Contour plots showing the effects of (A) τ₁ vs. μ_s (B) τ₂ vs. μ_s (C) μ_c vs. μ_m.

6. Conclusion

In this paper, a deterministic mathematical model of transmission dynamics of Schistosomiasis disease that incorporates acute and chronic infected humans is formulated. The main objective of the model is to study the impact of snail control on the dynamics of Schistosomiasis. The model analysis comprised of establishing the invariant region, the positivity of solutions, the existence of disease-free and endemic equilibrium, bifurcation analysis and computation of the basic reproduction number, R₀ using the next-generation approach. We proved that whenever R₀ < 1, the disease-free equilibrium is globally asymptotically stable, which means that the disease will fizzle out within the population over time. The endemic equilibrium was found to be globally asymptotically stable implying that if not controlled, Schistosomiasis will persist within the population. We also carried out a sensitivity analysis of the parameters governing the basic reproduction number, R₀, using the Latin Hypercube Sampling (LHS) scheme-Partial Rank Correlation Coefficient (PRCC) technique. The results of the sensitivity analysis revealed that the spread of schistosomiasis infection can be controlled and eradicated in the population by increasing the natural death rate of the snail population. It also established that there is need to increase the decay rates of Miracidia and Cercariae within the environment by applying a control measure. Therefore, the results of our study effectively conclude that the most effective method to control the transmission dynamics and spread of schistosomiasis disease is by rapidly killing of infected snails and applying control measures that will clear the parasite, miracidia and cercariae, in the environment.

One limitation of our model is that it was not fitted to real-life data. Instead, our model was parameterized using accurate parameter estimates from existing literature on Schistosomiasis dynamics with inclusive heterogeneity in the disease transmission. Secondly, the assumption of mammals could be made specific by considering only cattle in the modeling framework. To bypass this complexity, we carried out a global sensitivity analysis to vary the parameter values and establish the uncertainty in the model parameters. Notwithstanding this limitation, we maintain that our findings remain valid and applicable if followed by policymakers to help eradicate the spread of Schistosomiasis.

Future work may look at an optimal control analysis of the model presented in the presence of different control strategies. The model may also be extended to include infection dynamics in the presence of snail age stratification.

Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.

Author contributions

CM and IO contributed to conception and model formulation of the study. CM, ZC, IO, CC, and FF contributed in the analysis and write-up of sections of the manuscript. All authors contributed to manuscript revision, read, and approved the final version of the manuscript.

Acknowledgments

The authors would like to thank their respective universities and the reviewers for their support toward this manuscript.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. WHO. WHO Guideline on Control and Elimination of Human Schistosomiasis. WHO (2022).

Google Scholar

2. Ross AG, Vickers D, Olds GR, Shah SM, McManus DP. Katayama syndrome. Lancet Infect Dis. (2007) 7:218–24. doi: 10.1016/S1473-3099(07)70053-1

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Shiff C. The importance of definitive diagnosis in chronic schistosomiasis, with reference to Schistosoma haematobium. J Parasitol Res. (2012) 2012:761269. doi: 10.1155/2012/761269

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Costa de Limae M, Pocha R, Filhocoura P, Katz N. A 13-years follow-up of treatment and snail control in the endemic area for Schistosoma mansoni in Brazil, incidence of infection and re-infection. Bull World Health Organ. (1993) 71:197–205.

PubMed Abstract | Google Scholar

5. Guiro A, Ouaro S, Traore A. Stability analysis of a schistosomiasis model with delays. Adv Differ Equat. (2013) 2013:1–15. doi: 10.1186/1687-1847-2013-303

CrossRef Full Text | Google Scholar

6. Chen Z, Zou L, Shen D, Zhang W, Ruan S. Mathematical modelling and control of Schistosomiasis in Hubei Province, China. Acta Trop. (2010) 115:119–25. doi: 10.1016/j.actatropica.2010.02.012

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Sokolow SH, Wood CL, Jones IJ, Swartz SJ, Lopez M, Hsieh MH, et al. Global assessment of schistosomiasis control over the past century shows targeting the snail intermediate host works best. PLoS Neglect Trop Dis. (2016) 10:e0004794. doi: 10.1371/journal.pntd.0004794

PubMed Abstract | CrossRef Full Text | Google Scholar

8. McCullough F, Gayral P, Duncan J, Christie J. Molluscicides in schistosomiasis control. Bull World Health Organ. (1980) 58:681.

Google Scholar

9. Lardans V, Dissous C. Snail control strategies for reduction of schistosomiasis transmission. Parasitol Today. (1998) 14:413–7. doi: 10.1016/S0169-4758(98)01320-9

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Zhang H, Harvim P, Georgescu P. Preventing the spread of schistosomiasis in Ghana: possible outcomes of integrated optimal control strategies. J Biol Syst. (2017) 25:625–55. doi: 10.1142/S0218339017400058

CrossRef Full Text | Google Scholar

11. Zhang T, Zhao XQ. Mathematical modeling for schistosomiasis with seasonal influence: a case study in Hubei, China. SIAM J Appl Dyn Syst. (2020) 19:1438–71. doi: 10.1137/19M1280259

CrossRef Full Text | Google Scholar

12. Bakare E, Nwozo C. Mathematical analysis of malaria-schistosomiasis coinfection model. Epidemiol Res Int. (2016) 2016:3854902. doi: 10.1155/2016/3854902

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Akanni JO, Olaniyi S, Akinpelu FO. Global asymptotic dynamics of a nonlinear illicit drug use system. J Appl Math Comput. (2021) 66:39–60. doi: 10.1007/s12190-020-01423-7

CrossRef Full Text | Google Scholar

14. Mustapha UT, Musa SS, Lawan MA, Abba A, Hincal E, Mohammed MD, et al. Mathematical modeling and analysis of schistosomiasis transmission dynamics. Int J Model Simul Sci Comput. (2021) 12:2150021. doi: 10.1142/S1793962321500215

CrossRef Full Text | Google Scholar

15. Wu YF, Li MT, Sun GQ. Asymptotic analysis of schistosomiasis persistence in models with general functions. J Frankl Instit. (2016) 353:4772–84. doi: 10.1016/j.jfranklin.2016.09.012

CrossRef Full Text | Google Scholar

16. Ronoh M, Chirove F, Pedro SA, Tchamga MSS, Madubueze CE, Madubueze SC, et al. Modelling the spread of schistosomiasis in humans with environmental transmission. Appl Math Modell. (2021) 95:159–75. doi: 10.1016/j.apm.2021.01.046

CrossRef Full Text | Google Scholar

17. Qi L, Cui JA, Gao Y, Zhu H. Modeling the schistosomiasis on the islets in Nanjing. Int J Biomath. (2012) 5:1250037. doi: 10.1142/S1793524511001799

CrossRef Full Text | Google Scholar

18. Abimbade SF, Olaniyi S, Ajala OA. Recurrent malaria dynamics: insight from mathematical modelling. Eur Phys J Plus. (2022) 137:1–16. doi: 10.1140/epjp/s13360-022-02510-3

CrossRef Full Text | Google Scholar

19. Hethcote HW. The mathematics of infectious diseases. SIAM Rev. (2000) 42:599–653. doi: 10.1137/S0036144500371907

CrossRef Full Text | Google Scholar

20. Acheampong E, Okyere E, Iddi S, Bonney JH, Asamoah JKK, Wattis JA, et al. Mathematical modelling of earlier stages of COVID-19 transmission dynamics in Ghana. Results Phys. (2022) 34:105193. doi: 10.1016/j.rinp.2022.105193

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Birkhoff G, Rota GC. Ordinary Differential Equations. New York, NY: John Wiley & Sons (1978).

Google Scholar

22. Van den Driessche P, Watmough J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Math Biosci. (2002) 180:29–48. doi: 10.1016/S0025-5564(02)00108-6

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Diekmann O, Heesterbeek JAP, Metz JA. On the definition and the computation of the basic reproduction ratio R 0 in models for infectious diseases in heterogeneous populations. J Math Biol. (1990) 28:365–82. doi: 10.1007/BF00178324

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Kim JH, Su W, Song YJ. On stability of a polynomial. J Appl Math Inform. (2018) 36:231–6. doi: 10.14317/jami.2018.231

CrossRef Full Text | Google Scholar

25. Heffernan JM, Smith RJ, Wahl LM. Perspectives on the basic reproductive ratio. J R Soc Interface. (2005) 2:281–93. doi: 10.1098/rsif.2005.0042

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Chavez CC, Feng Z, Huang W. On the computation of R0 and its role on global stability. In: Castillo-Chavex C Blower S Driessche P Kirschner D Yakubu A-A, editors. Mathematical Approaches for Emerging and Re-emerging Infection Diseases: An Introduction, Vol. 125. New York, NY: Springer (2002) p. 31–65. doi: 10.1007/978-1-4757-3667-0_13

CrossRef Full Text | Google Scholar

27. Castillo-Chavez C, Song B. Dynamical models of tuberculosis and their applications. Math Biosci Eng. (2004) 1:361. doi: 10.3934/mbe.2004.1.361

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Shuai Z, van den Driessche P. Global stability of infectious disease models using Lyapunov functions. SIAM J Appl Math. (2013) 73:1513–32. doi: 10.1137/120876642

CrossRef Full Text | Google Scholar

29. Kalinda C, Mushayabasa S, Chimbari MJ, Mukaratirwa S. Optimal control applied to a temperature dependent schistosomiasis model. Biosystems. (2019) 175:47–56. doi: 10.1016/j.biosystems.2018.11.008

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Feng Z, Eppert A, Milner FA, Minchella D. Estimation of parameters governing the transmission dynamics of schistosomes. Appl Math Lett. (2004) 17:1105–12. doi: 10.1016/j.aml.2004.02.002

CrossRef Full Text | Google Scholar

31. Chiyaka ET, GARIRA W. Mathematical analysis of the transmission dynamics of schistosomiasis in the human-snail hosts. J Biol Syst. (2009) 17:397–423. doi: 10.1142/S0218339009002910

CrossRef Full Text | Google Scholar

32. Abokwara A, Madubueze CE. The role of non-pharmacological interventions on the dynamics of schistosomiasis. J Math Fund Sci. (2021) 53:243–260. doi: 10.5614/j.math.fund.sci.2021.53.2.6

CrossRef Full Text | Google Scholar

33. Kanyi E, Afolabi AS, Onyango NO. Mathematical modeling and analysis of transmission dynamics and control of schistosomiasis. J Appl Math. (2021) 2021:6653796. doi: 10.1155/2021/6653796

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Mangal TD, Paterson S, Fenton A. Predicting the impact of long-term temperature changes on the epidemiology and control of schistosomiasis: a mechanistic model. PLoS ONE. (2008) 3:e1438. doi: 10.1371/journal.pone.0001438

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Li Y, Teng Z, Ruan S, Li M, Feng X. A mathematical model for the seasonal transmission of schistosomiasis in the lake and marshland regions of China. Math Biosci Eng. (2017) 14:1279. doi: 10.3934/mbe.2017066

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Xiang J, Chen H, Ishikawa H. A mathematical model for the transmission of Schistosoma japonicum in consideration of seasonal water level fluctuations of Poyang Lake in Jiangxi, China. Parasitol Int. (2013) 62:118–26. doi: 10.1016/j.parint.2012.10.004

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Gao SJ, Cao HH, He YY, Liu YJ, Zhang XY, Yang GJ, et al. The basic reproductive ratio of Barbour's two-host schistosomiasis model with seasonal fluctuations. Parasites Vectors. (2017) 10:1–9. doi: 10.1186/s13071-017-1983-1

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Ding C, Qiu Z, Zhu H. Multi-host transmission dynamics of schistosomiasis and its optimal control. Math Biosci Eng. (2015) 12:983. doi: 10.3934/mbe.2015.12.983

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Blower SM, Dowlatabadi H. Sensitivity and uncertainty analysis of complex models of disease transmission: an HIV model, as an example. Int Stat Rev. (1994) 62:229–43. doi: 10.2307/1403510

CrossRef Full Text | Google Scholar

40. Herdicho FF, Chukwu W, Tasman H. An optimal control of malaria transmission model with mosquito seasonal factor. Results Phys. (2021) 25:104238. doi: 10.1016/j.rinp.2021.104238

CrossRef Full Text | Google Scholar

41. Chukwu C, Nyabadza F. A theoretical model of listeriosis driven by cross contamination of ready-to-eat food products. Int J Math Math Sci. (2020) 2020:9207403. doi: 10.1155/2020/9207403

CrossRef Full Text | Google Scholar

42. Chazuka Z, Moremedi GM. In-host dynamics of the human papillomavirus (HPV) in the presence of immune response. In:Mondaini RP, , editor. Trends in Biomathematics: Chaos and Control in Epidemics, Ecosystems, and Cells: Selected Works from the 20th BIOMAT Consortium Lectures, Rio de Janeiro, Brazil, 2020. Cham: Springer (2021). p. 79–97. doi: 10.1007/978-3-030-73241-7_6

CrossRef Full Text | Google Scholar

43. Chukwu CW, Mushanyu J, Juga ML, Fatmawathi. A mathematical model for co-dynamics of listeriosis and bacterial meningitis diseases. Commun Math Biol Neurosci. (2020) 2020:ID83.

Google Scholar

44. Madubueze CE, Chazuka Z. An optimal control model for the transmission dynamics of Lassa fever. Res Squ. [Preprint]. doi: 10.21203/rs.3.rs-1513399/v1

CrossRef Full Text | Google Scholar

45. Chukwu C, Juga M, Chazuka Z, Mushanyu J. Mathematical analysis and sensitivity assessment of HIV/AIDS-listeriosis co-infection dynamics. Int J Appl Comput Math. (2022) 8:1–21. doi: 10.1007/s40819-022-01458-3

PubMed Abstract | CrossRef Full Text | Google Scholar