1. Introduction
Schistosomiasis is one of the most prevalent parasitic diseases. In 2004, WHO suggested that there were 200 million individuals infected worldwide [1] . However, today, more than 207 million people, 85% of whom live in Africa, are infected with schistosomiasis [2] , and the estimated 700 million people are at risk of infection in 76 countries where the disease is considered endemic, as their agricultural work, domestic chores, and recreational activities expose them to infested water [2] [3] . Globally, 200,000 deaths are attributed to schistosomiasis annually [4] . Thus, controlling schistosomiasis is a long-term task in the developing nations, and mathematical modeling of Schistosoma japonicum transmission is beneficial for the development of new strategies for control.
After the pioneering work of Macdonald [5] , the dynamics properties (including stability, persistence, and oscillatory behavior) of the schistosomiasis models that have significant biological background have been one of the most active areas of research and have attracted great attention of many researchers. Many excellent and interesting results have been obtained (see [6] - [18] ). Besides, a simplification of the two-strain, vector-host model was proposed by Feng and Velasco-Hernández [19] for Dengue fever. The model couples a simple XYX model for the hosts with an SI model for the vectors. The four compartments correspond to infected hosts (Y), infected vectors (I), susceptible hosts (X) and susceptible vectors (S). Hosts are infected by contacts with infected vectors, and vectors are in turn infected by contacts with infected hosts. These infection rates are given by the two terms
and
. The model is written as follows:
(1.1)
where
and
are the disease transmission coefficients;
is the recovery rate of infected host. The birth and death rates have been scaled to
for the host and
for the vector.
Note that an infected snail cannot infect susceptible man (or an animal) directly and vice versa. Schistosomiasis has a complicated life cycle involving two free living stages, the miracidia and the cercariae; and two host populations, the human and the snail. The parasite eggs hatch into free-swimming larva called miracidia in water; the miracidium then penetrates an appropriate snail at suitable temperature. In the infected snail, the miracidium undergoes asexual multiplication through a series of stages called sporocyst; then thousands of free-swimming cercariae are released. Cercariae are shed from the infected snail and penetrate the skin of a definitive host (such as human) within a few minutes after exposure and transform into schistosomula, which travel through the bloodstream to the liver, where they mature into adults and start producing eggs [20] . The eggs infiltrate through the tissues and are passed in the feces. That finishes schistosomiasis life cycle. Besides, the model (1.1) can be used to describe the transmission of schistosomiasis since schistosomiasis is a snail-vector disease.
It is known that there are prepatent periods of schistosoma. In fact, it is about five weeks from the time of cercaria penetration through skins of human hosts to the time when eggs are discharged [21] . That is, a susceptible host becomes infection for some time and then excretes feces with parasite eggs. It is easy to see that the prepatent period of hosts is very important for Schistosome japonicum transmission. Hence, it is necessary to study the impact of the prepatent period on schistosomiasis transmission. The aim of this paper is to incorporate the prepatent period of infected hosts into (1.1), and estimate the impact of the prepatent period on the schistosomiasis transmission. In this paper, we incorporate effects of the prepatent period of infected hosts into the model (1.1) and propose a schistosomiasis model as follows:
(1.2)
with initial conditions
,
,
,
and
,
,
, where
, the space of continuous functions mapping the interval
into
, where
.
In system (1.2), A and
are the recruitment rates of hosts and snails, respectively. The constant
is the per capita rate of infection of hosts by cercaria released by a infected snail,
is the per capita rate of infection of snails by miracidia from the parasite eggs from a infected host. The constant v is the recovery rate of infected host. Constants
and
represent the natural death rate and disease inducing death rate of hosts and snails, respectively.
is the prepatent period in host. We assume that all parameters are positive constants.
In fact, the reciprocal of the death rate,
, is equivalent to the life expectancy of human. Assume that
, then
. This means that the survival rate
is infinitely close to 1. Therefore, we make a simplification common in system (1.2) with the survival rate and assume that the survival rate has negligible impact on dynamics. Thus system (1.2) can be written in the following form:
(1.3)
The effects of time delays on the dynamical behaviors of schistosomiasis have been investigated in the literatures [21] [22] [23] [24] . For example, Liang et al. [21] investigated the development period of worms in human hosts, they described temperature-dependent and precipitation-dependent effects on snail abundance and infection as well as seasonal aspects of local agricultural practice. In [23] , a fixed delay was inspired by the life history of schistosomes, they investigated the impact of the delay on the invasion and persistence of drug-resistant parasite strains as well as on multi-strain coexistence. The main purpose of this paper is to study the effects of the time delay on the dynamical behaviors of (1.3), and discuss the direction of the bifurcation and stability of the bifurcating periodic solutions.
The remainder of the paper is organized as follows. In Section 2, we obtain the stability of disease-free equilibrium and the existence of the endemic equilibrium. In Section 3, by analyzing the characteristic equation of the linearized system of system (1.3) at the endemic equilibrium, we discuss the stability of the endemic equilibrium and the existence of the Hopf bifurcations occurring at the endemic equilibrium. In Section 4, by using the normal form theory and the Center Manifold Theorem, the formulae determining the direction of the Hopf bifurcations and the stability of bifurcating periodic solutions are obtained. Some numerical simulations are presented to illustrate our theoretical results in Section 5. This paper ends with a brief conclusion.
2. Equilibrium Analysis
In this section, we discuss the existence of equilibria and the stability of the disease free equilibrium.
When the infective hosts and the infective snails do not exist, i.e.,
, then
and
. This is the infection free equilibrium
for schistosomiasis. The following theorem determines stability of
and existence of endemic equilibrium in terms of a threshold parameter
Theorem 2.1. If
, then system (1.3) has a unique equilibrium
and
is stable if
. If
, then system (1.3) has an endemic equilibrium
except the disease free equilibrium
, where
(2.1)
Proof Computing the nonnegative solutions of the following equations:
(2.2)
we can easily obtain the existence of two equilibria
and
.
Next, we show the stability behavior of equilibrium
by finding the eigenvalues of the corresponding Jacobian matrix obtained for system (1.3).
The Jacobian matrix for system (1.3) is as follows:
(2.3)
Let
be the Jacobian matrix J evaluated at the equilibrium
. From
, it is easy to calculate the associated characteristic equation of system (1.3) at
and obtain
(2.4)
It is obvious that
are negative characteristic roots of (2.4). Hence, we only need to discuss the roots of the following equation:
(2.5)
Note that all characteristic roots of (2.5) are negative if
. This yields that all roots of (2.4) are negative if
. We complete the proof.
3. Endemic Equilibrium and Hopf Bifurcation
In this section, We investigate the stability of the endemic equilibrium
and existence of the Hopf bifurcation occurring at
.
Similar to Section 2, let
be the Jacobian matrix (2.3) evaluated at the equilibrium
. Then we calculate the characteristic equation of system (1.3) at
and obtain
(3.1)
where
and
.
In the above expression of
and
,
and
are given as follows:
Taking
, we can rewrite (3.1) as
(3.2)
where
Now it is easy to see that
. Thus for the local stability of the endemic equilibrium
of system (1.3) without delay, we have the following result.
Theorem 3.1. When
, the endemic equilibrium
is locally asymptotically stable for
if the following conditions are satisfied:
where
and
are defined as above.
Then we turn to an investigation of local stability of the endemic equilibrium
.
We know all roots of characteristic Equation (3.1) have negative real parts at
when the conditions in Theorem 3.1 are satisfied. Next we will show that there is a unique pair of purely imaginary roots
for characteristic Equation (3.1).
Assume that for some
,
is a root of (3.1), which implies
(3.3)
Separating real and imaginary parts, we get the following equations:
(3.4)
Now squaring and adding Equations (3.4), we get
(3.5)
where
Substituting
in above Equation (3.5), we have
(3.6)
Now if the coefficients in
satisfy the conditions of the Routh-Hurwitz criterion, then Equation (3.6) will not have any positive real root, thus we may not get any positive value of
which satisfies the Equation (3.5). In this case the result may be written in the form of following theorem.
Theorem 3.2. Assume that the coefficients in
defined in (3.6) satisfy the conditions of the Routh-Hurwitz criterion, then the endemic equilibrium
of system (1.3) is asymptotically stable for all delay
if it is stable in the absence of delay.
Assuming contrary that the values of
in (3.6) do not satisfy the Routh-Hurwitz criterion. In this case a simple assumption for the existence of a positive root of Equation (3.6) is
, which implies
(3.7)
Now if condition (3.7) holds, then Equation (3.6) has a positive root
and thus Equation (3.5) has a pair of purely imaginary roots
. It follows from Equations (3.4) that
Then
corresponding to this positive value of
is given as follows:
By using Butler’s Lemma, we can say that the endemic equilibrium
remains stable for
.
Next we investigate whether there is a phenomenon of Hopf bifurcation as
increases through
. For this the following lemma is needed.
Lemma 3.1. The following transversality condition is satisfied:
(3.8)
provided that condition (3.7) holds.
Proof Differentiating Equation (3.1) with respect to
, we get
(3.9)
where
Therefore,
From the above argument, we know
, then
(3.10)
Here it may be noted that
if the condition (3.7) is satisfied. This proves the Lemma 3.1. Thus we have the following result:
Theorem 3.3. If
and the condition (3.7) hold, then the endemic equilibrium
of system (1.3) remains stable for all
and becomes unstable for
. System (1.3) with
undergoes a Hopf bifurcation.
Remark. It must be pointed out that Theorem 3.3 cannot determine the stability and direction of bifurcation periodic solutions. That is to say, the periodic solutions may exist for
near
. Next, in Section 5 the stability of bifurcating periodic solutions is investigated by analyzing higher order terms according to Hassard et al. [25] .
4. Stability and Direction of Hopf Bifurcation
In this section, in term of the center manifold and normal form theory due to Hassard et al. [25] , the direction of hopf bifurcation and the stability of periodic bifurcation solution are discussed.
Without loss of generality, and
. So,
is the Hopf value of system (1.3).
Let
,
,
,
,
, and dropping the bars for simplification of notations, system (1.3) becomes a functional differential equation in
as
(4.1)
where
, and
,
are given, respectively, by
(4.2)
and
(4.3)
By the Riesz representation theorem, there exists a function
of bounded variation for
such that
(4.4)
In fact, we can choose
(4.5)
where
denote the Dirac delta function:
For
, define
and
Then system (1.3) is equivalent to
(4.6)
where
for
.
For
, define
and a bilinear inner product
(4.7)
where
. Then
and
are adjoint operators. By the above discussion, we know that
are eigenvalues of
. Hence, they are also eigenvalues of
. We first need to compute the eigenvectors of
and
corresponding to
and
, respectively.
Supposed
is the eigenvectors of
corresponding to
, then
. Then from the definition of
and (4.2), (4.4) and (4.5), we have
For
, then we obtain
Similarly, we can obtain the eigenvector
of
corresponding to
, where
In order to assure
, we need to determine the value of D. By (4.7), we have
Therefore, we can choose D as
Next we will compute the coordinate to describe the center manifold
at
. Let
be the solution of (4.6) when
. Define
(4.8)
On the center manifold
, we have
(4.9)
where z and
are local coordinates for center manifold
in the direction of
and
. Note that W is real if
is real. We only consider real solutions. For solution
of (4.6), since
, we obtain
(4.10)
where
(4.11)
It follows from (4.8) and (4.9) that
and
, then
Then from the definition of
, we obtain
(4.12)
Comparing the coefficients with (4.11), we have
In order to assure
,
and
are needed to compute.
From (4.6), (4.8) and (4.10), we have
(4.13)
where
(4.14)
It follows from (4.13) and (4.14) that
(4.15)
From (4.9) and (4.10), we have
(4.16)
and
(4.17)
Substituting (4.16) and (4.17) into (4.15) and comparing the coefficients of
and
, we have
(4.18)
From (4.11) and (4.13), we know that for
(4.19)
Comparing the coefficients with (4.19) gives that
(4.20)
and
(4.21)
From the definition of
and (4.18) and (4.21), we obtain
For
, we have
(4.22)
where
is a constant vector. Similarly, from (4.18) and (4.22), we know
(4.23)
where
is a constant vector.
Finally, we will seek the values of
and
. From the definition of
and (4.18), we have
(4.24)
and
(4.25)
where
. By (4.13), we know when
That is
(4.26)
By (4.3), we have
By (4.8), we obtain
Then, we have
(4.27)
By (4.26) and (4.27), we have
(4.28)
and
(4.29)
For
is the eigenvalue of
and
is the corresponding eigenvector, we obtain
and
So, substituting (4.22) and (4.28) into (4.24), we obtain
That is
It follows that
Similarly, substituting (4.23) and (4.29) into (4.25), we also get
and
Thus, we can determine
and
from (4.22) and (4.23). Further, we can compute
. Thus we can compute the following values:
which determine the qualities of bifurcation periodic solution in the center manifold at the critical values
.
According to [25] , we can obtain the following result.
Theorem 4.1. Assume that
and the condition (3.7) hold, we have:
1) if
(
), then the Hopf bifurcation is supercritical (subcritical) and the bifurcating periodic solutions exist for
.
2) if
(
), then the bifurcating periodic solutions are stable (unstable).
3) if
(<0), then the period of the bifurcating periodic solutions increases (decreases).
From Theorem 4.1, we know that the value of
determines the directions of the Hopf bifurcation, the values of
and
determine the stability and the period of the bifurcating periodic solutions, respectively.
5. A Numerical Example
In this section, we implement numerical simulations to testify the above theoretical results. Let
,
,
,
,
,
,
,
,
. For
, using the software Matlab, we derive
. Thus the endemic equilibrium
is stable when
. Figures 1(a)-(f) show that the endemic equilibrium
is asymptotically stable when
. When
passes through the critical value
, the endemic equilibrium
loses its stability and a Hopf bifurcation occurs, that is, a family of periodic solutions bifurcate from the endemic equilibrium
. Figures 2(a)-(f) suggest that Hopf bifurcation occurs from the endemic equilibrium
when
.
6. Conclusion
In this paper, we have investigated a delayed schistomiasis model, and studied
the local stability of the equilibria and Hopf bifurcation. We have shown that if
, the disease-free equilibrium is locally asymptotically stable. Further, the sufficient conditions for the stability of the endemic equilibrium are obtained. That is, if
and the condition (3.7) hold, the endemic equilibrium
is asymptotically stable for all
. As
increases, the endemic equilibrium loses its stability and a sequence of Hopf bifurcations occurs at the endemic equilibrium; that is, urcates from the equilibrium. This shows that the density of thea family of periodic orbits bif susceptible human, snails and the infected human, snails may keep in an oscillatory mode near the endemic equilibrium. By the normal form theory and the Center Manifold Theorem, the direction of Hopf bifurcation and the stability of the bifurcating periodic orbits have been discussed.
Figure 2. (a)-(f) The dynamical behavior of system (1.3) with
. Hopf bifurcation occurs from the endemic equilibrium
.
Acknowledgements
This work is supported by The National Natural Science Foundation of China (No.11561004), The Science and Technology research project of Jiangxi Provincial Education Department (No.GJJ170815) and The bidding project of Gannan Normal University (No.16zb02).