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

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

\begin{document}
\title[\hfilneg EJDE-2013/09\hfil Stability and Hopf bifurcation]
{Stability and Hopf bifurcation in a symmetric
Lotka-Volterra predator-prey system with delays}

\author[Jing Xia, Zhixian Yu, Rong Yuan \hfil EJDE-2013/09\hfilneg]
{Jing Xia, Zhixian Yu, Rong Yuan}  % in alphabetical order

\address{Jing Xia \newline
School of Mathematical Sciences, 
Peking University, Beijing 100871, China}
\email{xiajing2005@mail.bnu.edu.cn}

\address{Zhixian Yu \newline
College of Science, University of Shanghai for Science and Technology,
Shanghai 200093, China}
\email{yuzx0902@yahoo.com.cn}

\address{Rong Yuan \newline
School of Mathematical Sciences,
Beijing Normal University,
Beijing 100875, China}
\email{ryuan@bnu.edu.cn}

\thanks{Submitted September 28, 2012. Published January 9, 2013.}
\subjclass[2000]{34K18, 37G05, 37G10, 92D25}
\keywords{Predator-prey system; delay; stability; Hopf
bifurcation; normal form}

\begin{abstract}
 This article concerns a symmetrical Lotka-Volterra predator-prey system
 with delays. By analyzing the associated characteristic equation of
 the original system at the positive equilibrium and choosing
 the delay as the bifurcation parameter, the local stability and Hopf
 bifurcation of the system are investigated. Using the normal form
 theory, we also establish the direction and stability of the Hopf bifurcation.
 Numerical simulations suggest an existence of Hopf bifurcation near
 a critical value of time delay.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
%\newtheorem{proposition}[theorem]{Proposition}
%\newtheorem{corollary}[theorem]{Corollary}
%\newtheorem{definition}[theorem]{Definition}
%\newtheorem{example}[theorem]{Example}
%\newtheorem{remark}[theorem]{Remark}
\allowdisplaybreaks

\def\C{\mathbb C}
\def\R{\mathbb R}
\def\X{\mathbb X}
\def\cA{\mathcal A}
\def\cT{\mathcal T}
\def\Z{\mathbb Z}
\def\Y{\mathbb Y}
\def\Z{\mathbb Z}
\def\N{\mathbb N}

\section{Introduction}


 Since the pioneering work of Volterra \cite{Volterra}, the
delayed predator-prey models of Lotka-Volterra type have
played an important role in the development of population dynamics.
For instance, the predator-prey system
\begin{equation}\label{1.10}
\begin{gathered}
\dot{x}(t)=x(t)[r_1-a_{11}x(t-\tau)-a_{12}y(t)],\\
\dot{y}(t)=y(t)[-r_2+a_{21}x(t)-a_{22}y(t)]
\end{gathered}
\end{equation}
was first proposed and discussed briefly by May \cite{May} in 1973.
The dynamics of system \eqref{1.10} or systems similar to (1.1) has
been widely investigated, concerning boundedness of solutions,
stability, permanence, persistence, periodic phenomena, bifurcation
and chaos (see, for example, \cite{Beretta,Kuang,Nakao,Ruan,Saito}
and their references).
In particular, research on the stability of positive
equilibrium and periodic solutions arising from the Hopf bifurcation
is very critical, which helps us to realize the law and trend for
species quantity. Also, extensive models of various fields (such as
biology, neural network, engineering systems, epidemic disease) have
been embodied these properties,
 (see \cite{Faria3,Liu,Wei,Sun}).

In mathematical biology, Lotka-Volterra type predator-prey with a
delay have been well studied by Song and Wei \cite{song}, Yan
 and Li \cite{yanLi}, Yan and Zhang \cite{YanZhang} and so on. In addition,
 models involving two discrete delays are increasingly applied to
 the study of a variety of situations. Faria \cite{Faria3}
 considered the predator-prey system of the form
\begin{equation}\label{1.0}
\begin{gathered}
\dot{x}(t)=x(t)[r_1-a_{11}x(t)-a_{12}y(t-\tau_1)],\\
\dot{y}(t)=y(t)[-r_2+a_{21}x(t-\tau_2)-a_{22}y(t)].
\end{gathered}
\end{equation}
The author took $\tau_2$ as a bifurcation parameter and obtained the
existence of local Hopf bifurcation. Moreover, Song et. al took the
sum of two delays $\tau_1+\tau_2$ as a parameter and showed that
when $\tau$ passed through some critical values, Hopf bifurcation
occurs.

 In the present article, we focus on the following symmetrical
Lotka-Volterra predator-prey system:
\begin{equation}\label{1.20}
\begin{gathered}
\dot{x}(t)=x(t)[r_1-ax(t)-bx(t-\tau_1)-cy(t-\tau_2)],\\
\dot{y}(t)=y(t)[r_2-ay(t)+cx(t-\tau_1)-by(t-\tau_2)],
\end{gathered}
\end{equation}
where $r_1$, $r_2$, $a$, $b$, $c$, $\tau_1$, $\tau_2$ are constants
with $a>0$, $b>0$, $c>0$ and $\tau_1\geq 0$, $\tau_2\geq 0$. The
variables $x(t)$ and $y(t)$ denote the population densities of prey
and predator at time t, respectively. The initial conditions for
system \eqref{1.20} are
\begin{gather*}
x(\theta)=\phi_1(\theta),\quad -\tau_1\leq\theta\leq 0,\; \phi_1(0)>0,\\
y(\theta)=\phi_2(\theta),\quad -\tau_2\leq\theta\leq 0,\; \phi_2(0)>0.
 \end{gather*}
Obviously, in the case $b=0$ and $a\neq b$, systems \eqref{1.20} is
reduced to system \eqref{1.0}. Systems \eqref{1.20} shows that, in
the absence of predator species, the prey species are governed by
general case of Hutchinson's equation
\begin{equation}\label{1.30}
\dot{x}(t)=x(t)[r_1-ax(t)-bx(t-\tau_1)].
\end{equation}
(There is an extensive literature on various aspects of
\eqref{1.30}, see, e.g., \cite{Gopalsamy, sunhan}). In the
presence of predators, there is a hunting term, $cy(t-\tau_2)$, with
a certain delay $\tau_2$, called the hunting delay. In the absence
of preys, the feedback $cx(t-\tau_1)$ has a nonnegative delay which
is the time of the predator maturation. For system \eqref{1.20},
Saito et al studied the permanence and the global asymptotic
stability of positive equilibrium including the cases $b\leq 0$. Our
main purpose is to study stability and Hopf bifurcation of
\eqref{1.20}.

This article is organized as follows. In Section 2, the local
stability of the positive equilibrium of system \eqref{1.20} is
addressed. In Section 3, we show that the system \eqref{1.20}
with $\tau_1=\tau_2$
undergoes Hopf bifurcation as the delays cross some critical values.
In Section 4, following the procedure of deriving normal form due to
Faria and Magalh\~{a}es \cite{Faria1,Faria2}, the direction and
stability of the Hopf bifurcation are determined. Finally, we give
some numerical simulations to justify the theoretical analysis.

\section{Local stability}

 The following lemma will be needed in analyzing the
 characteristic equation of the linearize system.

 \begin{lemma}[\cite{Cooke}] \label{lem0}
 Let $f(\lambda,\tau)=\lambda^2+a\lambda+b\lambda e^{-\tau \lambda}
+c+de^{-\tau \lambda}$, where $a$, $b$, $c$, $d$, $\tau$
are real numbers and $\tau\geq 0$. Then, as $\tau$ varies,
the sum of the multiplicities of zeros of $f$ in the open right
half-plane can change only if a zero appears on or crosses the imaginary axis.
\end{lemma}

 In the rest of this paper, we assume that
 system \eqref{1.20} has a unique positive equilibrium
$E (x^*,y^*)$:
\[
x^*=\frac{r_1(a+b)-r_2c}{(a+b)^2+c^2}>0, \quad
y^*=\frac{r_2(a+b)+r_1c}{(a+b)^2+c^2}>0,
\]
and $x^*\neq y^*$.

By using the transformation $\overline{x}(t)=x(t)-x^*$,
$\overline{y}(t)=y(t)-y^*$, and dropping the bars for simplicity
of notation, system \eqref{1.20} can be transformed into the
 system
\begin{equation}\label{1.40}
\begin{aligned}
 \dot{x}(t)
&=-ax^*x(t)-bx^*x(t-\tau_1)-cx^*y(t-\tau_2)-ax^2(t)\\
 &\quad -bx(t)x(t-\tau_1)-c x(t)y(t-\tau_2), \\
\dot{y}(t)&=-ay^*y(t)+cy^*x(t-\tau_1)-by^*y(t-\tau_2)-ay^2(t)\\
&\quad +c x(t-\tau_1)y(t)-by(t)y(t-\tau_2).
 \end{aligned}
