# A Schistosomiasis Model with Praziquantel Resistance

Authors: Longxing Qi; Jing-an Cui
Journal: Discrete Dynamics in Nature and Society
 (2013)

---

## Abstract

A compartmental model is established for schistosomiasis with praziquantel resistance. The model considers the impact of genetic resistance and drug treatment on the transmission of schistosomiasis. We calculate the basic reproductive number and discuss the existence and stability of disease-free equilibrium, boundary equilibrium, and coexistence equilibrium. Our analysis shows that regardless of whether drug treatment leads to the emergence of resistance, once the impact of genetic resistance is larger, the resistant strain will be dominant, which is detrimental to the control of schistosomiasis. In addition, once the proportion of human with drug-resistant strain produced by drug treatment is larger, the number of human and snails with resistant strain is larger. This is not a good result for drug treatment with praziquantel.


---

## Body

## 1. Introduction

Currently, treatment of human beings infected by schistosomiasis primarily focuses on chemotherapy with praziquantel (PZQ). PZQ appeared as a new schistosomicidal compound during the 1970s [1]. In recent years, PZQ has become the drug of choice in most endemic areas because of its efficacy, its ease of administration, its tolerable side-effects, and its cost [1]. Although the effectiveness of PZQ against schistosomiasis is well documented, the precise mode of action of the drug has not clearly defined [2]. It is reported that the chemotherapy of many helminth infections is complicated by the occurrence of drug resistance and drug tolerance (a natural resistance) to certain anthelmintics [2]. Not surprisingly, recent epidemiological evidence suggests the emergence of PZQ-resistant∖tolerant schistosomes [1, 3, 4]. Resistance is defined as a genetically transmitted loss of sensitivity in a parasite population that was previously sensitive to a given drug [2]. Tolerance is an innate insusceptibility of a parasite to a drug, with the caveat that the parasite must not have been previously exposed to the drug [2]. The first report of possible PZQ resistance came from an intensive focus in northern Senegal, where the drug had produced very low cure rates (18–39%) [5, 6]. And snails collected in the area carried schistosoma strains. When tested in the laboratory, those snails had a decreased susceptibility to PZQ [7, 8]. Additional PZQ-resistant evidence was collected in Egypt [4, 9]. Preliminary studies have begun on these isolates to identify genetic, physiological, and morphological characteristics associated with PZQ resistance, and some of these may find use as markers for monitoring whether or not resistance is developing in endemic areas, where the drug is used [10].Many papers have reported that drug treatment results in the emergence of schistosome resistance to PZQ [11–13]. Drug treatment can remove drug-susceptible parasites in infected human beings, while resistant parasites survive. However, many investigations find that traits of PZQ resistance of Schistosoma mansoni are dominant inheritance [14–16]. The resistant worms can reproduce and pass the resistant genes to the next generation. Furthermore, the resistance of Schistosoma mansoni to PZQ can be expressed in eggs, miracidia, cercariae, adults, and all stages of development [14–16]. In other words, some definitive hosts carrying resistant schistosomes can infect snails and make those snails carry resistant schistosomes. On the contrary, some snails carrying resistant schistosomes can infect definitive hosts and also make those definitive hosts carry resistant schistosomes. It goes on, the control of schistosomiasis will face enormous difficulties. Therefore, it is necessary to study the impact of this schistosome genetic resistance on the transmission of schistosomiasis.In previous schistosomiasis models, resistant problems have been studied by considering that the resistance of schistosomiasis is due to drug treatment [11–13]. In [11, 13], the authors proposed a multistrain schistosome model including sensitive and resistant parasite strains. Their goal was to infer the impact of drug treatment on the maintenance of schistosome genetic diversity. In their assumptions, the drug-sensitive parasite strain had an additional per capita death rate, σ, due to treatment. For a parasite strain that had developed drug resistance with a resistance level θ (θ>1), this treatment-related death rate was assumed to be reduced by the factor θ to σ/θ. Their results implied that higher treatment rate could allow for coexistence between sensitive and resistant parasite strains. In [12], the authors formulated a deterministic model with multiple strains of schistosomes in order to explore the role of drug treatment in the maintenance of a polymorphism of parasite strains that differed in their resistance levels. And snails infected by parasite strains were divided into multistrain subclasses according to the different level θi of parasite strains. Analysis of the model showed that the likelihood that resistant strains would increase in frequency depended on the interplay between their relative fitness, the cost of resistance, and the degree of selection pressure exerted by drug treatments.Motivated by [11–13], we establish a new model considering hosts with sensitive and resistant strains in this paper. Our purpose is mainly to study the impact of drug treatment and genetic resistance on the transmission of schistosomiasis.Our paper is organized as follows. In Section2, we establish a mathematical model with praziquantel resistance and obtain basic reproductive number and existence of equilibria. And then the stability of disease-free equilibrium is obtained in Section 3. Section 4 devotes to stability analysis of boundary equilibria. In Section 5, stability analysis of endemic equilibrium is performed.


## 2. Mathematical Model

