www.biomathforum.org/biomath/index.php/biomath ORIGINAL ARTICLE Stability Analysis of a Schistosomiasis Transmission Model with Control Strategies Mouhamadou Diaby Inria, Université de Lorraine, CNRS. ISGMP Bat. A, Ile du Saulcy, 57045 Metz Cedex 01, France. UMI-IRD-209 UMMISCO and LANI Université Gaston Berger, Saint-Louis, Sénégal. diabloss84@yahoo.fr Received: 24 November 2014, accepted: 16 April 2015, published: 20 May 2015 Abstract—We have established and rigorously analyzed a new mathematical model that describes the dynamics schistosomiasis infection. This model incorporates several realistic features including density-dependent births rate of snails and reduced fecundity in snail hosts. Our qualitative analysis of the deterministic model is made with respect to the stability of the disease free equilibrium and the unique endemic equilibrium. Some biological consequences and control strategies are discussed. We have derived the basic reproduction number above which the infection will be controlled under certain levels. We have shown that the disease free equilibrium is globally asymptotically stable when the basic reproduction number R0 is less than one. We have proved the existence and global asymptotic stability of an endemic steady state when R0 > 1. This mathematical analysis of the model gives in- sight about the effects of the reduced fecundity and intermediate host density-dependent birth rate. Finally, numerical simulations are performed to illustrate the main results. Keywords-Epidemic models; Nonlinear dynamical systems; Monotone systems; Global stability; Repro- duction number; Schistosomiasis. I. INTRODUCTION Mathematical modeling of the spread of infec- tious diseases is an important instrument in the comprehension of the dynamics of diseases and in decision making processes regarding intervention programs for controlling these diseases in many countries. The high prevalence of infectious diseases as schistosomiasis has prompted the mathematical modelers to deploy multiple infectious disease models in recent years, and different mechanisms have been suggested to explain their occurrence. Schistosomiasis is transmitted by human contact with contaminated fresh water (lakes and ponds, rivers, dams) inhabited by snails carrying the parasite. Many people die from schistosomiasis disease every year. Schistosomiasis is the most deadly NTD (Neglected Tropical Diseases), killing an estimated 280, 000 people each year [20]. It can result in liver damage and anemia, especially among children [22]. It is found mostly in rural areas in tropical and subtropical countries and infects humans and other vertebrates, using snails in most cases as intermediate hosts. Infection takes place when Citation: Mouhamadou Diaby, Stability Analysis of a Schistosomiasis Transmission Model with Control Strategies, Biomath 1 (2015), 1504161, http://dx.doi.org/10.11145/j.biomath.2015.04.161 Page 1 of 13 http://www.biomathforum.org/biomath/index.php/biomath http://dx.doi.org/10.11145/j.biomath.2015.04.161 M. Diaby, Stability Analysis of a Schistosomiasis Transmission Model ... larvae of the parasite released by the intermediate hosts, penetrate the skin of an individual when in contact with infested water. In the body, the larvae develop and cross over into adult schistosomes. These parasites live in the blood vessels, where the females lay their eggs. Some eggs out of the body via the feces or urine and the parasitic life cycle continues. Others are trapped in the tissues of the body, causing an immune reaction in active lesions and organs. A possible method of struggle against schisto- somiasis is to treat water bodies with a mollusci- cide to minimize the number of snail intermediate hosts, thus breaking the cycle of disease transmis- sion [23], [24]. Beside this one, there are other such that biological control, chemical mollusci- cides, chemotherapy, and more permanent methods such as the provision of safe water and sanitary facilities [1], [6]. Control strategies focused on intermediate hosts which are the sites of intense proliferation of this parasite are seen as a priority of the reduction of schistosomiasis transmission [28]. Many researchers studied the dynamics of schis- tosomiasis using systems similar to the one we are interested [2], [7], [13], [14], [15], [16], [21]. That systems depend on the epidemiological for- mulation, but also on the demographic process incorporated into the model under different control strategies. Almost all of the work to date on schistosomiasis transmission has assumed constant immigration in the snail population [25], [36] or linear birth rates [26], [1], [27] not allowing the possibility of the reduced fecundity of the infected population snails. A drawback of the models with birth and death rates proportional to the size of the population are that the population size decreases or increases exponentially, except in the special case where births exactly balance deaths. Moreover, researchers in laboratory reported the influence of a trematode on life history traits of adult Lymnaea elodes snails. It has been displayed re- duced fecundity relative to uninfected snails [29]. To our knowledge, due to various considerations of the factors related to the transmission of the disease, impaired fertility and density-dependent birth rate are not taken account simultaneously with chemotherapy and chemical molluscicides as control strategies. We shall introduce and analyze a compartmental model, which considers a host (human), interme- diate host (snail) as well as free-living cercariae and miracidia and their interaction. The model focuses essentially on the aquatic life stages of the parasite. A full mathematical analysis of the model is derived. The aim focus of this paper is to study the effect of chemotherapy combined with chemical molluscicide as control strategies on the dynamics of the model with density-dependent birth rate and impaired fertility in snail host. The paper is organized as follows. In Section II we start by defining the mathematical framework we use and focus on the different processes of transmission that might be appropriate to under- stand schistosomiasis. The main results are devel- oped in Section III, including the determination of R0 the basic reproductive number of the model and the analysis of its stability properties. We show that the disease-free equilibrium is global asymptotic stability when the basic reproductive number is less than one. When the basic reproduc- tive number is larger than one, we prove the global asymptotic stability of a unique endemic equilib- rium by using some properties of K-monotone systems (see [17]). Section IV is devoted to the numerical analysis and control strategies. Finally, in Section V, concluding remarks close the paper. II. MODEL DERIVATION We build an evolutionary outcomes model of interactions between a complex life-cycle parasite schistosoma and its hosts (humans and snails). The parasite populations at the free-living stages are modeled explicitly through miracidia and cercarie. The model consists of a system of ordinary differ- ential equations. Furthermore, it is admitted that the infected snails did not recover from schistosomiasis that their lifetimes are short compared to that of hu- mans. The model sub-divides the total human Biomath 1 (2015), 1504161, http://dx.doi.org/10.11145/j.biomath.2015.04.161 Page 2 of 13 http://dx.doi.org/10.11145/j.biomath.2015.04.161 M. Diaby, Stability Analysis of a Schistosomiasis Transmission Model ... population at time t, denoted by Nh(t), into the following sub-populations of susceptible individ- uals (Sh(t)) and individuals with schistosomiasis symptoms (Ih(t)). So that Nh(t) = Sh(t) + Ih(t). The total snails population at time t, denoted by Ns(t), is sub-divided into susceptible snails (Ss(t)) and infectious snails (Is(t)). Thus, Ns(t) = Ss(t) + Is(t). Let Wm and Wp be the population of miracidia and cercariae, respectively. This model assumes a constant per capita rate of exposure between hosts and sensitive parasites, where exposure and sus- ceptibility to infection are regrouped in the trans- mission coefficient. We considered the mass action transmission model. We follow of the available model for schistosomiasis [30] by incorporating the human interaction. Susceptible host snails increase through density- dependent births with maximum rate bs and com- petitive intensity c. Assume that individuals are born uninfected. Infected hosts suffer from re- duced fecundity (0 ≤ ρ < 1). Susceptible snails die at background death rate ds, due to parasite natural death rate ds� and the elimination rates θs of snails and become infected at the per capita infection rate βs, through contact with infectious miracidia. Infected snails die at an elevated death rate, ds, due to parasite natural death rate ds� and the elimination rates θs of snails at which is added the rate α due to parasite virulence. The susceptible human population is increased by the recruitment of individuals in the population (assumed susceptible), at a rate Λ. The population of individuals is further decreased by natural death (at a rate dh). It is also assumed that infected individuals have additional host mortality during the given short time considered at the rate µ. Susceptible humans become infected only through contact with free-living pathogen cercariae in in- fested water at the per capita infection rate βh and recover at the per capita rate η. Free-living miracidia are introduced into the aquatic environ- ment at a per capita rate k, but they are depleted during the infection process at the per capita rate δ and die naturally if they do not find snails to infect at the rate dm. Here, depletion of parasites through transmission depends on total host density. Infected snails will then free up second form of larvae called cercariae at a rate γ to be able to infect humans. Some cercariae also die at an elevated death rate dc due to the natural death rate dc� and cercariae elimination θc by the chemical molluscicide. The time evolution of the different populations is governed by the following system of equations:   dSs dt = bs (Ss + ρIs)(1 − c (Ss + Is)) − ds︷ ︸︸ ︷ (ds� + θs) Ss −βs Ss Wm, dIs dt = βs Ss Wm − (ds� + θs + α) Is, dWm dt = k Ih −δ (Ss + Is) Wm −dm Wm, dWc dt = γ Is− dc︷ ︸︸ ︷ (dc� + θc) Wc, dSh dt = Λ −βh Wc Sh −dh Sh + η Ih, dIh dt = βh Wc Sh − (dh + µ + η) Ih. (1) For convenience of simplicity, we denote ds = ds� + θs, dc = dc� + θs. The dynamics of the total snails population (Ns = Ss + Is) is governed by dNs dt = bs (Ss + ρIs) (1 − cNs) −ds Ns −αIs ≤ bs Ns − bs cN2s −ds Ns ≤ ( bs −ds bs c −Ns ) bs cNs. The dynamics of the total humans population (Nh = Sh + Ih) is governed by dNh dt = Λ −dh Nh −µIh ≤ Λ dh . Biomath 1 (2015), 1504161, http://dx.doi.org/10.11145/j.biomath.2015.04.161 Page 3 of 13 http://dx.doi.org/10.11145/j.biomath.2015.04.161 M. Diaby, Stability Analysis of a Schistosomiasis Transmission Model ... Thus the region D = {(Ss,Is,Wm,Wc,Sh,Ih) ∈ R6+ : Ns ≤ bs −ds cbs , Wm ≥ 0, Wc ≥ 0, Nh ≤ Λ dh } is a positively invariant set of system (1). Furthermore, the model (1) is well-posed epi- demiologically and we will consider dynamic be- havior of model (1) on D. III. MAIN RESULTS In this section, we firstly derive the basic re- production number for system (1) by the method of next generation matrix formulated in [9], [5]. The basic reproduction number R0 is used to assess the stability of the disease free equilibrium and the endemic equilibrium. Then we derive the asymptotic behavior of the system (1) by its limiting system. Discussions of the relation of the dynamics between (1) and its limiting system can be found in [32], [33], [34]. In particular, there holds the following significant result : Lemma III.1. Consider the following two systems dx dt = f(t,x), dy dt = g(y), where x,y ∈ Rn, f and g are continuous, satisfy a local Lipschitz condition in any compact set Ω ∈ Rn, and f(t,x) → g(x) as t → ∞, so that the second system is the limit system for the first system. Let Φ(t,t0,x0 and Φ(t,x) be solutions of theses systems, respectively. Suppose that p ∈ Ω is a locally asymptotically stable equilibrium of the limit system and its attractive region is W(p) = {y ∈ Ω|Φ(t,t0,y) → p,t → +∞}. Let WΦ be the omega limit set of Φ(t,t0,x0). If WΦ ∩W(p) 6= ∅, then limt→+∞Φ(t,t0,x0) = p. This limiting system is the key to understanding the global dynamics of system (1). By a method due to Castillo-Chavez et al [31], the stability of the disease free equilibrium of the proposed limit- ing system is discussed. Moreover, the stability of positive equilibrium of a proposed limiting system of model (1) is analyzed theoretically thanks to monotone dynamical system theory [11]. A. The reproduction number The reproduction number is the expected num- ber of secondary cases produced in a completely susceptible population by a typical infective in- dividual. It is easy to see that system (1) ad- mits always a disease-free equilibrium, E0 =( bs −ds bs c , 0, 0, 0, Λ µ ) . Let x = (Is,Ih,Wm,Wp,Ss,Sh)T . Then sys- tem (1) can be written as x′ = F(x) −V(x), where F =   βs Ss Wm βh Sh Wc 0 0 0 0   , and V =   (ds + α) Is (µ + dh + η) Ih −k Ih + δ (Ss + Is) Wm + dm Wm −γ Is + dc Wc −bs (Ss −ρIs)(1 − c (Ss + Is)) +ds Ss + βs Ss Wm −Λ + βh Wc Sh + µIh + η Ih   . We can get F =   0 0 βs (bs −ds) bs c 0 0 0 0 βh Λ µ 0 0 0 0 0 0 0 0   and V =   α + ds 0 0 0 0 η + µ 0 0 0 −k δ (bs −ds) bs c + dm 0 −γ 0 0 dc   F V −1 is the next generation matrix for model (1). It then follows that the spectral radius Biomath 1 (2015), 1504161, http://dx.doi.org/10.11145/j.biomath.2015.04.161 Page 4 of 13 http://dx.doi.org/10.11145/j.biomath.2015.04.161 M. Diaby, Stability Analysis of a Schistosomiasis Transmission Model ... of matrix F V −1 is ρ(F V −1). According to Theorem 2 in [18], the basic reproduction number of model (1) is R0 = ρ(F V −1) = √ β γ k Λ (bs −ds) βh dh dc (ds + α) (η + µ + dh) (bs (cdm + δ) −δds) B. The local stability of the DFE Before proceeding with the mathematical anal- ysis of the model (1), the following simplification would be made. To make the mathematical analy- sis more tractable, we will consider the case where the the disease-induced death and the parasite virulence parameters in the model will be set to zero. That is, from now on, it is assumed that α = µ = 0 (it should be mentioned that such assumption may not be appropriate in modeling the transmission dynamics in hardest hit areas). However, this assumption will be relaxed in the numerical simulations section. Without loss of generality (see [33], [32]), we assume that our population of humans has reached its limiting value, i.e, N∗h ≡ Λ dh (with α = µ = 0). By eliminating the equation for dSh dt , we get from (1) the equivalent limiting system :  dSs dt = bs (Ss + ρIs)(1 − c (Ss + Is)) −ds Ss −βs Ss Wm, dIs dt = βs Ss Wm −ds Is, dWm dt = k Ih −δ (Ss + Is) Wm −dm Wm, dWc dt = γ Is −dc Wc, dIh dt = βh Wc (N ∗ h − Ih) − (dh + η) Ih. (2) defined on D0 = {(Ss,Is,Wm,Wc,Ih) ∈ R5+ : Ns ≤ bs −ds cbs , Wm ≥ 0, Wc ≥ 0, Ih ≤ Λ dh } . We now study the local behavior of the disease free equilibrium for system (2), which also gives the analogous behavior for system (1). Proposition III.1. The DFE for limiting system (2) E0 is LAS (locally asymptotically stable) if R0 < 1 with ρ ≤ 1 −ds/bs and is unstable if R0 > 1. Proof: The Jacobian matrix of (2) at E0 is J0 = ( J11 J12 J21 J22 ) where, J11 = ( ds − bs (ρ + 1)ds − bs 0 −ds ) , J12 =   βs (ds − bs) bs c 0 0 βs (bs −ds) bs c 0 0   ,J21 =   0 00 γ 0 0   , J22 =   δ ds − bs (δ + cdm) bs c 0 k 0 −dc 0 0 Nhβh −η −dh   . Let us do the following transformation: Js0 = T.J0.T where, T =   1 0 0 0 0 0 −1 0 0 0 0 0 −1 0 0 0 0 0 −1 0 0 0 0 0 −1   . Then, Js0 and J0 are similar. Js0 is a Metzler matrix and we can write Js0 = F + V with F =   0 bs −ds βs (bs −ds) bs c 0 0 0 0 βs (bs −ds) bs c 0 0 0 0 0 0 k 0 γ 0 0 0 0 0 0 Nh βh 0   , Biomath 1 (2015), 1504161, http://dx.doi.org/10.11145/j.biomath.2015.04.161 Page 5 of 13 http://dx.doi.org/10.11145/j.biomath.2015.04.161 M. Diaby, Stability Analysis of a Schistosomiasis Transmission Model ... V=   ds − bs−ρds 0 0 0 0 −ds 0 0 0 0 0 δ ds−bs(δ+cdm) bs c 0 0 0 0 0 −dc 0 0 0 0 0 −η −dh   . We have F > 0 and V is Metzler stable, see [3], [12], [10], [4]. Thanks to Varga’s Theorem in [19]: s(J0) ≤ 0 iff ρ(−F V −1) ≤ 1. Since ρ(−F V −1) = R40, we then deduce that E0 is LAS if R0 < 1 and is unstable if R0 > 1. C. Global stability of the DFE The object of this subsection is to prove the GAS result using the second generation approach given in Castillo-Chavez et al [31]. Herein stated to be self-contained. Theorem III.1. If the system (2) can be written in the form dX dt = F(x,Z), dZ dt = G(X,Z),G(x,Z) = 0, where the components of X ∈ Rm denotes the number of uninfected individuals and the components of Z ∈ Rn the number of infected individuals. let E0 = (x∗, 0) be the disease free equilibrium of the system. Assume that • For dX dt = F(X, 0), X∗ is globally asymp- totically stable (GAS), • G(X,Z) = AZ − Ĝ(X,Z), Ĝ(X,Z) ≥ 0 for (X,Z) ∈ D, where A = DZ G(X∗, 0) is an M-Matrix (the off diagonal of A are nonnegative) and D the region where the model makes biological sense. Then the fixed E0 = (x∗, 0) is a globally asymptotically stable equilibrium of model system (2) provided that R0 < 1. Applying Theorem III.1 to system (2) gives Ĝ(X,Z) =   Wm βs (S ∗ −Ss) Wm δ (Is + Ss −S∗s ) 0 Wc Ih βh   . Since, S ≤ S∗ ≡ N∗h thus Ĝ(X,Z) ≥ 0, and by Theorem III.1, E0 is GAS. We summarize the result in the following theo- rem: Theorem III.2. If R0 < 1 then the DFE is GAS. D. Local Stability of the Endemic Equilibrium We proceed to investigate the local stability of the model. We use the following limiting system (3) to obtain the information for the whole system. Since the population size, Ns(t) (with α = µ = 0) also satisfies Ns → N∗s := bs −ds cbs as t →∞. Using results from Lemma III.1 , we can get analytical results by considering the limiting system of (3) in which the total population of snails and humans both have reached the limiting states. Then, we obtain the reduced limiting dynamical system:  dIs dt = βs (Ns − Is) Wm −ds Is, dWm dt = k Ih −δ N∗s Wm −dm Wm, dWc dt = γ Is −dc Wc, dIh dt = βh Wc (N ∗ h − Ih) − (dh + η) Ih. (3) The dynamics of system (1) can be focused on this restricted region D1 = {(Is,Wm,Wc,Ih) ∈ R4+ : Is ≤ bs −ds cbs , Wm ≥ 0, Wc ≥ 0, Ih ≤ Λ dh } Clearly, the system (1) is asymptotic to the limiting system (3). Thus, the dynamics of system (1) are qualitatively equivalent to the dynamics of its limiting system. The variation matrix of system (3) at E∗ an equilibrium point is J(E∗) =   J11(E∗) J12(E∗) J21(E ∗) J22(E ∗)   , where, J11(E ∗) = ( −W∗m βs −ds (N∗s − I∗s )βs 0 −N∗s δ −dm ) , Biomath 1 (2015), 1504161, http://dx.doi.org/10.11145/j.biomath.2015.04.161 Page 6 of 13 http://dx.doi.org/10.11145/j.biomath.2015.04.161 M. Diaby, Stability Analysis of a Schistosomiasis Transmission Model ... J12(E ∗) = ( 0 0 0 k ) , J21(E ∗) = ( γ 0 0 0 ) , J22(E ∗) = ( −dc 0 (N∗h − I ∗ h)βh −η −W ∗ c βh −dh ) Theorem III.3. If R0 > 1, then the positive endemic equilibrium E∗ of the limiting system (3) is locally asymptotically stable. Before proving Theorem III.3, we state the following useful lemma : Lemma III.2. (Kamkang [35], Proposition 3.1) Let M be a square Metzler matrix written in block form M = ( A B C D ) , with A and D square matrices. M is Metzler stable if only if matrices A and D −C A−1 B are metzler stable. Proof of Theorem III.3: J(E∗) is a Metzler matrix and we can write : J(E∗) = ( A B C D ) , where A=  −W∗m βs −ds (N∗s − I∗s )βs 0 −N∗s δ −dm   , D=   −dc 0 (N∗h − I ∗ h)βh −η −W ∗ c βh −dh   . B= ( 0 0 0 k ) ,C= ( γ 0 0 0 ) . Clearly, A is a stable Metzler matrix. Then, after some computations, we get D −C A−1 B =   M11 M12 M21 M22   . where, M11 = −dc, M12 = k ( N∗s − I∗h)β γ W∗m N ∗ s β δ + N ∗ s ds δ + W ∗ m β dm + ds dm M21 = (N ∗ h −I ∗ h)βh, M22 = −η−W ∗ c βh−dh. D −C A−1 B is a stable Metzler matrix iff χ:=β γ k βh(I ∗ h −N ∗ s )(N ∗ h − I ∗ h) + dc (dm + δ N ∗ s ) (ds + β W ∗ m) . (W∗c βh + η + dh) > 0. The endemic equilibrium satisfies  βs (Is −Ns) Wm = −ds Is, Wc = γ dc Is. k Ih = (δ Ns + dm) Wm. (4) Hence, χ=−γ k βh ds Wm Is(N ∗ h − I ∗ h) + dc (dm + δ N ∗ s ) . (ds + βs W ∗ m) (W ∗ c βh + η + dh) >−γ k βh ds Wm Is(N ∗ h − I ∗ h) +dc (dm + δ N ∗ s ) ds βh γ Is dc > ds βh γ Is [ − k Wm (N∗h − I ∗ h) + (dm + δ N ∗ s ) ] >ds βh γ Is [ − k Wm I∗h + (dm + δ N ∗ s ) ] >ds βh γ Is [ − Wm (δ Ns + dm) Wm I∗h + (dm + δ N ∗ s ) ] > 0, which implies J(E∗) is a Metzler stable matrix. Thus the unique endemic equilibrium is locally assymptoticaly stable. E. Global Stability of the Endemic Equilibrium In this section we will establish the global stability of the unique endemic equilibrium point when R0 > 1. We shall use the properties of K- monotone systems for the analysis of our system (see [17]). Biomath 1 (2015), 1504161, http://dx.doi.org/10.11145/j.biomath.2015.04.161 Page 7 of 13 http://dx.doi.org/10.11145/j.biomath.2015.04.161 M. Diaby, Stability Analysis of a Schistosomiasis Transmission Model ... Theorem III.4. If R0 > 1, then the positive en- demic equilibrium state E∗ of the limiting system (3) is globally asymptotically stable in the interior of the set D1. Proof: It is useful to be able to test mono- tonicity directly in terms of vector fields. We recall here a wider criterion to define monotonicity. A system X = f(X) is said to be monotone in the set K if there exists a diagonal matrix H = diag(�1,�2, · · · ,�n) where each �i is 1 or −1 such that the matrix product H Df(X) H has only non positive values outside the diagonal for any X (see Smith [17], Lemma 2.1). In the above statement, Df(X) is the Jacobian matrix of f and K is a convex set. The Jacobian of the system (3), J(E∗), is a Metzler matrix and is irreducible, which implies the strong monotonicity of the system under a usual order defined by the orthant K = {(Is,Wm,Wc,Ih) ∈ R4+ : Is, Wm, Wc, Ih > 0} The significant result of convergence proved by Hirsh in the late 1980s can be described in our framework as follows : given an autonomous system that is strongly monotone with respect to some proper cone and assuming that there is a unique equilibrium in an open set of points with compact orbit closures, every initial condition with bounded solution converges to the unique equilibrium; (see Hirsh [11], theorem 10.3). Thanks Hirsch’s theorem and the fact that we have only one endemic equilibrium E∗ in D̊1 which is locally asymptotically stable when R0 > 1 we can conclude that E∗ is globally asymptoti- cally stable in D̊1 when R0 > 1. IV. NUMERICAL ANALYSIS AND CONTROL STRATEGIES In this section, we present some numerical simulation results to confirm our analytical pre- dictions on the global dynamics of the schisto- somiasis models. The two selected control strate- gies identified in this study against the spread of schistosomiasis are snails reduction strategies and chemotherapy treatments. Snail reduction strate- gies include the elimination of snails as well as the elimination of free-living cercariae using appro- priate biological agents. The strategy of applying chemotherapy treatments requires an effective drug delivery as praziquantel to infected humans. In order to examine the effects of these two classes of anti-schistosomiasis strategies, the model (1) is simulated using the set of parameter values given in Table I. The effect of chemical molluscicide is incorpo- rated in our model by increasing the death rate parameter of the snails and free-living cercariae at the rate θs and θc, respectively. In other words, a higher values in θs and θc serve to reduce the number of snails and free-living cercariae, and thus can be used to investigate the effects of chemical molluscicide on the epidemic. Similarly, the effects of the drug administration method are modelled by increasing a recover rate η to the schistosomiasis. The basic reproduction number can be applied to measuring the control efforts needed to reduce or eliminate an infection. Furthermore, since R0 is a decreasing function of θs, θc and η, the use of any preventive strategy that can increase θs, θc or η results in a reduction of schistosomiasis infections. We, consequently, seek to compare the effects of chemical molluscicide and drug administration on the control of the spread of schistosomiasis in humans. We first consider the original no-control model by setting η = θs = θc = 0 in system (1). We set the population size of initial conditions as Ss(0) = 104; Is(0) = 5 × 104; Wm(0) = 5×103; Wc(0) = 9×103; Sh(0) = 2×103; Ih(0) = 550, and take the values of the parameters from table I. It should be mentioned that some of these parameters (e.g. θs, θc, η) are estimated in the con- venience use. We pick the following conclusions for the parameter values used. The results in Fig. 1 show a marked increase in the number of asymptotically infected snails and infected humans at steady state with no control strategies (θs ≈ θc ≈ η ≈ 0). With the same configurations, we now add control measures and simulate the control model (1). The simulation Biomath 1 (2015), 1504161, http://dx.doi.org/10.11145/j.biomath.2015.04.161 Page 8 of 13 http://dx.doi.org/10.11145/j.biomath.2015.04.161 M. Diaby, Stability Analysis of a Schistosomiasis Transmission Model ... TABLE I PARAMETER VALUES ESTIMED Parameter Definition Default value Literature source bs Maximum birth rathe of hosts 0.06 day−1 [30] ρ Relative fecundity of infected hosts 0.75 day−1 −− c Strenght of density dependence on host birth rate 0.025 L day−1 −− dh Nartural mortality rate of humans 0.0000384 day −1 −− ds Natural mortality rate of snails 0.003 day−1 −− µ Parasite dependent mortality of humans 0.0039 day−1 −− dm Loss rate of free-living miracidia 2.5 day−1 −− dc Loss rate of free-living cerscariae 2 day −1 −− βs Rate at which the miracidia successfully infects a snail 0.615 L mir−1 day−1 −− βh Rate at which susceptible humans get infected 0.406 L mir −1 day−1 [36] α Parasite virulence on survival 0.007 day−1 [30] γ Per capita production rate of cercariae by infected snails 100 cerc host−1 day−1 [36] k Per capita production rate of miracidia by infected humans 0.00232 day−1 [36] η Treatment rate of infected human 0.03 day−1 Estimated δ Per capita depletion rate of miracidia 0.0039 day−1 [30] θs, θc elimination rates of snails and cercariae, respectively 0.05 day−1 Estimated Λ Recruitment rate for humans 8000 day−1 [36] shown in Figure 2 demonstrate the effectiveness of the control strategies adopted. We observe that an epidemic outbreak occurs for some period of time. For instance, if human recover at the rate η = 0.9 and snails be eliminated at the rate η = 0.05, then the number of asymptomatically infected humans tend to zeros after a period of time t ≈ 1050 days. Mention should be made here of the fact that, in the last scenario, we have values for each of the following parameters η, θs, θc, so that R0 < 1. Fig. 3 and 4 also show that the efficacy of snail elimination can be important for such strategies to make a meaningful impact in combatting schistosomiasis in humans. V. SUMMARY AND CONCLUSIONS In this paper, we have presented a stability analysis of a schistosomiasis infection model that explicitly includes density dependent births rate and impaired fecundity in the snail host pop- ulation. The model also captures the effect of two control strategies: chemotherapy and snail elimination by molluscicides. Six sub population sizes were considered: human host susceptible and infected, snail intermediate host susceptible and infected, free-living miracidia and cercariae. It is shown that the combination of chemotherapy and snail elimination can be effective control strategy that is there will not be an epidemic. Mathematical properties of the model are ana- lyzed and used to reduce the dimension of the sys- tem under consideration. The reproductive number R0 is then analytically and explicitly computed. We proved that the disease-free steady state E0 is globally asymptotically stable if R0 ≤ 1. We have also established the global asymptotic stability of the endemic equilibrium E∗ when it exists i.e., when R0 > 1 using some properties of monotone systems. Biomath 1 (2015), 1504161, http://dx.doi.org/10.11145/j.biomath.2015.04.161 Page 9 of 13 http://dx.doi.org/10.11145/j.biomath.2015.04.161 M. Diaby, Stability Analysis of a Schistosomiasis Transmission Model ... 0 200 400 600 800 1000 1200 1400 1600 1800 2000 0 0.5 1 1.5 2 2.5 3 3.5 x 10 6 time In fe c te d s n a il s (a) 0 200 400 600 800 1000 1200 1400 1600 1800 2000 0 0.5 1 1.5 2 2.5 x 10 5 time In fe c te d h u m a n s (b) Fig. 1. Time evolution of the number of infected snails population (a) and humans population (b) without control strategies with parameter values defined in Table I: θc = 0, θs = 0, η = 0. These parameters correspond to R0 > 1. The initial condition is Ss = 10000, Is = 50000, Wm = 5000, Wc = 9000, Sh = 2000, Ih = 550. 0 50 100 150 200 250 300 350 400 450 0 1 2 3 4 5 6 x 10 4 time In fe c te d s n a il s (a) 0 200 400 600 800 1000 1200 1400 1600 0 2 4 6 8 10 12 14 16 18 x 10 4 time In fe c te d h u m a n s (b) Fig. 2. Effect of combined human treatment and snails elimination on snails population (a) and humans population behavior (b) with parameter values defined in Table I: θc = 10, θs = 0.05, η = 0.90. These parameters correspond to R0 < 1. The initial condition is Ss = 10000, Is = 50000, Wm = 5000, Wc = 9000, Sh = 2000, Ih = 550. Biomath 1 (2015), 1504161, http://dx.doi.org/10.11145/j.biomath.2015.04.161 Page 10 of 13 http://dx.doi.org/10.11145/j.biomath.2015.04.161 M. Diaby, Stability Analysis of a Schistosomiasis Transmission Model ... 0 50 100 150 200 250 300 350 400 450 0 0.5 1 1.5 2 2.5 3 x 10 6 time (days) In fe c te d s n a il s (a) 0 200 400 600 800 1000 1200 1400 1600 1800 2000 0 0.5 1 1.5 2 2.5 x 10 5 time (days) In fe c te d h u m a n s (b) Fig. 3. Effect of human treatment only, by drug administra- tion, on snails population (a) and humans population behavior (b) with parameter values defined in Table I: θc = 0, θs = 0, η = 0.90. The initial condition is Ss = 10000, Is = 50000, Wm = 5000, Wc = 9000, Sh = 2000, Ih = 550. 0 50 100 150 200 250 300 350 400 450 0 1 2 3 4 5 6 x 10 4 time (days) In fe c te d s n a il s (a) 0 200 400 600 800 1000 1200 1400 1600 1800 2000 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x 10 5 time (days) In fe c te d h u m a n s (b) Fig. 4. Effect of snails elimination only, by molluscicide, on snails population (a) and humans population behavior (b) with parameter values defined in Table I: θc = 10, θs = 0.05, η = 0. The initial condition is Ss = 10000, Is = 50000, Wm = 5000, Wc = 9000, Sh = 2000, Ih = 550. Biomath 1 (2015), 1504161, http://dx.doi.org/10.11145/j.biomath.2015.04.161 Page 11 of 13 http://dx.doi.org/10.11145/j.biomath.2015.04.161 M. Diaby, Stability Analysis of a Schistosomiasis Transmission Model ... REFERENCES [1] E. J. Allen and H. D. J. Victory. Modelling and simulation of a schistosomiasis infection with biological control. Acta Trop, 87(2):251–267, 2003. [2] R. Anderson and R. May. Helminth infections of humans: mathematical models, population dynamics, and control. Advances in parasitology, 24:1–101, 1985. [3] A. Berman and R. J. Plemmons. Nonnegative matrices in the mathematical sciences. SIAM, 1994. [4] F. Brauer and C. Castillo-Chávez. Mathematical models in population biology and epidemiology, volume 40 of Texts in Applied Mathematics. Springer-Verlag, New York, 2001. [5] C. Castillo-Chavez, Z. Feng, and W. Huang. On the computation of R0 and its role on global stability. In Mathematical approaches for emerging and reemerging infectious diseases: an introduction (Minneapolis, MN, 1999), volume 125 of IMA Vol. Math. Appl., pages 229– 250. Springer, New York, 2002. [6] Z. Chen, L. Zou, D. Shen, W. Zhang, and S. Ruan. Mathematical modelling and control of schistosomiasis in hubei province, china. Acta Tropica, 115(1–2):119– 125, 2010. [7] J. E. Cohen. Mathematical models of schistosomiasis. Annual Review of Ecology and Systematics, 8(1):209– 233, 1977. [8] Diaby, M., Iggidr, A., Sy, M., & Sene, A. (2014). Global analysis of a schistosomiasis infection model with bio- logical control. Applied Mathematics and Computation, 246, 731-742. [9] O. Diekmann, J. Heesterbeek, and J. Metz. On the defini- tion and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations. J. Math. Biol., 28(4):365–382, 1990. [10] O. Diekmann and J. Heesterbeek. Mathematical epi- demiology of infectious diseases. Model building, analy- sis and interpretation. Wiley Series in Mathematical and Computational Biology. Chichester: Wiley. 300 p. , 1999. [11] M. W. Hirsch. Stability and convergence in strongly monotone dynamical systems. J. Reine Angew. Math., 383:1–53, 1988. [12] J. A. Jacquez and C. P. Simon. Qualitative theory of compartmental systems. SIAM Rev., 35(1):43–79, 1993. [13] A. R. Kimbir. A mathematical model of the transmission dynamics of schistosomiasis. J. Nigerian Math. Soc., 16/17:39–63, 1997/98. [14] G. Macdonald. The dynamics of helminth infections, with special reference to schistosomes. Transactions of the Royal Society of Tropical Medicine and Hygiene, 59(5):489–506, 1965. [15] F. A. Milner and R. Zhao. A deterministic model of schistosomiasis with spatial structure. Math. Biosci. Eng., 5(3):505–522, 2008. http://dx.doi.org/10.3934/mbe.2008.5.505 [16] J. P. Pointier and J. Jourdane. Biological control of the snail hosts of schistosomiasis in areas of low transmis- sion: the example of the caribbean area. Acta Tropica, 77(1):53–60, 2000. [17] H. L. Smith. Systems of Ordinary Differential Equations Which Generate an Order Preserving Flow. A Survey of Results. SIAM Review, 30(1):87–113, 1988. [18] P. van den Driessche and J. Watmough. Reproduction numbers and sub-threshold endemic equilibria for com- partmental models of disease transmission. Math. Biosci., 180:29–48, 2002. John A. Jacquez memorial volume. [19] R. Varga. Matrix iterative analysis. 2. printing. (Prentice-Hall Series in Automatic Computation). En- glewood Cliffs, New Jersey:Prentice-Hall, Inc. XIII, 322 p. , 1962. [20] WHO, Schistosomiasis. http://www3.imperial.ac.uk/schisto/schistosomiasis (2014). [21] M. E. J. Woolhouse. On the application of mathematical models of schistosome transmission dynamics. i. natural transmission. Acta Tropica, 49(4):241–270, 1991. [22] C. N. Mascie-Taylor and E. Karim. The burden of chronic disease. Science, 302(5652):1921–1922, 2003. [23] F. McCullough, P. Gayral, J. Duncan, and J. Christie. Molluscicides in schistosomiasis control. Bulletin of the World Health Organization, 58(5):681, 1980. [24] G. Greer, P.-B. Tchounwou, I. Takougang, and A. Monkiedje. Field tests of a village-based mollus- ciciding programme for the control of snail hosts of human schistosomes in cameroon. Tropical Medicine & International Health, 1(3):320–327, 1996. [25] S. Gao, Y. Liu, Y. Luo, and D. Xie. Control problems of a mathematical model for schistosomiasis transmission dynamics. Nonlinear Dynamics, 63(3):503–512, 2011. [26] Z. Feng, J. Curtis, and D. J. Minchella. The influence of drug treatment on the maintenance of schistosome genetic diversity. Journal of Mathematical Biology, 43(1):52–68, 2001. [27] M. Diaby, A. Iggidr, M. Sy, A. Sène, et al. Global analysis of a shistosomiasis infection with biological control. 2012. [28] V. Lardans and C. Dissous. Snail control strategies for reduction of schistosomiasis transmission. Parasitology Today, 14(10):413–417, 1998. [29] R. Wilson and J. Denison. The parasitic castration and gigantism oflymnaea truncatula infected with the larval stages offasciola hepatica. Zeitschrift für Parasitenkunde, 61(2):109–119, 1980. [30] D. J. Civitello and J. R. Rohr. Disentangling the effects of exposure and susceptibility on transmission of the zoonotic parasite schistosoma mansoni. Journal of Animal Ecology, 2014. [31] A.-A. Yakubu and C. Castillo-Chavez. Interplay be- tween local dynamics and dispersal in discrete-time metapopulation models. Journal of theoretical biology, 218(3):273–288, 2002. [32] C. Castillo-Chevez and H. R. Thieme. Asymptotically autonomous epidemic models, volume 94. Mathematical Sciences Institute, Cornell University, 1994. Biomath 1 (2015), 1504161, http://dx.doi.org/10.11145/j.biomath.2015.04.161 Page 12 of 13 http://dx.doi.org/10.3934/mbe.2008.5.505 http://www3.imperial.ac.uk/schisto/schistosomiasis http://dx.doi.org/10.11145/j.biomath.2015.04.161 M. Diaby, Stability Analysis of a Schistosomiasis Transmission Model ... [33] H. R. Thieme. Convergence results and a poincaré- bendixson trichotomy for asymptotically autonomous dif- ferential equations. Journal of mathematical biology, 30(7):755–763, 1992. [34] C. Castillo-Chavez, W. Huang, and J. Li. Competitive exclusion in gonorrhea models and other sexually trans- mitted diseases. SIAM Journal on Applied Mathematics, 56(2):494–508, 1996. [35] J. C. Kamgang and G. Sallet. Global asymptotic stabil- ity for the disease free equilibrium for epidemiological models. Comptes Rendus Mathematique, 341(7):433– 438, 2005. [36] E. T. Chiyaka and W. GARIRA. Mathematical anal- ysis of the transmission dynamics of schistosomiasis in the human-snail hosts. Journal of Biological Systems, 17(03):397–423, 2009. Biomath 1 (2015), 1504161, http://dx.doi.org/10.11145/j.biomath.2015.04.161 Page 13 of 13 http://dx.doi.org/10.11145/j.biomath.2015.04.161 Introduction Model derivation Main results The reproduction number The local stability of the DFE Global stability of the DFE Local Stability of the Endemic Equilibrium Global Stability of the Endemic Equilibrium Numerical analysis and control strategies Summary and conclusions References