Stability and bifurcation analysis of two-neuron network with discrete and distributed delays

Author’s Accepted Manuscript Stability and Bifurcation Analysis of Two-Neuron Network with Discrete and Distributed Delays Esra Karaoğlu, Enes Yılmaz, Hüseyin Merdan

www.elsevier.com/locate/neucom

PII: DOI: Reference:

S0925-2312(15)01916-5 http://dx.doi.org/10.1016/j.neucom.2015.12.006 NEUCOM16507

To appear in: Neurocomputing Received date: 8 May 2015 Revised date: 6 October 2015 Accepted date: 8 December 2015 Cite this article as: Esra Karaoğlu, Enes Yılmaz and Hüseyin Merdan, Stability and Bifurcation Analysis of Two-Neuron Network with Discrete and Distributed Delays, Neurocomputing, http://dx.doi.org/10.1016/j.neucom.2015.12.006 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Stability and Bifurcation Analysis of Two-Neuron Network with Discrete and Distributed Delays Esra Karao˘glua,c,1, Enes Yılmazb,∗, H¨ useyin Merdanc a b

University of Pittsburgh, Department of Mathematics, Pittsburgh, PA 15260, USA Gazi University, Polatlı Faculty of Science and Arts, Department of Mathematics, 06900 Ankara, Turkey c TOBB University of Economics and Technology, Faculty of Science and Letters, Department of Mathematics, 06530, Ankara, Turkey

Abstract In this paper, we give a detailed Hopf bifurcation analysis of a recurrent neural network system involving both discrete and distributed delays. Choosing the sum of the discrete delay terms as a bifurcation parameter the existence of Hopf bifurcation is demonstrated. In particular, the formulae determining the direction of the bifurcations and the stability of the bifurcating periodic solutions are studied by using the normal form theory and center manifold theorem. Finally, numerical simulations supporting our theoretical results are also included. Keywords: Hopf bifurcation, stability, neural networks, delay differential equations, periodic solution 1. Introduction Recently, research of recurrent neural networks (RNNs) (especially Hopfield neural networks) has attracted the attention of great number of investigators [8, 10, 13, 20, 22, 31, 32]. Major contributions to the theory and design of neural networks dates back to the 1980s [13]. Hopfield [14] ∗

Corresponding author Email addresses: [email protected], [email protected] (Esra Karao˘ glu ), [email protected], [email protected] (Enes Yılmaz ), [email protected] (H¨ useyin Merdan) 1 E. Karao˘ glu was supported by TUBITAK (The Scientific and Technological Research Council of Turkey). Preprint submitted to Elsevier

December 14, 2015

has been a pioneer study in recurrent networks with symmetric synaptic connections using energy function. Based on the Hopfield network model, Babcock and Westervelt [2] have studied the dynamics of simple electronic neural networks and analyzed two neurons with discrete time delays. Marcus and Westervelt [18] have considered the stability of analog neural networks with delayed response and showed that continuous time analog networks can exhibit sustained oscillation when time delay is included. Afterwards, the dynamic properties of neural networks have been studied very much in the literature [3, 7, 13, 21, 23, 29]. RNNs which generate stable oscillations have been used to model certain biological phenomena (see [28] and the references therein). Physiological experiments suggest that brain has chaotic structure. If this chaotic behaviour changes, as in Alzheimer disease, the brain could be slower on rapid state transitions essential for information processing [9]. For epilepsy disease, it has been assumed that the stabilization of unstable patterns in the healthy chaotic brain dynamics is the cause of the increased rhythmicity observed in EEG recordings at the onset of epileptic seizures. In this sense, neural networks are important for controlling chaotic dynamical systems. If RNNs have periodic orbits, then such periodic orbits are meant to capture the idea that certain activities or motions are learned by repetition [28]. It is known that time delay may affect dynamical behaviors of neural networks. The delayed axonal signal transmissions in neural networks make the dynamic behaviors more complicated and may destabilize stable equilibria, and give rise to periodic oscillation, bifurcation and chaos ([10]). Therefore, stability analysis of neural networks with constant or time-varying delays has been an attractive subject of research. Although delay differential equations present richer dynamics than simple low dimensional differential equations, qualitative analysis of delay differential equations is much more complex [4, 5, 15]. Bifurcation is one of the most important dynamical phenomena for the nonlinear neural systems [26]. Changing control parameter may cause instability, losing equilibrium points or occurring periodic solutions. It is determined locally whether a system has periodic solutions via Hopf bifurcation theory. In bifurcation theory literature, choosing the bifurcation parameter as delay parameter is common to see the effect of delay term. Neural networks with delays have very rich and complex dynamics. There are some other methods to investigate limit cycles or periodic orbits in RNNs [28]. Using Hopf bifurcation theorem is one of the widespread techniques. 2

Since determination of the roots of the characteristic polynomial corresponding to linearized system and then applying normal form theory are rather arduous, researchers have paid attention to study the dynamics of small-scale neurons rather than large-scale networks to understand these complicated structures. As pointed out by Ruan and Filfil [21], neural networks usually have a spatial extent due to the presence of a multitude of parallel pathways with a variety of axon sizes and lengths. Hence, there will be a distribution of conduction velocities along these pathways and a distribution of propagation delays. In these circumstances the signal propagation is not instantaneous and cannot be modeled with discrete delays. Therefore, a more appropriate way is to incorporate continuously distributed delays in finite or/and infinite time. Also, in some cases, the entire history affects the current state, so considering infinite time delay is more general [16, 24, 30]. Moreover, a distributed delay becomes a discrete delay when the delay kernel is a Dirac delta function at a certain time. We refer to [6, 10, 17, 21, 34] and the references therein for related work on networks with distributed delays. Olien and B´elair [19] investigated the stability of a two-neuron system with discrete time delays, that is, 