\end{equation}
Then we obtain the linearized system
\begin{equation} \label{1.50}
\begin{gathered}
 \dot{x}(t)=-ax^*x(t)-bx^*x(t-\tau_1)-cx^*y(t-\tau_2),\\
 \dot{y}(t)=-ay^*y(t)+cy^*x(t-\tau_1)-by^*y(t-\tau_2).\\
 \end{gathered}
\end{equation}
The associated characteristic equation of system \eqref{1.50} is
\begin{equation}\label{1.60}
\begin{split}
&\lambda ^2+\lambda (ax^*+ay^*+bx^* e^{-\lambda
\tau_1}+by^*e^{-\lambda \tau_2})+a^2x^*y^*
+abx^*y^*e^{-\lambda \tau_2}\\
&+abx^*y^*e^{-\lambda \tau_1}
+(b^2x^*y^*+c^2x^*y^*)e^{-\lambda (\tau_1+\tau_2)}=0.
\end{split}
\end{equation}
Obviously, $\lambda=0$ is not a root of \eqref{1.60}. When
$\tau_1=\tau_2=0$, the characteristic equation \eqref{1.60} becomes
\begin{equation}\label{1.70}
 \lambda ^2+\lambda (ax^*+ay^*+bx^*+by^*)+a^2x^*y^*+2abx^*y^*
+b^2x^*y^*+c^2x^*y^*=0.
\end{equation}
Letting $\lambda_1$, $\lambda_2$ be the characteristic roots of
\eqref{1.70}, then $\lambda_1+\lambda_2=-(ax^*+ay^*+bx^*+by^*)<0$,
$\lambda_1\lambda_2=a^2x^*y^*+2abx^*y^* +b^2x^*y^*+c^2x^*y^*>0$.
Thus, all roots of \eqref{1.70} have negative real parts.

When $\tau_2=0$, Equation \eqref{1.60} becomes
\begin{equation}\label{1.80}
 \lambda ^2+p\lambda +q+s e^{-\lambda
\tau_1}+l\lambda e^{-\lambda \tau_1}=0,
\end{equation}
where $p=ax^*+ay^*+by^*$, $q=(a^2+ab)x^*y^*$,
$s=(ab+b^2+c^2)x^*y^*$, $l=bx^*$. We first study \eqref{1.80}
with the delay $\tau_1$. We consider the stable
intervals of $\tau_1$ by taking $\tau_1$ as a parameter and applying
the lemma in \cite{Cooke}. Secondly, we investigate \eqref{1.60} in
the stable interval of $\tau_1$. Choosing $\tau_2$ as a parameter
and using the lemma in \cite{Cooke} once more, we will find a stable
interval (depending on $\tau_1$) for $\tau_2$. This will be the
stable interval for \eqref{1.60}.

By putting $\tau_1=0$ in \eqref{1.80}, we obtain
\[
 \lambda ^2+\lambda(l+p)+s+q=0.
\]
Based on the above analysis, we know that two roots of the above
equation have negative real parts. In addition, $\lambda=i\sigma$ is
a root of \eqref{1.80} if and only if $\sigma$ satisfies:
\[
 -\sigma^2+i\sigma p+q+se^{-i\sigma
\tau_1}+li\sigma e^{-i\sigma \tau_1}=0.
\]
Separating the real and imaginary parts, we have
 \begin{gather*}
 \sigma^2-q=s\cos\sigma\tau_1+l\sigma\sin\sigma\tau_1, \\
p \sigma=s\sin\sigma\tau_1-l\sigma\cos\sigma\tau_1,\\
 \end{gather*}
which leads to
\[
\sigma^4+\sigma^2(p^2-2q-l^2)+q^2-s^2=0.
\]

In this article we will use the following hypotheses:
\begin{itemize}
\item[(H1)] $q^2<s^2$,

\item[(H2)] $q^2>s^2$, $2q+l^2-p^2>0$, $(2q+l^2-p^2)^2>4(q^2-s^2)$.

\end{itemize}
Then we easily obtain the following results.

\begin{lemma}\label{lem2.10}
For \eqref{1.80}, we have
\begin{itemize}
\item[1.] If {\rm (H1)} holds, then when $\tau_1=\tau_{1,k}^{(1)}$, Equation
\eqref{1.80} has a pair of purely imaginary roots $\pm i\sigma_+$;

\item[2.] If {\rm (H2)} holds, then when
 $\tau_1=\tau_{1,k}^{(1)}(\tau_1=\tau_{1,k}^{(2)})$, Equation
\eqref{1.80} has a pair of purely imaginary roots
$\pm i\sigma_+(\pm i\sigma_-)$;

\item[3.] If neither {\rm (H1)} nor {\rm (H2)} hold, then \eqref{1.80} has no
purely imaginary root,
\end{itemize}
where
 \begin{equation} \label{1.90}
 \begin{gathered}
 \sigma_{\pm}^2=\frac{2q+l^2-p^2\pm\sqrt{(2q+l^2-p^2)^2-4(q^2-s^2)}}{2},\\
\tau_{1,k}^{(1)}=\frac{1}{\sigma_+}\arccos\frac{\sigma_+^2s-qs
 -lp\sigma_+^2}{s^2+l^2\sigma_+^2}+ \frac{2k\pi}{\sigma_+},\\
\tau_{1,k}^{(2)}=\frac{1}{\sigma_-}\arccos\frac{\sigma_-^2s-qs
 -lp\sigma_-^2}{s^2+l^2\sigma_-^2}+ \frac{2k\pi}{\sigma_-}
\end{gathered}
 \end{equation}
\end{lemma}

Let
$\lambda_{n,k}=\mu_{n,k}(\tau_1)+i\sigma_{n,k}(\tau_1)$($n=1,2$) denote the
roots of \eqref{1.80} satisfying
\begin{gather*}
\mu_{1,k}(\tau_{1,k}^{(1)})=0,\quad\sigma_{1,k}(\tau_{1,k}^{(1)})=\sigma_+,\\
\mu_{2,k}(\tau_{1,k}^{(2)})=0,\quad\sigma_{2,k}(\tau_{1,k}^{(2)})=\sigma_-.
 \end{gather*}
In what follows, we can verify the following transversality conditions.

 \begin{lemma}\label{lem20}
Assume that either {\rm (H1)} or {\rm (H2)} holds, then
\begin{equation}\label{1.100}
\frac{d\operatorname{Re}\lambda_{1,k}(\tau_{1,k}^{(1)})}{d\tau_1}>0,\quad
\frac{d\operatorname{Re}\lambda_{2,k}(\tau_{1,k}^{(2)})}{d\tau_1}<0.
\end{equation}
\end{lemma}

\begin{proof}
 Differentiating the characteristic equation \eqref{1.80} with
respect to $\tau_1$, we obtain
\[
\big(\frac{d\lambda}{d\tau}\big)^{-1}
=-\frac{2\lambda+p}{\lambda(\lambda^2+p\lambda+q)}
+\frac{l}{\lambda^2l+\lambda s}-\frac{\tau_1}{\lambda}.
\]
Thus, if $\lambda=i\sigma_+$, then
\begin{align*}
&\operatorname{sign}\Big\{\operatorname{Re}(\frac{d\lambda}{d\tau_1})\big|_{\tau_1
 =\tau_{1,k}^{(1)}}\Big\}
=\operatorname{sign}\Big\{\operatorname{Re}(\frac{d\tau}{d\lambda})
 \big|_{\lambda=i\sigma_+}\Big\}
\\
&=\operatorname{sign}\Big \{
\frac{\sigma_+^4l^2+2\sigma_+^2s^2-2qs^2+p^2s^2-l^2q^2}{[(q-\sigma_+^2)^2
 +p^2\sigma_+^2](\sigma_+^2l^2+s^2)}\Big\}.
\end{align*}
Since
\begin{align*}
&\sigma_+^4l^2+2\sigma_+^2s^2-2qs^2+p^2s^2-l^2q^2\\
&=\frac{1}{2}l^2\sqrt{(2q+l^2-p^2)^2-4(q^2-s^2)}\\
&\quad\times [\sqrt{(2q+l^2-p^2)^2-4(q^2-s^2)}+
(2q+l^2-p^2)]\\
&\quad +s^2\sqrt{(2q+l^2-p^2)^2-4(q^2-s^2)},
\end{align*}
thus, either (H1) or (H2) implies that
$\operatorname{sign}
\big\{\operatorname{Re}(\frac{d\lambda}{d\tau_1})\big|_{\tau_1=\tau_{1,k}^{(1)}}
\big\}>0$.
If $\lambda=i\sigma_-$, we can similarly prove that
\[
\operatorname{sign}\Big\{\operatorname{Re}(\frac{d\lambda}{d\tau_1})\big|_{\tau_1=\tau_{1,k}^{(2)}}\Big\}
<0.
\]
The proof is complete.
\end{proof}


