\documentclass[reqno]{amsart} \usepackage{graphicx} \AtBeginDocument{{\noindent\small {\em Electronic Journal of Differential Equations}, Vol. 2007(2007), No. 77, pp. 1--17.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu (login: ftp)} \thanks{\copyright 2007 Texas State University - San Marcos.} \vspace{9mm}} \begin{document} \title[\hfilneg EJDE-2007/77\hfil Bifurcation analysis on SIS epidemic model] {Bifurcation analysis on a delayed SIS epidemic model with stage structure} \author[L. Liu, X. Li, K. Zhuang\hfil EJDE-2007/77\hfilneg] {Li Liu, Xiangao Li, Kejun Zhuang} \address{Li Liu \newline School of Mathematical Sciences, South China Normal University, Guangzhou 510631, China} \email{liuli\_0926@163.com} \address{Xiangao Li \newline School of Mathematical Sciences, South China Normal University, Guangzhou 510631, China} \email{lixg@scnu.edu.cn} \address{Kejun Zhuang \newline School of Mathematical Sciences, South China Normal University, Guangzhou 510631, China} \email{zhkj123@163.com} \thanks{Submitted February 6, 2007. Published May 22, 2007.} \thanks{Supported by grant 10571064 from the National Natural Science Foundation of China} \subjclass[2000]{34K13, 34K18, 34K20} \keywords{SIS model; delay; Hopf bifurcation; stability; periodic solution} \begin{abstract} In this paper, a delayed SIS (Susceptible Infectious Susceptible) model with stage structure is investigated. We study the Hopf bifurcations and stability of the model. Applying the normal form theory and the center manifold argument, we derive the explicit formulas determining the properties of the bifurcating periodic solutions. The conditions to guarantee the global existence of periodic solutions are established. Also some numerical simulations for supporting the theoretical are given. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{proposition}[theorem]{Proposition} \newtheorem{lemma}[theorem]{Lemma} \section{Introduction} In this paper, we study the bifurcation properties of Hopf branches in the stage-structured SIS (Susceptible Infectious Susceptible) model \begin{equation} \label{1.1} \begin{gathered} \dot{x}_1(t)=a_1x_2(t)-c_1x_2(t-\tau)-rx_1(t)-a_2x_1(t)y(t)+b_1y(t),\\ \dot{y}(t)=a_2x_1(t)-b_1y(t)-(r+\alpha)y(t),\\ \dot{x}_2(t)=c_1x_2(t-\tau)-c_2{x_2}^2(t), \end{gathered} \end{equation} with time delay $\tau>0$ as bifurcation parameter, where $x_1(t)$ and $y(t)$ denote the densities of immature susceptible and infectious population, respectively, $x_2(t)$ denotes the density of mature population which do not contact the disease. All coefficients but $\alpha$ are positive constant. System (\ref{1.1}) was proposed by Xiao et al. \cite{Xiao} who showed the stability properties of the system. Stage-structured models have long been and will continue to be of interest to both applied mathematicians and ecologists due to its universal existence and importance \cite{Xiao}. This is not only because they are much more simple than the models governed by partial differential equations but also because they can exhibit phenomena similar to those of partial differential models \cite{Bence}. Thus, the SIS stage-structured models have been widely studied in recent years \cite{Hethcote,Lefere,Liu}. It is well known that models with complexities such as delays can have periodic solutions while the corresponding ordinary differential equation model have globally asymptotically stable equilibrium points, one might expect a more complex behavior for models with stage structure. Thus it is important to carefully analyze epidemic model with stage structure. Our purpose of this paper is to analyze the local and global Hopf bifurcations on the SIS model to give more knowledge to the dynamics of the model. This paper is organized as follows. In the next section, we consider the local Hopf bifurcation and stability of the positive equilibrium. Section 3 presents the properties of Hopf bifurcation on the center manifold. The global existence of multiple periodic solutions is discussed in Section 4. Finally, some numerical simulations are given. \section{Stability of the positive equilibrium and local Hopf bifurcation} Substituting $z(t)$ by $x_1(t)+y(t)$ in system (\ref{1.1}), system (\ref{1.1}) becomes \begin{equation} \label{2.1} \begin{gathered} \dot{z}(t)=a_1x_2(t)-c_1x_2(t-\tau)-rz(t)-\alpha y(t),\\ \dot{y}(t)=y(t)[-(b_1+r+\alpha)-a_2y(t)+a_2z(t)],\\ \dot{x}_2(t)=c_1x_2(t-\tau)-c_2{x_2}^2(t). \end{gathered} \end{equation} Define \begin{equation} \label{C} R_0=\frac{a_2c_1(a_1-c_1)}{rc_2(b_1+r+\alpha)}>1,\quad a_1>c_1, \end{equation} where $R_0$ has been identified for \eqref{2.1} as the basic reproductive number in Xiao et al. \cite{Xiao}. It has been shown that $R_0$ represents the number of secondary infectious caused by an average infective during the course of the disease. If the condition \eqref{C} holds, then system \eqref{2.1} has a unique positive equilibrium $E_*=(z^*,y^*,x_2^*)$, where \begin{equation*} %\label{2.2} z^*=\frac{c_1(a_1-c_1)}{c_2(r+\alpha)} +\frac{\alpha(b_1+r+\alpha)}{a_2(r+\alpha)}, y^*=\frac{c_1(a_1-c_1)}{(r+\alpha)c_2} -\frac{r(b_1+r+\alpha)}{a_2(r+\alpha)}, x_2^*=\frac{c_1}{c_2} \end{equation*} By the linear transform $u_1(t)=z(t)-z^*$, $u_2(t)=y(t)-y^*$, $u_3(t)=x_2(t)-x_2^*$, system \eqref{2.1} becomes \begin{equation} \label{2.3} \begin{gathered} \dot{u}_1(t)=a_1u_3(t)-c_1u_3(t-\tau)-ru_1(t)-\alpha u_2(t),\\ \dot{u}_2(t)=a_2y^*u_1(t)-[-(b_1+r+\alpha)-a_2z^*+2a_2y^*]u_2(t)-a_2u_2^2(t)+a_2u_1(t)u_2(t),\\ \dot{u}_3(t)=c_1u_3(t-\tau)-2c_2{x_2}^*u_3(t)-c_1u_3^2(t). \end{gathered} \end{equation} The associated characteristic equation of system (\ref{2.3}) is \begin{equation} \label{2.4} \lambda^3+m_2\lambda^2+m_1\lambda+m_0+(n_2\lambda^2+n_1\lambda+n_0) e^{-\lambda\tau}=0, \end{equation} where $m_1=rd+a_2\alpha y^*+2c_2x_2^*(r+d)$, $m_2=r+d+2c_2x_2^*$, $m_0=2c_2x_2^*rd+2c_2x_2^*\alpha y*$, $n_2=-c_1,n_1=-c_1(r+d)$, $n_0=-c_1(rd+a_2\alpha y*)$, $d=-(b_1+r+\alpha)-a_2z^*+2a_2y^*$. Obviously, $i\omega$ ($\omega>0$) is a root of (\ref{2.4}) if and only if $\omega$ satisfies \[ - {\omega}^3 i-m_2{\omega}^2+m_1\omega i+m_0+(-n_2\omega^2+n_1\omega i+n_0)(\cos{\omega\tau}-i\sin{\omega\tau})=0, \] separating the real and imaginary parts, we have \begin{equation} \label{2.5} \begin{gathered} m_2\omega^2-m_0=(n_0-n_2\omega^2)\cos{\omega\tau}+n_1\omega\sin{\omega\tau},\\ -\omega^3+m_1\omega=(n_0-n_2\omega^2)\sin{\omega\tau}-n_1\cos{\omega\tau}, \end{gathered} \end{equation} which is equivalent to \begin{equation} \label{2.6} \omega^6+(m_2^2-n_2^2-2m_1)\omega^4+[m_1^2-2m_0m_2-n_1^2+2n_0n_2] \omega^2+m_0^2-n_0^2=0. \end{equation} Let $z=\omega^2$ and denote $p=m_2^2-n_2^2-2m_1$, $q=m_1^2-2m_0m_2-n_1^2+2n_0n_2$ and $R=m_0^2-n_0^2$. Then \eqref{2.6} becomes \begin{equation} \label{2.7} z^3+pz^2+qz+R=0. \end{equation} Now, we introduce the following results which was proved by Song and Yuan \cite{Song}. Denote \begin{equation} \label{2.8} h(z)=z^3+pz^2+qz+R, \end{equation} then \[ \frac{dh(z)}{dz}=3z^2+2pz+q. \] Thus equation \begin{equation} \label{2.9} 3z^2+2pz+q=0 \end{equation} has two real roots \begin{equation} \label{2.10} z_1^*=\frac{-p+\sqrt{\Delta}}{3},\quad z_2^*=\frac{-p-\sqrt{\Delta}}{3} \end{equation} \begin{lemma}[\cite{Song}] \label{lem2.1} For the polynomial equation \eqref{2.7}, we have the following results: \begin{itemize} \item[(i)] If $R<0$,then \eqref{2.7} has at least one positive root. \item[(ii)] If $R \geq 0$ and $\Delta=p^2-3q \leq 0$.then \eqref{2.7} has no positive root. \item[(iii)] If $R \geq 0$ and $\Delta > 0$,then \eqref{2.7} have positive roots if and only if $z_1^*=\frac{-p+\sqrt{\Delta}}{3}>0$ and $h(z^*) \leq 0$. \end{itemize} \end{lemma} \begin{lemma}[\cite{Song}] \label{lem2.2} Suppose that $z_k=\omega_k^2,n_0 \neq n_2\omega_k^2$, and $h'(z_k)\neq0$, then \[ \frac{d(\mathop{\rm Re}(\lambda(\tau_k^{(j)}))}{d\tau}\neq 0 \] and $d(\mathop{\rm Re}(\lambda(\tau_k^{(j)}))/d\tau$ and $h'(z_k)$ have the same sign. \end{lemma} Now we study the characteristic equation (\ref{2.4}) of system (\ref{2.3}). As we know $m_i,n_i$ ($i=0,1,2,3$), we can obtain $p=m_2^2-n_2^2-2m_1$, $q=m_1^2-2m_0m_2-n_1^2+2n_0n_2$, $R=m_0^2-n_0^2$. Thus \[ \Delta=p^2-3q ,\quad h(z)=z^3+pz^2+qz+R,\quad z_1^*=\frac{-p+\sqrt{\Delta}}{3}. \] Note that when $\tau=0$, \eqref{2.4} becomes \begin{equation} \label{2.11} \lambda^3+m_2\lambda^2+m_1\lambda+m_0=0 \end{equation} Applying Lemma \ref{lem2.1} to \eqref{2.4}, we obtain the following lemma. \begin{lemma} \label{lem2.3} For the third degree transcendental \eqref{2.4}, we have \begin{itemize} \item[(i)] if $R \geq 0$ and $\Delta \leq 0$, then all roots with positive real parts of \eqref{2.4} has the same sum to those of the polynomial \eqref{2.11} for all $\tau \geq 0$; \item[(ii)] if either $R<0$ or $R \geq 0,\Delta>0,z_1^*>0$ and $h(z_1^*)\leq 0$, then all roots with positive real parts of \eqref{2.4} have the same sum to those of the polynomial \eqref{2.11} for $\tau \in [0,\tau_0)$. \end{itemize} \end{lemma} To state the next lemma we define the hyothesis \begin{itemize} \item[(H)] $ m_2m_1-m_0>0$ and $m_0>0$. \end{itemize} \begin{lemma} \label{lem2.4} Suppose that $R<0$ and the condition holds. Then all roots of \eqref{2.4} have negative real parts when $\tau \in [0,\tau_0)$. \end{lemma} \begin{proof} From (\ref{2.11}), by the Routh-Huriwz criterion, the real parts of all roots of \eqref{2.11} are negative if and only if (H) is satisfied. Then using the results in Lemma \ref{lem2.3}, we complete the proof. \end{proof} Suppose that \eqref{2.7} have positive roots, without loss of generality,we assume that it has three positive roots $z_1,z_2,z_3$. Then \eqref{2.6} has three positive roots $\omega_1=\sqrt{z_1}$, $\omega_2=\sqrt{z_2}$, $\omega_3=\sqrt{z_3}$. From (\ref{2.5}), we have \[ \cos(\omega\tau)=\frac{n_1\omega^2(\omega^2-m_1) -(m_2\omega^2-m_0)(n_2\omega^2-n_0)}{(n_2\omega^2-n_0)+n_1^2\omega^2}. \] Thus, if we denote \begin{equation} \label{2.12} \tau_k^{(j)}=\frac{1}{\omega_k}\{\cos^{-1} (\frac{(n_1-m_2n_2)\omega_k^4+(m_0n_2+m_2n_0-m_1n_1)\omega_k^2-m_0n_0}{(n_2\omega_k^2-n_0)+n_1^2\omega_k^2})+2j\pi\}, \end{equation} where $k=1,2,3$; $j=0,1,\dots$, then $\pm\omega_k$ is a pair of purely imaginary roots of \eqref{2.4} with $\tau_k^{(j)}$. Define \begin{equation} \label{2.13} \tau_0=\tau_{k0}^{(0)}=min\{\tau_k^{(0)}\},\quad \omega_0=\omega_{k0} \end{equation} Let $\lambda(\tau)=\alpha(\tau)+i\omega(\tau)$ be the root of \eqref{2.4} near $\tau=\tau_k^{(j)}$ satisfying $\alpha(\tau_k^{(j)})=0,\omega(\tau_k^{(j)})=\omega_k$. Applying Lemmas \ref{lem2.1}--\ref{lem2.4} to \eqref{2.4},we have the following theorem about the stability of the positive equilibrium of system \eqref{2.1} and Hopf bifurcations. \begin{theorem} \label{thm2.5} Let $\tau_k^{(j)}$ and $\omega_0,\tau_0$ be defined by \eqref{2.12} and \eqref{2.13}, respectively. Suppose that the condition \eqref{C} and (H) hold. (i) if $R \geq 0$ and $\Delta \leq 0$, then the positive equilibrium $E_*$ of system \eqref{4.2} is asymptotically stable for all $\tau \geq 0$, that is \eqref{2.1} is absolutely stable; (ii) if either $R<0$ or $R \geq 0,\Delta>0,z_1^*>0$ and $h(z_1^*)\leq 0$, then the positive equilibrium $E_*$ of system \eqref{2.1} stable for $\tau \in [0,\tau_0)$, and unstable when $\tau>\tau_0$; (iii) if the conditions of (ii) are satisfied,and $p^2<3q$, then system \eqref{2.1} undergoes a Hopf bifurcation at the equilibrium $E_*$ when $\tau=\tau_k^{(j)}$, where $\tau_k^{(j)}$ is defined by \eqref{2.12}. \end{theorem} \section{Direction and stability of the Hopf bifurcation} From section 2, we obtained a set of conditions for system \eqref{2.1} to undergo Hopf bifurcation. In this section, we shall derive the explicit formulae determining the direction, stability and period of the bifurcating non-trivial periodic solutions at $\tau=\tau_k^{(j)}$. To accomplish this, we use the normal form and center manifold theory developed by Hassard et al. \cite{Hassaard}. Let $u_1(t)=z-z^*$, $u_2(t)=y-y^*$, $u_3(t)=x_2(t)-x_2^*$, $\bar{u}_i(t)=u_i(\tau t)$, $\tau=\tau_k+\mu$ and drop the bar for simplification, system \eqref{2.1} is transformed into an FDE in $C=C([-1,0],R^3)$ as \begin{equation}\label{3.1} \dot{u}(t)=L_\mu(u_t)+f(\mu,u_t), \end{equation} where $u(t)=(u_1(t),u_2(t),u_3(t))^T \in R^3$, and $L_\mu:C \to R$, $f:R \times C \to R$ are given respectively by $$ L_\mu\phi=(\tau_k+\mu)B_1\phi(0)+(\tau_k+\mu)B_2\phi(-1), $$ where $B_1$ and $B_2$ are defined as \[ B_1=\begin{pmatrix} -r & -\alpha & a_1 \\ a_2y^* & -d & 0\\0 & 0 & 2c_1 \end{pmatrix},\quad B_2=\begin{pmatrix} 0 & 0 & -c_1 \\ 0 & 0 & 0\\ 0 & 0 & c_1 \end{pmatrix} \] and \begin{equation} \label{3.2} f(\mu,\phi)=(\tau_k+\mu) \begin{pmatrix} 0 \\-a_2\phi_2^2(0)+a_2\phi_1(0)\phi_2(0) \\ -c_1^2\phi_3^2(0). \end{pmatrix} \end{equation} By the Riesz representation theorem, there exist a matrix whose components are bounded variation functions $\eta(\theta,\mu)$ in $\theta \in [-1,0]$ such that \begin{equation} \label{3.3} L_\mu\phi=\int_{-1}^0 d\eta(\theta,\mu)\phi(\theta) \quad\mbox{for} \quad \phi \in C. \end{equation} In fact, if we choose \begin{equation}\label{3.4} \eta(\theta,\mu)= \begin{cases} B_1\delta(\theta), & \theta=0, \\ -B_2\delta(\theta+1), & \theta \in [-1,0), \end{cases} \end{equation} where $\delta$ is a Dirac Delta function, then (\ref{3.3}) is satisfied. For $ \phi = ( \phi_1, \phi_2)^T \in C$, define: \begin{gather*} A(\mu)\phi= \begin{cases} \frac{d\phi(\theta)}{d\theta}, & \theta \in [-1,0), \\ \int_{-1}^0 d\eta(\mu,s)\phi(s), & \theta=0, \end{cases} \\ R(\mu)\phi= \begin{cases} 0, & \theta \in [-1,0),\\ f(\mu,\phi), & \theta=0. \end{cases} \end{gather*} Hence, we can rewrite (\ref{3.1}) as \begin{equation} \label{3.5} \dot{u}_t=A(\mu)u_t+R(\mu)u_t, \end{equation} where $u=(u_1,u_2,u_3)^T$, $u_t=u(t+\theta)$, for $\theta \in [-1,0]$. For $\psi \in C^1([0,1],(R^3)^*)$, define \begin{equation}\label{3.6} A^*\psi(s)= \begin{cases} -\frac{d\psi(s)}{ds}, & s \in (0,1], \\ \int_{-1}^0 d\eta(t,0)\psi(-t) ,& s=0. \end{cases} \end{equation} For $\phi \in C^1([-1,0],R^2)$, and $\psi \in C^1([0,1],(R^2)^*)$, define the bilinear form \begin{equation} \label{3.7} \langle \psi,\phi\rangle =\bar{\psi}(0)\phi(0) -\int_{-1}^0\int_{\xi=0}^\theta \bar{\psi}(\xi-\theta)d\eta(\theta)\phi(\xi)d\xi, \end{equation} where $\eta(\theta)=\eta(\theta,0)$. Obviously $A^*$ and $A$ are adjoint operators. By the results in section 2, we know that $\pm i\tilde{\tau}\omega$ are eigenvalues of $A(0)$. Thus, they are also eigenvalues of $A^*$. Let $q(\theta)=(1,\beta,\gamma)^Te^{i\theta\omega{\tau_k}}$ is the eigenvector of $A(0)$ corresponding to $i\tilde{\tau}\omega$ and $q^*(s)=D(1,\beta^*,\gamma^*)e^{is\omega{\tau_k}}$ is the eigenvector of $A^*$ corresponding to $-i\tilde{\tau}\omega$. Moreover, by a direct computation, we can show that \[ q(\theta)=\begin{pmatrix} 1\\ \frac{i\omega+d}{a_2y^*}\\ \frac{a_2y^*r+\alpha b_1+\alpha r+\alpha^2+2\alpha a_2y^*-\alpha a_2z^*+(\omega\alpha+\omega a_2y^*)i}{a_2y^*(a_1-c_1e^{i\omega \tau_k})} \end{pmatrix} e^{i\theta\omega{\tau_k}}, \] and \[ q^*(s)=D\begin{pmatrix} 1\\ \frac{a_2y^*(c_1e^{i\omega \tau_k}+2c_1)}{-r(c_1e^{-i\omega \tau_k}-2c_1)+i\omega}\\ \frac{-i\omega+dc_1e^{i\omega \tau_k}-2c_1d}{d(c_1e^{i\omega \tau_k}-a_1)}+\frac{\alpha a_2y^*(c_1e^{i\omega \tau_k}-2c_1)^2}{r(c_1e^{-i\omega \tau_k}-2c_1)-i\omega} \end{pmatrix} ^Te^{i\omega \tau_k s}. \] To assure $=1$, we need to determine the value of $D$. From (\ref{3.7}),we have \begin{align*} \langle q^*(s),q(\theta)\rangle &=\bar{D}(1,\beta^*,\gamma^*)(1,\beta,\gamma)^T\\ &\quad -\int_{-1}^0\int_{\xi=0}^\theta (1,\bar{\beta^*},\bar{\gamma^*}) e^{-i(\xi-\theta)\omega\tau_k}d\eta(\theta)(1,\beta,\gamma)^T e^{i\xi\omega\tau_k}d\xi\\ &=\bar{D}\{1 + \beta \beta^*+\gamma \gamma^* - \int_{-1}^0(1, \bar{\beta^*}, \bar{\gamma^*}) \theta e ^{i \theta \omega \tau_k} d\eta(\theta)(1,\beta,\gamma)^T \}\\ &=\bar{D}\{1+\beta\beta^*+\gamma\gamma^*-\tau_ke^{-i\omega\tau_k} (c_1r-c_1r\bar{r}^*)\}. \end{align*} Thus, we can choose $D$ as \[ D=\frac{1}{1+\beta\beta^*+\gamma\gamma^*-\tau_ke^{i\omega\tau_k} (c_1r-c_1r\bar{r}^*)}. \] Using the same notation as in Hassard et al. \cite{Hassaard}, we first compute the coordinates to describe the center manifold $\textbf{C}_0$ at $\mu=0$. Let $u_t$ be the solution of (\ref{3.1}) when $\mu=0$. Define \begin{equation} \label{3.8} z(t)=\langle q^*,u_t\rangle ,\quad w(t,\theta)=u_t-2\mathop{\rm Re}\{z(t)q(\theta)\}. \end{equation} On the center manifold $\textbf{C}_0$ we have \[ w(t,\theta)=w(z(t),\bar{z}(t),\theta), \] where \begin{equation} \label{3.9} w(z(t),\bar{z}(t),\theta)=w_{20}(\theta)\frac{z^2}{2}+w_{11}(\theta)z\bar{z} +w_{02}\frac{\bar{z}^2}{2}+\dots, \end{equation} $z$ and $\bar{z}$ are local coordinates for center manifold $\textbf{C}_0$ in the direction of $q^*$ and $\bar{q}^*$. Note that $w$ is real if $u_t$ is real.We consider only real solutions.For solution $u_t \in \textbf{C}_0$ of \eqref{3.1}, since $\mu=0$, \begin{equation} \label{3.10} \begin{aligned} \dot{z}(t) &=i\omega\tilde{\tau}z+\langle \bar{q^*}(\theta),f(0,w(z,\bar{z},\theta) +2\mathop{\rm Re}\{z(t)q(\theta)\})\rangle \\ &=i\omega\tilde{\tau}z+\bar{q^*}(0)f_0(z,\bar{z}); \end{aligned} \end{equation} that is, \begin{equation} \label{3.11} \dot{z}(t)=i\omega\tilde{\tau}z(t)+g(z,\bar{z}), \end{equation} where \begin{equation} \label{3.12} g(z,\bar{z})= \bar{q^*}(0)f_0(z,\bar{z}) = g_{20}\frac{z^2}{2}+g_{11}z\bar{z}+g_{02}\frac{\bar{z}^2}{2}+g_{21} \frac{z^2\bar{z}}{2}+\dots. \end{equation} Then from (\ref{3.12}), we have \begin{align*} g(z,\bar{z}) &=\bar{q}^*(0)f_0(z,\bar{z})=\bar{q}*(0)f(0,u_t)\\ &=\tau_k \bar{D}(1,\bar{\beta^*},\bar{\gamma^*}) \begin{pmatrix} 0 \\ -a_2u_{2t}(0)^2+a_2u_{1t}(0) \\ -c_1u_{3t}^2(0) \end{pmatrix} \\ &=\tau_k \bar{D}[-a_2u_{2t}^2(0)+a_2u_{1t}(0)u_{2t}(0)\bar{\beta^*} -\bar{\gamma^*}c_1u_{3t}^2(0)]\\ &=\tau_k \bar{D} \bar{\beta^*} \{-a_2[\beta z+\bar{\beta} \bar{z} +w_{20}^{(2)}(0) \frac{z^2}{2}+w_{11}^{2}(0) z \bar{z} +w_{02}^{2}(0)\frac{\bar{z}^2}{2}+o(|(z,\bar{z})|^3)]^2\\ &\quad+a_2[z+\bar{z}+w_{20}^{(1)}(0)\frac{z^2}{1}+w_{11}^{1}(0)z\bar{z} +w_{02}^{1}(0)\frac{\bar{z}^2}{2}+o(|(z,\bar{z})|^3)]\\ &\quad\times [\beta z + \bar{\beta} \bar{z}\ + w_{20}^{(2)}(0) \frac{z^2}{2} +w_{11}^{(2)}(0) z \bar{z}+w_{02}^{2}(0) \frac{\bar{z}^2}{2} + o(| (z,\bar{z}) |^3)]\}\\ & -\tau_k \bar{D} \bar{\gamma^*}c_1[rz+\bar{r} \bar{z}+w_{20}^{(3)}(0)\frac{z^2}{3}+w_{11}^{3}(0)z\bar{z} +w_{02}^{3}(0)\frac{\bar{z}^2}{2}+o(|(z,\bar{z})|^3]^2. \end{align*} Comparing the coefficients with (\ref{3.12}), we obtain \begin{equation} \label{3.13} \begin{aligned} g_{20}&=2\tau_k \bar{D}\bar{\beta^*}a_2(\beta-\beta^2)- 2\tau_k \bar{D} \bar{\gamma^*}\gamma^2,\\ g_{11}&=2\tau_k \bar{D}\bar{\beta^*}a_2(-{| \beta |}^2+2\mathop{\rm Re}\beta)-2\tau_k \bar{D}\bar{\gamma^*}c_1{|\bar{\gamma}^2|}^2,\\ g_{02}&=2\tau_k \bar{D}\bar{\beta^*}a_2(\bar{\beta}-\bar{\beta}^2) - 2\tau_k \bar{D}\bar{\gamma^*}\gamma^2,\\ g_{21}&=\tau_k \bar{D}\bar{\beta^*}a_2(-w_{20}^{(2)}(0)\bar{\beta} -w_{11}^{2}(0)\bar{\beta}+w_{11}^{2}(0)\\ &\quad +w_{20}^{(2)}(0)+w_{11}^{1}(0) \bar{\beta}+w_{11}^{1}(0)\beta) -\tau_k \bar{D}\bar{\gamma^*}c_1(\gamma\omega_{11}^{(3)} +\bar{\gamma}\omega_{20}^{(3)}). \end{aligned} \end{equation} We still need to compute $w_{20}(\theta)$ and $w_{11}(\theta)$. From (\ref{3.5}) and (\ref{3.8}) ,we have \begin{equation} \label{3.14} \begin{aligned} \dot{w}=\dot{u}_t-\dot{z}q-\dot{\bar{z}}\dot{\bar{z}} &=\begin{cases} Aw-2\mathop{\rm Re}{\bar{q}^*(0)f_0q(\theta)}, & \theta \in [-1,0),\\ Aw-2\mathop{\rm Re}{\bar{q}^*(0)f_0q(0)}+f_0, & \theta=0, \end{cases} \\ &=Aw+H(z,\bar{z},\theta), \end{aligned} \end{equation} where \begin{equation} \label{3.15} H(z,\bar{z},\theta)=H_{20}(\theta)\frac{z^2}{2}+H_{11}(\theta)z\bar{z} +H_{02}(\theta)\frac{\bar{z}^2}{2}. \end{equation} Expanding the above series and comparing the coefficients, we obtain \begin{equation} \label{3.16} (A-2i\omega\tau_k)w_{20}(\theta)=-H_{20}(\theta),\quad Aw_{11}(\theta)=-H_{11}(\theta). \end{equation} From (\ref{3.14}) we know that for $\theta \in [-1,0)$, \begin{equation} \label{3.17} H(z,\bar{z},\theta)=-\bar{q}^*f_0q(\theta)-q^*(0)\bar{f}_0\bar{q}(\theta) =-g(z,\bar{z})q(\theta)-\bar{g}(z,\bar{z})\bar{q}(\theta). \end{equation} Comparing the coefficients with (\ref{3.15}), we obtain \begin{gather} \label{3.18} H_{20}=-g_{20}q(\theta)-\bar{g}_{02}\bar{q}(\theta), \\ \label{3.19} H_{11}=-g_{11}q(\theta)-\bar{g}_{11}\bar{q}(\theta). \end{gather} From (\ref{3.16}) and (\ref{3.18}), we obtain \[ \dot{w}_{20}=2i\omega\tilde{\tau}w_{20}(\theta)+g_{20}q(\theta)+\bar{g}_{02}\bar{q}(\theta). \] Note that $q(\theta)=(1,\alpha,\beta)^Te^{i\theta\omega\tau_k}$, hence \begin{equation} \label{3.20} w_{20}(\theta)=\frac{ig_{20}}{\omega\tau_k}q(0)e^{i\omega\tau_k\theta} +\frac{i\bar{g}_{02}}{3\omega\tau_k}\bar{q}(0)e^{-i\omega\tau_k\theta}. +E_1e^{2i\omega\tau_k\theta}. \end{equation} Similarly, from (\ref{3.16}) and (\ref{3.19}), we have \begin{equation} \label{3.21} \begin{gathered} \dot{w}_{11}=g_{11}q(\theta)+\bar{q}_{11}\bar{q}(\theta), \\ \omega_{11}(\theta)=-\frac{ig_{20}}{\omega, \tau_k}q(0)e^{i\omega\tau_k\theta} +\frac{i\bar{g}_{11}}{\omega\tau_k}\bar{q}(0)e^{-i\omega\tau_k\theta}+E_2, \end{gathered} \end{equation} where $E_1$ and $E_2$ are both 3-dimensional vectors, and can be determined by setting $\theta=0$ in $H$. In fact, from the definition of $A$ and (\ref{3.16}) that \begin{gather} \label{3.22} \int_{-1}^0d\eta(\theta)w_{20}(\theta)=2i\omega\tau_kw_{20}(0)-H_{20}(0), \\ \label{3.23} \int_{-1}^0d\eta(\theta)w_{11}(\theta)=-H_{11}(0), \end{gather} where $\eta(\theta)=\eta(0,\theta)$. From (\ref{3.14}), we have \begin{gather} \label{3.24} H_{20}=-g_{20}q(0)-\bar{g}_{02}\bar{q}(0)+2\begin{pmatrix} 0 \\ -a_2\beta^2+a_2\beta\\ -c_1\gamma^2 \end{pmatrix}, \\ \label{3.25} H_{11}=-g_{11}q(0)-\bar{g}_{11}\bar{q}(0)+2\begin{pmatrix} 0 \\ -a_2|\beta|^2+a_2\mathop{\rm Re}\beta\\ -c_1|\gamma|^2 \end{pmatrix}. \end{gather} Substituting (\ref{3.20}) and (\ref{3.24}) into (\ref{3.22}), and noticing that \begin{gather*} (i \omega\tau_k I - \int_{-1}^0 e^{i\theta\omega\tau_k}d\eta(\theta))q(0)=0, \\ (-i \omega \tau_k I -\int_{-1}^0 e^{-i\theta\omega\tau_k}d\eta(\theta))\bar{q}(0)=0, \end{gather*} we obtain \begin{equation} \label{3.26} (2i\omega\tau_k I-\int_{-1}^0e^{2i\theta\omega\tau_k}d\eta(\theta))E_1=2\tau_k \begin{pmatrix} 0 \\ -a_2\beta^2+a_2\beta\\ -c_1\gamma^2 \end{pmatrix}; \end{equation} that is, \begin{equation} \label{3.27} \begin{pmatrix} 2i\omega+r & \alpha & -a+c_1e^{-2i\omega\tau_k}\\ -a_2y^* & 2i\omega+d & 0 \\ 0 & 0 & 2i\omega c_1-c_1e^{-2i\omega\tau_k} \end{pmatrix} E_1=2 \begin{pmatrix} 0 \\ -a_2\beta^2+a_2\beta\\ -c_1\gamma^2 \end{pmatrix}. \end{equation} Similarly, substituting (\ref{3.21}) and (\ref{3.25}) into (\ref{3.23}), we have \[ \int_{-1}^0d\eta(\theta)E_2=2\begin{pmatrix} 0 \\ -a_2|\beta|^2+a_2\mathop{\rm Re}\beta\\ -c_1|\gamma|^2 \end{pmatrix}, \] which leads to \begin{equation} \label{3.28} \begin{pmatrix} r & \alpha & -a+c_1\\ -a_2y^* & d & 0 \\ 0 & 0 & c_1 \end{pmatrix} E_2=2\begin{pmatrix} 0 \\ -a_2|\beta|^2+a_2\mathop{\rm Re}\beta\\ -c_1|\gamma|^2 \end{pmatrix}. \end{equation} Thus, $g_{21}$ can be computed. Then, we can compute the following values: \begin{equation}\label{3.29} \begin{gathered} c_1(0)=\frac{i}{2\omega\tilde{\tau}}(g_{11}g_{20}-2{| g_{11} |}^2 -\frac{{| g_{02} |}^2}{3})+\frac{g_{21}}{2},\\ \mu_2=\frac{\mathop{\rm Re}{c_1(0)}}{\mathop{\rm Re}{\lambda'_0(\tilde{\tau})}},\\ \beta_2 =2\mathop{\rm Re}{c_1(0)},\\ T_2 =-\frac{Im(c_1(0))+\mu_2Im(\lambda'_0(\tilde{\tau}))}{\omega}, \end{gathered} \end{equation} which determine the quantities of bifurcating periodic solutions in the center manifold at the critical value $\tau_k$. we know that (see Hassard et.al. \cite{Hassaard}) (i) $\mu_2$ determines the directions of the Hopf bifurcation:if $\mu_2>0$ ($<0$), then the Hopf bifurcation is supercritical(subcritical) meaning that bifurcated periodic solution exists for $\tau > \tau_k (\tau<\tau_k)$; (ii) $\beta_2$ determines the stability of the bifurcating periodic solutions:the bifurcating periodic solutions in the center manifold are stable(unstable)if $\beta_2<0(\beta_2>0)$; and $T_2$ determines the period of the bifurcating periodic solutions:the period increase(decrease)if $T_2>0$ ($T_2<0$). \section{Existence of global Hopf bifurcation} In this section, we will study the global existence of non-trivial periodic solution using global Hopf bifurcation theorem given by Wu \cite{Wu}. At first, we introduce some notation: Let$u_1(t)=z(t)$, $u_2(t)=y(t)$, $u_3(t)=x_2(t)$, then system \eqref{2.1} can be written as \begin{equation} \label{4.1} \dot{u}(t)=F(u_t,\tau,p), \end{equation} where $u(t)=(u_1(t),u_2(t),u_3(t))\in R^3$, denote $u^*=(u_1^*,u_2^*,u_3^*)$ as the positive equilibrium and $u_t=(u_{1t},u_{2t},u_{3t})$. Let $X=C([-\tau,0],R^3)$, $N=\{(\bar{u},\tau,p);F(\bar{u},\bar{\tau},\bar{p})=0\}$, $$ \Sigma=Cl\{(u,\tau,p)\in X\times R\times R^+;u \mbox{ is a $p$-periodic solution of \eqref{2.1}}\}, $$ and $\ell(u^\ast,\tau^{(j)}_k,\frac{2\pi}{\omega_0})$ denotes the connected component through $(u^\ast,\tau^{(j)}_k,\frac{2\pi}{\omega_0})$ in $\Sigma$. First we present a global Hopf bifurcation result from Wu \cite[Theorem 3.3]{Wu}. \begin{proposition} \label{prop4.1} Assume that (A1)-(A6) hold. Let $ \Sigma(F)=Cl\{(x,\alpha,p)\in X\times R\times R^+;u \mbox{is a $p$-periodic solution of \eqref{4.1}}\}$ which is a sub set of $X\times\mathbb{R}\times\mathbb{R}$. Let $N(F)=\{(\hat{x},\alpha,p); F(\hat{x},\alpha,p)=0\}$. Let $C(\hat{x}_0,\alpha_0,p_0)$ be the connected component in $\Sigma(F)$ of an isolated center $(\hat{x}_0,\alpha_0,p_0)$. Then either \begin{itemize} \item[(i)] $C(\hat{x}_0,\alpha_0,p_0)$ is unbounded, or \item[(ii)] $C(\hat{x}_0,\alpha_0,p_0)$ is bounded, $C(\hat{x}_0,\alpha_0,p_0)\cap N(F)$ is finite and \begin{equation} \label{4.2} \sum_{(\hat{x},\alpha,p)\in C(\hat{x}_0,\alpha_0,p_0)\cap N(F)}\gamma_m(\hat{x},\alpha,p)= 0 \end{equation} for all $m=1,2,\ldots$,where $\gamma_m(\hat{x},\alpha,p)$ is the mth crossing number of $(\hat{x},\alpha,p)$ if $m \in J(\hat{x}_0,\alpha_0,p_0)$ or it is zero if otherwise. \end{itemize} \end{proposition} Another technical issue when applying proposition \ref{prop4.1} is to prove that (\ref{4.1}) with $\tau=0$ has non-constant periodic solutions. This will be done by applying a high-dimensional Bendixson's criterion of Li and Muldowney \cite{Muldowney}, which we briefly summarize in the following. Consider a system of ordinary differential equations \begin{equation} \label{4.3} \dot{x}=f(x),\quad x \in \mathbb{R}^n,\quad f \in C^1 \end{equation} for any finite $n$. Denote \begin{equation} \label{4.4} z'(t)=\frac{\partial f^{[2]}}{\partial x}(x(t,x_0))z(t) \end{equation} as the second compound equation with respect to a solution $x(t,x_0)\in D$ to (\ref{4.5}) (see Fiedler \cite{Fiedler} and Muldowney \cite{Muldowney}). \begin{proposition} \label{prop4.2} Let $D \subset \mathbb{R}^n$ be a simply connected region. Assume that the family of linear systems \[ z^{\prime}(t)=\frac{\partial f^{[2]}}{\partial x}(x(t,x_0))z(t),\quad x_0\in D \] is equi-uniformly asymptotically stable. Then (a) $D$ contains no simple closed invariant curves including periodic orbits, homoclinic orbits,heteroclinic cycles; (b) each semi-orbit in D converges to a single equilibrium. In particular, if $D$ is positively invariant and contains an unique equilibrium $\bar{x}$, the $\bar{x}$ is globally asymptotically stable in $D$. \end{proposition} To investigate the global properties of the solution, now we present some basic results of system (\ref{4.1}). \begin{lemma} \label{lem4.1} Non-constant periodic solutions of (\ref{4.1}) are uniformly bounded. \end{lemma} \begin{proof} For periodic functions $u_1(t)$, $u_2(t)$, $u_3(t)$, we define \begin{equation} \label{4.5} \begin{gathered} u_1(\xi_1)=\min\{u_1(t)\},u_1(\eta_1)=\max\{u_1(t)\},\\ u_2(\xi_2)=\min\{u_2(t)\},u_2(\eta_2)=\max\{u_2(t)\},\\ u_3(\xi_3)=\min\{u_2(t)\},u_3(\eta_3)=\max\{u_3(t)\}, \end{gathered} \end{equation} Let $(u_1(t)$, $u_2(t)$, $u_3(t))$ be a nonconstant periodic solution of (\ref{4.1}). Then we know either $u_i(t)\equiv 0$ or $u_i(t)\neq 0$ ($i=1,2,3$). Thus, there are four cases to be considered. \noindent\textbf{Case 1.} When $u_3(t)>0$, we have, from the third equation of (\ref{4.1}), \[ 0 = c_1u_3(\eta_3-\tau)-c_2u_3^2(\eta_3) \leq c_1u_3(\eta_3)-c_2u_3^2(\eta_3), \] which leads to \[ u_3(\eta_3) \leq \frac{c_1}{c_2}. \] (1) $u_1(t)>0$, $u_2(t)<0$. From the second equation of (\ref{4.1}), \[ 0 = -(b_1+r+\alpha)-a_2u_2(\xi_2)+a_2u_1(\xi_2) > -(b_1+r+\alpha)-a_2u_2(\xi_2), \] which implies \[ u_2(\xi_2) > -\frac{b_1+r+\alpha}{a_2}. \] From the first equation of (\ref{4.1}), we get \[ 0 = a_1u_3(\eta_1)-c_1u_3(\eta_1-\tau)-ru_1(\eta_1)-\alpha u_2(\eta_1) \leq a_1 \cdot \frac{c_1}{c_2}-ru_1(\eta_1)+\alpha \cdot \frac{b_1+r+\alpha}{a_2} \] which induces \[ u_1(\eta_1) \leq \frac{a_1c_1}{rc_2}+\frac{\alpha(b_1+r+\alpha)}{ra_2}\,. \] \noindent (2) $u_1(t)>0, u_2(t)>0$, \[ \dot{u}_1(t)=a_1u_3(t)-c_1u_3(t-\tau)-ru_1(t)-\alpha u_2(t) \leq a_1\frac{c_1}{c_2}-ru_1(t). \] Consider the comparison equation \begin{equation} \label{4.6} v(t)= a_1\frac{c_1}{c_2}-rv(t), \end{equation} since ${a_1c_1}/{rc_2}$ is the unique positive equilibrium of \eqref{4.6}, and it is globally asymptotically stable. Let $u_1(0) \leq v(0)$, then by comparison principle, $u_1(t) \leq v(t)$ for $t \geq 0$. Furthermore, \[ \lim_{t \to \infty}\sup u_1(t) \leq \lim_{t \to \infty}\sup v(t) = \frac{a_1c_1}{rc_2} \] Hence, there exists a $M>0$, such that for $t>M$, \[ u_1(t) \leq \frac{a_1c_1}{rc_2}\,. \] For $t\in[0,M]$, there exists a $M' > 0$, such that $u_1(t)\leq M'$. Then, we can choose $M_1=\max\{M',\frac{a_1c_1}{rc_2}\}$, for $t \geq 0,u_1(t) \leq M_1$. Similarly, we can prove there exists a $M_2>0$, such that $u_2(t) \leq M_2$. \begin{align*} \dot{u}_2(t)&=u_2(t)[-(b_1+r+\alpha)-a_2u_2(t)+a_2u_1(t)]\\ &\leq [-(b_1+r+\alpha)+a_2M_1-a_2u_2(t)]u_2(t) \end{align*} Let $N=-(b_1+r+\alpha)+a_2M_1$, and the comparison equation is \[ \dot{v}(t)=v(t)(N-a_2v(t)). \] Then similar to the proof of $u_1(t)$, there exists a $M_2>0$, such that $u_2(t) \leq M_2$. This completes the proof. \noindent (3) $u_1(t)<0$, $u_2(t)>0$. From the second equation of (\ref{4.1}), \[ 0 = -(b_1+r+\alpha)-a_2u_2(\eta_2)+a_2u_1(\eta_2) < -(b_1+r+\alpha)-a_2u_2(\eta_2), \] which means \[ u_2(\eta_2) < -\frac{b_1+r+\alpha}{a_2} < 0. \] Obviously, it is a contradiction. Also, \[ 0 = a_1u_3(\xi_1)-c_1u_3(\xi_1-\tau)-ru_1(\xi_1)-\alpha u_2(\xi_1) \geq -c_1\frac{c_1}{c_2}-ru_1(\xi_1)+\frac{\alpha(b_1+r+\alpha)}{a_2}, \] which leads to \[ u_1(\xi_1) \geq \frac{\alpha(b_1+r+\alpha)}{ra_2}-\frac{c_1^2}{ra_2}. \] Clearly, it is also a contradiction. Thus, there are no nontrivial periodic solution to (\ref{4.1}) in this case. \noindent (4) $u_1(t)<0$, $u_2(t)<0$. From the first and second equation of (\ref{4.1}), we get \[ \dot{u}_1(t)=a_1u_3(t)-c_1u_3(t-\tau)-ru_1(t)-\alpha u_2(t) \geq -c_1\frac{c_1}{c_2}-ru_1(t), \] and \[ \dot{u}_2(t)=u_2(t)[-(b_1+r+\alpha)-a_2u_2(t)+a_2u_1(t)] > u_2(t)[-(b_1+r+\alpha)-a_2u_2(t)]. \] Then we can use the comparison principle similar to (2), we get the results. \noindent\textbf{Case 2.} When $u_3(t)<0$, from the third equation of (\ref{4.1}), we have \[ 0 = c_1u_3(\xi_3-\tau)-c_2u_3^2(\xi_3) \geq c_1u_3(\xi_3) - c_2u_3^2(\xi_3) \] which implies \[ u_3(\xi_3) \geq \frac{c_1}{c_2}. \] Under the assumption $u_3(t)<0$, by the similar analysis of Case 1, we can prove the nontrivial periodic solution of (\ref{4.1}) is uniformly bounded. \noindent\textbf{Case 3.} $u_i(t)\equiv 0,u_j(t)\neq 0$ ($i=1,2,3, j=1,2,3, i \neq j$). \noindent\textbf{Case 4.} $u_i(t) \neq 0,u_j(t)\equiv 0$ ($i=1,2,3, j=1,2,3, i \neq j$). The method for proving these two cases is similar to that of Cases 1 and 2. Mainly use the comparison principle and analyze the sign of $u_i(t)(i=1,2,3)$. We omit their proofs. \end{proof} For the next lemma, we define the hypothesis \begin{itemize} \item[(H')] There exist $N_1, N_2$ such that \begin{equation} \label{4.7} \begin{aligned} &\sup_{x\in R}\{[-r-(b_1+r+\alpha)]+N_1(c_1-a_1)+2a_2|u_2|,\\ &(-r+c_1)|u_2(t)|+2N_2c_2|u_2(t)|+N_2\alpha,\\ &(c_1-b_1-r-\alpha)+(a_2/N_2)|u_2(t)|+2c_2|u_3(t)|+2a_2|u_2(t)|\}<0. \end{aligned} \end{equation} \end{itemize} \begin{lemma} \label{lem4.2} When the condition (H') hold, system \eqref{4.1} has no nontrivial $\tau$-periodic solution. \end{lemma} \begin{proof} Denote $u=(u_1,u_2,u_3)^T$. Then system (\ref{4.1}) with $\tau=0$ becomes \begin{equation} \label{4.8} f(u_1,u_2,u_3)=\begin{pmatrix} \dot{u}_1 \\ \dot{u}_2 \\ \dot{u}_3 \end{pmatrix} =\begin{pmatrix} (a_1-c_1)u_3(t)-ru_1(t)-\alpha u_2(t)\\ u_2(t)[-(b_1+r+\alpha)-a_2u_2(t)+a_2u_1(t)] \\ c_1u_3(t)-c_2{u_3}^2(t) \end{pmatrix}. \end{equation} We have \[ \frac{\partial f}{\partial u}=\begin{pmatrix} -r & -\alpha & a_1-c_1 \\ a_2u_2(t) & -(b_1+r+\alpha)-2a_2u_2(t) & 0 \\ 0 & 0 & c_1-2c_2u_3(t) \end{pmatrix} \] and \[ \frac{\partial f^{[2]}}{\partial u}=\begin{pmatrix} a_{11} & 0 & c_1-a_1 \\ 0 & -r+c_1-2c_2u_3(t) & -\alpha \\ 0 & a_2u_2(t) & a_{23} \end{pmatrix} \] where $a_{11}=-r-(b_1+r+\alpha)-2a_2u_2(t),a_{23} = -(b_1+r+\alpha)-2a_2u_2(t)+c_1-2c_2u_3(t)$. The second compound system \[ \begin{pmatrix} \dot{u}_1 \\ \dot{u}_2 \\ \dot{u}_3 \end{pmatrix} =\frac{\partial f^{[2]}}{\partial u} \begin{pmatrix} u_1 \\ u_2 \\ u_3 \end{pmatrix} \] is \begin{equation} \label{4.9} \begin{gathered} \dot{u}_1(t)=[-r-(b_1+r+\alpha)-2a_2u_2(t)]u_1(t)+(c_1-a_1)u_3(t)\\ \dot{u}_2(t)=u_2(t)[-r+c_1-2c_2u_3(t)]-\alpha u_3(t)\\ \dot{u}_3(t)=a_2{u_2}^2(t)-[-(b_1+r+\alpha)-2a_2u_2(t)+c_1-2c_2u_3(t)]u_3(t) \end{gathered} \end{equation} Set \begin{equation} \label{4.10} W(u)=\max\{N_1|u_1|,N_2|u_2|,|u_3|\} \end{equation} where $N_1,N_2$ are constants. Then direct calculation leads to the equation \begin{align*} \frac{d^+}{dt}N_1|u_1(t)| &\leq[-r-(b_1+r+\alpha)]N1|u_1(t)|-2a_2N_1|u_1(t)||u_2(t)|\\ &\quad +N_1||u_1|(t)|+N_1(c_1-a_1)|u_3(t)|\\ \frac{d^+}{dt}N_2|u_2(t)| &\leq |u_2(t)|(-r+c_1)N_2|u_2(t)|-2N_2c_2|u_2(t)||u_3(t)|-N_2\alpha|u_3(t)|\\ \frac{d^+}{dt}|u_3(t) &|\leq \frac{a_2}{N_2}|u_2(t)|N_2|u_2(t)|+(c_1-b_1-r-\alpha)|u_3(t)|\\ &\quad -2a_2|u_2(t)||u_3(t)|-2c_2|u_3(t)|^2 \end{align*} where $\frac{d^+}{dt}$ denotes the right-hand derivative. Therefore, \[ \frac{d^+}{dt}W(u)\leq \mu(t)W(u(t)) \] with \begin{align*} \mu(t)&=\max\big\{[-r-(b_1+r+\alpha)]+N_1(c_1-a_1)+2a_2|u_2|,\\ &\quad (-r+c_1)|u_2(t)|+2N_2c_2|u_2(t)|+N_2\alpha,\\ &\quad (c_1-b_1-r-\alpha)+(a_2/N_2)|u_2(t)|+2c_2|u_3(t)|+2a_2|u_2(t)| \big\} \end{align*} Thus, under the hypothesis (H'), and by the boundedness of solution to (\ref{4.8}), there exists a $\delta>0$ such that $\mu(t)\leq -\delta < 0$, and \[ W(u(t)) \leq W(u(t))e^{-\delta(t-s)},\quad t \geq s > 0. \] This implies that the second compound system (\ref{4.8}) is equi-uniform asymptotically stable, and hence the system (\ref{4.1}) has no non-constant periodic solutions when $\tau=0$ follows from proposition \ref{prop4.2}. \end{proof} \begin{theorem} \label{thm4.2} Assume that \begin{itemize} \item[(i)] the condition \eqref{C}, (H) and (H') hold; \item[(ii)] the condition $(iii)$ of Theorem \ref{thm2.5} hold. \end{itemize} Then, for every $\tau> \tau_k^{(j)}$, system \eqref{4.1} has a $p$-periodic solution. \end{theorem} \begin{proof} We regard $(\tau,p)$ as parameters and apply proposition \ref{prop4.1}. It is sufficient to prove that the projection of $\ell(u^\ast,\tau_k^{(j)},\frac{2\pi}{\omega_0})$ onto $\tau$-space is unbounded. The characteristic matrix of (\ref{4.1}) at an equilibrium $\bar{u}=(\bar{u}^{(1)}, \bar{u}^{(2)}, \bar{u}^{(3)}) \in R^3$ takes the following form \begin{equation}\label{4.11} \Delta(\bar{u},\lambda,\tau)=\lambda Id - DF(\bar{u},\bar{\tau},\bar{p})(e^\lambda Id) \end{equation} By \eqref{2.1}, we can easily show that $(\bar{u}_i,\tau,p)(i=1,2)$ are the stationary solution of (\ref{4.1})(where $\bar{u}_1=(0,0,0), \bar{u}_2=(\frac{c_1(a_1-c_1)}{rc_2},0,\frac{c_1}{c_2})$) and the corresponding characteristic matrices are \begin{gather}\label{4.12} \Delta_{(\bar{u}_1,\tau,p)}(\lambda) =\begin{pmatrix} \lambda+r_1 & \alpha & -a_1+c_1e^{-\lambda\tau}\\ 0 & \lambda+(b_1+r+\alpha) & 0\\ 0 & 0 & \lambda-c_1e^{-\lambda\tau} \end{pmatrix},\\ \label{4.13} \Delta_{(\bar{u}_2,\tau,p)}(\lambda)= \begin{pmatrix} \lambda+r_1 & \alpha & -a_1+c_1e^{-\lambda\tau}\\ 0 & \lambda+(b_1+r+\alpha) & 0\\ 0 & 0 & \lambda-c_1e^{-\lambda\tau}+2c_1. \end{pmatrix} \end{gather} The point $(\bar{u},\bar{\tau},\bar{p})$ is called a center if $F(\bar{u},\bar{\tau},\bar{p})=0$ and $\det(\Delta_{(\bar{u},\bar{\tau},\bar{p})}(\lambda))(\frac{2\pi}{p}i)=0$. A center $(\bar{u},\bar{\tau},\bar{p})$ is said to be isolated if it is the only center in some neighborhood of $(\bar{u},\bar{\tau},\bar{p})$. It follows from (\ref{4.12}) and (\ref{4.13}) that \begin{gather} \label{4.14} \det(\Delta_{(\bar{u}_1,\tau,p)}(\lambda)) =(\lambda+r_1)(\lambda+b_1+r+\alpha)(\lambda-c_1e^{-\lambda\tau})=0 \\ \label{4.15} \det(\Delta_{(\bar{u}_1,\tau,p)}(\lambda))=(\lambda+r_1) (\lambda+b_1+r+\alpha)(\lambda-c_1e^{-\lambda\tau}+2c_1)=0. \end{gather} Obviously, \eqref{4.15} has no purely imaginary roots. Thus, we conclude that (\ref{4.1}) has no center of the form $(\bar{u}_2,\tau,p)$. For $\omega>0$, $i\omega$ is a root of \eqref{4.14} if and only if \[ i\omega-c_1(\cos\omega\tau-i\sin\omega\tau)=0 \] Separating the real and imaginary parts, we have \[ c_1\cos\omega\tau=0, \quad \sin\omega\tau=\omega \] which implies \[ \omega=c_1, \quad \tau_k=\frac{\pi}{2c_1}+\frac{2k\pi}{c_1}\,. \] Thus, when $\tau_k=\frac{\pi}{2c_1}+\frac{2k\pi}{c_1}$, Equation (\ref{4.14}) has a pair of simple imaginary roots $\pm ic_1$. By direct computation, we may obtain that \[ \mathop{\rm Re}\{\frac{d\lambda}{d\tau}|_{\tau=\tau_k}\} =\frac{c_1^2}{1+c_1^2\tau_k^2}>0 \] Therefore, we conclude that $(\bar{u}_1, \tau_k, \frac{2\pi}{c_1})$ is an isolated center stated as above. On the other hand, from the discussion about the local Hopf bifurcation in Section 2, it is easy to verify that $(u^*,\tau_k^{(j)},\frac{2\pi}{\omega_0})$ is also an isolated center, and there exist $\varepsilon>0,\delta>0$ and a smooth curve $\lambda:(\tau_k^{(j)}-\delta,\tau_k^{(j)}+\delta)\to C$ such that $\det(\Delta(\lambda(\tau)))=0$, $|\lambda(\tau)-i\omega_0|<\varepsilon$ for all $\tau\in[\tau_k^{(j)}-\delta,\tau_k^{(j)}+\delta]$, and \[ \lambda(\tau_k^{(j)})=i\omega_0,\ \ \frac{d}{d\tau}\mathop{\rm Re}\lambda(\tau)|_{\tau=\tau_k^{(j)}}\neq 0. \] Let \[ \Omega_\varepsilon=\{(v,p);0 \tau_k^{(j)}$, there exists an integer m such that $\tau/m < \frac{2\pi}{\omega_0} < \tau$. Since system (\ref{4.1}) has no periodic solution when $\tau=0$, it has no $\tau/n$-periodic solution for any integer $n$. Thus, the period $p$ of a periodic solution on the connected component $\ell(u^\ast,\tau_k^{(j)},\frac{2\pi}{\omega_0})$ satisfies $\tau/m < p < \tau$, through $(u^\ast,\tau_k^{(j)},\frac{2\pi}{\omega_0})$ is unbounded. Now we prove that the projection of $\ell(u^\ast,\tau_k^{(j)},\frac{2\pi}{\omega_0})$ onto $\tau-$space is unbounded. Lemma \ref{lem4.1} implies that the projection of $\ell(u^\ast,\tau_k^{(j)},\frac{2\pi}{\omega_0})$ onto the $u$-space is bounded. Aslo, note that Lemma \ref{lem4.2} shows that the system (\ref{4.1}) with $\tau=0$ has no non-constant periodic solution. Therefore, the projection of $\ell(u^\ast,\tau_k^{(j)},\frac{2\pi}{\omega_0})$ onto the $\tau$-space is bounded below. The definition of $\tau_k^{(j)}$ in \eqref{2.12} implies that $\frac{2\pi}{\omega_0} <\tau_k^{(j)},k=1,2,\ldots$ and applying Lemma \ref{lem4.2}, we know that for $\tau > \tau_k^{(j)}$, there exists an integer m such that $\tau/m < \frac{2\pi}{\omega_0} < \tau$. Since system (\ref{4.1}) has no periodic solution when $\tau=0$, it has no $\tau/n$-periodic solution for any integer $n$. Thus, the period p of a periodic solution on the connected component $\ell(u^\ast,\tau_k^{(j)},\frac{2\pi}{\omega_0})$ satisfies $\tau/m < p < \tau$. This shows that in order for $\ell(u^\ast,\tau_k^{(j)},\frac{2\pi}{\omega_0})$ to be unbounded, the projection of $\ell(u^\ast,\tau_k^{(j)},\frac{2\pi}{\omega_0})$ onto $\tau-$space must be unbounded. Consequently, the projection of $\ell(u^\ast,\tau_k^{(j)},\frac{2\pi}{\omega_0})$onto $\tau-$space must be an interval $[\alpha,\infty]$ with $0<\alpha<\tau_k^{(j)}$. This shows for each $\tau > \tau_k^{(j)}$, system (\ref{4.1}) has a p-periodic solution. This completes the proof. \end{proof} \section{An example} We consider the system \begin{equation}\label{5.1} \begin{gathered} \dot{x}(t)=2z(t)-1.2z(t-\tau)-0.8x(t)-0.16x(t)y(t)+0.75y(t),\\ \dot{y}(t)=0.16x(t)-0.75y(t)-0.98y(t),\\ \dot{z}(t)=1.2z(t-\tau)-0.03z^2(t). \end{gathered} \end{equation} which has a positive equilibrium $E_*=(27.9622,2.5861,40)$. It follows from section 2 that $z_1=1.5812$, $\tau_0=1.0441$, $\tau=1.0441+3.9736j$($j=0,1,2,\dots$) and $h(z_1)=-25.8025\neq 0$. By Theorem \ref{thm2.5}, we know that the positive equilibrium is stable when $\tau < \tau_0$. Figure 1 shows that the positive equilibrium is asymptotically stable when $\tau=0.4$. \begin{figure}[ht] \begin{center} \includegraphics[width=0.325\textwidth, height=0.23\textheight]{figure11} \includegraphics[width=0.325\textwidth, height=0.23\textheight]{figure12} \includegraphics[width=0.325\textwidth, height=0.23\textheight]{figure13} \end{center} \caption{} \end{figure} By section 2,we know that a Hopf bifurcation occurs when $\tau=\tau_0$, the positive equilibrium loses its stability and a periodic solution bifurcating from the positive equilibrium exists for $\tau=1.1>\tau_0$. The bifurcation is supercritical and the bifurcating periodic solution is orbitally asymptotically stable,as depicted in Figure 2. \begin{figure}[ht] \begin{center} \includegraphics[width=0.325\textwidth, height=0.23\textheight]{figure21} \includegraphics[width=0.325\textwidth, height=0.23\textheight]{figure22} \includegraphics[width=0.325\textwidth, height=0.23\textheight]{figure23} \end{center} \caption{} \end{figure} \begin{thebibliography}{00} \bibitem{Bence} J. R. Bence, R. M. Nisbet; \emph{Space-limited recruitment in open systems:The importance of time delays}, Ecology, (70), (1989) 1434-1441. \bibitem{Fiedler} M. Fiedler; \emph{Additive compound matrices and inequality for eigenvalues of stochastic matrices}, Czech, Math. J., (99), (1974) 392-402. \bibitem{Hassaard} B. D. Hassaard, N. D. Kazarinoff, Y. H. Wan; \emph{Theory and Application of Hopf Bifurcation, Cambrigdge University Press, Cambridge, 1981.} \bibitem{Hethcote} H. W. Hethcote; \emph{A Thousand and one epidemic models:Frontiers in Mathematical Biology}, Levin, S. A. Lecture Notes in biomathematics 100, Springer-Verlag, Berlin, Heideberg, New York, (1994), 304-305. \bibitem{Lefere} C. Lefere; \emph{Threshold behavior for a chain binomial SIS infectious disease}, J. Math. Biol, (24), (1986) 59-70. \bibitem{Li} M. Y. Li, J. Muldowney; \emph{On Bendixson's criterion}, J. Differential Equation, (106), (1994) 27-39. \bibitem{Liu} W. M. Liu, H. W. Hethcote, S. A. Levin; \emph{Dynamical behavior of epidemiological models with nonlinear incidence rates}, J. Math. Biol, (25), (1987) 359-386. \bibitem{Muldowney} J. S. Muldowney; \emph{Compound matrices and ordinary differential equations}, Rocky Mountain J. Math, (20), (1990) 857-871. \bibitem{Song} Y. L. Song, S. Yuan; \emph{Bifurcation analysis in a predator-prey system with time delay}, Nonlinear Analysis, (7) (2006) 265-284. \bibitem{Wu} J.Wu; \emph{Symmetric functional differential equations and neural networks with memory}, Trans. Amer. Math. Soc, (350), (1998) 4799-4838. \bibitem{Xiao} Y. N. Xiao, L. S. Chen; \emph{On an SIS epidemic model with stage structure}, J. Systems Science and Complexity, (16), (2003) 275-288. \end{thebibliography} \end{document}