\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2013 (2013), No. 198, pp. 1--12.\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/198\hfil Stability and bifurcation analysis]
{Stability and bifurcation analysis for a discrete-time bidirectional
ring neural network model with delay}

\author[Y.-K. Du, R. Xu, Q.-M. Liu \hfil EJDE-2013/198\hfilneg]
{Yan-Ke Du, Rui Xu, Qi-Ming Liu} 

\address{Yan-Ke Du \newline
Institute of Applied Mathematics,
Shijiazhuang Mechanical Engineering College,
Shijiazhuang. 050003, China}
\email{yankedu2011@163.com}

\address{Rui Xu \newline
Institute of Applied Mathematics,
Shijiazhuang Mechanical Engineering College,
Shijiazhuang. 050003, China}
\email{rxu88@163.com}

\address{Qi-Ming Liu \newline
Institute of Applied Mathematics,
Shijiazhuang Mechanical Engineering College,
Shijiazhuang. 050003, China}
\email{lqmmath@163.com}

\thanks{Submitted January 3, 2012. Published September 5, 2013.}
\subjclass[2000]{92B20, 34K18, 34K20, 37G05}
\keywords{Neural network; time delay; stability; bifurcation}

\begin{abstract}
 We study a class of discrete-time bidirectional ring neural
 network model with delay. We discuss the asymptotic stability 
 of the origin  and the existence of Neimark-Sacker bifurcations, 
 by analyzing the corresponding characteristic equation.
 Employing M-matrix theory and the Lyapunov functional method,
 global asymptotic stability of the origin is derived.
 Applying the normal form theory and the center manifold theorem,
 the direction of the Neimark-Sacker bifurcation and the stability
 of bifurcating periodic solutions are obtained. Numerical simulations
 are given to illustrate the main results.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{example}[theorem]{Example}
\newtheorem{remark}[theorem]{Remark}
\allowdisplaybreaks

\section{Introduction}

 Since Hopfiled's pioneering work
 \cite{s3,s4}, the dynamic  behavior (including stability,
periodic oscillatory and chaos) of continuous-time Hopfield
neural networks has received much attention
due to their applications in optimization, signal processing, image
processing, solving nonlinear algebraic equation, pattern
recognition, associative memories and so on (see, \cite{s9,s10,s13,s6} and
references therein).

It is well known that time delays in the information
processing of neurons exist. The delayed axonal signal transmissions in
the neural networks make the dynamic behaviors more
complicated, and may destabilize stable equilibria and give rise
to periodic oscillation, bifurcation and chaos (see
\cite{s11,s12,s13,s14}). Therefore, the delay is inevitable and cannot
be neglected.
For computer simulations, experimental or computational purposes, it
is common to discretize the continuous-time neural networks. In some
sense, the discrete-time model inherits the dynamical
characteristics of the continuous-time networks. We refer the reader to
\cite{s15,s8,s16,s7} for related discussions on the need and
importance of discrete-time analogues to reflect the dynamics of
their continuous-time counterparts.

In the field of neural networks, rings are studied to gain insight
into the mechanisms underlying the behavior of recurrent  networks.
 In \cite{s5}, Wang and Han investigated the following continuous-time 
bidirectional ring network model with delay
 \begin{equation}\label{a1}
\begin{gathered}
\dot{x}=-x+\alpha f(y(t-\tau))+\beta f(z(t-\tau)),\\
\dot{y}=-y+\alpha f(z(t-\tau))+\beta f(x(t-\tau)),\\
\dot{z}=-z+\alpha f(x(t-\tau))+\beta f(y(t-\tau)),
\end{gathered}
\end{equation}
where $\tau$ denotes the synaptic transmission delay, $\alpha$ and
$\beta$ are connection strengths, $f:$ $R\to R$ is the
activation function. In \cite{s5}, some conditions on the linear
stability of the trivial solution of system \eqref{a1} were given
and Hopf bifurcation, including its direction and stability, were
investigated.

Motivated by the work of Wang and Han \cite{s5} and the discussions
above, in the present paper, for simplicity, assuming the neurons in
the network be identical (see \cite{s17}), we are concerned
 with the stability and bifurcation analysis of the following discrete-time
bidirectional ring neural network model with delay
 \begin{equation}\label{a2}
\begin{gathered}
x(n+1)=ax(n)+\beta f(y(n-k))+\beta f(z(n-k)),\\
y(n+1)=ay(n)+\beta f(z(n-k))+\beta f(x(n-k)),\\
z(n+1)=az(n)+\beta f(x(n-k))+\beta f(y(n-k)),\\
\end{gathered}
\end{equation}
where $a\in(0,1)$ is the internal decay of the neurons, $\beta$
is the connection strength, $k\in N$ is the time delay.

This paper  contributes to understanding of neural networks as follows:

(1) There is a  large body of work discussing the stability and bifurcation 
of neural networks with delays, but most of them deal only with  
continuous-time neural network models,
or  discrete-time neural network models of two neurons with or without time 
delays (\cite{s8,s7}).
Here  we discuss the dynamic  behavior of a tri-neuron discrete-time 
bidirectional ring neural network with delay.
The characteristic equation of the neural network is a polynomial equation
 with high order terms. Using a new  approach, sufficient and necessary 
