Enhancement of synchronization in inter–intra-connected neuronal networks

JID:PLA

AID:23362 /DIS

Doctopic: Biological physics

[m5G; v1.160; Prn:7/10/2015; 10:34] P.1 (1-7)

Physics Letters A ••• (••••) •••–•••

Contents lists available at ScienceDirect

Physics Letters A www.elsevier.com/locate/pla

Discussion

Enhancement of synchronization in inter–intra-connected neuronal networks F.M. Moukam Kakmeni a,b,∗ , V.M. Nguemaha a,c a

Complex Systems and Theoretical Biology Group, Laboratory of Research on Advanced Materials and Nonlinear Science (LaRAMaNS), Department of Physics, Faculty of Science, University of Buea, P.O. Box 63, Buea, Cameroon b International Centre of Insect Physiology and Ecology, P.O. Box 30772-00100, Nairobi, Kenya c Department of Physics, Florida State University, 77 Chieftan way, Tallahassee, FL 32306, USA

a r t i c l e

i n f o

Article history: Received 5 September 2014 Received in revised form 4 August 2015 Accepted 4 August 2015 Available online xxxx Communicated by Z. Siwy Keywords: Neural synchrony Active control Amplitude death Chaos suppression

a b s t r a c t We study the enhancement of neural synchrony in a network of electrically coupled Hindmarsh Rose (HR) neurons. The behavior of the network under control by an external environment modeled by the Fitzhugh Nagumo (FN) is analyzed. Biologically, such a control system could mimic the modification of normal neuronal dynamics due to drugs or other chemical substances. We show that the environment could have as effect the suppression of chaos, enhancement of synchrony and favor interesting properties such as sub-threshold membrane oscillations, and oscillation death for relatively strong local coupling. Interestingly, we find that the electrical coupling between each two coupled HR and FN is less important to synchronization than the local coupling between the HR and the FN neurons. In other words, local interactions are found to play a stronger role in synchronization than long-range (global) interactions. © 2015 Elsevier B.V. All rights reserved.

1. Introduction A common approach to study synchronization is to consider a network in which the neurons are coupled either chemically or electrically. Synchronization is then analyzed by varying a single parameter that determines the coupling strength between the neurons. The simplistic nature of this architecture limits the level of control over the system and has motivated studies of the combined effect of chemical and electrical synapses on synchronization [1–12]. In this case, two parameters are required to model the network topology, allowing more flexibility in the control of synchronization properties of the system. In the present work, we consider an even more flexible neural network architecture in which each unit is a 3 dimensional (3D) system, which is dynamically coupled to a 2D control. We represent the 3D system with the HR neuron [13,14], whose dynamic is given by Eq. (1).

x˙ = y + ax2 − x3 − z + I e

*

Corresponding author at: Complex Systems and Theoretical Biology Group, Laboratory of Research on Advanced Materials and Nonlinear Science (LaRAMaNS), Department of Physics, Faculty of Science, University of Buea, P.O. Box 63, Buea, Cameroon. E-mail addresses: [email protected], [email protected] (F.M. Moukam Kakmeni). http://dx.doi.org/10.1016/j.physleta.2015.08.004 0375-9601/© 2015 Elsevier B.V. All rights reserved.

y˙ = 1 − bx2 − y z˙ = r [s(x − xo ) − z]

(1)

where x(t ) represents the membrane potential, which is written in dimensionless units. The spiking variables y (t ) and the bursting variable z(t ) describe the rates of ion transport through fast and slow ion channels respectively. The variable a controls the spiking activity of the HR neuron and allows the switching between bursting and spiking. r controls the speed of variation of the slow (bursting) variable and s is a recovery parameter. Following the literatures, values adopted here are a = 1, b = 5, r = 0.006, s = 4. The external input current I e mimics a membrane input excitation current. A minimum initial excitation current of I e ≈ 1.37 is required to set the HR into an active state. In this work, we use I e = 3.2, a value at which the HR is chaotic. The 2D control system is modeled by the spiking Fitzhugh Nagumo (FN) neuron [15], whose dynamics is given by Eq. (2).

u˙ = u − u 3 − v