u1 (t) = −u1 (t) + a11 f (u1(t − τ1 )) + a12 f (u2 (t − τ2 )), 

u2 (t) = −u2 (t) + a21 f (u1(t − τ1 )) + a22 f (u2 (t − τ2 )).

(1)

They studied several cases, including τ1 = τ2 , a11 = a22 . They found that system (1) may undergo Hopf bifurcations at certain values of delay. Following this work, properties of two-neuron networks with delays have been studied intensively [6, 11, 21, 23, 29, 32, 34]. In the case of multiple delays, the dynamics of systems could be more complex and interesting since the characteristic equation is transcendental. Recently, Li and Hu [17] studied the following differential equations with multiple delays:  t  x1 (t) = −x1 (t) + a11 f ( F (t − s)x1 (s)ds) + a12 f (x2 (t − τ )), −∞ (2)  t  F (t − s)x2 (s)ds). x2 (t) = −x2 (t) + a21 f (x1 (t − τ )) + a22 f ( −∞

First, they investigated the stability of the zero equilibrium using RouthHurwitz criterion when delay term τ = 0. And then taking discrete delay τ 3

as bifurcation parameter, they showed the existence of local Hopf bifurcation using Hopf bifurcation theorem conditions. In this paper, extending the idea above, we take discrete delays between neurons that appears in Eq.(2) differently and consider the following general two-neuron network with discrete and distributed delays, that is,  t  x1 (t) = −x1 (t) + a11 f11 ( F (t − s)x1 (s)ds) + a12 f12 (x2 (t − τ2 )), −∞ (3)  t  F (t − s)x2 (s)ds), x2 (t) = −x2 (t) + a21 f21 (x1 (t − τ1 )) + a22 f22 ( −∞



i where xi (t) = dx , xi (t) represents the state of the ith neuron at time t and dt aij > 0 (i = 1, 2 and j = 1, 2) are real constants. Here, F (.) is nonnegative bounded delay kernel defined on [0, ∞) which reflects the influence of the past states on the current dynamics. System (3) is reduced to system (2) if fij = f , i = 1, 2 and j = 1, 2 and τ1 = τ2 . Similarly, it is reduced to system (1) when the delay kernel is taken as Dirac delta function and fij = tanh (see [17]). The architecture of system (3) is illustrated in Figure 1. One should underline that system (3) can maintain a periodic orbit such that these periodic orbits present periodic pattern and have been used in learning theory, which are meant to capture the idea that certain activities or motions are learned by repetition [28].

Figure 1: Architecture of the model (3). Two neurons send signals to each other with a discrete delay, τj , j = 1, 2. One element receives one delayed input from itself with distributed delay, that has shown with dashed line.

Our aim in this paper is to give a detailed Hopf bifurcation analysis of Eq.(3) with second distributed delay is omitted as the first case. A complete 4

Hopf bifurcation analysis of Eq.(3) will be studied in a separate paper later since the dimension of the system is getting higher and characteristic equation is more complex. Choosing τ = τ1 + τ2 as a bifurcation parameter we study the stability of the zero solution and investigate the local Hopf bifurcation properties. This paper is organized as follows. In Section 2, stability of the equilibrium and the existence of Hopf bifurcation are investigated. In Section 3, the direction of Hopf bifurcation and the stability and period of bifurcating periodic solutions on the center manifold are determined. Finally, in Section 4, we consider an example and simulate it using MATLAB to support our theoretical results. 2. Stability Analysis and Hopf Bifurcation As we mentioned it in Introduction, we first study a simplified version of Eq.(3) in this paper. This model that we consider has the following form  t  F (t − s)x1 (s)ds) + a12 f12 (x2 (t − τ2 )) x1 (t) = −x1 (t) + a11 f11 ( (4) −∞  x2 (t) = −x2 (t) + a21 f21 (x1 (t − τ1 )) + a22 f22 (x2 (t)). In order  ∞ not to affect the equilibrium values, we normalize the kernel such that 0 F (s)ds = 1. In this paper, we consider only the weak kernel, that is, F (s) = αe−αs , α > 0, where α reflects the mean delay of the weak kernel. Now, it is necessary to make the following assumptions: (H1) fij ∈ C 3 , fij (0) = 0, (i = 1, 2 and j = 1, 2), (H2) τ = τ1 + τ2 . For convenience, we define a new variable as follows  t x3 (t) = F (t − s)x1 (s)ds. −∞

Then by the linear chain trick technique, system (4) can be transformed into the following system 

x1 (t) = −x1 (t) + a11 f11 (x3 (t)) + a12 f12 (x2 (t − τ2 )), 

x2 (t) = −x2 (t) + a21 f21 (x1 (t − τ1 )) + a22 f22 (x2 (t)), 

x3 (t) = −αx3 (t) + αx1 (t). 5

(5)

By the hypothesis (H1), it is easy to see that the origin (0, 0, 0) is an equilibrium point of system (5). Letting u1 (t) = x1 (t − τ1 ), u2 (t) = x2 (t), u3 (t) = x3 (t − τ1 ) and using the hypothesis (H2), system (5) can be rewritten as the following equivalent system 

u1 (t) = −u1 (t) + a11 f11 (u3 (t)) + a12 f12 (u2 (t − τ )), 

u2 (t) = −u2 (t) + a21 f21 (u1 (t)) + a22 f22 (u2 (t)),

(6)



u3 (t) = −αu3 (t) + αu1 (t). Under the hypothesis (H1), the linearization of system (6) at (0, 0, 0) is 

u1 (t) = −u1 (t) + α11 u3 (t) + α12 u2 (t − τ ), 

u2 (t) = −u2 (t) + α21 u1 (t) + α22 u2 (t),

(7)



u3 (t) = −αu3 (t) + αu1 (t), 

where αij = aij fij (0), (i = 1, 2 and j = 1, 2). The corresponding characteristic equation of (7) is λ3 + a2 λ2 + a1 λ + a0 + (b1 λ + b0 )e−λτ = 0,

(8)

where a2 = α − α22 + 2 a1 = 2α − α(α11 + α22 ) − α22 + 1 a0 = α − α(α11 + α22 ) + αα11 α22

b1 = −α12 α21 b0 = −αα12 α21 .

(9)

If τ = 0, that is, when there is no discrete time delay, Eq.(8) will be λ3 + a2 λ2 + (a1 + b1 )λ + (a0 + b0 ) = 0.

(10)

Now, it is necessary to investigate the distribution of roots of Eq.(10) in order to determine the stability of the origin. Using Routh-Hurwitz criteria for n = 3, all of the roots of the polynomial (10) are negative or have negative real parts if and only if the following conditions hold (i) a2 > 0, (ii) a0 + b0 > 0, (iii) a2 (a1 + b1 ) > a0 + b0 . 6

Now, let us take τ = 0. We shall investigate the roots of the transcendental equation (8) that lie in the left half of the complex plane. Suppose that λ = iw be a root of the characteristic equation (8) with w > 0. Substituting this into (8) and separating real and imaginary parts yield the following equations − a2 w 2 + a0 + b1 w sin wτ + b0 cos wτ = 0, −w 3 + a1 w + b1 w cos wτ − b0 sin wτ = 0.

(11) (12)

By taking square of (11) and (12) and then adding them up, one obtains w 6 + pw 4 + qw 2 + r = 0,

(13)

where p = a22 − 2a1 , q = a21 − 2a0 a2 − b21 and r = a20 − b20 . Setting z = w 2 , one can obtain z 3 + pz 2 + qz + r = 0. (14) Let us denote Eq.(14) as h(z) = z 3 + pz 2 + qz + r.

(15)

First, suppose that r < 0. Since limz→∞ h(z) = ∞, Eq.(15) has at least one positive root. On the other hand, suppose now r > 0. From Eq.(15), we have dh(z) = 3z 2 + 2pz + q. (16) dz If Δ = p2 − 3q ≤ 0, then Eq.(16) has no real roots. So, the function h(z) is monotone increasing on [0, ∞). Thus, when r > 0 and Δ ≤ 0, Eq.(15) has no positive real root. If Δ > 0, one can easily obtain that Eq.(16) has two real roots as follows √ √ −p + Δ −p − Δ ∗ ∗ z1 = , z2 = . 3 3 Now, we can give the following lemma without proof which is given in [25, 34] before. Lemma 1. For Eq.(14), we have the followings: (i) If r < 0, Eq.(14) has at least one positive root. (ii) If r ≥ 0 and Δ = p2 − 3q ≤ 0, then Eq.(14) has no positive root. 7

(iii) If r ≥ 0√and Δ > 0, then Eq.(14) has positive roots if and only if z1∗ = −p+3 Δ > 0 and h(z1∗ ) ≤ 0. (see Lemma 1 in [34]) Suppose that Eq.(14) has positive roots. Without loss of generality, we can assume that it has three positive roots, denoted by z1 , z2 and z3 , respectively. Then, Eq.(13) has three positive roots as follows √ √ √ w1 = z1 , w2 = z2 and w3 = z3 . For k = 1, 2, 3, there exists a sequence {τkj |j = 1, 2, 3, ...} such that Eq.(13) holds. One can easily obtain    b0 (a2 w 2 − a0 ) + b1 (w 4 − a1 w 2 ) 1  (j) + 2πj . (17) τk = arccos wk b20 + b21 w 2 (j)

Thus, ±iwk is a pair of purely imaginary roots of Eq.(8) with τ = τk . (0) (0) Define τ0 = τk0 = min{τk |k = 1, 2, 3}, w0 = wk0 . Let λ(τ ) = α(τ ) + iw(τ ) be the root of Eq.(8) near τ = τ0 satisfying α(τ0 ) = 0, w(τ0 ) = w0 . Then we have the following transversality condition. Lemma 2. Suppose that zk = wk2 , and h (zk ) = 0, where h(z) is defined by (15). Then  d(Reλ(τ )) = 0 (18) (j) dτ τ =τ k



and [(d(Reλ(τ ))/dτ )]τ =τ (j) and h (zk ) have the same sign. (see Lemma 4 in k [34]) Now, we can give our main theorem, which is one of the most important auxiliary results of the present paper. Theorem 1. The followings hold i) If r ≥ 0 and Δ = p2 − 3q ≤ 0, then the equilibrium point (0, 0, 0) of system (6) is asymptotically stable for all τ ≥ 0. ii) If either r < 0 or r ≥ 0, Δ > 0, z1∗ > 0 and h(z1∗ ) ≤ 0, then the equilibrium point (0, 0, 0) of system (6) is asymptotically stable for τ ∈ [0, τ0 ), iii) If the conditions of (ii) are satisfied, and h (zk ) = 0, then system (6) (j) undergoes a Hopf bifurcation at origin when τ = τk (j = 0, 1, 2, ...).

