DYNAMICS OF A DELAY SCHISTOSOMIASIS MODEL IN SNAIL INFECTIONS

In this paper we modify and study a system of delay differential equations model proposed by N̊asell and Hirsch (1973) for the transmission dynamics of schistosomiasis. The modified stochastic version of MacDonald’s model takes into account the time delay for the transmission of infection. We carry out bifurcation studies of the model. The saddle-node bifurcation of the model suggests that the transmission and spread of schistosomiasis is initial size dependent. The existence of a Hopf bifurcation due to the delay indicates that the transmission can be periodic.


1.
Introduction. Schistosomiasis is a parasitic disease affecting approximately 200 million people worldwide. There is an estimated 650 million people at risk of infection in 76 countries and territories [22]. Schistosomiasis transmission has re-emerged in 38 counties including seven provinces in China. Those counties, having reached a critical point of transmission, came under control and the increase was presented in focal areas along the Yangtze River. This was due to changes of environmental, ecological, societal and economic status, as well as on the forces of control [27]. The transmission of this infection depends on a complex cycle of events which includes humans (the definitive host), certain parasitic flatworms (the schistosome) and particular species of snails (the intermediate host). In this paper, we shall focus on attention to schistosoma mansoni and schistosoma haematobium. Particularly, in endemic foci of infection where human is the only significant vertebrate host.
The life cycle of the schistosome is described in details in [12,20]. There are two adult forms of schistosome, male and female which mate heterosexually. In the liver of the definitive host mating occurs. Some of the fertilized eggs succeed in passing through the wall of the blood vessel in which they are laid, hence through the wall of intestine or bladder until they fall into the lumen of organ and are voided with faeces or urine. Some of the eggs ultimately may be deposited in fresh water streams or lakes, where small ciliated larvae (miracidia) emerge. They swim actively and if one comes into contact with an appropriate molluscan host (snail), it rapidly penetrates the snail tissue. By a peculiar process of repeated asexual multiplication within the snail, thousands of a second larval form (cercariae) are produced. When mature, the cercariae are shed by the snail and enter a free-swimming state designed for invasion of the definitive host. On coming into contact with a definitive host a cercaria attaches itself to the skin and quite rapidly penetrates it, while at the same time sloughing larval structures to become a juvenile schistosome. This is followed by migration to the liver, maturation, mating, migration to the permanent vascular abode, copulation, and oviposition, which begins the cycle all over again. Viewed schematically the transmission dynamics of schistosomiasis can be seen to depend on two interrelated flows: one a drift of eggs from definitive host to snail, the other a stream of cercariae from snail to definitive host.
During the last three decades, mathematical models have been developed and analysis performed on the study of the transmission dynamics of such parasitic infection. The first schistosomiasis mathematical model, was proposed and studied by MacDonald [17]. A stochastic version of MacDonald's model was constructed by Nåsell and Hirsch [19]. In one of the early modeling studies of schistosomiasis in China, Wu et al. [24] built and analyzed a model of schistosomiasis with delay, and later Wu and Feng [23] generalized the model to multiple definitive hosts. Recently, there are several studies addressing the mating structure of schistosomiasis, see [2] and reference therein. A general model incorporating acquired immunity and contact pattern with infested water is developed by Yang in [25]. Mathematical models that include mass chemotherapy in human hosts with explicit snail dynamics were formulated and studied by Feng et al. in [8]. Zhang et al. [26] proposed a schistosomiasis model with an age-structure in human hosts and studied the treatment strategies. Some models describing the schistosomiasis dynamics involving two migrating human groups were showed in [9]. However most models ignore the fact that the transmission process requires a definite length of time.
In this paper, we consider a mathematical model of schistosomiasis which takes into account the time delay for the transmission of infection. The model is studied analytically in terms of steady states, stability and bifurcation. By analysing the transcendental characteristic equation we study the effect of the time delay on the local stability of the equilibrium. We also study the stability of the equilibria when it is non-hyperbolic, and a saddle-node bifurcation is found. We find the Hopf bifurcation which can occur in the model when the threshold value increases through certain values of the delay, and numerical simulations are carried out for the periodic solution. By means of Lyapunov functionals, sufficient conditions are derived for the global asymptotic stability of the disease-free equilibrium of the model. This paper is organized as follows. In section 2 we introduce the model with delay, and study the basic mathematical properties of the solutions of the model equations after having supplemented it by the appropriate initial conditions. In section 3 we study the local stability of the equilibria, and in some critical situation the bifurcation analysis is carried out in section 4. In section 5 we deals with the global stability of the disease-free equilibria. In section 6 we simulate the periodic solution bifurcated from one of the endemic equilibria in a special case that the initial distributions of the number of male and female worms in definitive host follow the poisson distribution. The conclusion and discussion are in the final section.
2. The model with delay. Our model with time delay is based on the early work of Nåsell and Hirsch [19]. The original model has the following form: The authors [19] considered an idealized focus of infection-a relative isolated community whose definitive hosts (humans) are equally exposed to the risk of infection and are not subjected to the process of birth, death, immigration, nor emigration. Birth and death, but not immigration or emigration, will be assumed to occur in the intermediate host population under the simplifying hypothesis that at the instant a snail dies and an uninfected snail is born. We denote by N 1 and N 2 , respectively, the sizes of the definitive and intermediate host populations. In Eq.(1), β(t) represents the expected worms burden per definitive host during the time interval [0, t) and survive to time t, y(t) represents the expected proportion of infected intermediate hosts at time t. The term 1 2 ψ(β(t)) can be interpreted as the expected worm pair burden in a definitive host initially uninfected and the term F (e −µ1t , β(t)) represents a transient perturbation due to the magnitude of the initial distribution of the infection at time t. The quantities T 1 = ν11ν12 µ1 N 2 and T 2 = ν21ν22 µ2 N 1 have suggestive interpretations. T 1 represents the potential of the intermediate host population to deliver live schistosome to a given definitive host and T 2 represents the potential of the definitive host population to deliver live miracidia to a given uninfected snail. T 1 and T 2 are transmission factors. µ 1 and µ 2 are the instantaneous death rate per schistosome and per infected snail respectively. ν 11 is the instantaneous cercarial shedding rate per infected snail. ν 12 is the probability that a given free-swimming cercaria infects a given human host. ν 21 is the instantaneous rate of oviposition of a paired female schistosome. ν 22 is the probability that a fertilized egg gives rise to a miracidium which infects a given uninfected snail. All parameters were defined in [19], here we summarize them in Table 1.
According to [19], the function ψ(x) was defined as ψ(x) = x[1 − e −x (I 0 (x) + I 1 (x))] where I n denotes the modified Bessel function of the first kind of order n (n = 0, 1). The function F (a, b) is given by the equation m are the initial distributions of the number of male and female worms respectively in definitive host k (k = 1, 2, · · · , N 1 ). We can refer to [19] for more details about the definition and the properties of function G(a, b, f, m), ψ(x) and F (a, b).
As in most other models for schistosoma, in model equations (1), it is tacitly assumed that the transmission process requires no time. However, the actual situation in nature is that the time required for the mutual infection can not be neglected [12,24]. Various factors such as the latent periods during which cercariae and schistosome develop to their mature forms, the time required for miracidium instantaneous death rate per schistosome ν 11 instantaneous cercarial shedding rate per infected snail ν 12 probability that a given free-swimming cercaria infects a given human host µ 2 instantaneous death rate per infected snail ν 21 instantaneous rate of oviposition of a paired female schistosome ν 22 probability that a fertilized egg gives rise to a miracidium which infects a uninfected snail to find a snail, the latent periods between the moment when a snail is infected by a miracidium and the moment it begins shedding cercariae and the time required for cercariae to find an appropriate definitive host, cause a time lag for the transmission of infection. Therefore the model should be a system of differential equations with delays. We assume that the transit time for the transmission process -a mated female schistosomiasis produces eggs, eggs are voided with feces or urine and deposited in fresh water, the miracidia developed from the eggs and a snail is infected by the miracidia -is τ . We also assume that the transit time for the transmission process -the snails shed cercariae, the cercariae find an appropriate definitive, the juvenile schistosome migrate to liver, mature and mate-isτ . Then a modified model from (1) reads where all parameters are the same as in the systems (1) except for the positive constants τ andτ which represents the related delays respectively. In this paper, we first study the case whenτ = 0, then the model to be considered becomes We will then come back to comment on the case ofτ > 0 in Section 6. We shall consider system (3) with the initial conditions Since the expected worm burden of the definitive host cannot exceed the potential of the snail population to deliver live schistosome [19], we assume First we summarize and state some of the basic properties of the model (3). We adopt the following notations: For b > a, we denote C([a, b], R n ) the Banach space of continuous functions mapping the interval [a, b] into R n with the topology Assume Ω is a subset of R × C, f : Ω → R n is a given function, and " · " represents the right-hand derivative; then any retarded delay differential equation can be written We will need the following lemmas [13,14,16].
We now give a proposition without proof.
According to Lemma 2.1, Proposition 1 implies the local existence and uniqueness of the solution of initial problem (3) and (4).
Proof. For interior points the lemma is obvious. For boundary points the result is an immediate consequence of the following relations given in the table: there is a unique solution (β, y) of initial value problem (3) and (4) defined for all t ≥ t 0 and satisfy the initial conditions Moreover, for all t ≥ t 0 , (β(t), y(t)) ∈ R.
Proof. Proposition 1 insures a unique local solution. As f : Ω → R n is completely continuous and Proposition 2 shows that for the initial condition (4) the solution remains in R, the solution can be continued for all t ≥ t 0 . Global uniqueness is a consequence of local uniqueness.
3. Equilibria and their stability. In this section, we discuss the existence and local stability of the equilibria of system (3). In stead of working on system (3) directly, we shall study the following equivalent autonomous system subject to initial conditions Let D denote the rectangular parallelepiped It follows from Theorem 2.3 that D is positively invariant with respect to the system (6). Accordingly we will restrict the stability analysis of the equilibrium on D.
The equilibria of (6) are solutions of: Here we first summarize the discussion about the existence of equilibria with (8) from [19]. For the discussion, we let R 0 = T2 T20(T1) . Based on the analysis in [19], we can find R 0 is a threshold parameter for the existence of the equilibrium. One can see that the role of the parameter R 0 is similar to the so-called basic-reproduction number defined in [6] for disease modelling. Regarding the existence of equilibria in terms of the threshold value R 0 , we have the following proposition [19].
is the only equilibria of (6).
Proposition 4. The disease-free equilibrium E 0 is locally asymptotically stable for all τ ≥ 0.

