Multi-host transmission dynamics of schistosomiasis and its optimal control.

In this paper we formulate a dynamical model to study the transmission dynamics of schistosomiasis in humans and snails. We also incorporate bovines in the model to study their impact on transmission and controlling the spread of Schistosoma japonicum in humans in China. The dynamics of the model is rigorously analyzed by using the theory of dynamical systems. The theoretical results show that the disease free equilibrium is globally asymptotically stable if R0 < 1, and if R0 > 1 the system has only one positive equilibrium. The local stability of the unique positive equilibrium is investigated and sufficient conditions are also provided for the global stability of the positive equilibrium. The optimal control theory are further applied to the model to study the corresponding optimal control problem. Both analytical and numerical results suggest that: (a) the infected bovines play an important role in the spread of schistosomiasis among humans, and killing the infected bovines will be useful to prevent transmission of schistosomiasis among humans; (b) optimal control strategy performs better than the constant controls in reducing the prevalence of the infected human and the cost for implementing optimal control is much less than that for constant controls; and


(Communicated by Zhilan Feng)
Abstract. In this paper we formulate a dynamical model to study the transmission dynamics of schistosomiasis in humans and snails. We also incorporate bovines in the model to study their impact on transmission and controlling the spread of Schistosoma japonicum in humans in China. The dynamics of the model is rigorously analyzed by using the theory of dynamical systems. The theoretical results show that the disease free equilibrium is globally asymptotically stable if R 0 < 1, and if R 0 > 1 the system has only one positive equilibrium. The local stability of the unique positive equilibrium is investigated and sufficient conditions are also provided for the global stability of the positive equilibrium. The optimal control theory are further applied to the model to study the corresponding optimal control problem. Both analytical and numerical results suggest that: (a) the infected bovines play an important role in the spread of schistosomiasis among humans, and killing the infected bovines will be useful to prevent transmission of schistosomiasis among humans; (b) optimal control strategy performs better than the constant controls in reducing the prevalence of the infected human and the cost for implementing optimal control is much less than that for constant controls; and (c) improving the treatment rate of infected humans, the killing rate of the infected bovines and the fishing rate of snails in the early stage of spread of schistosomiasis are very helpful to contain the prevalence of infected human case as well as minimize the total cost. parasitic disease. Due to the great harmfulness of schistosomiasis to public society, it is an important issue to increase our understanding of schistosomiasis transmission dynamics and further to take effective measures in containing or eliminating its transmission. Mathematical models are powerful tools for gaining insights into the transmission and control of infectious diseases, and can be used to address such an important issue. There have been a lot of mathematical models to investigate transmission dynamics of schistosomiasis ( [5,6,15,21,23] and references therein). These models have provided useful information for understanding the mechanics of schistosomiasis transmission.
Most of the schistosomiasis models mentioned above are the snail-schistosomehuman system which considered only one definitive host, i.e., human beings. However, Schistosome is a multi-host parasite, and the parasite is indirectly transmitted by snails, intermediate host, between definitive hosts which include not only humans but also many other mammals that may act as zoonotic reservoirs [14]. Perhaps it is partly because of the presence of zoonotic reservoirs that schistosomiasis still remains a major public health concern in China, despite the control program of schistosomiasis in China launched in 1950 and sustained over 60 years. Thus, in order to get a comprehensive understanding of schistosomiasis transmission dynamics, much more attentions should be paid to the potential role of zoonotic reservoirs. Assessing the potential role may help guide us to reduce and eliminate incidence of schistosomiasis in humans effectively. Over the past decades, studies on the epidemic dynamics between human and other mammals have been given abroad attentions ( [7,9,10,13,22] and the references therein).
Besides, most of the papers mentioned above assumed that the recruitment rates of the definitive host and intermediate host are constants. In fact, the different species have different recruitment rates. For example, the number of humans in an area can be regarded as a constant due to the balance of natural birth and death; the snails may grow exponentially since the natural capacity of snails is sufficiently large, and the bovines may migrate to different places due to the business transaction. Thus it is also necessary to incorporate the different growth rates for different hosts in order to make the model more reasonable.
In this paper, we will formulate a mathematical model to describe transmission dynamics of schistosomiasis in China. Because substantial literatures recognize that bovines play a major role in the transmission of Schistosoma japonicum to human in China [9], we will include two definitive hosts, humans and bovines, and intermediate host, snails. Then we will use the model to investigate the multi-host transmission dynamics of schistosomiasis, to assess the impact of bovines on human schistosomiasis transmission, and to evaluate the effects of various control measures against schistosomiasis. In recent years, optimal control theory has been often applied to epidemic models ( [1,3]). In this paper, we will also apply optimal control theory to the multi-host schistosomiasis model to further study the corresponding optimal control problems.
The paper is organized as follows. Section 2 introduces the epidemic model to describe the transmission dynamics of schistosomiasis which incorporates two definitive hosts, bovines and humans, with different recruitment rates. In section 3, we calculate the basic reproduction number and stability analysis of the 9 dimension schistosomiasis system. In Section 4, we study the local and global stability analysis of the model. Section 5 derives the existence and uniqueness for the optimal control strategies. We propose some numerical simulations, furthermore summarize some suggestions for schistosomiasis prevention and control in Section 6. In the last section, we conclude with a brief discussion of our theoretical and numerical results.
2. Model formulation. In this section we mainly formulate a mathematical model to describe transmission dynamics of Schistosomiasis japonica in China. Because a number of studies reported that bovines are the primary infection source of Schistosoma japonicum transmission in China [25], and the other mammals, such as dogs, pigs, rats and goats, are likely to contribute only minimally to overall transmission [9], in this paper we only consider the scenario that Schistosomiasis japonica are transmitted among two definitive hosts, humans and bovines, and one intermediate host, snails.
According to the transmission of Schistosoma japonicum, we assume that the total population of humans, N h (t), is divided into three epidemiological subclasses which are susceptible, infected and recovered, with sizes denoted by H s (t), H i (t) and H r (t), respectively. Similarly, the total populations of bovines, N b (t), and snails, N s (t), are divided into susceptible subclass and infected subclass, and their sizes are denoted by B s (t), B i (t), S s (t), S i (t), where B, S represent bovine and snail and s, i represent the subclasses of susceptible and infected. We use M (t) and C(t) to denote the mean spatial density of miracidia and cercaria.
The changes of human population are due mainly to immigrant and emigrant, and the changes of human population are relatively slow in a short duration, thus in this paper we assume that the recruitment rate of suspectable human is Λ h and per capita natural mortality rate of human is d h . Similarly, the changes of bovine population are also mainly due to immigrant and emigrant, and the changes of bovine population are also relatively slow in a short duration, thus in this paper we also assume that the recruitment rate of suspectable bovine is Λ b and per capita natural mortality rate of bovine is d b . The changes of snail population are mainly due to birth and death, thus in this paper we assume that the susceptible snails are recruited into the population at a rate b s (S s + qS i )(1 − Ns Ks ) because the schistosomiasis infection can reduce the reproduction of snails. The natural mortality rate of snail is assumed to be d h .
A transition diagram between these epidemiological classes is shown in Figure 1. Susceptible humans, bovines and snails acquire infection with schistosomiasis at the rates of β ch γ ch C, β cb γ cb C and β ms γ ms M , respectively, where β ch (β cb , β ms ) is the transmission probability from cercaria (cercaria, miracidia) to susceptible humans (bovines, snails), and γ ch (γ cb , γ ms ) is the contact rate between humans (bovines, snails) and cercaria (cercaria, miracidia). Infected humans are treated at the rate θ ir and recovered humans may lose immunity at the rate of θ rs . α h , α b , α s denote the per capita death rates of humans, bovines and snails induced by infection of schistosomiasis, respectively. Let b be the water surface area, σ be the average number of cercaria released by an infected snail per capita per unit time, and d c and d m be the natural death rates of cercaria and miracidia, respectively. The average number of eggs that are excreted by an infected human (bovine) per unit time is assumed to be the product of the amount of fecal output, g h (g b ), the average number of schistosomiasis eggs per gram of stool, h h (h b ), and hatch ability from the eggs into miracidia, α.
Based on the transition diagram in Figure 1, the model is described by the following system of 9 ordinary differential equations: Human equations : Bovine equations : Snail equations : (2.1m) 3. Reproduction number and steady states. Consider the case b s > d s , and let Then it can be seen that all solutions of system (1) starting in Γ remain in Γ for all t ≥ 0. In what follows, we always assume that the initial points lie in Γ. Straight forward computation yields that system (1) always has the boundary , 0, 0, 0, 0, 0), which is the snail free equilibrium (SFE).
After a simple calculation, we can easily prove that if b s ≤ d s the SFE E 00 is globally asymptotically stable in Γ. Thus, in this paper, we always assume that b s > d s . It is obvious that the SFE E 00 is globally asymptotically stable in Γ 00 = {(H s , 0, 0, B s , 0, 0, 0, 0, 0) ∈ Γ}, but it is unstable in Γ. Moreover, system (1) has the , 0, N s , 0, 0, 0), which is the disease free equilibrium (DFE). Similarly, it is easy to prove that E 0 is globally asymptotically stable in Γ 0 = {(H s , 0, 0, B s , 0, S s , 0, 0, 0) ∈ Γ}, but its local stability is completely determined by the basic reproduction number of system (1), which can be found from the next generation matrix for system (1).
Noting that system (1) has 5 infected populations, namely H i , B i , S i , C and M , it follows that, using the notation of Van den Driessche P. and Watmough J. [18], the matrices F and V, for the new infection terms and the remaining transfer terms, respectively, are given by A straightforward computation yields that the next generation matrix for schistosomiasis is After extensive algebraic calculations, the characteristic equation associated with FV −1 is Results in [18] imply that the basic reproduction number of system (1) is