8

3. Direction and Stability of Hopf Bifurcation In former section, we have shown system (6) undergoes a Hopf bifurcation (j) when τ passes through increasingly the critical values τk (j = 0, 1, 2, ...). In this section, we shall investigate the direction and stability of periodic solutions by using the normal form theory and center manifold reduction presented in [12]. For fixed k ∈ {1, 2, 3} and j ∈ {0, 1, 2, ...}, let us introduce the new (j) parameter μ = τ − τk . Normalizing the delay τ by the time scaling t → t/τ , Eq. (6) can be rewritten as ⎤ ⎤ ⎤ ⎡  ⎡ ⎡ x1 (t) x1 (t) x1 (t − 1) ⎣ x2 (t) ⎦ =(τ (j) + μ)A(τ ) ⎣ x2 (t) ⎦ + (τ (j) + μ)B(τ ) ⎣ x2 (t − 1) ⎦ x3 (t) x3 (t) x3 (t − 1) + (τ (j) + μ)f (x1 , x2 , x3 ),

(19)

where ⎡

⎤ ⎡ ⎤ 0 α12 0 −1 0 α11 A(τ ) = ⎣ α21 α22 − 1 0 ⎦, B(τ ) = ⎣ 0 0 0 ⎦, α 0 −α 0 0 0 ⎡ ⎤ f1 (j) ⎣ ⎦ f f (x1 , x2 , x3 ) = (τ + μ) 2 f3 where f1 = β11 x23 (t) + σ11 x33 (t) + β12 x22 (t − 1) + σ12 x32 (t − 1) + H.O.T., f2 = β21 x21 (t) + σ21 x31 (t) + β22 x22 (t) + σ22 x32 (t) + H.O.T., f3 = 0   in which βij = 12 aij fij (0) and σij = 16 aij fij (0) (i = 1, 2 and j = 1, 2) and H.O.T. denotes the higher order terms. Let u(t) = (x1 (t), x2 (t), x3 (t))T . The linearization of Eq.(19) around the origin is given by u (t) = (τ (j) + μ)A(τ )u(t) + (τ (j) + μ)B(τ )u(t − 1). For φ = (φ1 , φ2 , φ3 )T ∈ Ω = C([−1, 0], R3 ) we can define Lμ : Ω −→ R3 as follows Lμ (φ) = (τ (j) + μ)A(τ )φ(0) + (τ (j) + μ)B(τ )φ(−1).

