Dynamical Analysis of a Schistosomiasis japonicum Model with Time Delay

Fumin Zhang, Shujing Gao, Huahua Cao, Youquan Luo

Key Laboratory of Jiangxi Province for Numerical Simulation and Emulation Techniques, Gannan Normal University, Ganzhou, China.

DOI: 10.4236/jamp.2019.74064                

Abstract

In this paper, a Schistosomiasis japonicum model incorporating time delay is proposed which represents the developmental time from cercaria penetration through skins of human hosts to egg laying. By linearizing the system at the positive equilibrium and analyzing the associated characteristic equations, the local stability of the equilibria is investigated. And it proves that Hopf bifurcations occur when the time delay passes through a sequence of critical value. Furthermore, the explicit formulae for determining the stability and the direction of the Hopf bifurcation periodic solutions are derived by using techniques from the normal form theory and Center Manifold Theorem. Some numerical simulations which support our theoretical analysis are also conducted.

Keywords

Schistosomiasis japonicum, Time Delay, Stability, Hopf Bifurcation

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 ${\beta }_{1}XI$ and ${\beta }_{2}YS$. The model is written as follows:

$\{\begin{array}{l}\frac{\text{d}X}{\text{d}t}=b-{\beta }_{1}X\left(t\right)I\left(t\right)-bX\left(t\right)+vY\left(t\right),\hfill \\ \frac{\text{d}Y}{\text{d}t}={\beta }_{1}X\left(t\right)I\left(t\right)-\left(v+b\right)Y\left(t\right),\hfill \\ \frac{\text{d}S}{\text{d}t}=c-{\beta }_{2}Y\left(t\right)S\left(t\right)-cS\left(t\right),\hfill \\ \frac{\text{d}I}{\text{d}t}={\beta }_{2}Y\left(t\right)S\left(t\right)-cI\left(t\right),\hfill \end{array}$ (1.1)

where ${\beta }_1>0$ and ${\beta }_2>0$ are the disease transmission coefficients; $v>0$ is the recovery rate of infected host. The birth and death rates have been scaled to $b>0$ for the host and $c>0$ 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:

$\{\begin{array}{l}\frac{\text{d}X}{\text{d}t}=A-{\beta }_{1}X\left(t-\tau \right)I\left(t\right){\text{e}}^{-d_1\tau }-d_{1}X\left(t\right)+vY\left(t\right),\hfill \\ \frac{\text{d}Y}{\text{d}t}={\beta }_{1}X\left(t-\tau \right)I\left(t\right){\text{e}}^{-d_1\tau }-\left(v+{\alpha }_1+d_1\right)Y\left(t\right),\hfill \\ \frac{\text{d}S}{\text{d}t}=\Lambda -{\beta }_{2}Y\left(t\right)S\left(t\right)-d_{2}S\left(t\right),\hfill \\ \frac{\text{d}I}{\text{d}t}={\beta }_{2}Y\left(t\right)S\left(t\right)-\left({\alpha }_2+d_2\right)I\left(t\right),\hfill \end{array}$ (1.2)

with initial conditions $X\left(\theta \right)={\varphi }_1\left(\theta \right)$, $Y\left(\theta \right)={\varphi }_2\left(\theta \right)$, $S\left(\theta \right)={\varphi }_3\left(\theta \right)$, $I\left(\theta \right)={\varphi }_4\left(\theta \right)$ and ${\varphi }_i\left(\theta \right)\ge 0$, $\theta \in \left[-\tau \mathrm{,0}\right]$, ${\varphi }_i\left(0\right)>0\left(i=1,2,3,4\right)$, where $\left({\varphi }_1\left(\theta \right),{\varphi }_2\left(\theta \right),{\varphi }_3\left(\theta \right),{\varphi }_3\left(\theta \right)\right)\in C\left(\left[-\tau ,0\right],{ℝ}_{+0}^4\right)$, the space of continuous functions mapping the interval $\left[-\tau \mathrm{,0}\right]$ into ${ℝ}_{+0}^4$, where ${ℝ}_{+0}^4=\left\{\left(x_1,x_2,x_3,x_4\right):x_i\ge 0,i=1,2,3,4\right\}$.

In system (1.2), A and $\Lambda$ are the recruitment rates of hosts and snails, respectively. The constant ${\beta }_1$ is the per capita rate of infection of hosts by cercaria released by a infected snail, ${\beta }_2$ 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 $d_i$ and ${\alpha }_i\left(i=1,2\right)$ represent the natural death rate and disease inducing death rate of hosts and snails, respectively. $\tau$ is the prepatent period in host. We assume that all parameters are positive constants.

In fact, the reciprocal of the death rate, $1/d_1$, is equivalent to the life expectancy of human. Assume that $1/d_1=75\times 365\text{days}=27375\text{days}$, then $d_1=3.65\times {10}^{-5}$. This means that the survival rate ${\text{e}}^{-d_1\tau }$ 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:

$\{\begin{array}{l}\frac{\text{d}X}{\text{d}t}=A-{\beta }_{1}X\left(t-\tau \right)I\left(t\right)-d_{1}X\left(t\right)+vY\left(t\right),\hfill \\ \frac{\text{d}Y}{\text{d}t}={\beta }_{1}X\left(t-\tau \right)I\left(t\right)-\left(v+{\alpha }_1+d_1\right)Y\left(t\right),\hfill \\ \frac{\text{d}S}{\text{d}t}=\Lambda -{\beta }_{2}Y\left(t\right)S\left(t\right)-d_{2}S\left(t\right),\hfill \\ \frac{\text{d}I}{\text{d}t}={\beta }_{2}Y\left(t\right)S\left(t\right)-\left({\alpha }_2+d_2\right)I\left(t\right).\hfill \end{array}$ (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., $Y=I=0$, then $X=\frac{A}{d_1}$ and $S\mathrm{<mo>\; =\; </mo>}\frac{\Lambda }{d_2}$. This is the infection free equilibrium $E_0=\left(\frac{A}{d_1},0,\frac{\Lambda }{d_2},0\right)$ for schistosomiasis. The following theorem determines stability of $E_0$ and existence of endemic equilibrium in terms of a threshold parameter

$R_0=\frac{{\beta }_1{\beta }_2\Lambda A}{d_1{d}_2\left(v+{\alpha }_1+d_1\right)\left({\alpha }_2+d_2\right)}.$

Theorem 2.1. If $R_0<1$, then system (1.3) has a unique equilibrium $E_0=\left(\frac{A}{d_1},0,\frac{\Lambda }{d_2},0\right)$ and $E_0$ is stable if $R_0<1$. If $R_0>1$, then system (1.3) has an endemic equilibrium $E^{*}=\left(X^{*},Y^{*},S^{*},I^{*}\right)$ except the disease free equilibrium $E_0$, where

$\begin{array}{l}{X}^{*}=\frac{A}{d_1}-\frac{{\alpha }_1+d_1}{d_1}\cdot \frac{{\beta }_1{\beta }_2\Lambda A-d_1{d}_2\left(v+{\alpha }_1+d_1\right)\left({\alpha }_2+d_2\right)}{{\beta }_1{\beta }_2\Lambda \left({\alpha }_1+d_1\right)+d_1{\beta }_2\left(v+{\alpha }_1+d_1\right)\left({\alpha }_2+d_2\right)},\\ Y^{*}=\frac{{\beta }_1{\beta }_2\Lambda A-d_1{d}_2\left(v+{\alpha }_1+d_1\right)\left({\alpha }_2+d_2\right)}{{\beta }_1{\beta }_2\Lambda \left({\alpha }_1+d_1\right)+d{\beta }_2\left(v+{\alpha }_1+d_1\right)\left({\alpha }_2+d_2\right)},\\ S^{*}=\frac{\Lambda }{d_2}-\frac{{\alpha }_2+d_2}{d_2}\cdot \frac{{\beta }_1{\beta }_2\Lambda A-d_1{d}_2\left(v+{\alpha }_1+d_1\right)\left({\alpha }_2+d_2\right)}{{\beta }_1\left({\alpha }_2+d_2\right)\left[A{\beta }_2+d_2\left({\alpha }_1+d_1\right)\right]},\\ I^{*}=\frac{{\beta }_1{\beta }_2\Lambda A-d_1{d}_2\left(v+{\alpha }_1+d_1\right)\left({\alpha }_2+d_2\right)}{{\beta }_1\left({\alpha }_2+d_2\right)\left[A{\beta }_2+d_2\left({\alpha }_1+d_1\right)\right]}.\end{array}$ (2.1)

Proof Computing the nonnegative solutions of the following equations:

$\begin{array}{l}A-{\beta }_{1}X\left(t-\tau \right)I\left(t\right)-d_{1}X\left(t\right)+vY\left(t\right)=0,\\ {\beta }_{1}X\left(t-\tau \right)I\left(t\right)-\left(v+{\alpha }_1+d_1\right)Y\left(t\right)=0,\\ \Lambda -{\beta }_{2}Y\left(t\right)S\left(t\right)-d_{2}S\left(t\right)=0,\\ {\beta }_{2}Y\left(t\right)S\left(t\right)-\left({\alpha }_2+d_2\right)I\left(t\right)=0,\end{array}$ (2.2)

we can easily obtain the existence of two equilibria $E_0$ and $E^{\mathrm{<mo>\; *\; </mo>}}$.

Next, we show the stability behavior of equilibrium $E_0$ by finding the eigenvalues of the corresponding Jacobian matrix obtained for system (1.3).

The Jacobian matrix for system (1.3) is as follows:

$J=\left(\begin{array}{cccc}-d_1-{\beta }_{1}I\left(t\right){\text{e}}^{-\lambda \tau }& v& 0& -{\beta }_{1}X\left(t-\tau \right)\\ {\beta }_{1}I\left(t\right){\text{e}}^{-\lambda \tau }& -\left(v+{\alpha }_1+d_1\right)& 0& {\beta }_{1}X\left(t-\tau \right)\\ 0& -{\beta }_{2}S\left(t\right)& -{\beta }_{2}Y\left(t\right)-d_2& 0\\ 0& {\beta }_{2}S\left(t\right)& {\beta }_{2}Y\left(t\right)& -\left({\alpha }_2+d_2\right)\end{array}\right)$ (2.3)

Let $J\left(E_0\right)$ be the Jacobian matrix J evaluated at the equilibrium $E_0$. From $J\left(E_0\right)$, it is easy to calculate the associated characteristic equation of system (1.3) at $E_0$ and obtain

$\left(\lambda +d_1\right)\left(\lambda +d_2\right)\left[\left(\lambda +v+{\alpha }_1+d_1\right)\left(\lambda +{\alpha }_2+d_2\right)-{\beta }_1{\beta }_2\frac{A\Lambda }{d_1{d}_2}\right]=0.$ (2.4)

It is obvious that ${\lambda }_i=-d_i\left(i=1,2\right)$ are negative characteristic roots of (2.4). Hence, we only need to discuss the roots of the following equation:

$\left(\lambda +v+{\alpha }_1+d_1\right)\left(\lambda +{\alpha }_2+d_2\right)-{\beta }_1{\beta }_2\frac{A\Lambda }{d_1{d}_2}=0.$ (2.5)

Note that all characteristic roots of (2.5) are negative if $R_0<1$. This yields that all roots of (2.4) are negative if $R_0<1$. We complete the proof.

3. Endemic Equilibrium and Hopf Bifurcation

In this section, We investigate the stability of the endemic equilibrium $E^{\mathrm{<mo>\; *\; </mo>}}$ and existence of the Hopf bifurcation occurring at $E^{\mathrm{<mo>\; *\; </mo>}}$.

Similar to Section 2, let $J\left(E^{\mathrm{<mo>\; *\; </mo>}}\right)$ be the Jacobian matrix (2.3) evaluated at the equilibrium $E^{\mathrm{<mo>\; *\; </mo>}}$. Then we calculate the characteristic equation of system (1.3) at $E^{\mathrm{<mo>\; *\; </mo>}}$ and obtain

$P\left(\lambda \right)+Q\left(\lambda \right){\text{e}}^{-\lambda \tau }=0,$ (3.1)

where $P\left(\lambda \right)={\lambda }^4+A_1{\lambda }^3+A_2{\lambda }^2+A_3\lambda +A_4$ and $Q\left(\lambda \right)=B_1{\lambda }^3+B_2{\lambda }^2+B_3\lambda +B_4$.

In the above expression of $P\left(\lambda \right)$ and $Q\left(\lambda \right)$, $A_i\left(i=1,2,3,4\right)$ and $B_j\left(j=1,2,3,4\right)$ are given as follows:

$A_1={\beta }_2{Y}^{*}+2d_1+2d_2+{\alpha }_1+{\alpha }_2+v,$

$\begin{array}{c}{A}_2=\left({\beta }_2{Y}^{\mathrm{<mo>\; *\; </mo>}}+d_2\right)\left({\alpha }_2+d_2\right)+\left({\beta }_2{Y}^{\mathrm{<mo>\; *\; </mo>}}+2d_2+{\alpha }_2\right)\left(v+{\alpha }_1+2d_1\right)\\ +d_1\left(v+{\alpha }_1+d_1\right)-{\beta }_1{\beta }_2{X}^{\mathrm{<mo>\; *\; </mo>}}{S}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}\end{array}$

$\begin{array}{c}{A}_3=\left({\beta }_2{Y}^{\mathrm{<mo>\; *\; </mo>}}+d_2\right)\left({\alpha }_2+d_2\right)\left(v+{\alpha }_1+2d_1\right)\\ +d_1\left({\beta }_2{Y}^{\mathrm{<mo>\; *\; </mo>}}+2d_2+{\alpha }_2\right)\left(v+{\alpha }_1+d_1\right)-\left(d_1+d_2\right){\beta }_1{\beta }_2{X}^{\mathrm{<mo>\; *\; </mo>}}{S}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}\end{array}$

$A_4=d_1\left({\beta }_2{Y}^{*}+d_2\right)\left({\alpha }_2+d_2\right)\left(v+{\alpha }_1+d_1\right)-{\beta }_1{\beta }_2{d}_1{d}_2{X}^{*}{S}^{*},$

$B_1={\beta }_1{I}^{*},$

$B_2={\beta }_1{I}^{*}\left({\beta }_2{Y}^{*}+2d_2+d_1+{\alpha }_1+{\alpha }_2\right),$

$B_3={\beta }_1{I}^{*}\left[\left({\beta }_2{Y}^{*}+d_2\right)\left({\alpha }_2+d_2\right)+\left({\beta }_2{Y}^{*}+2d_2+{\alpha }_2\right)\left({\alpha }_1+d_1\right)\right],$

$B_4={\beta }_1{I}^{*}\left({\beta }_2{Y}^{*}+d_2\right)\left({\alpha }_1+d_1\right)\left({\alpha }_2+d_2\right).$

Taking $\tau =0$, we can rewrite (3.1) as

${\lambda }^4+a_1{\lambda }^3+a_2{\lambda }^2+a_3\lambda +a_4=0,$ (3.2)

where

$a_1={\beta }_2{Y}^{*}+2d_1+2d_2+{\alpha }_1+{\alpha }_2+v+{\beta }_1{I}^{*},$

$\begin{array}{c}{a}_2=\left({\beta }_2{Y}^{\mathrm{<mo>\; *\; </mo>}}+d_2\right)\left({\alpha }_2+d_2\right)+\left({\beta }_2{Y}^{\mathrm{<mo>\; *\; </mo>}}+2d_2+{\alpha }_2\right)\left(v+{\alpha }_1+2d_1+{\beta }_1{I}^{\mathrm{<mo>\; *\; </mo>}}\right)\\ +\left(d_1+{\beta }_1{I}^{\mathrm{<mo>\; *\; </mo>}}\right)\left(v+{\alpha }_1+d_1\right)-{\beta }_1{\beta }_2{X}^{\mathrm{<mo>\; *\; </mo>}}{S}^{\mathrm{<mo>\; *\; </mo>}}-v{\beta }_1{I}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}\end{array}$

$\begin{array}{c}{a}_3=\left({\beta }_2{Y}^{\mathrm{<mo>\; *\; </mo>}}+d_2\right)\left({\alpha }_2+d_2\right)\left(v+{\alpha }_1+2d_1+{\beta }_1{I}^{\mathrm{<mo>\; *\; </mo>}}\right)\\ +\left({\beta }_2{Y}^{\mathrm{<mo>\; *\; </mo>}}+2d_2+{\alpha }_2\right)\left[\left(d_1+{\beta }_1{I}^{\mathrm{<mo>\; *\; </mo>}}\right)\left(v+{\alpha }_1+d_1\right)-v{\beta }_1{I}^{\mathrm{<mo>\; *\; </mo>}}\right]\\ -\left(d_1+d_2\right){\beta }_1{\beta }_2{X}^{\mathrm{<mo>\; *\; </mo>}}{S}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}\end{array}$