From Lemmas \ref{lem2.10} and \ref{lem20}, we have the following result.

\begin{lemma} \label{lem2.3}
For Equation \eqref{1.80}, the following statements are true:

1. Assume that {\rm (H1)} holds, then for any $\tau_1\in[0,\tau_{1,0}^{(1)})$,
all roots of \eqref{1.80} have negative real parts, and for
 $\tau_1>\tau_{1,0}^{(1)}$,
Equation \eqref{1.80} has at least one root with positive real parts.

2. Assume that {\rm (H2)} holds, then there exists an integer $k_1>0$
such that there are $k_1$ switches from stability to instability to stability;
that is, when $\tau_1\in (\tau_{1,k}^{(2)},\tau_{1,k+1}^{(1)})$,
$k=-1,0,1,\dots,k_1-1$, all roots of \eqref{1.80} have negative real parts,
 where $\tau_{1,-1}^{(2)}=0$;
When $\tau_1\in (\tau_{1,k}^{(1)},\tau_{1,k}^{(2)})$ and
$\tau_1>\tau_{1,k_1}^{(1)}$,
$k=0,1,\dots,k_1-1$, Equation \eqref{1.80} has at least one root with
positive real parts.
\end{lemma}

In what follows, we consider \eqref{1.60} with $\tau_1$ in its stable interval.
We will take $\tau_2$
as a bifurcation parameter and have the following result.

\begin{lemma}\label{lem40}
Assume that all roots of \eqref{1.80} have negative real parts,
 then there exists $\tau_2(\tau_1)>0$ such that
when $\tau_2\in[0,\tau_2(\tau_1))$, all roots of \eqref{1.60} have negative
real parts.
\end{lemma}

\begin{proof}
Since all roots of \eqref{1.80} have negative real parts,
all roots of \eqref{1.60} with $\tau_2=0$
have negative real parts. Similar to the proof in \cite{Cooke},
 let $f(\lambda, \tau_1, \tau_2)$ be the left-hand side of \eqref{1.60}.
 Obviously, $f(\lambda, \tau_1, \tau_2)$ is analysis in $\lambda$ and $\tau_2$.
 We take $\tau_2$ as a bifurcation parameter,
 when $\tau_2$ varies, the sum of the multiplicities of
zero of $f(\lambda,\tau_1,\tau_2)$ in the open right half-plane can change
only if a zero appears on or crosses the imaginary axis.
Thus, there exists $\tau_2(\tau_1)>0$, such that when
$\tau_2\in \tau_2(\tau_1)$, all roots of \eqref{1.60}
have negative real parts. The proof is complete.
\end{proof}


Following Lemma \ref{lem40}, we can get sufficient conditions
for local stability of system \eqref{1.40}.

\begin{theorem}\label{thm10}
For system \eqref{1.40}, the following statements hold.

1. Suppose that {\rm (H1)} holds, then for any $\tau_1\in[0,\tau_{1,0}^{(1)})$,
there is a $\tau_2(\tau_1)>0$ such that
when $\tau_2\in[0,\tau_2(\tau_1))$, the zero solution of system \eqref{1.40}
 is locally asymptotically stable.


2. Suppose that {\rm (H2)} holds, then for any $\tau_1\in
(\tau_{1,k}^{(2)},\tau_{1,k+1}^{(1)})$, $k=-1,0,1,\dots,k_1-1$,
there is $\tau_2(\tau_1)>0$ such that
when $\tau_2\in[0,\tau_2(\tau_1))$, the zero solution of system \eqref{1.40}
 is locally asymptotically stable.

3. Suppose that neither {\rm (H1)} nor {\rm (H2)} holds, then for any
 $\tau_1\geq 0$, there is a $\tau_2(\tau_1)>0$, such that
when $\tau_2\in[0,\tau_2(\tau_1))$, the zero solution of system \eqref{1.40}
 is locally asymptotically stable.
\end{theorem}


\section{Hopf bifurcation}

Let $\tau_1=\tau_2=\tau$. System \eqref{1.20} becomes
\begin{equation}\label{90}
 \begin{gathered}
\dot{x}(t)=x(t)[r_1-ax(t)-bx(t-\tau)-cy(t-\tau)],\\
\dot{y}(t)=y(t)[r_2-ay(t)+cx(t-\tau)-by(t-\tau)].
 \end{gathered}
\end{equation}
and the characteristic equation \eqref{1.60} has the form
\begin{equation}\label{100}
 \lambda ^2+p_1\lambda +q_1+(s_1\lambda+l_1) e^{-\lambda
\tau}+m_1 e^{-2\lambda \tau}=0,
\end{equation}
where
\begin{equation}\label{105}
\begin{gathered}
p_1=ax^*+ay^*, \quad q_1=a^2x^*y^*,\quad l_1=2abx^*y^*, \\
 s_1=bx^*+by^*,\quad m_1=b^2x^*y^*+c^2x^*y^*.
\end{gathered}
\end{equation}

\begin{lemma}\label{lem50}
If $b^2+c^2>a^2$, then there exists $\sigma_+$ $(>0)$ such that
\eqref{100} with $\tau=\tau_k$ has a simple pair of conjugate purely
imaginary roots $\pm i\sigma_+$, where
\[
\tau_k=\frac{1}{\sigma_+}\arccos\frac{(l_1-s_1p_1)\sigma_+^2-l_1(q_1-m_1)}
{(-\sigma_+^2+q_1)^2-m_1^2+p_1^2\sigma_+^2}+\frac{2k\pi}{\sigma_+},
\quad k=0,1,\dots.
\]
\end{lemma}

\begin{proof}
Obviously, $\lambda=0$ is not a root of \eqref{100}.
If $\lambda=i\sigma$ is a root of \eqref{100}, then
\begin{equation}\label{110}
-\sigma^2+ip_1\sigma+q_1+l_1e^{-i\sigma\tau}+i\sigma
s_1e^{-i\sigma\tau}+m_1e^{-2i\sigma\tau}=0.
\end{equation}
Multiplying both sides of \eqref{110} by $e^{i\sigma\tau}$, we have
\begin{equation}\label{120}
(-\sigma^2+ip_1\sigma+q_1)e^{i\sigma\tau}+(l_1+i\sigma
s_1)+m_1e^{-i\sigma\tau}=0.
\end{equation}
Separating the real parts and imaginary parts of \eqref{120} leads to
\begin{equation}\label{130}
 \begin{gathered}
[(-\sigma^2+q_1)^2-m_1^2+p_1^2\sigma^2]\cos\sigma\tau
=(l_1-s_1p_1)\sigma^2-l_1(q_1-m_1), \\
[(-\sigma^2+q_1)^2-m_1^2+p_1^2\sigma^2]\sin\sigma\tau
=s_1\sigma^3+[l_1p_1-s_1(q_1+m_1)]\sigma,\\
\end{gathered}
\end{equation}
and thus,
\begin{align*} %\label{140}
&[(-\sigma^2+q_1)^2-m_1^2+p_1^2\sigma^2]^2\\
&=[(l_1-s_1p_1)\sigma^2-l_1(q_1-m_1)]^2
 +[s_1\sigma^3+[l_1p_1-s_1(q_1+m_1)]\sigma]^2.
\end{align*}
Let $y=\sigma^2$, define
\begin{equation}\label{150}
\begin{split}
 f(y)&=[(-y+q_1)^2-m_1^2+p_1^2y]^2-[(l_1-s_1p_1)y-l_1(q_1-m_1)]^2\\
 &\quad -y[s_1y+l_1p_1-s_1(q_1+m_1)]^2.
\end{split}
\end{equation}
Since $f(y)$ is a quartic function, we have $f(y)\to +\infty$ as $|y|\to
+\infty$. On the other hand,
$$
f(0)=(q_1^2-m_1^2)^2-l_1^2(q_1-m_1)^2=(q_1-m_1)^2(q_1+m_1+l_1)(q_1+m_1-l_1)>0,
$$
Define
\begin{gather*}
F(y)=[(-y+q_1)^2-m_1^2+p_1^2y]^2,\\
G(y)=-[(l_1-s_1p_1)y-l_1(q_1-m_1)]^2,\\
H(y)=-y[s_1y+l_1p_1-s_1(q_1+m_1)]^2,
\end{gather*}
Then $f(y)=F(y)+G(y)+H(y)$.
Assume that $g(y)=(-y+q_1)^2-m_1^2+p_1^2y$, then
$$
g(0)=q_1^2-m_1^2=(q_1+m_1)(q_1-m_1)=(q_1+m_1)x^*
y^*(a^2-b^2-c^2),
$$
In view of $b^2+c^2>a^2$, we have $g(0)<0$. On the other hand,
$g(y)\to +\infty$ as $y\to +\infty$.
Thus, there is a $z_0>0$ such that $g(z_0)=0$;
 i.e., $F(z_0)=0$. It is easily to see that
