Stability, bifurcations and synchronization in a delayed neural network model of n -identical neurons

Available online at www.sciencedirect.com

ScienceDirect Mathematics and Computers in Simulation 121 (2016) 12–33 www.elsevier.com/locate/matcom

Original articles

Stability, bifurcations and synchronization in a delayed neural network model of n-identical neurons✩ Amitava Kundu a , Pritha Das a,∗ , A.B. Roy b a Department of Mathematics, Indian Institute of Engineering Science and Technology, Shibpur, Howrah 711103, India b Department of Mathematics, Jadavpur University, Kolkata 700032, India

Received 31 May 2014; received in revised form 15 July 2015; accepted 25 July 2015 Available online 12 August 2015

Abstract This paper deals with dynamic behaviors of Hopfield type neural network model of n(≥3) identical neurons with two timedelayed connections coupled in a star configuration. Delay dependent as well as independent local stability conditions about trivial equilibrium is found. Considering synaptic weight and time delay as parameters Hopf-bifurcation, steady-state bifurcation and equivariant steady state bifurcation criteria are given. The criterion for the global stability of the system is presented by constructing a suitable Lyapunov functional. Also conditions for absolute synchronization about the trivial equilibrium are also shown. Using normal form method and the center manifold theory the direction of the Hopf-bifurcation, stability and the properties of Hopfbifurcating periodic solutions are determined. Numerical simulations are presented to verify the analytical results. The effect of synaptic weight and delay on different types of oscillations, e.g. in-phase, phase-locking, standing wave and oscillation death, has been shown numerically. c 2015 International Association for Mathematics and Computers in Simulation (IMACS). Published by Elsevier B.V. All rights ⃝ reserved. Keywords: Multiple-delay; Local stability; Global stability; Bifurcation; Synchronization

1. Introduction Neural networks usually incorporate time delays in the signal transmissions among neurons because of the finite propagation velocity of action potentials and non negligible time of a signal from a neuron to the receiving site of post-synaptic neuron. Time lags slow down the rate of information transmission among neurons. They can change the stability of neural network states. Recent years have witnessed a growing interest in the dynamics of interacting neurons with delayed feedback since Hopfield [19]. Stability, multi-stability, periodic phenomenon, bifurcations and chaos of delayed Hopfield neural networks have been studied by many researchers [1,4–6,12,23–27] to understand the dynamic behaviors of the identical neurons. The stability and Hopf-bifurcation analysis with respect to total delay has been done in Goodwin feedback model system also [20,34]. ✩ In memory of our mentor Prof. A.B. Roy. ∗ Corresponding author. Tel.: +91 9831022438.

E-mail addresses: kundoo [email protected] (A. Kundu), [email protected] (P. Das). http://dx.doi.org/10.1016/j.matcom.2015.07.006 c 2015 International Association for Mathematics and Computers in Simulation (IMACS). Published by Elsevier B.V. All rights 0378-4754/⃝ reserved.

A. Kundu et al. / Mathematics and Computers in Simulation 121 (2016) 12–33

13

Realistic modeling of neural networks inevitably requires careful design and variation of the connection topology. It is observed that synaptic coupling can change through learning [3] and a wide range of different behaviors can be established by varying the coupling strength and structure of neural networks. These observations provide us the motivation of this study to regard connection topology as the bifurcation parameter while time lags are fixed. In this paper, the generalization of Hopfield type arbitrary neural network model for n identical neurons with time delayed connections and bidirectional interconnections using delay-differential equations (DDE) is studied. The connectivity pattern is homogeneous in the sense that each neuron receives the same number of connections from other neurons with identical values of synaptic weights. This type of connectivity involves all permutation on n nodes (Sn symmetric group). Especially, Ashwin and Swift [2] and Hadley et al. [17] have studied a system of coupled oscillators with Sn symmetry. The presence of structural symmetry often causes purely imaginary eigen values with multiplicity in characteristic equation. As a result, several branches of phase-locked periodic solutions may bifurcate simultaneously from the trivial solution at some critical values of delay and weight. As neural oscillations and their synchronization support neural communication within and between cortical areas of brain and help in normal brain functioning [13], we analyze various types of synchronized or desynchronized (i.e. symmetric or symmetrybreaking) oscillations with inhibitory self-connection and excitatory interconnection and also for both inhibitory connections in our system. Xu et al. [32] studied the stability and Hopf bifurcation in a class of high-dimension neural network involving the discrete and distributed delays by introducing some virtual neurons to the original system. Campbell et al. [5] and Wu et al. [29] studied the stability and synchronization aspect of the three neurons homogeneous bidirectional associative memory (BAM) network model with single and multiple delays respectively. Using symmetric bifurcation theory both of them studied global stability analysis and discussed various boundaries of stability regions of trivial solution. Yuan [33] also studied a 4-neuron homogeneous BAM network model with different self and interconnecting synaptic weights. The synchronization aspect and asynchronous periodic waves arising from Hopf and equivariant Hopf-bifurcations are also investigated [11,8,9,7,10,30,31]. The plan for the article is as follows. In Section 2, construction of mathematical model is given. In Section 3, the delay-independent, delay-dependent stability conditions and instability criteria of the trivial equilibrium of the linearized system are shown. In Section 4, Hopf-bifurcations with respect to synaptic weight and time delay as parameters are discussed. Equivariant Hopf and equivariant steady-state bifurcation with respect to delay as parameter are also discussed. In Section 5, we consider the global stability criteria of the system using Lyapunov functional method. The synchronization property of the model is discussed in Section 6. Regions of absolute synchronization of single equilibrium, delay dependent de-synchronization and multiple equilibriums are plotted numerically. In Section 7, the direction, stability and period of the Hopf-bifurcating periodic solutions are determined by means of the normal form method and the center manifold theory. Example with numerical simulations is discussed to verify the analytical results obtained so far in Section 8. Finally, some concluding remarks are given in Section 9. 2. Mathematical model We consider an artificial n-neuron network model of Hopfield type with two-way (bidirectional) time delayed connections between the neurons and itself by the delay differential equations coupled with star configuration: x˙ j = −κ x j (t) + α f [x j (t − τ1 )] + β{g[x j−1 (t − τ2 )] + g[x j+1 (t − τ2 )] + · · · + g[x j+n−3 (t − τ2 )] + g[x j+n−2 (t − τ2 )]}

(1)

with j (mod n); where x j (t) is the activation state of neuron j ( j = 1, 2, . . . , n) at time t and x0 (t) = xn (t). κ > 0 is the decay rate of neurons, α(̸=0) is synaptic weight of self connection, β is the weight of synaptic connections from one neuron to another. τ1 and τ2 are the signal transmission delays for self and nearest-neighbor connections respectively. Throughout this paper, we use activation functions f, g : R → R satisfying the following assumptions for any x = (x1 , x2 , x3 , . . . , xn )T ∈ R n : (H1) (H2) (H3) (H4)

f (0) = g(0) = f ′′ (0) = g ′′ (0) = 0, f ′ (0) = g ′ (0) = 1, x f ′′ (x) < 0 and xg ′′ (x) < 0 for all x ̸= 0, −∞ < limx→±∞ f (x), g(x) < ∞, f ′′′ (0) and g ′′′ (0) are not zero simultaneously.

14

A. Kundu et al. / Mathematics and Computers in Simulation 121 (2016) 12–33

3. Local stability analysis Let h = max{τ1 , τ2 }. A natural phase space is the Banach space C = C([−h, 0]; R n ) of continuous mapping from [−h, 0] → R n equipped with usual super-norm ∥φ∥ = sup−h≤θ ≤0 |φ(0)| for φ ∈ C([−h, 0]; R n ). To investigate the linear stability analysis of the system, we consider the linearized stability of the trivial solution (0, 0, 0, . . . , 0) of (1). It is assumed that system (1) is provided with initial conditions: x j (s) = φ j (s) ̸= 0,

s ∈ [−h, 0], ( j = 1, 2, . . . , n),

(2)

where φi : [−h, 0] → R is assumed to be continuous. Linearizing (2) about (0, 0, 0, . . . , 0) it gives x˙ j = {−κ x j (t) + α[x j (t − τ1 )] + β[x j−1 (t − τ2 )] + β[x j+1 (t − τ2 )] + · · · + β[x j+n−3 (t − τ2 )] + β[x j+n−2 (t − τ2 )]}.

(3)