(20)

Now, system (6) can be written as a functional differential equation in Ω as u (t) = Lμ (ut ) + f (μ, ut),

9

(21)

where ut (θ) = u(t + θ) = (x1 (t + θ), x2 (t + θ), x3 (t + θ))T and f : R × Ω → R3 where ⎤ ⎡ β11 φ23 (0) + σ11 φ33 (0) + β12 φ22 (−1) + σ12 φ32 (−1) + H.O.T f (μ, φ) = (τ (j) +μ) ⎣ β21 φ21 (0) + σ21 φ31 (0) + β22 φ22 (0) + σ22 φ32 (0) + H.O.T ⎦ . 0 (22) By Riesz Representation theorem, there exists a 3×3 matrix-valued function η(θ, μ), θ ∈ [−1, 0], whose elements are of bounded variation such that  0 Lμ φ = dη(θ, μ)φ(θ). (23) −1

In fact we can choose η(θ, μ) = (τ (j) + μ)Aδ(θ) + (τ (j) + μ)Bδ(θ + 1), where δ is the Dirac delta function. For φ ∈ Ω, define  dφ(θ) , θ ∈ [−1, 0), dθ  A(μ)φ = 0 dη(ξ, μ)φ(ξ), θ = 0, −1 and

 R(μ)φ =

0, θ ∈ [−1, 0), . f (μ, φ), θ = 0.

(24)

(25)

(26)

Then the functional differential equation (21) is equivalent to the following abstract differential equation .

ut = A(μ)ut + R(μ)ut , where ut (θ) = u(t + θ) for θ ∈ [−1, 0). For ψ ∈ C([0, 1], (R3 )∗ ), define  − dψ(s) , s ∈ (0, 1], ∗ ds  A ψ(s) = 0 T dη (ξ, 0)ψ(−ξ), s = 0, −1

(27)

(28)

and a bilinear form 0 θ ψ(s), φ(θ) = ψ(0)φ(0) −

ψ(ξ − θ)dη(θ)φ(ξ)dξ, −1 ξ=0

10

(29)

where η(θ) = η(θ, 0). Then, A(0) and A∗ (0) are adjoint operators. Suppose that q(θ) and q ∗ (s) are eigenvectors of A(0) and A∗ (0) corresponding to λ = iω0 and λ = −iω0 , respectively. Let us take q(θ) = [1, q1 , q2 ]T eiω0 θ and q ∗ (s) = D1 [1, q1∗, q2∗ ]T eiω0 s . Using A(0)q(θ) = iω0 q(θ) and A∗ (0)q ∗ (θ) = −iω0 q ∗ (θ), one can easily obtain −τ (j) α21 , τ (j) (α22 − 1) − iω0 τ (j) α , q2 = (j) τ α + iω0

−τ (j) α12 eiω0 , τ (j) (α22 − 1) + iω0 τ (j) α11 . q2∗ = (j) τ α − iω0 q1∗ =

q1 =

Furthermore, using the relation q ∗ (s), q(θ) = 1 one can calculate D as follows D = 1 + q1∗ q1 + q2∗ q2 + τ (j) e−iω0 α12 q1 . In the remainder of this section, we use the same notations as in [12]. We first compute the coordinates to describe the center manifold C0 at μ = 0. Let ut be the solution of Eq.(27) with μ = 0. Define z(t) = q ∗ , ut ,

w(t, θ) = ut − 2Re{z(t)q(θ)}.

(30)

On the center manifold, we have w(t, θ) = w(z(t), z(t), θ) = w20 (θ)

z2 z2 + w11 (θ)zz + w02 (θ) + · · · , 2 2

(31)

where z and z are local coordinates for the center manifold C0 in the direction of q and q ∗ . For ut ∈ C0 we have .

.

z(t) = q ∗ , ut = q ∗ , Aut + Rut

= iω0 q ∗ , ut + q ∗ (0)f0 (z, z) ≡ iω0 z(t) + g(z, z), where g(z, z) = q ∗ (0)f0 (z, z) = g20