CHUNXIAO DING, ZHIPENG QIU AND HUAIPING ZHU
where ρ(FV −1 ) represents the spectral radius of the matrix FV −1 , and ; ; Now let us explain the biological interpretations of R 0 . If system (1) does not have any infected definitive host and intermediate host and the system is in balance, then the numbers of the susceptible humans, bovines and snails and N s , respectively. Under these conditions, the average increased spatial density of cercaria generated by infection of one infected snail can be expressed as represents the average increased density of cercaria that one infected snail will lead to per unit of time and 1 αs+ds is the mean lifespan of the infected snail. Similarly, the average number of the infected humans generated by per unit spatial density of cercaria can be defined by represents the average infected humans caused by per unit spatial density of cercaria in per unit of time and 1 β ch γ ch + β cb γ cb + d c is the mean lifespan of cercaria. Noting that R sh 0 = R sc R ch , it then follows that the parameter R sh 0 can be explained as the average number of the infected humans generated by one infected snail when the system does not have any infected definitive host and intermediate host and is in balance.
In the same way, the parameters R hs 0 , R sb 0 and R bs 0 can be explained as the average number of infected snails, bovines and snails generated by one infected human, snail and bovine, respectively, when the system does not have any infected definitive host and intermediate host and is in balance. Therefore, the total number of secondary infected snails generated by transmission of one infected snail can be defined by when the numbers of the susceptible humans, bovines and snails Now we are able to state the results on the existence of equilibria for system (1).
Theorem 3.1. Assume that b s > d s . The system (1) always has two boundary equilibria E 00 , E 0 , and 1) if R 0 ≤ 1, system (1) has no positive equilibrium; 2) if R 0 > 1, system (1) has a unique positive equilibrium Proof of Theorem 3.1. The existence of the boundary equilibria E 00 , E 0 has been illustrated at the beginning of this section. In the following, we only need to prove the results on the existence of the positive equilibrium. A positive equilibrium of system (1) is a positive solution to the equations: where From the first, third and sixth equations in (2), we have By substituting the expressions for N h , H r , C into the second equation in (2) we obtain Similarly, it follows from the fourth, fifth and sixth equations in (2) that we have Putting the expression for M into the eighth and ninth equations in (2), the positive equilibrium, if exists, is the interaction of the two curves With extensive algebraic manipulations, we have

