# Modelling the Role of Treatment, Public Health Education, and Chemical Control Strategies on Transmission Dynamics of Schistosomiasis

Authors: Timothy Kamara; Michael Byamukama; Martin Karuhanga
Journal: Journal of Mathematics
 (2022)

---

## Abstract

In this study, a mathematical model for the transmission dynamics and control of schistosomiasis is studied using a system of nonlinear ordinary differential equations. The basic reproduction numberℛ0 for the model is obtained, and its dependence on model parameters is discussed. Analytical results reveal that the disease-free equilibrium point is globally stable if and only if the basic reproduction number ℛ0<1, indicating that the disease would be wiped out of the community, and the endemic equilibrium point is globally stable if ℛ0>1 and the disease would persist at the endemic steady state. Numerical simulations reveal that a combination of treatment, public health education, and chemical control intervention strategies significantly increased the number of susceptible human population and susceptible snail population while significantly decreasing the number of infected humans, miracidia, infected snails, and cercariae. The results further indicate that a combination of treatment, public health education, and chemical control intervention strategies can effectively manage the transmission of schistosomiasis in endemic areas.


---

## Body

## 1. Introduction

Schistosomiasis, commonly known as bilharzia or snail fever, is a disease caused by parasitic flukes of genusSchistosoma. It is listed among the neglected tropical diseases (NTDs) and is endemic in many countries, particularly in tropical countries of Africa, South America, and Asia. Schistosomiasis is the second most devastating tropical disease after malaria [1], and its prevalence has continued to prevail amidst efforts to curb it.The human is the definitive host of the parasite and the snail is the intermediate host. The snail lives in fresh, shallow, and slow running water sources. The three main schistosomes that infect humans areSchistosoma japonicum, Schistosoma mansoni which cause intestinal and hepatosplenic disease, and Schistosoma haematobium that causes urogenital pathology. Each species of the schistosome is transmitted by a snail of different species; however, the species share the same pattern in the reproductive cycle [1–3].The infection cycle begins when an infected host eliminates eggs of the parasite through urine (Schistosoma haematobium) or feaces (Schistosoma mansoni and Schistosoma japonicum) into fresh water bodies [2]. The excreted eggs can stay up to seven days and hatch into miracidia that infect fresh water snails of appropriate species [4]. The miricidia then multiply through asexual reproduction within the body of the snail [1], and within 4–6 weeks after infection, the snails release cercariae into water which can then penetrate the skin of a new definitive host (human) and take about 4–6 weeks to mature in the blood vessels of the urinary or the digestive system. The cercariae continue to develop and mature to adult size and pair in sexual reproduction [5]. The period between cercariae penetration of the human skin and release of schistosome eggs in human waste takes about five weeks [6] and subsequent egg laying continues the life cycle of the parasite. Drinking water or eating food that has been washed with untreated water can also lead to infection.The disease is characterized by an itchy skin within the first few hours of penetration, fever, cough, headache, difficulty in breathing, loss of appetite and weight, diarrhea and blood in stool. In advanced cases, children develop thin legs, seizures and spinal cord inflammation which lead to physical and mental development problems [5, 7]. However, most of these conditions can be reversed and prevented through early and regular treatment with praziquantel [8]. Globally, 76.9 million people received treatment for schistosomiasis in 2020, and of these, 59.9 million were school aged children (55.2 million in African region) and 17 million adults (14.5 million in African region) [9].Studies to determine the effectiveness of the control strategies on schistosomiasis have been undertaken before. Gao et al. [10] studied the effect of treatment and snail control on schistosomiasis transmission, Das et al. [11] studied the optimal control of snail population and sanitation and treatment on schistosomiasis transmission, Nur et al. [12] determined the effect of health education and snail control on schistosomiasis transmission, and Kanyi et al. [13] determined the effect of treatment and sanitation on schistosomiasis transmission. Although there have been studies that address the role of reducing the intermediate host population (snail population) in the spread and control dynamics of schistosomiasis, the combined effect of treatment, public health education, and chemical control interventions on the disease transmission dynamics has not been explicitly established. Therefore, there is need to examine the effectiveness of control measures that include environmental modifications that separate humans from contaminated water sources through public health campaigns and snail reduction with chemical control measures to limit schistosomiasis transmission. In this study, a deterministic mathematical model describing the effect of treatment, public health education, and snail chemical control strategies on schistosomiasis transmission is formulated. The model considers control of schistosomiasis through human treatment using praziquantel and chemical control intervention by use of molluscicide (niclosamide) to reduce snail population and public health education which is intended to create change in human behavior, so as to maintain improved sanitation and hence reduce schistosome worm egg disposal in water by infected individuals.The rest of the study is structured as follows. In the next section, we present the schistosomiasis model. In Section3, we find and determine sufficient conditions for the stability of the disease-free and endemic equilibria and analyze the reproduction number of the disease, while Section 4 provides numerical results, and Section 5 concludes the study.


## 2. Model Formulation

### 2.1. Variables of the Model

The model consists of the human population divided into two epidemiological classes: the susceptible human populationSHt and the infected human population IHt such that the total human population is given by NHt=SHt+IHt. The snail host population is categorized into two epidemiological classes, namely, the susceptible snail population SSt and infected snail population ISt such that the total snail population is given by NSt=SSt+ISt. The model also includes the population of the free-living larva stages, the miracidia Mt and cercariae Ct.


### 2.2. Assumptions

The following assumptions are considered in the formulation of the model:(i)
There is no vertical transmission of the disease in humans(ii)
There is no immigration of infectious humans(iii)
Seasonal and weather variations do not affect snail populations(iv)
Susceptible host is infected only by contact with infected water(v)
The public health education awareness campaign covers the entire population(vi)
Treatment provides cure but does not prevent reinfection(vii)
The interaction between the miracidia and the susceptible snails, and the interaction between the cercariae and the susceptible humans has negligible effect on the pathogen population


### 2.3. Parameters of the Model

The parameters used in the model are described as follows:ΛH: Recruitment rate into the human populationΛS: Recruitment rate into the snail populationη: Recovery rate due to treatmentλH: Force of infection of the susceptible human populationλS: Force of infection of the susceptible snail populationϕ: Public health education efficacy parameterρ: Number of schistosome eggs excreted by infected humanγ1: Rate at which excreted eggs hatch into miracidiaγ2: Rate at which infected snails release cercariaeβ1: Rate at which susceptible humans get infectedβ2: Rate at which susceptible snails get infectedϵ: Limitation rate of growth velocityC0: Saturation constant for cercariaeℳ0: Saturation constant for miracidiaδH,δS: Disease induced mortality rates of the human and snail population, respectivelyθS,θM,θC: Elimination rates of snail, miracidia, and cercariae by molluscicide, respectivelyμH,μS: Natural death rates of the human and snail population, respectivelyμM,μC: Natural death rates of the miracidia and cercariae population, respectivelyThe parametersλH and λS are given by λH=β1C/C0+ϵC and λS=β2M/ℳ0+ϵM, respectively.


### 2.4. Equations of the Model

From the compartmental diagram in Figure1 and applying the assumptions, definitions of variables and parameters above, the following system of differential equations that describe the effect of integrated control strategies on schistosomiasis transmission is obtained:(1)dSHdt=ΛH+ηIH−1−ϕλHSH−μHSH,dIHdt=1−ϕλHSH−η+μH+δHIH,dMdt=1−ϕργ1IH−μM+θMM,dSSdt=ΛS−λSSS−μS+θSSS,dISdt=λSSS−μS+δS+θSIS,dCdt=γ2IS−μC+θCC.Figure 1
Compartmental diagram for schistosomiasis transmission.Letd1=η+μH+δH, d2=μM+θM, d3=μS+θS, d4=μS+δS+θS, and d5=μC+θC, then system (1) can be written as(2)dSHdt=ΛH+ηIH−1−ϕλHSH−μHSH,dIHdt=1−ϕλHSH−d1IH,dMdt=1−ϕργ1IH−d2M,dSSdt=ΛS−λSSS−d3SS,dISdt=λSSS−d4IS,dCdt=γ2IS−d5C.


## 2.1. Variables of the Model

The model consists of the human population divided into two epidemiological classes: the susceptible human populationSHt and the infected human population IHt such that the total human population is given by NHt=SHt+IHt. The snail host population is categorized into two epidemiological classes, namely, the susceptible snail population SSt and infected snail population ISt such that the total snail population is given by NSt=SSt+ISt. The model also includes the population of the free-living larva stages, the miracidia Mt and cercariae Ct.


