Table Of Content

Chimera patterns in two-dimensional networks of coupled neurons Alexander Schmidt,1 Theodoros Kasimatis,2 Johanne Hizanidis,2,3 Astero Provata,2 and Philipp Hövel1,4 1Institut fu¨r Theoretische Physik, Technische Universita¨t Berlin, Hardenbergstraße 36, 10623 Berlin, Germany 2Institute of Nanoscience and Nanotechnology, National Center for Scientific Research “Demokritos”, 15310 Athens, Greece 3Crete Center for Quantum Complexity and Nanotechnology, Department of Physics, University of Crete, 71003 Heraklion, Greece 4Bernstein Center for Computational Neuroscience Berlin, Humboldt-Universität zu Berlin, Philippstraße 13, 10115 Berlin, Germany (Dated: Received: January 3, 2017/ Revised version: date) We discuss synchronization patterns in networks of FitzHugh-Nagumo and Leaky Integrate-and- Fire oscillators coupled in a two-dimensional toroidal geometry. Common feature between the two modelsisthepresenceoffastandslowdynamics,atypicalcharacteristicofneurons. Earlierstudies have demonstrated that both models when coupled nonlocally in one-dimensional ring networks 7 produce chimera states for a large range of parameter values. In this study, we give evidence of 1 a plethora of two-dimensional chimera patterns of various shapes including spots, rings, stripes, 0 and grids, observed in both models, as well as additional patterns found mainly in the FitzHugh- 2 Nagumo system. Both systems exhibit multistability: For the same parameter values, different n initialconditionsgiverisetodifferentdynamicalstates. Transitionsoccurbetweenvariouspatterns a whentheparameters(couplingrange,couplingstrength,refractoryperiod,andcouplingphase)are J varied. Many patterns observe in the two models follow similar rules. For example the diameter of the rings grows linearly with the coupling radius. 2 ] PACSnumbers: 05.45.Xt,89.75.-k,87.19.lj O Keywords: synchronization,spatio-temporalpatternformation,nonlineardynamics,chimerastates A . I. INTRODUCTION couplingstrengthandrange. Thesemulti-chimerastates n i were shown to be robust when small inhomogeneities in l n the coupling topology (with identical elements) or in- Chimerasarehybridstatesthatemergespontaneously, [ homogeneous elements with regular nonlocal coupling combining both coherent and incoherent parts [1]. First were introduced [26]. For a constant coupling strength 1 foundinidenticalandsymmetricallycoupledphaseoscil- and given number of links, hierarchical connectivity was v lators[2],chimerastateshavebeenthefocusofextensive shown to induce nested multi-chimera patterns [26]. Co- 2 research for over a decade now. Both theoretical and ex- 4 existence of coherent and incoherent domains were also perimental works have shown that this counter-intuitive 3 observedinsystemsofexcitableFHNelementsunderthe collective phenomenon may arise in numerous systems 0 influence of noise. As Semenova et al. stress, the noise 0 includingmechanical,chemical,electro-chemical,electro- has often a constructive role: It shifts the dynamics of . optical, electronic, and superconducting coupled oscilla- 1 identical excitable elements into the oscillatory domain, tors [3–14]. 0 giving then rise to chimera states [27]. Moreover, it is 7 The phenomenon of chimera states has also been ad- possible to control the position of coherent and inco- 1 dressed in networks of biological neural oscillators. In herent domains of a multi-chimera state by introducing v: particular, Hindmarsh-Rose neural oscillators have been a block of excitable FHN elements in appropriate posi- i studied in networks with nonlocal [15] and nearest- tions, or to generate a chimera directly from the syn- X neighbor [16] coupling, as well as in modular networks chronous state [28]. This observation offers promising ar consisting of communities [17]. Potential relevance of ideasintermsofachievingdesiredstatesbyalocalmod- chimerastatesinthiscontextincludebumpstates[18,19] ification of the system parameters. For instance, it may andthephenomenonofunihemisphericsleepobservedin be possible in targeted medication to modify incoherent birds and dolphins [20], which sleep with one eye open, neurondynamicsbychangingonlythepotentialofafew meaning that half of the brain is synchronous with the neurons locally, leaving the coupling topology of the rest other half being asynchronous. Furthermore, it has been of the network untouched. In addition, Omelchenko et recently hypothesized that chimera states are the route al. proposed an alternative control scheme that extends of onset or termination of epileptic seizures [21–24]. the lifetime of chimeras, which are known to be chaotic Chimera states have also been reported in the one- transients [29], and at the same time reduces the erratic dimensional FitzHugh-Nagumo (FHN) model with non- drift present in small networks [30]. Controlling the po- local connectivity [25]. Patches of synchrony were ob- sitioncanalsobeachievedbyintroducinganasymmetric servedwithintheincoherentdomainsgivingrisetomulti- coupling strength [31]. chimera states, when the coupling constant increased Experimentally, the FHN model has been studied by above the weak limit. The multiplicity of the state, that Essaki Arumugam and Spano, in connection to synchro- is, the number of (in)coherent regions, depended on the nization phenomena associated with neurological disor- 2 ders such as epilepsy [32]. In their study, they imple- two systems. In Sec. II, the two models are briefly mented nine FHN neurons linked in a ring topology, via recapitulated. In Sec. III, main attributes of spot and discrete electronics. Introducing nonlocal coupling, a ring chimeras are presented for the two models, and chimera state appeared, while for local connectivity ei- their common and different properties are discussed. ther fully synchronous or asynchronous states were ob- Similarly, in Secs. IV and V stripe and grid chimeras served. Their results indicate that epilepsy can be un- are discussed, respectively. Additional miscellaneous derstood as a topological disease, strongly related to the chimerasaswellasotherpatternsaresummarizedinthe connectivity of the underlying network of neurons. appendices. Finally, the main results and open problems Chimera states were recently reported in the one- are discussed in a brief concluding section. dimensional Leaky Integrate-and-Fire (LIF) model, a system describing the spiking behavior of neuron cells. II. THE MODELS It was shown that identical LIF oscillators nonlocally linked in a one-dimensional ring geometry give rise to In this section, we summarize the two models con- multi-chimerastateswhosemultiplicitydependsbothon sidered in this article – FitzHugh-Nagumo (FHN) and the coupling strength and the refractory period of the Leaky Integrate-and-Fire (LIF) model – and introduce neuron cells [33]. In analogy with the FHN model, the the respective toroidal coupling schemes describing the introduction of a hierarchical topology in the coupling two-dimensional layout with nonlocal interactions. induced nested chimera states and also transitions between multi-chimera states with different multiplicities. Usinga differentgeometrical setupof twopopulationsof A. The FitzHugh-Nagumo model identical LIF oscillators coupled via excitatory coupling, Olmietal. studiedtheonsetofchimerasaswellasstates characterized by a different degree of synchronization in ThedynamicsofasingleFHNsystemaredescribedby the two populations [34]. Irregular synchronization phe- the following set of equations [40, 41]: nomena have been reported for LIF elements even in the case of all-to-all connectivity [35]. dx x3 (cid:15) =x− −y (1a) Chimera states have mainly been investigated in one- dt 3 dimensional systems. Recently, works involving two [36, dy =x+a, (1b) 37] and three-dimensional [38] oscillator arrays have re- dt vealed new types of chimera states, depending on the coupling function. These studies have mainly focused where x and y denote a fast activator and slow inhibitor on phase oscillators. In this paper, we go beyond the variable, respectively. The parameter (cid:15) determines the simple Kuramoto model and consider a two-dimensional timescaleseparationandisfixedinthisstudyat(cid:15)=0.05. network configuration using two different models related The threshold parameter a determines the oscillatory to neuronal spiking activity. This topology is motivated (|a| < 1) or excitable (|a| > 1) behavior in the system, frommedicalexperimentswherethinbrainslicesarecul- i.e., an unstable or stable fixed point at the intersection tured in Petri dishes and various electrical and chemical of the nullclines, and is set to a = 0.5 throughout this properties are recorded, see for example Ref. [39]. The study. Theinteractionofthetwosystemvariablesofeach two-dimensional nonlocal connectivity studied here can node in the network is realized via a rotational coupling be considered as an approximation of the acute brain matrix as proposed in Ref. [25] and detailed below. slices, whose connectivity is certainly more complex. WeconsiderthefollowingcouplingschemefortheFHN Our investigation follows a parallel presentation of model: A two-dimensional regular N ×N-network with common patterns in the two-dimensional FHN and N = 100 nodes and periodic, toroidal boundary condi- LIF systems highlighting the main features that are tions. Thissettingyieldsthefollowingsystemofnetwork responsible for the formation of each pattern in the equations: (cid:15)dxij =x − x3ij −y + σ (cid:88) [b (x −x )+b (y −y )] (2a) dt ij 3 ij N −1 xx ij mn xy ij mn r (m,n)∈Br(i,j) dyij =x +a+ σ (cid:88) [b (x −x )+b (y −y )], (2b) dt ij N −1 yx ij mn yy ij mn r (m,n)∈BFHN(i,j) r where we fix the coupling strength at σ = 0.1 in the following. Note that the double index of the dynamic 3 variables x and y with i,j = 1,...,N refers to the potential. Thevariableuevolvesaccordingtothefollow- ij ij position on the two-dimensional lattice. All oscillators ing equation: (x ,y )areidenticalandcoupleisotropicallytoallother ij ij du oscillators in a circular neighborhood given by: =µ−u. (6) dt BFHN(i,j)={(m,n):(m−i)2+(n−j)2 ≤r2}, (3) r Ifthemembranepotentialureachesathresholdu <µ, th where the differences need to be calculated in a man- it is reset to a resting potential: ner to account for a toroidal geometry. In other words, (m−i)2 +(n−j)2 is computed to obtain the shortest limu(t+ε)→u , when u≥u . (7) rest th distanceonthetorus. SeeFig.1(a)foraschematicdepic- ε→0 tion. Hence, the number of lattice points inside a circle Without loss of generality, the reset or resting potential of radius r is given by: can be set to u =0. rest Nr =1+4(cid:88)∞ (cid:18)(cid:22)4ir+2 1(cid:23)−(cid:22)4ir+2 3(cid:23)(cid:19), (4) µ −Th(eµ−soluu0t)ioen−towfitEhq.th(e6)iniistiatlhecnondgiivtieonn buy(0)u(=t) u=0 and is repeated after each reset starting again at u = i=0 0 u . The period between two resets T of a single el- where (cid:98)·(cid:99) is the floor function, which gives the largest rest s ement depends both on µ and u as follows: T = integer less than or equal to its argument. th s ln[(µ−u )/(µ−u )]. 0 th As an extension to the standard LIF model, we take a (cid:0)(cid:1)(cid:0)(cid:1) (cid:0)(cid:1) (cid:0)(cid:1)(cid:0)(cid:1) (cid:0)(cid:1) (cid:0)(cid:1) (cid:0)(cid:1) refractory period p into account via the following con- (cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1) (cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1) r j+r j+R dition: (cid:0)(cid:1)(cid:0)(cid:1) (cid:0)(cid:1) (cid:0)(cid:1)(cid:0)(cid:1) (cid:0)(cid:1) (cid:0)(cid:1) (cid:0)(cid:1) (cid:0)(cid:1)(cid:0)(cid:1) (cid:0)(cid:1) (cid:0)(cid:1)(cid:0)(cid:1) (cid:0)(cid:1) (cid:0)(cid:1) (cid:0)(cid:1) (cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1) (cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1) u(t)=urest,∀t:[n(Ts+pr)+Ts]≤t≤[(n+1)(Ts+pr)] j j (8) (cid:0)(cid:1)(cid:0)(cid:1) (cid:0)(cid:1) (cid:0)(cid:1)(cid:0)(cid:1) (cid:0)(cid:1) (cid:0)(cid:1) (cid:0)(cid:1) (cid:0)(cid:1)(cid:0)(cid:1) (cid:0)(cid:1) (cid:0)(cid:1)(cid:0)(cid:1) (cid:0)(cid:1) (cid:0)(cid:1) (cid:0)(cid:1) with n = 0,1,2,.... In other words, the LIF element is j−r (cid:0)(cid:1)(cid:0)(cid:1) (cid:0)(cid:1) (cid:0)(cid:1)(cid:0)(cid:1) j−R (cid:0)(cid:1) (cid:0)(cid:1) (cid:0)(cid:1) (cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1) (cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1) held at the rest potential during the refractory period. (cid:0)(cid:1)(cid:0)(cid:1) (cid:0)(cid:1) (cid:0)(cid:1)(cid:0)(cid:1) (cid:0)(cid:1) (cid:0)(cid:1) (cid:0)(cid:1) As a result, the total period T of the single element is (cid:0)(cid:1)(cid:0)(cid:1) (cid:0)(cid:1) (cid:0)(cid:1)(cid:0)(cid:1) (cid:0)(cid:1) (cid:0)(cid:1) (cid:0)(cid:1) (a) i−r i i+r (b) i−R i i+R now: T =T +p =ln[(µ−u )/(µ−u )]+p . s r 0 th r Similar to the case of the FHN model, we consider the FIG. 1: (Color online) Schematic depiction of the coupling LIFnetworkdynamicsonatwo-dimensionalregularN× topologyusedfortheFHNmodelofEqs.(2)inpanel(a)and N-network with N = 100 nodes and toroidal boundary fortheLIFmodelofEq.(9)inpanel(b)accordingtoEqs.(3) conditions. The LIF coupled network dynamics read: and (10), respectively. duij =µ−u + σ (cid:88) [u −u ], (9) Thecouplingbetweenthex-andy-variablesisrealized dt ij N −1 ij mn R by a rotational coupling matrix using the single param- (m,n)∈BRLIF(i,j) eter ϕ ∈ [0,2π), which allows diagonal or direct cou- where the neighborhood of the element u with i,j = ij pling (ϕ = 0 and ϕ = π), cross coupling (ϕ = π/2 and 1,...,N is given by ϕ = 3π/2) and mixed coupling scenarios. Accordingly, the coupling matrix is given by: BLIF(i,j)={(m,n):i−R≤m≤i+R∧j−R≤n≤j+R}, R (cid:18) (cid:19) (cid:18) (cid:19) (10) b b cosϕ sinϕ B = xx xy = . (5) which includes all elements on the grid within a square b b −sinϕ cosϕ yx yy of side 2R+1 as depicted in Fig. 1(b). The size of the For this type of coupling, the coupling phase ϕ corre- coupled region is then given by N =(2R+1)2. R sponds to the so-called phase frustration parameter α in Tocalculatecoherenceinbothmodelsweusethemean the Kuramoto system, which is known to be a crucial phase velocity ω of an oscillator at position (i,j) [25]. ij quantity for the occurrence of chimera states [42]. As If c (∆t) is the number of periods that the oscillator at ij derivedinRef.[25],themappingfromϕtoαisgivenvia position (i,j) has completed in a time interval ∆t, the a phase-reduction technique. mean phase velocity is defined as: ω =2πc (∆t)/∆t. ij ij In the following, we will use the coupling phase ϕ and Inthenextsections,weprovideaseriesofselectedpat- coupling radius r to explore the dynamical scenarios in terns that we observe for both the FHN and LIF model. the networked system. In order to explore the parameter space, we use random initial conditions. Whenever we observe a characteristic pattern for a specific choice of coupling parameters, we B. The Leaky Integrate-and-Fire model performacontinuationbyscanningthevicinityusingthe patternoftheprevioussimulationasaninitialcondition. TheLIFmodeldescribesaneuronviaasingledynam- Thisway,weareabletomapouttheregionsofexistence ical variable u that can be interpreted as its membrane of the different dynamical patterns. 4 III. SPOTS AND RING PATTERNS In this section, we focus on chimera states of spot and ring type in both FHN and LIF models. A. FitzHugh-Nagumo: Spots and ring patterns In the 100×100 FHN system described by Eqs. (2) and (3), we observe various chimera regions in the parameter space spanned by the coupling radius r ∈[1,49] and coupling phase ϕ ∈ [0,2π). A complete scan of the parameter space and several enlargements can be found in Appendix A. Here, we will focus on the range FIG. 2: (Color online) FHN system: Chimera state of (a) ϕ∈[1.3,1.65]. incoherent-spot and (b) coherent-spot type. The left panels An example for a single-headed incoherent spot sur- depict a snapshot of the activator variable x at t = 2000. ij roundedbycoherentoscillatorsisshowninthetoppanel Thecenterandrightpanelsshowthemeanphasevelocityω ij of Fig. 2(a). The left, center, and right panels depict a andahorizontalcutofω . Couplingparameters: (a)r=33, ij snapshot of xij, the mean phase velocity ωij, where the ϕ = π/2−0.2 and (b) r = 42, ϕ = π/2. Other parameters: average is computed over 1000 time units, and a hori- σ=0.1, a=0.5, and (cid:15)=0.05. zontal cut of ω through the center of the spot, respec- ij tively. The latter exhibits the typical arc-shape profile of a chimera state. When the spots deform to a square shape, the maximum of the ω ’s is less pronounced and ij becomes a plateau (not shown). The spots and the squares grow when the coupling radius r increases, as shown in Fig. 3(a) for all observed spots. The color code refers to different choices of ϕ. We observe two different domains of large and small coupling radii. These correspond to different areas in parameter space, which support spot chimeras, see Appendix A. In the range 20 ≤ r ≤ 30 the spot patterns lose their stability in space, wander around the grid in an unpredictable way, changetheirshapeovertimeorevensplitintopiecesand reemerge after some time. Anotherchimeraspotpattern–againfoundasasingle- headedandmulti-headedversion–isaspotwithacoher- ent center. An example is shown in the bottom panel of Fig. 2(b). As in the previous example, the ω -values are ij lower in the coherent area than in the incoherent region (middle panel). Furthermore, the section of the mean phase velocities (right panel) has a shape similar to incoherent spot type and many other chimera patterns in one dimension. This state can also be found with two coherent heads (not shown). The size of the coherent region is depicted in Fig. 3(b). We observe only a weak dependence on the coupling radius r, while the coupling phase ϕ has a stronger effect. Compared to the incoherent spot chimeras, the size increases with decreasing ϕ. In some cases, we observe that at the opposite side of the coherent region of the torus, a spiral can form, see FIG.3: (Coloronline)FHNsystem: Radiusρofthespot-like Appendix B for an example. The radius of the spiral chimera states, (a) incoherent spot and (b) coherent spot, increases with the coupling radius (not shown). in dependence on the coupling radius r for different coupling Between regions of spots (see Appendix A) and fully phasesϕgivenascolorcode. Theerrorbarsindicateanuncer- synchronoussolutions,weobservethatringchimerascan taintyof(a)∆ρ=0.3and(b)∆ρ=0.5inthemeasurements be found for ϕ < 1.35, see Fig. 4. The label “ring of the radii. Parameters as in Fig. 2. chimera” refers to the dip in the mean phase velocity ω in the center of the incoherent ring that can reach ij 5 the level of the coherent oscillators or have a different, Theregioninparameterspacediscussedaboveincludes larger ω -value as in the displayed example. We also many other patterns, like square- or cross-shaped inco- ij find that rings can occur as a single-headed state and a herent spot, stripe (cf. Sec. IV), alternating and more multi-headed state with 2 rings existing simultaneously complexchimeras,manyoftheminamultistableconfigu- atdifferentpositions(notshown). Similartotheincoher- ration. ExampleswillbediscussedinAppendixB.More- ent spots, the rings can either be round or have a square over, inthelimitr →0whichcorrespondstothecontin- shape. Figure 5 shows the dependence of the inner and uous reaction-diffusion FHN model, the well-known Tur- outer diameter of the ring chimeras on the coupling ra- ing patterns can be found. In addition, we find grid-like dius r using different levels of brightness to indicate the structuresandlinesofincoherentspotsinadifferentarea number of occurrences on the considered ϕ-range. We of the parameter space (cf. Sec. V). findthattheouterdiameterincreaseslinearwithr,while thesizeoftheinnerringdoesnotshowsuchastrongde- pendence. B. Leaky Integrate-and-Fire: Spots and ring patterns Spots and ring patterns are also observed in the LIF model given by Eqs. (9) and (10), in particular for small valuesofthecouplingstrengthσ. Thespotpatternexists for p =0, while for finite refractory periods the interior r of the spot synchronizes and ring patterns are obtained. As in the case of the FHN model we also use a N×N = 100×100latticewithtoroidalboundaryconditions,while FIG.4: (Coloronline)FHNsystem: Ringchimera,foundfor herethecontrolparametersarethecouplingradiusRand r = 33 and ϕ = π/2−0.24. A snapshot at t = 1900 of the the refractory period p . (color coded) activator variable x is shown on the left side r ij Figure 6 depicts an incoherent-spot chimera state for and the mean phase velocity and a section through the ring inthecenterandrightpanels,respectively. Otherparameters small values of the coupling constant and zero refractory as in Fig. 2. period. A typical snapshot of the LIF spot chimera is shown in the left panel, the corresponding mean phase velocityω isdepictedinthemiddleandahorizontalcut ij crossing the spot is depicted on the right. Single spots developintheLIFsystemforsmallvaluesofσ (hereσ = 0.1)whenthecouplingrangeRissmall. Asthecoupling range increases while keeping the refractory period at zerothespotbreaksintoseveralunequalsecondaryspots whosesizeanddistanceincreasewiththecouplingrange (not shown). FIG. 6: (Color online) LIF system: Chimera state of incoherent-spot type. Snapshot of u (left) and the corre- ij sponding mean phase velocity (middle and right) for small values of the coupling constant σ = 0.1 and for a coupling FIG. 5: (Color online) FHN system: Inner (red triangles) range R = 10. Other parameters: µ = 1.0, p = 0 and and outer (black squares) diameter D of the incoherent ring r N =100. in dependence on the coupling radius r. The error bars indicate an uncertainty of ∆D = 0.9 in the measurements of Ring patterns are produced, when the refractory pe- the diameters. The lines are linear fits as a guide to the eye. All ring patterns found for ϕ ∈ [π/2−0.24,π/2−0.2] are riod takes finite values pr > 0 and the interior of the shown. The inset shows the data with rescaled axes: Frac- spots synchronizes. The snapshot of uij, ωij and its cor- tion of incoherent elements N /N2 vs. number of coupled responding horizontal cut of a ring chimera is shown in inc neighbors Nr/N2. The brightness of the dots corresponds to Fig. 7. The outer border of the ring has a square shape the number of observed ring patterns for every (D,r)-pair in due to the square configuration of the connectivity ma- the considered ϕ-range. Other parameters as in Fig. 2. trixgivenbyEq.(10). Thispatterndiffersfromtheones 6 appearing in the one-dimensional configuration: From everypointoftheincoherentringonecanreachallother incoherent elements by moving continuously on the ring, and this region forms a set that is not simply connected and cannot exist in one dimension. As the size of the coupling range increases, the ring radius increases accordingly. When the coupling range approaches the system size, the single ring gives rise to multiple parallel lines of increasing size (not shown). In the limit of small coupling ranges (e.g. R = 1), we do not observe the formation of coherent and incoherent regions. Instead, Turing-like shapes are formed as in the case of the FHN model. u ω ij ij 100 1 100 1.52 1.52 80 0.8 80 FIG. 8: (Color online) LIF system: Inner and outer diam- 1.5 1.5 j 4600 00..46j 4600 j=15ωi,j ecoteurplDingofratnhgeeinRc.ohInerseent:t Rrienlgatpivaettneurnmsbaesr aoffiunnccothioernenotf etlhee- 1.47 1.47 mentsN /N2asafunctionoftherelativenumberofcoupled 20 0.2 20 inc neighbors N /N2. Parameters as in Fig. 7. 0 0 20 40 60 80 100 0 0 0 20 40 60 80 100 1.45 0 20 40 60 80 100 1.45 R i i i FIG.7: (Coloronline)LIFsystem: Ringchimera: Snapshot ofu (left)andthecorrespondingmeanphasevelocity(mid- ij dleandright)foracouplingrangeofR=10andp =0.22T . IV. STRIPES r s Other parameters as in Fig. 6. In this section, we focus on a different type of chimera states that exhibit one-dimensional structures extended Theinnerandouterdiameteroftheincoherentregions in the second spatial dimension: stripes. asafunctionofthecouplingrangeRisplottedinFig.8. Thesizeofbothdiametersincreasesproportionallytothe couplingrangeasthelinearfits(solidlines)indicate. The inset of the same figure depicts the relative number of incoherent elements N /N2 (where N is the number inc inc of oscillators in the incoherent domain) as a function of A. FitzHugh-Nagumo: Stripes the relative connectivity matrix size N /N2 = (2R + R 1)2/N2. From this figure we observe that the growth rate of Ninc/N2 is almost linear. We also find, that the Besides the two spot types and rings, we also iden- mean phase velocities are independent of the coupling tified stripe chimeras in the FHN model, as shown in range (not shown). As we will see in the next sections, Fig.9. Thesestatesconsistofanincoherentandacoher- the mean phase velocities are mainly influenced by the ent stripe region. The latter oscillates with lower mean coupling constant σ. phasevelocityω . Again,theω -profileofthestripeex- ij ij hibits an arc-shape. In the depicted example the stripe region is synchronous, but in most cases it contains a gradual phase shift, which means that waves travel in C. Comparison the direction of the stripe (see Appendix B). For this synchronized stripe type, we never observed a diagonal stripe pattern. We also found other stripe patterns that Comparing the FHN oscillator and LIF dynamics, we exhibit a wave-like dynamics within the stripe instead of identifysimilarspatialstructureswithspotorringshape. complete synchronization. Their size grows linearly with the coupling range. It is worth stressing that the ω -values in the incoherent re- Figure 10 depicts the width W of the coherent re- ij gions are larger than in the coherent ones. For small gion as a function of the coupling range r for the fully coupling ranges, i.e. the local coupling regime, we find synchronized stripe type. The color code refers to dif- Turing-like behavior. In the LIF model, the oscillators ferent coupling phases. We find that the width of the synchronize for large coupling ranges close to all-to-all coherent region increases for larger r and that smaller connectivity,andintheFHNmodeltheysynchronizefor phases result in larger values of W. Other multi-stripe coupling phases slightly smaller than in the area where chimeras are found outside the previously discussed re- spots and rings are observed. gion of ϕ∈[1.3,1.65] (see Appendix B). 7 not been observed. u ω ij ij 100 1 100 2.16 2.16 80 0.8 80 j 4600 00..46 j 4600 22..0182 22..0182j=35ωij 20 0.2 20 0 0 0 2.04 2.04 0 20 40 60 80 100 0 20 40 60 80 100 0 20 40 60 80 100 FIG. 9: (Color online) FHN system: Stripe chimera for i i i r = 39, ϕ = π/2−0.04, and t = 1900. The (color coded) FIG. 11: (Color online) LIF system: Snapshot of u (left), activator x (left panel) shows a synchronous stripe region, ij ij thecorrespondingmeanphasevelocity(middle)andasection whiletheotherpanelsshowthemeanphasevelocitydistribu- of it (right). Parameters are σ = 0.6, µ = 1.0, p = 0.6T , tion ω (middle panel) and the section of ω (right panel). r s ij ij N =100, and R=20. Other parameters as in Fig. 2. C. Comparison StripechimerashavebeenfoundinbothFHNandLIF models. They include single and multiple stripes. In the FHN system, we have observed an arc-shaped mean phase velocity profile. The width of the coherent region becomes smaller as the coupling phase increases. For multiple stripes in the LIF model, coherent regions emerge with different mean phase velocities, while ω ij was found to be equal for every coherent stripe of the multi-headed stripe chimeras in the FHN system. V. GRID PATTERNS FIG.10: (Coloronline)FHNsystem: WidthW ofthecoher- ent region in the fully synchronized stripes in dependence on the coupling range r for different coupling phases ϕ as color Grids are the most common patterns observed in both code. The size of the error bars refers to an uncertainty of FHNandLIFmodels. Thesestatesincludefeaturesfrom ∆W =0.9 in the measurements. Parameters as in Fig. 2. both spot and stripe chimeras. Examples are presented below. B. Leaky Integrate-and-Fire: Stripes A. FitzHugh-Nagumo: Grid patterns Stripesarealsopresent intheLIF model forlarge values of the coupling strength and the refractory period, Thegridchimeraisformedbyaregulargridofk×kin- and medium values of the coupling range. A stable con- coherent spots. The spots are either round or spiral-like figuration containing 6 coherent regions separated by 6 and have approximately the same size given the choice incoherent ones is presented in Fig. 11. In fact, the co- of r and ϕ. When approaching the direct-coupling case herentregionsareseparatedin2groups,onewithhigher with ϕ≈π they shrink to a size of only a few oscillators meanphasevelocitythantheother. Thisisevidentfrom and disappear. A complex coherent grid pattern with- the right panel of Fig. 11, where the section of the mean out incoherent cores remains, which consists of multiple phase velocity at j =35 demonstrates high values in the clusterswiththesamemeanphasevelocity,butnon-zero order of 2.15 for the first coherent group, intermediate phase lags. As seen in the enlargement of the param- valuesintheorderof2.10forthesecondcoherentgroup, eter space for the grid chimera states in Appendix A, and ω -values between 2.06 and 2.12 for the incoherent these cluster states partially surround the corresponding ij stripes. Other stripe multiplicities were also observed grid chimeras. Similar cluster states are also found for for different values of the coupling range R, in the in- the double spot lines and the twisted chimeras discussed terval 15 ≤ R ≤ 25. For smaller (larger) values of R below and also in Appendix B. a larger (smaller) number of stripes is supported with Figure 12(a) shows an example of spirals arranged on smaller (larger) width (not shown). Depending on the a grid with asynchronous, incoherent cores. The grids initialconditionsthestripesformparallelorperpendicu- canalsoappearasrotatedversionswithanon-zeroangle lar to the i-axis. In the LIF model diagonal stripes have comparedtotheunderlyinglatticestructure(notshown). 8 In this work, only the regular grids with the same orien- tation as the oscillator grid are presented. FIG. 12: (Color online) FHN system: (a) Grid state with 6×6 spots for r = 30 and ϕ = π/2+0.82. The k ×k- structure is visible in the (color coded) activator variable xij FIG. 13: (Color online) FHN system: Grid states: Depen- (left) at t=2000 and the (color coded) mean phase velocity dence of the number of asynchronous heads k in a single line (center). The right column shows a section of ωij. (b) Lines on the coupling radius r. The red line is an exponential fit. of incoherent spots for r=22, ϕ=π/2+0.96, and t=1900. The brightness of the dots indicates the number of occur- Other parameters in Fig. 2. rencesofk×k gridchimerasforϕ∈(1.71,4.65). Parameters as in Fig. 2. The number of spots in the baseline is found to be evenineverygridconsideredanddecreasesexponentially withincreasingcouplingradiusr,asshowninFig.13. In of an odd number of spot lines are observed as well, but case of an odd number of spots in a baseline, the spirals they are transients, that is, the lines tend to move, until cannot fulfill the boundary conditions. This leads to a two of them collide to a single line of spots. destructionofthestate. Nevertheless,thegridscanexist We also observe that the coherent wave pattern which as a k×l-state, where both k and l are even numbers. surroundstheincoherentspots[seeFigure12(b)],canbe Another grid pattern, which is observed commonly in orientedinspaceinotherwaysthatincludeananglewith the two-dimensional FHN model for a large parameter theunderlying grid-structureof the n×n oscillatorgrid, region, is a multi-headed, double-spotted line chimera whileleavingthestraightlinesofincoherentspotsintact. state. This state consists of two parallel straight lines We find in a distribution similar to Fig. 13, that for a of incoherent spots, while the surrounding area is filled given number of spots k the solutions that are aligned with a frequency-locked coherent pattern that exhibits withthegrid-structurehavebigcouplingradii,whilethe a gradual phase shift in space. An example for such a radius decreases when the angle increases (not shown). spotted chimera pattern is shown in Fig. 12(b). The spotted patterns can be arranged in a horizontal/vertical or diagonal configuration. When the cou- B. Leaky Integrate-and-Fire: Grid patterns pling radius increases, the number of spots in a single line decreases exponentially for both configurations (not In the LIF network, grid chimera patterns emerge shown), similar to the grid chimeras in Fig. 13. Again, when the coupling range takes relatively high values, they can be spiral shaped or round and shrink or disap- see Figs. 14 and 15. Before studying the effects of the pearcompletelynearthedirect-couplingcasewithϕ≈π, coupling range R on the form of chimera states, we whichleavesacoherentpattern. Forallcouplingphases, present the case of multistability observed for high val- every second spot in a line can be shifted off the line for ues of σ = 0.7 and intermediate values of R = 22. Fig- big coupling radii. ure 14 depicts snapshots of u (left panels), correspond- ij All spotted line states are multistable to other pat- ing mean phase velocities (middle panels) and section of terns of similar geometry and the grid states. Starting mean phase velocities (right panels) for LIF realizations from the basic patterns of two parallel lines of spots, the starting from two different random initial conditions. In change of parameters induces a transition to a different panels (a) we observe a 36-headed chimera arranged in a chimera pattern, wave patterns, or to synchronization. 6×6 grid, while panels (b) show a 12-headed chimera in One possibility is the collapse of the two spotted lines theconfigurationoftwohorizontaldouble-spottedbands into a single line, which then forms an incoherent stripe. with6incoherentregionseach. Wemuststressthatpat- Another possibility is the formation of more than two tern (b) is very rare, but stable in time and is present lines of incoherent spots, mostly 4 lines. Configurations only in a few cases around 22 ≤ R ≤ 24 and for spe- 9 cific initial conditions. Most initial conditions support A grouping phenomenon of mean phase velocities is the 36-headed chimeras. clearly visible in this 16-headed grid chimera pattern of Fig. 15(b). The coherent regions split into four groups, with each one of them having a distinct mean phase ve- (a) uij ωij locity, whiletheωij-levelsofthefourcoherentgroupsal- 100 1 100 2.83 2.83 ternateinspace. IntherightpanelofFig.15(b)onlytwo 80 0.8 80 ofthesegroupsarevisiblesincethesectionforj =25cuts j 4600 00..46 j 4600 2.8 2.8 j=15ωij throughthe2.785and2.79regions. Tracesofthisgroup- 20 0.2 20 ingphenomenonarealsovisibleinFigs.14(a)and15(a). 0 0 20 40 60 80 100 0 0 0 20 40 60 80 100 2.77 0 20 40 60 80 100 2.77 As the coupling range grows, the difference between the i i i mean phase velocities of the four groups increases (not (b) shown) and for large values of R the four groups are dis- 100 1 100 2.83 2.83 tinguishable. 80 0.8 80 Another observation worth studying is that for large j 4600 00..46 j 4600 2.8 2.8 j=25ωij values of σ the incoherent regions exhibit small ωij- values, whilethecoherentonescorrespondtolargeω ’s. 20 0.2 20 ij 0 0 0 2.77 2.77 Comparing the mean phase velocities between Figs. 7 0 20 40 60 80 100 0 20 40 60 80 100 0 20 40 60 80 100 i i i and 15, we observe that their values have been doubled. We find that it is the coupling strength that controls the FIG. 14: (Color online) LIF system: (a) Grid chimera, (b) magnitude of them. As mentioned in Sec. IIIB for small Double-spotted line. (a) and (b) patterns are formed starting from two different random initial conditions. Left pan- values of σ the ωij-values of the incoherent regions are els: snapshots of u ; middle panels: corresponding mean largerthanthecoherentones,whileweobserveherethat ij phase velocities; right panels: section of mean phase veloc- the opposite is true for large values of σ. This qualita- ities. Parameters: R = 22, σ = 0.7, µ = 1.0, pr = 0.22Ts, tive change indicates that there is an intermediate value and N =100. of the coupling constant where this qualitative change takes place. For the LIF model this critical point is lo- In Fig. 15 we present the snapshots of u (left pan- cated between the values 0.2 < σ < 0.6. It should ij crit els),correspondingmeanphasevelocities(middlepanels) correspond to full synchronization or traveling waves at and sections of mean phase velocities (right panels) for the cross-over ω = ω . Our simulations show that incoh coh LIFrealizationsfortwodifferentcouplingrangesR. The σ ∼ 0.3. For this value, we observe full synchroniza- crit 36-headed chimera persists as we increase the coupling tion, while for σ between approximately 0.4 and 0.5 the range, up to R ∼ 30. Above this value the system syn- systemstaysinthetransientregimeforalongtime. This chronizes (not shown). As the limit of all-to-all synchro- difficultytoattainastationarystateisanotherindication nizationisapproachedthereisanotherwindowofpartial that we are near a critical point. synchronization around R = 45, where a 16-headed grid Overall, the number of incoherent elements decreases chimera is observed, see Fig. 15(b). with the coupling strength σ as shown in Fig. 16. This figure is based only on the 36-headed chimera patterns, since the 16-headed is hard to detect in higher σ values. (a) u ω The exponential fit (solid line) describes closely the de- ij ij 100 1 100 2.84 2.84 crease of the size of the incoherent regions with R, with 80 0.8 80 an exponent β =−3.25 characterizing the decay. j 4600 00..46 j 4600 22..882 22..882j=15ωij staTtehseicnonthceluLsiIoFnmreoladteeldistothtahtetahtetrribatuitoesofofcothheercehnitmaenrda 20 0.2 20 0 0 0 2.78 2.78 incoherent elements and the mean phase velocities are 0 20 40 60 80 100 0 20 40 60 80 100 0 20 40 60 80 100 i i i mainly controlled by the coupling strength σ. (b) 100 1 100 80 0.8 80 2.79 2.79 C. Comparison j 4600 00..46 j 4600 2.78 2.78j=25ωij 20 0.2 20 2.77 2.77 Both models support the presence of incoherent spot 0 0 20 40 60 80 100 0 0 0 20 40 60 80 100 2.76 0 20 40 60 80 100 2.76 patterns arranged in grids. The spots are organized i i i equidistantly along straight lines. The grid consists of FIG. 15: (Color online) LIF system: The snapshots of u an even number of such spotted lines in the two spatial ij (left panels), corresponding mean phase velocities (middle directions, along which the number of lines is not neces- panels)andsectionsofmeanphasevelocities(rightpanels)in sarily equal. In the LIF model an odd number of lines is thelimitofstrongcouplingfortwodifferentcouplingranges: not observed for the considered parameter scans. In the (a) R = 23 and (b) R = 45. Parameters: σ = 0.7, µ = 1.0, FHN model spots arranged in odd numbers of lines may pr =0.22Ts, and N =100. arise, but they soon merge, zipping together and giving 10 models are described by different equations, both systems support hybrid states composed of coherent and incoherent regions when the elements are coupled. We have identified a number of common chimera patterns: spots, grids, rings, and stripes. Our simulations suggest that the coherent/incoherent pattern characteristics follow similar growth rules. For example, the diameter of the ring patterns grows linearly with the coupling range in both models. Other phenomena typical in nonlinear systems such as multistability and transitions between different patterns have been observed as well. The common behavior of the two models supports the universal occurrenceofthesepeculiardynamics: Chimerapatterns are persistent and independent of the specificity of the model, provided that the models retain the characteristics of spiky limit cycles. FIG.16: (Coloronline)LIFsystem: Theratioofincoherent Future investigations are needed to show whether elements in all observed 6×6-grids as a function of the cou- these chimera patterns will be persistent for a three- pling strength σ, using a fixed coupling range at R=25 and dimensional configuration and in more complex network other parameters as in Fig. 7. architectures, drawing from the current advances in the realistic recordings of the neuron connectivity in the different parts of the brain. Another aspect worth inves- risetostablepatternswithevennumberofspotsinboth tigating is the inhomogeneity of the single neurons. Bi- directions. ological experiments have shown that microscale inhomogeneities including different neuron bodies, dendritic structures, axonal fiber bundles etc. are common in VI. CONCLUSIONS AND OPEN PROBLEMS the brain structure and that they affect electrical sig- nalpropagationonamicroscale[43]. Ontheotherhand, Wehavestudiedthedynamicsoftwo-dimensionalnet- numerical experiments on the one-dimensional ring ar- works of neuronal oscillators with nonlocal coupling. chitecture support the view that at low levels of inho- This can be seen as an intermediate step extending mogeneity of the single neuron parameters the chimera the intensively studied one-dimensional ring geometries patterns persist [26]. It will be insightful to investigate towards a three-dimensional arrangement of the brain. thisphenomenoninthecaseofhigherspatialdimensions FitzHugh-Nagumo (FHN) and Leaky-Integrate-and-Fire with complex connectivity patterns inspired by the ones (LIF) models have been chosen to represent the single recorded in brain medical experiments. neuronal activity, as they are two of the most prominent paradigmatic descriptions. Finding common synchronization patterns in the two dynamical networks could VII. ACKNOWLEDGMENTS point the way toward identification of universal dynamical features present in brain activity. The aim of the current study was to concentrate on AS thanks Yu. Maistrenko for fruitful discussions. stable chimera patterns induced by the nonlocal connec- Funding for this study was provided by NINDS tivity and to identify patterns that are common in both R01-40596. PH acknowledges support by Deutsche models. In each model, we have considered a parameter Forschungsgemeinschaft (DFG) in the framework of the space spanned by three quantities. For the FHN model Collaborative Research Center 910. This work was par- thephaseconnectivityparameterandthecouplingrange tially supported by the SIEMENS research program on have been varied, while the coupling strength was kept “Establishing a Multidisciplinary and Effective Innova- fixed. FortheLIFmodelthecouplingstrength,coupling tion and Entrepreneurship Hub” and by computational range, and refractory period were varied. timegrantedfromtheGreekResearch&TechnologyNet- Our comparative study has demonstrated that, al- work (GRNET) in the National HPC facility - ARIS - though the dynamics of the single neurons in the two under project ID PA002002. [1] D.M.AbramsandS.H.Strogatz,Phys.Rev.Lett. 93, [3] M. J. Panaggio and D. Abrams, Nonlinearity 28, R67 174102 (2004). (2015). [2] Y.KuramotoandD.Battogtokh,NonlinearPhenomena [4] M.R.Tinsley,S.Nkomo,andK.Showalter,NaturePhys. in Complex Systems 5, 380 (2002). 8, 662 (2012).