and the fact that
was used in the last step. Thus the function N s = N s (S i ), S i ∈ (0, S i ) is increasing and concave, and the function is inverse. Let us denote its inverse function by , and it then follows from the prosperities of inverse function that the function is a convex and increasing function. Now let us consider two cases: We can easily verify that the condition N s (0) > N s implies that R 0 ≤ 1. In this case, since the intersection of the interval (0, N s ) and interval (N s (0), N s ( S i )] is empty, the two curves C 1 and C 2 do not meet each other (see the first diagram in Figure 2). Thus if R 0 ≤ 1 the system (1) has no positive equilibrium.
Case 2. N s (0) < N s . This condition implies that R 0 > 1. In this case, we have Direct mathematical analysis show that the two curves C 1 and C 2 meet at only one point due to the fact that the curve C 1 is concave and the curve C 2 is convex (see the second graph in Figure 2). Thus the equations (2) has only one positive solution if R 0 > 1, i.e., system (1) has a unique positive equilibrium. This completes the proof of Theorem 3.1. Figure 2. Graphs for the function C 1 and the inverse function of C 2 4. Disease dynamics. In this section, we mainly analyze the local and global stability of system (1). Using Theorem 2 in [18], we can easily obtain the following stability result.
In fact, the local stability of DFE E 0 implies the global stability, and we have the following result.
Proof of Theorem 4.2. If R 0 < 1, it follows Theorem 4.1 that the DFE E 0 is locally asymptotically stable. In the following, we only need to prove that E 0 is a global attractor.
Let us consider a positive solution of system (1). From the equations in system (1), it follows that