(2)

v˙ = u + ao − bo v , in which u is the membrane potential and v is the recovery variable representing the force that restores the resting state of the FN system. The parameters are chosen as ao = 0.05 and bo = 0.06, to allow for rapid spiking.

JID:PLA

AID:23362 /DIS

Doctopic: Biological physics

[m5G; v1.160; Prn:7/10/2015; 10:34] P.2 (1-7)

F.M. Moukam Kakmeni, V.M. Nguemaha / Physics Letters A ••• (••••) •••–•••

2

The dynamical properties of the FN and HR neurons are well known and several studies on the synchronization properties of the latter have been reported. However, it is uncertain how a local coupling of these two systems will affect synchronization when a network architecture is considered. Our goal is to study the effect of such local couplings on synchronization. The rest of this work is organized as follows: In Section 2, we analyze the synchronization behavior of a single system unit, consisting of a HR neuron coupled to an FN neuron. An analysis of the bifurcation route of the oscillator is presented. Section 3 deals with the stability of the synchronized system in a network architecture. We determine the contributions of local coupling (interactions within a single network unit) and global coupling (interactions between network units) to the stability of the synchronized network. Our conclusions and remarks appear in Section 4.

2.1. Fixed points and their stability In this section, the analysis of fixed point and their stability is performed. The fixed points are obtained by setting dx = dy = dt dt dz = du = dv = 0 in equation (3). The location of the fixed points dt dt dt (xe , y e , ze , u e , v e ) are determined by solving the coupled system:



2. The coupled HR and FN system The HR and FN models have independently been widely used in computational studies of neural networks. These systems are well known for exhibiting a rich variety of neuronal behaviors, including bursting and spiking. In the absence of electrical synapses, we represent the dynamics of a single unit of coupled HR and FN by

x˙ = y + ax2 − x3 − z + I e + 1 u y˙ = 1 − bx2 − y z˙ = r [s(x − xo ) − z]

Figs. 2(c) and (d) show wave pulses, suggesting a rhythmic interaction between the neurons and enhancement of synchronization. Finally, Figs. 2(e) and (f) ( = 8.7) reveal complete inhibition for the two neurons characterized by persistent sub-membrane oscillations. This suggests that any incoming stimulus has only a hyper polarizing (i.e. making the membrane potential to be more negative) effect on the receiving neuron. This phenomenon is known as oscillation death.

(3)

u˙ = u − u − v + 2 x 3

v˙ = u + ao − bo v . In Eq. (3), the parameters 1 and 2 are used to control the local coupling within a single network unit. For the rest of our analysis and without loss of generality, we will assume 1 = 2 =  . In engineering, a control could be represented by a mechanical system [16,17]. In biological neurons, a control could model the modification of normal neuronal dynamics, due to the action of drugs or other chemical substances introduced in the system. We note that communication between neurons entirely depends on the generation and propagation of the action potential. This is a chemical process that involves a change in ion concentration across the cell membrane. It is therefore not controversial to state that a chemical substance can trigger electrical synapses. The main difference between electrical and chemical synapses resides on the range of interaction. The former are short range while the latter are long range. The time series for the coupled systems for small values of  are plotted in Fig. 1. When the control variable is set to zero Fig. 1(a) and (b), we observe the normal spiking and bursting behaviors for the FN and the HR respectively. That is, the two neurons operate on different duty cycles (fraction of time during which the neuron is active). For relatively small  (weak coupling), we observe that bursting is damped in the HR while the duty cycle of the FN is reduced (Fig. 1(c) to (f)). This suggests that the two neurons exchange their properties towards the achievement of an equilibrium state. Interestingly, we find that the HR neuron is now active even for I e = 0, indicating the emergence of a self-regenerative excitation property. This implies that the system is less sensible to initial conditions and has a higher propensity to achieve synchronization under the action of the control system (FN). On the other hand, Fig. 2 shows time series of the HR and FN for large values of  (strong coupling). A large frequency mismatch (Fig. 2(a) and (b)) is observed, and the time series are patterns of spike-trains interrupted by periods of sub threshold membrane oscillations. This behavior has also been observed in the ION [18].

