Chimera patterns in three-dimensional locally coupled systems
CChimera patterns in three-dimensional locally coupled systems
Srilena Kundu , Bidesh K. Bera , , Dibakar Ghosh , ∗ and M. Lakshmanan Physics and Applied Mathematics Unit, Indian Statistical Institute, 203 B. T. Road, Kolkata-700108, India Department of Mathematics, Indian Institute of Technology Ropar, Punjab-140001, India Centre for Nonlinear Dynamics, School of Physics,Bharathidasan University, Tiruchirapalli-620024, India (Dated: February 5, 2019)The coexistence of coherent and incoherent domains, namely the appearance of chimera states, isbeing studied extensively in many contexts of science and technology since the past decade, thoughthe previous studies are mostly built on the framework of one-dimensional and two-dimensionalinteraction topologies. Recently, the emergence of such fascinating phenomena has been studiedin three-dimensional (3D) grid formation while considering only the nonlocal interaction. Here westudy the emergence and existence of chimera patterns in a three dimensional network of coupledStuart-Landau limit cycle oscillators and Hindmarsh-Rose neuronal oscillators with local (nearestneighbour) interaction topology. The emergence of different types of spatiotemporal chimera pat-terns is investigated by taking two distinct nonlinear interaction functions. We provide appropriateanalytical explanations in the 3D grid of network formation and the corresponding numerical jus-tifications are given. We extend our analysis on the basis of Ott-Antonsen reduction approachin the case of Stuart-Landau oscillators containing infinite number of oscillators. Particularly, inHindmarsh-Rose neuronal network the existence of non-stationary chimera states are characterizedby instantaneous strength of incoherence and instantaneous local order parameter. Besides, thecondition for achieving exact neuronal synchrony is obtained analytically through a linear stabilityanalysis. The different types of collective dynamics together with chimera states are mapped overa wide range of various parameter spaces.
PACS numbers: 87.10.-e, 05.45.Xt
I. INTRODUCTION
The study on the collective behaviour of coupled dy-namical systems has drawn a significant attention duringthe past two decades due to its immense applicabilityin several physical, chemical, biological and engineeringsystems. Among the varied collective dynamics, chimerastate [1–3] is one of the most studied research topics inrecent times and since its identification, the collectivedynamics of identical oscillators have been revitalized.Chimera state [4, 5] is a complex spatiotemporal patternthat consists of both coherent and incoherent domains si-multaneously and appear in a network of symmetricallycoupled identical oscillators. The existence and robust-ness of chimera states has been verified through a seriesof experiments in different systems, including chemicaloscillators [6], electrochemical oscillators [7], electroniccircuits [8], opto-elctronic systems [9], mechanical sys-tems [10] and frequency modulated time delay systems[11].Some aquatic mammals such as eared seal, dolphin,some migratory birds and manatees are engaged in theslow wave unihemispheric sleep [12, 13]. During theirsleeping state, half of their brain is shut down while theremaining portion of the hemisphere is awake to monitorwhat is happening in the environment. The neuronal os-cillations in the wake portion of the cerebral hemisphere ∗ Electronic address: [email protected] is desynchronized whereas oscillations in the sleepy partis very much synchronized. This type of neuronal processin the brain is closely related to Kuramoto’s observationof the coexistence of coherent and incoherent motionsin a network of nonlocally coupled identical phase os-cillators [4]. Later, Abrams et al. [5] provided sometheoretical explanations for the existence of such statesand named it as chimera state . The chimera-like fea-tures are also observed in many man-made systems andnatural phenomena such as power grid networks [14, 15],superconducting metamaterial [16], ventricular fibrilla-tion [17, 18] and different types of brain diseases [19, 20],like epileptic seizures, brain tumors and schizophrenia.The chimera state is a fascinating type of symmetry-breaking complex pattern and appears in symmetry pre-serving systems. After the detection of chimera statesin a network of phase oscillators [4, 5], the existence ofthese states has also been proved in networks of chaotic[21] and hyper chaotic oscillators [22], limit cycle oscilla-tors [23, 24], chaotic maps [25], time delayed oscillators[26] and even in neuronal oscillators [27] which exhibitdifferent types of bursting dynamics. At the beginningit was believed that nonlocal network configuration isthe essential requirement for the existence and emergenceof chimera states since such coupling topology producesmultistability in the systems and leads to the appear-ance of such a state. But recent results revealed thatproper design of the coupling functions and the nodaldynamics may generate such states in globally (all-to-all) [28–34] coupled networks and locally (nearest neigh-bour) coupled rings [35–38] and arrays [24] of oscillators. a r X i v : . [ n li n . C D ] F e b Besides these symmetric interacting configurations, theexistence of chimera states has also been uncovered inrandom and scale free networks [39], heterogeneous [40],multiplex networks [41–45], time varying networks [46],modular networks [47] and multiscale networks [48]. Theexistence of chimera states is also detected in two inter-acting different populations where each oscillator in eachpopulation interacts in nonlocal [49] and global fashion[50]. Recently, the emergence of alternating chimera in anetwork of identical neuronal systems induced by an ex-ternal electromagnetic field has been reported in [51]. De-pending on the spatiotemporal motion of the oscillatorsin the network, different types of non-stationary chimerapatterns are categorized which include imperfect chimera[52], travelling chimera [53], imperfect travelling chimera[54], and spiral wave chimera states [55, 56]. From theamplitude variation of each individual oscillators of thenetwork, various types of chimera states are classifiedsuch as chimera death [57], amplitude mediated chimera[58] and so on.In real world systems, the human brain is one of themost challenging complex systems due to the presenceof large number of neurons and myriads of interconnec-tions between them. Also the understanding of the neu-ronal communication and interaction patterns betweenthe neurons through the synapses is also a real chal-lenging issue. Mainly two structurally different types ofsynapses are identified in the interneuronal communica-tions, one is electrical gap junction and another is chem-ical synapse. In the electrical communication, the neu-ronal information are exchanged between two adjacentneurons by making the gap junction channel in the mem-brane potential, while through the chemical synapses thesignal is transferred chemically in the form of neurotrans-mitter molecules among the neurons. Normally, neuronsin the brain are situated and pass information in a twodimensional (2D) grid formation for their normal opera-tion. The existence and emergence of the chimera statesin the 2D grid of oscillations has been investigated inthe neuronal network by taking electrical and chemicalsynapses in nonlocal [59] and local [60] modes. Also,recently, the presence of chimera patterns have been re-ported in 3D network of neuron oscillators having nonlo-cal interaction topology [61].In this work, we systematically investigate the chimerapatterns in networks of locally coupled oscillators in 3Dgrid formation, which can be considered as a generaliza-tion of the 2D grid network. Here we consider the oscilla-tors in 3D grid to be interacting via a nonlinear couplingfunction along with local coupling topology. Previouslychimera state was observed to appear in 1D [38] and 2D[60] grid of locally coupled oscillators with nonlinear in-teraction function. Here also we observe that the pres-ence of nonlinearity in the coupling function plays a cru-cial role for the existence of chimera states in 3D grid oflocally coupled oscillators. To understand the mechanismof such emerging phenomena, we start with a network oflocally coupled Stuart-Landau oscillators in the 3D grid formation where the interaction takes place through anonlinear coupling function. First we numerically showthe existence and emergence of chimera states, whose sta-tionarity is justified by the long term spatiotemporal plot,and then by using Ott-Antonsen (OA) reduction method,we perform relevant calculations on the continuum limitapproach. The results obtained from this approach arefound to be in good agreement with the numerical find-ings. Next we extend our study on the realistic neu-ronal network where the Hindmarsh-Rose neurons areconsidered as local dynamical units of each node of thenetwork. Using nonlinear chemical synaptic interactionfunction in the local coupling topology, the emergenceof chimera states are articulated. In this case, our ob-served chimera pattern is of a non-stationary type whichmeans that the coherent and incoherent populations ofthe chimera states are changing erratically in space andtime. The existence of such spatiotemporal dynamics ischaracterized through instantaneous strength of incoher-ence and local order parameter. Moreover, the analyticalbound for exact synchronization in the case of synapti-cally coupled neuronal network in 3D local interactiontopology is derived using a linear stability analysis, andfurther verified by corresponding numerical results usingthe notion of Kuramoto order parameter. We further ex-plore the transitions between the several dynamical statesand map them over a wide range of different parameterspaces.The remaining parts of the paper are organized as fol-lows. In Sec. II, we present a general mathematicaldescription of the three dimensional grid of oscillators.The numerical and analytical investigations on the exis-tence and emergence of chimera patterns in a 3D grid ofStuart-Landau oscillators is provided in Sec. III. Thenthe subsequent Sec. IV focusses on the occurrence ofchimera patterns in the 3D grid of neuronal network andin the sub Sec. IV A, we derived the analytical conditionfor the neuronal synchrony threshold. The conclusions ofour findings are presented in Sec. V.
FIG. 1: Schematic diagram of the coupling scheme in 3D cubiclattice: the ( i, j, k )-th oscillator (red circle) is connected to itssix nearest-neighbour oscillators (green circles).
FIG. 2: Snapshots of the state variables x (real part of z ) in the (a) 3D grid of N oscillators, (b) 2D plane with N oscillatorskeeping k = 40 fixed, (c) in 1D with N oscillators keeping j = 20 , k = 40 fixed, showing the emergence of chimera states in thewhole network of Stuart-Landau oscillators (2). Here α = 1 . , β = − . , ˜ a = 1 . , (cid:15) = 0 .
09, and N = 40 are considered. II. GENERAL MATHEMATICAL FORM OF 3DCOUPLED OSCILLATORS
We here consider a network of N × N × N cubic latticeof locally coupled oscillators. The systematic interactionscheme is shown in Fig. 1. The local dynamics of an indi- vidual node of the network is given by ˙ X i,j,k = F ( X i,j,k ) , where X i,j,k represents an l -dimensional vector of the dy-namical state variables and F ( X i,j,k ) is the correspond-ing velocity field. The general mathematical equations oflocally coupled systems in a 3D cubic lattice of networkcan be described as˙ X i,j,k = F ( X i,j,k ) + K { H ( X i,j,k , X i − ,j,k ) + H ( X i,j,k , X i +1 ,j,k ) + H ( X i,j,k , X i,j − ,k ) + H ( X i,j,k , X i,j +1 ,k )+ H ( X i,j,k , X i,j,k − ) + H ( X i,j,k , X i,j,k +1 ) } , (1)where the subscripts ( i, j, k )( i, j, k = 1 , , ..., N ) in X i,j,k and F ( X i,j,k ) determine the position of the oscillatorin the 3D coupled network. The coupling function H : R l × R l → R describes the manner by which the( i, j, k )-th oscillator is connected with its nearest neigh-bouring oscillators. Here K = ( (cid:15) , (cid:15) , ..., (cid:15) l ) T is the cou-pling matrix, where T denotes transpose of a matrix.We use periodic boundary conditions as X ,j,k = X N,j,k , X N +1 ,j,k = X ,j,k and X i, ,k = X i,N,k , X i,N +1 ,k = X j, ,k and X i,j, = X i,j,N , X i,j,N +1 = X i,j, .In the subsequent sections, we provide adequate ev-idence for the existence and emergence of the chimerastates in 3D lattice of locally coupled oscillators by con-sidering two different dynamical systems, namely Stuart-Landau oscillators and Hindmarsh-Rose neuronal model.For Stuart-Landau oscillators, we consider “pull-push”type nonlinear coupling function and chemical synapticinteraction function is used in the Hindmarsh-Rose neu-ronal network. III. COUPLED STUART-LANDAU SYSTEM
We start our investigation by considering a network ofcoupled identical Stuart-Landau (SL) oscillators in thethree dimensional grid formation where each node of the network is interacting with the local neighbouring nodes.The dynamical evolution of the prescribed network is de-scribed by following set of coupled equations,˙ z i,j,k = (1 + ˆ iα ) z i,j,k − (1 + ˆ iβ ) | z i,j,k | z i,j,k + (cid:15) [ H ( z i − ,j,k ) + H ( z i +1 ,j,k ) + H ( z i,j − ,k ) + H ( z i,j +1 ,k + H ( z i,j,k − ) + H ( z i,j,k +1 ) − H ( z i,j,k )] , (2)for i, j, k = 1 , ..., N and N is the total number of oscil-lators in each direction of the three grid network. Thevariable z = x + ˆ iy is a complex variable with ˆ i = √− α , β are real parameters. The parameter (cid:15) denotesthe coupling strength and H ( z ) = ˜ a z − z | z | , representsthe nonlinear interaction function [62] with real constant˜ a . In contrast, if H is considered to be a linear cou-pling function i.e., simple diffusion (say), then the abovesystem (2) may not be able to produce chimera states,as observed in locally coupled 2D grid of oscillators [60].In absence of coupling strength, individual SL systemsexhibit limit cycle behaviour near the Hopf bifurcationpoint with intrinsic frequency α .In order to analyze system (2), we first numericallyshow the existence of chimera states in the network ofSL oscillators. To integrate the coupled system (2), weuse Euler integration scheme with time step ∆ t = 0 . x i,j,k (0) = 0 . N − ( i + j + k )) , y i,j,k (0) = 0 . N − ( i + j + k )) for i, j, k = 1 , , ..., N along with some additionalsmall random fluctuations. Here we considered N = 40in linear dimension of the 3D cubic lattice which accountsfor N = 64000 as the total number of oscillators in theentire 3D network. Figure 2(a) shows the snapshot ofthe state variables x i,j,k over the entire 3D lattice. Co-existence of coherence and incoherence and consequentlythe chimera pattern is easily identifiable from the figure.Snapshot of the state variables in 2D projection on theplane k = 40 is depicted in Fig. 2(b), whereas the 1Dsnapshot (with j = 20 , k = 40) of the oscillators distin-guishing the domain of coherent and incoherent popula-tions is shown in Fig. 2(c). The parameter values arefixed as α = 1 . , β = − . , ˜ a = 1 . , (cid:15) = 0 . x i,j,k = Re ( z i,j,k ) for the oscillators situ-ated along j = 20 and k = 40, while keeping the inter-action strength as (cid:15) = 0 .
09 fixed in the chimera region.The appearance of stationarity in the observed chimerapattern of 3D network of locally coupled SL oscillators isclearly visible from the figure. Moreover, to characterizethe ordination of coherence and incoherence among theoscillators, we use the measure of local order parameter [54] L i,j,k , which quantifies the degree of local orderingof the oscillators. The quantity is calculated using theusual formula L i,j (cid:48) ,k (cid:48) = (cid:12)(cid:12)(cid:12) δ (cid:80) | i − i (cid:48) |≤ δ e ˆ i Φ i (cid:48) (cid:12)(cid:12)(cid:12) , j (cid:48) = 20 , k (cid:48) = 40 . (3)Here Φ i = arctan ( y i,j (cid:48) ,k (cid:48) /x i,j (cid:48) ,k (cid:48) ) is the geometric phaseof the ( i, j (cid:48) , k (cid:48) )-th oscillator, δ is the number of oscil-lators in the neighborhood of ( i, j (cid:48) , k (cid:48) )-th oscillator andˆ i = √−
1. In Fig. 3(b), the spatiotemporal evolution ofthe local order parameter L i,j (cid:48) ,k (cid:48) has been depicted. Thevalue of the order parameter L i,j (cid:48) ,k (cid:48) (cid:39) i, j (cid:48) , k (cid:48) )-th oscillator belongs to the coherent part of thechimera state, whereas L i,j (cid:48) ,k (cid:48) < A. Analytical results: Ott-Antonsen Approach
To support the results obtained from our numericalexperiments and to validate them in the continuum limitwe go through the approach of Ott and Antonsen [63, 64].Following the procedure discussed in [60], we can directlyobtain the phase reduced model for the 3D coupled SLoscillators as˙ θ i,j,k = ω (cid:48) i,j,k − λ N (cid:88) l =1 N (cid:88) m =1 N (cid:88) n =1 A ijklmn sin( θ i,j,k − θ l,m,n + γ ) , (4)where ω (cid:48) i,j,k = ω i,j,k + (cid:15)β (˜ a − λ = (cid:15) (˜ a − (cid:112) β and γ = tan − β . Since all the systems are identical here weconsider ω i,j,k = ω, ∀ i, j, k = 1 , , ..., N with ω = α − β as in [60]. It should be noted that the above equation can berewritten in the continuous form as ∂θ ( x,y,z,t ) ∂t = ω (cid:48) − λ (cid:82) (cid:82) (cid:82) G ( x − x (cid:48) , y − y (cid:48) , z − z (cid:48) ) sin( θ ( x, y, z, t ) − θ ( x (cid:48) , y (cid:48) , z (cid:48) , t ) + γ ) dx (cid:48) dy (cid:48) dz (cid:48) . (5)Here the coupling kernel G can be written as G ( x − x (cid:48) , y − y (cid:48) , z − z (cid:48) ) = H [cos( (cid:112) ( x − x (cid:48) ) + ( y − y (cid:48) ) + ( z − z (cid:48) ) )2 π − cos(2 π/N )] and H ( . ) is the Heaviside step function.Considering the limiting case of an infinite number of oscillators, i.e. N → ∞ , the dynamics of the oscillators attime t can be described by a probability density function f ( x, y, z, θ, t ) which satisfies the continuity equation, ∂f∂t + ∂∂θ ( f v ) = 0 , (6)where v = dθdt = ω (cid:48) − i [ re ˆ iθ + ¯ re − ˆ iθ ] . (7)It is to be noted that here r is the order parameter given by, r ( x, y, z, t ) = λe ˆ iγ (cid:82) (cid:82) (cid:82) G ( x − x (cid:48) , y − y (cid:48) , z − z (cid:48) ) (cid:82) π e − ˆ iθ f ( x (cid:48) , y (cid:48) , z (cid:48) , θ, t ) dθdx (cid:48) dy (cid:48) dz (cid:48) . (8)Then the probability density function f ( x, y, z, θ, t ) can be expanded in terms of the Fourier series taking intoaccount the OA ansatz f n ( x, y, z, θ, t ) = h ( x, y, z, t ) n as f ( x, y, z, θ, t ) = π (cid:16) (cid:80) ∞ n =1 h ( x, y, z, t ) n e ˆ inθ + c.c. (cid:17) = π (cid:16) (cid:80) ∞ n =1 ( h n e ˆ inθ + ¯ h n e − ˆ inθ ) (cid:17) . (9) FIG. 3: Spatiotemporal evaluation of (a) the state variable x i,j,k , and (b) corresponding local order parameter L i,j,k of theSL oscillators. The plots are obtained along the cross-section j = 20 , k = 40, keeping the coupling strength (cid:15) = 0 .
09 fixed inthe chimera region.FIG. 4: (a) Snapshots of the local order parameter | h ( x, y, z ) | and (b) the local phase Ψ( x, y, z ) at a particular instant of time.(c) Spatiotemporal evolution of the phase Ψ( x, y, z, t ) for a long time interval, indicating stationarity of the observed chimerapattern. Substituting (9) into (6), we obtain the governing equation for the considered 3D network on the basis of OAreduction ∂h∂t = − ˆ iω (cid:48) h + (cid:0) ¯ rh + r (cid:1) , (10)where r ( x, y, z, t ) = λe ˆ iγ (cid:82) (cid:82) (cid:82) G ( x − x (cid:48) , y − y (cid:48) , z − z (cid:48) ) h ( x (cid:48) , y (cid:48) , z (cid:48) , t ) dx (cid:48) dy (cid:48) dz (cid:48) . (11)The modulus and argument of the complex OA ansatz h ( x, y, z, t ) = | h ( x, y, z, t ) | e − ˆ i Ψ( x,y,z,t ) respectively pro-vide information about the local order parameter and thelocal phase of the oscillators in the continuum limit. Theexistence of the appropriate complex OA ansatz func-tion h ( x, y, z, t ) is the condition for the stability of theobtained chimera state.Using Eqs. (10) and (11), we have enlightened theappearance of chimera states through the complex OAreduction approach in a 3D network consisting of infi-nite number of locally coupled SL oscillators. Figure 4presents the results obtained through numerical integra-tion of Eqs. (10) and (11). In Figs. 4(a) and 4(b) respec-tively, we have plotted the snapshots of the local orderparameter | h ( x, y, z ) | and the local phase Ψ( x, y, z ) of theoscillators located in the considered 3D grid at a particu-lar instant of time. Here in Fig. 4(a), | h ( x, y, z ) | = 1characterizes the domain of the phase locked oscilla-tors (i.e., the coherent group), whereas, | h ( x, y, z ) | < x, y, z ) to verify that the observed chimera pat-terns are stationary in time even in the continuum limitcase. From these results, we can conclude that the OAreduction technique on the basis of infinite number of os-cillators certainly predicts the same dynamics, which waspreviously realized through numerical analysis. IV. COUPLED HINDMARSH-ROSE NEURONMODEL
Next we consider the dynamics of the network of a 3Dcubic lattice of coupled Hindmarsh-Rose neuron model,consisting of N = 27000 neuronal oscillators. The localdynamics is modeled through Hindmarsh-Rose (HR) [65]neuron, F ( X ) = [ ax − x − y − z, ( a + α ) x − y, c ( bx − z + e )] T with coupling matrix K = diag ( (cid:15) , ,
0) withcoupling function H ( x (cid:48) , x ) = ( v s − x (cid:48) )Γ( x ), where Γ( x ) isthe chemical synaptic function. After local coupling, thecomplete set of equations becomes,˙ x i,j,k = ax i,j,k − x i,j,k − y i,j,k − z i,j,k + (cid:15) ( v s − x i,j,k ) { Γ( x i − ,j,k ) + Γ( x i +1 ,j,k )+ Γ( x i,j − ,k ) + Γ( x i,j +1 ,k ) + Γ( x i,j,k − ) + Γ( x i,j,k +1 ) } , ˙ y i,j,k = ( a + α ) x i,j,k − y i,j,k , ˙ z i,j,k = c ( bx i,j,k − z i,j,k + e ) , i, j, k = 1 , , ..., N. (12)The variable x j,j,k denotes the membrane potential of the( i, j, k )-th position of neuron in the 3D grid of neuronalnetwork. The other two state variables y i,j,k and z i,j,k areassociated with the ion transportation across the mem-brane potential in which y i,j,k represents the fast currentassociated with N a + or K + ions and z i,j,k is the slowcurrent associated with Ca + ions respectively. The pa-rameter c is controlled by the ratio of slow-fast dynamics.For the excitatory interaction, we choose v s = 2 so thatthe condition x i,j,k ( t ) < v s is always maintained for alltime t . The interaction function between the neuronsis taken as a chemical synaptic function Γ( x ) which isa nonlinear sigmoidal input-output function, describedby Γ( x ) = e − λ ( x − Θ s ) . The parameters λ and Θ s rep-resent the slope of the sigmoidal function and synapticfiring threshold, respectively. The constant (cid:15) denotesthe chemical synaptic interaction strength which sum-marizes how the information is delivered through the in-teraction among the neurons. In particular, if electricalsynapses are employed to interact between the neuronsin 3D grid then the desired chimera states cannot beproduced as verified in [60] for 2D grid. We fix the lo-cal system parameter values as a = 2 . , α = 1 . , b =9 . , c = 0 . , e = 5 . v s = 2 . , λ = 10 . , Θ s = − .
25. Next we study the be-havior of all the neurons in the network by varying thecoupling strength (cid:15) with N = 30.Now we study the emergence of chimera states in alocally connected 3D grid consisting of N = 30 HR os-cillators in each direction numerically. So the total num-ber of neurons corresponding to the linear dimensionin the 3D grid formation is N = 27000. For the nu-merical simulation, we again used the Euler integrationscheme with time-step size ∆ t = 0 .
01. The initial con-ditions for each node of the large neuronal network ischosen as x i,j,k (0) = 0 . N − ( i + j + k )) , y i,j,k (0) =0 . N − ( i + j + k )) , z i,j,k (0) = 0 . N − ( i + j + k ))for i, j, k = 1 , , ..., N with small added random fluctua-tions. Different dynamical states are portrayed in Fig. 5for various coupling strengths and the color bar denotesthe values of the membrane potential x i,j,k attained byeach neuron at a particular time t = 1495. The left col-umn of Fig. 5 represents the incoherent states on 3DHR network for the coupling strength (cid:15) = 0 .
12. Figs.5(a, d, g) represent the snapshots of the membrane po-tential x i,j,k in 3D space, 2D plane (taking cross-sectionalong k = 30) and 1D array (keeping j = 18 , k = 30fixed) respectively. In these figures, all the neurons are FIG. 5: Left, middle and right columns respectively represent (a, d, g) incoherent, (b, e, h) chimera and (c, f, i) coherentstates emerging in a 3D cubic lattice of HR neuron. Top, middle and bottom row respectively represent the (a, b, c) 3Dsnapshots plotted in ( i, j, k ) space, (d, e, f) 2D snapshots plotted in ( i, j ) plane keeping k = 30 fixed, and (g, h, i) 1D snapshotsplotted with respect to i keeping j = 18 , k = 30 fixed. The coupling strengths for incoherent, chimera and coherent states arerespectively (cid:15) = 0 . , . .
44. All the snapshots are taken at time t = 1495. randomly distributed which is a clear significance of in-coherent states. Similarly, the snapshots in the middlecolumn of Fig. 5 show the appearance of chimera statesin different dimensions for the synaptic coupling strength (cid:15) = 1 .
2. Here the group of synchronized neurons are sur-rounded by the population of desynchronized neurons. InFigs. 5(b) and 5(e), the blue regions correspond to thecoherent population of neurons and surrounding regionswith mixed colors indicate the incoherent group of neu-rons. This feature is also evident in the 1D snapshot (Fig.5(h)), where a group of neurons (14 ≤ i ≤
20) followsa smooth profile while the rest are randomly scattered.This dynamical behavior is a strong indication of the ex-istence of chimera patterns. Now for sufficiently highervalues of the synaptic interaction strength, (cid:15) = 1 .
44, allthe neurons in the coupled 3D network follow a coherentmotion and the corresponding snapshots are shown in theright column of Figs. 5 considering 3D, 2D and 1D cases.The continuous color variations (Figs. 5(c, f)) and thesmooth profile (Fig. 5(i)) of each neuron are bearing astrong resemblance of coherent states in the 3D neuronalnetwork.Next we move on to investigate the stationarity of the chimera patterns and we find that our observed chimerastate is not static, rather it varies in time. To prove theemergence of non-stationary chimera patterns, we plotthe spatiotemporal dynamics of the neuronal network inFig. 6(a) and its existence is characterized through theinstantaneous order parameter plot (using Eq. (3) with j (cid:48) = 18 , k (cid:48) = 30) in Fig. 6(b). Taking the cross-sectionalong j = 18 and k = 30, the space-time plot is drawnat coupling strength (cid:15) = 0 .
9. In the chimera pattern, theblue region is the coherent state and mixed colors cor-respond to the incoherent states in Fig 6(a). These twopopulations oscillate with respect to time and the instan-taneous local order parameter signifies the appearance ofthis feature since it takes the value 1 from color bar (darkred region of Fig. 6(b)) for the coherent states and takesother lower values for the incoherent states.To characterize the appearance of chimera states in 3Dlocally coupled network, we have considered the quantita-tive statistical measure based on the concept of strengthof incoherence (SI) as developed by Gopal et al. [22]. Forthis we begin by defining the transformations of the statevariables x i,j,k to new variables, w i,j,k = (cid:112) ( x i,j,k − x i +1 ,j,k ) + ( x i,j,k − x i,j +1 ,k ) + ( x i,j,k − x i,j,k +1 ) , i, j, k = 1 , , ..., N. (13) FIG. 6: Spatiotemporal evolution of (a) membrane potential x i,j,k , and (b) corresponding local order parameter L i,j,k of theoscillators. The plots are obtained along the cross-section j = 18 , k = 30, keeping the coupling strength (cid:15) = 0 . x of N neurons situated on the 3D lattice along the cross-section j = 18 , k = 30, (b)variation of SI with respect to time of all the oscillators in the network for a particular coupling strength (cid:15) = 0 .
9. (c) Scenarioof the strength of incoherence (SI) measure with respect to time and coupling strength.
Next we divide the whole three dimensional medium into M (= M h × M w × M b ) number of cubic bins each of equalvolume n (= N/M h × N/M w × N/M b ) and then calculate the local standard deviation σ ( p, q, r ) in each of these binsas σ ( p, q, r )( t ) = (cid:115) n n p (cid:80) i = n ( p − n q (cid:80) j = n ( q − n r (cid:80) k = n ( r − ( w i,j,k − (cid:104) w (cid:105) ) , (14)where (cid:104) w (cid:105) = N N (cid:80) i =1 N (cid:80) j =1 N (cid:80) k =1 w i,j,k , p = 1 , , ..., M h ; q = 1 , , ..., M w ; r = 1 , , ..., M b . Then the strength of incoherenceis defined as SI ( t ) = 1 − Mh (cid:80) p =1 Mw (cid:80) q =1 Mb (cid:80) r =1 s ( p,q,r ) M , s ( p, q, r ) = Θ( δ − σ ( p, q, r )) , (15)where Θ( . ) is the Heaviside step function and δ is thepredefined threshold. Consequently, the values SI ( t ) =1 , SI ( t ) = 0 and 0 < SI ( t ) < x i,j,k corresponding to Fig. 6(a)is plotted in Fig. 7(a). From this figure it can be ob-served that, all the neurons are in a quiescent states whilefollowing the coherent motion (say, 1660 < t < < t < SI ( t ) measure-ment, the coherence-incoherence patterns are character-ized. The variation of SI ( t ) is shown in Fig. 7(b) cor-responding to Fig. 6(a) and Fig. 7(a) for particular (cid:15) = 0 .
9. Here SI ( t ) takes the value 1 for desynchronizedstates (spiking region) and takes the value 0 for coherentstates (quiescent states) while intermediate values referto the chimera states. To reveal the coupling effect inthe emergence of chimera patterns, we draw a phase dia-gram using SI ( t ) measurement in Fig. 7(c) for varying (cid:15) by considering a long time interval. The color bar of Fig.7(c) shows the variation of SI ( t ) in which SI ( t ) = 1 and0 refer to the incoherent and coherent states respectively.The dark red and blue regions represent the incoherentand coherent states respectively while the intermediatecolors correspond to the chimera states. For very lowcoupling values, all the neurons are in incoherent statefor all times and chimera states (i.e., the co-existence ofcoherence and incoherence) appear after certain thresh-old value of (cid:15) > .
13. Further increase in the couplingstrength influences the dominating nature of the coherentstates, leading to a shrinkage of incoherent and chimeraregions.The effect of the network size on the emergence ofnon-stationary chimera patterns in the 3D neuronalnetwork is also discussed in Fig. 8. For the couplingstrengths (cid:15) = 0 . .
08, the dependency of N on chimera pattern with respect to time is plotted in Fig.8(a) and Fig. 8(b), respectively, and characterizedthrough instantaneous strength of incoherence. Fromthis figure it is clear that the imperfectness of non-stationary chimera patterns is reduced with an increasein the synaptic coupling strength and also the chimerapatterns are independent of size after certain networksize N = 9 . A. Synchronization threshold for 3D networkmodel
Following the approach of Belykh et al. [66] to calcu-late the synchronization threshold in the case of chem-ical synaptic neuronal networks, here we can also cal-culate the threshold value in the case of our considerednetwork with 3D lattice architecture. The important cri-terion to apply the mechanism proposed by Belykh etal. [66] for finding synchronization threshold in any ar-bitrary network topology is that the number of signalsreceived by each neuron should be equal to n i.e., the de-gree of each of the neurons should be equal. Motivatedby this, here we want to calculate the threshold for com-plete synchrony in case of our considered network (12)with 3D coupling topology, where the number of signalsreceived by each neuron is equal to six, i.e. n = 6. Onthe onset of complete synchronization each of the oscil-lators follow the single state x i,j,k ( t ) = x ( t ) , y i,j,k, ( t ) = y ( t ) , z i,j,k ( t ) = z ( t ) for i, j, k = 1 , , ..., N . We can saythat the whole 3D network of N oscillators synchronizes,when each of the N oscillators in each of the i th, j th and k th direction synchronizes. So for simplicity, we can cal-culate the threshold value by considering the oscillatorsalong 1D, say along i direction (for any arbitrary valuesof j = j (cid:48) , k = k (cid:48) ). Then if the N number of oscillatorsalong ( i, j (cid:48) , k (cid:48) ) synchronizes, the whole network of N os-cillators synchronizes. We can rewrite the equation of thenetwork model (12) along this specific ( i, j (cid:48) , k (cid:48) ) directionas0 FIG. 8: Variation of SI with respect to time obtained by varying the number of nodes in the network for two different couplingstrengths (a) (cid:15) = 0 . (cid:15) = 1 . ˙ x i,j (cid:48) ,k (cid:48) = ax i,j (cid:48) ,k (cid:48) − x i,j (cid:48) ,k (cid:48) − y i,j (cid:48) ,k (cid:48) − z i,j (cid:48) ,k (cid:48) + (cid:15) (cid:48) ( v s − x i,j (cid:48) ,k (cid:48) ) N (cid:80) i (cid:48)(cid:48) ,j (cid:48)(cid:48) ,k (cid:48)(cid:48) =1 C ii (cid:48)(cid:48) j (cid:48) j (cid:48)(cid:48) k (cid:48) k (cid:48)(cid:48) Γ( x i (cid:48)(cid:48) ,j (cid:48)(cid:48) ,k (cid:48)(cid:48) ) , ˙ y i,j (cid:48) ,k (cid:48) = ( a + α ) x i,j (cid:48) ,k (cid:48) − y i,j (cid:48) ,k (cid:48) , ˙ z i,j (cid:48) ,k (cid:48) = c ( bx i,j (cid:48) ,k (cid:48) − z i,j (cid:48) ,k (cid:48) + e ) , i = 1 , , ..., N. (16)Here C ii (cid:48)(cid:48) j (cid:48) j (cid:48)(cid:48) k (cid:48) k (cid:48)(cid:48) is the connectivity matrix where N (cid:80) i (cid:48)(cid:48) ,j (cid:48)(cid:48) ,k (cid:48)(cid:48) =1 C ii (cid:48)(cid:48) j (cid:48) j (cid:48)(cid:48) k (cid:48) k (cid:48)(cid:48) = n for each of the ( i, j (cid:48) , k (cid:48) )-th oscillatorsand (cid:15) (cid:48) = (cid:15) .Adding and subtracting an additional term (cid:15) (cid:48) ( v s − x i,j (cid:48) ,k (cid:48) ) N (cid:80) i (cid:48)(cid:48) ,j (cid:48)(cid:48) ,k (cid:48)(cid:48) =1 C ii (cid:48)(cid:48) j (cid:48) j (cid:48)(cid:48) k (cid:48) k (cid:48)(cid:48) Γ( x i,j (cid:48)(cid:48) ,k (cid:48)(cid:48) ) = n(cid:15) (cid:48) ( v s − x i,j (cid:48) ,k (cid:48) )Γ( x i,j (cid:48)(cid:48) ,k (cid:48)(cid:48) ) to the x equation of the above system, we get˙ x i,j (cid:48) ,k (cid:48) = ax i,j (cid:48) ,k (cid:48) − x i,j (cid:48) ,k (cid:48) − y i,j (cid:48) ,k (cid:48) − z i,j (cid:48) ,k (cid:48) + n(cid:15) (cid:48) ( v s − x i,j (cid:48) ,k (cid:48) )Γ( x i,j (cid:48) ,k (cid:48) )+ (cid:15) (cid:48) ( v s − x i,j (cid:48) ,k (cid:48) ) N (cid:80) i (cid:48)(cid:48) ,j (cid:48)(cid:48) ,k (cid:48)(cid:48) =1 C ii (cid:48)(cid:48) j (cid:48) j (cid:48)(cid:48) k (cid:48) k (cid:48)(cid:48) (Γ( x i (cid:48)(cid:48) ,j (cid:48)(cid:48) ,k (cid:48)(cid:48) ) − Γ( x i,j (cid:48) ,k (cid:48) ))˙ y i,j (cid:48) ,k (cid:48) = ( a + α ) x i,j (cid:48) ,k (cid:48) − y i,j (cid:48) ,k (cid:48) ˙ z i,j (cid:48) ,k (cid:48) = c ( bx i,j (cid:48) ,k (cid:48) − z i,j (cid:48) ,k (cid:48) + e ) , i = 1 , , ..., N. (17)Introducing the differences between the neural oscillator coordinates X ii (cid:48)(cid:48) ,j (cid:48) j (cid:48)(cid:48) ,k (cid:48) k (cid:48)(cid:48) = x i (cid:48)(cid:48) ,j (cid:48)(cid:48) ,k (cid:48)(cid:48) − x i,j (cid:48) ,k (cid:48) , Y ii (cid:48)(cid:48) ,j (cid:48) j (cid:48)(cid:48) ,k (cid:48) k (cid:48)(cid:48) = y i (cid:48)(cid:48) ,j (cid:48)(cid:48) ,k (cid:48)(cid:48) − y i,j (cid:48) ,k (cid:48) , Z ii (cid:48)(cid:48) ,j (cid:48) j (cid:48)(cid:48) ,k (cid:48) k (cid:48)(cid:48) = z i (cid:48)(cid:48) ,j (cid:48)(cid:48) ,k (cid:48)(cid:48) − z i,j (cid:48) ,k (cid:48) , we can derive the variational equations for transverse stability of the synchronization manifold as˙ X ii (cid:48)(cid:48) ,j (cid:48) j (cid:48)(cid:48) ,k (cid:48) k (cid:48)(cid:48) = (2 ax − x ) X ii (cid:48)(cid:48) ,j (cid:48) j (cid:48)(cid:48) ,k (cid:48) k (cid:48)(cid:48) − Y ii (cid:48)(cid:48) ,j (cid:48) j (cid:48)(cid:48) ,k (cid:48) k (cid:48)(cid:48) − Z ii (cid:48)(cid:48) ,j (cid:48) j (cid:48)(cid:48) ,k (cid:48) k (cid:48)(cid:48) − n(cid:15) (cid:48) Γ( x ) X ii (cid:48)(cid:48) ,j (cid:48) j (cid:48)(cid:48) ,k (cid:48) k (cid:48)(cid:48) + (cid:15) (cid:48) ( v s − x )Γ (cid:48) ( x )( nX ii (cid:48)(cid:48) ,j (cid:48) j (cid:48)(cid:48) ,k (cid:48) k (cid:48)(cid:48) + N (cid:80) p (cid:48) ,q (cid:48) ,r (cid:48) =1 { C i (cid:48)(cid:48) p (cid:48) j (cid:48)(cid:48) q (cid:48) k (cid:48)(cid:48) r (cid:48) X i (cid:48)(cid:48) p (cid:48) ,j (cid:48)(cid:48) q (cid:48) ,k (cid:48)(cid:48) r (cid:48) − C ip (cid:48) j (cid:48) q (cid:48) k (cid:48) r (cid:48) X ip (cid:48) ,j (cid:48) q (cid:48) ,k (cid:48) r (cid:48) } ) , ˙ Y ii (cid:48)(cid:48) ,j (cid:48) j (cid:48)(cid:48) ,k (cid:48) k (cid:48)(cid:48) = 2( a + α ) xX ii (cid:48)(cid:48) ,j (cid:48) j (cid:48)(cid:48) ,k (cid:48) k (cid:48)(cid:48) − Y ii (cid:48)(cid:48) ,j (cid:48) j (cid:48)(cid:48) ,k (cid:48) k (cid:48)(cid:48) , ˙ Z ii (cid:48)(cid:48) ,j (cid:48) j (cid:48)(cid:48) ,k (cid:48) k (cid:48)(cid:48) = c ( bX ii (cid:48)(cid:48) ,j (cid:48) j (cid:48)(cid:48) ,k (cid:48) k (cid:48)(cid:48) − Z ii (cid:48)(cid:48) ,j (cid:48) j (cid:48)(cid:48) ,k (cid:48) k (cid:48)(cid:48) ) . (18)The origin X ii (cid:48)(cid:48) ,j (cid:48) j (cid:48)(cid:48) ,k (cid:48) k (cid:48)(cid:48) = Y ii (cid:48)(cid:48) ,j (cid:48) j (cid:48)(cid:48) ,k (cid:48) k (cid:48)(cid:48) = Z ii (cid:48)(cid:48) ,j (cid:48) j (cid:48)(cid:48) ,k (cid:48) k (cid:48)(cid:48) = 0 is an equilibrium point of the system (18) whose stability determines the stabilityof the synchronization manifold. The first coupling1 FIG. 9: (a) Synchronization threshold (cid:15) ∗ in 3D grid of network with respect to network size, (b) neuronal synchronizationregion in the phase diagram of ( N, (cid:15) ) plane, characterized through global order parameter R . Here N = 30 is the number ofoscillators in one direction of the 3D lattice for which the network size is N . term S = − n(cid:15) (cid:48) Γ( x ) X ii (cid:48)(cid:48) ,j (cid:48) j (cid:48)(cid:48) ,k (cid:48) k (cid:48)(cid:48) depends on thenumber of inputs n , whereas the second couplingterm S = (cid:15) (cid:48) ( v s − x )Γ (cid:48) ( x )( . ) depends on the couplingconfiguration. In terms of the original variable x i,j (cid:48) ,k (cid:48) ,the coupling matrix G = C − n I is the negative of theLaplacian of the connected graph, which has one zeroeigenvalue λ and all other eigenvalues have non-positivereal parts. Then we can write the stability equation onthe synchronous manifold as˙ X = (2 ax − x ) X − Y − Z − Ω( x ) X, ˙ Y = 2( a + α ) xX − Y, ˙ Z = c ( bX − Z ) , (19)where Ω( x ) = n(cid:15) (cid:48) Γ( x ) − (cid:15) (cid:48) ( v s − x )Γ (cid:48) ( x )( n + λ ), λ beingthe eigenvalue with the largest real part. Equation (19)is the same as obtained in [66] which allows us to directlywrite the estimate for synchronization threshold as (cid:15) (cid:48)∗ = (cid:15) ∗ /n , where (cid:15) ∗ = 1 .
285 is the synchronization thresholdfor two synaptically coupled HR neurons. This showsthat the synchrony threshold does not depend on thelattice dimension of the interacting oscillators.Now we have analytically obtained the bound for fullycoherent region through the linear stability analysis. Toverify it with our numerical results, we have plotted thesynchronization threshold (cid:15) (cid:48)∗ against the network size N in Fig. 9(a), where dashed blue line and open red circlelines represent the theoretical and numerical results re-spectively. To numerically calculate the measure of thecoherence level in the 3D locally connected neurons, weuse Kuramoto order parameter R , defined as R = (cid:104) R ( t ) (cid:105) t = (cid:104)| N N (cid:80) i =1 N (cid:80) j =1 N (cid:80) k =1 e ˆ i Φ i,j,k ( t ) |(cid:105) t , (20)where Φ i,j,k is the phase of the ( i, j, k )-th neuron andˆ i = √− R (cid:39) R characterizes thedesynchronized motion of the coupled network. In Fig.9(a), the critical coupling strength (cid:15) (cid:48)∗ , for which R (cid:39) R in the( N, (cid:15) ) parameter space in Fig. 9(b). V. CONCLUSION
In this paper we have studied the existence and emer-gence of chimera patterns in three dimensional grids ofsymmetrically coupled identical oscillators. The nodes ofthe network are assumed to interact nonlinearly in thelocally connected network. In the first part, we haveshown the existence of stationary chimera states in thenetwork of identically coupled Stuart-Landau oscillators.Using the OA approach, we have investigated the chimerapatterns in the continuum limit corresponding to infinitenumber of oscillators and it is observed that the local or-der parameter and the local phase of the oscillators, ob-tained through the analysis of OA reduction technique,validate the numerical results. Then, by considering arealistic neuronal interacting medium as chemical synap-tic coupling function, we enunciated the emergence ofchimera patterns in the 3D grid of neuronal network. Themathematical form of the chemical synapse is of the sig-moidal type nonlinear function. Hindmarsh-Rose burst-ing neurons are considered as the local dynamics of theneuronal network. Here the important observation is theemergence of non-stationary chimera patterns. The com-bination of slow-fast dynamical features of each neuronand the nonlinear interaction function among them leads2to the appearance of such chimera patterns. Using dif-ferent characterizing quantities, such as strength of inco-herence and local order parameter, we have characterizedvarious collective dynamical states and articulated thetransition scenarios among them. We can say that thenonlinearity present in the interaction function acts as anessential condition for the existence of chimera states inour considered network. Further, we have also calculatedthe synchronization threshold for HR network interactingin 3D lattice architecture. Our finding reveals that thesynchrony threshold is independent of the lattice dimen-sion, which we consider to be an interesting observation.The central aim of the present work is to establish thegeneralization of the dynamics of one dimensional andtwo dimensional networks to three dimensional networks,truly inspired from the complex connectivity structure of the brain network. Our findings may be generalized inmore 3D complex structures which are related to con-nectivity formation in the brain. The present study maygive a better knowledge to understand the different neu-ronal disorder diseases. The effect of time-delay in thelocal coupling function with different 3D network topolo-gies can be extended for further studies. Also how doesthe increase in the lattice dimension affect the observednonlinear phenomena will be an interesting question todiscuss in the near future.
Acknowledgments:
DG was supported by DST-SERB (Department of Science and Technology), Gov-ernment of India (Project no. EMR/2016/001039). Thework of ML is supported by a DST-SERB DistinguishedFellowship program. [1] M. J. Panaggio, and D. M. Abrams,
Nonlinearity ,R67 (2015).[2] B. K. Bera, S. Majhi, D. Ghosh and M. Perc, Europhys.Letts. , 10001 (2017).[3] S. Majhi, B. K. Bera, D. Ghosh and M. Perc,
Phys. Lif.Rev. (doi:10.1016/j.plrev.2018.09.003), (2018).[4] Y. Kuramoto, and D. Battogtokh,
Nonlinear Phenom.Complex Syst. , 380 (2002).[5] D. M. Abrams, and S. H. Strogatz, Phys. Rev. Lett. ,174102 (2004).[6] M. R. Tinsley, S. Nkomo, and K. Showalter, Nat. Phys. , 662 (2012).[7] M. Wickramasinghe, I. Z. Kiss, PLoS ONE , e80586(2013).[8] L. V. Gambuzza, A. Buscarino, S. Chessari, L. Fortuna,R. Meucci, and M. Frasca, Phys. Rev. E , 032905(2014).[9] A. M. Hagerstrom, T. E. Murphy, R. Roy, P. H¨ovel, I.Omelchenko, and E. Sch¨oll, Nat. Phys. , 658 (2012).[10] E. A. Martens, S. Thutupalli, A. Fourri`ere, and O. Hal-latschek, Proc. Natl. Acad. Sci. USA , 10563 (2013).[11] L. Larger, B. Penkovsky, and Y. Maistrenko,
Phys. Rev.Lett. , 054103 (2013);
Nat. Commun. , 7752 (2015).[12] N. C. Rattenborg, C. J. Amlaner, and S. L. Lima, Neu-rosci. Biobehav. Rev. , 817 (2000).[13] N. C. Rattenborg, Naturwissenschaften , 413 (2006).[14] F. Dorfler, M. Chertkov, and F. Bullo, Proc. Natl. Acad.Sci. USA , 2005 (2013).[15] A. E. Motter, S. A. Myers, M. Anghel, and T. Nishikawa,
Nat. Phys. , 191 (2013).[16] N. Lazarides, G. Neofotistos, and G. P. Tsironis, Phys.Rev. B , 054303 (2015); J. Hizanidis, N. Lazarides,and G. P. Tsironis, Phys. Rev. E , 032219 (2016).[17] J. M. Davidenko, A. V. Pertsov, R. Salomonsz, W. Bax-ter, and J. Jalife, Nature , 349 (1992).[18] A. V. Panfilov,
Chaos , 57 (1998).[19] R. Levy, W. D. Hutchison, A. M. Lozano, and J. O.Dostrovsky, J. Neurosci. , 7766 (2000).[20] G. F. Ayala, M. Dichter, R. J. Gumnit, H. Matsumoto,and W. A. Spencer, Brain Res. , 1 (1973).[21] C. Gu, G. St-Yves, and J. Davidsen, Phys. Rev. Lett. , 134101 (2013). [22] R. Gopal, V. K. Chandrasekar, A. Venkatesan, and M.Lakshmanan,
Phys. Rev. E , 052914 (2014).[23] S. Ulonska, I. Omelchenko, A. Zakharova, and E. Sch¨oll, Chaos , 094825 (2016).[24] B. K. Bera, D. Ghosh, P. Parmananda, G. V. Osipov,and S. K. Dana, Chaos , 073108 (2017).[25] I. Omelchenko, Y. Maistrenko, P., H¨ovel, and E. Sch¨oll, Phys. Rev. Lett. , 234102 (2011).[26] S. Rakshit, B. K. Bera, M. Perc and D. Ghosh,
Sci. Rep. , 2412 (2017).[27] J. Hizanidis, V. Kanas, A. Bezerianos, and T. Bountis, Int. J. Bifurcat. Chaos , 1450030 (2014).[28] A. Yeldesbay, A. Pikovsky, and M. Rosenblum, Phys.Rev. Lett. , 144103 (2014).[29] V. K. Chandrasekar, R. Gopal, A. Venkatesan, and M.Lakshmanan,
Phys. Rev. E , 062913 (2014).[30] A. Mishra, C. Hens, M. Bose, P. K. Roy, and S. K. Dana, Phys. Rev. E , 062920 (2015).[31] G. C. Sethia, and A. Sen, Phys. Rev. Lett. , 144101(2014).[32] F. B¨ohm, A. Zakharova, E. Sch¨oll, and K. L¨udge,
Phys.Rev. E , 040901(R) (2015).[33] L. Schmidt, and K. Krischer, Phys. Rev. Lett. ,034101 (2015).[34] L. Schmidt, and K. Krischer,
Chaos , 064401 (2015).[35] C. R. Laing, Phys. Rev. E , 050904(R) (2015).[36] B. K. Bera, and D. Ghosh, Phys. Rev. E , 052223(2016).[37] J. Hizanidis, N. Lazarides, and G. P. Tsironis, Phys. Rev.E , 032219 (2014).[38] B. K. Bera, D. Ghosh, and M. Lakshmanan, Phys. Rev.E , 012205 (2016).[39] Y. Zhu, Z. Zheng, and J. Yang, Phys. Rev. E , 022914(2014).[40] J. Hizanidis, N. Lazarides, and G. P. Tsironis, Chaos ,013113 (2009).[41] S. Majhi, M. Perc, and D. Ghosh, Sci. Rep. , 39033(2016).[42] V. A. Maksimenko, V. V. Makarov, B. K. Bera, D.Ghosh, S. K. Dana, M. V. Goremyko, N. S. Frolov, A.A. Koronovskii, A. E. and Hramov, Phys. Rev. E ,052205 (2016). [43] S. Ghosh, and S. Jalan, Int. J. Bifur. Chaos , 1650120(2016).[44] S. Ghosh, A. Kumar, A. Zakharova, and S. Jalan, Euro-phys. Letts. , 60005 (2016).[45] S. Majhi, M. Perc, and D. Ghosh,
Chaos , 073109(2017).[46] A. Buscarino, M. Frasca, L. V. Gambuzza, and P. H¨ovel, Phys. Rev. E , 022817 (2015).[47] J. Hizanidis, N. E. Kouvaris, G. Zamora-L´opez, A. D´ıaz-Guilera, and C. G. Antonopoulos, Sci. Rep. , 19845(2016).[48] V. V. Makarov, S. Kundu, D. V. Kirsanov, N. S. Frolov,V. A. Maksimenko, D. Ghosh, S. K. Dana, and A. E.Hramov, Commun. Nonlinear Sci. Numer. Simul. ,118 (2019).[49] K. Premalatha, V. K. Chandrasekar, M. Senthilvelan,and M. Lakshmanan, Phys. Rev. E , 012311 (2016).[50] K. Premalatha, V. K. Chandrasekar, M. Senthilvelan,and M. Lakshmanan, Phys. Rev. E , 022208 (2017).[51] S. Majhi, and D. Ghosh, Chaos , 083113 (2018).[52] T. Kapitaniak, P. Kuzma, J. Wojewoda, K. Czolczynski,and Y. Maistrenko, Sci. Rep. , 6379 (2014).[53] J. Xie, E. Knobloch, and H. C. Kao, Phys. Rev. E ,022919 (2014).[54] B. K. Bera, D. Ghosh, and T. Banerjee, Phys. Rev. E , 012215 (2016).[55] B. W. Li, and H. Dierckx, Phys. Rev. E , 020202(R)(2016).[56] S. Kundu, S. Majhi, P. Muruganandam, D. Ghosh, Eur.Phys. J. Spec. Top. , 983 (2018).[57] A. Zakharova, M. Kapeller, and E. Sch¨oll,
Phys. Rev.Lett. , 154101 (2014).[58] G. C. Sethia, A. Sen, and G. L. Johnston,
Phys. Rev. E , 042917 (2013).[59] A. Schmidt, T. Kasimatis, J. Hizanidis, A. Provata, andP. H¨ovel, Phys. Rev. E , 032224 (2017).[60] S. Kundu, S. Majhi, B. K. Bera, D. Ghosh, and M. Lak-shmanan, Phys. Rev. E , 022201 (2018).[61] T. Kasimatis, J. Hizanidis, and A. Provata, Phys. Rev.E , 052213 (2018).[62] L. Janagal and P. Parmananda, Phys. Rev. E , 056213(2012).[63] E. Ott and T. M. Antonsen, Chaos , 037113 (2008).[64] E. Ott and T. M. Antonsen, Chaos , 023117 (2009).[65] J. L. Hindmarsh and M. Rose, Proc. R. Soc. London,Ser. B , 87 (1984).[66] I. Belykh, E. de Lange, M. Hasler,
Phys. Rev. Lett.94