The characteristic equation of (3) is −λτ1 ∆n−1 + βe−λτ2 ]n−1 [λ + κ − αe−λτ1 − (n − 1)βe−λτ2 ] = 0. 1 (λ)∆2 (λ) = [λ + κ − αe

(4)

Since ∆1 and ∆2 are continuous functions limλ→+∞ ∆i (λ) = +∞ (i = 1, 2), there exists at least one positive root in (4). It is well known that the trivial solution of (1) will be locally asymptotically stable if all the roots of the characteristic equation (4) have negative real parts and unstable if at least one root has a positive real part. We thus establish results about the local asymptotic stability of the trivial solution of (1) by analyzing the roots of (4). Clearly, λ is a zero of (4) if and only if ∆i (λ) = 0, i = 1, 2. The delay-independent condition for the local stability of the trivial equilibrium of the network in parametric region is given as follows. Theorem 1. If the parameters satisfy 0 < κ − |α| − (n − 1) |β| , τ1 ≥ 0 and τ2 ≥ 0, then all the roots of the characteristic equation (4) have negative real parts. Proof. Let λ = µ + iω be an arbitrary root of ∆2 (λ) = 0, then ∆2 (µ + iω) = R(µ, ω) + i I (µ, ω), where R(µ, ω) = µ + κ − αe−µτ1 cos(ωτ1 ) − (n − 1)βe−µτ2 cos(ωτ2 ), I (µ, ω) = ω + αe−µτ1 sin(ωτ1 ) + (n − 1)βe−µτ2 sin(ωτ2 ). Clearly, R(µ, ω) ≥ µ + κ − |α| e−µτ1 − (n − 1) |β| e−µτ2 = R1 (µ). Then R1 (0) = κ − |α| − (n − 1) |β| > 0 under the assumptions of the theorem. Furthermore, R ′ (µ) = 1 + |α| τ1 e−µτ1 + (n − 1) |β| τ2 e−µτ2 > 0. Hence, R ′ (µ) > 0 for all µ ≥ 0 and R(µ, ω) > 0 for all µ ≥ 0, ω ∈ R. But from the above discussion we can say that R(µ, ω) = 0 and I (µ, ω) = 0 implies that µ < 0. Similar conclusion can be drawn by taking first factor of (4) into consideration. Thus, all the roots of the characteristic equation (4) have negative real parts. This completes the proof. The conditions for delay-dependent linear stability of the trivial equilibrium with parametric region are given as follows: 1 , n > 2 and τ1 ≥ 0, then all the Theorem 2. If the parameters satisfy −κ < β < 0, α < |β| , 0 ≤ τ2 < − 2(n−1)β roots of the characteristic equation (4) have negative real parts.

Proof. To begin, let λ = µ + iω be an arbitrary root of second factor of (4), separating the real and imaginary parts, we obtain µ = −κ + (n − 1)βe−µτ2 cos(ωτ2 ) + αe−µτ1 cos(ωτ1 ),

(5)

ω = −(n − 1)βe

(6)

−µτ2

sin(ωτ2 ) − αe

−µτ1

sin(ωτ1 ).

We now assume that (5) and (6) have roots µ and ω, where ω ≥ 0. From (6), if µ ≥ 0 and using the conditions 1 imposed on α, we find that ω < −nβ. The condition 0 ≤ τ2 < − 2(n−1)β then implies that 0 ≤ ωτ2 < 1. Hence,

A. Kundu et al. / Mathematics and Computers in Simulation 121 (2016) 12–33

15

Fig. 1. Linear delay dependent stability region (shaded portion in black color) shown in (α, β, τ1 ) and (α, β, τ2 ) space corresponding to Eq. (3).

1/2 < cos(1) < cos(ωτ1 ) < sin(1) < 1. Isolating the last term in both (5) and (6), squaring and adding, we obtain the necessary condition (µ + κ)2 + ω2 − 2(n − 1)βe−µτ2 {(µ + κ) cos(ωτ2 ) − ω sin(ωτ2 )} + (n − 1)2 β 2 e−2µτ2 − α 2 e−2µτ1 = 0, (7) for a solution to exist. For fixed values of ω, τ1 and τ2 , let us suppose the left-hand side of (7) is equal to M(µ) and note that M(0) = κ 2 + ω2 − 2(n − 1)κβ cos(ωτ2 ) + 2(n − 1)βω sin(ωτ2 ) + (n − 1)2 β 2 − α 2 . 1 Since sin(ωτ1 ) < ωτ1 and τ2 < − 2(n−1)β , then ω2 + 2(n − 1)βω sin(ωτ2 ) ≥ ω2 {1 + 2(n − 1)βτ2 } > 0, yields M(0) > 0. Taking the derivative of M(µ) with respect to µ, we obtain dM = 2{α 2 τ1 e−2µτ1 − (n − 1)βτ2 sin(ωτ2 ) dµ + (µ + κ)[1 + (n − 1)βτ2 e−µτ2 cos(ωτ2 )] − (n − 1)βe−µτ2 [cos(ωτ2 ) + (n − 1)βτ2 e−µτ2 ]}.

(8)

Since β < 0, ω ≥ 0, τ1 ≥ 0, τ2 ≥ 0, µ ≥ 0 and sin(ωτ2 ) ≥ 0, the first two terms in the first line of expression (8) are nonnegative. We now consider the other two terms in turn. 1 and µ ≥ 0, we have 0 ≤ e−µτ2 ≤ 1. Combining this with cos(ωτ2 ) ≤ 1, we obtain I. From 0 ≤ τ2 < − 2(n−1)β   1 > 0. (µ + κ)[1 + (n − 1)βτ2 e−µτ2 cos(ωτ2 )] ≥ (µ + κ) 1 − 2

II. From β < 0,

1 < cos(1) < cos(ωτ2 ), τ2 < − 2(n−1)β , and 0 ≤ e−µτ2 ≤ 1,   1 −(n − 1)βe−µτ2 [cos(ωτ2 ) + (n − 1)βτ2 e−µτ2 ] > −(n − 1)β cos(1) − > 0. 2 1 2

M Thus, ddµ > 0 for µ ≥ 0. Since M(0) > 0, we conclude that M(µ) > 0 if µ ≥ 0. Thus if M(µ) = 0, then µ < 0; i.e., all roots of ∆2 (λ) have negative real part. It may be shown in a similar manner that all roots of ∆1 (λ) have negative real part. α 1 Theorem 3. If the parameters satisfy −κ < α < 0, |β| < − (n−1) , n > 2, 0 ≤ τ1 < − 2α and τ2 ≥ 0, then all the roots of the characteristic equation (4) have negative real parts.

Proof. The proof of the above theorem can be shown in a manner similar to that of [5]. Delay dependent local stability region are shown in Fig. 1. The conditions for the instability of the trivial equilibrium of the network are given as follows:

16

A. Kundu et al. / Mathematics and Computers in Simulation 121 (2016) 12–33

Theorem 4. If α > κ, then trivial equilibrium of (1) is unstable for all values of β, τ1 ≥ 0 and τ2 ≥ 0. If α = κ, then the trivial solution of (1) is unstable for all values of β ̸= 0, τ1 ≥ 0 and τ2 ≥ 0. Proof. From the characteristic equation (4) we can write ∆1 (λ) = (λ + κ − αe−λτ1 + βe−λτ2 ). Then, with α > κ and β ≤ 0, ∆1 (0) = (κ − α + β) < 0. And, for λ ∈ R, lim ∆1 (λ) = lim (λ + κ − αe−λτ1 + βe−λτ2 ) = +∞ for all β ≤ 0, τ1 ≥ 0 and τ2 ≥ 0.

λ→+∞

λ→+∞

Hence, as ∆1 : R → R is a continuous function, there exists a λ∗ > 0 such that ∆1 (λ∗ ) = 0 for any fixed values of τ1 ≥ 0, τ2 ≥ 0, β ≤ 0 and α > κ. Now consider ∆2 (λ) = {λ + κ − αe−λτ1 − (n − 1)βe−λτ2 }. For β ≥ 0, with α > κ, ∆2 (0) = {κ − α − (n − 1)β} < 0 and, for λ ∈ R, lim ∆2 (λ) = lim {λ + κ − αe−λτ1 − (n − 1)βe−λτ2 } = +∞.

λ→+∞

λ→+∞

Hence, as ∆2 : R → R is a continuous function, there exists a λ∗ > 0 such that ∆2 (λ∗ ) = 0 for any fixed values of τ1 ≥ 0, τ2 ≥ 0, β ≥ 0 and α > κ. When α = κ, similar arguments show that ∆1 (λ) has a positive real root for any fixed values of τ1 ≥ 0, τ2 ≥ 0 and β ≥ 0. Thus, the characteristic equation has a positive real root for all β and all τ1 ≥ 0, τ2 ≥ 0 when α > κ and for all β ̸= 0 and all τ1 ≥ 0, τ2 ≥ 0 when α = κ. Hence the trivial solution is unstable for these parameter values. Theorem 5. The trivial equilibrium is unstable for all τ1 , τ2 ≥ 0, if either κ < (α − β) or κ < {α + (n − 1)β}. Proof. At the trivial equilibrium point, we have in the characteristic equation ∆1 (0) = κ − α + β < 0

when κ < (α − β),

∆2 (0) = κ − α − (n − 1)β < 0

when κ < {α + (n − 1)β}.

Since ∆1 and ∆2 are continuous functions limλ→∞ ∆i = +∞ (i = 1, 2) there exists λi > 0, such that ∆i (λi ) > 0 (i = 1, 2) respectively. Therefore, there exists at least one positive real root in (4), implying that the trivial solution is unstable in the given regions of parameters κ, α and β. 4. Bifurcation analysis As the parameters are varied, stability may be lost if a real root of the characteristic equation (4) passing through zero or by a pair of complex conjugate roots passing through the imaginary axis. For our system (3), there are four such codimension one bifurcations: i. A standard steady state bifurcation, which can occur when the factor ∆2 (λ) of the characteristic equation (4) has a single zero root; ii. An equivariant steady state bifurcation, which can occur when the factor ∆1 (λ) of the characteristic equation (4) has a single zero root; iii. A standard Hopf-bifurcation, which can occur when the factor ∆2 (λ) has simple pair of pure imaginary eigen values; iv. An equivariant Hopf-bifurcation, which can occur when the factor ∆1 (λ) has a repeated two pairs of pure imaginary roots. So the standard steady state bifurcation occurs when β = when β = α − κ.

1 (n−1) (κ − α) and equivariant steady state bifurcation occurs

4.1. Hopf-bifurcation with respect to synaptic weight (β) as a parameter Now we are interested to analyze the bifurcation phenomenon considering synaptic weight β as bifurcation parameter.

17

A. Kundu et al. / Mathematics and Computers in Simulation 121 (2016) 12–33

Let at β = β0 , characteristic equation (4) has a simple pair of purely imaginary roots ±iw0 . That is ∆2 (iw0 ) = 0 in the second factor of (4). This yields the following equations  κ − α cos(w0 τ1 ) − (n − 1)β0 cos(w0 τ2 ) = 0, (9) w0 + α sin(w0 τ1 ) + (n − 1)β0 sin(w0 τ2 ) = 0. Squaring both sides and rearranging we get w02 = α 2 + (n − 1)2 β02 + 2(n − 1)αβ0 cos w0 (τ1 − τ2 ) − κ 2 .

(10)

Now differentiating both sides of characteristic equation with respect to β and substituting λ = iw0 for β = β0 , it is obtained as    (n − 1){(n − 1)β0 τ2 + cos(w0 τ2 ) + ατ1 cos w0 (τ1 − τ2 )} dλ  Re = , (11) dβ β=β0 k12 + k22 where 

k1 = 1 + ατ1 cos(w0 τ1 ) + (n − 1)βτ2 cos(w0 τ2 ), k2 = ατ1 sin(w0 τ1 ) + (n − 1)βτ2 cos(w0 τ2 ).

Hence, the usual transversality condition is met if and only if (n − 1)β0 τ2 + cos(w0 τ2 ) + ατ1 cos w0 (τ1 − τ2 ) ̸= 0. Therefore we have the following theorem: Theorem 6. For given τ1 , τ2 , κ, α if at β = β0 all the roots of characteristic equation (4) have negative real parts except the pair of purely imaginary ones ±iw0 , then system (3) undergoes a Hopf-bifurcation with respect to synaptic weight at β = β0 where β0 lies on the surface  α cos(w0 τ1 ) κ = α cos(w0 τ1 ) + (n − 1)β0 cos(w0 τ2 ) , β0 ̸= − , n ≥ 3. (12) w0 = −α sin(w0 τ1 ) − (n − 1)β0 sin(w0 τ2 ) (n − 1) cos(w0 τ2 ) 4.2. Hopf-Bifurcation with respect to time delay (τ2 ) as a parameter The characteristic equation has a simple pair of pure imaginary roots λ = iω for parameter values such that ∆2 (±iω) = 0. This occurs when α, β, τ1 , τ2 , ω satisfy  κ − α cos(ωτ1 ) = (n − 1)β cos(ωτ2 ), (13) ω + α sin(ωτ1 ) = −(n − 1)β sin(ωτ2 ). + 1,σ +1 − 2,σ +1 For fixed α and τ1 this occurs along the curves (β H , τ2 ) and (β H , τ2 ) where

1  2 κ + α 2 + ω2 + 2αω sin(ωτ1 ) − 2ακ cos(ωτ1 ), (n − 1)  2σ π     , κ − α cos(ωτ1 ) > 0, 1 −1 −ω − α sin(ωτ1 ) ω = tan +  ω κ − α cos(ωτ1 )  (2σ + 1)π , κ − α cos(ωτ1 ) < 0. ω  2σ π     , κ − α cos(ωτ1 ) < 0, 1 −1 −ω − α sin(ωτ1 ) ω = tan +  ω κ − α cos(ωτ1 )  (2σ + 1)π , κ − α cos(ωτ1 ) > 0. ω

± βH =±

τ21,σ +1

τ22,σ +1

Therefore, we can obtain the following result.

(14)

σ = 0, 1, 2 . . .

(15)

σ = 0, 1, 2 . . . .

(16)

18

A. Kundu et al. / Mathematics and Computers in Simulation 121 (2016) 12–33

± Lemma 1. If the parameters satisfy β = β H , τ2 = τ2u,σ +1 (u = 1, 2, σ = 0, 1, 2 . . .), then there exists a sequence of a simple pair of pure imaginary eigen values for any values of α and τ1 for the characteristic equation (4). ± From Lemma 1, we know that the system (3) may go through Hopf-bifurcation at β = β H , τ2 = τ2u,σ +1 (u = + 1, 2, σ = 0, 1, 2 . . .) (defined in (15) and (16)) if the transversality condition is hold. For example, when β = β H , j,1 Hopf-bifurcation takes place at τ2 = τ2 (defined in (15)) if    (n − 1)ωβ{k4 cos(ωτ2 ) − k3 sin(ωτ2 )} dλ  = ̸= 0 (17) Re dτ2 λ=iω k32 + k42

where k3 = 1 + ατ1 cos(ωτ1 ) + (n − 1)βτ2 cos(ωτ2 ), k4 = ατ1 sin(ωτ1 ) + (n − 1)βτ2 sin(ωτ2 ) with k3 , k4 not zero simultaneously. 4.3. Equivariant Hopf-bifurcation with respect to time delay (τ2 ) as a parameter As we notice that the characteristic equation (4) has repeated two pairs of pure imaginary roots λ = ±iω for parameter values such that ∆1 (±iω) = 0 (multiplicity of ±iω is four). This occurs when α, β, τ1 , τ2 , ω satisfy  κ − α cos(ωτ1 ) = −β cos(ωτ2 ), (18) ω + α sin(ωτ1 ) = β sin(ωτ2 ). 3,σ +1 4,σ +1 + − For fixed α and τ1 this occurs along the curves (β H ) and (β H ) where d , τ2 d , τ2  ± 2 2 2 βH d = ± κ + α + ω + 2αω sin(ωτ1 ) − 2ακ cos(ωτ1 ) ± = (n − 1)β H .  2σ π     , κ − α cos(ωτ1 ) < 0, 1 3,σ +1 −1 −ω − α sin(ωτ1 ) ω τ2 = tan σ = 0, 1, 2 . . . +  ω κ − α cos(ωτ1 )  (2σ + 1)π , κ − α cos(ωτ1 ) > 0. ω  2σ π     , κ − α cos(ωτ1 ) > 0, 1 4,σ +1 −1 −ω − α sin(ωτ1 ) ω τ2 = tan σ = 0, 1, 2 . . . . +  ω κ − α cos(ωτ1 )  (2σ + 1)π , κ − α cos(ωτ1 ) < 0. ω

(19)

(20)

(21)

Therefore, we can obtain the following result. ξ,σ +1

± Lemma 2. If the parameters satisfy β = β H (ξ = 3, 4, σ = 0, 1, 2 . . .), the system possess two pairs d , τ2 = τ2 of pure imaginary eigen values for any values of α and τ1 . ξ,n+1

± From Lemma 2, we know that the system (3) may go through equivariant Hopf-bifurcation at β = β H d , τ2 = τ2 + (ξ = 3, 4, σ = 0, 1, 2 . . .) (defined in (20) and (21)) if the transversality condition holds. For example, when β = β H d, ξ,1 Hopf-bifurcation takes place at τ2 = τ2 (defined in (21)) if    ωβ{k5 sin(ωτ2 ) − k6 cos(ωτ2 )} dλ  Re = ̸= 0, (22) dτ2 λ=iω k52 + k62

where k5 = 1 + ατ1 cos(ωτ1 ) − βτ2 cos(ωτ2 ), k6 = ατ1 sin(ωτ1 ) − βτ2 sin(ωτ2 ) with k5 , k6 not zero simultaneously.

A. Kundu et al. / Mathematics and Computers in Simulation 121 (2016) 12–33

19

5. Global stability analysis This section derives the criterion for the global asymptotic stability of the trivial equilibrium of the network by constructing a suitable Lyapunov functional. Theorem 7. If |α| + (n − 1) |β| < κ, then the trivial solution of (1) is globally, asymptotically stable for any τ1 ≥ 0, τ2 ≥ 0 and n > 2. Proof. We consider Lyapunov functional defined by n n  t n     2 2 V (t) = κ xm (t) + κ |α| f (xm (v))dv + κ |β| m=1 t−τ1

m=1  t

g (xm+1 (v))dv · · · + 2

+

t

t−τ2

g 2 (xm−1 (v))dv

t−τ2

m=1



t



g (xm+n−2 (v))dv . 2

t−τ2

Then, n n n    dV = 2κ xm (t)x˙m (t) + κ |α| [ f 2 (xm (t)) − f 2 (xm (t − τ1 ))] + κ |β| [g 2 (xm−1 (t)) dt m=1 m=1 m=1

− g 2 (xm−1 (t − τ2 )) + g 2 (xm+1 (t)) − g 2 (xm+1 (t − τ2 )) + · · · + g 2 (xm+n−2 (t)) − g 2 (xm+n−2 (t − τ2 ))], n  = 2κ xm {−κ xm (t) + α f [xm (t − τ1 )] + β{g[xm−1 (t − τ2 )] m=1

+ g[xm+1 (t − τ2 )] + g[xm+2 (t − τ2 )] + · · · + g[xm+n−2 (t − τ2 )]}} n n   + κ |α| [ f 2 (xm (t)) − f 2 (xm (t − τ1 ))] + κ |β| [g 2 (xm−1 (t)) − g 2 (xm−1 (t − τ2 )) m=1

m=1

+ g (xm+1 (t)) − g (xm+1 (t − τ2 )) + · · · + g (xm+n−2 (t)) − g 2 (xm+n−2 (t − τ2 ))] n n n n     ≤ −2κ 2 xm2 + κ |α| xm2 + (n − 1)κ |β| xm2 + κ |α| f 2 (xm (t)) 2

2

m=1 n 

+ κ |β|

2

m=1

m=1

m=1

[g 2 (xm−1 (t)) + g 2 (xm+1 (t)) + · · · + g 2 (xm+n−2 (t))].

m=1

From the hypothesis on f and g, we have the following expressions:  1 f (xm (t)) = pm (t)xm (t), pm (t) = f ′ (vxm (t))dv, g(xm (t)) = qm (t)xm (t),

qm (t) =



0 1

g ′ (vxm (t))dv and there exist p ∗ , q ∗ ∈ (0, 1] such that

0

pm (t) ≤ p ∗ ≤ 1, qm (t) ≤ q ∗ ≤ 1; (m = 1, 2, . . . , n).    Therefore nm=1 f 2 (xm (t)) ≤ nm=1 ( p ∗ )2 xm2 ≤ nm=1 xm2 . n n n n n 2 2 2 2 2 Similarly, m=1 g (x m±1 (t)) ≤ m=1 x m , m=1 g (x m+2 (t)) ≤ m=1 x m , . . . , m=1 g (x m+n−2 (t)) ≤ n 2 m=1 x m . Thus, n n n    dV ≤ (−2κ 2 + κ |α| + (n − 1)κ |β|) xm2 + κ |α| xm2 + (n − 1)κ |β| xm2 , dt m=1 m=1 m=1 ∴

n  dV ≤ −2κ(κ − |α| − (n − 1) |β|) xm2 . dt m=1

If |α| + (n − 1) |β| < κ,

dV dt

< 0 when x ̸= 0, hence the trivial solution of (1) is globally asymptotically stable.

20

A. Kundu et al. / Mathematics and Computers in Simulation 121 (2016) 12–33

6. Synchronization In this section apart from assumptions (H1)–(H4), we make additional assumption: (H5) 0 < g ′ (x) ≤ f ′ (x) ≤ 1. Our next focus, in the remaining part of this paper, is on the system (1), where j (mod n), f, g : R → R is C n -smooth and satisfies the conditions (H1)–(H5). A solution of (1) is asymptotically synchronous if the ω-limit set of the solution is contained in the set of synchronous phase points given by {Φ = (Φ1 , Φ2 , . . . , Φn )T ∈ C : Φ1 = Φ2 = · · · = Φn }.

(23)

If every solution of (1) is asymptotically synchronous for all fixed non-negative values of τ1 and τ2 , then the system is absolutely synchronous. Using the similar process as in [28], we can obtain the following result. Theorem 8. If |α| + |β| < κ, then the trivial solution of (1) is absolutely synchronous. Proof. Let τ1 , τ2 ≥ 0 be fixed and define h = max{τ1 , τ2 }. Consider a given solution x : [−h, ∞) → R n of (1) and let y(t) = x1 (t) − x2 (t). Then from (1) we find, for all t ≥ 0 y˙ (t) = −κ y(t) + α[ f (x1 (t − τ1 )) − f (x2 (t − τ1 ))] − β[g(x1 (t − τ2 )) − g(x2 (t − τ2 ))] = −κ y(t) + αp(t)y(t − τ1 ) − βq(t)y(t − τ2 ),

(24)

where 

1

p(t) =

f ′ (vx1 (t − τ1 ) + (1 − v)x2 (t − τ1 ))dv,

0

 q(t) =

1

(25) g (vx1 (t − τ1 ) + (1 − v)x2 (t − τ1 ))dv. ′

0

Due to the C 1 smoothness of f and g, the boundedness of x1 , x2 : [−h, ∞) → R, and the normalization and concavity conditions of f and g, we can find p ∗ , q ∗ ∈ (0, 1] such that p(t) ≤ p ∗ and q(t) ≤ q ∗ for all t ≥ 0.  t  t Let ν(t) = κ y 2 (t) + κ |α| p ∗ y 2 (θ )dθ + κ |β| q ∗ y 2 (θ )dθ, t ≥ 0, (26) t−τ1

∴

t−τ2

dν(t) = 2κ y(t)[−κ y(t) + αp(t)y(t − τ1 ) − βq(t)y(t − τ2 )] dt + κ |α| p ∗ [y 2 (t) − y 2 (t − τ1 )] + κ |β| q ∗ [y 2 (t) − y 2 (t − τ2 )] ≤ −2κ 2 y 2 (t) + κ |α| p ∗ [y 2 (t) + y 2 (t − τ1 )] + κ |β| q ∗ [y 2 (t) + y 2 (t − τ2 )] + κ |α| p ∗ [y 2 (t) − y 2 (t − τ1 )] + κ |β| q ∗ [y 2 (t) − y 2 (t − τ2 )] ≤ −2κ(κ − |α| − |β|)y 2 (t).

(27)

From which it follows that  t  y 2 (t) + |α| p ∗ y 2 (θ )dθ + |β| q ∗ t−τ1

≤ y 2 (0) + |α| p ∗

t

y 2 (θ )dθ + 2(κ − |α| − |β|)

0

−τ1

y 2 (θ )dθ + |β| q ∗



t

y 2 (θ )dθ

0

t−τ2





0

y 2 (θ )dθ.

−τ2

∞ Consequently, 0 y 2 (θ )dθ < ∞. Note that y is bounded on [−h, ∞) and thus y˙ is bounded on [0, ∞). This together with y ∈ L 2 ([0, ∞)), implies that limt→∞ y(t) = 0.

A. Kundu et al. / Mathematics and Computers in Simulation 121 (2016) 12–33

21

Therefore limt→∞ [x1 (t) − x2 (t)] = 0. Similarly we can show that lim [x2 (t) − x3 (t)] = lim [x3 (t) − x4 (t)] = · · · = lim [xn (t) − x1 (t)] = 0.

t→∞

t→∞

t→∞

The above result shows that if the sum of self-connection and nearest neighbor interaction weights satisfies |α| + |β| < κ, then the system (1) is absolutely synchronous. Besides, assuming f (x) = g(x), the dynamics of system (1) is completely characterized by z˙ (t) = −κz(t) + {α + (n − 1)β} f (z(t − h))

(28)

where z is the common component of a synchronous solution of system (1), h = max{τ1 , τ2 }. We shall say Ras = {(α, β) : |α| + |β| < κ} as the absolutely synchronous region. It is easy to utilize the same argument for Theorem 8 to show that every solution z of (28) converges to zero as t → ∞ if |α| + (n − 1) |β| < κ. Consequently we have the following result. Lemma 3. If |α| + |β| < κ and |α| + (n − 1) |β| < κ, then every solution of (1) converges to zero as t → ∞. In other words, if we write Rasse = {(α, β) : |α| + |β| < κ, |α| + (n − 1) |β| < κ}, then the region Rasse of (1) can store and recover information in the form of absolutely synchronous single stable equilibrium. We consider a region R01 outside of Rasse as R01 = {(α, β) : |α| + (n − 1) |β| > κ, α + (n − 1)β > κ}, then we have the following theorem. Theorem 9. If (α, β) ∈ R01 , the system (1) has at least three equilibriums: u + , 0, u − with u ± = (x± , x± , . . . , x± )T , 0 = (0, 0, . . . , 0)T . Moreover, the trivial equilibrium is unstable and the two non-trivial equilibrium u ± are asymptotically stable. Proof. The equilibrium solution u ± , 0, satisfies the equation κu = α f (u) + (n − 1)βg(u). Let F(u) = α f (u) + (n − 1)βg(u) − κu, then F(0) = 0, F ′ (0) = α f ′ (0) + (n − 1)βg ′ (0) − κ = α + (n − 1)β − κ > 0, so the trivial equilibrium point is unstable in R01 . If (α, β) ∈ R01 , from the assumptions given in (H1)–(H4), we have   F(u) f (u) g(u) lim = lim α + (n − 1)β − κ = α f ′ (0) + (n − 1)βg ′ (0) − κ = α + (n − 1)β − κ > 0. u u u→0+ u u→0+ Then there exists u 1 > 0 and close to zero such that F(u 1 ) > 0. Since f (u) g(u) limu→+∞ F(u) u = limu→+∞ [α u + (n − 1)β u − κ] = −κ < 0,there exists u + > 0, such that F(u + ) < 0 and F ′ (u + ) < 0. Similarly, there exists u − < 0, such that F(u − ) < 0. Therefore, if (α, β) ∈ R01 , there are at least three equilibriums u ± and 0. The nontrivial equilibriums u ± are asymptotically stable since F(u ± ) < 0. This completes the proof. In Fig. 2, regions of absolute synchronization single stable equilibrium (Rasse ), absolute synchronization (Ras ), multiple equilibrium (R01 ) and delay dependent de-synchronization (R02 ∪ R03 ∪ R04 ) are illustrated for system (1). We numerically simulated three tables taking points from R02 ∪ R03 ∪ R04 . 7. Direction and stability of the Hopf-bifurcation In the previous analysis, we obtained conditions for a Hopf-bifurcation and equivariant Hopf-bifurcation with respect to time-delay (τ2 ) as parameter. In this section, we will determine the direction, stability and period of Hopfbifurcating periodic solution for n = 5 using the normal form and center manifold theorem introduced by Hassard et al. [18]. By similar way we can present the analysis of direction and stability of equivariant Hopf-bifurcation also. Let C([−τ2 , 0], R 5 ) be the Banach space of continuous mapping from [−τ2 , 0] into R 5 equipped with the supremum norm ∥ϕ∥ = sup−τ2 ≤θ≤0 |ϕ(θ )| for ϕ ∈ C([−τ2 , 0], R 5 ). Without loss of generality, it can be assumed that τ2 ≥ τ1 . Eq. (1) for n = 5 can be rewritten as an operator differential equation in C([−τ2 , 0], R 5 ).

22

A. Kundu et al. / Mathematics and Computers in Simulation 121 (2016) 12–33

Fig. 2. A partition in (α, β) plane: absolutely synchronous single stable equilibrium region Rasse (parallelogram shaded in red), absolute synchronization region Ras (square region shaded in green), multiple equilibrium region R01 , delay dependent de-synchronization regions R02 ∪ R03 ∪ R04 . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

For simplicity, it is assumed that fl′′ (0) = 0, gl′′ (0) = 0, l = 1, 2, . . . , 5. Recasting (1) as the following form with truncated Taylor expansions, we get  x˙   1        x˙2     x˙ 3     x˙4         x˙5  

= −κ x1 (t) + αx1 (t − τ1 ) + βx2 (t − τ2 ) + βx3 (t − τ2 ) + βx4 (t − τ2 ) + βx5 (t − τ2 ) + α1 x13 (t − τ1 ) + β1 x23 (t − τ2 ) + β1 x33 (t − τ2 ) + β1 x43 (t − τ2 ) + β1 x53 (t − τ2 ), = −κ x2 (t) + βx1 (t − τ2 ) + αx2 (t − τ1 ) + βx3 (t − τ2 ) + βx4 (t − τ2 ) + βx5 (t − τ2 ) + β1 x23 (t − τ2 ) + α1 x23 (t − τ1 ) + β1 x33 (t − τ2 ) + β1 x43 (t − τ2 ) + β1 x53 (t − τ2 ), = −κ x3 (t) + βx1 (t − τ2 ) + βx2 (t − τ2 ) + αx3 (t − τ1 ) + βx4 (t − τ2 ) + βx5 (t − τ2 ) + β1 x13 (t − τ2 ) + β1 x23 (t − τ2 ) + α1 x33 (t − τ1 ) + β1 x43 (t − τ2 ) + β1 x53 (t − τ2 ), = −κ x4 (t) + βx1 (t − τ2 ) + βx2 (t − τ2 ) + β1 x3 (t − τ2 ) + αx4 (t − τ1 ) + β1 x5 (t − τ2 ) + β1 x13 (t − τ2 ) + β1 x23 (t − τ2 ) + β1 x33 (t − τ2 ) + α1 x43 (t − τ1 ) + β1 x53 (t − τ2 ), = −κ x5 (t) + βx1 (t − τ2 ) + βx2 (t − τ2 ) + βx3 (t − τ2 ) + βx4 (t − τ2 ) + αx5 (t − τ1 ) + β1 x13 (t − τ2 ) + β1 x23 (t − τ2 ) + β1 x33 (t − τ2 ) + β1 x43 (t − τ2 ) + α1 x53 (t − τ1 )

(29)

where α1 = α f ′′′ (0)/6, β1 = βg ′′′ (0)/6. For convenience, let τ2 = τc + µ, (µ ∈ R). Then µ = 0 is the bifurcation value for (1) when n = 5. By the Riesz representation theorem, there exists a matrix whose components are bounded variation functions η(θ, µ) 0 for θ ∈ [−τc , 0] → R 5 such that L µ ϕ = −τc dη(θ, µ)ϕ(θ ). In fact we can choose  M1 , θ = 0,  M2 δ(θ + τ1 ), θ ∈ [−τ1 , 0), η(θ, µ) =  −M3 δ(θ + τc ), θ ∈ [−τc , −τ1 ) where δ(θ ) is Dirac function. For ϕ(θ ) = [ϕ1 (θ ), ϕ2 (θ ), ϕ3 (θ ), ϕ4 (θ ), ϕ5 (θ )]T ∈ C([−τc , 0] , R 5 ), we define an operator L µ ϕ = M1 ϕ(0) + M2 ϕ(−τ1 ) + M3 ϕ(−τc ),

23

A. Kundu et al. / Mathematics and Computers in Simulation 121 (2016) 12–33

where −κ  0  M1 =  0  0 0

0 −κ 0 0 0

0 0 −κ 0 0

0 0 0 −κ 0

0  0   0 , 0  −κ

α 0  M2 =  0 0 0

0 α 0 0 0

0 0 α 0 0

0 0 0 α 0

0 0  0 , 0 α

0 β  M3 = β β β

β 0 β β β

So using the obtained results from Appendix, the following parameters can be calculated:     1 g21 i g20 g11 − 2 |g11 |2 − |g02 |2 + , µ2 = −Re {C1 (0)} /Re λ′ (0) , C1 (0) = 2ωc 3 2 β2 = 2 Re {C1 (0)} ,  . T2 = −[Im{C1 (0)} + µ2 Im{λ′ (0)}]/ωc , where λ′ (0) = λ′ (τc ) = λ′ (τ2 ) λ=iωc

β β 0 β β

β β β 0 β

β β  β . β 0 (30)

(31)

Now as shown in [18], we have the following theorem for system (29). Theorem 10. (i) The bifurcating periodic solution is supercritical (subcritical) if µ2 > 0 (µ2 < 0). (ii) The bifurcating periodic solutions are asymptotically stable (unstable) if β2 < 0(β2 > 0). (iii) The period of the bifurcating periodic solution increases (decreases) if T2 > 0(T2 < 0). 8. Numerical results and discussion In this section, the analytical results obtained above are verified by using the following numerical example given below. We used Matlab 7.10 for simulation. The special aspect in our model is that there are different time delays in the self connection and nearest-neighbor connections between neurons. Besides, in this section we assume n = 5 (see Fig. 3), f (x) = g(x) = tanh(x), the decay rate κ = 0.5 with inhibitory self-connections but excitatory and inhibitory interconnections. If n = 5, in (α, β) parametric plane, delay independent local stability region, delay dependent stability region, the region of instability, steady-state bifurcation (SB) and equivariant steady-state bifurcation (ESB) lines are illustrated in Fig. 4 Example. x˙1 = −κ x1 (t) + α tanh[x1 (t − τ1 )] + β tanh[x2 (t − τ2 )] + β tanh[x3 (t − τ2 )] + β tanh[x4 (t − τ2 )] + β tanh[x5 (t − τ2 )], x˙2 = −κ x2 (t) + β tanh[x1 (t − τ2 )] + α tanh[x2 (t − τ1 )] + β tanh[x3 (t − τ2 )] + β tanh[x4 (t − τ2 )] + β tanh[x5 (t − τ2 )], x˙3 = −κ x3 (t) + β tanh[x1 (t − τ2 )] + β tanh[x2 (t − τ2 )] + α tanh[x3 (t − τ1 )] + β tanh[x4 (t − τ2 )] + β tanh[x5 (t − τ2 )],

(32)

x˙4 = −κ x4 (t) + β tanh[x1 (t − τ2 )] + β tanh[x2 (t − τ2 )] + β tanh[x3 (t − τ2 )] + α tanh[x4 (t − τ1 )] + β tanh[x5 (t − τ2 )], x˙5 = −κ x5 (t) + β tanh[x1 (t − τ2 )] + β tanh[x2 (t − τ2 )] + β tanh[x3 (t − τ2 )] + β tanh[x4 (t − τ2 )] + α tanh[x5 (t − τ1 )]. The conditions for delay-dependent linear stability of the trivial equilibrium proved in Theorems 2 and 3 are verified using two sets of parameter values: (i) κ = 0.5, α = −0.3, β = 0.05, τ1 = 1, τ2 = 2; (ii) κ = 0.5, β = −0.3, α = 0.2, τ2 = 0.2, τ1 = 2 (see Fig. 5). To verify equivariant steady state bifurcation results, we choose parameter values α = −1.12, β = −1.62 such that α − β = 0.5. Equivariant steady state and phase locked oscillation (period-doubling) are observed in Fig. 6(a) and (b) when the system (32) has zero eigen value of multiplicity four.

24

A. Kundu et al. / Mathematics and Computers in Simulation 121 (2016) 12–33

Fig. 3. Architecture of the model of 5 neurons and 2 time delays.

Fig. 4. Steady-state bifurcation (SB) and equivariant steady-state bifurcation (ESB) lines are shown in (α, β) plane. The parallelogram S1 is the local stability region; outside this region, S2 is the delay dependent stability region (shaded in blue) and S3 is the region of instability (shaded in green) corresponding to Eq. (3). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 5. Time series showing asymptotic stable behavior of the system (32) satisfying delay dependent linear stability conditions.

Choosing τ1 = 3, τ2 = 1.5, κ = 0.5, α = −1.12, we observe that Theorem 6 is verified. Here, origin loses its stability and a Hopf-bifurcation occurs when β = β0 = 0.41 at w0 = −1.48. Bifurcation diagram with respect to synaptic weight β is illustrated in Fig. 7(a) which clearly indicates that for β > β0 equilibrium is asymptotically stable. Time series of the system (32) showing periodic behavior is depicted in Fig. 7(b) while phase portrait of the system is illustrated in Fig. 8(a) for τ1 = 3, τ2 = 1.5 and β = .08 < β0 . But when β = 0.6 > β0 the other parameters

A. Kundu et al. / Mathematics and Computers in Simulation 121 (2016) 12–33

25

Fig. 6. (a) Bifurcation diagram with respect to delay τ2 of the system (32) when α = −1.12, β = −1.62 and τ1 = 4; (b) Phase locked oscillation (period-doubling) of (32) occurs when α = −1.12, β = −1.62, τ1 = 4 and τ2 = 0.25.

Fig. 7. (a) Bifurcation diagram with respect to parameter β when α = −1.12, τ1 = 3, τ2 = 1.5; (b) Time series of the system (32) showing periodic behavior when α = −1.12, τ1 = 3, τ2 = 1.5 and β = 0.08.

Fig. 8. (a) Projection of phase spaces in (x2 , x3 , x5 ) when α = −1.12, τ1 = 3, τ2 = 1.5 and β = 0.08; (b) Time series of the system (32) showing oscillation death when α = −1.12, τ1 = 3, τ2 = 1.5 and β = 0.6.

same as before, origin is asymptotically stable (oscillation death) as shown in Fig. 8(b). Also choosing τ1 = 1.5, j,1 + κ = 0.5, α = −1.12 we observe that equivariant Hopf-bifurcation occurs when τ2 = τ2 = 0.38 for β = β H d = 0.3 and ω = 0.6. Similarly by choosing τ1 = 2, κ = 0.5, α = −1.12, equivariant Hopf-bifurcation occurs when j,1 − τ2 = τ2 = 0.77 for β = β H d = −0.4 and ω = 0.06. These occur when α, β, τ1 , τ2 and ω satisfy (18). Bifurcation diagrams with respect to delay τ2 are illustrated in Figs. 9 and 10 for different parametric values. Here, switching behavior from in-phase oscillation to death of oscillation and then back to it via mixed and death of oscillation are seen in Fig. 9 verifying Table 1 and Hopf-bifurcating periodical solution is subcritical [15].

26

A. Kundu et al. / Mathematics and Computers in Simulation 121 (2016) 12–33

Fig. 9. Bifurcation diagram with respect to parameter τ2 showing stability switch when κ = 0.5, α = −1.12, β = 0.3 and τ1 = 1.5.

Fig. 10. Bifurcation diagram with respect to parameter τ2 showing stability switch when κ = 0.5, α = −1.12, β = −0.4 and τ1 = 2. Table 1 Self-inhibitory, excitatory interconnection with varying delay (κ = 0.5, α = −1.12 and β = 0.3). Delay τ1

Delay τ2

Dynamical behavior

1.5

0.01–0.3 1.5–2.1 2.6–6

In-phase Mixed-type In-phase

2.5

0.01–1 1.1–1.4 1.5–2.6 2.7–3.7 3.8–6

In-phase Phase locked Mixed-type Phase locked In-phase

4

0.01–1.9 1.95–2.1 2.2–3.9 4–5.4 5.5–6

In-phase Phase locked Mixed-type Phase locked In-phase

Substituting the parameters into (31) gives µ2 = −0.144 < 0,

T2 = −1.96 < 0,

β2 = −0.08 < 0,

λ′ (0) = −0.561 − 0.051i.

(33)

Hence by Theorem 10, the Hopf-bifurcating periodic solution of (32) is subcritical as µ2 < 0, the bifurcating periodic solutions are asymptotically stable as β2 < 0 and the period of bifurcating periodic solution decreases. Switching behavior from mixed type oscillation to in-phase oscillation via oscillation death can also be seen in Fig. 10. To simulate synchronized and de-synchronized oscillations, we choose parameter values (α, β) ∈ R02 ∪ R03 and numerically calculated three tables.

A. Kundu et al. / Mathematics and Computers in Simulation 121 (2016) 12–33

27

Fig. 11. (a) Synchronous periodic solution (In-phase) when α = −1.12, β = 0.3, τ1 = 1.5 and τ2 = 3; (b) phase portrait in (x2 , x3 , x5 ) space with same parameter value of the system (32).

Fig. 12. Solution trajectory showing phase locked oscillation and phase diagram in (x1 , x2 , x3 ) space of the system (32) when α = −1.12, β = 0.3, τ1 = 2 and τ2 = 2.4.

For self-inhibitory (α < 0) and excitatory (β > 0) interconnection, in-phase synchronized oscillation can be observed for small value of inter-connection delay τ2 , independent of self-connection delay τ1 . Phase locked oscillation is observed for long range of τ1 starting from τ1 = 2.5 (see Table 1). When both connections are inhibitory (α < 0 and β < 0), for large value of τ2 , in-phase synchronized oscillation is always obtained. But de-synchronization occurs only for small values of τ2 , independent of τ1 (see Table 2). In-phase synchronized oscillation can be observed repetitively depending on synaptic weight β (excitatory) and τ2 , independent of τ1 (see Table 3). We observe that five identical neurons either in-phase (see Fig. 11) or show phase locked oscillation (periodic  oscillate  solutions satisfying x j (t) = x j+1 t ± r5p , t ∈ R, r ̸= 0 (mod 5), j (mod 5), where p > 0 is the period) or standing   wave oscillation (periodic solutions satisfying x5− j (t) = x j t − 2p , t ∈ R, r ̸= 0 (mod 5), j (mod 5), where p > 0 is the period) for self-inhibitory, excitatory interconnection for different values of time delays (see Figs. 12 and 13). We can notice the following different types of periodic oscillations (de-synchronized and synchronized) when all synaptic weights are inhibitory satisfying Table 2 (see Figs. 14 and 15). We can also observe that five neurons gradually separate into three groups (four in two pair and one single). Two of the neurons from each group move in same wave form with phase shifting half a period (anti-phase periodic solutions between two pair), while the fifth neuron moves with different amplitude and period from other two groups as depicted in Figs. 13(a) and 15(a) illustrating standing wave solution. Synchronous periodic solution (In-phase) can be observed when α = −1.12, β = −0.4, τ1 = 2, τ2 = 1 (same as Fig. 11). For verification the results of Table 3 we can see in-phase synchronized oscillation independent of selfconnection delay τ1 in Fig. 16.

28

A. Kundu et al. / Mathematics and Computers in Simulation 121 (2016) 12–33

Fig. 13. (a) Mixed-type (standing wave) oscillation occurs when α = −1.12, β = 0.3, τ1 = 1.5 and τ2 = 2; (b) phase portrait in (x2 , x3 , x5 ) space with same parameter value of the system (32). Table 2 Both self-inhibitory, inhibitory interconnection with varying delay (κ = 0.5, α = −1.12 and β = −0.4). Delay τ2

Delay τ1

Dynamical behavior

0.1

1.6–1.9 2–6

Mixed-type Phase locked

0.5

1.9–3.4 3.5–6

Mixed-type Phase locked

1

0.5–2.7 2.8–4 4.1–6

In-phase Mixed-type In-phase

2.5

0.4–6

In-phase

Table 3 Self-inhibitory, excitatory inter-connection with varying β (κ = 0.5, α = −1.12). Delay τ1

Delay τ2

Synaptic weight β

Dynamical behavior

1.5

1.5

0.3 0.4

Mixed-type Phase locked

1.5

2 3

0.3–0.4 0.2–0.4

Mixed-type In-phase

3

1 2,3 3.5,4.5 5

0.1–0.4

In-phase Mixed-type Phase locked In-phase

6

4 5 6,7 8

0.42–0.43

In-phase Mixed-type Phase locked In-phase

9. Conclusions This paper addresses both the global stability and local stability of artificial neural network model of n identical neurons with two delays. The results obtained enforce more restrictions on synaptic weight and delay parameters. It is observed that there is a region Ras in (α, β) parameter plane where every solution is absolutely synchronized. The synchronized stable pattern can be in the form of either trivial equilibrium (if |α| + (n − 1) |β| < κ) or multiple

A. Kundu et al. / Mathematics and Computers in Simulation 121 (2016) 12–33

29

Fig. 14. Solution trajectory showing phase locked oscillation of the inhibitory system and phase diagram in (x2 , x3 , x5 ) space of the system (32) when α = −1.12, β = −0.4, τ1 = 3 and τ2 = 0.1.

Fig. 15. (a) Mixed-type (standing wave) oscillation occurs when α = −1.12, β = −0.4, τ1 = 2, τ2 = 0.5; (b) phase portrait in (x2 , x3 , x5 ) space with same parameter value of the system (32).

Fig. 16. In-phase oscillation occurs when α = −1.12, β = 0.3, τ1 = 4 and τ2 = 1 of the system (32).

equilibrium (if |α| + (n − 1) |β| > κ). The coexistence of phase locked and synchronized periodic solutions suggest possible complicated structure due to Sn symmetry. Using numerical simulation of the system (32) we have obtained the following dynamic behavior for n = 5. We have mentioned standing wave pattern under mixed type oscillation, resulting from the presence of two or more waves of the same frequency with different directions which are usually seen in BAM neural networks. The behavior of the waves at the points of minimum and maximum vibrations contributes to the constructive interference [22].

30

A. Kundu et al. / Mathematics and Computers in Simulation 121 (2016) 12–33

Besides we observe some important features of the system such as scroll-like chaotic attractor near the equivariant Hopf-bifurcation with respect to delay τ2 in three dimensional phase portraits. There is significant evidence that rhythmic activity and synchronization (in-phase) of neuronal firing play a profound role in the information processing in many brain systems, including the olfactory bulb, hippocampus and thalamo-cortical system [21]. Altered oscillation or oscillation death associated with certain neuropsychiatric disorders involves dysfunction for many cognitive processes. To get prominent idea of objects, information needed to be linked up. This is viewed as ‘binding function’ of different information encoded in different brain areas. If there is no phase synchronization or if there is phase synchronization with an inappropriate phase lag, the action potentials will arrive at a phase in which sensitivity to neural input is non-optimal, and communication is likely to be blocked [14,16,28]. Generally both excitatory and inhibitory synaptic interaction may be responsible for synchronization of neuronal activity. Thus in-phase oscillation, phase-locked periodic oscillations, standing wave and oscillation death correlate neuronal processing and provide helpful measures for the assessment of normal and pathological brain functions. Appendix To describe the nonlinear terms in (29), let  α1 φ13 (−τ1 ) + β1 φ23 (−τc ) + β1 φ33 (−τc ) + β1 φ43 (−τc ) + β1 φ53 (−τc ) β φ 3 (−τ ) + α φ 3 (−τ ) + β φ 3 (−τ ) + β φ 3 (−τ ) + β φ 3 (−τ )  1 1 c c c c  1 2 1 1 3 1 4 1 5   3 3 3 3 3  G(µ, ϕ) =  β1 φ1 (−τc ) + β1 φ2 (−τc ) + α1 φ3 (−τ1 ) + β1 φ4 (−τc ) + β1 φ5 (−τc ) .   β1 φ13 (−τc ) + β1 φ23 (−τc ) + β1 φ33 (−τc ) + α1 φ43 (−τ1 ) + β1 φ53 (−τc ) 

(A.1)

β1 φ13 (−τc ) + β1 φ23 (−τc ) + β1 φ33 (−τc ) + β1 φ43 (−τc ) + α1 φ53 (−τ1 ) For ϕ ∈ C([−τc , 0] , R 5 ), we define two operators  dϕ(θ )   θ ∈ [−τc , 0),  dθ , A(µ)ϕ =  0   dη(θ, µ)ϕ(θ ), θ = 0. 

(A.2)

−τc

 R(µ)φ =

0, θ ∈ [−τc , 0), G(µ, ϕ), θ = 0.

(A.3)

In order to study Hopf-bifurcation problem, we transform system (29) into operator equation of the form: x˙t = A(µ)xt + R(µ)xt .

(A.4)

Eq. (A.4) is just the trivial equation d xt /dt = d xt /dθ when θ ∈ [−τc , 0), and it becomes (29) when θ = 0. Here, xt = x(t + θ ), x = (x1 , x2 , x3 , x4 , x5 )T . The adjoint operator A∗ of A for ψ ∈ C([0, τc ], R 5 ) is defined by  ∗ ∗ θ ∗ ∈ (0, τc ],  −dψ(θ )/dθ ,  0 ∗ A (µ)ψ = (A.5)  dη T (θ ∗ , µ)ψ(−θ ∗ ), θ ∗ = 0  −τc

η T (θ ∗ , µ)

where is the transpose of the matrix η. In order to normalize the eigen vectors of operator A and adjoint operator A∗ , the following bilinear form is needed to introduce: ⟨ψ, ϕ⟩ = ψ¯ T (0) · ϕ(0) − Here η(θ ) = η(θ, µ).



0

θ =−τc



θ

ξ =0

ψ¯ T (ξ − θ )dη(θ )ϕ(ξ )dξ.

(A.6)

A. Kundu et al. / Mathematics and Computers in Simulation 121 (2016) 12–33

31

For convenience in computation, it is beneficial to allow all functions to fall into C 5 instead of R 5 . We know that ±iωc are eigen values of A(0) and other eigen values have strictly negative real parts. Thus they are also eigen values of A∗ (0). Next, we calculate the eigen vector q(θ ) of A(0) belonging to the eigen value iωc and the eigen vector q ∗ (θ ) of A∗ (θ ) belonging to the eigen value (−iωc ). Then we obtain  A(0)q(0) = iωc q(θ ), for θ ∈ [−τc , 0), θ ∗ ∈ (0, τc ]. (A.7) A∗ (0)q ∗ (θ ) = −iωc q ∗ (θ ∗ ), Substituting (A.2) and (A.5) into (A.6) yields q(θ ) = q(0)eiωc θ ,

q ∗ (θ ∗ ) = q ∗ (0)eiωc θ , ∗

for θ ∈ [−τc , 0), θ ∗ ∈ (0, τc ]

(A.8)

and q(0) = [q1

q2

q3

q4

1]T ,

q ∗ (0) =

1 ∗ [q ρ 1

q2∗

q3∗

q4∗

1],

for θ = 0, θ ∗ = 0

(A.9)

where βe−iwc τc , iwc + k + αe−iwc τ1 − 3βe−iwc τc −βeiwc τc q1∗ = q2∗ = q3∗ = q4∗ = , −iwc + k + αeiwc τ1 + 3βeiwc τc q1 = q2 = q3 = q4 =

q(0) and q ∗ (0) have non zero solutions if and only if the following conditions hold: 0 0 ∗ det[iωc I − −τc dη(θ )eiωc θ ] = 0, det[iωc I + −τc dη T (θ ∗ )e−iωc θ ] = 0, where I is the identity matrix. The normalization conditions of q and q ∗ are  ∗    q , q = 1, q ∗ , q¯ = 0.

(A.10)

Substituting (A.8) and (A.9) into the first equation of (A.10) gives  0  θ  ∗  q , q = q¯ ∗ (0).q(0) − q¯ ∗T (ξ − θ )dη(θ )q(ξ )dξ, θ =−τc

ξ =0

1 = (q1 q¯1∗ + q2 q¯2∗ + q3 q¯3∗ + q4 q¯4∗ + 1) ρ¯  0  θ 1 ∗ ∗ ∗ ∗ [q1 , q2 , q3 , q4 , 1] × dη(θ )e−iωc θ [q1 , q2 , q3 , q4 , 1]T dξ, − ρ ¯ θ =−τc ξ =0 1 = (q1 q¯1∗ + q2 q¯2∗ + q3 q¯3∗ + q4 q¯4∗ + 1) + (ατ1 e−iωc τ1 + 3βe−iωc τc )(q¯1∗ + q¯2∗ + q¯3∗ + q¯4∗ ), ρ¯ = 1.

(A.11)

Then we can determine ρ in (A.9) from (A.11) as follows 4β 2 (iwc + k + αe−iwc τ1 − 3βe−iwc τ2 )(−iwc + k + αeiwc τ1 + 3βeiwc τ2 )   4βe−iwc τc − (4ατ1 e−iwc τ1 + 3βe−iwc τ2 ) . −iwc + k + αeiwc τ1 + 3βeiwc τ2

ρ¯ = 1 −

(A.12)

The second equation of (A.10) can be proved by simple calculation. Now the expressions of q(θ ) for θ ∈ [−τc , 0] and q ∗ (θ ∗ ) for θ ∗ ∈ [0, τc ] can be obtained. Next, we study the stability of bifurcating periodic solutions. Using Hassard et al. [18], the bifurcating periodic solutions Z (t, µ(ε)) has amplitude O(ε) and non-zero Floquet exponent β(ε) with β(0) = 0. Under the hypothesis,

32

A. Kundu et al. / Mathematics and Computers in Simulation 121 (2016) 12–33

µ, β are given by  µ = µ2 ε 2 + µ4 ε 4 + · · · , β = β2 ε 2 + β4 ε 4 + · · · . The sign of µ2 indicates the direction of bifurcation while that β2 determines the stability of Z (t, µ(ε)). Z (t, µ(ε)) is stable if β2 < 0 and unstable if β2 > 0. In the following, we will show how to derive the coefficients in this expressions, and we compute µ2 and β2 describing center manifold Ω0 at µ = 0. We define z(t) = ⟨q ∗ , xt ⟩, W (t, θ ) = xt − 2 Re [z(t)q(θ )] where xt is a solution of (A.4). On the center manifold Ω0 , we have W (t, θ ) = W (z, z¯ , θ ), = W20

z¯ 2 z2 + W11 (θ )z z¯ + W02 (θ ) + · · · . 2 2

(A.13)

Since µ = 0, z˙ (t) = ⟨q ∗ , x˙t ⟩, = ⟨q ∗ , Axt ⟩ + ⟨q ∗ , Rxt ⟩, = iwc z + q¯ ∗ (0)T · G(0, xt ). Let g(z, z¯ ) = g20

z2

(A.14) z¯ 2

z 2 z¯

+ g11 z z¯ + g02 + g21 + ···. (A.15) 2 2 2 In (A.14), G is an operator acting upon its argument taken as a function of θ, however, g is a function in z and z¯ , the dependence of θ is removed. From (A.1), (A.9) and (A.15) we may write   α1 x13 (−τ1 ) + β1 x23 (−τc ) + β1 x33 (−τc ) + β1 x43 (−τc ) + β1 x53 (−τc ) β x 3 (−τ ) + α x 3 (−τ ) + β x 3 (−τ ) + β x 3 (−τ ) + β x 3 (−τ )  1 1 c c c c  1 2 1 1 3 1 4 1 5    1 ∗ ∗ ∗ ∗ 3 3 3 3 3  q1 q2 q3 q4 1 β1 x1 (−τc ) + β1 x2 (−τc ) + α1 x3 (−τ1 ) + β1 x4 (−τc ) + β1 x5 (−τc ) g(z, z¯ ) = . ρ¯   3 3 3 3 3 β1 x1 (−τc ) + β1 x2 (−τc ) + β1 x3 (−τc ) + α1 x4 (−τ1 ) + β1 x5 (−τc ) β1 x13 (−τ1 ) + β1 x23 (−τc ) + β1 x33 (−τc ) + β1 x43 (−τc ) + α1 x53 (−τ1 ) (A.16) From the second equation in (A.1) and q(θ ) = [q1

q2

q3

xdt (−τ1 ) = W (γ ) (−τ1 ) + qγ ze−iωc τ1 + q¯γ z¯ eiωc τ1 , x1t (−τc ) = W

(1)

(−τc ) + (q2 + q3 + q4 + 1)ze

−iωc τc

q4

1]T eiωc θ , we obtain

γ = 1, . . . , 5, + (q¯2 + q¯3 + q¯4 + 1)¯z eiωc τc ,

x2t (−τc ) = W (2) (−τc ) + (q1 + q3 + q4 + 1)ze−iωc τc + (q¯1 + q¯3 + q¯4 + 1)¯z eiωc τc , x3t (−τc ) = W (3) (−τc ) + (q1 + q2 + q4 + 1)ze−iωc τc + (q¯1 + q¯2 + q¯4 + 1)¯z eiωc τc ,

(A.17)

x4t (−τc ) = W (4) (−τc ) + (q1 + q2 + q3 + 1)ze−iωc τc + (q¯1 + q¯2 + q¯3 + 1)¯z eiωc τc , x5t (−τc ) = W (5) (−τc ) + (q1 + q2 + q3 + q4 )ze−iωc τc + (q¯1 + q¯2 + q¯3 + q¯4 )¯z eiωc τc . where W (θ ) = [W 1 (θ ) W 2 (θ ) W 3 (θ ) W 4 (θ ) W 5 (θ )]T . By substituting (A.17) into (A.16), equating the result to the right hand side of (A.15) and making comparison of the coefficients of z 2 /2, z z¯ , z 2 z¯ /2 we get g20 = 0, g11 = 0, g02 = 0 and  1  ′′′  α f (0) 4 |q1 |4 + 1 e−iωc τ1 + βg ′′′ (0) g21 = ρ¯   × 4(3q¯1∗ + 1)(3q1 + 1)(9 |q1 |2 + 3(q1 + q¯1 ) + 1) + 256q¯1∗ |q1 |2 e−iωc τc . (A.18)

A. Kundu et al. / Mathematics and Computers in Simulation 121 (2016) 12–33

33

References [1] S. Arik, V. Tavsanoglu, Global asymptotic stability analysis of bidirectional associative memory neural networks with constant time delays, Neurocomputing 68 (2005) 161–176. [2] P. Ashwin, J.W. Swift, The dynamics of n weakly coupled identical oscillators, J. Nonlinear Sci. 2 (1992) 69–108. [3] F.M. Atay, J. Josty, A. Wende, Delays, connection topology, and synchronization of coupled chaotic maps, Phys. Rev. Lett. 92 (14) (2004) 144101. [4] J. B´elair, S.A. Campbell, P.V. Driessche, Frustration, stability, and delay-induced oscillations in a neural network model, SIAM J. Appl. Math. 56 (1996) 245–255. [5] S.A. Campbell, I. Ncube, J. Wu, Multi-stability and stable asynchronous periodic oscillations in a multiple delayed neural system, Physica D 214 (2006) 101–119. [6] J. Cao, Global stability conditions for delayed CNNs, IEEE Trans. Circuits Syst. I 48 (11) (2001) 1330–1333. [7] J. Cao, A. Alofi, A. Al-Mazrooei, A. Elaiw, Synchronization of switched interval networks and applications to chaotic neural networks, Abstr. Appl. Anal. 2013 (2013) 1–11. Article ID 940573. [8] J. Cao, G. Chen, P. Li, Global synchronization in an array of delayed neural networks with hybrid coupling, IEEE Trans. Syst. Man Cybern. B 38 (2) (2008) 488–498. [9] J. Cao, L. Li, Cluster synchronization in an array of hybrid coupled neural networks with delay, Neural Netw. 22 (2009) 335–342. [10] J. Cao, Y. Wang, Matrix measure strategies for stability and synchronization of inertial BAM neural network with time delays, Neural Netw. 53 (2014) 165–172. [11] J. Cao, M. Xiao, Stability and Hopf bifurcation in a simplified BAM neural network with two time delays, IEEE Trans. Neural Netw. 18 (2) (2007) 416–430. [12] A. Das, A.B. Roy, P. Das, Chaos in a three dimensional general model of neural network, Int. J. Bifurcation Chaos 12 (10) (2002) 2271–2281. [13] J. Fell, N. Axmacher, The role of phase synchronization in memory processes, Nat. Rev. Neurosci. 12 (2011) 105–118. [14] P. Fries, A mechanism for cognitive dynamics: neuronal communication through neuronal coherence, Trends Cogn. Sci. 9 (10) (2005) 474–480. [15] M. Golubitsky, W.F. Langford, Classification and unfoldings of degenerate Hopf bifurcations, J. Differential Equations 41 (3) (1981) 375–415. [16] C.M. Gray, Synchronous oscillations in neuronal systems: Mechanisms and functions, J. Comput. Neurosci. 1 (1994) 11–38. [17] P. Hadley, M.R. Beasley, K. Wiesenfeld, Phase locking of Josephson-junction series arrays, Phys. Rev. B 38 (1988) 8712–8719. [18] B.D. Hassard, N.D. Kazarinoff, Y.H. Wan, Theory and Applications of Hopf Bifurcation, Cambridge University Press, Cambridge, 1981. [19] J. Hopfield, Neural network and physical system with emergent collective computational abilities, Proc. Natl. Acad. Sci. USA 79 (1982) 2554–2558. [20] C. Huang, J. Cao, Hopf bifurcation in an n-dimensional Goodwin model via multiple delays feedback, Nonlinear Dynam. 79 (2015) 2541–2552. [21] E.M. Izhikevich, Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting, The MIT press, London, England, 2007. [22] C. King, The central enigma of consciousness, J. Conscious. Explor. Res. 2 (1) (2011) 013–044. [23] A. Kundu, P. Das, A.B. Roy, Complex dynamics of a four neuron network model having a pair of short-cut connections with multiple delays, Nonlinear Dynam. 72 (3) (2013) 643–662. [24] L. Olien, J. B´elair, Bifurcations, stability and monotonicity properties of a delayed neural network model, Physica D 102 (1997) 349–363. [25] 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. [26] Y. Song, T. Zhang, M.O. Tad´e, Stability switches, Hopf bifurcations, and spatio-temporal patterns in a delayed neural model with bidirectional coupling, J. Nonlinear Sci. 19 (2009) 597–632. [27] F. Tu, X. Liao, W. Zhang, Delay-dependent asymptotic stability of a two-neuron system with different time delays, Chaos Solitons Fractals 28 (2006) 437–447. [28] T. Womelsdorf, J. Schoffelen, R. Oostenveld, W. Singer, R. Desimone, A.K. Engel, P. Fries, Modulation of neuronal interactions through neuronal synchronization, Science 316 (2007) 1609–1612. [29] J. Wu, T. Faria, Y.S. Huang, Synchronization and stable phase-locking in a network of neurons with memory, Math. Comput. Modelling 30 (1999) 117–138. [30] M. Xiao, W.X. Zheng, J. Cao, Bifurcation and control in a neural network with small and large delays, Neural Netw. 44 (2013) 132–142. [31] M. Xiao, W.X. Zheng, J. Cao, Hopf bifurcation of an (n + 1) neuron BAM neural network model with delays, IEEE Trans. Neural Netw. Learn. Syst. 24 (2013) 118–132. [32] W. Xu, J. Cao, M. Xiao, D. Ho, G. Wen, A new framework for analysis on stability and bifurcation in a class of neural networks with discrete and distributed delays, IEEE Trans. Cybern. 99 (2014) 1–13. [33] Y. Yuan, Dynamics in a delayed neural network, Chaos Solitons Fractals 33 (2007) 443–454. [34] Y. Zhang, J. Cao, W. Xu, Stability and Hopf bifurcation of a Goodwin model with four different delays, Neurocomputing 165 (2015) 144–151.