u e3 + ( b1 − 1)u e + xe3

+ (b

0

− a)xe2

a0 b0

−  xe = 0 + sxe − sx0 − 1 − I −  u e = 0.

(4)

System (4) is a set of cubic polynomials in u e and xe respectively, with discriminants given by:

⎧  3  2 a0 1 ⎪ ⎪ = 4 − 1 + 27 −  x  u e ⎪ e b0 b0 ⎪ ⎪ ⎨ x = 4(b − a)3 (−sx0 − 1 − I −  u e ) e + 27(−sx0 − 1 − I −  u e )2 ⎪ ⎪ ⎪ ⎪ − 54(b − a)(−sx0 − 1 − I −  u e ) ⎪ ⎩ − 9(b − a)2 + 108.

(5)

Hence, each can admit 3 equilibrium points when u e < 0 (resp. xe < 0), 2 equilibrium points when ue = 0 (xe = 0) and a unique equilibrium point when u e > 0 (xe > 0). We note that equilibria of dynamical systems are not always stable. Since stable and unstable equilibria play quite different roles in the dynamics of a system, we find it useful to classify the equilibrium points (xe , y e , ze , u e , v e ) for our system based on their stability. Therefore, we study the stability of (xe , y e , ze , u e , v e ) by analyzing the eigenvalues of the characteristic equation of the jacobian matrix given by (6), obtained by linear approximation of (3).

⎡

⎤

2axe − 3xe2 1 −1  0 ⎢ −2bxe −1 0 0 0 ⎥ ⎢ ⎥ rs 0 −r 0 0 ⎥ J0 = ⎢ ⎢ ⎥ ⎣  0 0 1 − 3u e2 −1 ⎦ 0 0 0 1 −b0

(6)

It is easy to see from properties of matrix algebra that the eigenvalues λ1 = −1, λ2 = −r and λ3 = −b0 have negative real parts, leaving the remaining reduced matrix

 J (xe , u e ) =

−3xe2 − 2(a − b)xe − s





 1 − 3u e2

−

1 b0

(7)

The type of equilibrium point can be determined by the signs of the trace Tr[ J (xe , u e )] = [−3(xe2 + u e2 ) − 2(a − b)xe + (1 − s − b1 )] 0

and the determinant Det[ J (xe , u e )] = [(−3xe2 − 2(a − b)xe − s)(1 − 3u e2 − b1 ) −  2 ] [19]. The eigenvalues λ4/5 of the jacobian J (xe , u e ) 0 are given by

1

λ4/5 = {Tr[ J (xe , u e )] 2  ± {Tr[ J (xe , u e )]}2 − 4Det[ J (xe , u e )]}

(8)

Table 1 below gives a summary of possible types of equilibria according to the values of xe , u e , Tr− , Tr+ , Det− and Det+ . Tr− and Tr + are respectively the negative and the positive zeroes of

JID:PLA

AID:23362 /DIS

Doctopic: Biological physics

[m5G; v1.160; Prn:7/10/2015; 10:34] P.3 (1-7)

F.M. Moukam Kakmeni, V.M. Nguemaha / Physics Letters A ••• (••••) •••–•••

Fig. 1. (Color online.) Time series of membrane potentials for the Hindmarsh Rose (left) and Fitzhugh Nagumo (right) for weak control: (a)–(b)  = 1.5.

3

 = 0.0; (c)–(d)  = 0.5; (e)–(f)

JID:PLA

4

AID:23362 /DIS

Doctopic: Biological physics

[m5G; v1.160; Prn:7/10/2015; 10:34] P.4 (1-7)

F.M. Moukam Kakmeni, V.M. Nguemaha / Physics Letters A ••• (••••) •••–•••

Fig. 2. (Color online.) Time series for the membrane potential of the Hindmarsh Rose (left) and Fitzhugh Nagumo (right) for strong control: (a)–(b) (e)–(f)  = 8.7.

 = 8.0; (c)–(d)  = 8.5;

JID:PLA

AID:23362 /DIS

Doctopic: Biological physics