positive roots of $F$, $G$ and $H$ are mutually different when $x^*\neq y^*$.
 Consequently, $G(z_0)<0$, $H(z_0)<0$. Thus, $f(z_0)<0$, which together with
 $f(0)>0$ imply that there exists a $y_0>0$ such that $f(y_0)=0$.
Obviously, there exists another positive root because $f(y)\to+\infty$ as
 $y\to +\infty$. Furthermore, one of the two roots is a simple root at least.
Let $y_0$ be such a simple root and
$y_0=\sigma^2$. Substituting $y_0=\sigma^2$ into \eqref{130}, we have
\begin{equation}\label{160}
\tau_k=\frac{1}{\sigma_+}\arccos\frac{(l_1-s_1p_1)\sigma_+^2-l_1(q_1-m_1)}
{(-\sigma_+^2+q_1)^2-m_1^2+p_1^2\sigma^2}+\frac{2k\pi}{\sigma_+},\quad
k=0,1,\dots.
\end{equation}
The proof is complete.
\end{proof}

Let $\lambda(\tau)=\mu(\tau)+i\sigma(\tau)$ be the root of \eqref{100}
near $\tau=\tau_k$ satisfying $\mu(\tau_k)=0$, $\sigma(\tau_k)=\sigma_+$,
and $k=0,1,2,\dots$, where $\tau_k$ is defined by \eqref{160}.

\begin{lemma}[\cite{Saito}] \label{lem60}
Assume that $b^2+c^2>a^2$ and $a>b$, then
$$
\operatorname{sign} \Big\{\operatorname{Re}(\frac{d\lambda}{d\tau})\big
|_{\tau=\tau_k}\Big\}>0.
$$
\end{lemma}

The proof of the above lemma is similar to that in \cite[p. 549]{Saito}
and is omitted.
By Lemma \ref{lem50} and \ref{lem60}, we have the following result.

\begin{theorem}\label{th60}
Assume that $b^2+c^2>a^2$ and $a>b$. Then
\begin{itemize}
\item[1.] When $\tau\in [0,\tau_0)$, the equilibrium $E=(x^*,y^*)$
 of \eqref{90}is locally asymptotically stable;

\item[2.] When $\tau>\tau_0$, the equilibrium $E=(x^*,y^*)$ of \eqref{90}
 is unstable;

\item[3.] System \eqref{90} undergoes a Hopf bifurcation at $E=(x^*,y^*)$
and $\tau=\tau_k$.
\end{itemize}
\end{theorem}

\section{Direction and stability of the Hopf bifurcation}

In this section, we shall determine the
direction and stability of the bifurcating periodic solutions by
employing the normal form theory duo to Faria and
Magalh\~aes \cite{Faria1, Faria2}.

Let $z_1(t)=x(\tau t)-x^*$, $z_2(t)=y(\tau t)-y^*$.
Then system \eqref{90} can be
rewritten as a functional differential equation in
$C([-1,0],\mathbb{R}^2)$,
\begin{equation}\label{3.10}
\begin{gathered}
\begin{aligned}
 \dot{z}_1(t)&=x^*\tau[-az_1(t)-bz_1(t-1)-cz_2(t-1)]-\tau z_1(t)[az_1(t)\\
&\quad +bz_1(t-1)+cz_2(t-1)],
\end{aligned} \\
\begin{aligned}
\dot{z}_2(t)&=y^*\tau[cz_1(t-1)-az_2(t)-bz_2(t-1)]+\tau z_2(t)[cz_1(t-1)\\
&\quad -az_1(t)-bz_2(t-1)].
 \end{aligned}
\end{gathered}
\end{equation}
and its linearized system is
\begin{equation}\label{3.20}
 \begin{gathered}
 \dot{z}_1(t)=\tau[-ax^*z_1(t)-bx^*z_1(t-1)-cx^*z_2(t-1)], \\
\dot{z}_2(t)=\tau[cy^*z_1(t-1)-ay^*z_2(t)-by^*z_2(t-1)].\\
 \end{gathered}
\end{equation}
The characteristic equation associated with system \eqref{3.20} has
the form
\begin{equation}\label{3.30}
\lambda^2+\lambda \tau (ay^*+ax^*+bx^*e^{-\lambda}+by^*e^{-\lambda})
+\tau ^2x^*y^*(a^2+2abe^{-\lambda}+b^2e^{-2\lambda}+c^2e^{-2\lambda})=0.
\end{equation}
By the change of variable $\lambda=\xi
\tau$, Equation \eqref{3.30} becomes
\begin{equation}\label{3.40}
\xi^2+p_1\xi+q_1+(s_1\xi+l_1)e^{-\xi \tau}+m_1e^{-2\xi \tau}=0,
\end{equation}
where $p_1,q_1,s_1,l_1,m_1$ are defined by \eqref{105}.
Equation \eqref{3.40} is the same as \eqref{100}. Therefore, from
Section 3, we know that \eqref{3.40} has a simple pair of
conjugate complex roots $\xi(\tau)=\mu(\tau) \pm i\sigma(\tau)$
which satisfy
\[
\mu(\tau_k)=0,\quad \sigma(\tau_k)=\sigma_+,\quad \mu'(\tau_k)\neq
0,
\]
 where
$\tau_k$ is given in \eqref{160}. Thus, Eq. \eqref{3.30} has two
complex roots $\lambda(\tau)=\tau \mu(\tau)\pm i \tau \sigma(\tau)$,
which satisfy
\[
\frac{d\operatorname{Re}\lambda(\tau)}{d\tau}|_{\tau=\tau_k}= \tau
\frac{d\operatorname{Re}\xi(\tau)}{d\tau}|_{\tau=\tau_k}
+\operatorname{Re}\xi(\tau)|_{\tau=\tau_k}=\tau_k\mu'(\tau_k)\neq
0.
\]
 Letting $z=(z_1,z_2)^T$, then system \eqref{3.10} can further be written as
\begin{equation}\label{3.50}
\dot{z}(t)=L(\tau)z_t+F(z_t,\tau),
\end{equation}
where
\begin{gather*}
L(\tau)\varphi = \tau \begin{pmatrix}
 -ax^*\varphi_1(0)-bx^*\varphi_1(-1)-cx^*\varphi_2(-1)\\
 cy^*\varphi_1(-1)-ay^*\varphi_2(0)-by^*\varphi_2(-1)
 \end{pmatrix},
\\
F(\varphi,\tau)= \tau \begin{pmatrix}
 -a\varphi_1^2(0)-b\varphi_1(0)\varphi_1(-1)-c\varphi_1(0)\varphi_2(-1)\\
 -a\varphi_2^2(0)+c\varphi_1(-1)\varphi_2(0)-b\varphi_2(0)\varphi_2(-1)
\end{pmatrix},
\end{gather*}
 and $\varphi=(\varphi_1,\varphi_2)^T\in C([-1,0],\mathbb{R}^2)$.

Introducing a new parameter $\alpha=\tau-\tau_k$,
system \eqref{3.50} is equivalent to
\begin{equation}\label{3.70}
\dot{z}(t)=L(\tau_k)(z_t)+F_0(z_t,\alpha),
\end{equation}
where $F_0(z_t,\alpha)=L(\alpha)z_t+F(z_t,\tau_k+\alpha)$.

 Let $A_0$ be the infinitesimal generator corresponding
to $\dot{z}(t)=L(\tau_k)z_t$. Then $A_0$ has a pair of purely
imaginary characteristic roots
$\pm i\sigma_k(\sigma_k=\tau_k\sigma_+)$, which are simple, and no other
characteristic roots with zero real part. Define
$\Lambda=\{-i\sigma_k,i\sigma_k\}$ and denote by $P$ the invariant
space of $A_0$ corresponding to $\Lambda$, where the dimension of
$P$ equals to 2. The phase space is $C=C([-1,0],\mathbb{R}^2)$, so
we can decompose the space $C$ as $C=P \oplus Q$ by applying the
formal adjoint theory for functional differential equations in
\cite{jh77}. Consider complex coordinates and still write as
$C([-1,0],\mathbb{C}^2)$. Suppose that $\Phi_1,\Phi_2$ are
the eigenvectors of the operator $A_0$ associated with eigenvalues
$i\sigma_k,-i\sigma_k$, respectively. Thus, $\Phi=(\Phi_1,\Phi_2)$
is a basis of $P$, $\Phi_1(\theta)=e^{i\sigma_k\theta} v^T
=e^{i\sigma_k\theta}(1,v_2)^T$,
$\Phi_2(\theta)=\overline{\Phi_1(\theta)}$, $-1\leq \theta \leq 0,$
where $v=(v_1, v_2)^T$ is a vector in $\mathbb{C}^2$ and
$$
v_2=-\frac{i\sigma_k+\tau_k a x^*+\tau_kbx^*
e^{-i\sigma_k}}{\tau_kcx^*e^{-i\sigma_k}},
$$
which satisfies
\[
L(\tau_k)(\Phi_1)=i\sigma_k v^T.
\]
For $\Phi=(\Phi_1,\Phi_2)$, we note that
$\dot{\Phi}=\Phi B$, where $B$ is a diagonal matrix of the form
 $\operatorname{diag}(i\sigma_k,-i\sigma_k)$.
 Also, the formal adjoint