z2 z2 z2z + g11 zz + g02 + g21 + ···. 2 2 2

(32)

Here, f0 (z, z) denotes f (z, z) at μ = 0. Notice that ut (u1t (θ), u2t (θ), u3t (θ)) =

11

w(t, θ) + zq(θ) + zq(θ) and q(θ) = [1, q1 , q2 ]T eiω0 θ , so we have z2 z2 (1) (1) + w11 (0)zz + w02 (0) + O(|z, z|3 ), 2 2 2 z2 z (2) (2) (2) u2t (0) = zq1 + zq1 + w20 (0) + w11 (0)zz + w02 (0) + O(|z, z|3 ), 2 2 2 2 z z (3) (3) (3) u3t (0) = zq2 + zq2 + w20 (0) + w11 (0)zz + w02 (0) + O(|z, z|3 ), 2 2 2 z (2) (2) u2t (−1) = zq1 e−iω0 + zq1 e−iω0 + w20 (−1) + w11 (−1)zz 2 2 z (2) +w02 (−1) + O(|z, z|3 ). 2 (1)

u1t (0) = z + z + w20 (0)

Thus, it follows from the definition of f (μ, φ) and (32) that ⎡ 0 ⎤ f1 (j) τ ∗ ∗ ⎣ ∗ [1, q1 , q2 ] f20 ⎦ , g(z, z) = q (0)f0 (z, z) = D f0 3

where f10 = β11 u23t (0) + σ11 u33t (0) + β12 u22t (−1) + σ12 u32t (−1), f20 = β21 u21t (0) + σ21 u31t (0) + β22 u22t (0) + σ22 u32t (0), f30 = 0. Thus, we have g(z, z) = q ∗ (0)f0 (z, z)  τ (j)  β11 q22 + β12 q12 e−2iω0 + β21 q1∗ + β22 q1∗ q12 z 2 =  D + 2β11 q2 q2 + 2β12 q1 q1 e−2iω0 + 2β21 q1∗ + 2β22 q1∗ q1 q1 zz   + β11 q2 2 + β12 q1 2 e−2iω0 + β21 q1∗ + β22 q1∗ q1 2 z 2  3 3 + 2β11 q2 w11 (0) + β11 q2 w20 (0) + 3σ11 q22 q2 2 2 + 2β12 q1 e−iω0 w11 (−1) + β12 q1 e−iω0 w20 (−1) + 3σ12 q12 q1 e−3iω0 1 1 + 2β21 q1∗ w11 (0) + β21 q1∗ w20 (0) + 3σ21 q1∗

  2 2 + 2β22 q1 q1∗ w11 (0) + β22 q1 q1∗ w20 (0) + 3σ22 q1 q1∗ q12 z 2 z + H.O.T. 12

Comparing the coefficients in (32) one gets the coefficients as follows:  τ (j)  β11 q22 + β12 q12 e−2iω0 + β21 q1∗ + β22 q1∗ q12 , D (j)   τ g11 = 2β11 q2 q2 + 2β12 q1 q1 e−2iω0 + 2β21 q1∗ + 2β22 q1∗ q1 q1 , D  τ (j)  β11 q2 2 + β12 q1 2 e−2iω0 + β21 q1∗ + β22 q1∗ q1 2 , g02 = 2 D τ (j)  3 3 2β11 q2 w11 (0) + β11 q2 w20 (0) + 3σ11 q22 q2 g21 = 2 D 2 2 + 2β12 q1 e−iω0 w11 (−1) + β12 q1 e−iω0 w20 (−1) + 3σ12 q12 q1 e−3iω0 g20 = 2

1 1 + 2β21 q1∗ w11 (0) + β21 q1∗ w20 (0) + 3σ21 q1∗

 2 2 + 2β22 q1 q1∗ w11 (0) + β22 q1 q1∗ w20 (0) + 3σ22 q1 q1∗ q12 . In order to determine g21 , we need to compute w11 (θ) and w20 (θ). From (30) we can write .

.

.

w(t, θ) = xt − 2Re{z(t)q(θ)}  θ ∈ [−1, 0) Aw − 2Re {q ∗ (0)f0 q(θ)} , = ∗ Aw − 2Re {q (0)f0 q(θ)} + f0 , θ=0 : ≡ Aw + H(z, z, θ), where H(z, z, θ) = H20 On the other hand, one has .

z2 z2 + H11 zz + H02 + · · · . 2 2

(33)

.

.

w = wz z + wz z on the center manifold. Thus, comparing the coefficients one obtains that (A − 2iω0 ) w20 (θ) = −H20 (θ), Aw11 (θ) = −H11 (θ). For θ ∈ [−1, 0),

(34)

.

H(z, z, θ) = −2Re{z(t)q(θ)}.

Comparing the coefficients of (34) with those of (33) we obtain the following equalities H20 (θ) = − (q(θ)g20 + q(θ)g 02 ) , H11 (θ) = − (q(θ)g11 + q(θ)g 11 ) , H02 (θ) = − (q(θ)g02 + q(θ)g 20 ) . 13

From (34) and the definition of the A, we get 

w20 (θ) − 2iω0 w20 (θ) = q(θ)g20 + q(θ)g 02 . Then, since q(θ) = q(0)eiω0 θ , we have w20 (θ) = ⎡

(1)

i i g20 q(0)eiω0 θ + g q(0)e−iω0 θ + E1 e2iω0 θ , ω0 3ω0 02 ⎤

E1 ⎢ (2) ⎥ where E1 = ⎣ E1 ⎦ ∈ R3 is a constant vector. Similarly, (3) E1 w11 (θ) = ⎡ ⎢ where E2 = ⎣

(1) E2 (2) E2 (3) E2

⎤

−i i g11 q(0)eiω0 θ + g 11 q(0)e−iω0 θ + E2 , ω0 ω0