conditions are derived to ensure that all the roots of the characteristic 
equation stay inside or on the unit circle.

(2) We remove some restrictions on the conditions required by \cite{s5}, 
and in a sense, our results on the  asymptotic stability of the origin  
are less restrictive than those for the corresponding continuous system in
\cite{s5}.

(3) Employing M-matrix theory and the Lyapunov functional
method, global asymptotic stability of
the origin is derived, which was not taken into account in \cite{s5}.
 The stability criterion is simple and can be easily checked.

The rest of this paper is organized as follows. In Section 2, we
analyze the location of roots of a class of polynomial equation. In
Section 3, the local stability of the origin and the existence of
Neimark-Sacker bifurcations are discussed by analyzing the
corresponding characteristic equation, and global asymptotic
stability is derived using the method of M-matrix and Lyapunov
function. In Section 4, we discuss the stability and direction of
the Neimark-Sacker bifurcation by employing the normal form method and
the center manifold theorem. Some numerical simulations are carried
out in Section 5 to illustrate the main results. In Section 6, a
brief discussion is given to conclude the work of this paper.

\section{analysis of polynomial equations}

 In this section, we analyze the location of the roots
of the polynomial equation 
\begin{equation}\label{c1}
\lambda^{k+1}-a\lambda^k-b=0,~~a\in (0,1),~b\in R,
\end{equation}
which  will be used to determine the  asymptotic  stability of
system \eqref{a2}.

 Suppose that $\lambda=\mathrm{e}^{i\theta}$ is a root
of \eqref{c1}. Substituting it into \eqref{c1} and separating the
real and imaginary parts, we have
\begin{equation}\label{c2}
\begin{gathered}
\cos ((k+1)\theta)-a\cos (k\theta)=b,\\
\sin ((k+1)\theta)-a\sin (k\theta)=0.
\end{gathered}
\end{equation}
Using the identities
$\cos ((k+1)\theta)=\cos (k\theta)\cos \theta
-\sin (k\theta)\sin \theta$ and
$ \sin ((k+1)\theta)=\sin (k\theta)\cos \theta
+\cos (k\theta)\sin \theta$,
we rewrite \eqref{c2} as
\begin{gather*}
\sqrt{a^2+1-2a\cos \theta}
\big[
\frac{\cos \theta-a}{\sqrt{a^2+1-2a\cos \theta}}\cos (k\theta)
-\frac{\sin \theta}{\sqrt{a^2+1-2a\cos \theta}}\sin (k\theta)
\big]=b,\\
\sqrt{a^2+1-2a\cos \theta}\big[
\frac{\cos \theta-a}{\sqrt{a^2+1-2a\cos \theta}}\sin (k\theta)
+\frac{\sin \theta}{\sqrt{a^2+1-2a\cos \theta}}\cos (k\theta)
\big]=0.
\end{gather*}
It is easy to see that
if $\theta\in (0,\pi)$, \eqref{c2} is equivalent to the
 equations
\begin{equation}\label{c3}
\begin{gathered}
\sqrt{a^2+1-2a\cos \theta}\cdot\cos (h(\theta))=b,\\
\sin (h(\theta))=0,
\end{gathered}
\end{equation}
where
$$
h(\theta)=\operatorname{arccot}\frac{\cos \theta-a}{\sin \theta}+k\theta.
$$
Since
$$
h'(\theta)=\frac{1}{1+\big(\frac{\cos \theta-a}{\sin \theta}\big)^2}
\cdot \frac{1-a\cos \theta}{\sin ^2\theta}+k>0
$$
for $\theta\in(0,\pi)$ and
$$
\lim_{\theta\to 0^+}h(\theta)=0,\quad
\lim_{\theta\to \pi^-}h(\theta)=(k+1)\pi,
$$
we derive that $h(\theta):(0,\pi)\to(0,(k+1)\pi)$ is an
increasing bijective function.
From the second equation in \eqref{c3}, we know that 
$h(\theta)=j\pi$, $j=1,2,\dots, k$. Denote
$\theta_j=h^{-1}(j\pi)$, $j=1,2,\dots, k$. Then $\theta_j$ satisfies
the  equation
$$
j\pi=\operatorname{arccot}\frac{\cos \theta-a}{\sin \theta}+k\theta,
$$
which yields $f(\theta)=0$, where
$$
f(\theta)=\sin ((k+1)\theta)-a\sin (k\theta).
$$
Obviously,
$$
f(0)=0,~f'(0^+)=(1-a)k+1>0,\quad 
f(\frac{j\pi}{k+1})=\begin{cases}
a\sin \frac{j\pi}{k+1}>0,&\text{if $j$ is even}\\
-a\sin \frac{j\pi}{k+1}<0, &\text{if $j$ is odd}
\end{cases}
$$
$j=1,2,\dots,k$.
Therefore, we can deduce that
$\theta_j\in(\frac{(j-1)\pi}{k+1},\frac{j\pi}{k+1})$,
$j=1,2,\dots,k$.
From the first equation in \eqref{c3}, we get that 
\begin{equation}\label{c4}
b=b_j=(-1)^j\sqrt{a^2+1-2a\cos \theta_j},~j=1,2,\dots, k.
\end{equation}
If $\theta=0$, then $b=b_0=1-a>0$; if $\theta=\pi$,
then $b=b_{k+1}=(-1)^{k+1}(1+a)$.

 Obviously, if $\theta$ is a root of \eqref{c2}, $-\theta$ is also