According to different level of parasite strains, we divide infected hosts into sensitive and resistant strains. Considering resistance and inheritance of resistance, new resistant strains are composed of two parts. We classify definitive and intermediate hosts as susceptible, sensitive and resistant in the following:(i)
Xs(t), the population of susceptible human;(ii)
Xi(t), the population of human infected with sensitive parasite strain;(iii)
Xr(t), the population of human infected with resistant parasite strain;(iv)
Ys(t), the population of susceptible snail host;(v)
Yi(t), the population of snail host carrying sensitive parasites;(vi)
Yr(t), the population of snail host carrying resistant parasites.We follow some of the available models for schistosomiasis [17–19] and assume that the reproduction rate of hosts is constant, and we ignore the recovery class of host since the life span of infected is short in comparison to that of human [18, 20]. All parameters in the model are assumed to be nonnegative constants:(i)
Λ1, recruitment rate of human;(ii)
μ1, per capita natural death rate of human;(iii)
α1, per capita disease-induced death rate of human;(iv)
β1, per capita contact transmission rate from infected snails to susceptible human;(v)
Λ2, recruitment rate of snails;(vi)
μ2, per capita natural death rate of snails;(vii)
α2, per capita disease-induced death rate of snails;(viii)
β2, per capita contact transmission rate from infected human to susceptible snails;(ix)
θ, the rate of treatment for infected human;(x)
ρ, the proportion of human with drug-resistant strain produced by treatment;(xi)
k, represents the impact of inheritance and the cost of resistance on transmission rate, we assume that k≤1 since the transmission rate is reduced due to resistance [12];(xii)
k1, represents the impact of resistance on disease-induced death rate of human;(xiii)
k2, represents the impact of resistance on disease-induced death rate of snails.Then, we have a model with the form:(1)dXsdt=Λ1-β1XsYi-kβ1XsYr-μ1Xs+(1-ρ)σXi,dXidt=β1XsYi-μ1Xi-α1Xi-σXi,dXrdt=kβ1XsYr+ρσXi-μ1Xr-k1α1Xr,dYsdt=Λ2-β2YsXi-kβ2YsXr-μ2Ys,dYidt=β2YsXi-μ2Yi-α2Yi,dYrdt=kβ2YsXr-μ2Yr-k2α2Yr.
Using standard methods, it is easy to see that disease free equilibrium E0=(Λ1/μ1,0,0,Λ2/μ2,0,0) always exists. Let Xs0=Λ1/μ1, Ys0=Λ2/μ2, m1=μ1+α1+σ, m2=μ1+k1α1, m3=μ2+α2, and m4=μ2+k2α2. According to the concept of next generation matrix [21] and the formula of the basic reproductive number for ODE compartmental models [22], if we let
(2)F=(00β1Xs00000kβ1Xs0β2Ys00000kβ2Ys000),V=(m10000m20000m30000m4),
one can calculate that the eigenvalues of the next generation matrix FV-1 are given by
(3)R0i=β1β2Xs0Ys0m1m3,R0r=k2β1β2Xs0Ys0m2m4.
Then, it follows that the basic reproductive number for the system (1) is given by
(4)R0=max{R0i,R0r}.To obtain other equilibria, we let the right-hand side of (1) equal to zero and obtain
(5)Λ1-β1XsYi-kβ1XsYr-μ1Xs+(1-ρ)σXi=0,(6)β1XsYi-μ1Xi-α1Xi-σXi=0,(7)kβ1XsYr+ρσXi-μ1Xr-k1α1Xr=0,(8)Λ2-β2YsXi-kβ2YsXr-μ2Ys=0,(9)β2YsXi-μ2Yi-α2Yi=0,(10)kβ2YsXr-μ2Yr-k2α2Yr=0.
If Xi=0, we have Yi=0 and a formula of Xr as follows:
(11)(k2β1β2Λ2m2+μ1m2m4kβ2)Xr=k2β1β2Λ1Λ2-μ1μ2m2m4.Note thatR0r>1 equals to k2β1β2Λ1Λ2>μ1μ2m2m4. Hence, if R0r>1, there exists a boundary equilibrium with only resistant type, given as E0r=(Xsr,0,Xrr,Ysr,0,Yrr), where
(12)Xsr=Λ1(μ1m4R0r2+kβ1Λ2)(μ12m4+kβ1μ1Λ2)R0r2,Xrr=μ1μ2m4(R0r2-1)(μ1m4+kβ1Λ2)kβ2,Ysr=Λ2(μ1m4+kβ1Λ2)μ1μ2m4R0r2+kβ1Λ2μ2,Yrr=Λ2μ1(R0r2-1)μ1m4R0r2+kβ1Λ2.If drug treatment does not lead to drug resistance, that is,ρ=0, then when R0i>1, we can obtain the other boundary equilibrium with only sensitive type, given as E0i=(Xsi,Xii,0,Ysi,Yii,0), where
(13)Xsi=Λ1m1m3+(μ1+α1)μ2m1m3β1β2Λ2(μ1+α1)+μ1m1m3,Xii=μ1μ2m1m3(R0i2-1)β1β2Λ2(μ1+α1)+μ1m1m3,Ysi=β1β2Λ22(μ1+α1)+Λ2μ1m1m3β2μ1μ2m1m3(R0i2-1)+β1β2Λ2μ2(μ1+α1)+μ1μ2m1m3,Yii=Λ2β2μ1μ2m1(R0i2-1)β2μ1μ2m1m3(R0i2-1)+β1β2Λ2μ2(μ1+α1)+μ1μ2m1m3.Now, we study existence of coexistence equilibrium for the system (1). From (5)–(10), we obtain
(14)[β1β2Λ2(μ1+α1)+μ1m1m3β2]Xi+(μ1m1m3kβ2+β1β2Λ2m2)Xr=β1β2Λ1Λ2-μ1μ2m1m3.
Following (6), (7), (9), and (10), we have
(15)YiYr=km1Xim2Xr-ρσXi=m4Xikm3Xr.In the case thatρ=0, (15) leads to
(16)km1Xim2Xr=m4Xikm3Xr.
Note that R0i=R0r equals to km1/m2=m4/km3. It is easy to see that if R0i≠R0r, (16) cannot hold, which implies that there is not coexistence equilibrium. If R0i=R0r, (16) always holds. Hence, if R0r=R0i>1, that is, β1β2Λ1Λ2>μ1μ2m1m3, (14) represents the existence of coexistence equilibrium in the form of a line.In the case thatρ≠0, (15) leads to
(17)m4ρσXi=(m2m4-k2m1m3)Xr.
Note that R0i>R0r equals to m2m4-k2m1m3>0. Following (14) and (17), we can obtain that if R0i>R0r and R0i>1, then the unique coexistence equilibrium is given by Ec=(Xsc,Xic,Xrc,Ysc,Yic,Yrc), where
(18)Xrc=m4ρσ(β1β2Λ1Λ2-μ1μ2m1m3)Xrc=×([β1β2Λ2(μ1+α1)+μ1m1m3β2]Xrc=×(m2m4-k2m1m3)Xrc=+m4ρσ(μ1m1m3kβ2+β1β2Λ2m2))-1,Xic=m2m4-k2m1m3m4ρσXrc,Ysc=Λ2β2Xic+kβ2Xrc+μ2,Yrc=kΛ2β2Xrcm4(β2Xic+kβ2Xrc+μ2),Yic=β2YscXicμ2+α2,Xsc=Λ1+(1-ρ)σXicβ1Yic+kβ1Yrc+μ1.Summarizing above analyses, we have the following result.Lemma 1.
The existence of equilibria for the system (1) is as follows:(1)
the disease free equilibriumE0 always exists;(2)
ifR0r>1, there exists a boundary equilibrium with only resistant type E0r;(3)
ifρ=0 and R0i>1, there exists a boundary equilibrium with only sensitive type E0i;(4)
ifρ=0 and R0i=R0r>1, there exists coexistence equilibrium in the form of a line (14);(5)
ifρ≠0 and R0i>R0r and R0i>1, there exists a unique coexistence equilibrium Ec.The following section shows that the basic reproductive numberR0 provides a threshold condition for schistosoma extinction in (1).