## 2.2. Assumptions

The following assumptions are considered in the formulation of the model:(i)
There is no vertical transmission of the disease in humans(ii)
There is no immigration of infectious humans(iii)
Seasonal and weather variations do not affect snail populations(iv)
Susceptible host is infected only by contact with infected water(v)
The public health education awareness campaign covers the entire population(vi)
Treatment provides cure but does not prevent reinfection(vii)
The interaction between the miracidia and the susceptible snails, and the interaction between the cercariae and the susceptible humans has negligible effect on the pathogen population


## 2.3. Parameters of the Model

The parameters used in the model are described as follows:ΛH: Recruitment rate into the human populationΛS: Recruitment rate into the snail populationη: Recovery rate due to treatmentλH: Force of infection of the susceptible human populationλS: Force of infection of the susceptible snail populationϕ: Public health education efficacy parameterρ: Number of schistosome eggs excreted by infected humanγ1: Rate at which excreted eggs hatch into miracidiaγ2: Rate at which infected snails release cercariaeβ1: Rate at which susceptible humans get infectedβ2: Rate at which susceptible snails get infectedϵ: Limitation rate of growth velocityC0: Saturation constant for cercariaeℳ0: Saturation constant for miracidiaδH,δS: Disease induced mortality rates of the human and snail population, respectivelyθS,θM,θC: Elimination rates of snail, miracidia, and cercariae by molluscicide, respectivelyμH,μS: Natural death rates of the human and snail population, respectivelyμM,μC: Natural death rates of the miracidia and cercariae population, respectivelyThe parametersλH and λS are given by λH=β1C/C0+ϵC and λS=β2M/ℳ0+ϵM, respectively.


## 2.4. Equations of the Model

From the compartmental diagram in Figure1 and applying the assumptions, definitions of variables and parameters above, the following system of differential equations that describe the effect of integrated control strategies on schistosomiasis transmission is obtained:(1)dSHdt=ΛH+ηIH−1−ϕλHSH−μHSH,dIHdt=1−ϕλHSH−η+μH+δHIH,dMdt=1−ϕργ1IH−μM+θMM,dSSdt=ΛS−λSSS−μS+θSSS,dISdt=λSSS−μS+δS+θSIS,dCdt=γ2IS−μC+θCC.Figure 1
Compartmental diagram for schistosomiasis transmission.Letd1=η+μH+δH, d2=μM+θM, d3=μS+θS, d4=μS+δS+θS, and d5=μC+θC, then system (1) can be written as(2)dSHdt=ΛH+ηIH−1−ϕλHSH−μHSH,dIHdt=1−ϕλHSH−d1IH,dMdt=1−ϕργ1IH−d2M,dSSdt=ΛS−λSSS−d3SS,dISdt=λSSS−d4IS,dCdt=γ2IS−d5C.


## 3. Model Analysis

### 3.1. Positivity and Boundedness of the Solutions

Theorem 1.
Let the initial population be(3)SH0,IH0,M0,SS0,IS0,C0≥0∈Ω.Then, the solution setSHt,IHt,Mt,SSt,ISt,Ct of (2) is positive for all t≥0.Proof.
From the first equation of system (2), we get(4)ΛH+ηIH−1−ϕλH+μHSH≥−1−ϕλH+μHSH.
This implies that(5)dSHdt≥−1−ϕλH+μHSH.
By separating variables, (5) is integrated to give(6)lnSH≥−1−ϕλH+μHt+K,where K is the constant of integration. Applying the initial condition SH0=SH0 to (6) gives K=lnSH0.
Hence,(7)SHt≥SH0e−1−ϕλH+μHt.
Therefore,(8)SHt≥0.
Again from system (2), one gets(9)dIHdt≥−d1IH,dMdt≥−d2M,dSSdt≥−λSSS−d3SS,dISdt≥−d4IS,dCdt≥−d5C.
Now, separating the variables and integrating with the initial conditionsIH0=IH0,M0=M0,SS0=SS0,IS0=IS0,C0=C0 gives(10)IHt≥IH0e−d1t≥0,Mt≥M0e−d2t≥0,SSt≥SS0e−λS+d3t≥0,ISt≥IS0e−d4t≥0,Ct≥C0e−d5t≥0.
Thus, the solution set of model (2) is positive for all values of t≥0.Theorem 2.
The feasible solution set of the model system (2) is mathematically and epidemiologically well-posed in the domain as(11)Ω=SH,IH,M,SS,IS,C∈ℝ6+:NH≤ΛHμH,M≤1−ϕργ1ΛHμHd2,NS≤ΛSd4,C≤γ2ΛSd3d5.Proof.
FromNHt=SHt+IHt, it follows that(12)dNHdt=ΛH−μHN−δHIH≤ΛH−μHN,and integrating (12) with the initial condition NH0=NH0 gives(13)NHt≤ΛHμH−1μHΛH−μHNH0e−μHt,and as t⟶∞ in (13), we get(14)NHt≤ΛHμH.
For the snail host population, we get(15)NSt=SSt+ISt,dNSdt=ΛS−d3NS−δSIS≤ΛS−d3NS.
Integrating (15) with initial condition NS0=NS0 gives(16)NSt≤ΛSd3−1d3ΛS−d3NS0e−d3t.
Ast⟶∞ in (16),(17)NSt≤ΛSd3.
For the miracidia free-living larva stage subpopulation,(18)dMdt=1−ϕργ1IH−d2M≤1−ϕργ1ΛHμH−d2M,and integrating (18) with the initial condition M0=M0 gives(19)Mt≤1−ϕργ1ΛHμHd2−1d21−ϕργ1ΛHμH−d2M0e−d2t.
Ast⟶∞ in (19), we get(20)Mt≤1−ϕργ1ΛHμHd2.
For the cercariae free-living larva stage subpopulation,(21)dCdt=γ2IS−d5C≤γ2ΛSd3−d5C.
Integrating (21) with initial condition C0=C0 gives(22)Ct≤γ2ΛSd3d5−1d5γ2ΛSd3−d5C0e−d5t.
Ast⟶∞ in (22), we have(23)Ct≤γ2ΛSd3d5.
Thus, the regionΩ is positively invariant, and hence, the model is well posed and biologically meaningful.


### 3.2. The Disease-Free Equilibrium Point

The disease-free equilibrium point is obtained by setting the right-hand side of (2) to zero and is defined at the steady state where there is no disease. Therefore, in absence of the disease, we have IH=M=IS=C=0. Hence, the disease-free equilibrium point is given as(24)E0=SH,0,0,SS,0,0=ΛHμH,0,0,ΛSd3,0,0.


### 3.3. Basic Reproduction Number

The basic reproduction numberℛ0 is the average number of secondary infections caused by an infectious individual in a fully susceptible population. It is a threshold for predicting the disease as well as evaluating the effectiveness of the intervention strategies.The value ofℛ0 is obtained using the next generation operator approach described in [14, 15]. Let Fi be the matrix of the rates of the appearance of new infections in compartment i and Vi be matrix of the rates of the transfer of individuals into and out of compartment i. Thus, from (2), we have(25)Fi=1−ϕλHSH0λSSS0,Vi=d1IHd2M−1−ϕργ1λHd4ISd5C−γ2IS.LetF and V be the Jacobian matrices of Fi and Vi, respectively, evaluated at the disease-free equilibrium point (24). Then, one gets(26)F=0001−ϕβ1ΛHC0μH00000β2ΛSℳ0d3000000,V=d1000−1−ϕργ1d20000d4000−γ2d5.EvaluatingFV−1, we get(27)FV−1=00β1γ21−ϕΛHC0d4d5μH1−ϕβ1ΛHC0d5μH00001−ϕργ1β2ΛSℳ0d1d2d3β2ΛSℳ0d2d3000000.The eigenvalues of matrixFV−1 are obtained from(28)detFV−1−λI=0,where I is a 4 by 4 identity matrix. This gives(29)λ2λ2−1−ϕβ1ργ1ΛHC0d1d2μH1−ϕβ2ργ2ΛSℳ0d3d4d5=0,and the roots of the characteristic (29) are λ=0,β1β2ργ1γ2ΛHΛS1−ϕ2/C0ℳ0μHd1d2d3d4d5 and λ=−β1β2ργ1γ2ΛHΛS1−ϕ2/C0ℳ0μHd1d2d3d4d5.The basic reproduction number is the dominant eigenvalue of the matrixFV−1 and is thus given as(30)ℛ0=β1β2ργ1γ2ΛHΛS1−ϕ2C0ℳ0μHd1d2d3d4d5=1−ϕβ1ργ1ΛHC0d1d2μH1−ϕβ2ργ2ΛSℳ0d3d4d5=ℛ0Hℛ0S,where(31)ℛ0H=1−ϕβ1ργ1ΛHC0d1d2μH,ℛ0S=1−ϕβ2ργ2ΛSℳ0d3d4d5.Thus, the basic reproduction number,ℛ0, can be written in terms of ℛ0H and R0S, where ℛ0H describes the average number of secondary infections by an infected snail in a completely susceptible population (the average number of susceptible humans infected by the infectious snail) and ℛ0S describes the average number of secondary infections by an infected human in a completely susceptible population (the average number of susceptible snails infected by the infectious human).