a root of \eqref{c2}. Hence, we only need to consider the roots
$\lambda=\mathrm{e}^{i\theta}$ of \eqref{c1} in $[0,\pi]$.
Further, from \eqref{c4}, we
deduce that 
\begin{equation}\label{c5}
\dots<b_3<b_1<b_0<b_2<b_4<\dots.
\end{equation}

On the other hand, from \eqref{c1}, we have
$$
\frac{\mathrm{d}|\lambda|^2}{\mathrm{d}b}
=\frac{\mathrm{d}(\lambda\bar{\lambda})}{\mathrm{d}b}
=\lambda\frac{\mathrm{d}\bar{\lambda}}{\mathrm{d}b}+\bar{\lambda}\frac{\mathrm{d}\lambda}{\mathrm{d}b}
=\frac{2\operatorname{Re}(\lambda P(\lambda))}{|P(\lambda)|^2},
$$ 
where $P(\lambda)=\lambda^{k-1}[(k+1)\lambda-ak]$.
It follows that
$$
\frac{\mathrm{d}|\lambda|^2}{\mathrm{d}b}\big|_{b=b_j}
=\frac{2b_j}{|P(\mathrm{e}^{i\theta_j})|^2}
\Big(\frac{1-a\cos \theta_j}{a^2+1-2a\cos \theta_j}+k\Big).
$$
Hence, we have
$$
\operatorname{sign}\big\{\frac{\mathrm{d}|\lambda|^2}{\mathrm{d}b}\big|_{b=b_j}
\big\}=
\operatorname{sign}\{b_j\}=(-1)^j,\quad
j=0,1,\dots,k+1.
$$
If $b=0$, then system \eqref{c1} has two roots $a$ and $0$, which are
 inside the unit circle. As the parameter $b$ varies,
 there exist roots of \eqref{c1} which appear on or cross the unit circle.

On the other hand, it can be easily verified that both $\pi/2$ and
$2\pi/3$ are not roots of the second equation of \eqref{c2}, which
implies that $\mathrm{e}^{\pm is\theta_j}\neq 1$ for
$s=1,2,3,4$ and $j=1,2,\dots,k$.

From the discussions above, we have the following results.

\begin{theorem} \label{thm2.1}
 For equation \eqref{c1} with $a\in (0,1)$ and
$b\in R$, we have
\begin{itemize}
\item[(i)]  if $b\in(b_1,b_0)$, all the roots of \eqref{c1} are
inside the unit circle;

\item[(ii)]  if $b=b_0=1-a$, there is a simple root $\lambda=1$  of
\eqref{c1} on the unit circle, and all the other roots are inside the
unit circle;

\item[(iii) ] if $b=b_j$, $j=1,2,\dots,k$, there is a pair of complex
roots $\mathrm{e}^{\pm i\theta_j}$  on the unit
circle. Further, $\mathrm{e}^{\pm is\theta_j}\neq 1$ 
for $s=1,2,3,4$. Moreover, for the case of $b=b_1$, all other roots
are inside the unit circle;

\item[(iv)] if $b=b_{k+1}=(-1)^{k+1}(1+a)$, \eqref{c1} has a
simple root $\lambda=-1$ on the unit circle;

\item[(v)]  if $|b|>|b_{k+1}|=1+a$, all roots of \eqref{c1} are
outside the unit circle;

\item[(vi)] $\frac{\mathrm{d}|\lambda|^2}{\mathrm{d}b}\big|_{b=b_j}\neq 0$
 for $j=0,1,\dots,k+1$,
where $b_j=(-1)^j\sqrt{a^2+1-2a\cos \theta_j}$, and  $\theta_j$ is the 
unique solution in
$(\frac{(j-1)\pi}{k+1},\frac{j\pi}{k+1})$ of the equation
$\sin ((k+1)\theta)-a\sin (k\theta)=0$  for
$j=1,2,\dots, k$.
\end{itemize}
\end{theorem}

\section{Stability analysis and existence of bifurcations}

 Throughout this paper, we assume that
\begin{itemize}
\item[(H1)] $f(0)=0,~f(\cdot)\in C^3(R)$.
\end{itemize}
Denote $x_0(n)=x(n)$, $x_j(n+1)=x_{j-1}(n)$; $y_0(n)=y(n)$,
$y_j(n+1)=y_{j-1}(n)$; $z_0(n)=z(n)$, $z_j(n+1)=z_{j-1}(n)$, $j=1,2,\dots,k$.
Then we can transform system \eqref{a2}  into the
following system of $3k+3$ difference equations without delays
\begin{equation}\label{b1}
\begin{gathered}
x_0(n+1)=ax_0(n)+\beta f(y_k(n))+\beta f(z_k(n)),\\
y_0(n+1)=ay_0(n)+\beta f(z_k(n))+\beta f(x_k(n)),\\
z_0(n+1)=az_0(n)+\beta f(x_k(n))+\beta f(y_k(n)),\\
x_j(n+1)=x_{j-1}(n),\\
y_j(n+1)=y_{j-1}(n),\\
z_j(n+1)=z_{j-1}(n),\quad j=1,2,\dots,k.
\end{gathered}
\end{equation}
For convenience, we denote $c=\beta f'(0)$. The Jacobian matrix of
system \eqref{b1} at the equilibrium $E=(0,\dots,0)$ is as follows
\begin{equation}\label{b2}
A=\begin{bmatrix}
 B & 0 & 0 & c & 0 & c\\
 I_k & 0 & 0 & 0 & 0 & 0\\
 0 & c & B & 0 & 0 & c\\
 0 & 0 & I_k & 0 & 0 & 0\\
  0 & c & 0 & c & B & 0\\
 0 & 0 & 0 & 0 & I_k & 0