## 3. Stability Analysis of the Disease Free Equilibrium

In this section, we will analyze stability of the disease free equilibrium of the model (1). The stability of the disease free equilibrium determines whether schistosomiasis will be permanent in an uninfected population. The following result shows that schistosome will go extinct if R0<1.Theorem 2.
The disease free equilibriumE0 of the system (1) is locally asymptotically stable if R0<1 and unstable if R0>1.Proof.
The Jacobian matrix for the system (1) is given by(19)J=(-β1Yi-kβ1Yr-μ1(1-ρ)σ00-β1Xs-kβ1Xsβ1Yi-m100β1Xs0kβ1Yrρσ-m200kβ1Xs0-β2Ys-kβ2Ys-β2Xi-kβ2Xr-μ2000β2Ys0β2Xi-m3000kβ2Yskβ2Xr0-m4).Then, the eigenvalues of E0 are -μ1, -μ2 and roots of the following equations:
(20)λ2+(m1+m3)λ+m1m3-β1β2Xs0Ys0=0,λ2+(m2+m4)λ+m2m4-k2β1β2Xs0Ys0=0.
Note that R0<1 equals to R0i<1 and R0r<1, which leads to m1m3-β1β2Xs0Ys0>0 and m2m4-k2β1β2Xs0Ys0>0. Hence, if R0<1, then all the roots of (20) have negative real parts. Hence, using the Routh-Hurwitz criterion, we can obtain that the disease free equilibrium E0 of the system (1) is locally asymptotically stable if R0<1 and unstable if R0>1.Now, we turn to the study of the global stability of the disease free equilibrium of the model (1) by using Metzler matrix theory and the technique of Kamgang and Sallet [23].Consider systems of the following form:(21)dx1dt=f(x1,x2),dx2dt=g(x1,x2),
where x1∈R+n1, x2∈R+n2, and f and g are C1. We denote by x=(x1,x2) the state of the system and (x1*,0) is a disease free equilibrium on a positively invariant set Ω⊂R+n1+n2. Now rewrite (21) as
(22)dx1dt=A1(x)·(x1-x1*)+A12(x)·x2,dx2dt=A2(x)x2.
For the system (22), we make the following assumptions.(h1)
The system is defined on the positively invariant setΩ of the nonnegative orthant. The system is dissipative on Ω.(h2)
The subsystemdx1/dt=A1(x)·(x1-x1*) is globally asymptotically stable at the equilibrium x1* on the canonical projection of Ω on R+n1.(h3)
The matrixA2(x) is Metzler and irreducible for any given x∈Ω.(h4)
There exists an maximum matrixA-2, then for any x-∈Ω such that A-2=A2(x-), x-∈R+n1×{0}.(h5)
α(A-2)≤0, that is, the greatest real part of eigenvalues of A-2 is nonnegative.For convenience, we state two lemmas due to Kamgang and Sallet [23].Lemma 3.
If the above hypotheses, h1–h5, are satisfied, then the disease free equilibrium is globally asymptotically stable in Ω.Lemma 4.
If the same notations and hypotheses in Lemma3 hold and if, furthermore, we have A-2=A2(x1*,0), the disease free equilibrium is globally asymptotically stable if and only if α(A-2)≤0.Next, we discuss the global stability of the disease free equilibriumE0 of the system (1) using the above two Lemmas.From the system (1), we know
(23)d(Xs+Xi+Xr)dt=Λ1-μ1(Xs+Xi+Xr)-α1Xi-k2α1Xr≤Λ1-μ1(Xs+Xi+Xr),d(Ys+Yi+Yr)dt=Λ2-μ2(Ys+Yi+Yr)-α2Yi-k4α2Yr≤Λ2-μ2(Ys+Yi+Yr).
This proves that the set
(24)Ω={Λ2μ2(Xs,Xi,Xr,Ys,Yi,Yr)∈ℝ+5∣Xs+Xi+Xraaaa≤Λ1μ1,Ys+Yi+Yr≤Λ2μ2}
is a compact positively invariant absorbing set contained in the nonnegative orthant. Thus, the system (1) is dissipative on Ω because the trajectories of (1) are forward bounded. Now, we will study the system (1) on Ω.We set for system (1) x1=(Xs,Ys), x2=(Xi,Xr,Yi,Yr), and x1*=(Λ1/μ1,Λ2/μ2). As in [23], we express the subsystem as dx1/dt=A1(x1,0)(x1-x1*) and
(25)dXsdt=Λ1-μ1Xs,dYsdt=Λ2-μ2Ys.
This is a linear system, and its unique equilibrium (A1/μ1,A2/μ2) (corresponding to the disease free equilibrium of (1)) is globally asymptotically stable, hence the assumptions (h1) and (h2) are satisfied.The matrixA2(x) is given by
(26)A2(x)=(-m10β1Xs0ρσ-m20kβ1Xsβ2Ys0-m300kβ2Ys0-m4).
As required by hypothesis h3, for any x∈Ω, the matrix A2(x) is irreducible.Now, let us check (h4). There is a maximum which is uniquely realized in Ω if Xs=Λ1/μ1 and Ys=Λ2/μ2, which corresponds to the disease free equilibrium. This maximum matrix is then J2, the subblock of the Jacobian matrix at the disease free equilibrium, corresponding to the matrix A2(x). The matrix J2 is given by
(27)J2=(-m10β1Xs00ρσ-m20kβ1Xs0β2Ys00-m300kβ2Ys00-m4).
Therefore, we are in the situation of Lemma 4, where the maximum is attained at the disease free equilibrium.The hypothesis (h5) requires that α(J2)≤0. Writing J2 as a block matrix J2=(A2×2B2×2C2×2D2×2). Since A is already a Metzler stable matrix, the condition α(J2)≤0 is equivalent to the condition α(D-CA-1B)≤0 [23], where
(28)D-CA-1B=(β1β2Xs0Ys0m1-m30ρσkβ1β2Xs0Ys0m1m2k2β1β2Xs0Ys0m2-m4).
Then,
(29)α(D-CA-1B)=max{β1β2Xs0Ys0-m1m3m1,k2β1β2Xs0Ys0-m2m4m2}.
Hence, the condition α(D-CA-1B)≤0 is equivalent to R0≤1. We have seen that the hypotheses (h1), (h2), (h3), (h4), and (h5) are satisfied. Then, by Lemma 4, we have the following result.Theorem 5.
The disease free equilibriumE0 of the system (1) is globally asymptotically stable if R0≤1.