$a_4=\left({\beta }_2{Y}^{\mathrm{<mo>\; *\; </mo>}}+d_2\right)\left({\alpha }_2+d_2\right)\left[\left(d_1+{\beta }_1{I}^{\mathrm{<mo>\; *\; </mo>}}\right)\left(v+{\alpha }_1+d_1\right)-v{\beta }_1{I}^{\mathrm{<mo>\; *\; </mo>}}\right]-{\beta }_1{\beta }_2{d}_1{d}_2{X}^{\mathrm{<mo>\; *\; </mo>}}{S}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{.}$

Now it is easy to see that $a_1>0$. Thus for the local stability of the endemic equilibrium $E^{\mathrm{<mo>\; *\; </mo>}}\left(X^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{Y}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{S}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{I}^{\mathrm{<mo>\; *\; </mo>}}\right)$ of system (1.3) without delay, we have the following result.

Theorem 3.1. When $R_0>1$, the endemic equilibrium $E^{\mathrm{<mo>\; *\; </mo>}}\left(X^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{Y}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{S}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{I}^{\mathrm{<mo>\; *\; </mo>}}\right)$ is locally asymptotically stable for $\tau =0$ if the following conditions are satisfied:

$a_4>0,a_1{a}_2-a_3>0,a_3\left(a_1{a}_2-a_3\right)-a_1^2{a}_4>0,$

where $a_1\mathrm{,}{a}_2\mathrm{,}{a}_3$ and $a_4$ are defined as above.

Then we turn to an investigation of local stability of the endemic equilibrium $E^{\mathrm{<mo>\; *\; </mo>}}$.

We know all roots of characteristic Equation (3.1) have negative real parts at $\tau =0$ when the conditions in Theorem 3.1 are satisfied. Next we will show that there is a unique pair of purely imaginary roots $\pm i{\omega }_0\left({\omega }_0>0\right)$ for characteristic Equation (3.1).