\end{bmatrix},
\end{equation}
where $B=(a,0,\dots, 0)_{1\times k}$, $I_k$ is a $k\times k$
identity matrix, $0$ is a zero matrix of appropriate size.

 The associated characteristic equation of system \eqref{b1} is
\begin{equation}\label{b3}
(\lambda^{k+1}-a\lambda^k+c)^2(\lambda^{k+1}-a\lambda^k-2c)=0.
\end{equation}
Applying Theorem \ref{thm2.1} to \eqref{b3} and noting that
$d_1\leq -d_0$, we can obtain the following results.

\begin{theorem} \label{thm3.1} 
 Assume {\rm (H1)} and that $0<a<1$. Then we have 
\begin{itemize}
\item[(i)] if $\max \{-b_0,b_1/2\}<c<b_0/2$,   the origin of system
\eqref{a2} is asymptotically stable;

\item[(ii)]  if $c=-b_0$ or $b_0/2$, a fold bifurcation occurs at the
origin in system \eqref{a2};

\item[(iii)]  if $c=b_1/2$, a Neimark-Sacker bifurcation
occurs at the origin in system \eqref{a2},
where $b_0=1-a$, $b_1=-\sqrt{a^2+1-2a\cos \theta_1}$, 
in which $\theta_1$ is the unique solution in $(0,\frac{\pi}{k+1})$
 of the equation
$\sin ((k+1)\theta)-a\sin (k\theta)=0$.
\end{itemize}
\end{theorem} 

\begin{remark} \label{rmk 3.1} \rm
 As to the case of $c=b_1/2$, if $c=-b_j$ or
$b_j/2$, Neimark-Sacker bifurcations occur at the origin in system
\eqref{a2}, where $b_j=(-1)^j\sqrt{a^2+1-2a\cos \theta_j}$,
in which $\theta_j$ is the unique solution in
$(\frac{(j-1)\pi}{k+1},\frac{j\pi}{k+1})$ of the equation
$\sin ((k+1)\theta)-a\sin (k\theta)=0$ for $j=1,2,\dots,
k$. If $c=-b_{k+1}$ or $b_{k+1}/2$, system \eqref{a2} has a Flip
bifurcation at the origin, where $d_{k+1}=(-1)^{k+1}(1+a)$.
\end{remark}

In what follows, we investigate the global asymptotic stability of
system \eqref{a2}.

\begin{theorem} \label{thm3.2} 
 Under assumption {\rm  (H1)}, the origin of
 \eqref{a2} is globally asymptotically stable if the following
conditions hold.
\begin{itemize}
\item[(H2)]  There exists a constant $L>0$ such
that $|f'(\cdot)|\leq L$.

\item[(H3)] $1-a-2L|\beta|>0$.
\end{itemize}
\end{theorem}

\begin{proof}
 Since $1-a-2L|\beta|>0$, the matrix
$$
A=\begin{pmatrix}
 1-a & -L|\beta| & -L|\beta| \\
 -L|\beta| & 1-a & -L|\beta| \\
 -L|\beta| & -L|\beta| & 1-a
\end{pmatrix}
$$
is an M-matrix, and there exists a vector $p=(p_i)_{1\times 3}>0$
such that $pA>0$ (\cite{s18}); that is,
$$
p_i(1-a)-(\sum_{j=1}^{3}p_j-p_i)L|\beta|>0,~i=1,2,3.
$$
Hence, we can choose $\lambda>1$ such that 
\begin{equation} \label{b4}
p_i(1-a\lambda)-\Big(\sum_{j=1}^{3}p_j-p_i\Big)L|\beta|\lambda^{k+1}>0,
\quad i=1,2,3.
\end{equation}
Let $U_1(n)=\lambda^n|x(n)|$, $U_2(n)=\lambda^n|y(n)|$,
$U_3(n)=\lambda^n|z(n)|$. From \eqref{a2}, we have
\begin{equation}\label{b5}
 U_i(n+1)\leq a\lambda
