\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2009(2009), No. 50, pp. 1--13.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2009 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document}
\title[\hfilneg EJDE-2009/50\hfil Stability and bifurcation]
{Stability and bifurcation in a delayed predator-prey
system with stage structure and functional response}

\author[X. Zhang, R. Xu, Q. Gan\hfil EJDE-2009/50\hfilneg]
{Xiao Zhang, Rui Xu, Qintao Gan}  % not in alphabetical order

\address{Institute of Applied Mathematics,
Mechanical Engineering College, Shijiazhuang 050003, China}

\email[Xiao Zhang]{xzhang\_82@163.com} 
\email[Rui Xu]{rxu88@yahoo.com.cn} 
\email[Qintao Gan]{ganqintao@sina.com}

\thanks{Submitted September 16, 2008. Published April 7, 2009.}
\thanks{Supported by grant 10671209 from the National Natural Science
 Foundation of China}
\subjclass[2000]{34K18, 34K20, 34K60, 92D25}
\keywords{Time delay; stage structure; stability;  Hopf
bifurcation; \hfill\break\indent functional response}

\begin{abstract}
 In this paper, a delayed predator-prey system with
 stage structure and Holling type-II functional response is
 investigated. The local stability of a positive equilibrium and the
 existence of Hopf bifurcations are established. By using the normal
 form theory and center manifold reduction, explicit formulae
 determining the stability, direction of the bifurcating periodic
 solutions are derived. Finally, numerical simulations are carried
 out to illustrate the theoretical results.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{example}[theorem]{Example}

\section{Introduction}

In the natural world, almost all animals have the stage structure of
immature and mature. And the vital rates (rates of survival,
development and reproduction and so on) are often quite different in
these two stage. Hence, it is of ecological importance to
investigate the effects of such a subdivision on the interaction of
species. Chen \cite{chenl} introduced the following stage-structured
single-species population model without time delay
\begin{equation}\label{z1}
\begin{gathered}
\dot{N}_{i}(t)= B(t)-D_{i}(t)-W(t),\\
\dot{N}_{m}(t)= \alpha W(t)-D_{m}(t),
\end{gathered}
\end{equation}
where $N_{i}(t)$ and $N_{m}(t)$ denote the immature and mature
population densities at time t; $B(t)$ is the birth rate of the
immature population at time t; $D_{i}(t)$ and $D_{m}(t)$ are the
death rates of the immature and mature at time t; $W(t)$ represents
the transformation rate of the immature into the mature; $\alpha$ is
the probability of the successful transformation of the immature
into the mature. If the birth rate of model (\ref{z1}) obeys the
Malthus rule, i.e., $B(t)=aN_{m}$; and the death rates of the
immature and mature populations are logistic; i.e.,
\[
D_{i}(t)=r_{i}N_{i}(t)+b_{i}N_{i}^2(t), \quad
D_{m}(t)=r_{m}N_{m}(t)+b_{m}N_{m}^2(t),
\]
and the transformation rate of the immature into mature is
proportional the immature population, i.e., $W(t)=bN_{i}(t)$. Then
we rewrite model (\ref{z1}) as
\begin{equation}\label{z2}
\begin{aligned}
\dot{N}_{i}(t)= aN_{m}-r_{i}N_{i}(t)-b_{i}N_{i}^2(t)-bN_{i}(t),\\
\dot{N}_{m}(t)= bN_{i}(t)-r_{m}N_{m}(t)-b_{m}N_{m}^2(t).
\end{aligned}\end{equation}
Based on the idea above, many authors studied different kinds of
stage-structured models and a significant body of work has been
carried out (see, for example,
\cite{cuij,cuis,gm,gsj,mzp,rx,wang,zh}).

Xu and Ma \cite{xr2} considered  the following model, where prey
population have stage structure and immature individuals are
consumed by their predator
\begin{equation}\label{z3}
\begin{gathered}
\dot{x}_{1}(t)=ax_{2}-r_{1}x_{1}(t)-a_{11}x_{1}^2(t)-bx_{1}(t)-a_{1}x_{1}(t)y(t),\\
\dot{x}_{2}(t)=bx_{1}(t)-r_{2}x_{2}(t),\\
\dot{y}(t)=a_{2}x_{1}(t-\tau)y(t-\tau)-ry(t)-\beta y^2(t),
\end{gathered}
\end{equation}
where $a>0$ is the birth rate of immature population; $r_1>0$ is the
death rate of the immature population; $a_{11}>0$ is the
intra-specific competition rate of the immature population; $b>0$ is
the transformation rate from the immature individuals to mature
individuals; $a_1>0$ is the capturing rate of the predator
population; $r_2>0$ is the death rate of the mature population;
$a_2/a_1$ is the conversion rate of nutrients into the reproduction
of the predator; $r>0$ is the death rate of the predator; $\beta>0$
is the intra-specific competition rate of the predator, $\tau\geq0$
is a constant delay due to the gestation of the predator. In
\cite{xr2}, the global stability of the positive equilibrium and the
two boundary equilibria of the model (\ref{z3}) is discussed.

In system \eqref{z3}, the capture ability of the predator population
is proportional the number of prey population. That is to say if the
number of prey population is large, then predator will capture more
prey in unit time. Obviously, this is reasonless. The predator
always have a handling time. In 1965, Holling proposed three types
of functional response, proved that the functional response plays an
important role in predator-prey systems. Based on the work developed
in \cite{xr2}, in the present paper,  we are concerned with the
effect of functional response on the dynamics of a predator-prey
model. To this end, we consider the following delay differential
equations
 \begin{equation}\label{z4}
\begin{gathered}
\dot{x}_{1}(t)=ax_{2}-r_{1}x_{1}(t)-a_{11}x_{1}^2(t)-bx_{1}(t)-\frac{a_{1}x_{1}(t)y(t)}{1+mx_{1}(t)},\\
\dot{x}_{2}(t)=bx_{1}(t)-r_{2}x_{2}(t),\\
\dot{y}(t)=\frac{a_{2}x_{1}(t-\tau)y(t-\tau)}{1+mx_{1}(t-\tau)}-ry(t),
\end{gathered}
\end{equation}
where the positive $a,r_1,a_{11},b,r_2,r,a_1,a_2$ are  meanings to
same as system \eqref{z3}; the parameter $m$ stands for half
capturing saturation; $\frac{x_1(t)}{1+mx_1(t)}$ is Holling type-II
functional response, which reflects the capture ability of predator.
In the system \eqref{z4} we ignore the intra-specific competition
rate of the predator population.

 The initial conditions for system \eqref{z4} take the form
\begin{equation}\label{z5}
\begin{gathered}
x_1(\theta)=\phi_1(\theta), \quad x_2(\theta)=\phi_2(\theta),\quad
y(\theta)=\psi(\theta), \\
\phi_1(\theta)\geq 0, \quad \phi_2(\theta)\geq 0,\quad
\psi(\theta)\geq 0,\quad \theta
\in [-\tau, 0],\\
 \phi_1(0)>0, \quad \phi_2(0)>0,\quad \psi(0)>0,