CHUNXIAO DING, ZHIPENG QIU AND HUAIPING ZHU
Consider the following differential equations We can easily verify that system (5) is a cooperative irreducible linear system in R 5 + , and then the global stability of the origin of system (5) is completely determined by the stability of Jacobian matrix J = F − V. If R 0 < 1, Theorem 2 in [18] implies that the matrix J is stable. Then we have By the comparison principle it then follows that H i (t) → 0, B i (t) → 0, S i (t) → 0, C(t) → 0, M (t) → 0 as t → +∞. Substituting them into system (1), we can obtain the limiting system Solving the above equations yields that (6) is the limiting system of (1), it follows from Theorem 2.3 in paper [4] that the DFE E 0 is a global attractor of system (1). This completes the proof of Theorem 4.2. Now let us consider the local stability of the unique positive equilibrium E * when R 0 > 1. Linearizing system (1) around the positive equilibrium E * , we obtain the following Jacobian matrix , a 66 = β ch γ ch +β cb γ cb +d c , a 77 = β ms γ ms +d m , a 87 = β ms γ ms (N * S − S * i ), a 88 = β ms γ ms M * + (α s + d s ). After extensive algebraic calculations, its characteristic equation is given by where where ♦ is a positive constant. The calculation process for det(J(E * )) is detailed introduced in the Appendix A. Since A 9 (E * ) > 0, the relations between the roots and the polynomial coefficients imply 9 i=1 λ i = −A 9 (E * ) < 0, where λ i , i = 1, 2, · · · , 9, are the roots of the characteristic equation (7). Then it follows that the equation (7) has no zero roots and the number of roots with nonnegative real parts are even. Summarizing the above conclusions, we have Theorem 4.3. If R 0 > 1, the unique positive equilibrium E * of system (1) is either locally asymptotically stable, or the dimension of its unstable and center manifold are both even.
We are interested in the global behavior of system (1), and many numerical simulations show that the positive equilibrium is globally asymptotically stable so long as it exists, i.e., R 0 > 1. We conjecture that the simulation result holds for the system, but unfortunately we can not rigorously classify the global dynamics of system (1) if R 0 > 1. In the rest of this section, we will investigate the global behavior of system (1) for a special case that q = 1, i.e., the trematode infection completely inhibits snail reproductive activity. Sufficient conditions are provided in the following theorem for the special case.
In order to prove Theorem 4.4, let us consider the following ordinary differential equations: where f : U × ∆ → R n is continuous, where U ∈ R n and ∆ ∈ R k and D x f (x, λ) is continuous on U × ∆. We assume that solutions of initial value problems are unique and remain in U for all t > 0 and each λ ∈ ∆. We write x(t, z, λ) for the solution of (9) satisfying x(0) = z.
have a negative real part, and x 0 is globally attracting for solutions of (9) with λ = λ 0 . If there exists a compact set K ⊂ U such that for each λ ∈ ∆ and each z ∈ U , x(t, z, λ) ∈ K for all large t, then there exist > 0 and a unique pointx(λ) ∈ U for λ ∈ B ∆ (λ 0 , ) such that f (x(λ), λ) = 0 and x(t, z, λ) →x(λ) as t → +∞ for all z ∈ U . Now we are able to prove Theorem 4.4.