U_i(n)+L|\beta|\lambda^{k+1}\Big[\sum_{j=1}^{3}U_j(n-k)-U_i(n-k)\Big],\quad
i=1,2,3.
\end{equation}
Define a Lyapunov function
$$
V(n)=\sum_{i=1}^{3}p_iU_i(n)+\sum_{l=n-k}^{n-1}\sum_{i=1}^{3}\Big[
\Big(\sum_{j=1}^{3}p_j-p_i\Big)L|\beta|\lambda^{k+1}U_i(l)\Big].
$$
Then  from \eqref{b4} and \eqref{b5}, we deduce that
$$
\Delta V(n)=V(n+1)-V(n)
=-\sum_{i=1}^{3}\Big[p_i(1-a\lambda)-(\sum_{j=1}^{3}p_j-p_i)L|\beta|
\lambda^{k+1}\Big]U_i(n)
\leq 0,
$$
which implies that $V(n)\leq V(0)$.
Note that
\begin{gather*}
V(n)\geq m_0\lambda^n(|x(n)|+|y(n)|+|z(n)|),\\
V(0)=\sum_{i=1}^{3}p_iU_i(0)+\sum_{l=-k}^{-1}\sum_{i=1}^{3}
\Big[(\sum_{j=1}^{3}p_j-p_i)L|\beta|\lambda^{k+1}U_i(l)\Big]:=M_0,
\end{gather*}
where $m_0=\min_{i=1,2,3}\{p_i\}$, $M_0$ is a positive
constant.
Thus,
$$
|x(n)|+|y(n)|+|z(n)|\leq \frac{M_0}{m_0}\lambda^{-n}.
$$
Noting that $\lambda>1$, we get $\lim_{n\to +\infty}x(n)=0$,
$\lim_{n\to +\infty}y(n)=0$,
$\lim_{n\to +\infty}z(n)=0$. Then, the origin of system \eqref{a2} is globally
attractive. On the other hand, it is easy to verify that $|c|<b_0/2$
under conditions (H2) and (H3). Considering that
$b_1\leq b_0$, we have $\max \{-b_0,b_1/2\}<c<b_0/2$ if
$|c|<b_0/2$. Then the origin of system \eqref{a2} is asymptotically
stable. Consequently, the origin of system \eqref{a2} is globally
asymptotically stable.
\end{proof}

\section{Direction and stability of Neimark-Sacker bifurcation}

 In this section, employing the normal form theory and
the center manifold theorem for discrete-time system developed by
Kuznetsov \cite{s1}, we study the direction of Neimark-Sacker
bifurcation and the stability of periodic solutions bifurcating from
the origin of system \eqref{a2}.

 From the discussions in
Section 3, we know that if $c=b_1/2$, a Neimark-Sacker bifurcation
occurs at the origin in system \eqref{a2}. For convenience, we
denote $b_1/2$ by $\tilde{b}$, and the pair of imaginary roots of the
associated characteristic equation for system \eqref{a2} by
$\mathrm{e}^{\pm i\tilde{\theta}}$.

 Denote $\lambda=\mathrm{e}^{i\tilde{\theta}}$. Let
 $q$ be an eigenvector of $A$ corresponding to eigenvalue
 $\lambda$, and $p$ be an eigenvector of $A^\textrm{T}$ corresponding 
to eigenvalue  $\bar{\lambda}$, that is, $Aq=\lambda q$ and 
$A^\textrm{T}p=\bar{\lambda} p$. Obviously,
 $q$ and $p$ are $3(k+1)$-dimensional complex column vectors.

 By direct calculations, we obtain that
 \begin{equation}\label{d1}
\begin{gathered}
\begin{aligned}
q&=\Big(\frac{2c}{\lambda-a},\frac{2c}{\lambda(\lambda-a)},\dots,
 \frac{2c}{\lambda^{k-1}(\lambda-a)},\frac{2c}{\lambda^k(\lambda-a)},\\
&\quad \lambda^k,\lambda^{k-1},\dots,\lambda, 1,
\lambda^k,\lambda^{k-1},\dots,\lambda, 1 \Big)^\textrm{T},
\end{aligned}\\
\begin{aligned}
p&=E\Big(\frac{2c}{\bar{\lambda}^{2k-1}(\bar{\lambda}-a)^2},
 \frac{2c}{\bar{\lambda}^{2k-1}(\bar{\lambda}-a)},
 \frac{2c}{\bar{\lambda}^{2k-2}(\bar{\lambda}-a)},\dots,
 \frac{2c}{\bar{\lambda}^{k+1}(\bar{\lambda}-a)},\\
&\quad \frac{2c}{\bar{\lambda}^k(\bar{\lambda}-a)},
 \frac{1}{\bar{\lambda}^{k-1}(\bar{\lambda}-a)},\frac{1}{\bar{\lambda}^{k-1}},\frac{1}{\bar{\lambda}^{k-2}},\dots,
\frac{1}{\bar{\lambda}},1,\\
&\quad \frac{1}{\bar{\lambda}^{k-1}(\bar{\lambda}-a)},
 \frac{1}{\bar{\lambda}^{k-1}},\frac{1}{\bar{\lambda}^{k-2}},\dots,
\frac{1}{\bar{\lambda}},1 \Big)^\textrm{T},
\end{aligned}
\end{gathered}
\end{equation}
in which
$$
E=\frac{1}{2k+\frac{2\bar{\lambda}}{\bar{\lambda}-a}
 +\frac{4kc^2}{\bar{\lambda}^{2k}(\bar{\lambda}-a)^2}+
\frac{4c^2}{\bar{\lambda}^{2k-1}(\bar{\lambda}-a)^3}},
$$
where $p$ and $q$ satisfy $\langle p,q\rangle=1$.