\end{gathered}
\end{equation}
where $(\phi_1(\theta), \phi_2(\theta), \psi(\theta))\in C([-\tau,
0], \mathbb{R}^3_{+0})$, the Banach space of continuous functions mapping the
interval $[-\tau,0]$ into $\mathbb{R}^3_{+0}$, where
$ \mathbb{R}^3_{+0}=\{(x_1,x_2,x_3): x_{i}\geq 0,i=1,2,3\}$.

This paper is organized as follows. In Section 2, we discuss the
local stability of a positive equilibrium and the existence of Hopf
bifurcations of system \eqref{z4}. In Section 3, we study the
stability of the bifurcating periodic solutions and the direction of
the Hopf bifurcation in system \eqref{z4} by using the normal form
theory and center manifold reduction. Finally, some numerical
examples are given to illustrate the results above.


\section{Stability and Hopf bifurcations}

In this section, we discuss the local stability of a positive
equilibrium and the existence of Hopf bifurcations of system
\eqref{z4}.

It is easy to show that if the following holds:
\begin{itemize}
\item[(H1)] $a_2-rm>0$,
 $\frac{ab}{r_2}-\frac{a_{11}r}{a_2-rm}-r_1-b>0$,
\end{itemize}
system \eqref{z4} has a unique positive
equilibrium $E^*=(x_{1}^*,x_{2}^*,y^*)$, where
$$
x_{1}^*=\frac{r}{a_2-rm},\quad
x_{2}^*=\frac{br}{r_2(a_2-rm)}, \quad
y^*=\frac{a_2}{a_1(a_2-rm)}\big(\frac{ab}{r_2}
-\frac{a_{11}r}{a_2-rm}-r_1-b\big).
$$
Linearizing system \eqref{z4} at $E^*$, we derive
\begin{equation}\label{a1}
\begin{gathered}
\dot{x}_1(t)=p_1x_1(t)+p_2x_2(t)+p_3y(t),\\
\dot{x}_2(t)=p_4x_1(t)+p_5x_2(t),\\
\dot  y(t)=p_6x_1(t-\tau)+p_7y(t-\tau)+p_{8}y(t).\\
\end{gathered}
\end{equation}
where
\begin{gather*}
p_1=-r_1-b-2a_{11}x_1^{*}-\frac{a_1y^{*}}{(1+mx_1^{*})^2},\quad
p_2=a,\quad p_3=-\frac{a_1r}{a_2},\\
p_4=b,\quad p_5=-r_2,\quad p_6=\frac{a_2y^*}{(1+mx_1^{*})^2},\quad
p_7=r,\quad p_8=-r.
\end{gather*}
The characteristic equation of system \eqref{a1} takes the form
 \begin{equation}\label{a2}
\lambda^3+A\lambda^2+B\lambda+C+(D\lambda^2+E\lambda+F)e^{-\lambda\tau}=0,
\end{equation}
where
\begin{gather*}
A=-(p_1+p_5+p_8),\quad
B=p_1p_5+p_1p_8+p_5p_8-p_2p_4,\quad C=p_8(p_2p_4-p_1p_5),\\
D=-p_7,\quad E=p_1p_7+p_5p_7-p_3p_6,\quad
F=p_2p_4p_7+p_3p_5p_6-p_1p_5p_7.
\end{gather*}
If $  i\omega(\omega>0)$ is a
root of  \eqref{a2}, then
\begin{equation}\label{a3}
- i\omega^3-A\omega^2+B i\omega+C+(-D\omega^2+E i\omega+F)
e^{-i\omega\tau}=0.
\end{equation}
Separating the real and imaginary parts, we obtain
\begin{equation}\label{a4}
\begin{gathered}
C-\omega^2=D\omega^2\cos\omega\tau-E\omega\sin\omega\tau-F\cos\omega\tau,\\
-\omega^3+B\omega=-D\omega^2\sin\omega\tau-E\omega\cos\omega\tau+F\sin\omega\tau.
\end{gathered}
\end{equation}
It follows from \eqref{a4} that
\begin{equation}\label{a5}
\omega^6+Q_1\omega^4+Q_2\omega^2+Q_3=0,
\end{equation}
where
$$
Q_1=A^2-2B-D^2,\quad Q_2=B^2+2DF-2AC-E^2,\quad Q_3=C^2-F^2.
$$
 It is easy to show that
 \begin{gather*}
Q_1=A^2-D^2-2B=p_1^2+p_5^2+2p_2p_4>0,\\
Q_2=(p_1p_5-p_2p_4)^2+p_3p_6(2p_1p_7-p_3p_6).
 \end{gather*}
 Noting that if $Q_3>0$, we have $2(p_1p_5-p_2p_4)p_7>p_3p_5p_6$, it
 is easy to see $Q_2>0$. Hence, the stability of positive equilibrium $E^*$ of system \eqref{z4}
is unchanged for all $\tau>0$. When $\tau=0$, by Routh-Hurwitz
Criterion we obtain that all roots of  \eqref{a2} have negative
real parts only if
\begin{itemize}
\item[(H2)] $ p_1p_3p_6>(p_1+p_5)(p_1p_5-p_2p_4)$.
\end{itemize}
Therefore if (H2) holds, the positive equilibrium $E^*$ of system
\eqref{z4} is local asymptotically stable for all $\tau>0$; or else
the positive equilibrium $E^*$ of system \eqref{z4} is unstable for
all $\tau>0$.

If $Q_3<0$, we know that \eqref{a5} has only one unique positive
root $\omega_0$.
 Define
\begin{equation}\label{a6}
\tau_j=\frac{1}{\omega_0}\Big(\arccos\frac{(E-DA)\omega_0^4+(CD+AF-BE)\omega_0^2-FC}
{D^2\omega_0^4+(E^2-2FD)\omega_0^2+F^2}+2j\pi\Big),
\end{equation}
for $j=0,1,\dots$. Then $(\tau_j,\omega_0)$ solves  \eqref{a3}.
This means that when $\tau=\tau_j$, Equation \eqref{a2} has a pair of
purely imaginary roots $\pm  i\omega_0$.

