Synchronous dynamics of two coupled oscillators with inhibitory-to-inhibitory connection

Commun Nonlinear Sci Numer Simulat 15 (2010) 4131–4148

Contents lists available at ScienceDirect

Commun Nonlinear Sci Numer Simulat journal homepage: www.elsevier.com/locate/cnsns

Synchronous dynamics of two coupled oscillators with inhibitory-to-inhibitory connection Jian Peng a,b, Shangjiang Guo a,* a b

College of Mathematics and Econometrics, Hunan University, Changsha, Hunan 410082, People’s Republic of China College of Mechanics and Aerospace, Hunan University, Changsha, Hunan 410082, People’s Republic of China

a r t i c l e

i n f o

Article history: Received 25 March 2009 Accepted 7 February 2010 Available online 12 February 2010 Keywords: Coupled oscillators Synchronization Delay Bifurcation Stability Normal form

a b s t r a c t We investigate the behaviour of a neural network model consisting of two coupled oscillators with delays and inhibitory-to-inhibitory connections. We consider the absolute synchronization and show that the connection topology of the network plays a fundamental role in classifying the rich dynamics and bifurcation phenomena. Regarding eigenvalues of the connection matrix as bifurcation parameters, we obtain codimension one bifurcations (including fold bifurcation and Hopf bifurcation) and codimension two bifurcation (including fold-Hopf bifurcations and Hopf–Hopf bifurcations). Based on the normal form theory and center manifold reduction, we obtain detailed information about the bifurcation direction and stability of various bifurcated equilibria as well as periodic solutions with some kinds of spatio-temporal patterns. Numerical simulation is also given to support the obtained results. Ó 2010 Elsevier B.V. All rights reserved.

1. Introduction The main functional unit of oscillatory neural networks is a neural oscillator. Some neural oscillator models consider oscillations as an endogenous property of a pacemaker neuron, such as Van der Pol model [13], Hindmarsh–Rose model [10], and Hodgkin–Huxley model [8]. Another approach suggests that oscillations arise as a result of the interactions between neural populations, for example, between excitatory and inhibitory populations, such as Wilson–Cowan model [21], integrate and fire model [11], and McGregor model [15]. In this paper, the model of a single neural oscillator is represented by a system of two autonomous differential equations describing the dynamics of average activities of the excitatory and the inhibitory populations (measured as a portion of firing neurons in each population). Denoting these activities by E and I, respectively, the model reads



E0 ðtÞ ¼ EðtÞ þ af ðc1 Iðt  sÞÞ; I0 ðtÞ ¼ IðtÞ þ af ðc2 Eðt  sÞÞ;

ð1Þ

where a; c1 ; c2 ; s are positive constants and f : R ! R is a C1 -smooth function. Throughout this paper, we always assume that f ð0Þ ¼ 0 and f 0 ð0Þ ¼ 1. For a number of two neuron models and their linear stability analysis, we refer to the works of Marcus et al. [14], Babcock and Westervelt [1,2], Gopalsamy and Leung [4]. In order to understand complicate phenomena arising from large networks with delayed feedback, it is much interesting to consider two identical neural oscillators given by (1) and coupled using terms that may be interpreted as an additional external input. This system of coupled oscillators has the form * Corresponding author. Fax: +86 731 8823056. E-mail address: [email protected] (S. Guo). 1007-5704/$ - see front matter Ó 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.cnsns.2010.02.008

4132

J. Peng, S. Guo / Commun Nonlinear Sci Numer Simulat 15 (2010) 4131–4148

8 0 E1 ðtÞ ¼ E1 ðtÞ þ af ðc1 I1 ðt  sÞ þ P12 ðt  sÞÞ; > > > < I0 ðtÞ ¼ I ðtÞ þ af ðc E ðt  sÞ þ Q ðt  sÞÞ; 1 2 1 12 1 > E02 ðtÞ ¼ E2 ðtÞ þ af ðc1 I2 ðt  sÞ þ P21 ðt  sÞÞ; > > : 0 I2 ðtÞ ¼ I2 ðtÞ þ af ðc2 E2 ðt  sÞ þ Q 21 ðt  sÞÞ;

ð2Þ

where P 12 ¼ a1 E2  a2 I2 ; P21 ¼ a1 E1  a2 I1 ; Q 12 ¼ a3 E2  a4 I2 ; Q 21 ¼ a3 E1  a4 I1 . Here ðE1 ; I1 Þ and ðE2 ; I2 Þ describe activities of the first and the second oscillators, respectively. The terms P 12 ; P 21 ; Q 12 ; Q 21 describe connections between subpopulations related to different oscillators. Since we have symmetric coupling of identical oscillators, Eq. (2) has the reflection symmetry of interchange of two oscillators. A coupled network of simple networks like (2) can exhibit a range of interesting behaviour, qualitatively very different from their behaviour in isolation, such as synchronization, phase trapping, phase locking, and amplitude death. Among them, synchronization, which implies that the individual networks oscillate at a common frequency and phase when coupled, has been causing researchers’ wide focus in neuroscience field since the works by Pecora and Carroll [16,17]. There is significant evidence that synchronization of neuronal firing play a profound role in the information processing in many brain systems, including the olfactory bulb, hippocampus, and thalamo-cortical system [5]. Such synchronization is often responsible for the formation of certain beneficial biological functions, such as orchestrating the sleep/wake cycle (circadian rhythm) of mammals [18] and a more efficient data transmission in the brain [19]. In this paper, we consider the connections between the inhibitory populations: P 12 ¼ P21 ¼ 0; Q 12 ¼ aI2 ; Q 21 ¼ aI1 , where a is the control parameter. Namely, we consider the following model (see Fig. 1):

8 0 E1 ðtÞ ¼ E1 ðtÞ þ af ðc1 I1 ðt  sÞÞ; > > > 0 < I1 ðtÞ ¼ I1 ðtÞ þ af ðc2 E1 ðt  sÞ  aI2 ðt  sÞÞ; 0 > E > 2 ðtÞ ¼ E2 ðtÞ þ af ðc1 I2 ðt  sÞÞ; > : 0 I2 ðtÞ ¼ I2 ðtÞ þ af ðc2 E2 ðt  sÞ  aI1 ðt  sÞÞ:

ð3Þ

  A natural phase space is the Banach space C ½s; 0; R4 of all continuous mappings from ½s; 0 into R4 equipped with the   super-norm. Then, ut ðhÞ ¼ uðt þ hÞ for h 2 ½s; 0. Every initial state u 2 C ½s; 0; R4 defines uniquely a solution u u u T xu ðtÞ ¼ ðEu 1 ðtÞ; I 1 ðtÞ; E2 ðtÞ; I 2 ðtÞÞ

  of (3) for all t P s, and system (3) generates a semi-flow in C ½s; 0; R4 . A phase state (point) u ¼ ðu1 ; . . . ; u4 ÞT 2 C is said to be synchronous (respectively, anti-synchronous) if ðu1 ; u2 Þ ¼ ðu3 ; u4 Þ (respectively, ðu1 ; u2 Þ ¼ ðu3 ; u4 Þ). Due to the uniqueness of the Cauchy initial value problem of system (3), every synchronous (respectively, anti-synchronous) phase u u u point u gives a synchronous (respectively, anti-synchronous) solution xu of (3), that is, Eu 1 ðtÞ ¼ E2 ðtÞ and I 1 ðtÞ ¼ I 2 ðtÞ u u u u (respectively, E1 ðtÞ ¼ E2 ðtÞ and I1 ðtÞ ¼ I2 ðtÞ) for all t P s. These provide formal definitions of several concepts related to synchrony [20,23]. Another important class of solutions of (3) is a anti-phased periodic solutions which is of period p and   u u u u satisfies E1 ðtÞ ¼ E2 t þ 2p and I1 ðtÞ ¼ I2 t þ 2p for all t P s. Our focus in this paper is on the pattern formation of system (3) for different coupling strength. The study of such problems is important in various areas, for example, in the theory and applications of content addressable memories where a stable solution can be used as coded information of a memory of the system to be stored and retrieved. The rest of this paper is organized as follows. In Section 2, we consider the absolute synchronization of (3). We discuss the associated characteristic equation and obtain criteria ensuring the linear stability of the trivial solution in Section 3. Sections 4 and 5 are devoted to the existence, direction, and stability of the Hopf bifurcation and fold bifurcation analysis, respectively. Section 6 is devoted to codimension two bifurcations including fold-Hopf bifurcation and Hopf–Hopf bifurcations. Some numerical simulations

Fig. 1. Architecture of the model described by (3).

J. Peng, S. Guo / Commun Nonlinear Sci Numer Simulat 15 (2010) 4131–4148

4133

are performed to illustrate the analysis results in Section 7. Appendix A contains some of the details of the calculation of some critical coefficients for the normal forms of Hopf bifurcation. 2. Absolute synchronization In this section, we consider the absolutely synchronization of (3) in the case where f satisfies the normalization, monotonicity, concavity, and boundedness conditions that f 0 ðxÞ > 0 for all x 2 R; f 00 ðxÞx < 0 for all x–0; f 000 ð0Þ < 0; and  1 < limx!1 f ðxÞ < þ1. We now define asymptotical synchronization. Definition 1. A solution xu of (3) is said to be asymptotically synchronous if its x-limit set xðuÞ is contained in the set of synchronous phase given by

n o   w ¼ ðw1 ; . . . ; w4 ÞT 2 C ½s; 0; R4 : w1 ¼ w3 ; w2 ¼ w4 : We say that system (3) is absolutely synchronous if every solution of (3) is asymptotically synchronous, i.e.,

lim jE1 ðtÞ  E2 ðtÞj ¼ 0

t!1

and

lim jI1 ðtÞ  I2 ðtÞj ¼ 0

t!1

for every fixed nonnegative delay s. Theorem 1. If a < 0 and aðc1 þ c2 Þ < 2 or a > 0 and aðc1 þ c2 þ 2aÞ < 2, then every solution of (3) with arbitrarily given s is asymptotically synchronous. Proof. Consider a given solution x : ½s; 1Þ ! R4 of (3) and let yðtÞ ¼ E1 ðtÞ  E2 ðtÞ and zðtÞ ¼ I1 ðtÞ  I2 ðtÞ. Then from (3), we get for all t P 0,