System \eqref{b1} can be rewritten as
$$
u(n+1)=Au(n)+F(u(n)),
$$
where $A$ is defined by \eqref{b2}, and
\begin{gather*}
u(n)=(x_0(n),x_1(n),\dots,x_k(n),y_0(n),y_1(n),\dots,y_k(n),z_0(n),z_1(n),
\dots,z_k(n))^\textrm{T},\\
F(u(n))=(F_1(u(n)),0\dots,0,F_2(u(n)),0\dots,0,F_3(u(n)),0\dots,0)^\textrm{T},
\end{gather*}
in which $F_i$ $(i=1,2,3)$ can be expanded in the form
\begin{gather*}
F_1(\xi)=\frac{1}{2}\beta
f''(0)[\xi_{2k+2}^2+\xi_{3k+3}^2]+
\frac{1}{6}\beta f'''(0)[\xi_{2k+2}^3+\xi_{3k+3}^3]+\text{h.o.t.},\\
F_2(\xi)=\frac{1}{2}\beta
f''(0)[\xi_{k+1}^2+\xi_{3k+3}^2]+
\frac{1}{6}\beta f'''(0)[\xi_{k+1}^3+\xi_{3k+3}^3]+\text{h.o.t.},\\
F_3(\xi)=\frac{1}{2}\beta
f''(0)[\xi_{k+1}^2+\xi_{2k+2}^2]+ \frac{1}{6}\beta
f'''(0)[\xi_{k+1}^3+\xi_{2k+2}^3]+\text{h.o.t.}.
\end{gather*}
 Simple calculations yield 
\begin{gather}
B_1(x,y)=\sum_{j,k=1}^{3k+3}\frac{\partial^2F_1(\xi)}{\partial\xi_j\partial\xi_k}\big|_{\xi=0}x_jy_k
=\beta f''(0)[x_{2k+2}y_{2k+2}+x_{3k+3}y_{3k+3}],
\nonumber \\
B_2(x,y)=\sum_{j,k=1}^{3k+3}\frac{\partial^2F_2(\xi)}{\partial\xi_j\partial\xi_k}\big|_{\xi=0}x_jy_k
=\beta f''(0)[x_{k+1}y_{k+1}+x_{3k+3}y_{3k+3}],
\nonumber \\
B_3(x,y)=\sum_{j,k=1}^{3k+3}\frac{\partial^2F_3(\xi)}{\partial\xi_j\partial\xi_k}\big|_{\xi=0}x_jy_k
=\beta f''(0)[x_{k+1}y_{k+1}+x_{2k+2}y_{2k+2}],
\nonumber \\
\begin{aligned}
C_1(x,y,z)&=\sum_{j,k,l=1}^{3k+3}\frac{\partial^3F_1(\xi)}{\partial\xi_j\partial\xi_k\partial\xi_l}
\big|_{\xi=0}x_jy_kz_l
\\
&=\beta f'''(0)[x_{2k+2}y_{2k+2}z_{2k+2}+x_{3k+3}y_{3k+3}z_{3k+3}],
\end{aligned} \nonumber \\
\begin{aligned}
C_2(x,y,z)&=\sum_{j,k,l=1}^{3k+3}\frac{\partial^3F_2(\xi)}{\partial\xi_j\partial\xi_k\partial\xi_l}
\big|_{\xi=0}x_jy_kz_l
\\
&=\beta f'''(0)[x_{k+1}y_{k+1}z_{k+1}+x_{3k+3}y_{3k+3}z_{3k+3}],
\end{aligned} \nonumber \\
\begin{aligned}
C_3(x,y,z)&=\sum_{j,k,l=1}^{3k+3}\frac{\partial^3F_3(\xi)}
{\partial\xi_j\partial\xi_k\partial\xi_l}
\big|_{\xi=0}x_jy_kz_l
\\
&=\beta f'''(0)[x_{k+1}y_{k+1}z_{k+1}+x_{2k+2}y_{2k+2}z_{2k+2}].
\end{aligned}  \label{d2}
\end{gather}
If $b=\tilde{b}$, the restriction of system \eqref{b1} to its
two-dimensional center manifold at the critical parameter value can
be transformed into the normal form written in complex coordinates
(see \cite{s1,s2}):
$$
w\to\lambda w(1+\frac{1}{2}d|w|^2)+O(|w|^4),\quad w\in C,
$$
in which
$$
d=\bar{\lambda}\left\langle p,\hat{C}(q,q,\bar{q})+2\hat{B}
\big(q,(I-A)^{-1}\hat{B}(q,\bar{q})\big)+
\hat{B}\big(\bar{q},(\lambda^2I-A)^{-1}\hat{B}(q,q)\big)\right\rangle,
$$
with $p$ and $q$ defined by \eqref{d1} and
$\hat{B}=(B_1,0,\dots,0,B_2,0,\dots,0,B_3,0,\dots,0,)^\textrm{T}$,
$\hat{C}=(C_1,0,\dots,0,C_2,0,\dots,0,C_3,0,\dots,0,)^\textrm{T}$
with $B_i,~C_i$ $(i=1,2,3)$ defined by \eqref{d2}.

\begin{theorem}[\cite{s2}] \label{thm4.1} 
 The direction and stability
of the Neimark-Sacker bifurcation are determined by the sign of
 $\operatorname{Re}(d)$.  If $\operatorname{Re}(d)<0$,
 then the bifurcation is supercritical, that is, the closed invariant 