### 3.4. The Endemic Equilibrium Point

The endemic equilibrium point determines the persistence of the disease in the community and this is evaluated whenℛ0>1.Theorem 3.
System (2) has a unique endemic equilibrium point when R0>1.Proof.
Setting the right-hand side of system (2) to zero, we have(32)ΛH+ηIH−1−ϕλHSH−μHSH=0,1−ϕλHSH−d1IH=0,1−ϕργ1IH−d2M=0,ΛS−λSSS−d3SS=0,λSSS−d4IS=0,γ2IS−d5C=0.
Adding the first two equations of (15) and rearranging gives(33)SH=ΛH−μH+δHIHμH.
From the third equation of (15), M=1−ϕργ1/d2IH. Substituting this into the fourth equation gives(34)SS=ΛSℳ0d2+ϵ1−ϕργ1IHℳ0d2d3+β2+ϵd31−ϕργ1IH.
Putting (33) and (34) into the fifth equation of (32) gives(35)IS=ΛSβ21−ϕργ1IHℳ0d2d3d4+d4β2+ϵd31−ϕργ1IH.
From the sixth equation of (32), we have(36)C=γ2d6IS.
Consequently, (35) and (36) yield(37)C=ΛSγ2β21−ϕργ1IHd5ℳ0d2d3d4+d4β2+ϵd31−ϕργ1IH.
Now, using (34) and (36) in the second equation of (32) and then simplifying gives(38)IHAIH−BℛO2−1=0,where(39)A=ργ11−ϕC0μHd1d4d5ϵd3+β2+ϵμHΛSγ2β2+ΛSβ1β2γ21−ϕμH+δH,B=C0ℳ0d1d2d3d4d5.
From (38), IH=0 corresponds to the disease-free equilibrium and IH=B/AℛO2−1forA≠0 indicates endemicity of the disease. Thus, system (2) has a unique endemic equilibrium point E1SH∗,IH∗,M∗,SS∗,IS∗,C∗ when ℛ0>1, where(40)SH∗=ΛH−μH+δHIH∗μH,IH∗=C0ℳ0d1d2d3d4d5ργ11−ϕC0μHd1d4d5ϵd3+β2+ϵμHΛSγ2β2+ΛSβ1β2γ21−ϕμH+δHRO2−1,M∗=1−ϕργ1d2IH∗,SS∗=ΛSℳ0d2+ϵ1−ϕργ1IH∗ℳ0d2d3+β2+ϵd31−ϕργ1IH∗,IS∗=ΛSβ21−ϕργ1IH∗ℳ0d2d3d4+d4β2+ϵd31−ϕργ1IH∗,C∗=ΛSγ2β21−ϕργ1IH∗d5ℳ0d2d3d4+d4β2+ϵd31−ϕργ1IH∗.


### 3.5. Local and Global Stability of the Disease-Free Equilibrium Point

Theorem 4.
The disease-free equilibrium pointE0 of system (2) is locally asymptotically stable whenever R0<1 and unstable whenever R0>1.Proof.
The Jacobian matrix of system (2) evaluated at the disease-free equilibrium point E0 is given by(41)JE0=−μH0000−1−ϕβ1ΛHC0μH0−d10001−ϕβ1ΛHC0μH01−ϕργ1−d200000−β2ΛSℳ0d3−d30000β2ΛSℳ0d30−d400000γ2−d5.
The disease-free equilibrium pointE0 is locally asymptotically stable if the real parts of the eigenvalues of JE0 are negative; otherwise, it is unstable.
The eigenvalues ofJE0 are obtained by setting the characteristic polynomial.
P6λ=detJE0−λI equals to zero. Thus, we have(42)P6λ=−μH−λ−d3−λ−d1−λ00−1−ϕβ1ΛHC0μH1−ϕργ1−d2−λ000−β2ΛSℳ0d3−d4−λ000γ2−d5−λ=0.
From (42), it is clear that λ1=−μH, λ2=−d3. The other remaining eigenvalues are obtained by considering the matrix(43)−d1−λ00−1−ϕβ1ΛHC0μH1−ϕργ1−d2−λ000−β2ΛSℳ0d3−d4−λ000γ2−d5−λ,whose characteristic equation is(44)λ4+a1λ3+a2λ2+a3λ+a4=0,where a1=d1+d2+d4+d5,a2=d1d2+d1d4+d1d5+d2d4+d2d5+d4d5,a3=d1d2d4+d1d2d5+d1d4d5+d2d4d5,a4=d1d2d4d51−ℛ02.
According to [16], the Routh–Hurwitz conditions for the stability of the polynomial (44) are as follows:(i)
a1>0(ii)
a1a2−a3>0(iii)
a3a1a2−a3−a12a4>0(iv)
a4>0
Clearly (i) is satisfied and (iv) holds true wheneverℛ0<1. It remains to check conditions (ii) and (iii):
and(45)a1a2−a3=d12d2+d12d4+d12d5+d1d22+2d1d2d4+2d1d2d5+d22d4+d22d5+d1d42+2d1d4d5+d2d42+2d2d4d5+d42d5+d2d52+d4d52>0,and(46)a3a1a2−a3−a12a4=d1d2d4+d1d2d5+d1d4d5+d2d4d5d1+d2+d4+d5d1d2+d1d4+d1d5+d2d4+d2d5+d4d5−d1d2d4+d1d2d5+d1d4d5+d2d4d5−d1+d2+d4+d52d1d2d4d51−R02>0,ifℛ0<1.
Thus, the Routh–Hurwitz conditions hold when all the eigenvalues have negative real parts makingE0 locally asymptotically stable whenever ℛ0<1.
To study the global asymptotic stability of the disease-free equilibrium, the method described by Castillo-Chavez et al. [17] and in [18] has been employed. The model system (2) is expressed in the form:(47)dX1dt=FX1,X2,dX2dt=GX1,X2,GX1,0=0,where X1=SH,SS comprises of the uninfected subpopulation and X2=IH,M,IS,C comprises of the infected subpopulation. Let U0=X1∗,0 be the disease-free equilibrium point of system (47). Then, U0 is globally asymptotically stable provided that ℛ0<1 and that the following conditions hold true:H1: For dX1/dt=FX1,0,X1∗=SH0,SS0 is globally asymptotically stableH2: GX1,X2=AX2−G^X1,X2,G^X1,X2≥0 for X1,X2∈Ω, where A=DX2GX1∗,0 is an M-matrix (the off-diagonal elements are non-negative) and Ω is the region where the model makes biological senseTheorem 5.
The disease-free equilibrium pointE0 for system (2) is globally asymptotically stable provided ℛ0<1 and conditions H1 and H2 hold true.Proof.
For system (2), we have(48)dX1dt=FX1,X2=ΛH+ηIH−1−ϕλHSH−μHSHΛS−λSSS−d3SS,dX2dt=GX1,X2=1−ϕλHSH−d1IH1−ϕργ1IH−d2MλSSS−d4ISγ2IS−d5C,FX1,0=ΛH−μHSHΛS−d3SS.
Ast⟶∞, we observe that X1,0⟶E0∈Ω, that is, SH⟶ΛH/μH and SS⟶ΛS/d3. Thus, H1 holds true:(49)A=∂GX1,0∂X2|U0=−d1001−ϕβ1SH0C01−ϕργ2−d2000β2SH0ℳ0−d4000γ2−d5,and(50)G^X1,X2=AX2−GX1,X2=β11−ϕCSH0C0−SHC0+ϵC0β2MSS0ℳ0−SSℳ0+ϵM0.
From (49), clearly A is an M-matrix (all the off-diagonal elements are non-negative), and from (50), G^X1,X2≥0 since SH0/C0≥SH/C0+ϵC and SS0/ℳ0≥SS/ℳ0+ϵM. Thus, H2 holds true, and the proof is complete.