## 4. Stability Analysis of the Boundary Equilibria

In this section, we turn to study stability of the two boundary equilibria. From Lemma1, we know that if R0r>1, there exists a boundary equilibrium with only resistant type E0r. Through calculations, we can obtain the characteristic equation as following:
(30){(λ+m1)(λ+m3)-β1β2XsrYsr}×{k2β1β2XsrYsr(λ+μ1)(λ+μ2)(λ+m2)(λ+m4)(λ+μ1+kβ1Yrr)aaaaaaa×(λ+μ2+kβ2Xrr)aaaaaaa-k2β1β2XsrYsr(λ+μ1)(λ+μ2)}=0.
Hence, the eigenvalue of E0r are roots of the following equations:
(31)λ2+(m1+m3)λ+m1m3-β1β2XsrYsr=0,(32)λ4+a1λ3+a2λ2+a3λ+a4=0.
Here,
(33)a1=m2+m4+Λ1Xsr+Λ2Ysr>0,a2=m2m4+(m2+m4)(Λ1Xsr+Λ2Ysr)a2=+Λ1XsrΛ2Ysr-k2β1β2XsrYsr,a3=m2m4(Λ1Xsr+Λ2Ysr)+(m2+m4)Λ1XsrΛ2Ysra3=-k2β1β2XsrYsr(μ1+μ2),a4=m2m4Λ1XsrΛ2Ysr-k2β1β2XsrYsrμ1μ2.From (12), we can obtain XsrYsr=Xs0Ys0/R0r2=m2m4/k2β1β2. Then,
(34)m1m3-β1β2XsrYsr=m1m3-m2m4k2.
It is easy to see that R0i2/R0r2=(m2m4/k2)/m1m3. Then, R0i<R0r equals to m1m3>β1β2XsrYsr, which implies that the roots of (31) have negative real parts if R0i<R0r.Based on(Λ1/Xsr)(Λ2/Ysr)=μ1μ2R0r2 and k2β1β2XsrYsr=m2m4, we can obtain
(35)a2=(m2+m4)(Λ1Xsr+Λ2Ysr)+Λ1XsrΛ2Ysr>0,a3=m2m4[kβ1μ1Λ2(R0r2-1)μ1m4R0r2+kβ1Λ2+μ1μ2m4(R0r2-1)μ1m4+kβ1Λ2]a3=+(m2+m4)μ1μ2R0r2,a4=m2m4μ1μ2(R0r2-1).
Here, a3>0 and a4>0 if R0r>1.IfR0r>1, let kβ1μ1Λ2(R0r2-1)/(μ1m4R0r2+kβ1Λ2)=p1>0 and μ1μ2m4(R0r2-1)/(μ1m4+kβ1Λ2)=p2>0, then Λ1/Xsr=p1+μ1>μ1 and Λ2/Ysr=p2+μ2>μ2. Note that m2=μ1+k1α1>μ1 and m4=μ2+k2α2>μ2, then m2+m4>μ1+μ2 and m2m4>μ1μ2. Hence, we can obtain
(36)H1=a1>0,H2=|a1a31a2|=a1a2-a3=(m22+m42+m2m4+μ1μ2R0r2)(Λ1Xsr+Λ2Ysr)+(m2+m4)(Λ1Xsr+Λ2Ysr)2+m2m4(μ1+μ2)>0,H3=|a1a301a200a1a3|=a1a2a3-a32-a12a4=(m2+m4)2m2m4(Λ1Xsr+Λ2Ysr)2+(m2+m4)3(Λ1Xsr+Λ2Ysr)μ1μ2R0r2+μ1μ2R0r2(m2+m4)m2m4(Λ1Xsr+Λ2Ysr)+(m2+m4)m2m4(Λ1Xsr+Λ2Ysr)3+(m2+m4)2(Λ1Xsr+Λ2Ysr)2μ1μ2R0r2+(m2+m4)(Λ1Xsr+Λ2Ysr)μ12μ22R0r4+2(μ1+μ2)m22m42(Λ1Xsr+Λ2Ysr)+(m2+m4)μ1μ2R0r2m2m4(μ1+μ2)+(m2+m4)2μ1μ2m2m4+2(m2+m4)(Λ1Xsr+Λ2Ysr)μ1μ2m2m4+(Λ1Xsr+Λ2Ysr)2μ1μ2m2m4-(m2+m4)2m2m4(Λ1Xsr+Λ2Ysr)(μ1+μ2)-(m2+m4)m2m4(Λ1Xsr+Λ2Ysr)2(μ1+μ2)-μ1μ2R0r2(μ1+μ2)m2m4(Λ1Xsr+Λ2Ysr)-(Λ1Xsr+Λ2Ysr)2m22m42-m22m42(μ1+μ2)2-2m2m4(m2+m4)(Λ1Xsr+Λ2Ysr)μ1μ2R0r2-m2m4(m2+m4)2μ1μ2R0r2>(m22+m42)2m2m4(p12+p22+p1p2+p1μ1+p2μ2)+m22m42(p12+p22+2p1μ1+2p2μ2)>0,H4=|a1a3001a2a400a1a3001a2a4|=a4H3>0.
It follows from Routh-Hurwitz criterion that all roots of (32) have negative real parts if R0r>1. Summering above analyses, we have the following result.Theorem 6.
The boundary equilibriumE0r of the system (1) is locally asymptotically stable if R0r>1 and R0i<R0r.Now, we study the global stability of the boundary equilibriumE0r. Consider the Lyapunov function
(37)V1=d1Xi+d2Yi,
where d1=(m3+β2(Λ2/μ2))/m1m3(R0i+1) and d2=(m1+β1(Λ1/μ1))/m1m3(R0i+1). The Lyapunov derivative is
(38)dV1dt=d1(β1XsYi-m1Xi)+d2(β2YsXi-m3Yi)≤(d1β1Λ1μ1-d2m3)Yi+(d2β2Λ2μ2-d1m1)Xi≤(R0i-1)(Yi+Xi)≤0forR0i≤1.
This implies that the sensitive type dies out if R0i≤1. Then the largest compact invariant set of the system (1) in the set {dV1/dt=0} is Ω1=(Xs,Xi,Xr,Ys,Yi,Yr)∈Ω∣Xi=Yi=0. Using the LaSalle-Lyapunov theorem, we know that all trajectories in Ω eventually tend to Ω1 as t→∞. Then, we only need to study the dynamical behavior of (1) in Ω1. At this time, (1) reduces to the following system
(39)dXsdt=Λ1-kβ1XsYr-μ1Xs,dXrdt=kβ1XsYr-m2Xr,dYsdt=Λ2-kβ2YsXr-μ2Ys,dYrdt=kβ2YsXr-m4Yr.
To show that all trajectories of (39) in the interior of Ω1 approach the point (Xsr,Xrr,Ysr,Yrr) corresponding to the boundary equilibrium E0r, consider the Lyapunov function
(40)V2=d-1{[Xs-Xsr-Xsrln(XsXsr)]aaaa+[Xr-Xrr-Xrrln(XrXrr)]}+d-2{[Ys-Ysr-Ysrln(YsYsr)]aaaaa+[Yr-Yrr-Yrrln(YrYrr)]},
where positive constants d-1 and d-2 are defined in the following. It is easy to see that V2≥0 for (Xs,Xr,Ys,Yr)∈Ω1, and V2=0⇔(Xs,Xr,Ys,Yr)=(Xsr,Xrr,Ysr,Yrr). Hence, the function V2 is positive definite with respect to the point (Xsr,Xrr,Ysr,Yrr).Computing the derivative ofV2 along solutions of system (39), we have
(41)dV2dt=d-1{[dXsdt-dXsdt(XsrXs)]+[dXrdt-dXrdt(XrrXr)]}+d-2{[dYsdt-dYsdt(YsrYs)]+[dYrdt-dYrdt(YrrYr)]}=d-1[Λ1-μ1Xs-XsrXsΛ1+kβ1XsrYr+μ1Xsraaaa-m2Xr-kβ1XsYrXrrXr+m2Xrr]+d-2[Λ2-μ2Ys-YsrYsΛ2+kβ2YsrXr+μ2Ysraaaaa-m4Yr-kβ2YsXrYrrYr+m4Yrr].
Substituting Λ1=kβ1XsrYrr+μ1Xsr, Λ2=kβ2YsrXrr+μ2Ysr, m2Xrr=kβ1XsrYrr and m4Yrr=kβ2YsrXrr into dV2/dt, we obtain
(42)dV2dt=d-1{(Xsr)2Xsμ1Xsr(2-XsXsr-XsrXs)+(kβ1XsrYr-m2Xr)+[2kβ1XsrYrr-kβ1Yrr(Xsr)2Xs-kβ1XsYrXrrXr]}+d-2{(Ysr)2Ysμ2Ysr(2-YsYsr-YsrYs)+(kβ2YsrXr-m4Yr)+[2kβ2YsrXrr-kβ2Xrr(Ysr)2Ys-kβ2YsXrYrrYr]}.
Let d-1=kβ2YsrXrr and d-2=kβ1XsrYrr. Following from m2Xrr=kβ1XsrYrr and m4Yrr=kβ2YsrXrr, we have
(43)d-1(kβ1XsrYr-m2Xr)+d-2(kβ2YsrXr-m4Yr)=0.
Hence, for all (Xs,Xr,Ys,Yr)∈Ω1,
(44)dV2dt=d-1{μ1Xsr(2-XsXsr-XsrXs)aaa+kβ1XsrYrr[2-XsrXs-XsYrXrrXsrYrrXr]}+d-2{μ2Ysr(2-YsYsr-YsrYs)aaaaa+kβ2YsrXrr[2-YsrYs-YsXrYrrYsrXrrYr]}=d-1μ1Xsr(2-XsXsr-XsrXs)+d-2μ2Ysr(2-YsYsr-YsrYs)+d-1d-2(4-XsrXs-XsYrXrrXsrYrrXr-YsrYs-YsXrYrrYsrXrrYr)≤0.
It is easy to see that
(45)dV2dt=0⟺Xs=Xsr,Xr=aXrr,Ys=Ysr,Yr=aYrr,
where a is an arbitrary positive number. Substituting Xs=Xsr, Xr=aXrr, Ys=Ysr, and Yr=aYrr into the first equation of system (39), we obtain
(46)0=Λ1-akβ1XsrYrr-μ1Xsr,
and then the above formula holds if and only if a=1. Therefore, the only compact invariant subset of the set where dV2/dt=0 is the point (Xsr,Xrr,Ysr,Yrr), corresponding to the boundary equilibrium E0r. By LaSalle’s Invariance Principle, E0r is globally asymptotically stable if R0r>1 and R0i≤1. Summering above analyses, we have the following result.Theorem 7.
The boundary equilibriumE0r of the system (1) is globally asymptotically stable if R0r>1≥R0i.Now, we turn to the other boundary equilibrium. From Lemma1, we know that if ρ=0 and R0i>1, the boundary equilibrium with only sensitive type E0i exists. Through calculations, we can obtain the characteristic equation as following:
(47){(λ+m2)(λ+m4)-k2β1β2XsiYsi}×{(λ+m1)(λ+m3)(λ+μ1)(λ+μ2+β2Xii)aaaaaaaa+β1Yii(λ+m3)(λ+μ2+β2Xii)(λ+m1-σ)aaaaaaaa-β1β2XsiYsi(λ+μ1)(λ+μ2)}=0.
Hence, the eigenvalue of E0i are roots of the following equations:
(48)λ2+(m2+m4)λ+m2m4-k2β1β2XsiYsi=0,(49)λ4+b1λ3+b2λ2+b3λ+b4=0.
Here,
(50)b1=m1+m3+μ1+μ2+β2Xii+β1Yii>0,b2=μ1(m1+m3+μ2+β2Xii)+m1(m3+μ2+β2Xii)b2=+m3(μ2+β2Xii)+β1Yii(m1-σ+m3+μ2+β2Xii)b2=-m1m3b2=μ1(m1+m3+μ2+β2Xii)+m1(μ2+β2Xii)b2=+m3(μ2+β2Xii)+β1Yii(μ1+α+m3+μ2+β2Xii)b2>0,b3=μ1m1(m3+μ2+β2Xii)+m1m3(μ2+β2Xii)b3=+β1Yii[m3(μ2+β2Xii)+(m1-σ)(m3+μ2+β2Xii)]b3=-m1m3(μ1+μ2)b3=μ1m1(μ2+β2Xii)+m1m3β2Xiib3=+β1Yii[m3(μ2+β2Xii)+(μ1+α)(m3+μ2+β2Xii)]b3>0,b4=μ1m1m3(μ2+β2Xii)+β1Yiim3(μ2+β2Xii)(m1-σ)b4=-m1m3μ1μ2b4=μ1m1m3β2Xii+β1Yiim3(μ2+β2Xii)(μ1+α)>0.From (13), we can obtain that β1β2XsiYsi=m2m4. Then, it is easy to see that the roots of (48) have negative real parts if R0i>R0r. Similarly to the case of E0r, using Routh-Hurtwitz criterion, we can obtain that all roots of (49) have negative real parts if R0i>1. Summering above analysis, we have the following result.Theorem 8.
Whenρ=0, the boundary equilibrium E0i of the system (1) is locally asymptotically stable if R0i>1 and R0i>R0r.Now, we study the global stability of the boundary equilibriumE0i. Consider the Lyapunov function
(51)L1=l1Xr+l2Yr,
where l1=(m4+kβ2(Λ2/μ2))/m2m4(R0r+1) and l2=(m2+kβ1(Λ1/μ1))/m2m4(R0r+1). The Lyapunov derivative is
(52)dL1dt=l1(kβ1XsYr-m2Xr)+l2(kβ2YsXr-m3Yr)≤(l1kβ1Λ1μ1-l2m4)Yr+(l2kβ2Λ2μ2-l1m2)Xr≤(R0r-1)(Yr+Xr)≤0forR0r≤1.
This implies that the resistant type dies out if R0r≤1. Then, the largest compact invariant set of the system (1) in the set {dL1/dt=0} is Ω2={(Xs,Xi,Xr,Ys,Yi,Yr)∈Ω∣Xr=Yr=0}. Using the LaSalle-Lyapunov theorem, we know that all trajectories in Ω eventually tend to Ω2 as t→∞. Then, we only need to study the dynamical behavior of (1) in Ω2. At this time, (1) reduces to the following system:
(53)dXsdt=Λ1-β1XsYi-μ1Xs+σXi,dXidt=β1XsYi-m1Xi,dYsdt=Λ2-β2YsXi-μ2Ys,dYidt=β2YsXi-m3Yi.
To show that all trajectories of (53) in the interior of Ω2 approach the point (Xsi,Xii,Ysi,Yii) corresponding to the boundary equilibrium E0i, consider the Lyapunov function
(54)L2=l-1{[Xs-Xsi-Xsiln(XsXsi)]aaaa+[Xi-Xii-Xiiln(XiXii)]}+l-2{[Ys-Ysi-Ysiln(YsYsi)]aaaaa+[Yi-Yii-Yiiln(YiYii)]},
where positive constants l-1 and l-2 are defined in the following. It is easy to see that L2≥0 for (Xs,Xi,Ys,Yi)∈Ω1 and L2=0⇔(Xs,Xi,Ys,Yi)=(Xsi,Xii,Ysi,Yii). Hence, the function L2 is positive definite with respect to the point (Xsi,Xii,Ysi,Yii).Computing the derivative ofL2 along solutions of system (39), we have
(55)dL2dt=l-1{[dXsdt-dXsdt(XsiXs)]+[dXidt-dXidt(XiiXi)]}+l-2{[dYsdt-dYsdt(YsiYs)]+[dYidt-dYidt(YiiYi)]}=l-1[Λ1-μ1Xs-XsiXsΛ1+β1XsiYi+μ1Xsi-m1Xiaaa-β1XsYiXiiXr+m1Xii+σXi(1-XsiXs)]+l-2[Λ2-μ2Ys-YsiYsΛ2+β2YsiXi+μ2Ysi-m3Yiaaaaa-β2YsXiYiiYi+m3Yii].
Substituting Λ1=β1XsiYii+μ1Xsi-σXii, Λ2=β2YsiXii+μ2Ysi, m1Xii=β1XsiYii, and m3Yii=β2YsiXii into dL2/dt, we obtain
(56)dL2dt=l-1{(Xsi)2Xsμ1Xsi(2-XsXsi-XsiXs)+(β1XsiYi-m1Xi)+σ(Xi-Xii)(1-XsiXs)+[2β1XsiYii-β1Yii(Xsi)2Xs-β1XsYiXiiXi]}+l-2{(Ysi)2Ysμ2Ysi(2-YsYsi-YsiYs)+(β2YsiXi-m3Yi)+[2β2YsiXii-β2Xii(Ysi)2Ys-β2YsXiYiiYi]}.
Note that for the limiting system Xs+Xi=Xsi+Xii=Λ1/μ1, then (Xi-Xii)(Xs-Xsi)≤0. Let l-1=β2YsiXii and l-2=β1XsiYii. Following from m1Xii=β1XsiYii and m3Yii=β2YsiXii, we have
(57)l-1(β1XsiYi-m1Xi)+l-2(β2YsiXi-m3Yi)=0.
Hence, for all (Xs,Xi,Ys,Yi)∈Ω2,
(58)dL2dt≤l-1{μ1Xsi(2-XsXsi-XsiXs)aaa+β1XsiYii[2-XsiXs-XsYiXiiXsiYiiXi]}+l-2{μ2Ysi(2-YsYsi-YsiYs)aaaaa+β2YsiXii[2-YsiYs-YsXiYiiYsiXiiYi]}≤l-1μ1Xsi(2-XsXsi-XsiXs)+l-2μ2Ysi(2-YsYsi-YsiYs)+l-1l-2(4-XsiXs-XsYiXiiXsiYiiXi-YsiYs-YsXiYiiYsiXiiYi)≤0.
It is easy to see that
(59)dL2dt=0⟺Xs=Xsi,Xi=Xii,Ys=Ysi,Yi=Yii.
Therefore, the only compact invariant subset of the set where dL2/dt=0 is the point (Xsi,Xii,Ysi,Yii), corresponding to the boundary equilibrium E0i. By LaSalle’s Invariance Principle, E0i is globally asymptotically stable if R0i>1 and R0r≤1. Summering above analysis, we have the following result.Theorem 9.
Whenρ=0, the boundary equilibrium E0i of the system (1) is globally asymptotically stable if R0i>1≥R0r.