curve bifurcating from the origin is asymptotically stable. 
If $\operatorname{Re}(d)>0$,  then the bifurcation is subcritical;
that is, the closed invariant curve bifurcating from the origin is unstable.
\end{theorem}

\section{Numerical simulations}

 In this section, we give two examples to illustrate
the results derived in Sections 3 and 4. In system \eqref{a2}, 
we choose the activation function as the type of inverse tangent 
function or hyperbolic tangent function; i.e., $f(v)=\tanh(v)$, 
then $f'(0)=1$, $f''(0)=0$, $f'''(0)=-2$, and $|f'(x)|\leq 1$. 
In addition, in the following simulations, the numerical results of
 $y(n)$ and $z(n)$ are similar to those
of $x(n)$, so they are omitted.

\begin{example} \label{examp1}\rm
 For system \eqref{a2}, If $a=0.5$, $k=2$, then
$b_0=0.5$, $b_1\approx -0.7808$. Choose $\beta=-0.38$ and $0.24$,
respectively, we have
$\max \{-b_0,b_1/2\}=-0.3904<c<0.25=b_0/2$. By Theorem \ref{thm3.1},
the origin of system \eqref{a2} is asymptotically stable 
(see Figure \ref{fig1}). If $\beta=-0.3904\approx b_1/2$, a Neimark-Sacker 
bifurcation occurs at the origin in system \eqref{a2}. Furthermore, from the
formulae in Sections 3 and 4, and by direct computations, we obtain
$$\tilde{\theta}\approx 0.8758,~d\approx -0.6168+0.1248i.$$
Therefore, $\operatorname{Re}(d)\approx -0.6168<0$, and from 
Theorem \ref{thm4.1},
the Neimark-Sacker bifurcation is supercritical and the bifurcating
periodic solution is asymptotically stable (see Figure \ref{fig2}).
\end{example}

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.49\textwidth]{fig1a} % a1.eps, a2.eps
\includegraphics[width=0.49\textwidth]{fig1b}
\end{center}
\caption{(a) Dynamic behavior of system \eqref{a2} with $a=0.5$, 
$k=2$, $\beta=-0.38$. (b) Dynamic behavior of system \eqref{a2} 
with $a=0.5$, $k=2$, $\beta=0.24$. The origin of system \eqref{a2} 
is asymptotically stable} \label{fig1}
\end{figure}

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.49\textwidth]{fig2a} % b1.eps, b2.eps
\includegraphics[width=0.49\textwidth]{fig2b}
\end{center}
\caption{Dynamic behavior of system \eqref{a2} with $a=0.5$, $k=2$, 
$\beta=-0.3904$. A Neimark-Sacker bifurcation occurs at the origin in system
\eqref{a2}. (a) Diagram in $(n,x)$-plane. (b) Phase portrait}
\label{fig2}
\end{figure}

\begin{example} \label{examp2}\rm
Consider the following generalized system with
different connection weights through the neurons
 \begin{equation}\label{e1}
\begin{gathered}
x(n+1)=ax(n)+\alpha f(y(n-k))+\beta f(z(n-k)),\\
y(n+1)=ay(n)+\alpha f(z(n-k))+\beta f(x(n-k)),\\
z(n+1)=az(n)+\alpha f(x(n-k))+\beta f(y(n-k)).
\end{gathered}
\end{equation}
Let $a=0.5$, $\alpha=-0.4$, and choose $\beta=0.25,0.3,0.5,1.2$,
respectively, numerical simulations show that Neimark-Sacker
bifurcations occur at the origin in system \eqref{e1}. We further
find that with the increasing of $\beta$, the velocity of the
trajectory going to the bifurcating periodic solution is increasing
(see Figure \ref{fig3}). In Figure \ref{fig4}, we exhibit a typical
bifurcation and chaos
diagram when we fix other parameters and choose $\beta$ as a
bifurcation parameter $(-4<\beta<4)$. It clearly shows that system
\eqref{e1} admits rich dynamics including period-doubling
bifurcation and chaos.
\end{example}

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.49\textwidth]{fig3a} %c1.eps
\includegraphics[width=0.49\textwidth]{fig3b} \\
\includegraphics[width=0.49\textwidth]{fig3c}
\includegraphics[width=0.49\textwidth]{fig3d}
\end{center}
\caption{Phase plot in space $(x,y,z)$ for system \eqref{e1} with
(a) $c=0.25$, (b) $c=0.3$, (c) $c=0.5$ and (d) $c=1.2$}
\label{fig3}
\end{figure}

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.6\textwidth]{fig4} %d1.eps
\end{center}
\caption{Bifurcation diagram in $(c,x)$-plane for system \eqref{e1} with
step size of 0.01 for $c$. It shows the effect of parameter $c$ on
the dynamic behavior} \label{fig4}
\end{figure}

\subsection*{Conclusion}

  In this article, a class of discrete-time
bidirectional ring neural network model with delay was investigated.
Analyzing the corresponding characteristic equation, the
asymptotic stability of the origin was discussed, and choosing
$c=\beta f'(0)$ as a bifurcation parameter, we showed that system
\eqref{a2} undergoes Neimark-Sacker bifurcations at the origin. Using the
M-matrix method and the Lyapunov function, global
asymptotic stability of the origin was derived. Applying the
normal form theory and the center manifold theorem, the direction of
the Neimark-Sacker bifurcation and the stability of bifurcating
periodic solution were obtained.