y0 ðtÞ ¼ yðtÞ þ a½f ðc1 I1 ðt  sÞÞ  f ðc1 I2 ðt  sÞÞ ¼ yðtÞ  ac1 p1 ðtÞzðt  sÞ; z0 ðtÞ ¼ zðtÞ þ a½f ðc2 E1 ðt  sÞ  aI2 ðt  sÞÞ  f ðc2 E2 ðt  sÞ  aI1 ðt  sÞÞ ¼ zðtÞ þ ap2 ðtÞðc2 yðt  sÞ þ azðt  sÞÞ; where

p1 ðtÞ ¼

Z

1

f 0 ðsðc1 I1 ðt  sÞÞ þ ð1  sÞðc1 I2 ðt  sÞÞÞ ds 0

and

p2 ðtÞ ¼

Z

1

f 0 ðsðc2 E1 ðt  sÞ  aI2 ðt  sÞÞ þ ð1  sÞðc2 E2 ðt  sÞ  aI1 ðt  sÞÞÞ ds: 0

Due to the properties of f, there exists p 2 ð0; 1 such that p1;2 ðtÞ 6 p for all t P 0. Firstly, if a < 0, then let

VðtÞ ¼ y2 ðtÞ þ z2 ðtÞ þ aðc1  aÞp

Z

t

z2 ðsÞ ds þ ac2 p

Z

ts

t

y2 ðsÞ ds:

ts

Then

d VðtÞ ¼ 2yðtÞðyðtÞ  ac1 p1 ðtÞzðt  sÞÞ þ 2zðtÞðzðtÞ þ ap2 ðtÞðc2 yðt  sÞ þ azðt  sÞÞÞ þ aðc1  aÞp ðz2 ðtÞ  z2 ðt  sÞÞ dt þ ac2 p ðy2 ðtÞ  y2 ðt  sÞÞ 6 2y2 ðtÞ þ ac1 p y2 ðtÞ þ ac1 p z2 ðt  sÞ  2z2 ðtÞ þ ac2 p ðz2 ðtÞ þ y2 ðt  sÞÞ þ aap z2 ðtÞ  aap z2 ðt  sÞ þ aðc1  aÞp ðz2 ðtÞ  z2 ðt  sÞÞ þ ac2 p ðy2 ðtÞ  y2 ðt  sÞÞ ¼ ð2 þ aðc1 þ c2 Þp Þðy2 ðtÞ þ z2 ðtÞÞ: On the other hand, if a > 0, we take

VðtÞ ¼ y2 ðtÞ þ z2 ðtÞ þ aðc1 þ aÞp

Z

t

ts

z2 ðsÞds þ ac2 p

Z

t

y2 ðsÞ ds:

ts

Then

d VðtÞ 6 ð2 þ aðc1 þ c2 Þp Þy2 ðtÞ þ ð2 þ aðc1 þ c2 þ 2aÞp Þz2 ðtÞ: dt

4134

J. Peng, S. Guo / Commun Nonlinear Sci Numer Simulat 15 (2010) 4131–4148

Note that y; z are bounded on ½s; 1Þ, we see that y0 ; z0 are bounded on ½0; 1Þ. This, together with y; z 2 L2 ð½0; 1ÞÞ, implies that limt!þ1 yðtÞ ¼ limt!þ1 zðtÞ ¼ 0. This completes the proof. h 3. Linear stability analysis   As usual, in Banach space C ½s; 0; R4 , if r 2 R; A P 0 and u : ½r  s; r þ A ! R2 is continuous mapping, then ut 2 C for t 2 ½r; r þ A is defined by ut ðhÞ ¼ uðt þ hÞ for s 6 h 6 0. Linearizing system (3) at trivial solution leads to the following linear system:

8 0 E1 ðtÞ ¼ E1 ðtÞ  ac1 I1 ðt  sÞ; > > > < I0 ðtÞ ¼ I ðtÞ þ ac E ðt  sÞ  aaI ðt  sÞ; 1 2 1 2 1 > E02 ðtÞ ¼ E2 ðtÞ  ac1 I2 ðt  sÞ; > > : 0 I2 ðtÞ ¼ I2 ðtÞ þ ac2 E2 ðt  sÞ  aaI1 ðt  sÞ: Let g ¼ 2a a and b ¼ 2a

ð4Þ

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ja2  4c1 c2 j, for reasons explained below. The characteristic matrix for system (4) is

Dðg; b; kÞ ¼ ðk þ 1ÞId  Beks ; where Id is the identity matrix of size 4, and

2

0 6 ac 6 2 B¼6 4 0 0

ac1

0

0

0

0

0

aa

ac2

0

3

aa 7 7 7: ac1 5

ð5Þ

0

Thus, the characteristic equation is det Dðg; b; kÞ ¼ 0, that is,

h i2  2 ðk þ 1Þ2 þ a2 c1 c2 e2ks  ðk þ 1Þaaeks ¼ 0:

ð6Þ

Obviously, det Dðg; b; kÞ can be decomposed as

det Dðg; b; kÞ ¼ Dþ ðg; b; kÞ  D ðg; b; kÞ; where D ðg; b; kÞ ¼ ½k þ 1  ðg þ bÞeks ½k þ 1  ðg  bÞeks  if a2  4c1 c2 > 0, and D ðg; b; kÞ ¼ ½k þ 1  ðg þ ibÞeks ½k þ 1 ðg  ibÞeks  if a2  4c1 c2 < 0. Observe that ðg  bÞ (if a2  4c1 c2 > 0) or ðg  ibÞ (if a2  4c1 c2 < 0) are four eigenvalues of the matrix B. It follows that it is natural to regard g and b as the bifurcation parameters. As a result, it is necessary to determine when the infinitesimal generator Aðg; bÞ of the C 0 -semigroup generated by the linear system (4) has eigenvalues lying on the imaginary axis. For this purpose, we first consider

Pz ðkÞ ¼ k þ 1  zeks ;

ð7Þ

where z 2 C. Motivated by Guo et al. [6], we define a curve R with the following parametric equations:



uðtÞ ¼ cos st  t sin st;

v ðtÞ ¼ t cos st þ sin st;

t 2 R:

ð8Þ

Similar to the analysis in [6], the curve R intersects with the u-axis at ðcn ; 0Þ; n 2 N0 :¼ f0; 1; 2; . . .g, where qffiffiffiffiffiffiffiffiffiffiffiffiffi cn ¼ ð1Þn 1 þ t2n and ftn g1 n¼0 are the monotonic increasing sequence of the nonnegative zeros of v ðtÞ. The curve R is schematically illustrated in Fig. 2. In the sequel, we will identify R with fuðtÞ þ iv ðtÞ : t 2 Rg 2 C. The following lemma plays an important role in analyzing the distribution of zeros of (6). Lemma 1. [6] (i) Pz ðkÞ has a purely imaginary zero if and only if z 2 R. Moreover, if z ¼ uðhÞ þ iv ðhÞ then the purely imaginary zero is ih except that there is a pair of conjugate purely imaginary zeros itn if z ¼ cn for n 2 N. (ii) For each fixed z0 ¼ uðh0 Þ þ iv ðh0 Þ 2 R, there exists an open d-neighborhood of z0 in the complex plane, denoted by Bðz0 ; dÞ, and an analytical function k : Bðz0 ; dÞ ! C such that kðz0 Þ ¼ ih0 and kðzÞ is a zero of Pz ðkÞ for all z 2 Bðz0 ; dÞ. Moreover, along the outward-pointing normal vector to the curve R at z0 , the directional derivative of RekðzÞ at z0 is positive. (iii) Pz ðkÞ has only zeros with strictly negative real parts if and only if z is inside the curve R0 ; exactly j 2 N zeros with positive real parts if z is between Rj1 and Rj . In particular, if z 2 R0 ; Pz ðkÞ has either a simple real zero 0 (if z ¼ 1) or a simple purely imaginary zero (if ImðzÞ–0), or a pair of simple purely imaginary zeros (if z ¼ c1 ), and all other zeros have strictly negative real parts.

4135

J. Peng, S. Guo / Commun Nonlinear Sci Numer Simulat 15 (2010) 4131–4148

10 for t≥0 for t<0

8 Σ2

6 Σ

4

1

v(t)

2 0

Σ

0

γ

γ

γ

3

γ

0

1

2

−2 −4 −6 −8 −10 −10

−5

0

5

10

u(t) Fig. 2. The parametric curve R.

We introduce the following point set in R2 to simplify our presentation:

 R2þ ¼ ðg; bÞ 2 R2 : b > 0 ; n o Xj ¼ ðg; bÞ 2 R2þ : g  b ¼ cj ; g  b–cj and g  b–  ck for all k 2 N0 ; n o Cj ¼ ðg; bÞ 2 R2þ : g  b ¼ cj ; g  b–cj and g  b–  ck for all k 2 N0  for j 2 N0 . We call X j ; Cj ; j 2 N0 , critical lines.

Corollary 1. Suppose a2  4c1 c2 > 0, then,

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi (i) All zeros of det Dðg; b; kÞ have negative real parts if and only if a a þ a2  4c1 c2 < 2. S S S (ii) If and only if ðg; bÞ 2 Xþ Xn Cþn Cn for some n 2 N; det Dðg; b; kÞ has a pair of simple conjugate purely imaginary n zeros, which are it n . S S þS  X0 C0 C0 ; det Dðg; b; kÞ has a simple zero k ¼ 0. Moreover, if ðg; bÞ 2 Xþ0 , then all zeros but (iii) If and only if ðg; bÞ 2 Xþ 0 k ¼ 0 of det Dðg; b; kÞ have strictly negative real parts.