### 3.6. Local Stability of the Endemic Equilibrium

The center manifold theory by Castillo-Chavez [19] is used to establish the local asymptotic stability of the endemic equilibrium point.Theorem 6.
Consider the following general system of ordinary differential equations with a parameterϑ:(51)dxdt=fx,ϑ,f:ℝn×ℝ⟶ℝn,f∈ℂ2ℝn×ℝ,and without loss of generality, it is assumed that 0 is an equilibrium point of the system (51) for all values of the parameter ϑ, that is,(52)f0,ϑ≡0for allϑ.Assume the following points:A1. A=Dxf0,0=∂fi/∂xj0,0 is the linearization matrix of system (51) around the equilibrium 0 with ϑ evaluated at 0. Zero is a simple eigenvalue of A, and other eigenvalues have negative real parts.A2. The matrix A has a non-negative right eigenvector w and a left eigenvector v each corresponding to the zero eigenvalue.Letfk be the kth component of f, and then, we get(53)a=∑k,i,j=1nvkwiwj∂2fk∂xi∂xj0,0,b=∑k,i=1nvkwi∂2fk∂xi∂ϑ0,0.Then, the local dynamics of the system around 0 are completely determined by the signs ofa and b that can be given as follows:(i)
a>0,b>0. When ϑ<0 with ϑ≪1, 0 is locally asymptotically stable and there exists a positive unstable equilibrium; when 0<ϑ≪1, 0 is unstable, and there exists a negative and locally asymptotically stable equilibrium.(ii)
a<0,b<0. When ϑ<0 with ϑ≪1, 0 is unstable; when 0<ϑ≪1, 0 is locally asymptotically stable, and there exists a positive unstable equilibrium.(iii)
a>0,b<0. When ϑ<0 with ϑ≪1, 0 is unstable, and there exists a locally asymptotically stable negative equilibrium; when 0<ϑ≪1, 0 is stable, and a positive unstable equilibrium appears.(iv)
a<0,b>0. When ϑ changes from negative to positive, 0 changes its stability from stable to unstable. Correspondingly, a negative unstable equilibrium becomes positive and locally asymptotically stable.Theorem 7.
The unique endemic equilibrium point,E1 of system (2) is locally asymptotically stable if R0>1 but close to 1.Proof.
We apply the change of variable on system (2) as follows. Let SH=x1,IH=x2,M=x3,SS=x4,IS=x5,C=x6. The model system (2) can be written as dX/dt=F, where X=x1,x2,x3,x4,x5,x6T and F=f1,f2,f3,f4,f5,f6:(54)dx1dt=f1=ΛH+ηx2−1−ϕβ1x6C0+ϵx6x1−μHx1,dx2dt=f2=1−ϕβ1x6C0+ϵx6x1−d1x2,dx3dt=f3=1−ϕργ1x2−d2x3,dx4dt=f4=ΛS−β2x3ℳ0+ϵx3x4−d3x4,dx5dt=f5=β2x3ℳ0+ϵx3x4−d4x5,dx6dt=f6=γ2x5−d5x6.
Takingβ1 as the bifurcation parameter while considering the case ℛ0=1 gives(55)β1=β1∗=C0ℳ0μHd1d2d3d4d5β2ργ1γ2ΛHΛS1−ϕ2.
The linearization around the disease-free equilibrium point evaluated atβ1∗ gives(56)A=−μH0000−1−ϕβ1∗ΛHC0μH0−d10001−ϕβ1∗ΛHC0μH01−ϕργ1−d200000−β2ΛSℳ0d3−d30000β2ΛSℳ0d30−d400000γ2−d5.
MatrixA has zero as an eigenvalue, and this can be deduced from (21)–(23) with ℛ0=1, which yields(57)λλ3+a1λ2+a2λ+a3=0.
The right eigenvectorw=w1,w2,w3,w4,w5,w6T corresponding to the zero eigenvalue is computed from Aw=0, which yields(58)w1=−β1∗1−ϕΛHμH2C0w6,w2=β1∗1−ϕΛHC0μHd1w6,w3=ργ1β1∗ΛH1−ϕ2C0μHd1d2w6,w4=−ργ1β1∗β2ΛSΛH1−ϕ2ℳ0C0μHd1d2d32w6,w5=d5γ2w6,w6>0.
The left eigenvectorv=v1,v2,v3,v4,v5,v6T is computed from vA=0 which yields(59)v1=v4=0,v2=C0μHd5β1∗1−ϕΛHv6,v3=C0μHd1d5ργ1β1∗ΛH1−ϕ2v6,v5=γ2d4v6,v6>0.
Lettingw6=1, it can be seen from (58) that any other right eigenvector of matrix A corresponding to the zero eigenvalue is a constant multiple of w (the same is true for the left eigenvector in (59) by letting v6=1). Hence, zero is a simple eigenvalue.
In the computation ofa and b, only nonzero partial derivatives of system (54) are used. Since v1=v4=0, the nonzero partial derivatives computed at the disease-free equilibrium are from f2,f3,f5, and f6. These are obtained as(60)∂2f2∂x1∂x6=β1∗C01−ϕ,∂2f2∂2x6=−2ϵβ1∗C021−ϕx1∗,∂2f5∂x3∂x4=β2ℳ0,∂2f5/∂2x3=−2ϵβ2x5∗/M02 and ∂2f2/∂x6∂β1∗=1−ϕ/C0x1∗. Thus, it follows that(61)a=v2w1w6∂2f2∂x1∂x6+v2w62∂2f2∂2x6+v5w3w4∂2f5∂x3∂x4+v5w32∂2f5∂2x3,a=−v6w62β1∗1−ϕd5C0μH+2ϵd5C0+d4d52ΛSγ2+2ϵd3d4d52ΛSβ2γ22<0,b=d6β1∗v7w7>0.
Sincea<0 and b>0, then by condition iv of Theorem 5, system (54) exhibits a forward bifurcation at ℛ0=1, implying that the unique endemic equilibrium is locally asymptotically stable when ℛ0>1 and close to 1.


### 3.7. Global Stability of the Endemic Equilibrium Point

Theorem 8.
Ifℛ0≥1, the endemic equilibrium E1 of system (2) is globally asymptotically stable.Proof.
To analyze the global stability of endemic equilibrium, we consider the candidate Lyapunov functionLSH,IH,M,SS,IS,C given by(62)L=SH−SH∗+SH∗lnSH∗SH+IH−IH∗+IH∗lnIH∗IH+M−M∗+M∗lnM∗M+SS−SS∗+SS∗lnSS∗SS+IS−IS∗+IS∗lnIS∗IS+C−C∗+C∗lnC∗C.
Clearly,LSH∗,IH∗,M∗.SS∗,IS∗,C∗=0, and L>0, otherwise. Moreover, L is radially unbounded in R+6. So, it remains to determine the sign of L′. Differentiating (62) with respect to t gives(63)dLdt=SH−SH∗SHdSHdt+IH−IH∗IHdIHdt+M−M∗MdMdt+SS−SS∗SSdSSdt+IS−IS∗ISdISdt+C−C∗CdCdt.
Substituting for values ofdSH/dt, dIH/dt, dM/dt, dSS/dt, dIS/dt, and dC/dt in (63) yields(64)dLdt=SH−SH∗SHΛH+ηIH−1−ϕλHSH−μHSH+IH−IH∗IH1−ϕλHSH−d1IH+M−M∗M1−ϕργ1IH−d2M+SS−SS∗SSΛS−λSSS−d3SS+IS−IS∗ISλSSS−d4IS+C−C∗Cγ2IS−d5C.
Rearranging positive and negative terms in (64) gives(65)dLdt=P−Q,where(66)P=ΛH+ηIH+1−ϕλH+μH+λS+d3SH∗+1−ϕλHSH+λSSS+d1IH∗+η+1−ϕργ1IH+d2M∗+ΛS+λSSS+d4IS∗+γ2IS+d5C∗,Q=1−ϕλH+μHSH+ΛH+ηIHSH∗SH+d1IH+1−ϕλHSHIH∗IH+d2M+1−ϕργ1IHM∗M+λS+d3SS+ΛSSS∗SS+d4IS+λSSSIS∗IS+d5C+γ2ISC∗C.
IfP≤Q, then dL/dt≤0 noting that dL/dt=0 if and only if SH=SH∗, IH=IH∗, M=M∗, SS=SS∗, IS=IS∗, and C=C∗. Therefore, the largest compact invariant set in SH,IH,M,SS,IS,C∈Ω:dL/dt=0 is the singleton E1, where E1 is the endemic equilibrium point of system (2). Thus, by Lasalle’s invariance principle [20], it implies that E1 is globally asymptotically stable in Ω if P≤Q; which holds if and only if ℛ0≥1.


