Schistosomiasis Model Incorporating Snail Predator as Biological Control Agent

: Schistosomiasis is a parasitic disease caused by the schistosoma worm. A snail can act as the intermediate host for the parasite. Snail-population control is considered to be an effective way to control schistosomiasis spread. In this paper, we discuss the schistosomiasis model incorporating a snail predator as a biological control agent. We prove that the solutions of the model are non-negative and bounded. The existence condition of equilibrium points is investigated. We determine the basic reproduction number when the predator goes to extinction and when the predator survives. The local stability condition of disease-free equilibrium point is proved using linearization, and the Lienard–Chipart and Routh–Hurwitz criteria. We use center-manifold theory to prove the local stability condition of the endemic equilibrium points. Furthermore, we constructed a Lyapunov function to investigate the global stability condition of the disease-free equilibrium points. To support the analytical results, we presented some numerical simulation results. Our ﬁndings suggest that a snail predator as a biological control agent can reduce schistosomiasis prevalence. Moreover, the snail-predator birth rate plays an essential role in controlling schistosomiasis spread.


Introduction
Schistosomiasis is a parasitic disease caused by schistosoma-worm infection [1]. The parasite has a complicated life cycle. Some types of snail can act as intermediate hosts, while humans and mammals such as cows, pigs, and mice can serve as reservoir hosts [2]. Consequently, the parasite can infect mammals as a human substitute to complete its life cycle. This is one of the several factors that make schistosomiasis very difficult to eradicate. Controlling the snail population is the most effective way to control the spread of schistosomiasis. Some researchers recommended the use of molluscicides to manage the snail host population [1,3,4]. However, molluscicides have some negative environmental effects [5]. Another method that can be used to manage the snail population is an intervention with snail predators or competitors [6][7][8]. Sokolow et al. [7] stated that releasing river prawn, which is a snail predator, into a water contact site can reduce the snail population and schistosomiasis transmission.
Mathematical modeling is used to study the dynamics of the spread of the disease. The first mathematical model that is related to schistosomiasis is discussed in [9]. Schistosomiasis models considering parasite density in the environment are discussed in [10][11][12][13][14]. The authors divided the parasite into two compartments, i.e., miracidiae and cercariae, which can infect snails and humans or mammals, respectively. In 2009, Chiyaka et al. [10] investigated the host-parasite dynamics of schistosomiasis. They proposed a schistosomiasis model that consists of six first-order differential equations. They found that control interventions that target transmission to humans are more effective than control strategies that aim the transmission to snails are. Gao et al. [11] developed a schistosomiasis model to investigate some control strategies, i.e., health education, cercaria control, snail control, and drug treatment. They found that killing snails is the most effective method to manage schistosomiasis spread. Nur et al. [13] proposed a schistosomiasis model incorporating health education and molluscicide intervention. They found that the most effective way to control schistosomiasis prevalence is molluscicide intervention. Moreover, schistosomiasis cannot be eradicated if the only intervention is health education [13]. Diaby et al. [15] proposed a schistosomiasis model with biological control, i.e., competitor-resistant snails. They found that competitor-resistant snails can be used to manage the spread of schistosomiasis. Okamoto et al. [16,17] proposed mathematical models of vectorborne diseases. They found that biological control agents that can be used to control vectorborne diseases are natural predators, natural competitors, or parasites of the vector. Schistosomiasis is a waterborne disease. Hence, we study the dynamics of the spread of schistosomiasis when a snail predator is used as a biological control agent of snails. Different from the model discussed in [10][11][12][13][15][16][17], we propose a schistosomiasis model with a biological control agent of snails, namely, a snail predator.

Mathematical Preliminaries
In this section, we present a theorem that can be used to investigate the existence of backward and forward bifurcation. The theorem is very useful, especially in epidemic models. Theorem 1 can be used to investigate the local stability condition of endemic equilibrium [10,11,18]. The proof of Theorem 1 can be found in [19]. [19]) Consider the following general system of ordinary differential equations with a parameter ω.