The asymptotic equilibrium values for
The limiting system of system (1) is We can easily check that the Jacobian matrix of system (10) admits positive offdiagonal elements, and then it follows that system (10) is a cooperative irreducible system in R 5 + . The system (10) has only one boundary equilibrium E 0 (0, 0, 0, 0, 0), and the Jacobian matrix of system (10) at E 0 can be expressed asĴ(E 0 ) = F − V, where the matrices F and V are defined in Section 3. From the proof of Theorem 2 in [18] it follows that if R 0 > 1, the boundary equilibrium E 0 of system (10) is unstable. Following Smith [16], system (10) is strongly concave. Then it follows that if R 0 > 1 system (10) has an equilibrium E * , which is globally asymptotically stable. Since system (10) is the limiting system of (1), it follows from Theorem 2.3 in paper [4] that the unique positive equilibrium E * is a globally asymptotically stable equilibrium of system (1) when α h = 0, α b = 0, α s = 0,µ 2 = 0 and θ ir = 0. Now we extend the above conclusion by using perturbation arguments to show that α h = 0, α b = 0, α s = 0,µ 2 = 0 and θ ir = 0 can be replaced by nearly zero. Using a similar argument as in the proof of Theorem 2.3 in [19], we can conclude that if R 0 > 1 system (1) is uniformly persistent, i.e., there exists a constant δ > 0 such that This implies that we can take for the compact set K in Lemma 4.5 the set It follows from the above results that E * is global attracting for the case that α h = 0, α b = 0, α s = 0 and θ ir = 0. Furthermore, straightforward computation yields that all eigenvalues of the Jacobian matrix of system (1) at E * have negative real parts when α h = 0, α b = 0, α s = 0 and θ ir = 0. It follows from Lemma 4.5 that if R 0 > 1 for sufficiently small α h , α b , α s and θ ir the unique positive equilibrium E * of system (1) is globally asymptotically stable. This completes the proof of Theorem 4.4.

5.
Optimal control. In this section, we try to implement anti-Schistosomiasis control to protect humans and bovines while minimizing the total cost. In fact, the problem is a typical optimal control problems.
In this paper, we assume the control variables in the set in which all control variables are bounded and Lebesgue measurable and U i , i = 1, 2, 3 is denoted to the upper bound of the control variables. In our controls, the control function µ 1 (t) represents the enhanced drug treatment rate of infected humans, and the additional recovered rate of infected human due to enhancing treatment could be assumed to µ 1 (t) in this paper. µ 2 (t) indicates the killing rate of infected bovines. As well as the control function µ 3 (t) is the fishing rate of snails.
Consequently, additional reduced rates of infected bovines and snails due to control are represented by µ 2 (t) and µ 3 (t), respectively. In this paper, we choose the objective (cost) function where, a i , i = 1, 2 are positive constants representing the weight of infected humans and infected bovines, respectively. b i , i = 1, 2 are weight constants for human and bovines inspection. c i , i = 1, 2, 3 are weight constants represent for the costs of drug treatment of infected human, killing bovines and snails fishing. We assume that the costs are proportional to quadratic form of their corresponding control functions. Moreover, the coefficients in objective function, a i , b i , i = 1, 2, c i , i = 1, 2, 3, not only represent the weights, but also balance the different units in objective function because the magnitudes of the populations of humans, bovines and the control functions absolutely are on different scales. Our purpose is to find an optimal control pair (u * 1 (t), u * 2 (t), u * 3 (t)) in order to seek the minimum value of objective function J(µ * 1 (t), µ * 2 (t), µ * 3 (t)), such that subject to the the system given by (12). Now we derive the necessary conditions that a pair control and corresponding states must satisfy. By using the same method in [8], the existence of optimal control can be proved. In the above minimizing problem, we can easily verify that the objective function is convex on the closed, convex control set U . The optimal system is bounded which can determines the compactness needed for the existence of the optimal control.
In order to find the optimal solution of system (12), first let us define Hamiltonian functions H for the optimal control system (12) as With the existence of optimal control system, we now present and discuss the adjoint system and the characterizations of the optimal control system. For simplicity, we denote and λ = (λ 1 , λ 2 , · · · , λ 9 ).
Furthermore, an optimal control could be obtained Proof of Theorem 5.1. The Pontryagin's Maximum Principle [11] is used to find an optimal solution.
Applying the adjoint conditions to the Hamiltonian (13) with X = X * , that is