### 3.8. Sensitivity Analysis

In this subsection, sensitivity analysis of the reproduction numberℛ0 on different parameters is carried out in order to determine their relative importance to the disease transmission in the community. The approach by Chitnis et al. [21] has been employed to obtain the relative change in a state variable when the parameter changes. An expression for the sensitivity of ℛ0 with respect to p is given by(67)γpℛ0=∂ℛ0∂ppℛ0,where p is any sensitive parameter to be considered in relation to ℛ0.From (67), expressions for the sensitivity indices are as follows:(68)γβ1ℛ0=γβ2ℛ0=γρℛ0=γβ1ℛ0=γγ1ℛ0=γγ2ℛ0=γλHℛ0=γλSℛ0=12,γμCℛ0=−12μCd5,γθMℛ0=−12θMd2,γμMℛ0=−12μMd1,γδHℛ0=−12δHd1=γδSℛ0=−12δSd4,γθCℛ0=−12θCd5,γθSℛ0=−12θSd3+d4d3d4,γμSℛ0=−12μSd3+d4d3d4.The parameter values used areμC=0.004,μS=0.000569,δS=0.0004012,μM=0.008,θS=0.1,η=0.0181, and μH=0.0000384.Table1 shows the sensitivity indices of ℛ0 with respect to different parameters from the most sensitive to the least sensitive. It is noted that the value of ℛ0 increases when parameter values β1, β2, ρ, γ1, γ2, ΛH, and ΛS increase, while the other parameters values are kept constant. Since these parameters have positive indices, the endemicity of the disease is increased. When parameter μH, θS, θC, μC, δH, μS, θM, and δS are increased while the other parameters values are kept constant, the value of ℛ0 decreases implying that the endemicity of the disease is decreased.Table 1
Numerical values of sensitivity indices ofℛ0.
Parameter symbolSensitivity indexμH−0.5009β1+0.5000β2+0.5000ρ+0.5000γ1+0.5000γ2+0.5000ΛH+0.5000ΛS+0.5000μM−0.4956θS−0.4952θC−0.3333μC−0.1667δH−0.0885μS−0.0056θM−0.0044δS−0.0020


## 3.1. Positivity and Boundedness of the Solutions

Theorem 1.
Let the initial population be(3)SH0,IH0,M0,SS0,IS0,C0≥0∈Ω.Then, the solution setSHt,IHt,Mt,SSt,ISt,Ct of (2) is positive for all t≥0.Proof.
From the first equation of system (2), we get(4)ΛH+ηIH−1−ϕλH+μHSH≥−1−ϕλH+μHSH.
This implies that(5)dSHdt≥−1−ϕλH+μHSH.
By separating variables, (5) is integrated to give(6)lnSH≥−1−ϕλH+μHt+K,where K is the constant of integration. Applying the initial condition SH0=SH0 to (6) gives K=lnSH0.
Hence,(7)SHt≥SH0e−1−ϕλH+μHt.
Therefore,(8)SHt≥0.
Again from system (2), one gets(9)dIHdt≥−d1IH,dMdt≥−d2M,dSSdt≥−λSSS−d3SS,dISdt≥−d4IS,dCdt≥−d5C.
Now, separating the variables and integrating with the initial conditionsIH0=IH0,M0=M0,SS0=SS0,IS0=IS0,C0=C0 gives(10)IHt≥IH0e−d1t≥0,Mt≥M0e−d2t≥0,SSt≥SS0e−λS+d3t≥0,ISt≥IS0e−d4t≥0,Ct≥C0e−d5t≥0.
Thus, the solution set of model (2) is positive for all values of t≥0.Theorem 2.
The feasible solution set of the model system (2) is mathematically and epidemiologically well-posed in the domain as(11)Ω=SH,IH,M,SS,IS,C∈ℝ6+:NH≤ΛHμH,M≤1−ϕργ1ΛHμHd2,NS≤ΛSd4,C≤γ2ΛSd3d5.Proof.
FromNHt=SHt+IHt, it follows that(12)dNHdt=ΛH−μHN−δHIH≤ΛH−μHN,and integrating (12) with the initial condition NH0=NH0 gives(13)NHt≤ΛHμH−1μHΛH−μHNH0e−μHt,and as t⟶∞ in (13), we get(14)NHt≤ΛHμH.
For the snail host population, we get(15)NSt=SSt+ISt,dNSdt=ΛS−d3NS−δSIS≤ΛS−d3NS.
Integrating (15) with initial condition NS0=NS0 gives(16)NSt≤ΛSd3−1d3ΛS−d3NS0e−d3t.
Ast⟶∞ in (16),(17)NSt≤ΛSd3.
For the miracidia free-living larva stage subpopulation,(18)dMdt=1−ϕργ1IH−d2M≤1−ϕργ1ΛHμH−d2M,and integrating (18) with the initial condition M0=M0 gives(19)Mt≤1−ϕργ1ΛHμHd2−1d21−ϕργ1ΛHμH−d2M0e−d2t.
Ast⟶∞ in (19), we get(20)Mt≤1−ϕργ1ΛHμHd2.
For the cercariae free-living larva stage subpopulation,(21)dCdt=γ2IS−d5C≤γ2ΛSd3−d5C.
Integrating (21) with initial condition C0=C0 gives(22)Ct≤γ2ΛSd3d5−1d5γ2ΛSd3−d5C0e−d5t.
Ast⟶∞ in (22), we have(23)Ct≤γ2ΛSd3d5.
Thus, the regionΩ is positively invariant, and hence, the model is well posed and biologically meaningful.


## 3.2. The Disease-Free Equilibrium Point

The disease-free equilibrium point is obtained by setting the right-hand side of (2) to zero and is defined at the steady state where there is no disease. Therefore, in absence of the disease, we have IH=M=IS=C=0. Hence, the disease-free equilibrium point is given as(24)E0=SH,0,0,SS,0,0=ΛHμH,0,0,ΛSd3,0,0.


## 3.3. Basic Reproduction Number

The basic reproduction numberℛ0 is the average number of secondary infections caused by an infectious individual in a fully susceptible population. It is a threshold for predicting the disease as well as evaluating the effectiveness of the intervention strategies.The value ofℛ0 is obtained using the next generation operator approach described in [14, 15]. Let Fi be the matrix of the rates of the appearance of new infections in compartment i and Vi be matrix of the rates of the transfer of individuals into and out of compartment i. Thus, from (2), we have(25)Fi=1−ϕλHSH0λSSS0,Vi=d1IHd2M−1−ϕργ1λHd4ISd5C−γ2IS.LetF and V be the Jacobian matrices of Fi and Vi, respectively, evaluated at the disease-free equilibrium point (24). Then, one gets(26)F=0001−ϕβ1ΛHC0μH00000β2ΛSℳ0d3000000,V=d1000−1−ϕργ1d20000d4000−γ2d5.EvaluatingFV−1, we get(27)FV−1=00β1γ21−ϕΛHC0d4d5μH1−ϕβ1ΛHC0d5μH00001−ϕργ1β2ΛSℳ0d1d2d3β2ΛSℳ0d2d3000000.The eigenvalues of matrixFV−1 are obtained from(28)detFV−1−λI=0,where I is a 4 by 4 identity matrix. This gives(29)λ2λ2−1−ϕβ1ργ1ΛHC0d1d2μH1−ϕβ2ργ2ΛSℳ0d3d4d5=0,and the roots of the characteristic (29) are λ=0,β1β2ργ1γ2ΛHΛS1−ϕ2/C0ℳ0μHd1d2d3d4d5 and λ=−β1β2ργ1γ2ΛHΛS1−ϕ2/C0ℳ0μHd1d2d3d4d5.The basic reproduction number is the dominant eigenvalue of the matrixFV−1 and is thus given as(30)ℛ0=β1β2ργ1γ2ΛHΛS1−ϕ2C0ℳ0μHd1d2d3d4d5=1−ϕβ1ργ1ΛHC0d1d2μH1−ϕβ2ργ2ΛSℳ0d3d4d5=ℛ0Hℛ0S,where(31)ℛ0H=1−ϕβ1ργ1ΛHC0d1d2μH,ℛ0S=1−ϕβ2ργ2ΛSℳ0d3d4d5.Thus, the basic reproduction number,ℛ0, can be written in terms of ℛ0H and R0S, where ℛ0H describes the average number of secondary infections by an infected snail in a completely susceptible population (the average number of susceptible humans infected by the infectious snail) and ℛ0S describes the average number of secondary infections by an infected human in a completely susceptible population (the average number of susceptible snails infected by the infectious human).