Assume that for some $\tau >0$, $i\omega \left(\omega >0\right)$ is a root of (3.1), which implies

${\omega }^4-A_1{\omega }^{3}i-A_2{\omega }^2+A_3\omega i+A_4+\left(-B_1{\omega }^{3}i-B_2{\omega }^2+B_3\omega i+B_4\right){\text{e}}^{-i\omega \tau }=0.$ (3.3)

Separating real and imaginary parts, we get the following equations:

$\{\begin{array}{l}{\omega }^4-A_2{\omega }^2+A_4=\left(B_1{\omega }^3-B_3\omega \right)\sin \omega \tau +\left(B_2{\omega }^2-B_4\right)\cos \omega \tau \mathrm{,}\\ -A_1{\omega }^3+A_3\omega =\left(B_1{\omega }^3-B_3\omega \right)\cos \omega \tau -\left(B_2{\omega }^2-B_4\right)\sin \omega \tau \mathrm{.}\end{array}$ (3.4)

Now squaring and adding Equations (3.4), we get

${\omega }^8+D_1{\omega }^6+D_2{\omega }^4+D_3{\omega }^2+D_4=0,$ (3.5)

where

$D_1=A_1^2-2A_2-B_1^2,$

$D_2=A_2^2+2A_4-2A_1{A}_3+2B_1{B}_3-B_2^2,$

$D_3=A_3^2-2A_2{A}_4-B_3^2+2B_2{B}_4,$

$D_4=A_4^2-B_4^2.$

Substituting ${\omega }^2=\eta$ in above Equation (3.5), we have

$f\left(\eta \right)={\eta }^4+D_1{\eta }^3+D_2{\eta }^2+D_3\eta +D_4=0.$ (3.6)

Now if the coefficients in $f\left(\eta \right)$ 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 $\omega$ 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 $f\left(\eta \right)$ defined in (3.6) satisfy the conditions of the Routh-Hurwitz criterion, then the endemic equilibrium $E^{\mathrm{<mo>\; *\; </mo>}}$ of system (1.3) is asymptotically stable for all delay $\tau >0$ if it is stable in the absence of delay.

Assuming contrary that the values of $D_i\left(i=1,2,3,4\right)$ 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 $D_4<0$, which implies

$A_4^2-B_4^2<0.$ (3.7)

Now if condition (3.7) holds, then Equation (3.6) has a positive root ${\eta }_0$ and thus Equation (3.5) has a pair of purely imaginary roots $\pm i{\omega }_0$. It follows from Equations (3.4) that

$\cos \omega \tau =\frac{\left({\omega }_0^4-A_2{\omega }_0^2+A_4\right)\left(B_2{\omega }_0^2-B_4\right)+\left(A_3{\omega }_0-A_1{\omega }_0^3\right)\left(B_1{\omega }_0^2-B_3{\omega }_0\right)}{{\left(B_2{\omega }_0^2-B_4\right)}^2+{\left(B_1{\omega }_0^3-B_3{\omega }_0\right)}^2}.$

Then ${\tau }_k$ corresponding to this positive value of ${\omega }_0$ is given as follows:

$\begin{array}{c}{\tau }_k=\frac{1}{{\omega }_0}\arccos [\frac{\left({\omega }_0^4-A_2{\omega }_0^2+A_4\right)\left(B_2{\omega }_0^2-B_4\right)}{{\left(B_2{\omega }_0^2-B_4\right)}^2+{\left(B_1{\omega }_0^3-B_3{\omega }_0\right)}^2}\\ +\frac{\left(A_3{\omega }_0-A_1{\omega }_0^3\right)\left(B_1{\omega }_0^2-B_3{\omega }_0\right)}{{\left(B_2{\omega }_0^2-B_4\right)}^2+{\left(B_1{\omega }_0^3-B_3{\omega }_0\right)}^2}]+\frac{2k\text{π}}{{\omega }_0},\left(k=0,1,2,3,\cdots \right).\end{array}$

By using Butler’s Lemma, we can say that the endemic equilibrium $E^{\mathrm{<mo>\; *\; </mo>}}$ remains stable for $\tau <{\tau }_0$.

Next we investigate whether there is a phenomenon of Hopf bifurcation as $\tau$ increases through ${\tau }_0$. For this the following lemma is needed.

Lemma 3.1. The following transversality condition is satisfied:

${\text{sign}\frac{\text{d}\left(\mathrm{Re}\left(\lambda \right)\right)}{\text{d}\tau }|}_{\tau ={\tau }_0}>0,$ (3.8)

provided that condition (3.7) holds.

Proof Differentiating Equation (3.1) with respect to $\tau$, we get

${\left(\frac{\text{d}\lambda }{\text{d}\tau }\right)}^{-1}=\frac{P_1+P_2{\text{e}}^{-\lambda \tau }}{\lambda P_3{\text{e}}^{-\lambda \tau }}-\frac{\tau }{\lambda }=\frac{P_1}{\lambda P_3{\text{e}}^{-\lambda \tau }}+\frac{P_2}{\lambda P_3}-\frac{\tau }{\lambda },$ (3.9)

where

$P_1=4{\lambda }^3+3A_1{\lambda }^2+2A_2\lambda +A_3,$

$P_2=3B_1{\lambda }^2+2B_2\lambda +B_3,$

$P_3=B_1{\lambda }^3+B_2{\lambda }^2+B_3\lambda +B_4.$

Therefore,

$\begin{array}{l}\text{sign}{\left\{\frac{\text{d}\left(\mathrm{Re}\lambda \right)}{\text{d}\tau }\right\}}_{\tau ={\tau }_k}=\text{sign}{\left\{\mathrm{Re}{\left(\frac{\text{d}\lambda }{\text{d}\tau }\right)}^{-1}\right\}}_{\lambda =i{\omega }_0}\\ =\text{sign}\left\{\mathrm{Re}{\left[\frac{P_1}{\lambda P_3{\text{e}}^{-\lambda \tau }}\right]}_{\lambda =i{\omega }_0}+\mathrm{Re}{\left[\frac{P_2}{\lambda P_3}\right]}_{\lambda =i{\omega }_0}\right\}\\ =\text{sign}\{\frac{\left(A_3-3A_1{\omega }_0^2\right)\left(-A_1{\omega }_0^4+A_3{\omega }_0^2\right)}{{\left(B_1{\omega }_0^4-B_2{\omega }_0^2\right)}^2+{\left(B_4{\omega }_0-B_2{\omega }_0^3\right)}^2}\\ -\frac{\left(2A_2{\omega }_0-4{\omega }_0^3\right)\left({\omega }_0^5-A_2{\omega }_0^3+A_4{\omega }_0\right)}{{\left(B_1{\omega }_0^4-B_2{\omega }_0^2\right)}^2+{\left(B_4{\omega }_0-B_2{\omega }_0^3\right)}^2}\\ +\frac{2B_2{\omega }_0\left(B_4{\omega }_0-B_2{\omega }_0^3\right)+\left(B_3-3B_1{\omega }_0^2\right)\left(B_1{\omega }_0^4-B_3{\omega }_0^2\right)}{{\left(B_1{\omega }_0^4-B_2{\omega }_0^2\right)}^2+{\left(B_4{\omega }_0-B_2{\omega }_0^3\right)}^2}\}\end{array}$

$\begin{array}{l}=\text{sign}\{\frac{4{\omega }_0^6+\left(3A_1^2-6A_2-3B_1^2\right){\omega }_0^4}{{\left(B_1{\omega }_0^3-B_2{\omega }_0\right)}^2+{\left(B_4-B_2{\omega }_0^2\right)}^2}\\ +\frac{\left(-4A_1{A}_3+2A_2^2+4A_4-2B_2^2+4B_1{B}_3\right){\omega }_0^2}{{\left(B_1{\omega }_0^3-B_2{\omega }_0\right)}^2+{\left(B_4-B_2{\omega }_0^2\right)}^2}\\ +\frac{A_3^2-2A_2{A}_4+2B_2{B}_4-B_3^2}{{\left(B_1{\omega }_0^3-B_2{\omega }_0\right)}^2+{\left(B_4-B_2{\omega }_0^2\right)}^2}\}\\ =\text{sign}\left\{\frac{4{\omega }_0^6+3D_1{\omega }_0^4+2D_2{\omega }_0^2+D_3}{{\left(B_1{\omega }_0^3-B_2{\omega }_0\right)}^2+{\left(B_4-B_2{\omega }_0^2\right)}^2}\right\}.\end{array}$

From the above argument, we know ${\eta }_0={\omega }_0^2$, then

$\begin{array}{l}\text{sign}{\left\{\frac{\text{d}\left(\mathrm{Re}\lambda \right)}{\text{d}\tau }\right\}}_{\tau ={\tau }_k}=\text{sign}{\left\{\mathrm{Re}{\left(\frac{\text{d}\lambda }{\text{d}\tau }\right)}^{-1}\right\}}_{\lambda =i{\omega }_0}\\ =\frac{4{\eta }_0^3+3D_1{\eta }_0^2+2D_2{\eta }_0+D_3}{{\left(B_1{\omega }_0^3-B_2{\omega }_0\right)}^2+{\left(B_4-B_2{\omega }_0^2\right)}^2}\\ =\frac{f^{\prime }\left({\eta }_0\right)}{{\left(B_1{\omega }_0^3-B_2{\omega }_0\right)}^2+{\left(B_4-B_2{\omega }_0^2\right)}^2}.\end{array}$ (3.10)

Here it may be noted that $f^{\prime }\left({\eta }_0\right)>0$ if the condition (3.7) is satisfied. This proves the Lemma 3.1. Thus we have the following result:

Theorem 3.3. If $R_0>1$ and the condition (3.7) hold, then the endemic equilibrium $E^{\mathrm{<mo>\; *\; </mo>}}$ of system (1.3) remains stable for all $\tau \in \left[\mathrm{0,}{\tau }_0\right)$ and becomes unstable for $\tau >{\tau }_0$. System (1.3) with $\tau ={\tau }_0$ 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 $\tau >{\tau }_0$ near ${\tau }_0$. 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 $\tau =\mu +{\tau }_k$. So, $\mu =0$ is the Hopf value of system (1.3).

Let $u_1=X-X^{*}$, $u_2=Y-Y^{*}$, $u_3=S-S^{*}$, $u_4=I-I^{*}$, ${\bar{u}}_i\left(t\right)=u_i\left(\tau t\right)$, and dropping the bars for simplification of notations, system (1.3) becomes a functional differential equation in $C=C\left(\left[-1,0\right],{ℝ}^4\right)$ as

$\stackrel{\dot{}}{u}\left(t\right)=L_{\mu }\left(u_t\right)+f\left(\mu ,u_t\right),$ (4.1)

where $u\left(t\right)={\left(u_1\left(t\right)\mathrm{,}{u}_2\left(t\right)\mathrm{,}{u}_3\left(t\right)\mathrm{,}{u}_4\left(t\right)\right)}^{\text{T}}\in {ℝ}^4$, and $L_{\mu }\mathrm{<mo>\; :\; </mo>}C\to {ℝ}^4$, $f\mathrm{<mo>\; :\; </mo>}ℝ\times C\to {ℝ}^4$ are given, respectively, by

$\begin{array}{c}{L}_{\mu }\left(\varphi \right)=\left({\tau }_k+\mu \right)\left(\begin{array}{cccc}-d_1& v& 0& -{\beta }_1{X}^{\mathrm{<mo>\; *\; </mo>}}\\ 0& -\left(v+{\alpha }_1+d_1\right)& 0& {\beta }_1{X}^{\mathrm{<mo>\; *\; </mo>}}\\ 0& -{\beta }_2{S}^{\mathrm{<mo>\; *\; </mo>}}& -{\beta }_2{Y}^{\mathrm{<mo>\; *\; </mo>}}-d_2& 0\\ 0& {\beta }_2{S}^{\mathrm{<mo>\; *\; </mo>}}& -{\beta }_2{Y}^{\mathrm{<mo>\; *\; </mo>}}& -\left({\alpha }_2+d_2\right)\end{array}\right)\left(\begin{array}{c}{\varphi }_1\left(0\right)\\ {\varphi }_2\left(0\right)\\ {\varphi }_3\left(0\right)\\ {\varphi }_4\left(0\right)\end{array}\right)\\ +\left({\tau }_k+\mu \right)\left(\begin{array}{cccc}-{\beta }_1{I}^{\mathrm{<mo>\; *\; </mo>}}& 0& 0& 0\\ {\beta }_1{I}^{\mathrm{<mo>\; *\; </mo>}}& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right)\left(\begin{array}{c}{\varphi }_1\left(-1\right)\\ {\varphi }_2\left(-1\right)\\ {\varphi }_3\left(-1\right)\\ {\varphi }_4\left(-1\right)\end{array}\right)\mathrm{,}\end{array}$ (4.2)

and

$f\left(\mu \mathrm{,}\varphi \right)=\left({\tau }_k+\mu \right)\left(\begin{array}{c}-{\beta }_1{\varphi }_4\left(0\right){\varphi }_1\left(-1\right)\\ {\beta }_1{\varphi }_4\left(0\right){\varphi }_1\left(-1\right)\\ 0\\ 0\end{array}\right)\mathrm{.}$ (4.3)

By the Riesz representation theorem, there exists a function $\eta \left(\theta \mathrm{,}\mu \right)$ of bounded variation for $\theta \in \left[-\mathrm{1,0}\right]$ such that

$L_{\mu }\left(\varphi \right)={\displaystyle {\int }_{-1}^0}\text{d}\eta \left(\theta ,\mu \right)\varphi \left(\theta \right)\text{for}\varphi \in C.$ (4.4)

In fact, we can choose

$\begin{array}{c}\eta \left(\theta ,\mu \right)=\left({\tau }_k+\mu \right)\left(\begin{array}{cccc}-d_1& v& 0& -{\beta }_1{X}^{*}\\ 0& -\left(v+{\alpha }_1+d_1\right)& 0& {\beta }_1{X}^{*}\\ 0& -{\beta }_2{S}^{*}& -{\beta }_2{Y}^{*}-d_2& 0\\ 0& {\beta }_2{S}^{*}& -{\beta }_2{Y}^{*}& -\left({\alpha }_2+d_2\right)\end{array}\right)\delta \left(\theta \right)\\ +\left({\tau }_k+\mu \right)\left(\begin{array}{cccc}-{\beta }_1{I}^{*}& 0& 0& 0\\ {\beta }_1{I}^{*}& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right)\delta \left(\theta +1\right),\end{array}$ (4.5)

where $\delta$ denote the Dirac delta function:

$\delta \left(\theta \right)=\{\begin{array}{l}0,\theta \ne 0,\\ 1,\theta =0.\end{array}$

For $\varphi \in C\left(\left[-\mathrm{1,0}\right]\mathrm{,}{ℝ}^4\right)$, define

$A\left(\mu \right)\left(\varphi \right)=\{\begin{array}{l}\frac{\text{d}\varphi \left(\theta \right)}{\text{d}\theta },\theta \in \left[-1,0\right),\\ {\displaystyle {\int }_{-1}^0}\text{d}\eta \left(s,\mu \right)\varphi \left(s\right),\theta =0,\end{array}$

and

$R\left(\mu \right)\left(\varphi \right)=\{\begin{array}{l}0,\theta \in \left[-1,0\right),\\ f\left(\mu ,\varphi \right),\theta =0.\end{array}$

Then system (1.3) is equivalent to

$\stackrel{\dot{}}{u}\left(t\right)=A\left(\mu \right)u_t+R\left(\mu \right)u_t\mathrm{,}$ (4.6)

where $u_t\left(\theta \right)=u\left(t+\theta \right)$ for $\theta \in \left[-\mathrm{1,0}\right]$.

For $\psi \in C^1\left(\left[\mathrm{0,1}\right]\mathrm{,}{\left({ℝ}^4\right)}^{\mathrm{<mo>\; *\; </mo>}}\right)$, define

$A^{\mathrm{<mo>\; *\; </mo>}}\psi \left(s\right)=\{\begin{array}{l}-\frac{\text{d}\psi \left(s\right)}{\text{d}s}\mathrm{,}s\in \left[-\mathrm{1,0}\right)\mathrm{,}\\ {\displaystyle {\int }_{-1}^0}\psi \left(-t\right)\text{d}\eta \left(t\mathrm{,0}\right)\mathrm{,}s\mathrm{<mo>\; =\; </mo>0,}\end{array}$

and a bilinear inner product

$\langle \psi \left(s\right),\varphi \left(\theta \right)\rangle =\bar{\psi }\left(0\right)\varphi \left(0\right)-{\displaystyle {\int }_{-1}^0}{\displaystyle {\int }_{\xi =0}^{\theta }}\bar{\psi }\left(\xi -\theta \right)\text{d}\eta \left(\theta \right)\varphi \left(\xi \right)\text{d}\xi ,$ (4.7)

where $\eta \left(\theta \right)=\eta \left(\theta \mathrm{,0}\right)$. Then $A\left(0\right)$ and $A^{\mathrm{<mo>\; *\; </mo>}}$ are adjoint operators. By the above discussion, we know that $\pm i\omega {\tau }_k$ are eigenvalues of $A\left(0\right)$. Hence, they are also eigenvalues of $A^{\mathrm{<mo>\; *\; </mo>}}$. We first need to compute the eigenvectors of $A\left(0\right)$ and $A^{\mathrm{<mo>\; *\; </mo>}}$ corresponding to $i\omega {\tau }_k$ and $-i\omega {\tau }_k$, respectively.

Supposed $q\left(\theta \right)={\left(\mathrm{1,}{q}_1\mathrm{,}{q}_2\mathrm{,}{q}_3\right)}^{\text{T}}{\text{e}}^{i\omega {\tau }_k\theta }$ is the eigenvectors of $A\left(0\right)$ corresponding to $i\omega {\tau }_k$, then $A\left(0\right)q\left(\theta \right)=i\omega {\tau }_{k}q\left(\theta \right)$. Then from the definition of $A\left(0\right)$ and (4.2), (4.4) and (4.5), we have

$\begin{array}{l}\left(\begin{array}{cccc}-d_1& v& 0& -{\beta }_1{X}^{*}\\ 0& -\left(v+{\alpha }_1+d_1\right)& 0& {\beta }_1{X}^{*}\\ 0& -{\beta }_2{S}^{*}& -{\beta }_2{Y}^{*}-d_2& 0\\ 0& {\beta }_2{S}^{*}& -{\beta }_2{Y}^{*}& -\left({\alpha }_2+d_2\right)\end{array}\right)q\left(0\right)\\ +\left(\begin{array}{cccc}-{\beta }_1{I}^{*}& 0& 0& 0\\ {\beta }_1{I}^{*}& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right)q\left(-1\right)=i{\omega }_{0}q\left(0\right).\end{array}$

For $q\left(-1\right)={\left(1,q_1,q_2,q_3\right)}^{\text{T}}{\text{e}}^{-i{\omega }_0{\tau }_k}=q\left(0\right){\text{e}}^{-i{\omega }_0{\tau }_k}$, then we obtain

$\{\begin{array}{l}{q}_1=\frac{-i{\omega }_0-d_1}{{\alpha }_1+d_1+i{\omega }_0},\hfill \\ q_2=\frac{{\beta }_2{S}^{*}\left(i{\omega }_0+d_1\right)}{\left({\alpha }_1+d_1+i{\omega }_0\right)\left(d_2+{\beta }_2{Y}^{*}+i{\omega }_0\right)},\hfill \\ q_3=\frac{-{\beta }_2{S}^{*}\left(i{\omega }_0+d_1\right)\left(i{\omega }_0+d_2\right)}{\left({\alpha }_1+d_1+i{\omega }_0\right)\left({\alpha }_2+d_2+i{\omega }_0\right)\left(d_2+{\beta }_2{Y}^{*}+i{\omega }_0\right)}.\hfill \end{array}$

Similarly, we can obtain the eigenvector $q^{\mathrm{<mo>\; *\; </mo>}}\left(s\right)=D\left(\mathrm{1,}{q}_1^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{q}_2^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{q}_3^{\mathrm{<mo>\; *\; </mo>}}\right){\text{e}}^{i{\omega }_0{\tau }_{k}s}$ of $A^{\mathrm{<mo>\; *\; </mo>}}$ corresponding to $-i{\omega }_0{\tau }_k$, where

$\{\begin{array}{l}{q}_1^{*}=\frac{d_1+{\beta }_1{I}^{*}{\text{e}}^{-i{\omega }_0{\tau }_k}-i{\omega }_0}{{\beta }_1{I}^{*}{\text{e}}^{-i{\omega }_0{\tau }_k}},\hfill \\ q_2^{*}=\frac{{\beta }_1{\beta }_2{X}^{*}{Y}^{*}\left(d_1-i{\omega }_0\right)}{{\beta }_1{I}^{*}{\text{e}}^{-i{\omega }_0{\tau }_k}\left({\alpha }_2+d_2-i{\omega }_0\right)\left(d_2+{\beta }_2{Y}^{*}-i{\omega }_0\right)},\hfill \\ q_3^{*}=\frac{{\beta }_1{X}^{*}\left(d_1-i{\omega }_0\right)}{{\beta }_1{I}^{*}{\text{e}}^{-i{\omega }_0{\tau }_k}\left({\alpha }_2+d_2-i{\omega }_0\right)}.\hfill \end{array}$

In order to assure $\langle q^{*}\left(s\right),q\left(\theta \right)\rangle =1$, we need to determine the value of D. By (4.7), we have

$\begin{array}{l}\langle q^{*}\left(s\right),q\left(\theta \right)\rangle =\bar{D}\left(1,{\bar{q}}_1^{*},{\bar{q}}_2^{*},{\bar{q}}_3^{*}\right){\left(1,q_1,q_2,q_3\right)}^{\text{T}}\\ -{\displaystyle {\int }_{-1}^0}{\displaystyle {\int }_{\xi =0}^{\theta }}\bar{D}\left(1,{\bar{q}}_1^{*},{\bar{q}}_2^{*},{\bar{q}}_3^{*}\right){\text{e}}^{-i{\omega }_0{\tau }_k\left(\xi -\theta \right)}\text{d}\eta \left(\theta \right){\left(1,q_1,q_2,q_3\right)}^{\text{T}}{\text{e}}^{i{\omega }_0{\tau }_k\xi }\text{d}\xi \\ =\bar{D}\left\{1+q_1{\bar{q}}_1^{*}+q_2{\bar{q}}_2^{*}+q_3{\bar{q}}_3^{*}-{\displaystyle {\int }_{-1}^0}\left(1,{\bar{q}}_1^{*},{\bar{q}}_2^{*},{\bar{q}}_3^{*}\right)\theta {\text{e}}^{i{\omega }_0{\tau }_k\theta }\text{d}\eta \left(\theta \right){\left(1,q_1,q_2,q_3\right)}^{\text{T}}\right\}\\ =\bar{D}\left\{1+q_1{\bar{q}}_1^{*}+q_2{\bar{q}}_2^{*}+q_3{\bar{q}}_3^{*}+{\tau }_k{\beta }_1{I}^{*}{\text{e}}^{-i{\omega }_0{\tau }_k}\left(-1+{\bar{q}}_1^{*}\right)\right\}.\end{array}$

Therefore, we can choose D as

$D=\frac{1}{1+{\bar{q}}_1{q}_1^{*}+{\bar{q}}_2{q}_2^{*}+{\bar{q}}_3{q}_3^{*}+{\tau }_k{\beta }_1{I}^{*}{\text{e}}^{i{\omega }_0{\tau }_k}\left(-1+q_1^{*}\right)}.$

Next we will compute the coordinate to describe the center manifold $C_0$ at $\mu =0$. Let $u_t$ be the solution of (4.6) when $\mu =0$. Define

$z\left(t\right)=\langle q^{*}\left(s\right),u_t\left(\theta \right)\rangle ,W\left(t,\theta \right)=u_t\left(\theta \right)-2\mathrm{Re}\left\{z\left(t\right)q\left(\theta \right)\right\}.$ (4.8)

On the center manifold $C_0$, we have

$W\left(t,\theta \right)=W\left(z\left(t\right),\bar{z}\left(t\right),\theta \right)=W_{20}\left(\theta \right)\frac{z^2}{2}+W_{11}\left(\theta \right)z\bar{z}+W_{02}\left(\theta \right)\frac{{\bar{z}}^2}{2}+\cdots ,$ (4.9)

where z and $\bar{z}$ are local coordinates for center manifold $C_0$ in the direction of $q^{\mathrm{<mo>\; *\; </mo>}}$ and ${\bar{q}}^{\mathrm{<mo>\; *\; </mo>}}$. Note that W is real if $u_t$ is real. We only consider real solutions. For solution $u_t\in C_0$ of (4.6), since $\mu =0$, we obtain

$\begin{array}{c}{z}^{\prime }\left(t\right)=\langle q^{*},{u^{\prime }}_t\rangle =i{\omega }_0{\tau }_{k}z+{\bar{q}}^{*}\left(0\right)f\left(0,W\left(z,\bar{z},0\right)+2\mathrm{Re}\left\{zq\left(0\right)\right\}\right)\\ \triangleq i{\omega }_0{\tau }_{k}z+{\bar{q}}^{*}\left(0\right)f_0\left(z,\bar{z}\right)=i{\omega }_0{\tau }_{k}z+g\left(z,\bar{z}\right),\end{array}$ (4.10)

where

$g\left(z,\bar{z}\right)={\bar{q}}^{*}\left(0\right)f_0\left(z,\bar{z}\right)=g_{20}\left(\theta \right)\frac{z^2}{2}+g_{11}\left(\theta \right)z\bar{z}+g_{02}\left(\theta \right)\frac{{\bar{z}}^2}{2}+g_{21}\frac{z^2\bar{z}}{2}+\cdots .$ (4.11)

It follows from (4.8) and (4.9) that $u_t\left(\theta \right)=\left(u_{1t}\left(\theta \right),u_{2t}\left(\theta \right),u_{3t}\left(\theta \right),u_{4t}\left(\theta \right)\right)=W\left(t,\theta \right)+zq\left(\theta \right)+\overline{zq\left(\theta \right)}$ and $q\left(\theta \right)={\left(1,q_1,q_2,q_3\right)}^{\text{T}}{\text{e}}^{i{\omega }_0{\tau }_0\theta }$, then

$u_{4t}\left(0\right)=W_{20}^{\left(4\right)}\left(0\right)\frac{z^2}{2}+W_{11}^{\left(4\right)}\left(0\right)z\bar{z}+W_{02}^{\left(4\right)}\left(0\right)\frac{{\bar{z}}^2}{2}+q_{3}z+\overline{q_{3}z}+o\left({\left|\left(z,\bar{z}\right)\right|}^3\right),$

$\begin{array}{c}{u}_{1t}\left(-1\right)=W_{20}^{\left(1\right)}\left(-1\right)\frac{z^2}{2}+W_{11}^{\left(1\right)}\left(-1\right)z\bar{z}+W_{02}^{\left(1\right)}\left(-1\right)\frac{{\bar{z}}^2}{2}\\ +z{\text{e}}^{-i{\omega }_0{\tau }_k}+\bar{z}{\text{e}}^{i{\omega }_0{\tau }_k}+o\left({\left|\left(z,\bar{z}\right)\right|}^3\right).\end{array}$

Then from the definition of $f\left(\mu \mathrm{,}{u}_t\right)$, we obtain

$\begin{array}{l}g\left(z,\bar{z}\right)={\bar{q}}^{*}\left(0\right)f_0\left(z,\bar{z}\right)={\bar{q}}^{*}\left(0\right)f\left(0,u_t\right)\\ ={\tau }_k\bar{D}\left(1,{\bar{q}}_1^{*},{\bar{q}}_2^{*},{\bar{q}}_3^{*}\right)\left(\begin{array}{c}-{\beta }_1{u}_{4t}\left(0\right)u_{1t}\left(-1\right)\\ {\beta }_1{u}_{4t}\left(0\right)u_{1t}\left(-1\right)\\ 0\\ 0\end{array}\right)\\ ={\tau }_k\bar{D}{\beta }_1\left(-1+{\bar{q}}_1^{*}\right)\{W_{20}^{\left(4\right)}\left(0\right)\frac{z^2}{2}+W_{11}^{\left(4\right)}\left(0\right)z\bar{z}+W_{02}^{\left(4\right)}\left(0\right)\frac{{\bar{z}}^2}{2}+q_{3}z\\ \underset{}{\overset{}{}}+\overline{q_{3}z}+o\left({\left|\left(z,\bar{z}\right)\right|}^3\right)\}\{W_{20}^{\left(1\right)}\left(-1\right)\frac{z^2}{2}+W_{11}^{\left(1\right)}\left(-1\right)z\bar{z}\\ +W_{02}^{\left(1\right)}\left(-1\right)\frac{{\bar{z}}^2}{2}+z{\text{e}}^{-i{\omega }_0{\tau }_0}+\bar{z}{\text{e}}^{i{\omega }_0{\tau }_0}+o\left({\left|\left(z,\bar{z}\right)\right|}^3\right)\}\end{array}$ (4.12)

Comparing the coefficients with (4.11), we have

$g_{20}=2{\tau }_0\bar{D}{\beta }_1\left(-1+{\bar{q}}_1^{*}\right)q_3{\text{e}}^{-i{\omega }_0{\tau }_k},$

$g_{11}=2{\tau }_0\bar{D}{\beta }_1\left(-1+{\bar{q}}_1^{*}\right)\mathrm{Re}\left\{q_3{\text{e}}^{i{\omega }_0{\tau }_k}\right\},$

$g_{02}=2{\tau }_0\bar{D}{\beta }_1\left(-1+{\bar{q}}_1^{*}\right){\bar{q}}_3{\text{e}}^{i{\omega }_0{\tau }_k},$

$g_{21}=2{\tau }_0\bar{D}{\beta }_1\left(-1+{\bar{q}}_1^{*}\right)\left[W_{20}^{\left(4\right)}{\text{e}}^{i{\omega }_0{\tau }_k}+W_{20}^{\left(1\right)}\left(-1\right){\bar{q}}_3+2q_3{W}_{11}^{\left(1\right)}\left(-1\right)+2W_{11}^{\left(4\right)}\left(0\right){\text{e}}^{-i{\omega }_0{\tau }_k}\right].$

In order to assure $g_{21}$, $W_{20}\left(\theta \right)$ and $W_{11}\left(\theta \right)$ are needed to compute.

From (4.6), (4.8) and (4.10), we have

$\begin{array}{c}\stackrel{\dot{}}{W}={\stackrel{\dot{}}{u}}_t-z^{\prime }q-\stackrel{\dot{}}{\bar{z}}\bar{q}\\ =\{\begin{array}{l}A\left(0\right)W-2\mathrm{Re}\left\{{\bar{q}}^{*}\left(0\right)f_0\left(z,\bar{z}\right)q\left(\theta \right)\right\},\theta \in \left[-1,0\right),\\ A\left(0\right)W-2\mathrm{Re}\left\{{\bar{q}}^{*}\left(0\right)f_0\left(z,\bar{z}\right)q\left(\theta \right)\right\}+f_0\left(z,\bar{z}\right),\theta =0,\end{array}\\ =A\left(0\right)W+H\left(z,\bar{z},\theta \right),\end{array}$ (4.13)

where

$H\left(z\mathrm{,}\bar{z}\mathrm{,}\theta \right)=H_{20}\left(\theta \right)\frac{z^2}{2}+H_{11}\left(\theta \right)z\bar{z}+H_{02}\left(\theta \right)\frac{{\bar{z}}^2}{2}+\cdots \mathrm{.}$ (4.14)

It follows from (4.13) and (4.14) that

$A\left(0\right)W-\stackrel{\dot{}}{W}=-H\left(z,\bar{z},\theta \right)=-H_{20}\left(\theta \right)\frac{z^2}{2}-H_{11}\left(\theta \right)z\bar{z}-H_{02}\left(\theta \right)\frac{{\bar{z}}^2}{2}-\cdots .$ (4.15)

From (4.9) and (4.10), we have

$A\left(0\right)W=A\left(0\right)W_{20}\left(\theta \right)\frac{z^2}{2}+A\left(0\right)W_{11}\left(\theta \right)z\bar{z}+A\left(0\right)W_{02}\left(\theta \right)\frac{{\bar{z}}^2}{2}+\cdots$ (4.16)

and

$\begin{array}{l}\stackrel{\dot{}}{W}=W_z\stackrel{\dot{}}{z}+W_{\bar{z}}\stackrel{\dot{}}{\bar{z}}\\ =\left(W_{20}\left(\theta \right)z+W_{11}\left(\theta \right)\bar{z}\right)\left(i{\omega }_0{\tau }_{k}z+g_{20}\left(\theta \right)\frac{z^2}{2}+g_{11}\left(\theta \right)z\bar{z}+g_{02}\left(\theta \right)\frac{{\bar{z}}^2}{2}+\cdots \right)\\ +\left(W_{11}\left(\theta \right)z+W_{02}\left(\theta \right)\bar{z}\right)\left(-i{\omega }_0{\tau }_k\bar{z}+g_{20}\left(\theta \right)\frac{{\bar{z}}^2}{2}+g_{11}\left(\theta \right)z\bar{z}+g_{02}\left(\theta \right)\frac{z^2}{2}+\cdots \right)+\cdots .\end{array}$ (4.17)

Substituting (4.16) and (4.17) into (4.15) and comparing the coefficients of $z^2$ and $z\bar{z}$, we have

$\left(A\left(0\right)-2i{\omega }_0{\tau }_{k}I\right)W_{20}\left(\theta \right)=-H_{20}\left(\theta \right)\mathrm{,}A\left(0\right)W_{11}\left(\theta \right)=-H_{11}\left(\theta \right)\mathrm{.}$ (4.18)

From (4.11) and (4.13), we know that for $\theta \in \left[-\mathrm{1,0}\right)$

$\begin{array}{c}H\left(z,\bar{z},\theta \right)=-2\mathrm{Re}{\bar{q}}^{*}\left(0\right)f_{0}q\left(\theta \right)=-g\left(z,\bar{z}\right)q\left(\theta \right)-\bar{g}\left(z,\bar{z}\right)\bar{q}\left(\theta \right)\\ =-\frac{1}{2}\left(g_{20}q\left(\theta \right)+{\bar{g}}_{20}\bar{q}\left(\theta \right)\right)z^2-\left(g_{11}q\left(\theta \right)+{\bar{g}}_{11}\bar{q}\left(\theta \right)\right)z\bar{z}-\cdots .\end{array}$ (4.19)

Comparing the coefficients with (4.19) gives that

$H_{20}\left(\theta \right)=-g_{20}q\left(\theta \right)-{\bar{g}}_{02}\bar{q}\left(\theta \right),$ (4.20)

and

$H_{11}\left(\theta \right)=-g_{11}q\left(\theta \right)-{\bar{g}}_{11}\bar{q}\left(\theta \right).$ (4.21)

From the definition of $A\left(0\right)$ and (4.18) and (4.21), we obtain

${\stackrel{\dot{}}{W}}_{20}\left(\theta \right)=2i{\omega }_0{\tau }_k{W}_{\left(20\right)}\left(\theta \right)+g_{20}q\left(\theta \right)+{\bar{g}}_{20}\bar{q}\left(\theta \right).$

For $q\left(\theta \right)=q\left(0\right){\text{e}}^{i{\omega }_0{\tau }_k\theta }$, we have

$W_{20}\left(\theta \right)=\frac{i{g}_{20}}{{\omega }_0{\tau }_k}q\left(0\right){\text{e}}^{i{\omega }_0{\tau }_k\theta }+\frac{i{\bar{g}}_{02}}{3{\omega }_0{\tau }_k}\bar{q}\left(0\right){\text{e}}^{-i{\omega }_0{\tau }_k\theta }+E_1{\text{e}}^{2i{\omega }_0{\tau }_k\theta },$ (4.22)

where $E_1={\left(E_1^{\left(1\right)}\mathrm{,}{E}_1^{\left(2\right)}\mathrm{,}{E}_1^{\left(3\right)}\mathrm{,}{E}_1^{\left(4\right)}\right)}^{\text{T}}$ is a constant vector. Similarly, from (4.18) and (4.22), we know

$W_{11}\left(\theta \right)=-\frac{i{g}_{11}}{{\omega }_0{\tau }_k}q\left(0\right){\text{e}}^{i{\omega }_0{\tau }_k\theta }+\frac{i{\bar{g}}_{11}}{{\omega }_0{\tau }_k}\bar{q}\left(0\right){\text{e}}^{-i{\omega }_0{\tau }_k\theta }+E_2,$ (4.23)

where $E_2={\left(E_2^{\left(1\right)}\mathrm{,}{E}_2^{\left(2\right)}\mathrm{,}{E}_2^{\left(3\right)}\mathrm{,}{E}_2^{\left(4\right)}\right)}^{\text{T}}$ is a constant vector.

Finally, we will seek the values of $E_1$ and $E_2$. From the definition of $A\left(0\right)$ and (4.18), we have

${\displaystyle {\int }_{-1}^0}\text{d}\eta \left(\theta \right)W_{20}\left(\theta \right)=2i{\omega }_0{\tau }_k{W}_{20}\left(0\right)-H_{20}\left(0\right)$ (4.24)

and

${\displaystyle {\int }_{-1}^0}\text{d}\eta \left(\theta \right)W_{11}\left(\theta \right)=-H_{11}\left(0\right),$ (4.25)

where $\eta \left(\theta \right)=\eta \left(0,\theta \right)$. By (4.13), we know when $\theta =0$

$H\left(z,\bar{z},0\right)=f_0-g\left(z,\bar{z}\right)q\left(0\right)-\bar{g}\left(z,\bar{z}\right)\bar{q}\left(0\right).$

That is

$\begin{array}{l}{H}_{20}\left(\theta \right)\frac{z^2}{2}+H_{11}\left(\theta \right)z\bar{z}+H_{02}\left(\theta \right)\frac{{\bar{z}}^2}{2}+\cdots \\ =-q\left(0\right)\left(g_{20}\frac{z^2}{2}+g_{11}z\bar{z}+g_{02}\frac{{\bar{z}}^2}{2}+\cdots \right)\\ -\bar{q}\left(0\right)\left({\bar{g}}_{20}\frac{{\bar{z}}^2}{2}+{\bar{g}}_{11}z\bar{z}+{\bar{g}}_{02}\frac{z^2}{2}+\cdots \right)+f_0.\end{array}$ (4.26)

By (4.3), we have

$f_0={\tau }_k\left(\begin{array}{c}-{\beta }_1{u}_{4t}\left(0\right)u_{1t}\left(-1\right)\\ {\beta }_1{u}_{4t}\left(0\right)u_{1t}\left(-1\right)\\ 0\\ 0\end{array}\right)$

By (4.8), we obtain

$\begin{array}{c}{u}_t\left(\theta \right)=W\left(t\mathrm{,}\theta \right)+2\mathrm{Re}\left\{z\left(t\right)q\left(\theta \right)\right\}=W\left(t\mathrm{,}\theta \right)+z\left(t\right)q\left(\theta \right)+\bar{z}\left(t\right)\bar{q}\left(\theta \right)\\ =W_{20}\left(\theta \right)\frac{z^2}{2}+W_{11}\left(\theta \right)z\bar{z}+z\left(t\right)q\left(\theta \right)+\bar{z}\left(t\right)\bar{q}\left(\theta \right)+\cdots \mathrm{.}\end{array}$

Then, we have

$f_0={\tau }_k\left(\begin{array}{c}-{\beta }_1{q}_3{\text{e}}^{-i{\omega }_0{\tau }_k}\\ {\beta }_1{q}_3{\text{e}}^{-i{\omega }_0{\tau }_k}\\ 0\\ 0\end{array}\right)z^2+2{\tau }_k\mathrm{Re}\left\{q_3{\text{e}}^{i{\omega }_0{\tau }_k}\right\}\left(\begin{array}{c}-{\beta }_1\\ {\beta }_1\\ 0\\ 0\end{array}\right)z\bar{z}+\cdots .$ (4.27)

By (4.26) and (4.27), we have

$H_{20}\left(0\right)=-g_{20}q\left(0\right)-{\bar{g}}_{02}\bar{q}\left(0\right)+2{\tau }_k{q}_3{\text{e}}^{-i{\omega }_0{\tau }_k}\left(\begin{array}{c}-{\beta }_1\\ {\beta }_1\\ 0\\ 0\end{array}\right)$ (4.28)

and

$H_{11}\left(0\right)=-g_{11}q(0)-{\bar{g}}_{11}\bar{q}\left(0\right)+2{\tau }_k\mathrm{Re}\left\{q_3{\text{e}}^{i{\omega }_0{\tau }_k}\right\}\left(\begin{array}{c}-{\beta }_1\\ {\beta }_1\\ 0\\ 0\end{array}\right).$ (4.29)

For $i{\omega }_0{\tau }_0$ is the eigenvalue of $A\left(0\right)$ and $q\left(0\right)$ is the corresponding eigenvector, we obtain

$\left(i{\omega }_0{\tau }_{k}I-{\displaystyle {\int }_{-1}^0}{\text{e}}^{i{\omega }_0{\tau }_k\theta }\text{d}\eta \left(\theta \right)\right)q\left(0\right)=0$

and

$\left(-i{\omega }_0{\tau }_{k}I-{\displaystyle {\int }_{-1}^0}{\text{e}}^{-i{\omega }_0{\tau }_k\theta }\text{d}\eta \left(\theta \right)\right)\bar{q}\left(0\right)=0.$

So, substituting (4.22) and (4.28) into (4.24), we obtain

$\left(2i{\omega }_0{\tau }_{k}I-{\displaystyle {\int }_{-1}^0}{\text{e}}^{2i{\omega }_0{\tau }_k\theta }\text{d}\eta \left(\theta \right)\right)E_1=2{\tau }_k{q}_3{\text{e}}^{-i{\omega }_0{\tau }_k}\left(\begin{array}{c}-{\beta }_1\\ {\beta }_1\\ 0\\ 0\end{array}\right).$

That is

$\begin{array}{l}\left(\begin{array}{c}\begin{array}{cccc}2i{\omega }_0+d_1+{\beta }_1{I}^{*}{\text{e}}^{-2i{\omega }_0{\tau }_k}& -v& 0& {\beta }_1{X}^{*}\\ -{\beta }_1{I}^{*}{\text{e}}^{-2i{\omega }_0{\tau }_k}& 2i{\omega }_0+\left(v+{\alpha }_1+d_1\right)& 0& -{\beta }_1{X}^{*}\\ 0& {\beta }_2{S}^{*}& 2i{\omega }_0+{\beta }_2{Y}^{*}+d_2& 0\\ 0& -{\beta }_2{S}^{*}& -{\beta }_2{Y}^{*}& 2i{\omega }_0+{\alpha }_2+d_2\end{array}\end{array}\right)E_1\\ =2q_3{\text{e}}^{-i{\omega }_0{\tau }_k}\left(\begin{array}{c}-{\beta }_1\\ {\beta }_1\\ 0\\ 0\end{array}\right).\end{array}$

It follows that

$E_1=2q_3{\text{e}}^{-i{\omega }_0{\tau }_k}{\left(\begin{array}{c}\begin{array}{cccc}2i{\omega }_0+d_1+{\beta }_1{I}^{*}{\text{e}}^{-2i{\omega }_0{\tau }_k}& -v& 0& {\beta }_1{X}^{*}\\ -{\beta }_1{I}^{*}{\text{e}}^{-2i{\omega }_0{\tau }_k}& 2i{\omega }_0+\left(v+{\alpha }_1+d_1\right)& 0& -{\beta }_1{X}^{*}\\ 0& {\beta }_2{S}^{*}& 2i{\omega }_0+{\beta }_2{Y}^{*}+d_2& 0\\ 0& -{\beta }_2{S}^{*}& -{\beta }_2{Y}^{*}& 2i{\omega }_0+{\alpha }_2+d_2\end{array}\end{array}\right)}^{-1}\left(\begin{array}{c}-{\beta }_1\\ {\beta }_1\\ 0\\ 0\end{array}\right).$

Similarly, substituting (4.23) and (4.29) into (4.25), we also get

$-{\displaystyle {\int }_{-1}^0}\text{d}\eta \left(\theta \right)E_2=2{\tau }_k\mathrm{Re}\left\{q_3{\text{e}}^{i{\omega }_0{\tau }_k}\right\}\left(\begin{array}{c}-{\beta }_1\\ {\beta }_1\\ 0\\ 0\end{array}\right),$

and

$E_2=2\mathrm{Re}\left\{q_3{\text{e}}^{i{\omega }_0{\tau }_k}\right\}{\left(\begin{array}{cccc}{d}_1& -v& 0& {\beta }_1{X}^{*}\\ 0& v+{\alpha }_1+d_1& 0& -{\beta }_1{X}^{*}\\ 0& {\beta }_2{S}^{*}& {\beta }_2{Y}^{*}+d_2& 0\\ 0& -{\beta }_2{S}^{*}& {\beta }_2{Y}^{*}& {\alpha }_2+d_2\end{array}\right)}^{-1}\left(\begin{array}{c}-{\beta }_1\\ {\beta }_1\\ 0\\ 0\end{array}\right).$

Thus, we can determine $W_{20}\left(\theta \right)$ and $W_{11}\left(\theta \right)$ from (4.22) and (4.23). Further, we can compute $g_{21}$. Thus we can compute the following values:

$c_1\left(0\right)=\frac{i}{2{\omega }_0{\tau }_k}\left(g_{20}{g}_{11}-2{\left|g_{11}\right|}^2-\frac{{\left|g_{02}\right|}^2}{3}\right)+\frac{g_{21}}{2},$

${\mu }_2=-\frac{\mathrm{Re}{c}_1\left(0\right)}{\mathrm{Re}\left\{\frac{\text{d}\lambda \left({\tau }_k\right)}{\text{d}\tau }\right\}},$

${\beta }_2=2\mathrm{Re}\left(c_1\left(0\right)\right),$

$T_2=-\frac{\mathrm{Im}\left\{c_1\left(0\right)\right\}+{\mu }_2\mathrm{Im}\left\{\frac{\text{d}\lambda \left({\tau }_k\right)}{\text{d}\tau }\right\}}{{\omega }_0{\tau }_k},k=0,1,2,\cdots .$

which determine the qualities of bifurcation periodic solution in the center manifold at the critical values ${\tau }_k$.

According to [25] , we can obtain the following result.

Theorem 4.1. Assume that $R_0>1$ and the condition (3.7) hold, we have:

1) if ${\mu }_2>0$ ( ${\mu }_2<0$ ), then the Hopf bifurcation is supercritical (subcritical) and the bifurcating periodic solutions exist for ${\tau }_0$.

2) if ${\beta }_2<0$ ( ${\beta }_2>0$ ), then the bifurcating periodic solutions are stable (unstable).

3) if $T_2>0$ (<0), then the period of the bifurcating periodic solutions increases (decreases).

From Theorem 4.1, we know that the value of ${\mu }_2$ determines the directions of the Hopf bifurcation, the values of ${\beta }_2$ and $T_2$ 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 ${\alpha }_1=0.7$, ${\alpha }_2=0.30476509$, ${\beta }_1=0.040395636$, ${\beta }_2=0.0082258094$, $A=37.557699$, $v=0.5$, $\Lambda =49.721841$, $d_1=0.5$, $d_2=0.5$. For $k=0$, using the software Matlab, we derive ${\tau }_0\approx 5.476$. Thus the endemic equilibrium $E^{\mathrm{<mo>\; *\; </mo>}}$ is stable when $\tau <{\tau }_0\approx 5.476$. Figures 1(a)-(f) show that the endemic equilibrium $E^{\mathrm{<mo>\; *\; </mo>}}$ is asymptotically stable when $\tau =5.2<{\tau }_0\approx 5.476$. When $\tau$ passes through the critical value ${\tau }_0\approx 5.476$, the endemic equilibrium $E^{\mathrm{<mo>\; *\; </mo>}}$ loses its stability and a Hopf bifurcation occurs, that is, a family of periodic solutions bifurcate from the endemic equilibrium $E^{\mathrm{<mo>\; *\; </mo>}}$. Figures 2(a)-(f) suggest that Hopf bifurcation occurs from the endemic equilibrium $E^{\mathrm{<mo>\; *\; </mo>}}$ when $\tau =6>{\tau }_0\approx 5.476$.

6. Conclusion

In this paper, we have investigated a delayed schistomiasis model, and studied

(a) (b)(c) (d)(e) (f)

Figure 1. (a)-(f) The dynamical behavior of system (1.3) with $\tau =5.2<{\tau }_0\approx 5.476$. The endemic equilibrium $E^{\mathrm{<mo>\; *\; </mo>}}$ is asymptotically stable.

the local stability of the equilibria and Hopf bifurcation. We have shown that if $R_0<1$, the disease-free equilibrium is locally asymptotically stable. Further, the sufficient conditions for the stability of the endemic equilibrium are obtained. That is, if $R_0>1$ and the condition (3.7) hold, the endemic equilibrium $E^{\mathrm{<mo>\; *\; </mo>}}$ is asymptotically stable for all $\tau \in \left[\mathrm{0,}{\tau }_0\right)$. As $\tau$ 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 $\tau =6>{\tau }_0\approx 5.476$. Hopf bifurcation occurs from the endemic equilibrium $E^{\mathrm{<mo>\; *\; </mo>}}$.

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).

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Feng, Z., Eppert, A. and Milner, F. (2004) Estimation of Parameters Governing the Transmission Dynamics of Schistosomes. Applied Mathematics Letters, 17, 1105-1112. https://doi.org/10.1016/j.aml.2004.02.002
[2]	Malallah, H., Al-Onaizi, T. and Shuaib, A. (2017) Schistosomiasis as a Cause of Acute Appendicitis in Non-Endemic Areas. International Journal of Surgery, 47, S15-S108. https://doi.org/10.1016/j.ijsu.2017.08.137
[3]	Angeles, J., Goto, Y. and Kirinoki, M. (2011) Human Antibody Response to Thioredoxin Peroxidase-1 and Tandem Repeat Proteins as Immunodiagnostic Antigen Candidates for Schistosoma Japonicum Infection. The American Journal of Tropical Medicine and Hygiene, 85, 674-679. https://doi.org/10.4269/ajtmh.2011.11-0245
[4]	Chistulo, L., Loverde, P. and Engels, D. (2004) Focus: Schistosomiasis. Nature Reviews Microbiology, 2, 2-12.
[5]	Williams, G., Sleigh, A. and Li, Y. (2002) Mathematical Modelling of Schistosomiasis Japonica: Comparison of Control Strategies in the People’s Republic of China. Acta Tropica, 82, 253-262. https://doi.org/10.1016/S0001-706X(02)00017-7
[6]	Wu, Y., Li, M. and Sun, G. (2016) Asymptotic Analysis of Schistosomiasis Persistence in Models with General Functions. Journal of the Franklin Institute, 353, 4772-4784. https://doi.org/10.1016/j.jfranklin.2016.09.012
[7]	Anderson, R. and May, R. (1985) Helminth Infections of Humans: Mathematical Model, Population Dynamics, and Control. Advances in Parasitology, 24, 1-101. https://doi.org/10.1016/S0065-308X(08)60561-8
[8]	Barbour, A. (1996) Modeling the Transmission of Schistosomiasis: An Introductory View. The American Journal of Tropical Medicine and Hygiene, 55, 135-143. https://doi.org/10.4269/ajtmh.1996.55.135
[9]	Feng, Z., Li, C. and Milner, F. (2002) Schistosomiasis Models with Density Dependence and Age of Infection in Snail Dynamics. Mathematical Biosciences, 177-178, 271-286. https://doi.org/10.1016/S0025-5564(01)00115-8
[10]	Qi, L. and Cui, J. (2013) A Schistosomiasis Model with Mating Structure. Abstract and Applied Analysis, 2, 205-215. https://doi.org/10.1155/2013/741386
[11]	Gao, S., He, Y. and Liu, Y. (2013) Field Transmission Intensity of Schistosoma Japonicum Measured by Basic Reproduction Ratio from Modified Barbour’s Model. Parasites & Vectors, 6, 1-10. https://doi.org/10.1186/1756-3305-6-141
[12]	Zhang, P. and Milner, F. (2007) A Schistosomiasis Model with An Age-Structure in Human Hosts and Its Application to Treatment Strategies. Mathematical Biosciences, 205, 83-107. https://doi.org/10.1016/j.mbs.2006.06.006
[13]	Chiyaka, E. and Garira, W. (2009) Mathematical Analysis of the Transmission Dynamics of Schistosomiasis in the Human-Snail Host. Journal of Biological Systems, 17, 397-423. https://doi.org/10.1142/S0218339009002910
[14]	Zhang, X., Gao, S. and Cao, H. (2014) Threshold Dynamics for a Nonautonomous Schistosomiasis Model in a Periodic Environment. Journal of Applied Mathematics and Computing, 46, 305-319. https://doi.org/10.1007/s12190-013-0750-5
[15]	Bhunu, C. and Tchuenche, J. (2010) Modeling the Effects of Schistosomiasis on the Transmission Dynamics of HIV/AIDS. Journal of Biological Systems, 18, 277-297. https://doi.org/10.1142/S0218339010003196
[16]	Chen, Z., Zou, L. and Shen, D. (2010) Mathematical Modelling and Control of Schistosomiasis in Hubei Province, China. Acta Tropica, 115, 119-125. https://doi.org/10.1016/j.actatropica.2010.02.012
[17]	Qi, L. and Cui, J. (2012) Modeling the Schistosomiasis on the Islets in Nanjing. International Journal of Biomathematics, 5, 189-205. https://doi.org/10.1142/S1793524511001799
[18]	Mouhamadou, D., Abderrahman, I. and Mamadou, S. (2014) Global Analysis of a Shistosomiasis Infection with Biological Control. Applied Mathematics and Computation, 246, 731-742. https://doi.org/10.1016/j.amc.2014.08.061
[19]	Feng, Z. and Velasco-Hernández, X. (1997) Competitive Exclusion in a Vector-Host Model for the Dengue Fever. Journal of Mathematical Biology, 35, 523-544. https://doi.org/10.1007/s002850050064
[20]	Chiyaka, E., Magombedze, G. and Mutimbu, L. (2010) Modelling within Host Parasite Dynamics of Schistosomiasis. Computational and Mathematical Methods in Medicine, 11, 255-280. https://doi.org/10.1080/17486701003614336
[21]	Liang, S., Maszle, D. and Spear, R. (2002) A Quantitative Framework for a Multi-Group Model of Schistosomiasis Japonicum Transmission Dynamics and Control in Sichuan, China. Acta Tropica, 82, 263-277. https://doi.org/10.1016/S0001-706X(02)00018-9
[22]	Ruan, S., Xiao, D. and Beier, J. (2008) On the Delayed Ross-Macdonald Model for Malaria Transmission. Bulletin of Mathematical Biology, 70, 1098-1114. https://doi.org/10.1007/s11538-007-9292-z
[23]	Castillo-Chavez, C., Feng, Z. and Xu, D. (2008) A Schistosomiasis Model with Mating Structure and Time Delay. Mathematical Biosciences, 211, 333-341. https://doi.org/10.1016/j.mbs.2007.11.001
[24]	Qi, L. and Cui, J. (2012) The Delayed Barbour’s Model for Schistosomiasis. International Journal of Biomathematics, 5, 117-137. https://doi.org/10.1142/S1793524511001660
[25]	Hassard, B., Kazarino, D. and Wan, Y. (1981) Theory and Application of Hopf Bifurcation. Canbrige Universtiy Press, London, 5. http://booklikes.com/theory-and-application-of-hopf-bifurcation-b-d-hassard/book,3844179