[m5G; v1.160; Prn:7/10/2015; 10:34] P.5 (1-7)

F.M. Moukam Kakmeni, V.M. Nguemaha / Physics Letters A ••• (••••) •••–•••

5

Table 1 Different types of fixed points of the coupled system. Regions

Values of xe

Values of u e

Sign of Tr

Sign of Det

Type

I II III

xe < Tr+ Tr+ < xe < Tr− xe > Tr−

u e < Det− Det− < u e < Det+ u e > Det+

+ + −

+ − +

focal saddle focal

Fig. 3. Bifurcation diagram of the Hindmarsh Rose (HR) neurons, coupled to the Fitzhugh Nagumo (FN) as a function of the coupling strength ε . The parameters are the same as in Figs. 1 and 2.

Tr[ J (xe , u e )], while Det− and Det+ are zeroes of Det[ J (xe , u e )] for the respective values of Tr − and Tr+ . Based on the above analysis, only unstable focal (I), stable focal (III) and saddle-nodes (II) are likely to occur. It follows that saddle-nodes, symmetry-breaking or period-doubling bifurcations are expected in the model studied here. It is well known that the eigenvalues of the Jacobian matrix at the fixed points plays a crucial role in determining the stability of fixed points. Lyapunov exponents and Lyapunov numbers are related to the eigenvalues and they play an important role respectively in continuous and discrete systems as well. For continuous systems, governed by differential equations like in our case, the Lyapunov exponents are defined as the real parts of the eigenvalues at the equilibrium point. They indicate the rate of contraction or expansion near the equilibrium point when t tends to infinity. The amplitude dead arise in the system when the maximal Lyapunov exponent of the system have a near-zero value within the periodic region of the coupled HR and FN systems (λ2 = −0.006 in our case). This would signify a saddlenode bifurcation of the limit cycle and arrival of the system at a fixed point, indicating oscillation dead [21,26,30,31]. 2.2. The bifurcation diagram In order to visualize the type of bifurcation route to oscillation death, a bifurcation diagram of the system given by Eq. (3) has been plotted in Fig. 3 as a function of the coupling parameter  . For relatively small values (i.e.  < 2.86), the two oscillators are mutually coupled. The chaotic dynamics of the HR model is transmitted to the periodic oscillations of the FN model. For  between 2.86 and 3.48, a window of periodic motion is observed. The periodic attractor abruptly change to double and triple attractor. Further increase of the value of  take the system to a single attractor characterized by oscillation death, after a small window of multi periodicity. Contrary to what is observed in many cases, the route to oscillation death does not follow the common inverse period doubling. In this case, oscillation death occurs abruptly as the coupling between the two system changes. Recently, oscillation death was found to occur under the condition that two oscillators are mutually coupled [20–26]. The oscillation death phenomenon in coupled non-autonomous systems with parametrical excitation has been theoretically identified with two

coupled Duffing oscillators [20]. The existence of diverse routes of transition from amplitude death (AD) to oscillation death (OD) has also been reported [21–25]. However, these studies focused on limit cycles systems first, and then were later extended to coupled complex identical systems. This phenomenon has also been observed to coupled systems with large parameter mismatch [27]. In their work, Ryu et al. [26] analyzed the phenomenon of oscillation death in different oscillators. This implies that even though two different chaotic oscillators are coupled with each other, they exhibit oscillation quenching, just as coupled oscillators with a large parameter mismatch. That is the case here where the two neurons first undergo frequency synchronization, then evolve to amplitude death. The phenomenon observed here is not an isolated case. In fact, it was shown [20] on the basis of the analysis of the equations of Hodgkin and Huxley [28] that the local increase in the density of the protein trans-membrane sodium ion channels leads to a decrease in the excitation threshold of the action potential. More recently, it was proven that the redistribution of ion channels leads to the appearance of regions with high density of channels and areas with ion channel depletion. These areas can block the propagation of the action potential along the axon [29] and result in inactivity of the neuron. In the following, we give a global analysis of the stability of the system and the action of the coupling to a network of coupled HR–FN. 3. Networks analysis of the synchronized system In practice, neurons do not operate as single independent units. They always exist as part of a larger network. In this section, we investigate the synchronization properties of a network of HR–FN units described in Section 2. The units are made to interact through electrical coupling according to the equation