Denote $ \lambda=\lambda(\tau)$,  from  \eqref{a2} we have
$$
\Big(\frac{\mathrm{d}\lambda(\tau)}{\mathrm{d}\tau}\Big)^{-1}=-\frac{3\lambda^2+2A\lambda+B}
{\lambda(\lambda^3+A\lambda^2+B\lambda+C)}+\frac{2D\lambda+E}
{\lambda(D\lambda^2+E\lambda+F)}-\frac{\tau}{\lambda},
$$
which leads to   \begin{equation}\label{a7}
\begin{aligned}
&\mathop{\rm sign}\big\{\mathop{\rm Re}\big(\frac{\mathrm{d}\lambda}
{\mathrm{d}\tau}\big)_{\tau=\tau_j}\big\}\\
&=\mathop{\rm sign}\big\{\mathop{\rm Re}
\big(\frac{\mathrm{d}\lambda}{\mathrm{d}\tau}\big)_{\tau=\tau_j}^{-1}\big\}\\
&=\mathop{\rm sign}\big\{3\omega_0^4+(2A^2-4B-2D^2)\omega_0^2+B^2
+2DF-2AC-E^2\big\}\\
&=\mathop{\rm sign}\big\{3\omega_0^4+2Q_1\omega_0^2+Q_2\big\}.
\end{aligned}
\end{equation}
Because $Q_1>0$, so if the following holds
\begin{itemize}
\item[(H3)] $Q_2>0$,
\end{itemize}
then $3\omega_0^4+2Q_1\omega_0^2+Q_2>0$,
$\mathop{\rm sign}\big\{\mathop{\rm Re}
\big(\frac{\mathrm{d}\lambda}{\mathrm{d}\tau}\big)_{\tau=\tau_j}\big\}>0$.
And system \eqref{a1} undergoes a Hopf bifurcation.

Applying  \cite[Theorem 11.1]{hale}, we obtain the following results.

\begin{theorem} \label{thm2.1}
Let {\rm (H1)} hold, and
\begin{itemize}
\item[(i)] If $Q_3>0$ and {\rm (H2)} holds, then the positive
equilibrium $E^*$ of system
\eqref{z4} is local asymptotically stable for all $\tau>0$; or else
the positive equilibrium $E^*$ of system \eqref{z4} is unstable for
all $\tau>0$.

\item[(ii)] If $Q_3<0$ and {\rm (H2), (H3)} hold, the positive
equilibrium $E^*$ of system \eqref{z4} is asymptotically stable
for $\tau\in[0,\tau_0)$, and
$E^*$ is unstable for $\tau>\tau_0$; and the system \eqref{z4}
undergoes a Hopf Bifurcation at the positive equilibrium $E^*$ when
$\tau=\tau_j (j=0,1,\dots)$.
\end{itemize}
\end{theorem}

Xu and Ma \cite{xr2} obtained the following results:
\begin{itemize}
\item[(i)] If
\begin{align*}
&2a_2a_{11}{\beta[ab-r_2(r_1+b)]+a_1r_2r}\\
&+(a_{11}\beta-a_1a_2)a_2[ab-r_2(r_1+b)]-a_{11}r_2r>0,
\end{align*}
the positive equilibrium $E^*$ of system \eqref{z3} is local
asymptotically stable for all $\tau>0$.

\item[(ii)] If
\begin{align*}
&2a_2a_{11}{\beta[ab-r_2(r_1+b)]+a_1r_2r}\\
&+(a_{11}\beta-a_1a_2){a_2[ab-r_2(r_1+b)]-a_{11}r_2r}<0,
\end{align*}
the positive equilibrium $E^*$ of system \eqref{z3} is
asymptotically stable for $\tau\in[0,\tau_0)$, and $E^*$ is unstable
for $\tau>\tau_0$.
\end{itemize}

In system \eqref{z3}, we let $\beta=0$, then we know that if
$3a_{11}rr_2>a_2(ab-r_2(r_1+b))$, the system \eqref{z3} is local
asymptotically stable for all $\tau>0$; if
$3a_{11}rr_2<a_2(ab-r_2(r_1+b))$, the system \eqref{z3} is
asymptotically stable for $\tau\in[0,\tau_0)$, and unstable for
$\tau>\tau_0$.

If we let $m=0$ in system \eqref{z4}, then $Q_3>0$ means that
$3a_{11}rr_2>a_2(ab-r_2(r_1+b))$, same to the results of system
\eqref{z3}. So we know that system \eqref{z3}(when $\beta=0$) is the
special situation of system \eqref{z4}.


\section{Direction of Hopf bifurcation}

In this section, we derive explicit formulae to determine the
properties of the Hopf bifurcation at critical values $\tau_j$ by
using the normal form theory and center manifold reduction  (see,
for example, Hassard et al. \cite{hasd}).

Without loss of generality, denote the critical values $\tau_j$ by
$\tilde{\tau}$, and set $\tau=\tilde{\tau}+\mu$. Then $\mu=0$ is a
Hopf bifurcation value of system \eqref{z4}. Thus, we can work in
the phase space $C=C([-\tilde{\tau},0],\mathbb{R}^3)$.
Let
$$
u_1(t)=x_1(t)-x_1^*,\quad u_2(t)=x_2(t)-x_2^*,\quad
u_3(t)=y(t)-y^*.
$$
Then system \eqref{z4} is transformed into
\begin{equation}\label{c1}
\begin{gathered}
\dot u_1(t)=p_1u_1(t)+p_2u_2(t)+p_3u_3(t)+
\sum_{i+j+l\geq2}\frac{1}{i!j!l!}f_{ijl}^{(1)}u_1^i(t)u_2^j(t)u_3^l(t),\\
\dot u_2(t)=p_4u_1(t)+p_5u_2(t),\\
\begin{aligned}
\dot u_3(t)=&p_6u_1(t-\tau)+p_7u_3(t-\tau)+p_8u_3(t)\\
&+\sum_{i+j+l\geq2}\frac{1}{i!j!l!}f_{ijl}^{(3)}u_1^i(t-\tau)u_3^j(t-\tau)u_3^l(t),
\end{aligned}
\end{gathered}
\end{equation}
where
\[
f_{ijl}^{(1)}=\frac{\partial^{i+j+l}f^{(1)}}{\partial x_1^i\partial
x_2^j\partial y^l}\Big|_{(x_1^*,x_2^*,y^*)},\quad
f_{ijl}^{(3)}=\frac{\partial^{i+j+l}f^{(3)}}{\partial
x_1^i(t-\tau)\partial y^j(t-\tau)\partial
y^l}\Big|_{(x_1^*,y^*,y^*)},
\]
for $i,j,l\geq 0$, and
\begin{gather*}
f^{(1)}=ax_{2}-r_{1}x_{1}(t)-a_{11}x_{1}^2(t)-bx_{1}(t)
-\frac{a_{1}x_{1}(t)y(t)}{1+mx_{1}(t)},
\\
f^{(3)}=\frac{a_{2}x_{1}(t-\tau)y(t-\tau)}{1+mx_{1}(t-\tau)}-ry(t).
\end{gather*}
We rewrite \eqref{c1} as
\begin{equation}\label{c2}
\dot u(t) = L_\mu(u_t)+f(u_t,\mu),
\end{equation}
where $u(t)=(u_1(t),u_2(t),u_3(t))^T\in\mathbb{R}^3$,
 $u_t(\theta)\in C$ is defined by $u_t(\theta)=u(t+\theta)$,