## 5. Stability Analysis of the Coexistence Equilibrium

In this section, we turn to study the local stability of the coexistence equilibriumEc in the limiting system of (1) by using Krasnoselskii sublinearity trick [24], as in [25, 26]. In detail, if x′=f(x) is a system of differential equations and x* is an equilibrium point, then to prove the local asymptotical stability of x* is to prove that the linearized equation Z′=f′(x*)Z has no solutions of the form
(60)Z(t)=Z0eωt,
with Z0∈Cn, ω∈C and, Reω≥0. This implies that the eigenvalues of the characteristic polynomial associated with the linearized equations have negative real part, that is, Reω<0. Then, the coexistence equilibrium Ec is locally asymptotically stable.Considering the limiting system(61)dXidt=β1(N1-Xi-Xr)Yi-m1Xi,dXrdt=kβ1(N1-Xi-Xr)Yr+ρσXi-m2Xr,dYidt=β2(N2-Yi-Yr)Xi-m3Yi,dYrdt=kβ2(N2-Yi-Yr)Xr-m4Yr,
where N1=Λ1/μ1 and N2=Λ2/μ2. In this way, let Z0=(Z1,Z2,Z3,Z4), Zi∈C. Substituting a solution of the form (60) into the linearized system (61) of the coexistence equilibrium Ec, we obtain the following linear equations:
(62)ωZ1=-(β1Yic+m1)Z1-β1YicZ2ωZ1=+β1(N1-Xic-Xrc)Z3,ωZ2=(ρσ-kβ1Yrc)Z1-(kβ1Yrc+m2)Z2ωZ2=+kβ1(N1-Xic-Xrc)Z4,ωZ3=β2(N2-Yic-Yrc)Z1-(β2Xic+m3)-β2XicZ4,ωZ4=kβ2(N2-Yic-Yrc)Z2-kβ2XrcZ3ωZ4=-(kβ2Xrc+m4)Z4,
which is equivalent to the system
(63)[1+1m1(ω+β1Yic)]Z1=-β1Yicm1Z2+β1(N1-Xic-Xrc)m1Z3,[1+1m2(ω+kβ1Yrc)]Z2=(ρσ-kβ1Yrc)m2Z1+kβ1(N1-Xic-Xrc)m2Z4,[1+1m3(ω+β2Xic)]Z3=β2(N2-Yic-Yrc)m3Z1-β2Xicm3Z4,[1+1m4(ω+kβ2Xrc)]Z4=kβ2(N2-Yic-Yrc)m4Z2-kβ2Xrcm4Z3.
Moving all the negative terms to the left-hand side, after some manipulations we obtain the system
(64)[1+F1(ω)]Z1+[1+F2(ω)]Z2=(HZ)1+(HZ)2,[1+F3(ω)]Z3+[1+F4(ω)]Z4=(HZ)3+(HZ)4,
where
(65)F1(ω)=ω+β1Yicm1+kβ1Yrcm2,F2(ω)=ω+kβ1Yrcm2+β1Yicm1,F3(ω)=ω+β2Xicm3+kβ2Xrcm4,F4(ω)=ω+kβ2Xrcm4+β2Xicm3,H=(00β1(N1-Xic-Xrc)m10ρσm200kβ1(N1-Xic-Xrc)β2(N2-Yic-Yrc)m3000000kβ2(N2-Yic-Yrc)m4).
Note that the matrix H has nonnegative entries, and Ec=(Xic,Xrc,Yic,Yrc) satisfies
(66)Ec=HEc.To show thatReω<0, we distinguish two cases: ω=0 and ω≠0. In the first case, (62) is a homogeneous linear system. Through calculations, we have the determinant of (62) is
(67)Δ=|-(β1Yic+m1)-β1Yicβ1Xsc0ρσ-kβ1Yrc-(kβ1Yrc+m2)0kβ1Xscβ2Ysc0-(β2Xic+m3)-β2Xic0kβ2Ysc-kβ2Xrc-(kβ2Xrc+m4)|=-(β1Yic+m1)A1+β1YicA2+β1XscA3,
where
(68)A1=m1m3k2β2Xic+m1m32k2-m4kβ1β2XicYrcA1=-m3k2β1β2XrcYrc-m3m4kβ1Yrc-m2m4β2XicA1=-m2m3kβ2Xrc-m2m3m4,A2=ρσ(m3kβ2Xrc+m4β2Xic+m3m4)-m4kβ1β2XicYrcA2=-m3k2β1β2XrcYrc-m3m4kβ1Yrc-m1m3k2β2Xrc,A3=k2β1β22XscYsc2+ρσkβ22XicYsc-k2β1β22XicYrcYscA3=-β2Ysc(k2β1β2XrcYrc+m4kβ1YrcA3=aaaaaaa+m2kβ2Xrc+m2m4).
Then, we can obtain
(69)Δ=ρσ(m3kβ1β2XrcYic+m4β1β2XicYicaaaa+m3m4β1Yic+m1m3kβ2Xic)+(β1β2XicYic+m1β2Xic+m3β2Yic)×(m2m4-m1m3k2)+β1β2XrcYicm3k(m2-m1k)+β1β2XicYrcm1k(m4-m3k).
Note that the coexistence equilibrium Ec exists under the condition R0i>R0r, that is, m2m4>m1m3k2. Then Δ>0 if m2-m1k≥0 and m4-m3k≥0. Hence, for the case that ω=0, (62) has only the trivial solution if m2-m1k≥0 and m4-m3k≥0. This implies that ω≠0.Now assume thatω≠0 and Reω≥0. Let F(ω)=min{|1+Fi(ω)|,i=1,2,3,4}. It is easy to prove that |1+Fi(ω)|>1 for all i, and therefore F(ω)>1. In the other hand, since the coordinates of Ec are positive, if Z¯ is any solution of (64), then there exists a minimal positive real number s such that
(70)|Z¯|≤sEc,
where |Z¯|=(|Z1|,|Z2|,|Z3|,|Z4|), and |·| is the norm in C. Note that s is also the minimal positive real number such that
(71)|Z1|+|Z2|≤s(Xic+Xrc),|Z3|+|Z4|≤s(Yic+Yrc).
Taking norms on both side of (64) and using (66), (71), and the fact that H is nonnegative, we obtain
(72)F(ω)(|Z1|+|Z2|)≤H(|Z|)1+H(|Z|)2≤s(Xic+Xrc),
that is,
(73)|Z1|+|Z2|≤sF(ω)(Xic+Xrc)<s(Xic+Xrc),
which contradicts the minimality of s. Hence, Reω<0.Summering above analysis, we have the following result.Theorem 10.
Whenρ≠0, the coexistence equilibrium Ec of the system (1) is locally asymptotically stable if k≤min{m2/m1,m4/m3} and R0i>1.The existence and stability of equilibria can be summered in Table1.Table 1
Existence and stability of equilibria. The sign “∃  !” means the existence and uniqueness. LAS means locally asymptotically stable, and GAS means globally asymptotically stable.
ρ
=
0
ρ
≠
0
Existence
Stability
Existence
Stability
E
0
∃  !
IfR0<1, LAS; If R0≤1, GAS
∃  !
IfR0<1, LAS; If R0≤1, GAS
E
0
i
IfR0i>1, ∃  !
IfR0i>1 and R0i>R0r, LAS; If R0i>1≥R0r, GAS
Not exist
/
E
0
r
IfR0r>1, ∃  !
IfR0r>1 and R0r>R0i, LAS; If R0r>1≥R0i, GAS
IfR0r>1, ∃  !
IfR0r>1 and R0r>R0i, LAS; If R0r>1≥R0i, GAS
E
c
IfR0i=R0r>1, a line
/
IfR0i>1 and R0i>R0r, ∃  !
IfR0i>1 and k≤min{m2/m1,m4/m3}, LAS