X˙ i = F ( X i ) + K

N 

G i j H ( X j ),

(9)

j =1

where N is the number of units in the network and i = 1, 2, . . . , N. X˙ = F ( X ) is the 5-dimensional differential equation which describes the local coupling (through the parameter  ) of a single HR neuron with the FN as in Eq. (3). K represents the strength of the electrical (global) coupling between nearest neurons. G = (G i j ) is the coupling matrix, and H is a coupling function. Here, G ii = −1 and G i j = 1 if i = j. Because the rows of the coupling matrix G have zero sum, the network allows for a synchronization manifold X 1 = X 2 = . . . = X N = X s , where X s is the synchronous state and X˙s = F ( X s ). If the network is synchronizable, then it will approach the synchronous state asymptotically given a random initial condition. The topology described by Eq. (9) is a ring of electrically coupled neurons, subject to the boundary condition X N +1 = X 1 , X N −1 = X 0 = X N . Each unit interacts directly only with its two closest neighbors. Though this topology is not ideal (in practice, a single neuron in the brain could be linked to about 104 neighbors) it has been used in many interesting studies. For instance, a large number of oscillators given by Eq. (9) can be used to model intestinal signals through synchronization studies [10]. The ring model has also been used to study colorectal myoelectrical activities in humans [11]. In addition, it has been established [4] that

JID:PLA

AID:23362 /DIS

Doctopic: Biological physics

[m5G; v1.160; Prn:7/10/2015; 10:34] P.6 (1-7)

F.M. Moukam Kakmeni, V.M. Nguemaha / Physics Letters A ••• (••••) •••–•••

6

in a network of pulse-coupled neurons, the synchronization behavior of the network depends only on the number of pulses received by each neuron. In this case, the number of neighbors is of little significance. In order to study the stability of the synchronization, we make use of the Master Stability Function (MSF) approach. The MSF measures the exponential rate at which an infinitesimal perturbation of the synchronization manifold grows. For dynamical systems, the MSF is the maximum transverse Lyapunov exponent of the synchronization manifold. A necessary condition for stable synchronization in a network is the negativity of the MSF. The MSF approach is powerful because it allows the local properties the system to be separated from the coupling matrix characterizing the network topology. That is, the MSF can be obtained independently of the network topology, even for networks containing a very large number of neurons [32–36]. The stability of the synchronous state therefore reduces to the study of the system’s dynamical properties along directions that are transverse to the synchronization manifold. Linearizing the system around X s , and following the procedure in [6,32], one finally obtains the following set of N variational equations

η˙ k = [JF( X s ) + K γk JH( X s )]ηk

(10)

where k = 0, 1, . . . , N − 1 and the eigen values of the coupling matrix G for the diffusive nearest neighbor coupling are given by [6, 32–35]

γk = 4K sin2 (

πk N

)

in which k = 1, 2, . . . , N − 1 counts the number of eigen modes for the system. Each equation in Eq. (10) corresponds to a set of 5 j Conditional Lyapunov Exponents λk ( j = 1, . . . , 5) along the eigenmode corresponding to the specific eigenvalue γk . The case k = 0 represents the synchronization manifold (γ0 = 0) and the corresponding maximum conditional Lyapunov Exponent λ10 represents a single isolated HR neuron coupled to a FN neuron. The remaining variations ηk , k = 1, 2, . . . , N − 1 are transverse to the synchronization manifold X s and characterize the system’s response to small perturbations. Complete synchronization then exists if λk1 < 0, for k ≥ 1. The MSF has been plotted in Fig. 4 as function of K and N for 2 values of  . As a reminder, the global behavior of the network is tuned by the parameter K , while the local behavior within a single network unit is controlled by  . As shown in Fig. 4(a), we find that for weak local interactions ( = 0.5), the MSF is largely positive. This suggests that global electrical connections between network units has very little stabilizing effect on the synchronization of the network, even for high values of K . On the other hand, we find interestingly (Fig. 4(b)) that the stability window (region of negative MSF) is largely broadened for  = 1.5. This is true even for weak electrical coupling (small K values). These observations agree with the results of Beligkh et al. [37] which suggest that a network of bursting neurons with strong desynchronizing connections can be synchronized by an input from an external neuron whose duty cycle is sufficiently long. The key point of this analysis is that in a network topology described by Eq. (10), local interactions within units are more important for stable synchronization than global electrical interactions between network units. 4. Conclusion and discussions We have studied the effect a control system modeled by the Fitzhugh Nagumo (FN) on synchronization in a network of electrically coupled Hindmarsh Rose (HR) neurons. The FN and HR systems are two different but complementary artificial models of