and $L_\mu:C\to \mathbb{R}^3, f:R\times C\in \mathbb{R}^3$ are given by
\begin{equation}\label{c3}
L_\mu\phi=\begin{pmatrix}
p_1&p_2&p_3\\p_4&p_5&0\\0&0&p_8
\end{pmatrix} \phi(0) +
\begin{pmatrix} 0&0&0\\0&0&0\\p_6&0&p_7
\end{pmatrix} \phi(-\tau),
\end{equation}
and
\begin{equation}\label{c4}
f(\phi,\mu)=\begin{pmatrix}
\sum_{i+j+l\geq2}\frac{1}{i!j!l!}f_{ijl}^{(1)}\phi_1^i(0)\phi_2^j(0)
\phi_3^l(0)\\
0\\
\sum_{i+j+l\geq2}\frac{1}{i!j!l!}f_{ijl}^{(3)}
\phi_1^i(-\tau)\phi_3^j(-\tau)\phi_3^l(0)
\end{pmatrix}
\end{equation}
respectively. By Riesz representation theorem, there exists a
function $\eta(\theta,\mu)$ of bounded variation for
$\theta\in[-\tau,0]$ such that
\begin{equation}\label{c5}
L_\mu\phi=\int_{-\tilde{\tau}}^0\,\mathrm{d}\eta(\theta,0)
\phi(\theta)\quad\text{for }\phi\in C.
\end{equation}
In fact, we can choose
\begin{equation}\label{c6}
\eta(\theta,\mu)= \begin{pmatrix}
p_1&p_2&p_3\\p_4&p_5&0\\0&0&p_8
\end{pmatrix} \delta(\theta)- \begin{pmatrix}
0&0&0\\0&0&0\\p_6&0&p_7\end{pmatrix} \delta(\theta+\tau),
\end{equation}
where $\delta$ is the Dirac delta function. For
$\phi\in C^1([-\tau,0],\mathbb{R}^3)$, define
$$
A(\mu)\phi= \begin{cases}
\frac{\mathrm{d}\phi(\theta)}{\mathrm{d}\theta},&
\theta\in[-\tau,0),\\
\int_{-\tilde{\tau}}^0\,\mathrm{d}\eta(\mu,s)\phi(s), &\theta=0.
\end{cases}
$$
and
$$
R(\mu)\phi= \begin{cases}
0, &\theta\in[-\tau,0),\\
f(\phi), &\theta=0.
\end{cases}
$$
Then system \eqref{c2} is equivalent to
\begin{equation}\label{c7}
\dot u_t = A(\mu)u_t+R(\mu)u_t .
\end{equation}
For $\psi\in C^1([0,\tau],(\mathbb{R}^2)^*)$, define
$$
A^*\psi(s)= \begin{cases}
-\frac{\mathrm{d}\psi(s)}{\mathrm{d}s},& s\in(0,\tau],\\
\int_{-\tilde{\tau}}^0\,\mathrm{d}\eta^T(t,0)\psi(-t),& s=0,
\end{cases}
$$
and a bilinear inner product
\begin{equation}\label{c8}
\langle\psi(s),\phi(\theta)\rangle=\bar{\psi}(0)\phi(0)-
\int_{-\tilde{\tau}}^0\int_{\xi=0}^{\theta}\bar{\psi}(\xi-\theta)\,
\mathrm{d}\eta(\theta)\phi(\xi)d\xi,
\end{equation}
where $\eta(\theta)=\eta(\theta,0)$. Then $A(0)$ and $A^*$ are
adjoint operators. By discussions in Section 2 and foregoing
assumption, we know that $\pm i\omega_0$ are eigenvalues of
$A(0)$. Thus, they are also eigenvalues of $A^*$. We first need to
compute the eigenvector of $A(0)$ and $A^*$ corresponding to
$ i\omega_0$ and $- i\omega_0$, respectively.

Suppose that $q(\theta)=(1,\alpha,\beta)^Te^{ i\omega_0\theta}$
is the eigenvector of $A(0)$ corresponding to $ i\omega_0$. Then
$A(0)q(\theta)= i\omega_0q(\theta)$. It follows from the
definition of $A(0)$, \eqref{c5} and \eqref{c6} that
$$
\begin{pmatrix}
p_1- i\omega_0&p_2&p_3\\
p_4&p_5-i\omega_0&0\\
p_6e^{- i\omega_0\tau}&0&p_8+p_7e^{- i\omega_0\tau}- i\omega_0
\end{pmatrix} q(0) =
\begin{pmatrix} 0\\0\\0
\end{pmatrix}.
$$
We therefore derive that
$$
q(0)=(1,\alpha,\beta)^T=
\Big(1,-\frac{p_4}{p_5- i\omega_0},
-\frac{p_6e^{- i\omega_0\tau}}{p_8+p_7e^{- i\omega_0\tau}
- i\omega_0}\Big)^T.
$$
On the other hand, suppose that $q^*(s)=D(1,\alpha^*,\beta^*)
e^{ i\omega_0s}$ is the eigenvector of $A^*$ corresponding to
$- i\omega_0$. From the definition of $A^*$, \eqref{c5} and \eqref{c6}
we have
$$
\begin{pmatrix}
p_1+ i\omega_0&p_2&p_3\\
p_4&p_5+ i\omega_0&0\\
p_6e^{ i\omega_0\tau}&0&p_8+p_7e^{ i\omega_0\tau}+ i\omega_0
\end{pmatrix}^T(q^*(0))^T =
\begin{pmatrix} 0\\0\\0
\end{pmatrix},
$$
then
$$
q^*(0)=D(1,\alpha^*,\beta^*)=D\Big(1,-\frac{p_2}{p_5+ i\omega_0},-\frac{p_3}{p_8+p_7e^{ i\omega_0\tau}+ i\omega_0}\Big).
$$
To assure that $\langle q^*(s),q(\theta)\rangle=1$, we need to
determine the value of $D$. From \eqref{c8}, we can choose
$$
D=\frac{1}{1+\bar{\alpha}\alpha^*+\bar{\beta}\beta^*+(p_6\beta^*+p_7\bar{\beta}\beta^*)\tau
e^{ i\omega_0\tau}}.
$$