## 3.4. The Endemic Equilibrium Point

The endemic equilibrium point determines the persistence of the disease in the community and this is evaluated whenℛ0>1.Theorem 3.
System (2) has a unique endemic equilibrium point when R0>1.Proof.
Setting the right-hand side of system (2) to zero, we have(32)ΛH+ηIH−1−ϕλHSH−μHSH=0,1−ϕλHSH−d1IH=0,1−ϕργ1IH−d2M=0,ΛS−λSSS−d3SS=0,λSSS−d4IS=0,γ2IS−d5C=0.
Adding the first two equations of (15) and rearranging gives(33)SH=ΛH−μH+δHIHμH.
From the third equation of (15), M=1−ϕργ1/d2IH. Substituting this into the fourth equation gives(34)SS=ΛSℳ0d2+ϵ1−ϕργ1IHℳ0d2d3+β2+ϵd31−ϕργ1IH.
Putting (33) and (34) into the fifth equation of (32) gives(35)IS=ΛSβ21−ϕργ1IHℳ0d2d3d4+d4β2+ϵd31−ϕργ1IH.
From the sixth equation of (32), we have(36)C=γ2d6IS.
Consequently, (35) and (36) yield(37)C=ΛSγ2β21−ϕργ1IHd5ℳ0d2d3d4+d4β2+ϵd31−ϕργ1IH.
Now, using (34) and (36) in the second equation of (32) and then simplifying gives(38)IHAIH−BℛO2−1=0,where(39)A=ργ11−ϕC0μHd1d4d5ϵd3+β2+ϵμHΛSγ2β2+ΛSβ1β2γ21−ϕμH+δH,B=C0ℳ0d1d2d3d4d5.
From (38), IH=0 corresponds to the disease-free equilibrium and IH=B/AℛO2−1forA≠0 indicates endemicity of the disease. Thus, system (2) has a unique endemic equilibrium point E1SH∗,IH∗,M∗,SS∗,IS∗,C∗ when ℛ0>1, where(40)SH∗=ΛH−μH+δHIH∗μH,IH∗=C0ℳ0d1d2d3d4d5ργ11−ϕC0μHd1d4d5ϵd3+β2+ϵμHΛSγ2β2+ΛSβ1β2γ21−ϕμH+δHRO2−1,M∗=1−ϕργ1d2IH∗,SS∗=ΛSℳ0d2+ϵ1−ϕργ1IH∗ℳ0d2d3+β2+ϵd31−ϕργ1IH∗,IS∗=ΛSβ21−ϕργ1IH∗ℳ0d2d3d4+d4β2+ϵd31−ϕργ1IH∗,C∗=ΛSγ2β21−ϕργ1IH∗d5ℳ0d2d3d4+d4β2+ϵd31−ϕργ1IH∗.


## 3.5. Local and Global Stability of the Disease-Free Equilibrium Point

Theorem 4.
The disease-free equilibrium pointE0 of system (2) is locally asymptotically stable whenever R0<1 and unstable whenever R0>1.Proof.
The Jacobian matrix of system (2) evaluated at the disease-free equilibrium point E0 is given by(41)JE0=−μH0000−1−ϕβ1ΛHC0μH0−d10001−ϕβ1ΛHC0μH01−ϕργ1−d200000−β2ΛSℳ0d3−d30000β2ΛSℳ0d30−d400000γ2−d5.
The disease-free equilibrium pointE0 is locally asymptotically stable if the real parts of the eigenvalues of JE0 are negative; otherwise, it is unstable.
The eigenvalues ofJE0 are obtained by setting the characteristic polynomial.
P6λ=detJE0−λI equals to zero. Thus, we have(42)P6λ=−μH−λ−d3−λ−d1−λ00−1−ϕβ1ΛHC0μH1−ϕργ1−d2−λ000−β2ΛSℳ0d3−d4−λ000γ2−d5−λ=0.
From (42), it is clear that λ1=−μH, λ2=−d3. The other remaining eigenvalues are obtained by considering the matrix(43)−d1−λ00−1−ϕβ1ΛHC0μH1−ϕργ1−d2−λ000−β2ΛSℳ0d3−d4−λ000γ2−d5−λ,whose characteristic equation is(44)λ4+a1λ3+a2λ2+a3λ+a4=0,where a1=d1+d2+d4+d5,a2=d1d2+d1d4+d1d5+d2d4+d2d5+d4d5,a3=d1d2d4+d1d2d5+d1d4d5+d2d4d5,a4=d1d2d4d51−ℛ02.
According to [16], the Routh–Hurwitz conditions for the stability of the polynomial (44) are as follows:(i)
a1>0(ii)
a1a2−a3>0(iii)
a3a1a2−a3−a12a4>0(iv)
a4>0
Clearly (i) is satisfied and (iv) holds true wheneverℛ0<1. It remains to check conditions (ii) and (iii):
and(45)a1a2−a3=d12d2+d12d4+d12d5+d1d22+2d1d2d4+2d1d2d5+d22d4+d22d5+d1d42+2d1d4d5+d2d42+2d2d4d5+d42d5+d2d52+d4d52>0,and(46)a3a1a2−a3−a12a4=d1d2d4+d1d2d5+d1d4d5+d2d4d5d1+d2+d4+d5d1d2+d1d4+d1d5+d2d4+d2d5+d4d5−d1d2d4+d1d2d5+d1d4d5+d2d4d5−d1+d2+d4+d52d1d2d4d51−R02>0,ifℛ0<1.
Thus, the Routh–Hurwitz conditions hold when all the eigenvalues have negative real parts makingE0 locally asymptotically stable whenever ℛ0<1.
To study the global asymptotic stability of the disease-free equilibrium, the method described by Castillo-Chavez et al. [17] and in [18] has been employed. The model system (2) is expressed in the form:(47)dX1dt=FX1,X2,dX2dt=GX1,X2,GX1,0=0,where X1=SH,SS comprises of the uninfected subpopulation and X2=IH,M,IS,C comprises of the infected subpopulation. Let U0=X1∗,0 be the disease-free equilibrium point of system (47). Then, U0 is globally asymptotically stable provided that ℛ0<1 and that the following conditions hold true:H1: For dX1/dt=FX1,0,X1∗=SH0,SS0 is globally asymptotically stableH2: GX1,X2=AX2−G^X1,X2,G^X1,X2≥0 for X1,X2∈Ω, where A=DX2GX1∗,0 is an M-matrix (the off-diagonal elements are non-negative) and Ω is the region where the model makes biological senseTheorem 5.
The disease-free equilibrium pointE0 for system (2) is globally asymptotically stable provided ℛ0<1 and conditions H1 and H2 hold true.Proof.
For system (2), we have(48)dX1dt=FX1,X2=ΛH+ηIH−1−ϕλHSH−μHSHΛS−λSSS−d3SS,dX2dt=GX1,X2=1−ϕλHSH−d1IH1−ϕργ1IH−d2MλSSS−d4ISγ2IS−d5C,FX1,0=ΛH−μHSHΛS−d3SS.
Ast⟶∞, we observe that X1,0⟶E0∈Ω, that is, SH⟶ΛH/μH and SS⟶ΛS/d3. Thus, H1 holds true:(49)A=∂GX1,0∂X2|U0=−d1001−ϕβ1SH0C01−ϕργ2−d2000β2SH0ℳ0−d4000γ2−d5,and(50)G^X1,X2=AX2−GX1,X2=β11−ϕCSH0C0−SHC0+ϵC0β2MSS0ℳ0−SSℳ0+ϵM0.
From (49), clearly A is an M-matrix (all the off-diagonal elements are non-negative), and from (50), G^X1,X2≥0 since SH0/C0≥SH/C0+ϵC and SS0/ℳ0≥SS/ℳ0+ϵM. Thus, H2 holds true, and the proof is complete.