operator $A_0^*$ of $A_0$ has eigenvectors $\Psi_1$, $\Psi_2$
corresponding to $i\sigma_k$, $-i\sigma_k$.
$\Psi(s)=\operatorname{col}(\Psi_1(s),\Psi_2(s))$ constructs
a basis of the adjoint
space $P^*$ of $P$ and $\Psi_1(s)=e^{-i\sigma_k s}(u_1,u_2)$,
$\Psi_2(s)=\overline{\Psi_1(s)}$, $0\leq s\leq 1$. We know
$-\dot{\Psi}_1(0)=\int_{-1}^0\Psi_1(-\theta)d\eta(\theta)$,
therefore, we can choose
$$
u_2=\frac{i\sigma_k+\tau_k a x^*+\tau_kbx^*
e^{-i\sigma_k}}{\tau_kcy^*e^{-i\sigma_k}}u_1.
$$
 In order to have $(\Psi,\Phi)=((\Psi_j,\Phi_i),i,j=1,2)=I_2$; i.e., second
order identical matrix, where $(^.,^.)$ is a bilinear inner product
form
\[
(\psi,\varphi)=\psi(0)\varphi(0)-\int_{-1}^0\int_0^\theta
\psi(\xi-\theta)d \eta(\theta)\varphi(\xi)d\xi,
\]
 we have