In the remainder of this section, we use the same notations as in
Hassard et al. \cite{hasd}. We first compute the coordinates to
describe the center manifold $\mathbf{C}_0$ at $\mu=0$. Let $u_t$ be
the solution of  \eqref{c2} with $\mu=0$. Define
\begin{equation}\label{c9}
z(t)=\langle q^*,u_t\rangle,\quad
W(t,\theta)=u_t(\theta)-2\mathop{\rm Re}\{{z}(t){q}(\theta)\}.
\end{equation}
On the center manifold $\mathbf{C}_0$ we have
$$
W(t,\theta)=W(z(t),\bar{z}(t),\theta)=W_{20}(\theta)\frac{z^2}{2}+
W_{11}(\theta)z\bar{z}+W_{02}(\theta)\frac{\bar{z}^2}{2}+\dots,
$$
$z$ and $\bar{z}$ are local coordinates for center manifold
$\mathbf{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
the solution $u_t\in \mathbf{C}_0$ of \eqref{c2}, since $\mu=0$, we
have \begin{equation}\label{a}
\begin{aligned}
\dot z
=&\langle q^*, \dot{u_t} \rangle =\langle q^*, A(\mu)u_t\rangle + \langle q^*, R(\mu)u_t \rangle\\
=& i\omega_0z+\langle\bar{q}^*(\theta),f(0,W(z,\bar{z},
\theta)+2\mathop{\rm Re}\{{zq}(\theta)\})\rangle\\
=& i\omega_0z+\bar{q}^*(\theta)f(0,W(z,\bar{z},
\theta)+2\mathop{\rm Re}\{{zq}(\theta)\})\\
=& i\omega_0z+\bar{q}^*(0)f(0,W(z,\bar{z},0)+2\mathop{\rm Re}\{{zq}(0)\})\\
:=& i\omega_0z+\bar{q}^*(0)f_0(z,\bar{z}).
\end{aligned}
\end{equation}
We rewrite \eqref{a} as
$\dot z= i\omega_0z+g(z,\bar{z})$,
with \begin{equation}\label{c10}
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}
Noting that
$u_t(\theta)=(u_{1t}(\theta),u_{2t}(\theta),u_{3t}(\theta))
=W(t,\theta)+zq(\theta)+\bar{z}\bar{q}(\theta)$
and $q(\theta)=(1,\alpha,\beta)^Te^{ i\omega_0\theta}$, we have
\begin{gather*}
u_{1t}(0)=z+\bar{z}+W_{20}^{(1)}(0)\frac{z^2}{2}+W_{11}^{(1)}(0)z\bar{z}
+W_{02}^{(1)}(0)\frac{\bar{z}^2}{2}+\dots,\\
u_{2t}(0)=\alpha
z+\bar{\alpha}\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}+\dots,\\
u_{3t}(0)=\beta
z+\bar{\beta}\bar{z}+W_{20}^{(3)}(0)\frac{z^2}{2}+W_{11}^{(3)}(0)z\bar{z}
+W_{02}^{(3)}(0)\frac{\bar{z}^2}{2}+\dots,\\
u_{1t}(-\tilde{\tau})=e^{- i\omega_0\tilde{\tau}}z+e^{ i\omega_0\tilde{\tau}}\bar{z}+
W_{20}^{(1)}(-\tilde{\tau})\frac{z^2}{2}+W_{11}^{(1)}(-\tilde{\tau})z\bar{z}+
W_{02}^{(1)}(-\tilde{\tau})\frac{\bar{z}^2}{2}+\dots,\\
u_{3t}(-\tilde{\tau})=\beta e^{- i\omega_0\tilde{\tau}}z+\bar{\beta}e^{ i\omega_0\tilde{\tau}}\bar{z}+
W_{20}^{(3)}(-\tilde{\tau})\frac{z^2}{2}+W_{11}^{(3)}(-\tilde{\tau})z\bar{z}+
W_{02}^{(3)}(-\tilde{\tau})\frac{\bar{z}^2}{2}+\dots.
\end{gather*}
Thus, it follows from \eqref{c4} and \eqref{c10} that
\begin{align*}
g(z,\bar{z})&=\bar{q}^*(0)f_0(z,\bar{z})\\
&=\bar{D}(1,\bar{\alpha}^*,\bar{\beta}^*) \begin{pmatrix}
\sum_{i+j+l\geq2}\frac{1}{i!j!l!}f_{ijl}^{(1)}u_{1t}^i(0)
u_{2t}^j(0)u_{3t}^l(0)\\
0\\
\sum_{i+j+l\geq2}\frac{1}{i!j!l!}f_{ijl}^{(3)}u_{1t}^i(-\tilde{\tau})
u_{3t}^j(-\tilde{\tau})u_{3t}^l(0)
\end{pmatrix}
\\
&=\bar{D}\Big\{z^2\Big(\frac{1}{2}f_{200}^{(1)}+f_{101}^{(1)}\beta\Big)+\bar{\beta}^*z^2
\Big(\frac{1}{2}f_{200}^{(3)}e^{-2 i\omega_0\tilde{\tau}}
+f_{110}^{(3)}\beta e^{-2 i\omega_0\tilde{\tau}}\Big)\\
&\quad +z\bar{z}\Big(f_{200}^{(1)}+2f_{101}^{(1)}\mathop{\rm Re}\{\beta\}\Big)
+\bar{\beta}^*z\bar{z}\Big(f_{200}^{(3)}+2f_{110}^{(3)}\mathop{\rm Re}\{\beta\}\Big)\\
&\quad +\bar{z}^2\Big(\frac{1}{2}f_{200}^{(1)}+f_{101}^{(1)}\bar{\beta}\Big)
+\bar{\beta}^*\bar{z}^2\Big(\frac{1}{2}f_{200}^{(3)}e^{2 i\omega_0\tilde{\tau}} +f_{110}^{(3)}\bar{\beta}e^{2 i\omega_0\tilde{\tau}}\Big)\\
&\quad +z^2\bar{z}\Big(\frac{1}{2}f_{300}^{(1)}+\frac{1}{2}f_{201}^{(1)}(\bar{\beta}+2\beta)
+\frac{1}{2}f_{200}^{(1)}(W_{20}^{(1)}(0)+2W_{11}^{(1)}(0))\\
&\quad +f_{101}^{(1)}\big(W_{11}^{(3)}(0)+\frac{1}{2}W_{20}^{(1)}(0)\bar{\alpha}
+\frac{1}{2}W_{20}^{(3)}(0)+W_{11}^{(1)}(0)\beta\big)\Big)\\
&\quad +\bar{\beta}^*z^2\bar{z}\Big(\frac{1}{2}f_{300}^{(3)}e^{- i\omega_0\tilde{\tau}} +\frac{1}{2}f_{210}^{(3)}\big(2\beta e^{- i\omega_0\tilde{\tau}}+\bar{\beta}e^{- i\omega_0\tilde{\tau}}\big)\\
&\quad +f_{110}^{(3)}\big(\frac{1}{2}W_{20}^{(1)}(-\tilde{\tau})\bar{\beta}e^{ i\omega_0\tilde{\tau}}+ \frac{1}{2}W_{20}^{(3)}(-\tilde{\tau})e^{ i\omega_0\tilde{\tau}}\\
&\quad +W_{11}^{(1)}(-\tilde{\tau})\beta e^{- i\omega_0\tilde{\tau}}
+W_{11}^{(3)}(-\tilde{\tau})e^{- i\omega_0\tilde{\tau}}\big)\Big)+\dots \Big\}.
\end{align*}
Comparing the coefficients in \eqref{c10}, we get
\begin{gather*}
g_{20}=\bar{D}\Big(f_{200}^{(1)}+2f_{101}^{(1)}\beta
+\bar{\beta}^*\Big(f_{200}^{(3)}e^{-2 i\omega_0\tilde{\tau}}
+2f_{110}^{(3)}\beta e^{-2 i\omega_0\tilde{\tau}}\Big)\Big),\\
g_{11}=\bar{D}\Big(f_{200}^{(1)}+2f_{101}^{(1)}\mathop{\rm Re}\{\beta\}
+\bar{\beta}^*\Big(f_{200}^{(3)}+2f_{110}^{(3)}\mathop{\rm Re}\{\beta\}\Big)\Big),\\
g_{02}=\bar{D}\Big(f_{200}^{(1)}+2f_{101}^{(1)}\bar{\beta}
+\bar{\beta}^*\Big(f_{200}^{(3)}e^{2 i\omega_0\tilde{\tau}}+2f_{110}^{(3)}\bar{\beta}
e^{2 i\omega_0\tilde{\tau}}\Big)\Big),
\end{gather*}
\begin{equation}\label{c11}
\begin{aligned}
g_{21}&=\bar{D}\Big(f_{300}^{(1)}+f_{201}^{(1)}\big(\bar{\beta}
+2\beta\big)+f_{200}^{(1)}\big(W_{20}^{(1)}(0)+2W_{11}^{(1)}(0)\big)\\
&\quad +f_{101}^{(1)}\big(2W_{11}^{(3)}(0)+W_{20}^{(1)}(0))\bar{\alpha}
+W_{20}^{(3)}(0)+2W_{11}^{(1)}(0)\beta\big)\\
&\quad +\bar{\beta}^*\Big(f_{300}^{(3)}e^{- i\omega_0\tilde{\tau}}
+f_{210}^{(3)}\big(2\beta e^{- i\omega_0\tilde{\tau}}+\bar{\beta}e^{- i\omega_0\tilde{\tau}}\big)\\
&\quad +f_{110}^{(3)}\big(W_{20}^{(1)}(-\tilde{\tau})\bar{\beta}e^{ i\omega_0\tilde{\tau}}+ W_{20}^{(3)}(-\tilde{\tau})e^{ i\omega_0\tilde{\tau}}\\
&\quad +2W_{11}^{(1)}(-\tilde{\tau})\beta e^{- i\omega_0\tilde{\tau}}+2W_{11}^{(3)}(-\tilde{\tau})e^{- i\omega_0\tilde{\tau}}\big)\Big)\Big).
\end{aligned}
\end{equation}
We now compute $W_{20}(\theta)$ and $W_{11}(\theta)$.  It follows
from \eqref{c7} and \eqref{c9} that
\begin{equation}\label{c12}
\begin{aligned}
\dot W&=\dot u_t-\dot zq-\dot{\bar{z}}\bar{q}\\
&=\begin{cases}
AW-2\mathop{\rm Re}\{\bar{{q}}^*(0){f}_0{q}(\theta)\},
& \theta\in(0,\tilde{\tau}],\\
AW-2\mathop{\rm Re}\{\bar{{q}}^*(0){f}_0{q}(0)\}+{f}_0,
& \theta=0,
\end{cases}
&:=AW+H(z,\bar{z},\theta),
\end{aligned}
\end{equation}
where
\begin{equation}\label{c13}
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}+\dots.
\end{equation}
On the other hand, on
$\mathbf{C}_0$ near the origin
\begin{equation}\label{c14}
\dot W=W_z\dot z+W_{\bar{z}}\dot{\bar{z}}.
\end{equation}
We derive from \eqref{c12}-\eqref{c14} that
\begin{equation}\label{c15}
(A-2 i\omega_0)W_{20}(\theta)=-H_{20}(\theta),\quad
AW_{11}(\theta)=-H_{11}(\theta), \; \dots.
\end{equation}
It follows from \eqref{c10} and \eqref{c12} that for
$\theta\in[-\tilde{\tau},0)$,
\begin{equation}\label{c16}
H(z,\bar{z},\theta)=-\bar{q}^*(0)f_0q(\theta)
-q^*(0)\bar{f}_0\bar{q}(\theta)=
-gq(\theta)-\bar{g}\bar{q}(\theta).\end{equation}
Comparing the coefficients
in \eqref{c13} gives that for $\theta\in[-\tilde{\tau},0)$,
\begin{gather} \label{c17}
H_{20}(\theta)=-g_{20}q(\theta)-\bar{g}_{02}\bar{q}(\theta),\\
\label{c18}
H_{11}(\theta)=-g_{11}q(\theta)-\bar{g}_{11}\bar{q}(\theta).
\end{gather}
 We derive from \eqref{c15}, \eqref{c17} and the definition of $A$ that