Corollary 2. Suppose a2  4c1 c2 < 0, then, (i) det Dðg; b; kÞ has a purely imaginary zero if and only if ðg; bÞ or  ðg; bÞ lies on the curve R. The purely imaginary zero is give by ih, where h satisfies uðhÞ ¼ g and v ðhÞ ¼ b or uðhÞ ¼ g and v ðhÞ ¼ b. (ii) All zeros of det Dðg; b; kÞ have strictly negative real parts if and only if a2 c1 c2 < 1. pffiffiffiffiffiffiffiffiffiffiffiffiffiffi (iii) If 1 < a2 c1 c2 < 1 þ n2 , then det Dðg; b; kÞ has exactly a pair of simple purely imaginary zeroes and all other zeros has pffiffiffiffiffiffiffiffiffiffiffiffiffiffi strictly negative real parts; If a2 c1 c2 > 1 þ n2 , then det Dðg; b; kÞ has exactly a pair of simple purely imaginary zeroes and at least two zeros with positive real parts, where n is the minimal positive solution to equation n tanðsnÞ ¼ 1. In view of Corollaries 1 and 2, there are stability switches at the trivial solution as ðg; bÞ crosses the critical lines

Xj and Cj ðj 2 N0 Þ or the curves R. Then, system (3) may undergo fold bifurcations or Hopf bifurcations. In fact, we have Theorem 2. Suppose that a2  4c1 c2 > 0, then

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi (i) The trivial solution of system (3) is asymptotically stable if and only if a a þ a2  4c1 c2 < 2.  S S S þS  Xn Cn Cn a Hopf bifurcation occurs in system (3) , i.e., a unique branch of periodic solu(ii) Near each ðg; bÞ 2 n2N Xþ n tions bifurcates from the origin. S S þS  X0 C0 C0 a fold bifurcation (or saddle-node bifurcation) occurs in system (3) , i.e., a unique (iii) Near each ðg; bÞ 2 Xþ 0 branch of equilibria bifurcates from the

origin. 

 1þc 1c 1cj 1þc ;  2 j for j 2 N, a fold-Hopf bifurcation occurs in system (3). (iv) Near ðg; bÞ equal to one of the pairs 2 j ;  2 j ; 2

4136

J. Peng, S. Guo / Commun Nonlinear Sci Numer Simulat 15 (2010) 4131–4148



(v) Near ðg; bÞ equal to either cation occurs in system (3).

cn þcj cn cj 2

;



2

or



cn cj cn þcj 2

;

2



for some two different natural numbers n and j, a Hopf–Hopf bifur-

Theorem 3. Suppose that a2  4c1 c2 < 0, then (i) The trivial solution of system (3) is asymptotically stable if and only if a2 c1 c2 < 1. (ii) If either ðg; bÞ or ðg; bÞ lies on the curve R, then near ðg; bÞ a Hopf bifurcation occurs in system (3) , i.e., a unique branch of periodic solutions bifurcates from the origin.

Remark 1. Under assumption that a2  4c1 c2 < 0, b > 0 and hence g þ ib cannot be equal to ±1. Thus, fold bifurcations and hence fold-Hopf bifurcations cannot occur in system (3) when a2  4c1 c2 < 0. In the subsequently three sections, we consider two kinds of bifurcations: codimension one and codimension two bifurcations. 4. Hopf bifurcation In the previous section, we obtained some conditions on existence of Hopf bifurcation in system (3). It is much interesting to consider the spatio-temporal patterns of bifurcating periodic solutions. For this purpose, we need to explore the possible (spatial) symmetry of the system (3). Consider the following action of Z2 on R4 as:

q  ðE1 ; I1 ; E2 ; I2 Þ ¼ ðE2 ; I2 ; E1 ; I1 Þ:

ð9Þ

Obviously, system (3) is Z2 -equivariant. In view of the general symmetric Hopf bifurcation theorem has been developed by Wu [22] and Guo and Lamb [7], we have the following result. S  Theorem 4. Suppose that a2  4c1 c2 > 0, then the trivial solution of system (3) undergoes a Hopf bifurcation at ðg; bÞ 2 Cþ Cn n S (respectively, ðg; bÞ 2 Xþ X n n ) for some n 2 N, giving rise to one branch of synchronous (respectively, anti-phased) periodic solutions. Proof. On one hand, if ðg; bÞ 2 Cþ n ðhÞ, where spanned by qðhÞ and q

qðhÞ ¼ ð1; d; 1; dÞT eitn h ;

S

Cn for some n 2 N, then the center space U of (4) associated to eigenvalues itn is

h 2 ½s; 0;

ð10Þ 1

1

where d ¼ cn =ðac1 Þ. Obviously, U is Z2  S -isomorphic to C. It follows that the representation of Z2  S on C is given by

qz¼z

h  z ¼ eitn h z:

and

ð11Þ

Namely, the isotropic subgroup of Z2  S1 is only Z2 ðqÞ, which implies that the bifurcated periodic solutions are synchronous, i.e., taking the form xðtÞ ¼ ðuðtÞ; v ðtÞ; uðtÞ; v ðtÞÞT with u(t) and v(t) are periodic functions with the same minimal positive period. S  Xn for some n 2 N, then the center space U of (4) associated to eigenvalues itn is On the other hand, if ðg; bÞ 2 Xþ n ðhÞ, where spanned by qðhÞ and q

qðhÞ ¼ ð1; d; 1; dÞT eitn h ;

h 2 ½s; 0:

1

ð12Þ 1

Obviously, U is Z2  S -isomorphic to C. It follows that the representation of Z2  S on C is given by

q  z ¼ z

and

h  z ¼ eitn h z:

ð13Þ 1

Namely, the isotropic subgroup of Z2  S is only Z2 ðq; pÞ, which implies that the bifurcated periodic solutions are antiphased, i.e., taking the form



p

pT xðtÞ ¼ uðtÞ; v ðtÞ; u t þ ; v t þ ; 2 2 where p is a period of x(t).

h

By using a similar argument, we can verify the following result. Theorem 5. Suppose that a2  4c1 c2 < 0, then the trivial solution of system (3) undergoes a Hopf bifurcation at ðg; bÞ 2 R (respectively, ðg; bÞ 2 R), giving rise to one branch of anti-phased (respectively, synchronous) periodic solutions. In what follows, we consider the stability and bifurcation direction. We first focus on the following two cases: S  Cn for some n 2 N. Case (i): a2  4c1 c2 > 0 and ðg0 ; b0 Þ 2 Cþ n Case (ii): a2  4c1 c2 < 0 and ðg0 ; b0 Þ 2 R.

J. Peng, S. Guo / Commun Nonlinear Sci Numer Simulat 15 (2010) 4131–4148

4137

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi In Case (i), either g0 þ b0 ¼ cn or g0  b0 ¼ cn . In Case (ii), g20 þ b20 ¼ 1 þ x20 , where x0 ¼ a2 c1 c2  1. For the sake of simplification, let l ¼ g  b and l0 ¼ cn for Case (i), and l ¼ ðg; bÞ and l0 ¼ ðg0 ; b0 Þ for Case (ii). Let kðlÞ be the analytical function satisfying kðl0 Þ ¼ x0 and kðlÞ þ 1 þ lekðlÞs 0 for all sufficiently small jlj, where x0 is equal to tj in Case (ii). d Refkðl0 Þg < 0 in Case (ii), where # is the outObviously, it follows from Lemma 1(ii) that cn Rek0 ðl0 Þ < 0 in Case (i) and d# ward-pointing normal vector to the curve R at l0 . Thus, in both cases, k ¼ ix0 are the imaginary zeros of det Dðg0 ; b0 ; Þ. Based on the center manifold reduction and normal form approach introduced by Hassard et al. [9], we will compute the reduced system on the center manifold with the pair of conjugate complex, purely imaginary eigenvalues ix0 . By this reduction we can determine the Hopf bifurcation direction, i.e., to answer the question of whether the bifurcation branch of periodic solution exists locally for supercritical bifurcation or subcritical bifurcation. For this purpose, we further assume that f 2 C 3 ðR; RÞ and xf ðxÞ–0 when x–0. Letting uðtÞ ¼ ðE1 ðtÞ; I1 ðtÞ; E2 ðtÞ; I2 ðtÞÞT and ut ðhÞ ¼ uðt þ hÞ for h 2 ½s; 0, we can rewrite (3) as

_ uðtÞ ¼ Ll ut þ Gðut ; lÞ

ð14Þ

with Ll u ¼ uð0Þ þ BuðsÞ and

3 3 2 c21 u22 ðsÞ c31 u32 ðsÞ 7 7 6 6 6 ½c2 u1 ðsÞ  au4 ðsÞ2 7 a 000 6 ½c2 u1 ðsÞ  au4 ðsÞ3 7 a 7 þ f ð0Þ6 7 þ Oðjuj4 Þ: Gðu; lÞ ¼ f 00 ð0Þ6 7 6 7 6 6 2 c21 u24 ðsÞ c31 u34 ðsÞ 5 5 4 4 2

½c2 u3 ðsÞ  au2 ðsÞ2

½c2 u3 ðsÞ  au2 ðsÞ3

By the Riesz representation theorem, there exists a matrix whose entries are bounded variation functions Nðh; lÞ in h 2 ½s; 0 such that

Ll u ¼

Z

0

  dNðh; lÞuðhÞ for u 2 C ½s; 0; R4 :

s

In fact, we can choose

Nðh; lÞ ¼



h ¼ 0;

Id;

Bdðh þ sÞ; h 2 ½s; 0Þ;

  where dðhÞ is the Dirac delta function. Next, we define for u 2 C 1 ½s; 0; R4 ,