$$
u_1=\frac{y^*}{(1+i\sigma_k)(y^*-x^*v_2^2)+\tau_kax^*y^*(1-v_2^2)}.
$$
Using that the enlarged phase space
$BC:\{\varphi:[-1,0]\to \mathbb{C}^2| \varphi$
is continuous on $[-1,0)$, it exists $\lim_{\theta\to 0^-} \varphi(\theta)$.
The elements of
$BC$ have the form $\psi=\varphi+X_0\alpha$, where $\varphi\in C$,
$\alpha\in \mathbb{C}^2$, $X_0=X_0(\theta)$ is given by
\[
X_0(\theta)=\begin{cases}
I, & \theta=0,\\
0, & -1\leq \theta <0. \end{cases}
\]
The projection of $C$ onto $P$, $\varphi \mapsto
\Phi(\Psi,\varphi)$, associated with the decomposition
$C=P\oplus Q$ is replaced by $\pi :BC\mapsto P$ such that
\[
\pi(\varphi +X_0 \alpha)=\Phi[(\Psi,\varphi)+\Psi(0)\alpha].
\]
Thus, we have the decomposition $BC=P\oplus \ker\pi$.
 Using the decomposition
$z_t=\Phi x(t)+y_t$, $x(t)\in \mathbb{C}^2$,
$y_t \in \ker\pi\cap C^1=: Q^1$, the equation \eqref{3.70}
is equivalent to the system
\begin{equation}\label{3.80}
 \begin{gathered}
 \dot{x}=Bx+\Psi(0)F_0(\Phi x+y,\alpha), \\
\frac{d}{dt}y=A_{Q^1}y+(I-\pi)X_0 F_0(\Phi x+y,\alpha).
 \end{gathered}
\end{equation}
Consider the Taylor expansion as follows:
\begin{equation}\label{3.90}
\Psi(0)F_0(\Phi x+y,\alpha)=
\frac{1}{2}f_2^1(x,y,\alpha)+\frac{1}{3!}f_3^1(x,y,\alpha)+h.o.t.,
\end{equation}
\[
(I-\pi)X_0 F_0(\Phi x+y,\alpha)=
\frac{1}{2}f_2^2(x,y,\alpha)+\frac{1}{3!}f_3^2(x,y,\alpha)+h.o.t.,
\]
where $f_j^1(x,y,\alpha)$ and $f_j^2(x,y,\alpha)$ are homogeneous
polynomials in $(x,y,\alpha)$ of degree $j$, $j=2,3$, with
coefficients in $\mathbb{C}^2$ and $\ker\pi$,
 respectively, h.o.t. stands for higher order terms.
Thus, the normal form method gives a normal form on the
 center manifold of the origin for \eqref{3.80} as
\begin{equation}\label{3.100}
\dot{x}=Bx+\frac{1}{2}g_2^1(x,0,\alpha)+\frac{1}{3!}g_3^1(x,0,\alpha)+h.o.t,
\end{equation}
where $g_2^1$, $g_3^1$ are the second and third order terms in
$(x,\alpha)$, respectively. Always following \cite{Faria2}, $V_j^{m+p}(X)$
denotes the linear space of the homogeneous
polynomials of degree $j$ in $m+p$ real variables,
 $x=(x_1,\ldots,x_m),\alpha=(\alpha_1,\ldots,\alpha_p)$ with
coefficients in $X$; that is,
\[
V_j^{m+p}(X)=\Big\{ \sum_{|(q,l)|=j} c_{(q,l)}x^q\alpha^l:(q,l)\in
N_0^{m+p},c_{(q,l)}\in X\Big\}.
\]
 The operators $M_j^1$ are defined by
 \[
(M_j^1p)(x,\alpha)=D_xp(x,\alpha)Bx-Bp(x,\alpha),\quad j\geq 2,
\]
 and
$M_j^1$ act on
$V_j^3(\mathbb{C}^2)=\operatorname{Im}(M_j^1)\oplus
\operatorname{\ker}(M_j^1)$, and
\[
\ker(M_j^1)=\operatorname{span}\{x^q\alpha^le_k:
(q,\bar{\lambda})=\lambda_k, k=1,2, q\in \mathbb{N}_0^2, l\in
\mathbb{N}_0, |(q,l)|=j\},
\]
 where ${e_1,e_2}$ is the
canonical basis of $\mathbb{C}^2$. Hence,
\begin{gather*}
\ker(M_2^1)=\operatorname{span}\Big\{
 \begin{pmatrix}
 x_1\alpha \\
 0 \end{pmatrix},
 \begin{pmatrix}
 0 \\
 x_2\alpha
 \end{pmatrix}\Big\},
\\
\ker(M_3^1)=\operatorname{span}\Big\{
 \begin{pmatrix}
 x_1^2x_2 \\
 0 \end{pmatrix},
 \begin{pmatrix}
 x_1\alpha^2 \\
 0
 \end{pmatrix},
\begin{pmatrix}
 0 \\
 x_1x_2^2
 \end{pmatrix},
 \begin{pmatrix}
 0 \\
 x_2\alpha^2
 \end{pmatrix}\Big\}.
\end{gather*}
 From \eqref{3.90}, we have
 \begin{equation}\label{3.110}
f_2^1(x,y,\alpha)=\Psi(0)[2L(\alpha)(\Phi x+y)+F_2(\Phi x+y,\tau_k)],
\end{equation}
 where
\[
F_2(\varphi,\tau_k) = 2\tau_k
 \begin{pmatrix}
 -a\varphi_1^2(0)-b\varphi_1(0)\varphi_1(-1)-c\varphi_1(0)\varphi_2(-1)\\
 -a\varphi_2^2(0)+c\varphi_1(-1)\varphi_2(0)-b\varphi_2(0)\varphi_2(-1)
 \end{pmatrix},
\]
 Therefore, we  obtain
\begin{equation}\label{3.120}
f_2^1(x,0,\alpha)=
 \begin{pmatrix}
 2A_1x_1\alpha+ 2A_2x_2\alpha+a_{20}x_1^2+a_{02}x_2^2\\
 2\bar{A}_2x_1\alpha+ 2\bar{A}_1x_2\alpha+\bar{a}_{02}x_1^2+
 \bar{a}_{20}x_2^2
 \end{pmatrix},
\end{equation}
here
\begin{equation} \label{3.125}
\begin{gathered}
A_1=\frac{i\sigma_k}{\tau_k}(u_1+u_2v_2),\quad
A_2=\frac{-i\sigma_k}{\tau_k}(u_1+u_2\bar{v}_2),\\
a_{20}=2i\sigma_k(\frac{u_1}{x^*}+\frac{u_2v_2^2}{y^*}),\quad
a_{02}=-2i\sigma_k(\frac{u_1}{x^*}+\frac{u_2\overline{v}_2^2}{y^*}).
\end{gathered}
\end{equation}
 Since the second order terms
in $(\alpha,x)$ of the normal form on the center
manifold are given by
\[
\frac{1}{2}g_2^1(x,0,\alpha)=\frac{1}{2}\operatorname{Proj}_{\ker(M_2^1)}
f_2^1(x,0,\alpha),
\]
it follows that
\begin{equation}\label{3.130}
\frac{1}{2}g_2^1(x,0,\alpha)=
 \begin{pmatrix}
 A_1x_1\alpha \\
 \bar{A}_1x_2\alpha
 \end{pmatrix},
\end{equation}
where $A_1=\frac{i\sigma_k}{\tau_k}(u_1+u_2v_2)$. We notice that
\[
g_3^1(x,0,\alpha)\in \ker(M_3^1)=\operatorname{span}\Big\{
 \begin{pmatrix}
 x_1^2x_2 \\
 0
 \end{pmatrix},
 \begin{pmatrix}
 x_1\alpha^2 \\
 0
 \end{pmatrix},
  \begin{pmatrix}
 0 \\
 x_1x_2^2
 \end{pmatrix},
  \begin{pmatrix}
 0 \\
 x_2\alpha^2
 \end{pmatrix}\Big\}.
\]
However, the terms $O(|x|\alpha^2)$ are irrelevant to determine the
generic Hopf bifurcation. Thus, we only need to compute the coefficients of
$ \begin{pmatrix}
 x_1^2x_2 \\
 0  \end{pmatrix}$
 and $ \begin{pmatrix}
 0 \\
 x_1x_2^2  \end{pmatrix}$.
Note that
\[
\frac{1}{3!}g_3^1(x,0,\alpha)
=\frac{1}{3!}\operatorname{Proj}_{\ker(M_3^1)}\tilde{f}_3^1(x,0,\alpha)
=\frac{1}{3!}\operatorname{Proj}_s\tilde{f}_3^1(x,0,0)+O(|x|\alpha^2),
\]
where
\[
s:=\operatorname{span}\Big\{
 \begin{pmatrix}
 x_1^2x_2 \\
 0  \end{pmatrix},
 \begin{pmatrix}
 0 \\
 x_1x_2^2
 \end{pmatrix}\Big\}
 \]
and
\[
\tilde{f}_3^1(x,0,0)=f_3^1(x,0,0)+
\frac{3}{2}[(D_xf_2^1)U_2^1-(D_xU_2^1)g_2^1]_{(x,0,0)}
+\frac{3}{2}[(D_yf_2^1)h]_{(x,0,0)},
\]
here $\tilde{f}_3^1(x,0,0)$ is the third order terms of the equation
which is obtained after computing
the second terms of the normal form.

Now we compute $\frac{1}{3!}g_3^1(x,0,\alpha)$ step by step.

We first compute $\operatorname{Proj}_s[(D_xf_2^1)U_2^1]_{(x,0,0)}$.
Following \cite{Faria2}, we know that
\[
U_2^1(x,0)=(M_2^1)^{-1}P_{I,2}^1f_2^1(x,0,0).
\]
From \eqref{3.120}, we can easily see that
\[
f_2^1(x,0,0)= \begin{pmatrix}
 a_{20}x_1^2+a_{02}x_2^2 \\
 \bar{a}_{02}x_1^2+\bar{a}_{20}x_2^2
 \end{pmatrix}.
\]
According to the definition of $M_2^1$, the equation
$M_2^1U_2^1(x,0)=f_2^1(x,0,0)$ can be written as the following
differential equations:
\begin{gather*}
x_1\frac{\partial u_1}{\partial x_1}-x_2\frac{\partial
u_1}{\partial
x_2}-u_1=\frac{1}{i\sigma_k}(a_{20}x_1^2+a_{02}x_2^2),\\
x_1\frac{\partial u_2}{\partial x_1}-x_2\frac{\partial u_2}{\partial
x_2}+u_2=\frac{1}{i\sigma_k}(\bar{a}_{02}x_1^2+\bar{a}_{20}x_2^2).
\end{gather*}
For the equation, we can easily obtain
\[
U_2^1(x,0)= \begin{pmatrix}
 \frac{1}{i\sigma_k}(a_{20}x_1^2-\frac{1}{3}a_{02}x_2^2) \\
 \frac{1}{i\sigma_k}(\frac{1}{3}\bar{a}_{02}x_1^2-\bar{a}_{20}x_2^2)
 \end{pmatrix}.
\]
Thus,
\[
\operatorname{Proj}_s[(D_xf_2^1)U_2^1]_{(x,0,0)}
= \begin{pmatrix}
 \frac{2}{3i\sigma_k}|a_{02}|^2x_1^2x_2 \\
 \frac{-2}{3i\sigma_k}|a_{02}|^2x_1x_2^2
 \end{pmatrix}.
\]
From \eqref{3.130}, we know $g_2^1(x,0,0)=0$. Thus,
 $[(D_xU_2^1)g_2^1]_{(x,0,0)}=0$.

Next, we compute
$\operatorname{Proj}_s[(D_yf_2^1)h]_{(x,0,0)}$, where
$h=(h^1,h^2)^T$ is a second order homogeneous polynomial in
$(x_1,x_2,\alpha)$ with coefficients in $Q^1$, which has the form
\[
h=h(x_1,x_2,\alpha)
=h_{110}x_1x_2+h_{101}x_1\alpha+h_{011}x_2\alpha+h_{200}x_1^2+
h_{020}x_2^2+h_{002}\alpha^2
\]
and $h=h(x_1,x_2,\alpha)$ is the unique solution in $V_2^3(Q^1)$ of
the equation
\[
(M_2^2h)(x,\alpha)=(I-\pi)X_0[2L(\alpha)(\Phi x)+F_2(\Phi
x,\tau_k)].
\]
 As in \cite{Faria2}, we know that
\begin{align*}
(M_2^2h)(x,\alpha)
&=D_xh(x,\alpha)Bx-A_{Q^1}(h(x,\alpha))\\
&=D_xh(x,\alpha)Bx-\dot{h}(x,\alpha)-X_0[L(\tau_k)(h(x,\alpha))
 -\dot{h}(x,\alpha)(0)]\\
&=(I-\pi)X_0[2L(\alpha)(\Phi x)+F_2(\Phi x,\tau_k)].
\end{align*}
thus, $h(x,\theta)$ is also evaluated by the system
\begin{gather}\label{3.140}
\dot{h}(x)-D_xh(x)Bx=\Phi \Psi(0)[F_2(\Phi x,\tau_k)],\\
\label{3.150}
\dot{h}(x)(0)-L(\tau_k)(h(x))=F_2(\Phi x,\tau_k),
\end{gather}
where $\dot{h}$ denotes the derivative of $h(x)(\theta)$ relative to
$\theta$. According to \eqref{3.110}, it follows that
\begin{align*}
f_2^1(x,y,0)
&=\Psi(0)F_2(\Phi x+y,\tau_k)\\
&= \begin{pmatrix}
 u_1 &u_2\\
 \bar{u}_1 & \bar{u}_2
 \end{pmatrix}
 \tau_k
 \begin{pmatrix}
 -ad_1^2-bd_1d_2-cd_1n_2 \\
 -an_1^2+cd_2n_1-bn_1n_2
 \end{pmatrix},
\end{align*}
where
\begin{gather*}
d_1=x_1+x_2+y_1(0),\quad
d_2=e^{-i\sigma_k}x_1+e^{i\sigma_k}x_2+y_1(-1),\\
n_1=v_2x_1+\overline{v}_2x_2+y_2(0),\quad
l_2=e^{-i\sigma_k}v_2x_1+e^{i\sigma_k}\overline{v}_2x_2+y_2(-1).
\end{gather*}
Then we have
\begin{align*}
&[(D_yf_2^1)h]_{(x,0,0)}\\
&= 2\tau_k
 \begin{pmatrix}
 u^T \begin{pmatrix}
 -2ad_1'h^1(0)-bd_1'h^1(-1)-bd_2'h^1(0)-cd_1'h^2(-1)-cn_2'h^1(0) \\
 -2an_1'h^2(0)+cd_2'h^2(0)+cn_1'h^1(-1)-bn_1'h^2(-1)-bn_2'h^2(0)
 \end{pmatrix}
  \\
 \bar{u}^T \begin{pmatrix}
 -2ad_1'h^1(0)-bd_1'h^1(-1)-bd_2'h^1(0)-cd_1'h^2(-1)-cn_2'h^1(0) \\
 -2an_1'h^2(0)+cd_2'h^2(0)+cn_1'h^1(-1)-bn_1'h^2(-1)-bn_2'h^2(0)
 \end{pmatrix}
  \\
 \end{pmatrix},
 \end{align*}
where
\begin{equation}  \label{3.160}
\begin{gathered}
d_1'=x_1+x_2,\quad d_2'=e^{-i\sigma_k}x_1+e^{i\sigma_k}x_2,\\
n_1'=v_2x_1+\overline{v}_2x_2,\quad
n_2'=e^{-i\sigma_k}v_2x_1+e^{i\sigma_k}\overline{v}_2x_2.
\end{gathered}
\end{equation}
Thus,
\[
\operatorname{Proj}_s[(D_yf_2^1)h)]_{(x,0,0)}
= \begin{pmatrix}
 2c_3x_1^2x_2 \\
 2\bar{c}_3x_1x_2^2 \\
 \end{pmatrix},
\]
where
\begin{align*}
&c_3=u^T\tau_k \\
&\times\begin{pmatrix}
 Ah_{110}^1(0)+\overline{A}h_{200}^1(0)+Bh_{110}^1(-1)
 +\overline{B}h_{200}^1(-1)+Ch_{110}^2(-1)+\overline{C}h_{200}^2(-1) \\
 Dh_{110}^1(-1)+\overline{D}h_{200}^1(-1)+Eh_{110}^2(0)
 +\overline{E}h_{200}^2(0)
 +Fh_{110}^2(-1)+\overline{F}h_{200}^2(-1)
 \end{pmatrix}