Fig. 4. (Color online.) Variation of the Master Stability Function for the controlled system: (a)  = 0.5 The MSF is largely positive indicating unstable synchronization, even for large K; (b)  = 1.5 stable synchronization even for small K values.

real biological neurons. We have demonstrated that in the absence of electrical interaction between the neurons, each HR–FN unit can exhibit a rich variety of dynamical behaviors, including spike-trains, sub threshold membrane oscillations, and even persistent quiescence (oscillation death). We find contrary to what is observed in many cases that the route to oscillation death does not follow the common inverse period doubling. Rather, oscillation death occurs abruptly when the coupling between the two neurons reaches a critical value. This is a special case of synchronization where the systems synchronize to a steady state that may be different from the steady state of each uncoupled system. When the neurons are electrically coupled in a ring-like network, we show that synchronization can be achieved even in the weak coupling regime, due a damping effect from the control system. A Master Stability Function analysis, reveals that in a network setting, local interactions within network units play a stronger role than global electrical interactions in maintaining the stability of the synchronized state. The present results bring more insights to recent experimental findings which suggest the feasibility of closing the loop between seizure prediction and deep brain stimulation for a better control of seizures [38]. The active-control like formalism adopted here could potentially contribute to a deeper understanding of current and future seizure control schemes, as well as other processes that involve synchronization of neurons. Also, the study of hubs that are interconnected is currently of great interest to researchers, mainly because this topology is found in many realistic network models. The method proposed here is appealing because it could allow for the possibility of synchroniza-

JID:PLA

AID:23362 /DIS

Doctopic: Biological physics

[m5G; v1.160; Prn:7/10/2015; 10:34] P.7 (1-7)

F.M. Moukam Kakmeni, V.M. Nguemaha / Physics Letters A ••• (••••) •••–•••

tion without major changes or destructive effects on the network, provided the control is suitably tuned. Acknowledgements F.M. Moukam Kakmeni would like to express his gratitude to the Deutscher Akademischer Austausch Dienst (DAAD), Germany, under the Research Grant section ST 32 Africa, ICN: 91584646, and the International Centre of Insect Physiology and Ecology (ICIPE), for financial and material supports. The authors would like to thank the collaborators and the anonymous reviewer for their valuable comments and suggestions to improve the manuscript. References [1] M.S. Baptista, F.M. Moukam Kakmeni, C. Grebogi, Phys. Rev. E 82 (2010) 036203. [2] P.J. Uhlhaas, W. Singer, Neural synchrony in brain disorders: relevance for cognitive dysfunctions and pathophysiology, Neuron 52 (2006) 155–168. [3] R. Erichsen Jr, L.G. Brunnet, Phys. Rev. E 78 (2008) 061917. [4] Igor Belykh, Enno de Lange, Martin Hasler, Phys. Rev. Lett. 94 (2005) 188101. [5] V. Resmi, G. Ambika, R.E. Amritkar, http://arxiv.org/PS_cache/arxiv/pdf/1011/ 1011.4143v1, 2010. [6] F.M. Moukam Kakmeni, M.S. Baptista, Pramana J. Phys. 70 (2008) 1063. [7] S. Rajesh, Sudeshna Sinha, Somdata Sinha, Phys. Rev. E 75 (2007) 011906. [8] T. Pereira, M.S. Baptista, J. Kurths, Phys. J. 146 (2007) 155. [9] G.X. Qi, H.B. Huang, L. Chen, H.J. Wang, C.K. Shen, Europhys. Lett. 82 (2008) 38003. [10] B. Robertson-Dunn, D.A. Linkens, J. Med. Biol. Eng. 12 (1976) 750. [11] D.A. Linkens, I. Taylor, H.L. Duthie, IEEE Trans. Biomed. Eng. 23 (1876) 101. [12] M.L. Rosa, M.I. Rabinovich, R. Huerta, et al., Phys. Lett. A 266 (2000) 88. [13] J.L. Hindmarsh, R.M. Rose, Nature (London) 296 (1982) 162.