( du dh

Al u ¼

R0

if h 2 ½s; 0Þ;

;

s

dNðh; lÞuðhÞ ¼ Ll u; if h ¼ 0;

and

Rl u ¼ Since

dut dh



if h 2 ½s; 0Þ;

0;

Gðu; lÞ; if h ¼ 0:

t ¼ du , we can rewrite (14) in the form: dt

u_ t ¼ Al ut þ Rl ut :

ð15Þ

By direct computation, we can choose

qðhÞ ¼ ð1; d; 1; dÞT eix0 h ;

h 2 ½s; 0;

where d ¼ ð1 þ ix0 Þeix0 s =ðac1 Þ is an eigenvector of Al0 associated with the eigenvalue ix0 . i.e.,

Al0 qðhÞ ¼ ix0 qðhÞ: The adjoint operator Al0 is defined by



 Al0 w ðnÞ ¼

(

 dw ; if n 2 ð0; s; dn R0 wðtÞ dNðt; 0Þ; if n ¼ 0: s

ð16Þ

4138

J. Peng, S. Guo / Commun Nonlinear Sci Numer Simulat 15 (2010) 4131–4148

Note that the domains of Al0 and Al0 are C 1 ð½s; 0; C4 Þ and C 1 ð½0; s; C4 Þ, where C4 is the space of four-dimensional complex row vectors. It follows from (16) that ix0 is an eigenvalue of Al0 , and

Al0 q ðnÞ ¼ ix0 q ðnÞ

ð17Þ

for some nonzero row-vector function q ðnÞ; n 2 ½0; s. In order to construct coordinates to describe the center manifold Cl near the origin, we use the bilinear form

 uð0Þ  hw; ui ¼ wð0Þ

Z

0

h¼s

Z

h

  hÞ dNðh; l ÞuðnÞ dn wðn 0

ð18Þ

n¼0

E



D   for w 2 Cð½0; s; C4 Þ and u 2 Cð½s; 0; C4 Þ. Then, as usual, w; Al0 u ¼ Al0 w; u for ðu; wÞ 2 Dom Al0  Dom Al0 . We i ¼ 0. By direct computation, we obtain that normalize q and q by the condition hq ; qi ¼ 1 and hq ; q

    ix n q ðnÞ ¼ D b; 1; b; 1 e 0 ;

n 2 ½0; s 1

where b ¼ ac2 eix0 s =ð1 þ ix0 Þ and D ¼ f2ðb þ dÞ½1 þ sð1 þ ix0 Þg .   For each l 2 Dom Al0 , we may then associate with it the pair (z, w), where z ¼ hq ; ui and w ¼ u  2Refzqg. For a solution u(t) of (15) at l ¼ l0 . We define zðtÞ ¼ hq ; ut i and wðt; hÞ ¼ ut ðhÞ  2RefzðtÞqðhÞg. In fact, z and z are local coordinates for Cl0 in the directions of q and q . It is easy to see that hq ; wi ¼ 0. Now, for solutions ut 2 Cl0 of (16), hq ; u_ t i ¼ q ; Al0 ut þ Rl0 ut i. Then, when l ¼ l0 ,

z_ ðtÞ ¼ ix0 zðtÞ þ g ðz; zÞ;

ð19Þ

 ð0Þf0 ðz; zÞ and f 0 ðz; zÞ ¼ Gðwðz; z; hÞ þ 2RefzqðhÞg; l0 Þ. Let where g ðz; zÞ ¼ q

g ðz; zÞ ¼ g 20

z2 z2 z2 z þ : þ g 11 zz þ g 02 þ g 21 6 2 2

We leave the detailed derivations of the coefficients of g ðz; zÞ to Appendix A. Based on the above analysis, we can see that each g ij is determined by the parameters in system (3). Thus, we can compute the following quantities:

 1 1 g 20 g 11  2jg 11 j2  jg 02 j2 þ g 21 ; 2x0 3 2 RefC 1 ðl0 Þg ; r¼ Refk0 ðl0 Þg C 1 ð l0 Þ ¼

i



d ¼ 2RefC 1 ðl0 Þg: Therefore, following Hassard et al. [9], we obtain the following statements: r determines the directions of the Hopf bifurcation: if r > 0 (respectively, <0), then the Hopf bifurcation is supercritical (respectively, subcritical); d determines the stability of the bifurcating periodic solutions: the bifurcating periodic solutions are stable (respectively, unstable) if d < 0 (respectively, >0). If we further assume that: 00

(H) f : R ! R is a sufficiently smooth odd function satisfying xf ðxÞ < 0 for all x–0. Then, we have the following result. Theorem 6. Under condition (H) and assumptions of Case (i), let

 1 h  2 i aa h1 ¼ c1 c2 þ a2 c1 c4n 2 þ

cn

system (3) undergoes a Hopf bifurcation. The direction of Hopf bifurcation and stability of bifurcated synchronous periodic solutions are determined by the sign of h1 f 000 ð0Þ. More precisely, if h1 f 000 ð0Þ < 0 (respectively, >0), then the bifurcated periodic solutions have the same stability as the trivial solution had before the bifurcation (respectively, unstable), and Hopf bifurcation is supercritical (respectively, subcritical), i.e., the bifurcated periodic solutions exist for ðg; bÞ near ðg0 ; b0 Þ satisfying jg  bj > jcn j (respectively, < jcn j). Proof. Under assumption (H), we have g 20 ¼ g 11 ¼ g 02 ¼ 0 and

h  2  i 000 g 21 ðl0 Þ ¼ 2Df ð0Þa MN2 Nc31 eix0 s þ c2 eix0 s  aN c2 eix0 s  aN ;

J. Peng, S. Guo / Commun Nonlinear Sci Numer Simulat 15 (2010) 4131–4148

4139

where M ¼ ac2 =ð1 þ ix0 Þ and N ¼ ð1 þ ix0 Þ=ðac1 Þ. In view of Dþ ðg0 ; b0 ; ix0 Þ ¼ 0 and 1 þ ix0 þ ðg  bÞeix0 s ¼ 0, we have

  h  i f 000 ð0Þ 2RefC 1 ð0Þg ¼ Refg 21 g ¼ Re ac21 c2 jNj2 eix0 s þ jc2 eix0 s  aNj2 ac2 eix0 s  aaN ðb þ dÞ½1 þ sð1 þ ix0 Þ ( ) 1 i 000 f ð0Þ½1 þ sð1 þ ix0 Þ h 2 3 2 2ix0 s 2 2 ix0 s ¼ Re a c1 c2 jNj e þ jc2 e  aNj ð1 þ ix0 Þ 2ð1 þ ix0 Þ þ aaeix0 s ( ) 1 i f 000 ð0Þ½1 þ sð1 þ ix0 Þ h 2 2 2 ix0 s ix0 s ¼ Re c1 jNj ð1 þ ix0 þ aae Þ þ jc2 e  aNj ð1 þ ix0 Þ 2 þ aaeix0 s ð1 þ ix0 Þ1 (     1 ) f 000 ð0Þð1 þ s þ sx20 þ x0 iÞ 2 2 aa aa 2 ix0 s ¼ Re c jNj 1 þ e  a Nj þ jc 2 þ 2 1 cn cn ð1 þ sÞ2 þ s2 x20 ( "    2 # 1 ) 000 2 2 f ð0Þð1 þ s þ sx0 þ x0 iÞ cn aa ac aa 2þ ¼ Re 1þ þ c2 þ n cn cn a2 ac1 ð1 þ sÞ2 þ s2 x20 ( 1 ) i  2 2 2 f 000 ð0Þð1 þ s þ sx20 þ x0 iÞ h 2 2 aa 2 2þ ¼ Re a ðcn þ aacn Þ þ a c1 ða c1 c2 þ aacn Þ cn ð1 þ sÞ2 þ s2 x20 ( )   1  2 i f 000 ð0Þð1 þ s þ sx20 þ x0 iÞ h aa : ¼ Re c1 c2 þ a2 c1 c4n 2 þ 2 2 2 cn ð1 þ sÞ þ s x0

Some of the above equalities follows from the fact that c2n þ aacn þ a2 c1 c2 ¼ 0. Hence, RefC 1 ð0Þg has the same sign as that of h1 f 000 ð0Þ. Namely, if h1 f 000 ð0Þ > 0 (respectively, <0), then d ¼ 2RefC 1 ð0Þg > 0 (respectively, <0). In addition, it follows from Lemma 1(ii) that cn Rek0 ðl0 Þ < 0. Thus, if cn h1 f 000 ð0Þ > 0 (respectively, <0), then r ¼ RefC 1 ð0Þg=Refk0 ðl0 Þg > 0 (respectively, <0). The conclusions follow immediately the sign of cn . h In Case (ii), 1 þ ix0 þ ðg  ibÞeix0 s ¼ 0, then

  h  i f 000 ð0Þ ac21 c2 jNj2 eix0 s þ jc2 eix0 s  aNj2 ac2 eix0 s  aaN Refg 21 g ¼ Re ðb þ dÞ½1 þ sð1 þ ix0 Þ   h i f 000 ð0Þ ac1 eix0 s 2 ix0 s 2 2 ix0 s ix0 s ¼ Re c jNj e þ jc e  a Nj ac e  a a N ac 2 2 1 2 1 þ sð1 þ ix0 Þ aaeix0 s þ 2ð1 þ ix0 Þ  h 000 f ð0Þ 1 a2 c31 c2 jNj2 e2ix0 s þ jc2 eix0 s  aNj2 ¼ Re 1 þ sð1 þ ix0 Þ aaeix0 s þ 2ð1 þ ix0 Þ    a2 c1 c2 e2ix0 s þ aað1 þ ix0 Þeix0 s  h f 000 ð0Þ 1 ¼ Re c2 jNj2 ðaað1 þ ix0 Þeix0 s i x s 1 þ sð1 þ ix0 Þ aae 0 þ 2ð1 þ ix0 Þ 1

io þ 1 þ ix0 Þ2 Þ þ jc2 eix0 s  aNj2 ð1 þ ix0 Þ2  h

 f 000 ð0Þ 1 ¼ Re c21 jNj2 aað1 þ ix0 Þ2 =ðg  ibÞ þ ð1 þ ix0 Þ2 1 þ sð1 þ ix0 Þ aað1 þ ix0 Þ=ðg  ibÞ þ 2ð1 þ ix0 Þ  io  þ c2 eix0 s  aNj2 ð1 þ ix0 Þ2  h i f 000 ð0Þ 1 þ ix0 ¼ Re c21 jNj2 ðaa=ðg  ibÞ þ 1Þ þ jc2 eix0 s  aNj2 1 þ sð1 þ ix0 Þ aa=ðg  ibÞ þ 2  000     f ð0Þð1 þ ix0 Þ i  ¼ Re c1 c2 þ c22 g  ib c1 c2 þ c22 1 þ sð1 þ ix0 Þ 2b ( )   i f 000 ð0Þ½x0 þ ð1 þ s þ sx20 Þi haa  2 2 c  c1 c2  ib c1 c2 þ c2 ¼ Re 2 2 2b½ð1 þ sÞ2 þ s2 x20       1   f 000 ð0Þ i b 1 þ s þ sx20 c1 c2 þ c22 þ aax0 c1 c2  c22 : ¼ h 2 2b ð1 þ sÞ2 þ s2 x20 Then, we have the following result. Theorem 7. Under condition (H) and assumptions of Case (ii), let

   1   h2 ¼ b 1 þ s þ sx20 c1 c2 þ c22 þ aax0 c1 c2  c22 2 system (3) undergoes a Hopf bifurcation. The direction of Hopf bifurcation and stability of bifurcated synchronous periodic solutions are determined by the sign of h2 f 000 ð0Þ. More precisely, if h2 f 000 ð0Þ > 0 (respectively, <0), then the bifurcated periodic solutions have

4140

J. Peng, S. Guo / Commun Nonlinear Sci Numer Simulat 15 (2010) 4131–4148

the same stability as the trivial solution had before the bifurcation (respectively, unstable), and Hopf bifurcation is supercritical (respectively, subcritical), i.e., the bifurcated periodic solutions exist for g2 þ b2 > a2 c1 c2 (respectively, < a2 c1 c2 ). We end this section with the discussion focus on the other two cases. Case (iii): a2  4c1 c2 > 0 and ðg0 ; b0 Þ 2 Xþ n Case (iv): a2  4c1 c2 < 0 and ðg0 ; b0 Þ 2 R.

S

Xn for some n 2 N.

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi In Case (iii), either g0 þ b0 ¼ cn or g0  b0 ¼ cn . In Case (iv), g20 þ b20 ¼ 1 þ x20 , where x0 ¼ a2 c1 c2  1. We can discuss it   1; b;  1eix0 n ; n 2 ½0; s, where similarly as above by choosing qðhÞ ¼ ð1; d; 1; dÞT eix0 h ; h 2 ½s; 0 and q ðnÞ ¼ D1 b; 1 D1 ¼ f2ðb þ dÞ½1 þ sð1 þ ix0 Þg . Noting that

8 E ðt þ hÞ ¼ zðtÞeix0 h þ zðtÞeix0 h þ wð1Þ ðt; hÞ; > > > 1 > < I ðt þ hÞ ¼ zðtÞNeix0 ðhþsÞ þ zðtÞNeix0 ðhþsÞ þ wð2Þ ðt; hÞ; 1 > E2 ðt þ hÞ ¼ zðtÞeix0 h  zðtÞeix0 h þ wð3Þ ðt; hÞ; > > > : I2 ðt þ hÞ ¼ zðtÞNeix0 ðhþsÞ  zðtÞNeix0 ðhþsÞ þ wð4Þ ðt; hÞ:

we have g 20 ¼ g 11 ¼ g 02 ¼ 0 and

h  2  i g 21 ¼ 2D1 f 000 ð0Þa c31 MN2 Neix0 s  c2 eix0 s þ aN c2 eix0 s þ aN :

In view of D ðg0 ; b0 ; ix0 Þ ¼ 0, we have

  h  i f 000 ð0Þ 2RefC 1 ð0Þg ¼ Refg 21 g ¼ Re ac21 c2 jNj2 eix0 s  jc2 eix0 s þ aNj2 ac2 eix0 s þ aaN ðb þ dÞ½1 þ sð1 þ ix0 Þ ( ) 1 i f 000 ð0Þ½1 þ sð1 þ ix0 Þ h 2 3 2 2 2ix0 s 2 ix0 s ¼ Re a c1 c2 jNj e  jc2 e þ aNj ð1 þ ix0 Þ 2ð1 þ ix0 Þ þ aaeix0 s ( ) 1 h i f 000 ð0Þ½1 þ sð1 þ ix0 Þ 2 2 2 ix0 s ix0 s : ¼ Re c jNj ða a e  1  i x Þ  jc e þ a Nj ð1 þ i x Þ 0 2 0 1 2 þ aaeix0 s ð1 þ ix0 Þ1 In Case (iii), 1 þ ix0  ðg  bÞeix0 s ¼ 0, then

Refg 21 g ¼ Re

(  f 000 ð0Þð1 þ s þ sx20 þ x0 iÞ

   1 ) aa aa c21 jNj2 1  þ jc2 eix0 s þ aNj2 2 

cn cn ð1 þ sÞ2 þ s2 x20 ( "    2 # 1 ) 000 2 2 f ð0Þð1 þ s þ sx0 þ x0 iÞ cn aa ac aa 2 ¼ Re 1 þ c2  n cn cn a2 ac1 ð1 þ sÞ2 þ s2 x20 ( 1 ) h i  2 2 2 f 000 ð0Þð1 þ s þ sx20 þ x0 iÞ 2 2 aa 2 2 ¼ Re a ðcn  aacn Þ þ a c1 ða c1 c2  aacn Þ cn ð1 þ sÞ2 þ s2 x20 ( )   1  2 i f 000 ð0Þð1 þ s þ sx20 þ x0 iÞ h aa : ¼ Re c1 c2 þ a2 c1 c4n 2  2 2 2 cn ð1 þ sÞ þ s x0

Some of the above equalities follows from the fact that c2n  aacn þ a2 c1 c2 ¼ 0. In Case (iv), 1 þ ix0  ðg þ ibÞeix0 s ¼ 0, then

(

)  1  2 Refg 21 g ¼ Re c1 c2 ðaa  g  ibÞ  c2 ðg þ ibÞ 2bi ð1 þ sÞ2 þ s2 x20 ( ) f 000 ð0Þ½ð1 þ s þ sx20 Þi  x0  2 ¼ Re ½ðc1 c2 ðg  ibÞ  c2 ðg þ ibÞ 2b½ð1 þ sÞ2 þ s2 x20      1   f 000 ð0Þ 2 2 2 s þ s x Þ c c þ c ax c c  c bð1 þ a  : ¼ 1 2 0 1 2 0 2 2 2 2b½ð1 þ sÞ2 þ s2 x20  f 000 ð0Þð1 þ s þ sx20 þ x0 iÞ

Similarly, we have the following results. Theorem 8. Under condition (H) and assumptions of Case (iii), let

 1 h  2 i aa k1 ¼ c1 c2 þ a2 c1 c4n 2 

cn

J. Peng, S. Guo / Commun Nonlinear Sci Numer Simulat 15 (2010) 4131–4148

4141

system (3) undergoes a Hopf bifurcation. The direction of Hopf bifurcation and stability of bifurcated anti-phased periodic solutions are determined by the sign of k1 f 000 ð0Þ. More precisely, if k1 f 000 ð0Þ < 0 (respectively, >0), then the bifurcated periodic solutions have the same stability as the trivial solution had before the bifurcation (respectively, unstable), and Hopf bifurcation is supercritical (subcritical), i.e., the bifurcated periodic solutions exist for ðg; bÞ near ðg0 ; b0 Þ satisfying jg  bj > jcn j (respectively, < jcn j). Theorem 9. Under condition (H) and assumptions of Case (iv), let

   1   k2 ¼ b 1 þ s þ sx20 c1 c2 þ c22  aax0 c1 c2  c22 2 system (3) undergoes a Hopf bifurcation. The direction of Hopf bifurcation and stability of bifurcated anti-phased periodic solutions are determined by the sign of k2 f 000 ð0Þ. More precisely, if k2 f 000 ð0Þ > 0 (respectively, <0), then the bifurcated periodic solutions have the same stability as the trivial solution had before the bifurcation (respectively, unstable), and Hopf bifurcation is supercritical (subcritical), i.e., the bifurcated periodic solution exist for g2 þ b2 > a2 c1 c2 (respectively, < a2 c1 c2 ). 5. Fold bifurcation S S þS  X0 C0 C0 , which as usual is As stated in Section 3, when a2  4c1 c2 > 0, a fold bifurcation may occur at ðg; bÞ 2 Xþ 0 associated with the occurrence of new branches of nontrivial equilibria of (3), emanating from the trivial branch at ðg; bÞ. S  C0 . Then, either g þ b ¼ 1 or g  b ¼ 1. For the sake of convenience, let First, we focus on the case where ðg; bÞ 2 Cþ 0 l ¼ g  b þ 1. Similarly to Section 4, we rewrite (3) as (14), define associated operators Ll ; Al ; and Al , and also introduce the adjoint bilinear form h; i as (18). Denote by kðlÞ the simple real eigenvalue of Al satisfying kð0Þ ¼ 0. Let q and p be the eigenvectors for A0 and A0 associated with 0, respectively. In fact, we can choose qðhÞ ¼ ð1; d0 ; 1; d0 ÞT for h 2 ½s; 0 and pðnÞ ¼ D0 ðb0 ; 1; b0 ; 1Þ for n 2 ½0; s, where d0 ¼ 1=ðac1 Þ; b0 ¼ ac2 and D0 ¼ ½2ðb0 þ d0 Þð1 þ sÞ1 . Obviously, hp; qi ¼ 1. We associate each X 2 DomðAl Þ with the pair (x,w), where x ¼ hp; X i and w ¼ X  xq. For a solution X t of (15), we define x ¼ hp; X t i and wðx; lÞ ¼ X t  xq. It is easy to see that hp; wi ¼ 0. Now, on the center manifold Cl , we have

x_ ¼ gðx; lÞ;

ð20Þ



P where gðx; lÞ ¼ p; Al ðwðx; lÞ þ xqÞ þ Rl ðwðx; lÞ þ xqÞ . Let gðx; lÞ ¼ jP1 1j! g j ðlÞxj . By a direct computation, we have

g 1 ðlÞ ¼ lk0 ðlÞ þ Oðjlj2 Þ; g 2 ð0Þ ¼ f 00 ð0Þ½ac1 ðaa þ 2Þð1 þ sÞ1 ða þ 2ac1 c2 a þ a2 c21 c22 þ ac1 c2 Þ; and

g 3 ð0Þ ¼



i 3ac1 f 00 ð0Þ h ð2Þ ð4Þ ð1Þ ð4Þ ð3Þ ð2Þ c1 c2 w2 þ w2 þ ðc2  ad0 Þ c2 w2  aw2 þ c2 w2  aw2 4ðaa þ 2Þð1 þ sÞ   1 c1 c2 þ f 000 ð0Þ½ðaa þ 2Þð1 þ sÞ1 5 2  : a a c1

We still need to compute w2 ðhÞ for h 2 ½s; 0, that is similar to the computation in [6]. Next, we assume that f 00 ð0Þ–0, then the truncated form of (20) is

1 x_ ¼ lk0 ð0Þx þ g 2 ð0Þx2 : 2

ð21Þ

Therefore, (21) is equivalent to the following equation:

  1 x_ ¼ lk0 ð0Þx  f 00 ð0Þ½ac1 ðaa þ 2Þð1 þ sÞ1 a þ 2ac1 c2 a þ a2 c21 c22 þ ac1 c2 x2 : 2

ð22Þ

0

lk ð0Þ It is easy to see that (22) has two equilibria: x1 ¼ 0 and x2 ¼  f 002ð0Þg . Moreover, if 2 ð0Þ

lk0 ð0Þ > 0 then x1 is unstable and x2 is stable. If lk0 ð0Þ < 0 then x1 is stable and x2 is unstable. These two equilibria coalesce at l ¼ 0. Thus, we obtain a transcritical bifurcation of equilibria of system (3). Theorem 10. If f 00 ð0Þ–0, then near ðg0 ; b0 Þ 2 C 0 system (3) undergoes a transcritical bifurcation. Namely, besides the trivial solution, system (3) has a nonzero equilibrium, which continuously depends on ðg; bÞ for all sufficiently small jg  b þ 1j. Moreover, this nonzero equilibrium is stable if g  b > 1 and jg0  b0 j < g0  b0 ¼ 1, and unstable otherwise. However, under the condition (H), the truncated form of (20) is

1 x_ ¼ ðg  b þ 1Þð1 þ sÞ1 x þ g 3 ð0Þx3 ; 6

ð23Þ

4142

J. Peng, S. Guo / Commun Nonlinear Sci Numer Simulat 15 (2010) 4131–4148

i  c1ac2 . It is easy to see that (23) has only one equilibrium x1 ¼ 0 if pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðg  b þ 1Þf 000 ð0Þ < 0, and three equilibria x1 ¼ 0 and x2;3 ¼  6ð1  g  bÞ=½g 3 ð0Þð1 þ sÞ otherwise. Namely, there exists 000 a pitchfork bifurcation at l ¼ 0. More precisely, if f ð0Þ > 0, then for g þ b < 1, (23) has the stable equilibrium x1 ; for g  b > 1; x1 is still an equilibrium, but two new equilibria x2 and x3 appear. In this process, x1 becomes unstable for g  b > 1 while the other two equilibria are stable. Therefore, we have the following conclusion: where g 3 ð0Þ ¼ f 000 ð0Þ½ðaa þ 2Þð1 þ sÞ1

h

1 a5 c21

Theorem 11. Under the condition (H), near ðg0 ; b0 Þ 2 C 0 , system (3) undergoes a pitchfork bifurcation near 0. More precisely, (i) If f 000 ð0Þ < 0, two nontrivial equilibria exist for ðg; bÞ satisfying g  b > 1 (which are stable if g  b < 1 and unstable otherwise) and only the trivial equilibrium continues existing for g  b < 1. Moreover, the two nontrivial equilibria coalesce into zero as g  b ! 1. (ii) If f 000 ð0Þ > 0, two nontrivial equilibria exist for ðg; bÞ satisfying g  b < 1 (which are unstable) and only the trivial equilibrium continues existing for g  b > 1. Moreover, the two nontrivial equilibria coalesce into zero as g  b ! 1. S  Finally, the bifurcation case that ðg; bÞ 2 Xþ X0 can be discussed analogously by choosing qðhÞ ¼ ð1; d0 ; 1; d0 ÞT for 0 h 2 ½s; 0 and pðnÞ ¼ D2 ðb0 ; 1; b0 ; 1Þ for n 2 ½0; s, where D2 ¼ ½2ðb0 þ d0 Þð1 þ sÞ1 . Then, by computation, we have

g 1 ðlÞ ¼ lk0 ðlÞ þ Oðjlj2 Þ;

g 2 ð0Þ ¼ 0;

and

g 3 ð0Þ ¼



i 3ac1 f 00 ð0Þ h ð2Þ ð4Þ ð1Þ ð4Þ ð3Þ ð2Þ c1 c2 w2 þ w2  ðc2 þ ad0 Þ c2 w2  aw2 þ c2 w2  aw2 4ðaa  2Þð1 þ sÞ  1  2f 000 ð0Þ ðaa  2Þð1 þ sÞa5 c21 :

Similarly, we have the following results. Theorem 12. Under the condition (H), near ðg0 ; b0 Þ 2 X 0 , system (3) undergoes a pitchfork bifurcation near 0. More precisely, (i) If f 000 ð0Þ > 0, two nontrivial equilibria exist for ðg; bÞ satisfying g  b > 1 (which are stable if g  b < 1 and unstable otherwise) and only the trivial equilibrium continues existing for g  b < 1. Moreover, the two nontrivial equilibria coalesce into zero as g  b ! 1. (ii) If f 000 ð0Þ < 0, two nontrivial equilibria exist for ðg; bÞ satisfying g  b < 1 (which are unstable) and only the trivial equilibrium continues existing for g  b > 1. Moreover, the two nontrivial equilibria coalesce into zero as g  b ! 1.

Remark 2. Under the conditions of Theorem 12, further assume that f : R ! R is a sufficiently smooth odd function satis00 fying xf ðxÞ < 0 for all x–0. Then the bifurcated two nontrivial equilibria stated in Theorem 12 are both anti-synchronous. 6. Codimension two bifurcation In this section, we list the possible codimension two bifurcations (including a fold-Hopf bifurcation and a Hopf–Hopf bifurcation) in system (3). 6.1. Fold-Hopf (zero-pair) bifurcation It follows from Theorem 2(iv), near ðg; bÞ equal to one of the pairs



1þcn 2



 ;  12cn ; 12cn ;  1þ2cn for n 2 N, a fold-Hopf bifur-

cation occurs in system (3). Here, we only focus on the fold-Hopf bifurcation occurs when g ¼ 1þ2cn and b ¼ 12cn . For the sake

 of convenience, let l ¼ ðl1 ; l2 Þ ¼ g  1þ2cn ; b  12cn . Similarly to Section 4, we rewrite (3) as (14), define associated operators Ll ; Al ; and Al , and also introduce the bilinear form h; i as (18). For all sufficiently small jlj, let real-valued k0 ðlÞ and complex-valued k1 ðlÞ be eigenvalues of AðlÞ, satisfying k0 ð0Þ ¼ 0 and k1 ð0Þ ¼ it n , with associated eigenvectors   l l l l Q0 ; Q1 2 C ½s; 0; R4 , respectively. In view of the discussion in the previous sections, we can choose Q0 and Q1 such that

Q00 ðhÞ ¼ ð1; d0 ; 1; d0 ÞT ;

Q01 ðhÞ ¼ ð1; d1 ; 1; d1 ÞT eitn h ;

h 2 ½s; 0:

Moreover, let p0 and p1 be the eigenvectors of the adjoint operator A0 associated to eigenvalues 0 and it n . Normalize the







eigenvectors such that p0 ; Q00 ¼ p1 ; Q01 ¼ 1 and p1 ; Q00 ¼ p0 ; Q01 ¼ 0, then we have

p0 ðhÞ ¼ D0 ðb0 ; 1; b0 ; 1Þ;

    itn n ; p1 ðhÞ ¼ D1 b 1 ; 1; b1 ; 1 e

n 2 ½0; s; 1

where dk ¼ ck =ðac1 Þ; bk ¼ ac2 ck and Dk ¼ f2ðbk þ dk Þ½1 þ sð1 þ it k Þg

for k ¼ 0; 1.

4143

J. Peng, S. Guo / Commun Nonlinear Sci Numer Simulat 15 (2010) 4131–4148

 For a solution ut of (15) at l, we define xðtÞ ¼ hp0 ; ut i; zðtÞ ¼ hp1 ; ut i; and wðx; z; z; lÞ ¼ X t  xQ00 ðlÞ  2Re zQ01 ðlÞ . In the coordinates x 2 R and z 2 C, system (15) reads

_ xðtÞ ¼ g ðx; z; z; lÞ; z_ ðtÞ ¼ it n zðtÞ þ hðx; z; z; lÞ;

ð24Þ

where the smooth functions g and h are given as follows:





g ðx; z; z; lÞ ¼ p0 ; Al ðwðx; z; z; lÞ þ xq0 þ 2Refzq1 gÞ þ p0 ; Rl ðwðx; z; z; lÞ þ xq0 þ 2Refzq1 gÞ ;



hðx; z; z; lÞ ¼ p1 ; Al ðwðx; z; z; lÞ þ xq0 þ 2Refzq1 gÞ  it n z þ p1 ; Rl ðwðx; z; z; lÞ þ xq0 þ 2Refzq1 gÞ : System (24) is a system of ODEs and has a fold-Hopf singularity at center manifold is as follows:

l ¼ 0. Thus, near l ¼ 0 the critical normal form on the

 4  y_ ¼ e10 ðlÞy þ e11 ðlÞy2 þ e12 ðlÞjfj þ O ðy; f; fÞ ;

 4  f_ ¼ it n f þ e20 ðlÞf þ e21 ðlÞyf þ e22 ðlÞy2 f þ O ðy; f; fÞ ; where the expressions for eij come from Kuznentsov [12] and Guo et al. [6]. 6.2. Hopf–Hopf bifurcation In this section, we consider another type of codimension two bifurcation: Hopf–Hopf bifurcation. It follows from Theorem



 c þc c c c c c þc 2(v) that, near ðg; bÞ equal to either n 2 j ; n 2 j or n 2 j ; n 2 j for some two different natural numbers n and j, a Hopf–Hopf c þc

c c

bifurcation occurs in system (3). Here, we only focus on the Hopf–Hopf bifurcation occurs when g ¼ n 2 j and b ¼ n 2 j . For

 c þc c c the sake of convenience, let l ¼ ðl1 ; l2 Þ ¼ g  n 2 j ; b  n 2 j . Similarly to Section 4, we rewrite (3) as (14), define associated operators Ll ; Al ; and Al , and also introduce the bilinear form h; i as (18). For all sufficiently small jlj, Al has simple complex-valued eigenvalues k1 ðlÞ and k2 ðlÞ, satisfying kj ð0Þ ¼ ixj ; j ¼ 1; 2, with associated eigenvectors   l l l Qj 2 C ½s; 0; R4 ; j ¼ 1; 2, respectively. In view of the discussion in the previous sections, we can choose Q1 and Q2 such that

Q0j ðhÞ ¼ ð1; dj ; 1; dj ÞT eixj h ;

h 2 ½s; 0;

j ¼ 1; 2:

Moreover, let p1 and p2 be the eigenvectors of the adjoint operator A0 associated to eigenvalues ix1 and ix2 . Normalize the







eigenvectors such that p1 ; Q01 ¼ p2 ; Q02 ¼ 1 and p1 ; Q02 ¼ p2 ; Q01 ¼ 0, then we have

    ixj n ; pj ðhÞ ¼ Dj b j ; 1; bj ; 1 e

n 2 ½0; s;

where dj ¼ cj =ðac1 Þ; bj ¼ ac2 cj and Dj ¼ f2ðbj þ dj Þ½1 þ sð1 þ ixj Þg1 for j ¼ 1; 2.

0.8

E I

0.6

1

1

E

0.4

I

2

2

0.2 0 −0.2 −0.4 −0.6 −0.8

0

10

20

30

40

50

t Fig. 3. System (26) is asymptotically synchronous when a ¼ c2 ¼ 1; c1 ¼ 0:9; a ¼ 2.

4144

J. Peng, S. Guo / Commun Nonlinear Sci Numer Simulat 15 (2010) 4131–4148

We associate each X 2 DomðAðlÞÞ with ðz1 ; z1 ; z1 ; z1 ; wÞ. For a solution ut of (15) at l, we define zj ðtÞ ¼ 

pj ; X t ; j ¼ 1; 2; and wðx; z; z; lÞ ¼ X t  2Re z1 Q01 ðlÞ þ z2 Q02 ðlÞ . In these coordinates, system (15) takes the form

z_ 1 ðtÞ ¼ ix1 z1 ðtÞ þ g ðz1 ; z1 ; z2 ; z2 ; lÞ; z_ 2 ðtÞ ¼ ix2 z2 ðtÞ þ hðz1 ; z1 ; z2 ; z2 ; lÞ;

ð25Þ

where

D

E g ðz1 ; z1 ; z2 ; z2 ; lÞ ¼ p1 ; Rl z1 Q01 þ z1 Q01 þ z2 Q02 þ z2 Q02 ; l ; D

E hðz1 ; z1 ; z2 ; z2 ; lÞ ¼ p2 ; Rl z1 Q01 þ z1 Q01 þ z2 Q02 þ z2 Q02 ; l : The functions g and h are complex-value smooth functions of their arguments and have formal Taylor expansions with respect to the first four arguments. Thus, near l ¼ 0 the critical normal form on the center manifold is as follows:

0.6

E

1

I1

0.4

E I

2

2

0.2 0 −0.2 −0.4 −0.6

Fig. 4. The trivial solution of (26) is asymptotically stable for all

0.8

30

t

s P 0 when a a þ

E

I1

1

I1

0.6

E2

I

2

0.2 0

50

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a2  4c1 c2 < 2. Here a ¼ 0:5; a ¼ 2:1; c1 ¼ c2 ¼ 1; and

1

1

0

0

−1 −0.5

I2

0.4

40

2

20

E

10

0 E1

−1 −0.5

0.5

1

1

0

0

2

0

E

−0.8

−1 −0.5

−0.4

1

1

0

0

−0.6 −0.8

−1 −0.5 0

50

100

150

of

(26)

is

asymptotically

stable

for

all

−1 −0.5

0.5

2

0

I

1

t Fig. 5. The trivial solution 0:8; c2 ¼ 1:2; and s ¼ 0:3927.

0 E1

I

I2

−0.2

0.5

1

−1 −1

0 E1

0

−0.5

I1

E

s ¼ 0:3927.

0.5

0.5

1

0

0.5

2

s P 0 when a2 c1 c2 < 1 and a2 < 4c1 c2 . Here a ¼ 1; a ¼ 1:75; c1 ¼

4145

J. Peng, S. Guo / Commun Nonlinear Sci Numer Simulat 15 (2010) 4131–4148









m_ 1 ¼ ix1 m1 þ e10 ðlÞm1 þ e11 ðlÞjm1 j2 m1 þ e12 ðlÞjm2 j2 m1 þ O kmk5 ; m_ 2 ¼ ix2 m2 þ e20 ðlÞm2 þ e21 ðlÞjm2 j2 m2 þ e22 ðlÞjm1 j2 m2 þ O kmk5 ; where the expressions for eij come from Guo et al. [6] and Elphick et al. [3]. 7. Numerical simulation Let us now give some numerical simulations to illustrate the above results. For this purpose, we must specify the activation function. Throughout this section we always take f ðuÞ ¼ tanhðuÞ, i.e., we consider following system:

8 0 E1 ðtÞ ¼ E1 ðtÞ þ a tanhðc1 I1 ðt  sÞÞ; > > > 0 < I1 ðtÞ ¼ I1 ðtÞ þ a tanhðc2 E1 ðt  sÞ  aI2 ðt  sÞÞ; 0 > E > 2 ðtÞ ¼ E2 ðtÞ þ a tanhðc1 I2 ðt  sÞÞ; > : 0 I2 ðtÞ ¼ I2 ðtÞ þ a tanhðc2 E2 ðt  sÞ  aI1 ðt  sÞÞ;

0.6

I1

1

I1

E2

0.4

1

1

0

0

E2

E

−1 −1

I2

0.2

I2

0 −0.2

−0.5

1

0

0

−0.5

−0.4

0

0

150

200

I

1

0.4

E2

1

1

0

0

−1 −0.4

I2

0.2

I2

0 −0.2

I2

−0.4 −0.6

Fig. 7. When ðg0 ; b0 Þ ¼ ð0:81; 2:0756Þ and

200

0 E1

0.2

−1 −0.4

0.4

1

0

0

−0.2

0 E1

0.2

−1 −0.2

0.4

0.2

0.2

0

0

−0.2 −0.2 150

−0.2

0.2

−0.2 −0.4

100 t

−1 −1

1

E2

1

E1

50

0.5

−0.5

0

E1

0.5

−0.5

0 I1

0.5

1

−0.5

0 E2

0.5

1

s ¼ 0:3927, a branch of anti-phased periodic solutions is bifurcated from the trivial solution.

0.6

0

0 I1

2

100 t

−0.5

E

50

Fig. 6. When ðg0 ; b0 Þ ¼ ð0:81; 2:0756Þ and

−0.8

1

I2

0

−1 −1

0.5

1

−1 −1

I

−1

E1

0

I2

2

I

−0.8

−1 −1

0.5

1

−1 −1

−0.6

E1

0

E2

0.8

ð26Þ

0

0.2 I 1

0.4

0.6

−0.2 −1

−0.2

0 E1

0.2

0.4

0

0.2 I1

0.4

0.6

−0.5

E

0

2

s ¼ 0:3927, a branch of synchronous periodic solutions is bifurcated from the trivial solution.

0.5

4146

J. Peng, S. Guo / Commun Nonlinear Sci Numer Simulat 15 (2010) 4131–4148 000

where a; c1 ; c2 ; s and a are positive constants. Notice that f 000 ð0Þ ¼ tanh ð0Þ ¼ 2 < 0. Obviously, the function tanh x satisfies the normalization, concavity, boundedness conditions, from Theorem 1, we know that, when a ¼ c2 ¼ 1; c1 ¼ 0:9; a ¼ 2, every solution of (26) with arbitrarily given time delay s is asymptotically synchronous, which can be shown by Fig. 3.

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Taking a ¼ 0:5; a ¼ 2:1; c1 ¼ c2 ¼ 1, by a simple calculation, we have a a þ a2  4c1 c2 < 2, which together with Theorem 2, implies that the origin is delay-independently asymptotically stable (see Fig. 4). Taking a ¼ 1; a ¼ 1:75; c1 ¼ 0:8; c2 ¼ 1:2, it is easy to verify that a2 c1 c2 < 1 and a2 < 4c1 c2 . Hence, it follows from Theorem 3 that the origin is delay-independently asymptotically stable (see Fig. 5). Note that when either of ðg; bÞ and ðg; bÞ lies on the curve R, system (26) undergoes a Hopf bifurcation. So taking c1 ¼ c2 ¼ 1; a ¼ 1:2; a ¼ 1:35 (or 1.35), by Theorems 7 and 9, we know that the corresponding waveform and phase plot for s ¼ p=8 are shown in Figs. 6 and 7. All the numeric simulations presented here were performed by the package of MATLAB. Acknowledgements This work was partially supported by National Natural Science Foundation of China (Grant No.10971057), by the Scientific Research Foundation for the Returned Overseas Chinese Scholars (Grant No. [2008]890) and the Key Project (Grant No. [2009]41) of Education Ministry of China. Appendix A. Derivation of normal form (19) The purpose of this appendix is to derive the coefficients of g ðz; zÞ. By (15) and (19), we have

(

¼ _ ¼ u_ t  z_ q  z_ q w

Al0 w  2Refq ð0Þf0 ðz; zqÞg; h 2 ½s; 0; Al0 w  2Refq ð0Þf0 ðz; zqÞg þ f0 ðz; zÞ; h ¼ 0:

ðA:1Þ

We rewrite (A.1) as

_ ¼ Al0 w þ Hðz; z; hÞ; w

ðA:2Þ

where

z2 z2 z3 þ w11 ðhÞzz þ w02 ðhÞ þ w30 ðhÞ þ    ; 2 2 6 z2 z2 z3 Hðz; z; hÞ ¼ H20 ðhÞ þ H11 ðhÞzz þ H02 ðhÞ þ H30 ðhÞ þ    2 2 6 wðz; z; hÞ ¼ w20 ðhÞ

Expanding the above series and comparing the coefficients, we obtain



 Al0  2ix0 w20 ðhÞ ¼ H20 ðhÞ;

Al0 w11 ðhÞ ¼ H11 ðhÞ;   Al0 þ 2ix0 w02 ðhÞ ¼ H02 ðhÞ;   

ðA:3Þ

It is easy to see that

8 E ðt þ hÞ ¼ zðtÞeix0 h þ zðtÞeix0 h þ wð1Þ ðt; hÞ; > > > 1 > < I ðt þ hÞ ¼ zðtÞNeix0 ðhþsÞ þ zðtÞNeix0 ðhþsÞ þ wð2Þ ðt; hÞ; 1 > E2 ðt þ hÞ ¼ zðtÞeix0 h þ zðtÞeix0 h þ wð3Þ ðt; hÞ; > > > : I2 ðt þ hÞ ¼ zðtÞNeix0 ðhþsÞ þ zðtÞNeix0 ðhþsÞ þ wð4Þ ðt; hÞ;

ðA:4Þ

where wðiÞ ðt; hÞ denotes the ith component of wðt; hÞ, i.e., ðiÞ

wðiÞ ðt; hÞ ¼ w20 ðhÞ

z2 ðtÞ z2 ðtÞ ðiÞ ðiÞ þ w11 ðhÞzðtÞzðtÞ þ w02 ðhÞ þ ; 2 2

i ¼ 1; 2; 3; 4:

 ð0ÞGðwðz; z; hÞ þ 2RefzqðhÞg; 0Þ, then we have Note that g ðz; zÞ ¼ q

h i 00 g 20 ¼ 2Df ð0Þa MN2 c21 eix0 s þ ðc2 eix0 s  aNÞ2 ;    00 g 11 ¼ 2Df ð0Þa MNNc21 eix0 s þ ðc2 eix0 s  aNÞ c2 eix0 s  aN ; h i  2 00 g 02 ¼ 2Df ð0Þa MN2 c21 eix0 s þ c2 eix0 s  aN ; and

h

   

00 ð2Þ ð4Þ ð2Þ ð4Þ ð1Þ ð2Þ ð3Þ ð4Þ g 21 ¼ Df ð0Þa Mc21 eix0 s Nw11 þ Nw11 þ Nw20 þ Nw20 þ c2 eix0 s  aN c2 w11  aw11 þ c2 w11  aw11 i h  

 2  i 000 ð1Þ ð2Þ ð3Þ ð4Þ þ c2 eix0 s  aN c2 w20  aw20 þ c2 w20  aw20 þ 2Df ð0Þa c31 MN2 Neix0 s þ c2 eix0 s  aN c2 eix0 s  aN :

J. Peng, S. Guo / Commun Nonlinear Sci Numer Simulat 15 (2010) 4131–4148

4147

We still need to compute w11 and w20 for h 2 ½s; 0Þ. Indeed, we have

ðhÞ Hðz; z; hÞ ¼ 2Refq ð0Þf0 ðz; zÞqðhÞg ¼ g ðz; zÞqðhÞ  gðz; zÞq     2 2 2 z z2 z z z z2 ðhÞ: þ    qðhÞ  g20 þ g11 zz þ g02 þ    q ¼  g 20 þ g 11 zz þ g 02 þ g 21 6 2 2 2 2 Comparing the coefficients with (A.2) gives

ðhÞ; H20 ðhÞ ¼ g 20 qðhÞ  g02 q

ðhÞ: H11 ðhÞ ¼ g 11 qðhÞ  g11 q

It follows from (A.3) that

ð0Þeix0 h : _ 20 ðhÞ ¼ 2ix0 w20 ðhÞ þ g 20 qð0Þeix0 h þ g02 q w Solving for w20 ðhÞ, we obtain

w20 ðhÞ ¼ 

ð0Þ ix h g 20 qð0Þ ix0 h g02 q e  e 0 þ e1 e2ix0 h ; 3ix0 i x0

ðA:5Þ

and similarly, we have

w11 ðhÞ ¼ where

ð0Þ ix h g 11 qð0Þ ix0 h g11 q e  e 0 þ e2 ; ix0 ix0

ðA:6Þ

e1 and e2 are both four-dimensional vectors, and can be determined by setting h ¼ 0 in Hðz; z; hÞ. In fact, from   z2 z2 z2 z Hðz; z; 0Þ ¼  2Refq ð0Þf0 ðz; zÞqð0Þg þ f0 ðz; zÞ ¼  g 20 þ g 11 zz þ g 02 þ g 21 þ    qð0Þ 6 2 2 3 2  2 2 c1 zN þ zN þ wð2Þ 6     2 7 7 6   7 6 c2 zeix0 s þ zeix0 s þ wð1Þ  a zN þ zN þ wð4Þ z2 z2 00 7 þ ; ð0Þ þ af ð0Þ6  g20 þ g11 zz þ g02 þ    q   7 6 2 2 2 7 6 c21 zN þ zN þ wð4Þ 5 4   ix s    2 ix0 s ð3Þ ð2Þ 0    a zN þ zN þ w c2 ze þ ze þw

we have 00

ð0Þ þ af ð0Þn20 H20 ¼ g 20 qð0Þ  g02 q and 00

ð0Þ þ af ð0Þn11 ; H11 ¼ g 11 qð0Þ  g11 q where

2

c21 N2

3

7 6 6 ðc2 eix0 s  aNÞ2 7 7; n20 ¼ 6 7 6 c21 N2 5 4

3 c21 NN  7 6 6 ðc2 eix0 s  aNÞ c2 eix0 s  aN 7 7: ¼6 7 6 c21 NN 5 4  ix s   ix s  c2 e 0  aN c2 e 0  aN 2

n11

ðc2 eix0 s  aNÞ2

From (A.3) and the definition of Al0 , we have

w20 ð0Þ þ Bw20 ðsÞ ¼ 2ix0 w20 ð0Þ  H20 ð0Þ

ðA:7Þ

and

w11 ð0Þ þ Bw20 ðsÞ ¼ H11 ð0Þ: Substituting (A.5) into (A.7) and noticing that Dðl; ix0 Þ ¼ 0, we have e2 ¼ af 00 ð0ÞD1 ð0; 0Þn11 .

ðA:8Þ 00

1

e1 ¼ af ð0ÞD ð0; 2ix0 Þn20 . Similarly, we have

References [1] Babcock KL, Westervelt RM. Complex dynamics in simple neural circuits. In: Denker JS, editor. Neural networks for computing. New York: American Institute of Physics; 1986. p. 288–93. [2] Babcock KL, Westervelt RM. Dynamics of simple electronic neural networks with added inertia. Physica D 1987;23:464–9. [3] Elphick C, Tirapegui E, Brachet ME, Coullet P, Iooss G. A simple global characterization for normal forms of singular vector fields. Physica D 1987;29:95–127. [4] Gopalsamy K, Leung I. Delay induced periodicity in a neural netlet of excitation and inhibition. Physica D 1996;89:395–426. [5] Gray CM. Synchronous oscillations in neuronal systems: Mechanism and functions. J Comput Neurosci 1994;1:11–38. [6] Guo S, Chen Y, Wu J. Two-parameter bifurcations in a network of two neurons with multiple delays. J Diff Eqns 2008;244:444–86. [7] Guo S, Lamb JSW. Equivariant Hopf bifurcation for neutral functional differential equations. Proc Am Math Soc 2008;136:2031–41.

4148

J. Peng, S. Guo / Commun Nonlinear Sci Numer Simulat 15 (2010) 4131–4148

[8] Hansel D, Mato G, Meunier C. Phase dynamics for weakly coupled Hodgkin–Huxley neurons. Europhys Lett 1993;23:367–72. [9] Hassard B, Kazarinoff N, Wan YH. Theory of applications of Hopf bifurcation. London math, society lecture notes, series, vol. 41. Cambridge: Cambridge University Press; 1981. [10] Hindmarsh JL, Rose RM. A model of the nerve impulse using two first-order differential equations. Nature 1982;296:162–4. [11] Kryukov VI, Borisyuk GN, Borisyuk RM, Kirillov AB, Kovalenko EL. Metastable and unstable states in the brain. In: Dobrushin RL, Kryukov VI, Toom AL, editors. Stochastic cellular systems: ergodicity, memory, morphogenesis. Manchester: Manchester University Press; 1990. p. 225–357. [12] Kuznentsov YA. Elements of applied bifurcation theory. 2nd ed. New York: Springer; 1998. [13] Kawato M, Sokabe M, Suzuki R. Synergism and antagonism of neurons caused by an electrical synapse. Biol Cybern 1979;34:81–9. [14] Marcus CM, Waugh FR, Westervelt RM. Nonliear dynamics and stability of analog neural networks. Physica D 1991;51:234–47. [15] McGregor RJ. Neural and brain modelling. New York: Academic Press; 1987. [16] Pecora LM, Carroll TL. Synchronization in chaotic systems. Phys Rev Lett 1990;64:821–4. [17] Pecora LM, Carroll TL. Driving systems with chaotic signals. Phys Rev A 1991;44:2374–83. [18] Reppert SM, Weaver DR. Coordination of circadian timing in mammals. Nature 2002;418:935–41. [19] Samonds JM, Allison JD, Brown HA, Bonds AB. Cooperative synchronized assemblies enhance orientation discrimination. Proc Natl Acad Sci USA 2004;101:6722–7. [20] Stewart I, Golubitsky M, Pivato M. Symmetry groupoids and patterns of synchrony in coupled cell networks. SIAM J Appl Dyn Syst 2003;2:609–46. [21] Wilson HR, Cowan JD. Excitatory and inhibitory interactions in localized populations of model neurons. Biophys J 1972;12:1–24. [22] Wu J. Symmetric functional differential equations and neural networks with memory. Trans Am Math Soc 1999;350:4799–838. [23] Wu J, Faria T, Huang YS. Synchronization and stable phase-locking in a network of neurons with memory. Math Comput Model 1999;30:117–38.