## 3.6. Local Stability of the Endemic Equilibrium

The center manifold theory by Castillo-Chavez [19] is used to establish the local asymptotic stability of the endemic equilibrium point.Theorem 6.
Consider the following general system of ordinary differential equations with a parameterϑ:(51)dxdt=fx,ϑ,f:ℝn×ℝ⟶ℝn,f∈ℂ2ℝn×ℝ,and without loss of generality, it is assumed that 0 is an equilibrium point of the system (51) for all values of the parameter ϑ, that is,(52)f0,ϑ≡0for allϑ.Assume the following points:A1. A=Dxf0,0=∂fi/∂xj0,0 is the linearization matrix of system (51) around the equilibrium 0 with ϑ evaluated at 0. Zero is a simple eigenvalue of A, and other eigenvalues have negative real parts.A2. The matrix A has a non-negative right eigenvector w and a left eigenvector v each corresponding to the zero eigenvalue.Letfk be the kth component of f, and then, we get(53)a=∑k,i,j=1nvkwiwj∂2fk∂xi∂xj0,0,b=∑k,i=1nvkwi∂2fk∂xi∂ϑ0,0.Then, the local dynamics of the system around 0 are completely determined by the signs ofa and b that can be given as follows:(i)
a>0,b>0. When ϑ<0 with ϑ≪1, 0 is locally asymptotically stable and there exists a positive unstable equilibrium; when 0<ϑ≪1, 0 is unstable, and there exists a negative and locally asymptotically stable equilibrium.(ii)
a<0,b<0. When ϑ<0 with ϑ≪1, 0 is unstable; when 0<ϑ≪1, 0 is locally asymptotically stable, and there exists a positive unstable equilibrium.(iii)
a>0,b<0. When ϑ<0 with ϑ≪1, 0 is unstable, and there exists a locally asymptotically stable negative equilibrium; when 0<ϑ≪1, 0 is stable, and a positive unstable equilibrium appears.(iv)
a<0,b>0. When ϑ changes from negative to positive, 0 changes its stability from stable to unstable. Correspondingly, a negative unstable equilibrium becomes positive and locally asymptotically stable.Theorem 7.
The unique endemic equilibrium point,E1 of system (2) is locally asymptotically stable if R0>1 but close to 1.Proof.
We apply the change of variable on system (2) as follows. Let SH=x1,IH=x2,M=x3,SS=x4,IS=x5,C=x6. The model system (2) can be written as dX/dt=F, where X=x1,x2,x3,x4,x5,x6T and F=f1,f2,f3,f4,f5,f6:(54)dx1dt=f1=ΛH+ηx2−1−ϕβ1x6C0+ϵx6x1−μHx1,dx2dt=f2=1−ϕβ1x6C0+ϵx6x1−d1x2,dx3dt=f3=1−ϕργ1x2−d2x3,dx4dt=f4=ΛS−β2x3ℳ0+ϵx3x4−d3x4,dx5dt=f5=β2x3ℳ0+ϵx3x4−d4x5,dx6dt=f6=γ2x5−d5x6.
Takingβ1 as the bifurcation parameter while considering the case ℛ0=1 gives(55)β1=β1∗=C0ℳ0μHd1d2d3d4d5β2ργ1γ2ΛHΛS1−ϕ2.
The linearization around the disease-free equilibrium point evaluated atβ1∗ gives(56)A=−μH0000−1−ϕβ1∗ΛHC0μH0−d10001−ϕβ1∗ΛHC0μH01−ϕργ1−d200000−β2ΛSℳ0d3−d30000β2ΛSℳ0d30−d400000γ2−d5.
MatrixA has zero as an eigenvalue, and this can be deduced from (21)–(23) with ℛ0=1, which yields(57)λλ3+a1λ2+a2λ+a3=0.
The right eigenvectorw=w1,w2,w3,w4,w5,w6T corresponding to the zero eigenvalue is computed from Aw=0, which yields(58)w1=−β1∗1−ϕΛHμH2C0w6,w2=β1∗1−ϕΛHC0μHd1w6,w3=ργ1β1∗ΛH1−ϕ2C0μHd1d2w6,w4=−ργ1β1∗β2ΛSΛH1−ϕ2ℳ0C0μHd1d2d32w6,w5=d5γ2w6,w6>0.
The left eigenvectorv=v1,v2,v3,v4,v5,v6T is computed from vA=0 which yields(59)v1=v4=0,v2=C0μHd5β1∗1−ϕΛHv6,v3=C0μHd1d5ργ1β1∗ΛH1−ϕ2v6,v5=γ2d4v6,v6>0.
Lettingw6=1, it can be seen from (58) that any other right eigenvector of matrix A corresponding to the zero eigenvalue is a constant multiple of w (the same is true for the left eigenvector in (59) by letting v6=1). Hence, zero is a simple eigenvalue.
In the computation ofa and b, only nonzero partial derivatives of system (54) are used. Since v1=v4=0, the nonzero partial derivatives computed at the disease-free equilibrium are from f2,f3,f5, and f6. These are obtained as(60)∂2f2∂x1∂x6=β1∗C01−ϕ,∂2f2∂2x6=−2ϵβ1∗C021−ϕx1∗,∂2f5∂x3∂x4=β2ℳ0,∂2f5/∂2x3=−2ϵβ2x5∗/M02 and ∂2f2/∂x6∂β1∗=1−ϕ/C0x1∗. Thus, it follows that(61)a=v2w1w6∂2f2∂x1∂x6+v2w62∂2f2∂2x6+v5w3w4∂2f5∂x3∂x4+v5w32∂2f5∂2x3,a=−v6w62β1∗1−ϕd5C0μH+2ϵd5C0+d4d52ΛSγ2+2ϵd3d4d52ΛSβ2γ22<0,b=d6β1∗v7w7>0.
Sincea<0 and b>0, then by condition iv of Theorem 5, system (54) exhibits a forward bifurcation at ℛ0=1, implying that the unique endemic equilibrium is locally asymptotically stable when ℛ0>1 and close to 1.


## 3.7. Global Stability of the Endemic Equilibrium Point

Theorem 8.
Ifℛ0≥1, the endemic equilibrium E1 of system (2) is globally asymptotically stable.Proof.
To analyze the global stability of endemic equilibrium, we consider the candidate Lyapunov functionLSH,IH,M,SS,IS,C given by(62)L=SH−SH∗+SH∗lnSH∗SH+IH−IH∗+IH∗lnIH∗IH+M−M∗+M∗lnM∗M+SS−SS∗+SS∗lnSS∗SS+IS−IS∗+IS∗lnIS∗IS+C−C∗+C∗lnC∗C.
Clearly,LSH∗,IH∗,M∗.SS∗,IS∗,C∗=0, and L>0, otherwise. Moreover, L is radially unbounded in R+6. So, it remains to determine the sign of L′. Differentiating (62) with respect to t gives(63)dLdt=SH−SH∗SHdSHdt+IH−IH∗IHdIHdt+M−M∗MdMdt+SS−SS∗SSdSSdt+IS−IS∗ISdISdt+C−C∗CdCdt.
Substituting for values ofdSH/dt, dIH/dt, dM/dt, dSS/dt, dIS/dt, and dC/dt in (63) yields(64)dLdt=SH−SH∗SHΛH+ηIH−1−ϕλHSH−μHSH+IH−IH∗IH1−ϕλHSH−d1IH+M−M∗M1−ϕργ1IH−d2M+SS−SS∗SSΛS−λSSS−d3SS+IS−IS∗ISλSSS−d4IS+C−C∗Cγ2IS−d5C.
Rearranging positive and negative terms in (64) gives(65)dLdt=P−Q,where(66)P=ΛH+ηIH+1−ϕλH+μH+λS+d3SH∗+1−ϕλHSH+λSSS+d1IH∗+η+1−ϕργ1IH+d2M∗+ΛS+λSSS+d4IS∗+γ2IS+d5C∗,Q=1−ϕλH+μHSH+ΛH+ηIHSH∗SH+d1IH+1−ϕλHSHIH∗IH+d2M+1−ϕργ1IHM∗M+λS+d3SS+ΛSSS∗SS+d4IS+λSSSIS∗IS+d5C+γ2ISC∗C.
IfP≤Q, then dL/dt≤0 noting that dL/dt=0 if and only if SH=SH∗, IH=IH∗, M=M∗, SS=SS∗, IS=IS∗, and C=C∗. Therefore, the largest compact invariant set in SH,IH,M,SS,IS,C∈Ω:dL/dt=0 is the singleton E1, where E1 is the endemic equilibrium point of system (2). Thus, by Lasalle’s invariance principle [20], it implies that E1 is globally asymptotically stable in Ω if P≤Q; which holds if and only if ℛ0≥1.