⎥ ⎦ ∈ R3 is a constant vector. Let us find the values of E1

and E2 . If we take θ = 0 at (34), then one obtains that 0 dη(θ)w20 (θ) = 2iω0 w20 (0) − H20 (0),

(35)

dη(θ)w11 (θ) = −H11 (0).

(36)

−1

0 −1

Also, for θ = 0 ⎤ β11 q22 + β12 q12 e−2iω0 ⎦, β21 + β22 q12 H20 (0) = −g20 q(0) − g 02 q(0) + 2τ (j) ⎣ 0

(37)

⎤ 2β11 q2 q2 + 2β12 q1 q1 e−2iω0 ⎦. 2β21 + 2β22 q1 q1 H11 (0) = −g11 q(0) − g 11 q(0) + τ (j) ⎣ 0

(38)

⎡

and

⎡

14

On the other hand, since A(0)q(0) = iw0 q(0) and A(0)q(0) = iw0 q(0) we can write ⎡ ⎤ 0 ⎣iw0 I − eiw0 θ dη(θ)⎦ q(0) = 0, (39) ⎡

−1

⎣−iw0 I −

0

⎤ e−iw0 θ dη(θ)⎦ q(0) = 0.

(40)

−1

Substituting (37) into (35) and using (39) we obtain ⎡ ⎤ ⎤ ⎡ 0 β11 q22 + β12 q12 e−2iω0 ⎣2iw0 I − e2iw0 θ dη(θ)⎦ E1 = 2τ (j) ⎣ ⎦ β21 + β22 q12 0 −1 which is equal to ⎤ ⎡ ⎡ ⎤ −τ (j) α12 −τ (j) α11 2iw0 + τ (j) β11 q22 + β12 q12 e−2iω0 ⎦. ⎣ −τ (j) α21 2iw0 − τ (j) (α22 − 1) ⎦ E1 = 2τ (j) ⎣ β21 + β22 q12 0 (j) (j) 0 0 2iw0 + τ α −τ α Now, if one solves this system for E1 one obtains  2 2 −2iω0  −τ (j) α12 −τ (j) α11 (j)  β11 q2 + β12 q1 e 2τ  (1) 2 (j) β21 + β22 q1 2iw0 − τ (α22 − 1) 0 E1 = A1  0 0 2iw0 + τ (j) α

(2)

E1 =

(3)

E1 = where

2τ A1

(j)

2τ A1

(j)

  2iw0 + τ (j) β11 q22 + β12 q12 e−2iω0 −τ (j) α11   −τ (j) α21 β21 + β22 q12 0  (j)  −τ α 0 2iw0 + τ (j) α

   ,  

  2iw0 + τ (j) −τ (j) α12 β11 q22 + β12 q12 e−2iω0   −τ (j) α21 2iw0 − τ (j) (α22 − 1) β21 + β22 q12   −τ (j) α 0 0

  2iw0 + τ (j) −τ (j) α12 −τ (j) α11  (j) (j) 0 A1 =  −τ α21 2iw0 − τ (α22 − 1)  −τ (j) α 0 2iw0 + τ (j) α 15

   .  

   ,  

   ,  

Similarly, substituting (38) into (36) and utilizing (40) we can easily get ⎤ ⎤ ⎡ ⎡ α11 −2β11 q2 q2 − 2β12 q1 q1 e−2iω0 −1 α12 ⎦. ⎣ α21 α22 − 1 0 ⎦ E2 = ⎣ −2β21 − 2β22 q1 q1 α 0 −α 0 Solving this system, we have    −2β11 q2 q2 − 2β12 q1 q1 e−2iω0 α12 α11   1  (1) −2β21 − 2β22 q1 q1 α22 − 1 0  , E2 = A2  0 0 −α     −1 −2β11 q2 q2 − 2β12 q1 q1 e−2iω0 α11    1  (2) , α −2β − 2β q q 0 E2 = 21 21 22 1 1  A2  α 0 −α    −2iω0   −1 α −2β q q − 2β q q e 12 11 2 2 12 1 1  1  (3) , α α − 1 −2β − 2β q q E2 = 21 22 21 22 1 1  A2   α 0 0   −1 α12 α11  A2 =  α21 α22 − 1 0  α 0 −α

where

   .  

Finally, we substitute E1 and E2 into w11 (θ) and w20 (θ) and find the coefficients of g(z, z) to determine the following formulae to investigate the properties of bifurcating periodic solution on the center manifold at the crit(j) ical value τk . The formulae have the following forms c1 (0) =

i  |g02 |2  g21 + , g11 g20 − 2 |g11 |2 − 2w0 3 2 μ2 = −

Re{c1 (0)} (j)

Re{λ (τk )}

,

β2 = 2Re{c1 (0)}, (j)

T2 = −

Im{c1 (0)} + μ2 Im{λ (τk )} (j)

τk w0

.

These are the quantities that determine the properties of bifurcating periodic (j) solutions on the center manifold at τk . Here, μ2 determines the direction 16

of Hopf bifurcation, β2 determines the stability of the bifurcating periodic solution and T2 determines the period of the bifurcating solution. Hence, we have the following result. Theorem 2. μ2 determines the direction of Hopf bifurcation; • If μ2 > 0, then the Hopf bifurcation is supercritical and the bifurcating periodic solutions exist for τ > τ0 . • If μ2 < 0, then the Hopf bifurcation is subcritical and the bifurcating periodic solutions exist for τ < τ0 . β2 determines the stability of the bifurcating periodic solution; • If β2 < 0, bifurcating periodic solutions are stable. • If β2 > 0, bifurcating periodic solutions are unstable. T2 determines the period of the bifurcating solution; • If T2 > 0, the period increases. • If T2 < 0, the period decreases. 4. Numerical Simulations In this section, we present some numerical simulations to support our results in Lemma 1-2 and Theorem 1. As an example, we simulate system (6) with a11 = −0.5, a12 = −1.8, a21 = 1.3, a22 = 1.7 and α = 1, i.e., we consider the following system 