Theorem 1. (Castillo-Chavez and Song
where 0 is an equilibrium for System (1), such that f (0, ω) ≡ 0 for all ω. Assume The dynamics of System (1) around 0 is totally determined by the signs of A and B.
On the basis of Theorem 1, forward bifurcation occurs at ω = 0 if A < 0 and B > 0. Moreover, backward bifurcation occurs at ω = 0 if A > 0 and B > 0.

Model Formulation
On the basis of the life cycle of schistosoma worms [2, 20,21], it is clear that there is a latent period. Hence, the human population is divided into three compartments, i.e., susceptible humans (S h ), exposed humans (E h ), and infectious humans (I h ). In this model, we assumed that there was no latent period in snail population. Thus, the snail population was only divided into two compartments, i.e., susceptible snails (S v ) and infectious snails (I v ). We only consider two stages of schistosoma-worm development, i.e., the cercaria and miracidia stages. Therefore, the parasite is divided into two compartments, i.e., miracidiae (M) and cercariae (C). On the basis of the epidemiology of schistosomiasis [20][21][22], miracidiae can infect snails and cercariae can infect humans. It was assumed that there was only one type of snail predator in the environment, and the intermediate host (snail) was the only food for the predator. A reduction in parasites in the environment due to direct interaction with humans and snails was neglected because infectious humans can excrete large amounts of eggs that can hatch and release miracidiae, while infectious snails can excrete large amounts of cercariae [20]. We assumed that there was no recovery for infectious snails and no disease-related death in the snail and human populations. The transition and interaction between compartments are shown in Figure 1. The description of all parameters is given in Table 1.  On the basis of Figure 1, we constructed a schistosomiasis model as follows: where φ ∈ [0, 1), φ e ∈ [0, 1) and the other parameters are positive.

Invariant Region
In this subsection, we prove that the solutions of System (2) with non-negative initial value are non-negative and bounded.
Theorem 2. Solutions of System (2) are non-negative for all non-negative initial conditions. Proof. We prove this theorem by using a similar method as that used in [18,23]. From System (2), we have On the basis of Lemma 2 in [24], R 8 +0 is an invariant region of System (2). Hence, the solutions of System (2) with initial values in R 8 +0 remain in R 8 +0 . This completes the proof. Proof. On the basis of the assumptions that are used when formulating the model, we have Here, N h and N v are the total number of humans and total number of snails, respectively. On the basis of System (2), we have First, we prove that N h is bounded. It is clear that the solution of It easy to see that 0 We prove that W is bounded. The last two equations of (3) give dW dt where ξ ≥ τ and µ = min{(µ v + µ r ), µ p }. On the basis of Gronwall's lemma, we obtain According to Gronwall's lemma, we obtain C ≤ σΠ v µµ c and M ≤ αh h g h Π h µ h µ m . Thus, the solutions of System (2) are bounded. Therefore, System (2) is well-posed with invariant region Θ +

Equilibrium Points and Basic Reproduction Number
System (2) has four equilibrium points.
• Disease-free equilibrium point E a 0 .
E a 0 always exists in R 8 +0 .
• Endemic equilibrium point E a 1 .
The basic reproduction numbers are determined by using a next-generation matrix [25]. Here, we regard E h , I h , I v , C, and M as the infected compartments. Thus, we have F and V represent the new infection terms and transition terms, respectively. The Jacobian matrices of F and V at arbitrary equilibrium point (S * h , 0, 0, S * v , 0, 0, 0, P * ) are, respectively, given by The basic reproduction number is the spectral radius of FV −1 .
The characteristic polynomial of (6) is Thus, the spectral radius of FV −1 is By substituting E a 0 into (7), we obtain the basic reproduction number when the predator becomes extinct.
It is easy to see that R a 0 is completely independent from parameters that are related to the snail predator. This makes sense because E a 0 and E a 1 describe the condition when the predator becomes extinct.
After substituting E b 0 into (7), we obtain the basic reproduction number when the predator survives, namely, 0 decreases as τ increases. The higher the natural death rate of the snail predator is, the higher R b 0 is.

Stability Analysis
In this section, we investigate the stability condition of all equilibrium points. The Jacobian matrix of System (2) at E a 0 is The eigenvalues of (8) are zeroes of L(λ) where Clearly, l i for i = 1, 2, 3, 4, 5 are elementary symmetric functions [26]. It is obvious that (8) has two negative eigenvalues, i.e., The other eigenvalues are zeroes of L 1 (λ). It is easy to see that is zero. Following [27,28], we use a Routh-Hurwitz array shown in Table 2 to determine the Hurwitz determinant. Table 2. Routh-Hurwitz array associated with characteristic polynomial L 1 (λ).
Since D i > 0 for i = 1, 2, 3, 4, 5, it is clear that ∆ 2 > 0. We also proved that Ξ > 0 (see Appendix A). Thus, ∆ 4 > 0 if R a 0 < 1. Therefore, on the basis of the Liénard-Chipart criterion [29], all zeroes of L 1 (λ) have a negative real part if R a 0 < 1. Hence, all eigenvalues of J(E a 0 ) have a negative real part if R a 0 < 1 and τΠ v µ p (µ v +µ r ) < 1. Therefore, E a 0 is asymptotically stable if R a 0 < 1 and τΠ v µ p (µ v +µ r ) < 1. If R a 0 > 1, then there is one sign change in the sequence of L 1 (λ) coefficients. On the basis of Descartes' sign rule [26], there is exactly one positive real root if R a 0 > 1. Hence, E a 0 is unstable if R a 0 > 1.
Now, we present the local stability condition of E a 1 .
Proof. We use Theorem 1 to prove this theorem. Consider β ch as bifurcation parameter. We determine the bifurcation point when R a 0 = 1. The bifurcation point that is obtained . When β ch = β * ch , the characteristic polynomial (9) becomes where It is clear that characteristic polynomial (11) has simple zero root (λ 1 = 0) and two negative roots, i.e., λ 2 = −µ h and λ 3 Using the Routh-Hurwitz array in Tabel 3, we determine the conditions that guarantee that the other roots of (11) have a negative real part. Table 3. Routh-Hurwitz array associated with characteristic polynomial (12).

Column 1
Column 2 Column 3 Column 4 . Hence, on the basis of the Routh-Hurwitz criterion [27], all roots of (12) have a negative real part. These results imply that, if then J E a 0 , β * ch has one zero eigenvalue, and the other eigenvalues have a negative real part. Thus, Assumption A1 in Theorem 1 is satisfied if where v a 6 is arbitrarily positive. It is clear that v a 1 < 0 and v a 4 < 0. The left eigenvector of J E a 0 , β * ch , corresponding to a zero eigenvalue is where w a 6 is determined, such that w a . v a = 1. It is straightforward to show that w a 6 > 0. Let Now, we compute A and B that are defined in Theorem 1. It is clear that the only nonzero terms of A and B are Hence, we obtain According to Theorem 1, forward bifurcation occurs at R a 0 = 1. Consequently, endemic equilibrium point E a 1 , which exists when R a 0 > 1, is asymptotically stable if R a 0 > 1 (near 1) and τΠ v µ p (µ v +µ r ) < 1.

Theorem 6. Disease-free equilibrium point E a
0 is globally asymptotically stable if R a 0 ≤ 1 and Proof. Consider the following Lyapunov function: Thus [30], E a 0 is globally asymptotically stable if R a 0 ≤ 1 and τS a * v µ p < 1.

Now, we determine the local stability condition of
The eigenvalues of (13) are the solutions of H(λ) = 0. where It is clear that (13) has one negative eigenvalue, i.e., λ 1 = −µ h . The other eigenvalues are zeros of In the following, we apply a Routh-Hurwitz array shown in Table 4 to investigate the local stability condition of E b 0 . Table 4. Routh-Hurwitz array associated to characteristic polynomial H 1 (λ).

Column 1 Column 2 Column 3 Column 4 Column 5
It is clear that h 1 > 0. On the basis of the Routh-Hurwitz criterion [27,28], all roots of H 1 (λ) have a negative real part if the other entries in Column 1 are also positive. k 1 is always positive. It is obvious that h 7 > 0 if R b 0 < 1. Hence, all roots of H(λ) have a negative real part, which implies that Now, we investigate the local stability condition of E b 1 . We use the method that is presented in [19] to study the stability condition of E b 1 . Consider β mv as a bifurcation parameter. We determine the bifurcation point that is equivalent to R b 0 = 1. We ob- . If β mv = β * mv , characteristic polynomial (14) becomes where and It is clear that characteristic polynomial (15) has one zero root and one negative root, i.e., λ 1 = 0 and λ 2 = −µ h . Using the Routh-Hurwitz array in Table 5, we determine the conditions that guarantee that the other roots of (15) have a negative real part. Table 5. Routh-Hurwitz array associated to characteristic polynomial (16).

Column 1
Column 2 Column 3 Column 4 Column 5 On the basis of the Routh-Hurwitz criterion [27], all roots of (16) have a negative real part if h 1 , kz 1 , kz 3 , kz 5 , kz 6 , h 6 > 0. We recognize that h 1 , kz 1 , h 6 > 0. These results imply that if kz 3 , kz 5 , kz 6 > 0, then J E b 0 , β * mv has one zero eigenvalue, and the other eigenvalues have a negative real part. Hence, Assumption A1 in Theorem 1 is satisfied if kz 3 , kz 5 , kz 6 > 0. The right eigenvector of J E b 0 , β * mv corresponding to zero eigenvalue is where v b 7 is arbitrarily positive. It is clear that v b 1 < 0 and v b 4 < 0. The left eigenvector of J E b 0 , β * mv corresponding to a zero eigenvalue is Hence, we obtain According to Theorem 1, forward bifurcation occurs at R b 0 = 1. Consequently, endemic equilibrium point E b 1 that exists when R b 0 > 1 is asymptotically stable if kz 3 , kz 5 , kz 6 > 0, and R b 0 > 1 (close to 1).
Proof. Consider a Lyapunov function as follows: The time derivative of Q is dQ dt Arithmetic mean is always greater than the geometric mean. Hence, We proved the global stability condition of the disease-free equilibrium points, E a 0 and E b 0 , by formulating suitable Lyapunov functions. Since we have difficulty in finding a suitable Lyapunov function for the endemic equilibrium point, we only investigate the local stability condition of the endemic equilibrium points (E a 1 and E b 1 ). When the local stability condition of the endemic equilibrium point is met, the endemic equilibrium point is asymptotically stable, but may only attract a very small part of the state space. Therefore, studying the size of the basin of attraction of the endemic equilibrium point is relevant [31]. However, we do not discuss the basin of attraction of the endemic equilibrium point in this article. A method for numerical estimates of the size and shape of the basin of attraction, as well as the systems' return time to the attractor can be seen in [31,32].

Predator Becomes Extinct
In this subsection, we show the numerical results using τ = 0.00001 and β mv = 0.0004 365 . Here, we have τΠ v (µ v +µ r )µ p = 0.22054 < 1; thus, the existence condition of the equilibrium points says that the predator goes extinct. We can also show basic reproduction number R a 0 = 1 corresponds to the critical value of cercaria infection rate on susceptible humans, i.e., β * ch = 4.71859 × 10 −7 . If we set β ch = 1.71859 × 10 −7 < β * ch , then we obtain R a 0 = 0.60350 < 1. On the basis of Theorems 4 and 6, disease-free equilibrium point (E a 0 ) is asymptotically stable. This situation was confirmed by our simulation shown in Figure 2a where I h , I v , and P were convergent to zero. For the second simulation, we set β ch = 4.71859 × 10 −6 > β * ch which led to R a 0 = 3.16227 > 1. Theorem 5 states that endemic equilibrium point (E a 1 ) is asymptotically stable. The stability of (E a 1 ) is clearly shown in Figure 2b, namely, I h and I v converge to a positive equilibrium, whereas P goes to zero. To more clearly see the effect of the cercaria infection rate on humans, we performed numerical simulations by varying β ch from 2.35929 × 10 −7 to 7.07789 × 10 −7 , which corresponds to R a 0 between 0.70710 to 1.22474. Figure 2c indicates the occurrence of forward bifurcation driven by R a 0 or β ch . If R a 0 < 1, then E a 0 is asymptotically stable, while E a 1 does not exist. Conversely, if R a 0 > 1, then E a 0 is unstable, and there appears E a 1 , which is asymptotically stable. Bifurcation point R a 0 = 1 is achieved when β ch = β * ch = 4.71859 × 10 −7 .
(c) Figure 2. Dynamics of I h , I v , P for (a) β ch = 1.71859 × 10 −7 that corresponds to (R a 0 < 1), and (b) β ch = 4.71859 × 10 −6 that corresponds to (R a 0 > 1). (a,b) Reducing cercaria infection rate on humans by limiting the number of contacts between susceptible humans and contaminated areas or by using self-protective gear when entering contaminated areas can decrease schistosomiasis prevalence. (c) Forward bifurcation diagram of System (2) in the case of the snail predator becoming extinct, showing that schistosomiasis is eradicated if R a 0 < 1, but it becomes endemic if R a 0 > 1.

Predator Survives
We next performed numerical simulations using τ = 0.0001 and β ch = 1.914 × 10 −5 , which implies that τΠ v (µ v +µ r )µ p = 2.20543 > 1. On the basis of the existence condition of the equilibrium points, the predator survives. In this case, the basic reproduction number is unity (R b 0 = 1) if the miracidia transmission rate on snails reaches its critical value, i.e., β * mv = 1.31409 × 10 −7 . Hence, if we set β mv = 0.31409 × 10 −7 < β * mv , we obtain k 1 = 2.43742 > 0, k 4 = 0.37218 > 0, k 6 = 0.01988 > 0, k 8 = 0.00045 > 0, k 9 = 3.1 × 10 −6 > 0, and R b 0 = 0.48889 < 1. On the basis of Theorems 7 and 9, diseasefree equilibrium point (E b 0 ) is asymptotically stable. This behavior was confirmed by our simulation, where the numerical solution was convergent to (E b 0 ); see Figure 3a. On the other hand, if we take β mv = 4.31409 × 10 −7 > β * mv , such that R b 0 = 1.81188 > 1, then Theorem 8 shows that endemic equilibrium (E b 1 ) is asymptotically stable. The stability of (E b 1 ) is confirmed by Figure 3b, which shows that the numerical solution converges to (E b 1 ). When the predator goes extinct, schistosomiasis prevalence is higher than that when the predator survives. This means that the existence of a snail predator in the environment can reduce schistosomiasis prevalence.
(c) Figure 3. Dynamics of I h , I v , P for (a) β mv = 0.31409 × 10 −7 , which corresponds to R b 0 < 1; and (b) β mv = 4.31409 × 10 −7 , which corresponds to R b 0 > 1. (a,b) Reducing miracidia transmission rate on snails by improving irrigation facilities (e.g., lining canals with cement) can reduce schistosomiasis prevalence. (c) Forward bifurcation diagram of System (2) in the case of snail predator surviving, showing that schistosomiasis is eradicated if R b 0 < 1, but it becomes endemic if R b 0 > 1. Figure 3a,b indicate that there is an exchange of stability between the disease-free equilibrium point and the disease endemic equilibrium point when the miracidia transmission rate on snails is changed from β mv = 0.31409 × 10 −7 < β * mv to β mv = 4.31409 × 10 −7 > β * mv . To more clearly see the effect of β mv , we plotted in Figure 3c the steady state of I h against R b 0 . R b 0 between 0.70710 to 1.22474 corresponding to β mv from 0.65704 × 10 −7 to 1.97114 × 10 −7 ; R b 0 = 1 is obtained when β mv = β * mv = 1.31409 × 10 −7 . Figure 3c shows that System (2) experiences forward bifurcation. Thus, E b 0 is asymptotically stable, and E b 1 does not exist when R b 0 < 1. If R b 0 > 1, E b 0 becomes unstable, and E b 1 exists and it is asymptotically stable.

Impact of Snail Predator as Biological Control Agent
We now investigate the impact of a snail predator as biological control agent. Here, we study the effect of predation rate (ξ) by using parameter values as in Table 6: τ = 0.0001, β ch = 0.1 and β mv = 4.31409 × 10 −7 . According to [33], the predation rate is the proportion of prey killed per predator per time. It suggests that the predation rate is related to the effectiveness of the predator in hunting and killing snails. In this simulation, the basic reproduction number does not depend on the predation rate and is given by R b 0 = 18.11886. Figure 4 shows that all numerical solutions using ξ = 0.1, ξ = 0.15, and ξ = 0.3 are convergent to the same endemic equilibrium point. However, detailed observation shows that the schistosomiasis prevalence in the beginning of the intervention decreased as predation rate (ξ) increased. Several studies that are related to satiation-based predation models show that the predation rate is related to the gut capacity of the predator [34]. Hence, we should use snail predators that can effectively hunt and kill snails, and have high gut capacity. We next investigate the influence of conversion rate τ by performing simulations using the same parameter values as those in Figure 4, but with varying τ and fixed ξ = 0.01. According to [35], the conversion rate is related to the efficiency in turning predation into new predators. It implies that the conversion rate is related to the birth rate of snail predator. By choosing τ = 0.0001, τ = 0.00015, and τ = 0.0003, we get the basic reproduction number (R b 0 ) = 18.118869, (R b 0 ) = 12.07924, (R b 0 ) = 6.03962, respectively. Hence, the basic reproduction number decreases as τ increases. The dynamics of I h , I v , C, and M for τ = 0.0001, 0.00015, and 0.0003 is shown in Figure 5. Values of the steady state of I h , I v , C, and M were smaller for larger values of τ. Thus, schistosomiasis prevalence decreases as τ increases. So, the snail-predator birth rate plays an essential role in controlling schistosomiasis spread. Figures 4 and 5 show that using a natural predator of snails as a biological control agent can decrease schistosomiasis prevalence. These results are in line with the results in [7].  showing that the conversion rate of the snail predator that is related to its birth rate plays an essential role in reducing schistosomiasis prevalence.

Conclusions
In this work, a deterministic schistosomiasis model incorporating a snail predator as a biological control agent was discussed. The existence and stability conditions of all equilibrium points were investigated. Our findings suggest that a snail predator as a biological control agent can reduce the prevalence of schistosomiasis. Moreover, our results showed that the snail-predator birth rate plays an important role in controlling schistosomiasis spread. Despite the contributions of our work, it has some limitations. Due to the lack of data related to the parameters of our model, we only chose parameter values by considering the results of qualitative analysis, i.e., the existence and stability conditions of the equilibrium points. For further research, sensitivity and cost-effectiveness analysis of the interventions discussed in our proposed model will be considered.

Acknowledgments:
The authors thank the editors and the reviewers for their critical suggestions, contributions, and comments that greatly improved this manuscript.

Conflicts of Interest:
The authors declare no conflict of interest in this paper.