## 6. Discussion

In this paper we established a new schistosomiasis model. In contrast to previous schistosomiasis models with drug resistance, the model established in this study consider many aspects. First, snail is considered as a variable in the model since the resistance of schistosoma to PZQ can be expressed in snails [14–16]. Second, previous models considered that resistance was caused by drug treatment, while a large number of the literature show that some cases are due to inheritance [1, 3, 4]. Therefore, in this study, we consider the reasons for resistance are drug treatment and genetic. Last, in previous models a resistance level (θ) was used to discount the treatment rate. But the resistance level could not be measured, and its value was assumed. This model consider what percentage of infected human after treatment will recover and what percentage will emerge drug resistance. In medicine, from the occurrence of cases the value of this ratio can be identified. For example, this ratio was given in [27, 28] (1-ρ=0.28~0.609). Therefore, it is easy to operate. In addition, we separate the discussion of the case that treatment will cause resistance (ρ≠0) and treatment will not cause resistance (ρ=0).The reproductive numberR0i and R0r of the sensitive and resistant strains are given, respectively. It is easy to see that R0i is a decreasing function of σ, and R0r is an increasing function of k. If the basic reproductive number R0=max{R0i,R0r} of the model (1) is less than 1, one can prove the stability of the disease free equilibrium. This means that the spread of schistosomiasis can be effectively controlled. When the basic reproductive number R0 is greater than 1, we first consider the case that ρ=0. If R0i=R0r>1, the two strains can coexist. There is a line of coexistence equilibria in this case. The infected human will evolve to one of them with higher reproductive number (see Table 1).When drug treatment can not cause resistance (ρ=0), that is, the new resistant strain is due to the inheritance of resistance, the sensitive strain will dominate if the treatment rate is smaller and the impact of the inheritance of resistance is smaller such that R0i>1≥R0r. This result accords to the results of previous models. On the other hand, if the treatment rate and the impact of the inheritance of resistance are both larger such that R0r>1≥R0i, the resistant strain will dominate. This shows that although the assumption is that drug treatment does not result in the emergence of drug resistance, once the treatment rate is greater than a value, and the impact of genetic resistance is larger, there will still be the emergence of resistant strain, and the resistant strain is dominant. This further shows that genetic resistance has a great impact on the system.When drug treatment can cause resistance (ρ≠0), we can show that the sensitive strain either does not appear or coexist with the resistant strain under certain condition. If R0r>1, there is only resistant strain. It can be seen that, regardless of whether drug treatment leads to the emergence of resistance, once the impact of genetic resistance is larger, resistant strain will be dominant, which is detrimental to the control of schistosomiasis.Finally, from the formula of the coexistence equilibriumEc, it is easy to see that the value of the resistant strain is increased with the value of ρ. This means once the proportion of human with drug-resistant strain produced by drug treatment is larger, the number of human and snails with resistant strain is larger. This is not a good result for drug treatment with praziquantel.Hence, for poor treatment, there are two possible reasons: drug therapy and genetic. An important priority in developing new control strategies is to search new drug targets, in combination with selection of genetic methods such as that viable vaccine candidates. And there is already a need for alternative drugs to treat PZQ-resistant schistosomiasis, such as already exists in northern Senegal [10].



---

*Source: 945767-2013-03-06.xml*