u1 (t) = −u1 (t) − 0.5tanh(u3 (t)) − 1.8tanh(u2 (t − τ )), 

u2 (t) = −u2 (t) + 1.3tanh(u1 (t)) + 1.7tanh(u2 (t)),

(41)



u3 (t) = −u3 (t) + u1 (t). For simplicity, we take fij (.) = tanh(.) for i = 1, 2 and j = 1, 2. Then, we have w 6 + (1.49)w 4 − (2.7356)w 2 − 4.3731 = 0.

17

(42)

1

0.4

0.8

0.2

0.6

0

0.4 u2

u1

0.6

−0.2

0.2

−0.4

0

−0.6

−0.2

−0.8

−0.4

−1

0

100

300

200

400

−0.6

500

0

100

200

300

400

500

Time t

Time t

0.6

0.4

0.6 0.4

0.2

0.2

u3(t)

u3

0

−0.2

0 −0.2 −0.4

−0.4

−0.6 −0.6 −0.8 −1 −0.8

0

100

200

300

400

−1 −0.5

0

500

Time t

0 0.5

1

1

u1(t)

u2(t)

Figure 2: Graph of solution of system (6) with τ1 = 0.1, τ2 = 0.15, τ1 + τ2 = 0.25 < τ0 . The origin is asymptotically stable.

Eq.(42) has only one positive root, that is, w0 ≈ 1.2969. Also, one can easily obtain τ0 ≈ 0.2692 from (17). Using the results in Section 3, we can compute the values μ2 = 0.1030,

β2 = −0.2172,

T2 = 0.2687.

(43)

Choose τ1 = 0.1 and τ2 = 0.15, then τ = τ1 + τ2 = 0.25 < τ0 . Thus, from Theorem 1, the equilibrium (0, 0, 0) is stable when τ < τ0 which is supported in Figure 2. Since μ2 > 0, when τ passes through the critical value τ0 , the equilibrium loses its stability and a Hopf bifurcation occurs, i.e., a family of periodic solutions bifurcates from the origin. The periodic solutions are depicted in Fig 3 when τ = τ1 + τ2 = 0.3 ≈ τ0 . Since T2 > 0 and β2 < 0, the period of the periodic solutions increases as τ increases and periodic orbits are stable. If we choose τ1 = 0.2 and τ2 = 0.2, then τ = τ1 + τ2 = 0.4 > τ0 . When τ = 0.4 > τ0 in Figure 4, the corresponding periodic solutions have larger period than in Figure 3.

18

0.8

1

0.6

0.8

0.4

0.6

0.2 0.4 u2

u1

0 0.2

−0.2 0 −0.4 −0.2

−0.6

−0.4

−0.8 −1

0

100

200

300

400

−0.6

500

0

100

200

300

400

500

Time t

Time t

0.6

0.4 1 0.2

0.5

u3

u3(t)

0

−0.2

0

−0.5

−0.4

−1 1 0.5

−0.6

1 0.5

0 −0.8

0

−0.5 0

100

200

300

400

500

u2(t)

Time t

−0.5 −1

−1

u1(t)

Figure 3: The bifurcating periodic solution when τ1 = 0.1, τ2 = 0.2, τ1 + τ2 = 0.3 > τ0 .

5. Conclusion In this paper, we investigate local stability of the equilibrium (0, 0, 0) and local Hopf bifurcation in a two-neuron network model with multiple discrete and distributed delays. We show that the equilibrium is asymptotically stable when τ ∈ [0, τ0 ) and unstable for τ > τ0 . Also, we paid attention to the direction and the stability of the bifurcating periodic solutions by applying the normal form theory and the center manifold theorem. In particular, if we choose the kernel as delta function instead of weak kernel, that is, F (s) = δ(s − τi ), i = 1, 2, respectively, then system (3) reduces the model (1) that was studied in [19]. Specifically, employing the Routh-Hurwitz criterion and the results on distribution of the zeros of transcendental functions, we get a set of conditions to describe the stability of the fixed point of model (4) and the existence of Hopf bifurcations. Furthermore, we have shown the direction of the Hopf bifurcation and stability of the bifurcating periodic solutions by using the normal form method and center manifold theorem introduced by [12]. One should emphasize that in a real world problem, considering the neural net19

0.8

0.4

0.6

0.2

0.4

0

0.2 u2

1

0.6

u1

0.8

−0.2

0

−0.4

−0.2

−0.6

−0.4

−0.8 −1

−0.6

0

100

200

300

400

−0.8

500

0

100

200

Time t

300

400

500

Time t

0.6

0.4 0.5 0.2 0 u3

u3(t)

0

−0.5

−0.2

−0.4

−1 1 0.5

−0.6

1 0.5

0 −0.8

0

−0.5 0

100

200

300

400

500

u2(t)

Time t

−0.5 −1

−1

u1(t)

Figure 4: The bifurcating periodic solution when τ1 = 0.2, τ2 = 0.2, τ1 + τ2 = 0.4 > τ0 .

works with discrete and distributed delays is not only more general but also more realistic. Finally, we have performed some numerical simulations in order to support and extend the analysis results that we have obtained. References [1] F. M. Atay, Complex Time-Delay Systems: Theory and Applications, Springer, 2010. [2] K. L. Babcock, R. M. Westervelt, Dynamics of simple electronic neural networks, Physica D 28 (1987) 305–316. [3] P. Baldi, A. Atiya, How delays affect neural dynamics and learning, IEEE Trans. Neural Networks 5 (1994) 612–621. [4] B. Balachandran, T. K. Nagy, D. E. Gilsinn, Delay Differential Equations: Recent Advances and New Directions, Springer, 2009. [5] R. Bellman, K. L. Cooke, Differential-Difference Equations, Academic Press, 1963. 20