\end{align*}
and
\begin{gather*}
A=-2a-be^{-i\sigma_k}-ce^{-i\sigma_k}v_2,\quad
B=-b, \quad C=-c,\\
D=cv_2,\quad
E=-2av_2+ce^{-i\sigma_k}-be^{-i\sigma_k}v_2, \quad F=-bv_2.
\end{gather*}
Thus, we only need to compute $h_{110}(\theta)$ and
$h_{200}(\theta)$. From \eqref{3.140} and \eqref{3.150}, we know
that $h_{110}=(h_{110}^1,h_{110}^2)^T$ and
$h_{200}=(h_{200}^1,h_{200}^2)^T$ are respectively the solution of
the following two systems
\begin{equation}\label{3.170}
\begin{gathered}
\dot{h}_{110}= \begin{pmatrix}
 0 \\
 0
 \end{pmatrix},
\\
\dot{h}_{110}(0)-L(\tau_k)(h_{110})=2\tau_k
 \begin{pmatrix}
 a_1 \\
 b_1
 \end{pmatrix},
\end{gathered}
\end{equation}
and
\begin{equation}\label{3.180}
\begin{gathered}
\dot{h}_{200}-2i\sigma_kh_{200}=\begin{pmatrix}\Phi_1  \Phi_2\end{pmatrix}
 \begin{pmatrix}
 a_{20} \\
 \bar{a}_{02}
 \end{pmatrix},\\
\dot{h}_{200}(0)-L(\tau_k)(h_{200})=2\tau_k
 \begin{pmatrix}
 a_2 \\
 b_2
 \end{pmatrix},
\end{gathered}
\end{equation}
where
\begin{gather*}
a_1=-2a-be^{i\sigma_k}-be^{-i\sigma_k}-ce^{i\sigma_k}\bar{v}_2
 -ce^{-i\sigma_k}v_2=0,\\
b_1=-2av_2\bar{v}_2+ce^{-i\sigma_k}\bar{v}_2+ce^{i\sigma_k}v_2
 -be^{i\sigma_k}v_2\bar{v}_2
-be^{-i\sigma_k}v_2\bar{v}_2=0,\\
a_2=-a-be^{-i\sigma_k}-cv_2e^{-i\sigma_k},\\
b_2=-a v_2^2+ce^{-i\sigma_k}v_2-bv_2^2e^{-i\sigma_k}.
\end{gather*}
Solving \eqref{3.170} and \eqref{3.180}, we obtain
\begin{equation}\label{3.200}
h_{110}(\theta)=c_1,
\end{equation}
and
\begin{equation}\label{3.210}
h_{200}(\theta)=-\frac{1}{i\sigma_k}(a_{20}e^{i\sigma_k\theta}
v+\frac{1}{3}\bar{a}_{02}e^{-i\sigma_k\theta}\bar{v})+e^{2i\sigma_k\theta}c_2,
\end{equation}
where $c_1=(c_1^{(1)},c_1^{(2)})^T$, $c_2=(c_2^{(1)},c_2^{(2)})^T$
and
\begin{gather*}
c_1^{(1)}=\frac{2a_1(a+b)y^*-2b_1cx^*}{x^*y^*[(a+b)^2+c^2]},\quad
c_1^{(2)}=\frac{2a_1cy^*+2b_1(a+b)x^*}{x^*y^*[(a+b)^2+c^2]},
\\
\begin{aligned}
c_2^{(1)}&=\Big(2\tau_ka_2(2i\sigma_k+\tau_kay^*+\tau_kby^*e^{-2i\sigma_k})
-2\tau_k^2b_2cx^*e^{-2i\sigma_k}\Big)\\
&\quad\div \Big(
\tau_k^2c^2x^*y^*e^{-4i\sigma_k}+(2i\sigma_k+\tau_kax^*\\
&\quad +\tau_kbx^*e^{-2i\sigma_k})(2i\sigma_k+\tau_kay^*
 +\tau_kby^*e^{-2i\sigma_k})\Big),
\end{aligned}
\\
\begin{aligned}
c_2^{(2)}&=\Big(2\tau_kb_2(2i\sigma_k+\tau_kax^*
+\tau_kbx^*e^{-2i\sigma_k})+2\tau_k^2a_2cy^*e^{-2i\sigma_k}\Big)\\
&\quad\div
\Big(\tau_k^2c^2x^*y^*e^{-4i\sigma_k}+(2i\sigma_k+\tau_kax^*\\
&\quad +\tau_kbx^*e^{-2i\sigma_k})(2i\sigma_k+\tau_kay^*
 +\tau_kby^*e^{-2i\sigma_k})\Big).
\end{aligned}
\end{gather*}
Finally, we compute $\operatorname{Proj}_sf_3^1(x,0,0)$. From
 \eqref{3.90}, we can see that
$f_3^1(x,0,0)=0$.

Summarizing the above steps, we obtain
\[
\frac{1}{3!}g_3^1(x,0,0)=
 \begin{pmatrix}
 A_3x_1^2x_2 \\
 \bar{A}_3x_1x_2^2
 \end{pmatrix},
\]
where
\begin{equation}\label{3.215}
A_3=\frac{1}{6i\sigma_k}|a_{02}|^2+\frac{1}{2}c_3.
\end{equation}
Accordingly, the normal form of \eqref{3.100} has the form
\begin{align*}
\dot{x}
&= Bx+\frac{1}{2}g_2^1(x,0,\alpha)+\frac{1}{3!}g_3^1(x,0,\alpha)+h.o.t.\\
&= Bx+\begin{pmatrix}
 A_1x_1\alpha \\
 \bar{A}_1x_2\alpha
 \end{pmatrix}
  +  \begin{pmatrix}
 A_3x_1^2x_2 \\
 \bar{A}_3x_1x_2^2
 \end{pmatrix}
 +O(|x|\alpha^2+|x|^4).
\end{align*}
The normal form \eqref{3.100} relative to $P$ can be written in real
coordinates $(w_1,w_2)$ through the change of variables
$x_1=w_1-iw_2$, $x_2=w_1+iw_2$. Followed by the use of polar
coordinates $(\rho,\xi)$, $w_1=\rho\cos\xi$, $w_2=\rho\sin\xi$,
 this formal form becomes
\begin{equation}\label{3.220}
 \begin{gathered}