3.1.
Stability of E s . Now we consider the case R 0 > 1. To study the stability of the endemic equilibrium E * (E s and E u ), we start with the characteristic equation of system (6) at the point E * .
When τ = 0, we obtain the characteristic equation of the corresponding ODE A straightforward calculation from (12) yields that a 1 > 0 and a 2 + a 3 > 0.
As a result of Hurwitz criterion, all eigenvalues of the characteristic equation at E s have negative real parts, then the endemic equilibrium E s is locally asymptotically stable when τ = 0. For τ > 0, according to Lemma 3.1, as τ > 0 varies, the sum of the multiplicities of roots of (11) in the open right half-plane can change only if a root appears on or crosses the imaginary axis. Due to condition (3.1), λ = 0 is not a root of (11), then the imaginary axis cannot be crossed by λ(τ ) = 0 for any τ > 0, and due to the Lemma 3.1 a stability switch (or a cross of the imaginary axis) necessarily occurs with λ = ±iω with ω > 0 (As P (−iω) = P (−iω) Q(−iω) = Q(iω), if λ = iω for some real w > 0 is a characteristic root of (11), then also λ = −iω is a characteristic root).
Without loss of generality we assume that λ = iω, ω > 0 is a root of (11). This is the case if and only if ω satisfies the following equation: Separating the real and imaginary parts, we have the following system for ω: a 3 cos(ωτ ) = ω 2 − a 2 , a 3 sin(ωτ ) = a 1 ω.
Eliminating the trigonometric functions from (16), we obtain the following forth degree equation for ω: Since this equation contains only even powers of ω, we can reduce the order by letting z = ω 2 . Then Eq.(17) becomes a second order equation of z: One can verify that Then (18) has no positive real roots. This implies that there is no ω such that iω is an eigenvalue of the characteristic equation (11). In other words, there are no roots of (11) appearing on or crossing the imaginary axis when τ is increased. Therefore, the real parts of all the eigenvalues of (11) are negative for all τ ≥ 0. Hence there are no stability switches as τ varies. Summarizing the above analysis, we have the following theorem: Theorem 3.2. Suppose that R 0 > 1 is satisfied, then the endemic equilibrium E s of the delay model (6) is locally asymptotically stable for all τ ≥ 0.