7

[14] J.L. Hindmarsh, R.M. Rose, Proc. R. Soc. Lond. B, Biol. Sci. 221 (1984) 87; R.M. Rose, J.L. Hindmarsh, Proc. R. Soc. Lond. B, Biol. Sci. 161 (1985) 87. [15] R. FitzHugh, Biophys. J. 1 (1961) 445. [16] Rene Yamapi, Stefano Boccaletti, Phys. Lett. A 371 (2007) 48. [17] G.S. Mbouna Ngueuteu, R. Yamapi, P. Woafo, J. Sound Vib. 318 (2008) 1119. [18] Nicolas Schweighofer, Kenji Doya, Mitsuo Kawato, J. Neurophysiol. 82 (1999) 804. [19] J.H. Argyris, G. Fauts, M. Haase, Exploration of Chaos: An Introduction for Natural Scientists and Engineers, Elsevier, Amsterdam, 1994. [20] A.N. Pisarchik, Phys. Lett. A 318 (2003) 65. [21] C.R. Hens, Pinaki Pal, Sourav K. Bhowmick, Prodyot K. Roy, Abhijit Sen, Syamal K. Dana, Phys. Rev. E 89 (2014) 032901. [22] W. Zou, D.V. Senthilkumar, J. Duan, J. Kurths, Phys. Rev. E 90 (2014) 032906. [23] R. Banerjee, E. Padmanaban, S.K. Dana, Pramana 84 (2015) 203. [24] D. Ghosh, T. Banerjee, Phys. Rev. E 90 (2014) 062908. [25] M. Nandan, C.R. Hens, P. Pal, S.K. Dana, Chaos 24 (2014) 043103. [26] J.W. Ryu, D.S. Lee, Y.J. Park, C.M. Kim, J. Korean Phys. Soc. 55 (2009) 395. [27] G.V. Osipov, A.S. Pikovsky, M.G. Rosenblum, J. Kurths, Phys. Rev. E 55 (1997) 2353. [28] P.F.F. Almeida, W.L.C. Vaz, Lateral diffusion in membranes, in: R. Lipowsky, E. Sackmann (Eds.), Handbook of Biological Physics, vol. 1, 1995, p. 305. [29] M.N. Shneider, M. Pekker, Phys. Rev. E 89 (2014) 052713. [30] M.P. Bora, D. Sarmah, Pramana J. Phys. 81 (2013) 677. [31] A. Koseska, E. Volkov, J. Kurths, Phys. Rep. 531 (4) (2013) 173–199. [32] L.M. Pecora, T.L. Carroll, Phys. Rev. Lett. 64 (1990) 821; L.M. Pecora, T.L. Carroll, Phys. Rev. Lett. 80 (1998) 2109. [33] L. Huang, Q. Chen, Y.-C. Lai, L.M. Pecora, Phys. Rev. E 80 (2009) 036204. [34] L.M. Pecora, T.L. Carroll, G.A. Johnson, D.J. Mar, J.F. Heagy, Chaos 7 (1997) 520. [35] S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, D.-U. Hwang, Phys. Rep. 424 (2006). [36] S. Boccaletti, J. Kurths, G. Osipov, D.L. Valladares, C.S. Zhou, Phys. Rep. 366 (2002). [37] I. Belykh, A. Shilnikov, Phys. Rev. Lett. 101 (2008) 078102. [38] L.B. Good, S. Sabesan, S.T. Marsh, K. Tsakalis, D. Treiman, L. Iasemidis, Int. J. Neural Syst. 19 (2009) 173.