$$
\dot W_{20}(\theta)=2 i\omega_0W_{20}(\theta)
+g_{20}q(\theta)+\bar{g}_{02}\bar{q}(\theta).
$$
Noting that $q(\theta)=q(0)e^{ i\omega_0\theta}$, it follows
that
\begin{equation}\label{c19}
W_{20}(\theta)=\frac{ ig_{20}}{\omega_0}q(0)e^{ i\omega_0\theta}+
\frac{ i\bar{g}_{02}}{3\omega_0}\bar{q}(0)e^{- i\omega_0\theta}
+E_1e^{2 i\omega_0\theta},
\end{equation}
where $E_1=\big(E_1^{(1)},E_1^{(2)},E_1^{(3)}\big)\in \mathbb{R}^3$ is a
constant vector.
Similarly, from \eqref{c15} and \eqref{c18}, we  obtain
\begin{equation}\label{c20}
W_{11}(\theta)=-\frac{ ig_{11}}{\omega_0}q(0)e^{ i\omega_0\theta}+
\frac{ i\bar{g}_{11}}{\omega_0}\bar{q}(0)e^{- i\omega_0\theta}+E_2,
\end{equation}
where $E_2=\big(E_2^{(1)},E_2^{(2)},E_2^{(3)}\big)\in \mathbb{R}^3$
is also a constant vector.