Stability of E u and Hopf bifurcation.
In the following, we analyze the local stability of the endemic equilibrium E u . Let τ = 0 in (11), we again obtain the characteristic equation (14) of the corresponding ODE model. By calculating the coefficients of (14) at the endemic equilibrium E u , if we still use the same notations for the coefficients a j (j = 1, 2, 3) respectively as for E s , we obtain: Consequently, (14) has a positive real solution and the endemic equilibrium E u is unstable when τ = 0. Now consider the case when τ > 0. If we once again assume that λ = iω with ω > 0 a root of (11), we will also get (16). Due to the fact that a 2 1 − 2a 2 > 0, and from (19) it is readily verified that a 2 2 − a 2 3 < 0. Thus (18) has only one positive root, say z 0 , and it is a simple root. Therefore, (17) has only one positive root, denoted by ω 0 . From (16), we define Then the corresponding characteristic equation of E u (11) has a pair of purely imaginary roots ±iω 0 with τ = τ j , j = 1, 2, . . .. Furthermore, one can prove the following proposition.

Proposition 5.
Let h(z) = z 2 + (a 2 1 − 2a 2 )z + a 2 2 − a 2 3 , then dh(z) dz | z=ω 2 0 > 0. If we consider the root λ of (11) as a function of τ , then we can write where η(τ j ) = 0 and ω(τ j ) = ω 0 ( j = 1, 2 . . .). We try to determine the direction of motion of λ as τ is increased. As such, we need to determine the sign Since the left side of (11) is an analytic function of λ and τ , a root λ will be a differentiable function of τ , except at points where the root is multiple. Then we may consider λ = λ(τ ) to be a differentiable function, and if we differentiate (11) with respect to τ , we get
The above analysis can be summarized into the following theorem.