998
CHUNXIAO DING, ZHIPENG QIU AND HUAIPING ZHU (14) is obtained. The optimal conditions at U * could be calculated as follows.
Using the bounds on the controls, we can obtain the optimal control solutions of system (16) . The proof is completed.
6. Numerical simulations. In this section, we mainly present some numerical simulation results, which assess the impact of bovines on human Schistosoma japonicum transmission and explore the effect of various control measures against schistosomiasis.
For the numerical simulations, we choose the values for the parameters based on available information in the literature as follows: • , h h = 100, g b = 10000, h b = 240. Firstly, we use two measures to assess the impact of bovines on human Schistosoma japonicum transmission. One is the direct effect of bovine on human Schistosoma japonicum transmission, measured by the population attributable fraction (PAF) [2] of schistosomiasis infected human cases which transmitted from the bovines indirectly, defined as P AF (t) = Incidence of schistosomiasis infected human cases which transmitted from the bovines directly at time t Total incidence of schistosomiasis infected human cases at time t .
This represents the effect of bovines on schistosomiasis infected human cases at time t, the fraction of new schistosomiasis infected human cases that are transmitted from the bovines directly. The other one is the excess prevalence of schistosomiasis infected human cases, which is defined as the difference between human schistosomiasis prevalence when the system includes the bovines and the counterfactual scenario in which the system does not include the bovines. Furthermore, in order to assess the effect of killing bovines on human Schistosoma japonicum transmission, we also measure human schistosomiasis prevalence for different killing rates of infected bovines. Figure 3 illustrates the population attributable fraction (PAF) of schistosomiasis human incidences. Simulation results indicate that most of the new infected human cases are attributed indirectly to the infected bovines, and, at the endemic steady state, roughly 90% of the new human cases infected with schistosomiasis may be attributed to the infected bovines. Figure 4 shows the human schistosomiasis prevalence for the system with bovines and without bovines. Simulation results   show that the prevalence of infected human cases can reach 38% if the system has the bovines but the prevalence of infected human cases will keep at around 9% if the system does not have the bovines. The excess prevalence of schistosomiasis infected human cases indicates that the infected bovines may play an important role in the spread of schistosomiasis among humans. Figure 5 shows the human schistosomiasis prevalence for different killing rates of bovines. It follows from Figure 5 that if killing the infected bovines keep at the rate of 0.05, the human schistosomiasis prevalence will be reduced from 45% to 20%, and if killing the infected bovines keep at the rate of 0.5, the human schistosomiasis prevalence will fall to below 10%. These results imply that killing the infected will be helpful to prevent transmission of schistosomiasis among humans.  Next we examine the effect of various control measures on human Schistosoma japonicum control. Since the most important thing is to reduce the number of infected human cases, and the next important is to reduce the number of infected bovines, it is reasonable to assume that a 1 > a 2 . Based on these assumptions, in this paper the weights in the objective function are chosen for illustration purpose   Figure 6 illustrates the effectiveness of optimal control in comparison to two constant controls, µ 1 = µ 2 = µ 3 = 0.01 and µ 1 = µ 2 = µ 3 = 0.05. We can conclude that the optimal control strategy is a more beneficial choice in fighting the outbreak of schistosomiasis. The schistosomiasis infected human prevalence decreases rapidly under the optimal control strategy and then the prevalence will keep at the very low level. The values of objective functions for the optimal control strategy, µ 1 = µ 2 = µ 3 = 0.01 and µ 1 = µ 2 = µ 3 = 0.05 are 1.0514 × 10 6 , 1.4541 × 10 6 and 1.2070 × 10 6 respectively. We can easily see that the cost for implementing optimal control is much less than that for constant controls. The time dependent optimal control law is shown in Figure 7. From Figure 7 we can conclude that in the early stage of spread of schistosomiasis, improving the treatment rate of infected human µ 1 , the killing rate of the infected bovines µ 2 and the fishing rate of snails µ 3 , are very helpful to reduce the prevalence of infected human cases. At about 65 days, the control of intensive therapy for infected humans µ 1 begins to decrease with the reduce of Lagrangian function. Then the killing rate of the infected bovines µ 2 glides down around 120 days, and the Lagrangian function declines correspondingly. By comparing Figure 7 with Figure 8, we can see that, the controls implementing on humans and bovines are more effective than that on snails with respect to the decrease of cost function, and the arrows in these two figures appear in pair at every turning point. With the transversality condition λ(T ) = 0, all the control laws jump back to 0 in the end.
7. Discussion. In this paper, a multi-host dynamic model is formulated to describe the transmission dynamics of schistosomiasis, and the model includes two definitive hosts, humans and bovines, and intermediate host, snails. By using the method of next generation matrix [18], the expression for the basic reproduction number R 0 is derived, and the interpretation of R 0 is also explained biologically. Then in terms of the basic reproduction R 0 the dynamical behavior of the 9-dimensional nonlinear system is rigorously analyzed by using the theory of dynamical systems. The theoretical results show that if R 0 ≤ 1 the disease free equilibrium is globally asymptotically stable and if R 0 > 1 the system has only one positive equilibrium. In the case that R 0 > 1 the local stability of the unique positive equilibrium is investigated which shows that the dimension of its unstable and center manifold are both even. A great many of simulations show that the unique positive equilibrium is global asymptotically stable as long as it exists, but in this paper we only provide sufficient conditions for the global stability of the positive equilibrium.
Furthermore, the optimal control theory are applied to the multiple host schistosomiasis model to study the corresponding optimal control problem. By using the Pontryagin's Maximum Principle necessary conditions are provided for the existence of the optimal solution to the optimal control problem. Finally, numerical simulations are presented to verify the theoretical results, to assess the impact of bovines on human schistosomiasis transmission and to examine the effect of various control measures against schistosomiasis. From the numerical results, we can conclude that: 1) The infected bovines play an important role in the spread of schistosomiasis among humans, and killing the infected bovines will be helpful to contain the transmission of schistosomiasis among humans.
2) Optimal control strategy performs better than the constant controls in reducing the prevalence of the infected human and the cost for implementing optimal control is much less than that for constant controls.
3) Improving the treatment rate of infected human, the killing rate of the infected bovines and the fishing rate of snail in the early stage of schistosomiasis spread are very helpful to contain the prevalence of infected human cases and also minimize the total cost.
In this paper, the model presented here only involves a single strain of schistosomiasis in the host populations that includes two definitive hosts, human and bovines, and intermediate host, snails. It is interesting to model multi-strains of schistosomiasis in host population and to study the evolutionary questions on schistosomiasis. Our model is described by ordinary differential equations that do not include age structure of the human host population. However, lots of epidemiological evidences show that the human age structure certainly affects the transmission dynamics of schistosomiasis. Consequently, we should incorporate age structure into our modeling of the transmission dynamics of schistosomiasis. Finally, since parasites exhibit a wide degree of variability between their hosts, we also should incorporate the parasite distribution pattern into our model and analyze the consequences of such pattern on transmission dynamics of schistosomiasis. We keep these consideration for our future work.
By differentiating the function Υ(S i , N s ) with respect to S i , we have On the other hand, we have obtained continuously differential functions M = M (S i )S i in (3) Then it follows that since the function M (S i ) is a decreasing function of S i . The fact that β ms γ ms (N * s − S * i )M (S * i ) = µ 3 + α s + d s was used in the last step. From (19) and (20) we can easily see that det A > 0 since D 88 < 0.
Since Υ(Si,Ns) ∂Si | (S * i ,N * s ) < 0, it then follows from the implicit function theorem we conclude that there exist continuously differential function S i = S i (N s ) defined on a neighborhood Θ of N * s such that (1) .
As in the previous discussion, on the one hand, substituting the function S i = S i (N s ), N s ∈ Θ, into the expression f 9 (S i , N s ) yields that Ξ(N s ) := f 9 (S i (N s ), N s ). Differentiating the function Ξ(N s ) with respect to N s , we have