## 3.8. Sensitivity Analysis

In this subsection, sensitivity analysis of the reproduction numberℛ0 on different parameters is carried out in order to determine their relative importance to the disease transmission in the community. The approach by Chitnis et al. [21] has been employed to obtain the relative change in a state variable when the parameter changes. An expression for the sensitivity of ℛ0 with respect to p is given by(67)γpℛ0=∂ℛ0∂ppℛ0,where p is any sensitive parameter to be considered in relation to ℛ0.From (67), expressions for the sensitivity indices are as follows:(68)γβ1ℛ0=γβ2ℛ0=γρℛ0=γβ1ℛ0=γγ1ℛ0=γγ2ℛ0=γλHℛ0=γλSℛ0=12,γμCℛ0=−12μCd5,γθMℛ0=−12θMd2,γμMℛ0=−12μMd1,γδHℛ0=−12δHd1=γδSℛ0=−12δSd4,γθCℛ0=−12θCd5,γθSℛ0=−12θSd3+d4d3d4,γμSℛ0=−12μSd3+d4d3d4.The parameter values used areμC=0.004,μS=0.000569,δS=0.0004012,μM=0.008,θS=0.1,η=0.0181, and μH=0.0000384.Table1 shows the sensitivity indices of ℛ0 with respect to different parameters from the most sensitive to the least sensitive. It is noted that the value of ℛ0 increases when parameter values β1, β2, ρ, γ1, γ2, ΛH, and ΛS increase, while the other parameters values are kept constant. Since these parameters have positive indices, the endemicity of the disease is increased. When parameter μH, θS, θC, μC, δH, μS, θM, and δS are increased while the other parameters values are kept constant, the value of ℛ0 decreases implying that the endemicity of the disease is decreased.Table 1
Numerical values of sensitivity indices ofℛ0.
Parameter symbolSensitivity indexμH−0.5009β1+0.5000β2+0.5000ρ+0.5000γ1+0.5000γ2+0.5000ΛH+0.5000ΛS+0.5000μM−0.4956θS−0.4952θC−0.3333μC−0.1667δH−0.0885μS−0.0056θM−0.0044δS−0.0020


## 4. Numerical Simulations

In this section, we numerically simulate the combined effect of treatment with praziquantel, public health education and chemical control strategies on schistosomiasis transmission. The model was simulated using ODE solvers coded in MATLAB computer software. The parameter values used are shown in Table2 with the initial conditions such as SH0=120000, IH0=15000, M0=200, SS0=400, IS0=200, and C0=0.Table 2
Parameter values used in numerical simulations of the schistosomiasis model.
ParameterValueSourceμH0.0000384day−1[22]β10.406×10−8day−1[10]β20.615L−1day−1[23]ρ300[6]γ10.00232day−1[23]γ22.6day−1[23]ΛH300EstimatedΛS200[24]μM0.9day−1[23]θS0.01day−1[11]θC0.008day−1[11]C09×107[24]M01×106[10]ϕ0<ϕ<1Estimatedϵ0.3[9]μC0.004day−1[15]η0.0181day−1[26]δS0.0004012day−1[24]θM0.008day−1[11]μS0.000569day−1[23]δH0.0039day−1[27]Figures2 and 3 show numerical simulation for a mathematical model of schistosomiasis with the treatment strategy only, where η=0.0181, θS=0, θM=0, θC=0, and ϕ=0.Figure 2
Graph of population against time with treatment intervention only.Figure 3
Graph of human and snail subpopulations with treatment intervention only.The results indicate that with the treatment strategy only, the subpopulations of infected humans,IH, and miracidia, M, decrease in the early stages of treatment and begin to increase again. The subpopulations of susceptible humans and susceptible snails decrease, while the infected snails and cercariae increase. It is noted that with treatment intervention only, schistosomiasis persists within the population.Figure4 shows the numerical simulation of schistosomiasis model with treatment and public health education control strategy when ϕ=0.2, 0.5, and 0.9. The results reveal that the susceptible humans and snail subpopulations increase with increase in public health awareness from ϕ=0.2 to ϕ=0.9. The subpopulations of infected humans, infected snails, miracidia, and cercariae decrease as ϕ increases. This is because of increased awareness among the human population to maintain proper sanitation as well as reduced contact with contaminated water.Figure 4
Graph of population against time with treatment intervention and public health education control strategy whenϕ=0.2, 0.5, and 0.9.Figure5 shows numerical simulation of schistosomiasis model with the treatment and chemical control intervention strategy only with parameter values η=0.0181, θS=0.01, θM=0.008, θC=0.008, and ϕ=0. From the results, it is noted that the subpopulations of susceptible humans and susceptible snails increase relatively higher than that illustrated in Figures 2 and 3 which considers treatment only. The subpopulations of infected humans, infected snails, miracidia, and cercariae also decrease further as compared to those in Figures 2 and 3.Figure 5
Graph of population against time with treatment intervention and chemical control intervention strategy.Figure6 shows numerical simulation of schistosomiasis model with treatment, public health education, and chemical control intervention strategies with parameter values η=0.0181, θS=0.01, θM=0.008, θC=0.008, and ϕ=0.7. The results indicate that the subpopulations of infected humans, miracidia, infected snails, and cercariae decrease very significantly compared to other strategies illustrated in Figures 2–5, whereas the subpopulations of susceptible humans and susceptible snail increase significantly.Figure 6
Graph of population against time with treatment intervention, public health education, and chemical control intervention strategy.


## 5. Discussion and Conclusion

In this study, a deterministic mathematical model describing the effect of treatment, public health education, and chemical control intervention strategies on schistosomiasis transmission was formulated and analyzed. The model was found to be well posed and biologically meaningful.Analysis of the model showed existence of two equilibria, namely, the disease-free equilibrium and the endemic equilibrium. The basic reproductive numberℛ0 of the model was determined, and its biological interpretation explained. ℛ0 was further decomposed into individual reproductive numbers of the transmission agents of schistosomiasis; that is, ℛ0H is the average number of secondary infections by an infected snail in a completely susceptible population and ℛ0S is the average number of secondary infections by an infected human in a completely susceptible population. The disease-free steady state was found to be locally stable whenever ℛ0<1. The disease-free equilibrium point was also globally stable implying that all initial conditions would converge to it over time; hence, treatment can decrease the disease prevalence. The endemic equilibrium point of the model was found to be locally stable if and only if ℛ0>1. This meant that interventions by policy makers are necessary to eradicate schistosomiasis over time. In order to make endemic equilibrium point unstable, control measures such as mass drug administration in endemic areas are necessary.Sensitivity analysis identifies the rate at which susceptible humans get infected,β1 and rate at which susceptible snails get infected, β2 as the main parameters fueling the spread, while the θM and θS are the elimination rates of miracidia and snails by molluscicide (chemical control) as inhibitors of the disease. Thus, from the results, it is noted that a combination of treatment, public health education, and chemical control intervention strategies can effectively manage the transmission of schistosomiasis in endemic areas.Numerical simulations reveal that a combination of treatment, public health education, and chemical control intervention strategies significantly increase the number of susceptible human and susceptible snail populations while significantly decreasing the number of infected humans, miracidia, infected snails, and cercariae as shown in Figure6. The study shows that treatment only in Figures 2 and 3 is less effective and the disease persists within the population, which agrees with the results of Gurarie et al. The study also reveals that a combination of both treatment and chemical control intervention, illustrated in Figure 5 has more impact compared to that of treatment only in Figures 2 and 3 which agrees with the results of Gao et al. [10].In conclusion, in order to effectively manage the prevalence of schistosomiasis transmission, there should be a combined control intervention of treatment, public health education, and chemical control intervention strategies.



---

*Source: 2094979-2022-12-19.xml*