4.
Saddle-node bifurcation. In this section, we focus on the local stability and bifurcation analysis of the endemic equilibrium E c .
Proposition 6. For the endemic equilibrium E c , λ = 0 is a simple eigenvalue, and all other eigenvalues have negative real part for all τ > 0.
We argue Proposition 6 by contradiction and the proof is omitted. One can carry out the bifurcation and local stability of the endemic equilibrium E c using the center manifold theory introduced by Carr [1] and the normal form method due to Faria [7], but based on the Proposition 3 and Proposition 6, we can easily obtain that the saddle-node bifurcation occurs with the parameter R 0 = 1. Instead of the standard analysis of the center mainfold and normal form calculation, taking advantage of the stability of the equilibrium of E s , E u and the property of saddle-node bifurcation, we can obtain the stability of equilibrium of E c immediately.   Figure 1. Saddle-node bifurcation diagram using R 0 as parameter.
Recall that T 2 represents the potential of the definitive host population to deliver live miracidia to a given uninfected snail. From the definition of T 20 (·), it can be interpreted as the minimum level of the capacity of the definitive host population delivering live miracidia to a given uninfected snail such that the disease will prevail for a give T 1 , then we can see R 0 reflect the endemic level of the infection to some extent. Now we regard β 1 (T 1 , T 2 ) and β 2 (T 1 , T 2 ) as functions of T 1 and T 2 . A straightforward calculation yields Hence if we fix T 1 , then R 0 as a function of T 2 , there holds According to the above analysis, as shown in Figure 1, we plot the illustrative bifurcation curve in the R 0 −β plane. In the figure the arrows along the vertical lines represent the flow generated by (6) along the β-direction. For R 0 > 0, β = 0 always exists and is always stable for any values of R 0 . For R 0 > 1, there are two endemic equilibria, E s and E u bifurcating from E c . The stability analysis shows that E s corresponding to β = β 2 (T 1 , T 2 ), is stable (represented by the solid curve), and the other equilibrium E u corresponding to β = β 1 (T 1 , T 2 ), is unstable (represented by the broken curve).
Note that if T 2 is fixed and if we regard R 0 as a function of T 1 , we can have the analogous results. 5. Global stability analysis. In this section, we discuss the global asymptotic stability of the disease-free equilibrium E 0 of system (6). The strategy of the proof is to construct an appropriate Lyapunov functional.
Proof. We show that the equilibrium E 0 is global attractive with respect to domain Consider the positive definite Lyapunov functional where c > 0 will be chosen later. Differentiating along solutions of (6) yields where ξ ∈ (0, α) and M = max 0≤α≤1,0≤β≤T1 F α (α, β). Hence we obtaiṅ The inequality (21) implies that Choose c ≥ T2µ2M  [16] the global asymptotic stability of the disease-free equilibrium E 0 in D follows. This proves the theorem.
The expression 1 2 ψ(T 1 ) represents the maximum level of the expected worm pair burden in a definitive host initially uninfected, and condition 1 2 ψ(T 1 )T 2 ≤ 1 can be interpreted as that the maximum level of the two interrelated flows of the transmission of schistosomiasis should be under a certain level. 6. Simulation. Due to the complexity and abstraction of the function F (a, b), we consider a relative easier case [19]. As we know, p  F (a, b) are the initial distributions of the number of male and female worms respectively in definitive host k (k = 1, 2, · · · , N 1 ), so we let the initial distribution of the numbers of male and female worms follow the Poisson distribution with the parameter 1 2 ω 0 , i.e., p (k) n q (k) n have the following form Let then the function ω(t) and y(t) satisfy the differential equations where ω(t) represents the expected worm burden per definitive host at time t, y(t) represents the expected proportion of infected intermediate hosts at time t. Refer to [19] for more details. We choose the parameters as the following: N 1 = 1000, N 2 = 5000, µ 1 = 0.05, ν 11 = 0.1, ν 12 = 0.0001, µ 2 = 0.1, ν 21 = 0.3, ν 22 = 0.01. Then we have T 1 = 5, T 2 = 30, R 0 = 30/1.3743 > 1, and the system (22) has two positive equilibrium E u = (0.0271776, 0.0054355), E s = (4.8792352, 0.9794578). ±iω 0 = ±0.0643i are corresponding purely imaginary roots of the characteristic equation at the positive equilibrium E u when τ j = −23.0042 + 97.7795j(j = 1, 2, . . .) which crosses the imaginary axis from left to right transversally.According to Theorem 3.3, if R 0 > 1 system (22) undergoes Hopf bifurcation when τ = τ j (j = 1, 2, · · · ), and the bifurcated periodic solution from E u is unstable. Due to its instability, it is difficult to catch this limit cycle for the delay differential equations by simulation, but on the other hand, as the solution of system (22) is the continuously dependent on the initial value, we could observe the oscillation in the neighborhood of the periodic solution. See Figure 2.
In fact from the above simulations one can see that the solution undergoing transient oscillation with initial value (ω(t), y(t)) = (0.0267302, 0.0066452), t ∈ [−74.7853, 0] goes up rapidly after t = 4700, and will converge to the stable equilibrium E s = (4.8792352, 0.9794578) . Furthermore, if we change the initial value a little with (ω(t), y(t)) = (0.0267301, 0.0066452), t ∈ [−74.7853, 0], the solution will eventually converge to another stable equilibrium E 0 = (0, 0) which also goes through a similar transient oscillation as the previous case. 7. Discussion. Based on the modeling studies by Nåsell and Hirsch [19], in this paper we modify and propose a system of delay differential equations (3) to study the transmission dynamics of schistosomiasis. The new model takes into account the time lags for the transmission of infection. Mathematical properties of the model are analyzed in terms of the stability of the steady states. We also carry out the bifurcation analysis, including Hopf bifurcation and saddle-node bifurcation. The threshold parameter R 0 is introduced and shown to determine the existence of the steady states. Here the role of R 0 is similar to the so called basic reproduction number [6] reflecting the infection power, a popular concept in the modeling studies of infectious diseases. Usually, one would expect that the disease free equilibrium be local unstable when the infection power, or correspondingly the reproduction number is great than one [6,21]. In this model, we have just the opposite phenomena, the "disease free equilibrium", the boundary equilibrium is always stable no matter the threshold value R 0 is smaller or greater than one. One possible explanation for this to happen is due to the fact that the mating heterosexually of two adult forms of schistosome in the liver of the definitive host.
Indeed, in a study by Cosgrove and Southgate [5], the authors revealed that schistosome mating structure could be very complicated, and one species can have the competitive advantages over another one. A recent modeling study involving the mating structure of schistosomiasis was carried by Castillo-Chavez et al [2]. In the study the authors proposed a system of homogeneous equations with a time delay to model the population dynamics of schistosomes. Though the parasites mating structure is included, the saddle-node bifurcation was not discovered in the model. We realize that in our model, also in the model in [19], the age-specific and sex-specific infection rate in the definitive host population were ignored. It remains an open problem to figure out the reason for the saddle-node bifurcation through a comparison study for the models with mating structure.
We also observe that a Hopf bifurcation occurs in the model when the threshold value increases through one at certain values of the delay. The existence of Hopf bifurcation suggests that the periodic phenomena, or the recurrence of the schistosomiasis is possible when R 0 > 1 and the initial infection density is in certain range.
Hopf bifurcation is also found in some other schistosomiasis models [8,26], and usually the bifurcated periodic solutions are stable. Contrary to those cases, unstable periodic solutions are found in our model through which we can interpret some phenomenon in the transmission of schistosomiasis to some extent. As the endemic equilibrium E u corresponding to the case that the disease is prevalent in low level is unstable, the disease will ultimately die out or break out when the initial state is closed to the E u . On the other hand, due to the continuous dependence of solutions on the initial value, it follows that any solution that starts in a neighborhood of these unstable periodic solutions would have transient oscillatory behavior. This may account for the recurrent phenomenon that happened in China or other places: the disease is controlled to a low level or to the instance where it nearly dies out by effective measures. The disease prevails in a very small scale and might even present oscillation in a relatively long period such as thirty or forty years, and finally the resurgence of the schistosomiasis outbreak in a large scale as time goes through [27]. For the delay model (2) withτ = 0, one can verify from the characteristic equation (9) that the local analysis are the same, but the global dynamics could be different. We keep it as a future project. Our results indicate that the stability of the equilibrium of the model is not sensitive to the delay. This coincides with the theoretical analysis of Wu [24] and Nåsell [18] by introducing a transition rate to include the latency in snail infection in MacDonald's model.
The threshold parameter R 0 = T2 T20(T1) allows us to gain some insight into the transmission dynamics of infection and to determine the role of various parameters in the disease dynamics. According to the result of our model, there exists a one-dimensional unstable manifold for R 0 > 1 and we still can not determine the attractive region of the disease free equilibrium. From the viewpoint of biology it is difficult to control the disease in this case, So effective intervention towards disease control can be achieved by reducing R 0 . Since all the environmental and biological factors are included in T 1 and T 2 , one may try to reduce T 1 and T 2 such that R 0 is below one and 1 2 ψ(T 1 )T 2 ≤ 1, therefore the infection will die out. In this paper, we only consider an idealized focus of infection: we restrict our discussion to schistosoma mansoni and schistosoma haematobium in endemic foci of infection where human is the only significant vertebrate host. We ignore age-specific and sex-specific infection rate in the definitive host population and the effect of human acquired immunity. More insights may be obtained if these assumptions can be relaxed for which we leave for future studies.