We removed the restrictions of the
conditions that $f'(0) = 1$ and $xf (x) \neq 0$ for $x \neq 0$ in \cite{s5}, 
and derived conditions under which the origin is asymptotically stable. 
Noting that $−b_1 > 1 − a$, the hypotheses of Theorem \ref{thm3.1}
are less restrictive than those for the corresponding continuous system
($\alpha=\beta$) in \cite{s5}.
Moreover, when the connection weights through the neurons
in a bidirectional ring neural network model are different,
numerical simulations show that the  corresponding system still
undergoes Neimark-Sacker bifurcations at the origin. We leave
for future work the study of  \eqref{a2} with different connection
weights, activation functions and time delays.


\subsection*{Acknowledgements} 
This work was supported by the National Natural Science Foundation of 
China (No. 11071254), the Scientific
Research Foundation for the Returned Overseas Chinese Scholars, 
State Education Ministry and the Natural Science Foundation of 
Young Scientist of Hebei Province (No. A2013506012).

\begin{thebibliography}{00}

\bibitem{s9} J. Cao, W. Yu, Y. Qu;
 A new complex network model and convergence dynamics for reputation
computation in virtual organizations, {\it Phys. Lett. A} 
\textbf{356} (2006) 414-425.

\bibitem{s11} Q. Gan, R. Xu, W. Hu, P. Yang;
 Bifurcation analysis for a tri-neuron discrete-time BAM neural
network with delays, {\it Chaos Solitons Fractals} \textbf{42}
(2009) 2502-2511.

\bibitem{s15} S. Guo, Y. Chen;
 Stability and bifurcation of a discrete-time
three-neuron system with delays, {\it Int. J. Appl. Math. Eng. Sci.}
\textbf{1} (2007) 103-115.

\bibitem{s12} S. Guo, X. Tang, L. Huang; Bifurcation
analysis in a discrete-time single-directional network with delays,
{\it Neurocomputing} \textbf{71} (2008) 1422-1435.

\bibitem{s3} J. Hopfield;
 Neural networks and physical systems with
emergent collective computational abilities, {\it Proc. Natl. Acad.
Sci.} \textbf{79} (1982) 2554-2558.

\bibitem{s8} E. Kaslik, St. Balint;
 Bifurcation analysis for a discrete-time Hopfield
neural network of two neurons with two delays and self-connections,
{\it Chaos Solitons Fractals} \textbf{39} (2009) 83-91.

\bibitem{s1} Y. A. Kuznetsov; 
Elements of applied bifurcation theory,
 vol. 112 of Applied Mathematical Sciences, Springer, New York,
USA, 1998.

\bibitem{s2} Y. A. Kuznetsov, H. G. E. Meijer;
 Numerical normal forms for codim 2 bifurcations of fixed points with at most
two critical eigenvalues, {\it SIAM J. Sci. Comput.} 
\textbf{26} (2005) 1932-1954.

\bibitem{s10} S. M. Lee, O. M. Kwonb, Ju H. Park;
 A novel delay-dependent criterion for delayed neural
networks of neutral type, {\it Phys. Lett. A} \textbf{374} (2010)
1843-1848.

\bibitem{s13} X. Li, J. Cao;
 Delay-dependent stability of neural networks of neutral
 type with time delay in the leakage term, 
{\it Nonlinearity} \textbf{23} (2010) 1709-1726.

\bibitem{s16} A. Stuart, A. Humphries; 
Dynamical systems and numerical analysis,
Cambridge University Press, Cambridge, 1996.

\bibitem{s4}  D. Tank, J. Hopfield; 
Simple neural optimization networks: an A/D converter,
signal decision circuit and a linear programming circuit, {\it IEEE
Trans. Circ. Syst.} \textbf{33} (1986) 533-541.

\bibitem{s18} R. S. Varga;
 Matrix iterative analysis, vol.
27 of Springer Series in Computational Mathematics, Springer,
Berlin, Germany, 2000.

\bibitem{s14} B. Wang, J. Jian;
 Stability and Hopf bifurcation analysis on a four-neuron BAM neural
network with distributed delays, {\it Commun. Nonlinear Sci. Numer.
Simul.} \textbf{15} (2010) 189-204.

\bibitem{s5} L. Wang, X. Han;
 Stability and Hopf bifurcation
analysis in bidirectional ring network model, {\it Commun. Nonlinear
Sci. Numer. Simul.} \textbf{16} (2011) 3684-3695.

\bibitem{s7} H. Zhao, L. Wang, C. Ma;
 Hopf bifurcation and stability analysis on discrete-time Hopfield
neural network with delay, {\it Nonlinear Anal. Real World Appl.} \textbf{9} (2008)
103-113.

\bibitem{s6} J. Zhou, S. Li, Z. Yang;
 Global exponential stability of Hopfield neural networks
with distributed delays, {\it Appl. Math. Model.} \textbf{33}
(2009) 1513-1520.

\bibitem{s17} H. Zhu, L. Huang; 
Stability and bifurcation in a tri-neuron network model
with discrete and distributed delays, {\it Appl. Math. Comput.} 
\textbf{188} (2007) 1742-1756.

\end{thebibliography}

\end{document}