\dot{\rho}=k_1\alpha\rho+k_2\rho^3+O(\alpha^2\rho+|(\rho,\alpha)|^4), \\
\dot{\xi}=-\sigma_k+O(|(\rho,\alpha|),
 \end{gathered}
\end{equation}
where $k_1=\operatorname{Re}A_1,\quad k_2=\operatorname{Re}A_3$.
Following \cite{Chow}, we know that the sign of $k_1k_2$ determines the
direction of the bifurcation and the sign of $k_2$ determines the
stability of the nontrivial periodic solution bifurcating from Hopf
bifurcation. Therefore, we have the following theorem

\begin{theorem}\label{thm 3.10}
The flow of \eqref{3.70} on the center manifold of the origin at
$\tau=\tau_k$ is given by \eqref{3.220}. Hopf bifurcation is
supercritical if $k_1k_2<0$ and subcritical if $k_1k_2>0$. Moreover,
the nontrivial periodic solution is stable if $k_2<0$ and unstable
if $k_2>0$.
\end{theorem}

We remark that by the same method and calculation,
for the nonsymmetric system, we can also obtain a similar conclusion.
Thus, the symmetry does not affect the emergence of periodic solutions.


\section{Numerical simulations}

 In the section, we will give numerical simulations to illustrate
the conclusions obtained. Consider system \eqref{90} with
$r_1=3/2$, $r_2=1$, $a=1/2$,
$b=1/3$, $c=1$; that is,
\begin{equation}\label{4.10}
\begin{gathered}
\dot{x}(t)=x(t)[\frac{3}{2}-\frac{1}{2}x(t)-\frac{1}{3}x(t-\tau)-y(t-\tau)],\\
\dot{y}(t)=y(t)[1-\frac{1}{2}y(t)+x(t-\tau)-\frac{1}{3}y(t-\tau)].
 \end{gathered}
\end{equation}
In this case, system \eqref{4.10} has unique positive
 equilibrium $(x^*,y^*)=(\frac{9}{61},\frac{84}{61})$.
 By means of the software
 Maple, we know that $y_0=0.0599$ is a simple root of \eqref{150}.
By computation, we obtain
 $\sigma_+=0.2448, \tau_0=4.4898, \sigma_0=\tau_0\sigma_+=1.0989$.
 Applying Theorems \ref{th60} and \ref{thm 3.10}, we have the following theorem.

\begin{theorem}
The characteristic equation of \eqref{4.10} at
$(x^*,y^*)=(\frac{9}{61},\frac{84}{91})$ has a simple pair of purely
imaginary roots $\pm i\sigma_0$ for $\tau=\tau_0$ and the other
roots has non-zero real parts. System \eqref{4.10} with
$\tau_0=4.4898$ has a supercritical Hopf bifurcation and the
nontrivial periodic solution bifurcating from Hopf bifurcation of
\eqref{4.10} is stable on the center manifold.
\end{theorem}

\begin{proof}
 From Theorem \ref{th60}, we know that the equilibrium
$(x^*,y^*)$ is locally asymptotically stable when
$\tau \in[0,4.4898)$ and is unstable when $\tau >4.4898$.
In the following, we will consider the direction
and stability of Hopf bifurcation at $\tau =\tau_0=4.4898$.
By means of software Maple, we  obtain the following values:
$$
v= \begin{pmatrix}
 v_1 \\
 v_2  \end{pmatrix}
=\begin{pmatrix}
 1 \\
 0.9197-1.1994i  \end{pmatrix},
\quad u= \begin{pmatrix}
 u_1 \\
 u_2  \end{pmatrix}
 = \begin{pmatrix}
 0.2107-0.3369i\\
 0.0226+0.0602i
 \end{pmatrix},
$$
and
\begin{gather*}
a_{20} = 5.1560+3.3283i,\quad
 a_{02}=-4.9973-2.9057i,\\
 c_3 = -24.0701-5.0921i.
\end{gather*}
According to \eqref{3.125} and \eqref{3.215}, we obtain
\[
A_1=0.0756+0.0744i,\quad A_3=-12.0351-7.6141i.
\]
Thus, $k_1=0.0756$, $k_2=-12.0351$. Applying Theorem \ref{thm 3.10},
we obtain that system \eqref{4.10}
has a supercritical Hopf bifurcation at $\tau=\tau_0$ and the
nontrivial periodic solution
associated with Hopf bifurcation at $\tau=\tau_0$ is stable on the
center manifold. The proof is complete.
\end{proof}

The numerical simulations for $\tau=4.4898$ are shown in
Figure \ref{fig1}.

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.48\textwidth]{fig1a} %xiaa1.eps
\includegraphics[width=0.48\textwidth]{fig1b}\\ %xiab1.eps
\includegraphics[width=0.48\textwidth]{fig1c} % xiaa2.eps
\includegraphics[width=0.48\textwidth]{fig1d} % xiab2.eps
\includegraphics[width=0.48\textwidth]{fig1e} % xiaa3.eps
\includegraphics[width=0.48\textwidth]{fig1f} % xiab3.eps
\end{center}
\caption{Trajectories for system \eqref{4.10} with
$\tau=4.41$ (left column) and $\tau=4.49$ (right column),
on the t-x plane (top row), t-y plane (middle row), and
x-y plane (bottom row)}
\label{fig1}
\end{figure}


\subsection*{Acknowledgements}
 The authors would like to thank the anonymous referees for their
carefully reading of the original manuscript, and for their valuable
comments and suggestions for improving the results as well as the
exposition of this article.

Z. Yu was supported by grant 11101282 from the National Natural
Science Foundation of China.
R. Yuan was supported by National Natural
Science Foundation of China and RFDP.


\begin{thebibliography}{99}

\bibitem{Beretta} E. Beretta, Y. Kuang;
Convergence results in a well-known delayed
predator-prey systems, \emph{J. Math. Anal. Appl.}, \textbf{204} (1996),
840-853.

\bibitem{Chow} S. N. Chow, J. K. Hale;
Methods of Bifurcation Theory, Springer-Verlag, New York, 1982.

\bibitem{Cooke} K. L. Cooke, Z. Grossman;
Discrete delay, distributed delay and
stability switches, \emph{J. Math. Anal. Appl.}, \textbf{86} (1982),
592-627.

\bibitem{Faria1} T. Faria, L. T. Magalh\~{a}es;
Normal forms for retarded functional differential equations and
applications to Bogdanov-Takens singularity,
\emph{J. Differential Equations}, \textbf{122} (1995), 201-224.

\bibitem{Faria2} T. Faria, L. T. Magalh\~{a}es;
Normal forms for retarded functional differential equations with
parameters and applications to Hopf bifurcation,
\emph{J. Differential Equations}, \textbf{122} (1995), 181-20.

\bibitem{Faria3} T. Faria;
Stability and bifurcation for a delayed predator-prey
model and the effect of diffusion,
 \emph{J. Math. Anal. Appl.}, \textbf{254} (2001), 433-463.

\bibitem{Gopalsamy} K. Gopalsamy;
 \emph{Stability and oscillation in delayed differential
equations of population dynamics}, Dordrecht:Kluwer Academic Press,
1992.

\bibitem{jh77} J. K. Hale;
\emph{Theory of Functional Differential Equations},
Springer-Verlag, New York, 1977.

\bibitem{Kuang} Y. Kuang;
\emph{Delay Differential Equations with Application in
Population Dynamics}, Academic Press, New York, 1993.

\bibitem{Liu} Z. H. Liu, R. Yuan;
Stability and bifurcation in a harmonic
oscillator with delays, \emph{Chaos, Solitons and Fractals}, \textbf{23}
(2005), 551-562.

\bibitem{May} R. M. May;
Time delay versus stability in population models with two
and three trophic levels, \emph{Ecology}, \textbf{4} (1973), 315-325.

\bibitem{Nakao} S. Nakaoka, Y. Saito, Y. Takeuchi;
Stability, delay, and chaotic
behavior in a Lotka-Volterra predator-prey system, \emph{Math.
Biosci. Engineer.}, \textbf{3} (2006), 173-187.


\bibitem{Ruan} S. G. Ruan;
 Absolute stability, conditional stability and
bifurcation in Kolmogorov-type predator-prey systems with discrete
delays, \emph{Quarterly of Applied Mathematics}, \textbf{59} (2001),
159-173.

\bibitem{Saito} Y. Saito;
Permanence and global stability for general
Lotka-Volterra predator-prey systems with distributed delays,
\emph{Nonlinear Anal.,}, \textbf{47} (2001), 6157-6168.


\bibitem{song} Y. L. Song, J. J. Wei;
Local Hopf bifurcation and global
periodic solutions in a delayed predarot-prey systems,
\emph{J. Math. Anal. Appl.}, \textbf{301} (2005), 1-21.



\bibitem{sunhan} C. J. Sun, M. A. Han, Y. P. Yang;
Analysis of stability and Hopf bifurcation for a delayed logistic equation,
\emph{Chaos, Solitons and Fractals}, \textbf{31} (2007), 672-682.

\bibitem{Sun} C. Sun, Y. Lin, M. Han ;
Stability and Hopf bifurcation for an epidemic model with delay,
\emph{Chaos, Solitons and Fractals}, \textbf{30} (2006), 204-216.

\bibitem{Volterra} V. Volterra;
\emph{Lecons sur la th\'{e}orie mathematique de la lutte
pour la vie}, Gauthier-Villars, Paris, 1931.

\bibitem{Wei} J. J. Wei, S. G. Ruan;
 Stability and bifurcation in a neural network model with two delays,
\emph{Physica D}, \textbf{130} (1990), 255-272.

\bibitem{yanLi}  X. P. Yan, W.T. Li;
 Hopf bifurcation and global periodic
solutions in a delayed predator-prey system, \emph{Applied
Mathematics and Computation}, \textbf{177} (2006), 427-445.

\bibitem{YanZhang} X. P. Yan, C. H. Zhang;
 Hopf bifurcation in a delayed Lotka-Volterra predator-prey system,
 \emph{Nonlinear Anal.: Real World Appl.}, \textbf{9} (2008), 114-127.

\end{thebibliography}

\end{document}