[6] P. Bi, Z. Hu, Hopf bifurcation and stability for a neural network model with mixed delays, Applied Mathematics and Computation 218 (2012) 6748–6761, doi:10.1016/j.amc.2011.12.042. [7] S. A. Campbell, Stability and bifurcation of a simple neural network with multiple time delays, Fields Inst Commun 21 (1999) 65–79. [8] J. Cao, M. Xiao, Stability and Hopf bifurcation in a simplified BAM neural network with two time delays, IEEE Transactions on Neural Networks 18 (2007) 416–430, doi:10.1109/TNN.2006.886358. [9] P. Das, A. Kundu, Bifurcation and chaos in delayed cellular neural network model, Journal of Applied Mathematics and Physics 2 (2014) 219– 224. [10] Y. Du, R. Xu, Q. Liu, Stability and bifurcation analysis for a discrete time bidirectional ring neural network model with delay, Electronic Journal of Differential Equations 2013 (2013) 1–12. [11] S. Guo, L. Huang, L. Wang, Linear stability and Hopf bifurcation in a two-neuron network with three delays, International Journal of Bifurcation and Chaos 14 (2004) 2799–2810. [12] N. D. Hassard, Y. H. Kazarinoff, Theory and Applications of Hopf Bifurcation, Cambridge University Press, 1981. [13] S. Haykin, Neural Networks: A comprehensive Foundations, 2nd Edition, Prentice Hall, 1999. [14] J. Hopfield, Neurons with graded response have collective computational properties like those of two-state neurons, Proc. Nat. Acad. Sci. 81 (1984) 3088–3092. [15] Y. Kuang, Delay Differential Equations with Application in Population Dynamics, Academic Press, 1993. [16] T. Li, S. Fei, Stability analysis of Cohen-Grossberg neural networks with time-varying and distributed delays, Neurocomputing 71 (2008) 10691081, doi: 10.1016/j.neucom.2007.09.006.

21

[17] X. Li, G. Hu, Stability and hopf bifurcation on a neuron network with discrete and distributed delays, Applied Mathematical Sciences 5 (2011) 2077–2084. [18] C. Marcus, R. Westervelt, Stability of analog neural networks with delay, Phys. Rev. A 39 (1989) 347–359. [19] L. Olien, J. B´elair, Bifurcations, stability, and monotonicity properties of a delayed neural network model, Physica D 102 (1997) 349–363. [20] R. Rojas, Neural Networks: A Systematic Introduction, Springer-Verlag, 1996. [21] S. Ruan, R. S. Filfil, Dynamics of a two-neuron system with discrete and distributed delays, Physica D 191 (2004) 323–342. [22] A. Ruiz, D. H. Owens, S. Townley, Existence, learning, and replication of periodic motions in recurrent neural networks, IEEE Trans Neural Networks 9 (1998) 651–661. [23] L. P. Shayer, S. A. Campbell, Stability, bifurcation and multistability in a system of two coupled neurons with multiple time delays, SIAM J. Appl. Math. 61 (2000) 673–700. [24] Q. Song, J. Cao, Stability analysis of cohen-grossberg neural networks with both time-varying and continuously distributed delays, Journal of Computational and Applied Mathematics 197 (2006) 188203, doi:10.1016/j.cam.2005.10.029. [25] Y. Song, S. Yuan, Bifurcation analysis in a predator-prey system with time delay, Nonlinear Analysis: Real World Applications 7 (2006) 265– 284, doi:10.1016/j.nonrwa.2005.03.002. [26] Z. Song, J. Xu, Stability switches and double hopf bifurcation in a twoneural network system with multiple delays, Cogn Neurodyn 7 (2013) 505–521, doi:10.1007/s11571–013–9254–0. [27] A. K. O. Tiba, A. F. R. Araujo, M. N. Rabelo, Hopf bifurcation in a chaotic associative memory, Neurocomputing 152 (2015) 109–120, doi:10.1016/j.neucom.2014.11.013.

22

[28] S. Townley, A. Ilchmann, M. G. Weib, W. Mcclements, A. C. Ruiz, D. H. Owens, D. Pratzel-Wolters, Existence and learning of oscillations in recurrent neural networks, IEEE Transactions on Neural Networks 11 (2000) 205–214, doi:10.1109/72.822523. [29] J. Wei, S. Ruan, Stability and bifurcation in a neural network model with two delays, Physica D 130 (1999) 255–272, doi:10.1016/S0167– 2789(99)00009–3. [30] H. Xiang, J. Cao, Almost periodic solutions of recurrent neural networks with continuously distributed delays, Nonlinear Analysis: Theory, Methods and Applications 71 (2009) 60976108, doi:10.1016/j.na.2009.05.079. [31] E. Yilmaz, Neural networks with piecewise constant argument and impact activation, Ph.D. thesis, Middle East Technical University (2011). [32] W. Yu, J. Cao, Stability and hopf bifurcation on a two-neuron system with time delay in the frequency domain, International Journal of Bifurcation and Chaos 17 (2007) 1355–1366. [33] H. Zhang, Y. Wang, D. Liu, A comprehensive review of stability analysis of continuous-time recurrent neural networks, IEEE Transactions on Neural Networks and Learning Systems 25 (2014) 1229–1262, doi:10.1109/TNNLS.2014.2317880. [34] X. Zhou, Y. Wu, Y. Li, X. Yao, Stability and Hopf bifurcation analysis on a two-neuron network with discrete and distributed delays, Chaos, Solitons & Fractals 40 (2009) 1493–1505, doi:10.1016/j.chaos.2007.09.034.

23