In what follows, we seek appropriate $E_1$ and $E_2$. From the
definition of $A$ and \eqref{c15}, we obtain
\begin{gather} \label{c21}
\int_{-\tilde{\tau}}^0\,\mathrm{d}\eta(\theta)W_{20}(\theta)
=2 i\omega_0W_{20}(0)-H_{20}(0),\\
\label{c22}
\int_{-\tilde{\tau}}^0\,\mathrm{d}\eta(\theta)W_{11}(\theta)=-H_{11}(0),
\end{gather}
where $\eta(\theta)=\eta(0,\theta)$. From \eqref{c12}, it follows
that
\begin{gather} \label{c23}
H_{20}(0)=-g_{20}q(0)-\bar{g}_{02}\bar{q}(0)+\begin{pmatrix}
f_{200}^{(1)}+2f_{101}^{(1)}\beta\\
0\\
f_{200}^{(3)}e^{-2 i\omega_0\tilde{\tau}} +2f_{110}^{(3)}
\beta e^{-2 i\omega_0\tilde{\tau}}
\end{pmatrix},
\\
\label{c24}
H_{11}(0)=-g_{11}q(0)-\bar{g}_{11}\bar{q}(0)+\begin{pmatrix}
f_{200}^{(1)}+2f_{101}^{(1)}\mathop{\rm Re}\{\beta\}\\0\\
f_{200}^{(3)}+2f_{110}^{(3)}\mathop{\rm Re}\{\beta\}
\end{pmatrix}.
\end{gather}
Substituting \eqref{c19} and \eqref{c23} into \eqref{c21} and
noticing that
\begin{gather*}
\Big( i\omega_0I-\int_{-\tilde{\tau}}^0e^{ i\omega_0\theta}\,
\mathrm{d}\eta(\theta)\Big)q(0)=0, \\
\Big(- i\omega_0I-\int_{-\tilde{\tau}}^0e^{- i\omega_0\theta}\,
\mathrm{d}\eta(\theta)\Big)\bar{q}(0)=0,
\end{gather*}
we obtain
$$
\Big(2 i\omega_0I-\int_{-\tilde{\tau}}^0e^{2 i\omega_0\theta}\,
\mathrm{d}\eta(\theta)\Big)E_1=\begin{pmatrix}
f_{200}^{(1)}+2f_{101}^{(1)}\beta\\
0\\
f_{200}^{(3)}e^{-2 i\omega_0\tilde{\tau}}
+2f_{110}^{(3)}\beta e^{-2 i\omega_0\tilde{\tau}}
\end{pmatrix},
$$
which leads to
\begin{align*}
&\begin{pmatrix}
2 i\omega_0-p_1&-p_2&-p_3\\
-p_4&2 i\omega_0-p_5&0\\
-p_6e^{-2 i\omega_0\tilde{\tau}}&0
&2 i\omega_0-p_8-p_7e^{-2 i\omega_0\tilde{\tau}}
\end{pmatrix} E_1 \\
&=\begin{pmatrix}
f_{200}^{(1)}+2f_{101}^{(1)}\beta\\
0\\
f_{200}^{(3)}e^{-2 i\omega_0\tilde{\tau}}
+2f_{110}^{(3)}\beta e^{-2 i\omega_0\tilde{\tau}}
\end{pmatrix}.
\end{align*}
It follows that
\begin{align*}
&E_1^{(1)}\\
&=\frac{1}{A}\det\begin{pmatrix}
f_{200}^{(1)}+2f_{101}^{(1)}\beta&-p_2&-p_3\\0&2 i\omega_0-p_5&0\\
f_{200}^{(3)}e^{-2 i\omega_0\tilde{\tau}} +2f_{110}^{(3)}\beta
e^{-2 i\omega_0\tilde{\tau}}&0&2 i\omega_0-p_8
-p_7e^{-2 i\omega_0\tilde{\tau}}
\end{pmatrix},
\end{align*}
\begin{align*}
&E_1^{(2)}\\
&=\frac{1}{A}\det\begin{pmatrix} 2 i\omega_0
-p_1&f_{200}^{(1)}+2f_{101}^{(1)}\beta&-p_3\\-p_4&0&0\\
-p_6e^{-2 i\omega_0\tilde{\tau}}&f_{200}^{(3)}
e^{-2 i\omega_0\tilde{\tau}} +2f_{110}^{(3)}\beta
e^{-2 i\omega_0\tilde{\tau}}&2 i\omega_0-p_8-p_7e^{-2 i\omega_0\tilde{\tau}}
\end{pmatrix},
\end{align*}
\[
E_1^{(3)}=\frac{1}{A}\det\begin{pmatrix} 2 i\omega_0-p_1&
-p_2&f_{200}^{(1)}+2f_{101}^{(1)}\beta
\\-p_4&2 i\omega_0-p_5&0\\
-p_6e^{-2 i\omega_0\tilde{\tau}}&0&f_{200}^{(3)}
e^{-2 i\omega_0\tilde{\tau}} +2f_{110}^{(3)}\beta
e^{-2 i\omega_0\tilde{\tau}}
\end{pmatrix}.
\]
where
$$
A=\det\begin{pmatrix} 2 i\omega_0-p_1&-p_2&-p_3\\-p_4&2 i\omega_0-p_5&0\\
-p_6e^{-2 i\omega_0\tilde{\tau}}&0&2 i\omega_0-p_8-p_7e^{-2 i\omega_0\tilde{\tau}}
\end{pmatrix}.
$$
Similarly, we can get
$$
\begin{pmatrix}
-p_1&-p_2&-p_3\\-p_4&-p_5&0\\-p_6&0&0
\end{pmatrix} E_2 =
\begin{pmatrix}
f_{200}^{(1)}+2f_{101}^{(1)}\mathop{\rm Re}\{\beta\}\\0\\
f_{200}^{(3)}+2f_{110}^{(3)}\mathop{\rm Re}\{\beta\}
\end{pmatrix},
$$
and hence,
\begin{gather*}
E_2^{(1)}=\frac{f_{200}^{(3)}+2f_{110}^{(3)}\mathop{\rm Re}\{\beta\}}{p_6},\quad
E_2^{(2)}=\frac{p_4(f_{200}^{(3)}+2f_{110}^{(3)}\mathop{\rm Re}\{\beta\})}{p_5p_6},\\
E_2^{(3)}=\frac{1}{p_3p_5p_6}\det\begin{pmatrix}
-p_1&-p_2&f_{200}^{(1)}+2f_{101}^{(1)}\mathop{\rm Re}\{\beta\}\\
-p_4&-p_5&0\\
-p_6&0&f_{200}^{(3)}+2f_{110}^{(3)}\mathop{\rm Re}\{\beta\}
\end{pmatrix}.
\end{gather*}
Thus, we can determine $W_{20}(\theta)$ and $W_{11}(\theta)$ from
\eqref{c19} and \eqref{c20}. Furthermore, we can determine $g_{21}$.
Therefore, each $g_{ij}$ in \eqref{c11} is determined by the
parameters and delay in system \eqref{c1}. Thus, we can compute the
following values:
\begin{equation}\label{c25}
\begin{gathered}
c_1(0)=\frac{ i}{2\omega_0}\Big(g_{11}g_{20}-2|g_{11}|^2
-\frac{|g_{02}|^2}{3}\Big)+\frac{g_{21}}{2},\\
\mu_2=-\frac{\mathop{\rm Re}\{{c}_1(0)\}}{\mathop{\rm Re}
\{\lambda'(\tilde{\tau})\}},\quad
\beta_2=2\mathop{\rm Re}\{c_1(0)\},\\
T_2=-\frac{\mathop{\rm Im}\{{c}_1(0)\}+\mu_2\mathop{\rm Im}
\{\lambda'(\tilde{\tau})\}}{\omega_0},
\end{gathered}
\end{equation}
which determine the quantities of bifurcating periodic solutions in
the center manifold at the critical value $\tilde{\tau}$, i.e.,
$\mu_2$ determines the direction of the Hopf bifurcation: if
$\mu_2>0 (\mu_2<0)$, then the Hopf bifurcation is supercritical
(subcritical); $\beta_2$ determines the stability of the bifurcating
periodic solutions: the bifurcating periodic solutions 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{Numerical examples}

In this section, we give some examples to illustrate the results
above.

\begin{example} \label{exa1}\rm
 In system \eqref{z4}, we let
$a=4,a_1=a_2=2,a_{11}=3,b=3,r=r_1=r_2=1,m=1$,
 then system \eqref{z4} has a positive
equilibrium $E^*=(1, 3, 5)$.
 It is easy to show that
$\omega_0=0.2371, \tau_0=2.4144$. By Theorem \ref{thm2.1}, we see that the
positive equilibrium $E^*$ is stable when $\tau<\tau_0$ (see, Fig.
1); when $\tau>\tau_0$ $E^*$ is unstable (see, Fig. 2); and system
\eqref{z4} undergoes a Hopf bifurcation at $\tau_j$. When
$\tau=\tau_0$, $c_1(0)=-0.0710-0.2882 i$. It follows from
\eqref{c25} that $\mu_2>0$ and $\beta_2<0$. Therefore,  the Hopf
bifurcation of system \eqref{z4} is supercritical, and the
bifurcating periodic solutions are stable.
\end{example}

\begin{figure}[ht]
  \begin{center}
\includegraphics[width=0.4\textwidth]{fig1a} % hstable.eps
\includegraphics[width=0.4\textwidth]{fig1b} % hstable3.eps
 \end{center}
\caption{When $\tau=1<\tau_0$, the positive
equilibrium $E^*$ of system \eqref{z4} is asymptotically stable. The
initial value is $(0.5,2,4)$.}\label{F1}
\end{figure}

\begin{figure}[ht]
  \begin{center}
\includegraphics[width=0.32\textwidth]{fig2a} % 2.5.eps
\includegraphics[width=0.32\textwidth]{fig2b} % 3.eps
\includegraphics[width=0.32\textwidth]{fig2c} % hustable.eps
  \end{center}
  \caption{When $\tau=2.5,3,4>\tau_0$, bifurcating periodic
solutions from $E^*$ occur. The initial value is
$(0.5,2,4)$.}\label{F2}
\end{figure}

\begin{example} \label{exa2} \rm
 We consider the responding system without
functional response. In system \eqref{z3}, we let
$a=4,a_1=a_2=2,a_{11}=3,b=3,r=r_1=r_2=1,\beta=0$.
 It is easy to show that $E^*=(0.5, 1.5, 3.25)$,
$\omega_0=0.3722, \tau_0=2.5160$. We know that the positive
equilibrium $E^*$ of system \eqref{z3} is stable when $\tau<\tau_0$
(see, Fig. 3); when $\tau>\tau_0$, $E^*$ is unstable (see, Fig. 4),
and system \eqref{z3} undergoes a Hopf bifurcation at $\tau_j$.
\end{example}

\begin{figure}[ht]
  \begin{center}
\includegraphics[width=0.4\textwidth]{fig3a} % stable.eps
\includegraphics[width=0.4\textwidth]{fig3b} % stable3.eps
\caption{When $\tau=1<\tau_0$, the positive
equilibrium $E^*$ of system \eqref{z3} is asymptotically stable. The
initial value is $(0.5,2,4)$.}\label{F3}
  \end{center}
\end{figure}

\begin{figure}[ht]
  \begin{center}
\includegraphics[width=0.46\textwidth]{fig4a} % ustable.eps
\includegraphics[width=0.46\textwidth]{fig4b} % ustable3.eps
\end{center}
\caption{When $\tau=4>\tau_0$, bifurcating periodic
solutions from $E^*$ occur. The initial value is
$(0.5,2,4)$.}\label{F4}
\end{figure}


\subsection*{Discussion}

 In this article, we considered  a delayed predator-prey system
with stage structure and Holling type-II functional response. The
conditions of the local stability of system \eqref{z4} are obtained,
and the direction of the bifurcating periodic solutions are derived.
From Theorem \ref{thm2.1}, we know that when $Q_3>0$ the delay is harmless,
and the local stability of system \eqref{z4} doesn't change; when
$Q_3<0$, system \eqref{z4} will lose stability and a Hopf
bifurcation can occur as the delay $\tau$ increases. From fig 2, we
can see that the bifurcating periodic solutions of system \eqref{z4}
are different. The oscillatory extent of the bifurcating periodic
solution of system \eqref{z4} becomes more and more large as the
delay $\tau$ increases.

From the figs, we known that the rapidity of convergence to
equilibrium of system \eqref{z4} is slower than corresponding system
without functional response. Furthermore, the value of $\tau_0$
obtained for the model \eqref{z4} is smaller than the corresponding
value for a similar model without functional response. At the same
time stage-structured system with time delay and functional response
has a similar asymptotic behavior to that in the delayed system
without functional response.


\begin{thebibliography}{00}

\bibitem{chenl}L. Chen, X. Song, Z. Lu;
Mathematical Models and Methods in Ecology,
\emph{Sichuan Science and Technology Press}, 2003.

\bibitem{cuij}J. Cui, L. Chen, W. Wang;
 The effect of dispersal on population growth with stage-structure,
  \emph{Appl. Math. Comput.},  39 (2000), 91-102.

\bibitem{cuis}J. Cui, X. Song;
Permanence of predator-prey system with stage structure,
\emph{Discrete Contin. Dyn. Syst. Ser. B,}  4 (2004), 547-554.

\bibitem{gm}M. Gamez, C. Martinez;
Persistence and global stability in a
predator- prey system with delay, \emph{International Journal of
Bifurcation and Chaos}, 16 (2006), 2915-2922.

\bibitem{gsj} S. Gao, L. Chen, Z. Teng;
Hopf bifurcation and global stability for
a delayed predator-prey system with stage structure for predator,
\emph{Appl. Math. Comput.},  202 (2008), 721-729.

\bibitem{hale}J. Hale, D. Lunel;
Introduction to Functional Differential Equations, \emph{
Springer-Verlag}, New York, 1993.

\bibitem{hasd}B. Hassard, D. Kazarinoff;
Theory and Application of Hopf Bifurcation, \emph{
Cambridge University Press}, 1981.

\bibitem{mzp} Z. Ma, H. Huo, C. Liu;
Stability and Hopf bifurcation analysis on a predator-prey system
with discrete and distributed delays, \emph{Nonlinear Analysis: Real
World Applications}, 10 (2009), 1160-1172.

\bibitem{wang}J. Wang, K. Wang;
The optimal harvesting problems of a stage-structured population,
 \emph{Appl. Math. Comput.}, 148 (2004), 235-247.
\bibitem{rx}R. Xu, M. A. J. Chaplain, F. A.Davidson;
Global stability of a Lotka-Volterra type predator-prey
model with stage structure and time delay, \emph{Appl. Math.
Comput.},  159 (2004), 863-880.

\bibitem{xr2}R. Xu, Z. Ma;
The effect of stage-structure on the permanence
of a predator-prey system with time delay,
 \emph{Appl. Math. Comput.},  189 (2007), 1164-1177.

\bibitem{zh} H. Zhang, L. Chen, R. Zhu;
 Permanence and extinction of a
periodic predator-prey delay system with functional response and
stage structure for prey, \emph{Appl. Math. Comput.}, 184 (2007)
931-944.

\end{thebibliography}
\end{document}