Clustered Chimera States in Systems of Type-I Excitability

Chimera is a fascinating phenomenon of coexisting synchronized and desynchronized behaviour that was discovered in networks of nonlocally coupled identical phase oscillators over ten years ago. Since then, chimeras were found in numerous theoretical and experimental studies and more recently in models of neuronal dynamics as well. In this work, we consider a generic model for a saddle-node bifurcation on a limit cycle representative for neural excitability type I. We obtain chimera states with multiple coherent regions (clustered chimeras/multi-chimeras) depending on the distance from the excitability threshold, the range of nonlocal coupling as well as the coupling strength. A detailed stability diagram for these chimera states as well as other interesting coexisting patterns like traveling waves are presented.


Introduction
About ten years ago, a peculiar synchronization phenomenon was reported: in a system of nonlocally coupled oscillators, a state was discovered where synchronous and asynchronous oscillators coexist, even though the oscillators are identical and the interaction symmetric and translational invariant [1]. This phenomenon was termed the name "chimera" after the monstrous fire-breathing creature of Greek mythology composed of the parts of three different animals, a lion, a snake and a goat [2]. From the perspective of nonlinear dynamics, this surprising break of symmetry is observed by the coexistence of incongruent states of spatial coherence and disorder.
Real-world examples that exhibit a chimera state include electric-power grids, which rely on synchronized generators to avoid blackouts in power transmission. Also, certain patterns of intense heart-tissue contraction known as "spiral waves" in certain types of heart attacks have been observed in simulations of chimera states. Forms of chimera state may also be connected to large-scale synchronization patterns of neurons that have been observed during seizures. For a comprehensive review refer to [3] and references therein.
Chimera states were first reported by Kuramoto and Battogtokh in a model of densely and uniformly distributed oscillators, described by the complex Ginzburg-Landau equation in one spatial dimension, with nonlocal coupling of exponential form [1]. This seminal work was followed by the work of Abrams and Strogatz [2], who observed this phenomenon in a 1-dimensional ring continuum of phase oscillators assuming nonlocal coupling with a cosine kernel and coined the word "chimera" for it. The same authors also found chimera states in networks of identical, symmetrically coupled Kuramoto phase oscillators [4] by considering two subnetworks with all-toall coupling both within and between subnetworks, assuming strong coupling within each subnetwork and weaker coupling between them. This coupling scenario was also employed by C. G. Laing who demonstrated the presence of chimeras in coupled Stuart-Landau oscillators [5]. More recently, the same coupling scheme was used in a system of pendulum-like elements represented by phase oscillators with a second derivative term, where chimera states were also investigated [6]. Furthermore, Stuart-Landau oscillators have also been investigated related to amplitude-mediated chimera [7] and for symmetry-breaking coupling. The latter leads to a combination of chimeras and oscillation suppression, termed chimera death [8]. Chimeras have also been observed in many other systems, including coupled chaotic logistic maps and Rössler models [9,10]. Together with numerical, the theoretical studies of chimera states have been recently provided, such as general bifurcation analysis for chimeras with one and multiple incoherent domains in the system of nonlocally coupled phase oscillators [11].
The first experimental evidence of chimera states was found in populations of coupled chemical oscillators as well as in optical coupled-map lattices realized by liquid-crystal light modulators [12,13]. Recently, Martens and coauthors [14] showed that chimeras emerge naturally from a competition between two antagonistic synchronization patterns in a mechanical experiment involving two subpopulations of identical metronomes coupled in a hierarchical network. Furthermore, chimeras were experimentally realized using electrochemical oscillators [15] as well as electronic nonlinear delay oscillator [16].
The importance of chimera states is also very relevant for brain dynamics, since it is believed that they could potentially explain the so-called "bumps" of neuronal activity (proposed in mechanisms of visual orientation tuning, the rat head direction system, and working memory [17]) as well as the phenomenon of unihemispheric sleep [18] observed in dolphins and other animals which sleep with one eye open, suggesting that one hemisphere of the brain is synchronous the other being asynchronous. For this reason, it is particularly interesting that such states were recently observed in leaky integrate-and-fire neurons with excitatory coupling [19], as well as in networks of FitzHugh-Nagumo [20] and Hindmarsh-Rose [21] oscillators.
Excitability is an important feature of neuronal dynamics [22] as it determines the mechanism of the generation of action potentials (spikes) through which neurons communicate. There are two types of excitability: type I yields a response of finite amplitude and infinite period through a global bifurcation, and type II gives rise to zero-amplitude and finite period spikes via a Hopf bifurcation. Type-II excitability is often modeled by the FitzHugh-Nagumo system for which "multi-chimera" (or "clustered chimera" [23]) states, which consist of multiple coherent regions, were recently found slightly above the excitability threshold [20]. The Hindmarsh-Rose model which is representative for both type-I and type-II excitability, exhibits very complex behaviour including spiking, regular and chaotic bursting for which chimera states and other collective dynamics were identified [21].
In this work, we will focus on a generic model for type-I excitability and we will focus on the fundamental dynamics by performing a systematic analysis as far as chimera states are concerned. The system under consideration is representative for a global bifurcation, namely a saddle-node bifurcation on a limit cycle also known as Saddle-Node Infinite PERiod (SNIPER) bifurcation, which is also known as Saddle-Node bifurcation on an Invariant Circle (SNIC). It is defined by the following equations [24,25,26,27]: with the state variables x(t) and y(t), and b is the bifurcation parameter. For b < 1, there are three fixed points: an unstable focus at the origin and a pair of a saddle-point and a stable node on the unit circle with coordinates (b, , respectively. The latter two collide for b c = 1 at (x * , y * ) = (1, 0) and a limit cycle with constant radius ρ c = x 2 − y 2 = 1 is born. Above but close to the bifurcation, the frequency f of this limit cycle obeys a characteristic square-root scaling law f ∼ √ b 2 − 1. In the following, we choose b > b c so that the system operates in the oscillatory regime. The system oscillates with constant amplitude ρ = 1 and the period T 0 is given by 2π/ √ b 2 − 1. In figure 1 the numerical solution of x and y is shown for one period. For b = 1.05 (figure 1(a)), the dense region (the so-called "ghost") where the system slows down marks the collision point of the saddle and the node, i.e. (x * , y * ) = (1, 0). For this parameter value, the system remembers the collision point because it is close to the critical value b c . The phase velocity converges to a constant value as soon as b becomes large enough ( figure 1(b)). The rest of this paper is organized as follows: In Sec. 2, we introduce the coupling topology and describe the main features of the observed dynamics. In Sec. 3, we scan the parameter plane spanned by the bifurcation parameter and coupling range. Section 4 focuses on coexistence of chimeras and other patterns and in Sec. 5, we address the role of the coupling strength. Finally, we conclude with a summary in Sec. 6.

The model
We consider N nonlocally coupled SNIPER oscillators given by Eqs. (1) arranged on a ring:ẋ where k = 1, 2, . . . , N , σ > 0 is the coupling strength, and R ∈ [1, N/2] is the number of nearest neighbors of each oscillator on either side. The limit cases R = 1 and R = N/2 correspond to nearest-neighbour and all-to-all coupling, respectively. It is convenient to scale this parameter by the system size, which defines a coupling radius r = R/N ∈ [1/N, 0.5]. The coefficients b lm , where l, m ∈ {x, y}, are given by the elements of the rotational matrix: where φ ∈ [−π, π]. The matrix B allows for direct (xx)-and (yy)-coupling as well as cross coupling between x and y as in [20]. The diffusive coupling in Eqs. (2) is motivated by the electrical synapses (gap junctions) linking real neurons. Neuronal networks have a considerably higher amount of strong short-range connections rather than long-ranged links [28,29,30,31]. This property is implemented in our model by means of R-nearest-neighbour coupling in both directions. Recently, chimera states have also been reported for global coupling involving a mean-field via a nonlinear or linear coupling function as well as time delays [23,7,32,33]. The coupling phase φ parameterizing the matrix B can be related to the so-called phase lag parameter, which is as essential for the existence of chimera states as is the nonlocal coupling [1,2,34]. Figure 2(a) shows a snapshot of the variables x k at a fixed time, providing evidence of a classical chimera state: One group of neighboring oscillators on the ring is spatially coherent (blue dots) while the remaining elements form a a second, spatially incoherent group (black dots). These two domains of coherent and incoherent oscillators can be distinguished from each other through the mean phase velocity of each oscillator ω k = 2πM k /∆T , where M k is the number of periods of the the kth oscillator during a sufficiently long time interval ∆T [20]. Figure 2(b) shows the characteristic profile for the mean phase velocities ω k corresponding to the chimera state of figure 2(a). The oscillators in the coherent domain (blue) rotate along the unit circle at a constant speed ω coh , whereas the incoherent oscillators (black) have different mean phase velocities ω incoh with a maximum value denoted by ω ext incoh . If the difference defined as is sufficiently larger that a certain threshold value, we can ensure the existence of a chimera state. Note that, for the particular chimera state of figure 2(a), it holds that ω ext incoh > ω coh . Figure 2(c) shows the corresponding space-time plot for the variables x k . For weak coupling, which is the case here, the period of the oscillators converges to the period T 0 of the uncoupled system. Investigations of space-time plots for extended simulation times reveal that the (in)coherent domains are stationary, i.e. there is no "drift" on the ring. Finally, figure 2(d) shows the state of each oscillator at a certain time t in phase space (the blue dots mark the coherent oscillators while the black dots the incoherent ones).
In the following sections, we will systematically investigate the effect of the bifurcation parameter b as well as the coupling parameters R and σ on the chimera state. We will compare our results with findings of previous studies on chimera states in neuronal networks and shed light on new dynamical features. For the numerical integration of the Eqs. (2) we used the Euler method with step size dt = 0.01. The initial conditions for x k and y k are randomly distributed on the unit circle and we discard transients of 1000 time steps. For the mean phase velocities ω k , we average over a time interval ∆T = 10.000.
3. Impact of the bifurcation parameter and coupling range on the dynamics.
A stability diagram for the chimera states is displayed in figure 3 where the dependence of the modulus of ∆ω (equation (4)) is plotted with respect to the bifurcation parameter b and the coupling radius r = R/N .
Starting from the values b = 9, r = 0.43 and a certain set of initial conditions as described above, we perform a continuation on the direction of smaller r-values down to r = 0.06 and calculate ∆ω for each coupling radius. Subsequently, for values of r ∈ [0.04, 0.46] we perform a continuation in b-direction from b = 9 down to b = 0.1 starting again at r = 0.43. The coupling strength is fixed at a constant value σ = 0.1.
From figure 3 it is clear that |∆ω| has a non-monotonous behaviour in the (b, r)plane. Each "bump" in the 3D surface corresponds to a different type of chimera state associated to a different number of (in)coherent domains, marked in the square/curly brackets. Some of these states are explicitly shown below in figure 4 for certain combinations of b and r. Stability diagram in the (b, r)-plane: Modulus of the difference |∆ω| between the mean phase velocities of the coherent and incoherent oscillators (equation (4)) as a function of the bifurcation parameter b and the coupling radius r. The numbers in the brackets denote the number of the (in)coherent domains of the corresponding chimera state. Square and curly brackets refer to "normal" and "flipped" ω-profile, respectively. Parameters: σ = 0.1, φ = π/2 − 0.1, and N = 1000.
For large values of the bifurcation parameter (red-colored "bumps" in figure 3 and figure 4(a')) a classical chimera state with one group of (in)coherent oscillators exists. By decreasing r, which physically means removing more and more long-range connections, the number of clustered (in)coherent oscillators increases. In the red-colored "bumps" of figure 3 these so-called "multi-chimera" states exhibit the characteristic feature that ω ext incoh > ω coh (i.e. ∆ω > 0), shown in the corresponding mean phase velocity profiles in figure 4(b')-(d'). We denote these chimera states, for which ∆ω > 0, by the number of their (in)coherent domains in square brackets [1], [2], [4], and [6].
For lower values of b (blue-colored "bumps" in figure 3 and figures 4(a-e)), we exclusively find multi-chimera states. As in the case of larger b, the number of clustered (in)coherent oscillators increases with decreasing coupling radius r. However, there is a significant difference: The mean phase velocities of the incoherent oscillators is smaller than the velocity of the coherent ones, i.e. ∆ω < 0. Hence, there exists a critical value of the bifurcation parameter (found to be around b = 4), where ∆ω changes its sign, resulting in a "flip" in the mean phase velocity profile. The chimera states with a "flipped" ω-profile are denoted by the number of (in)coherent domains in curly brackets The characteristic form of the average phase velocities profile is commonly considered as a criteria to distinguish chimera states in the systems of coupled oscillators. The most often observed in the variety of systems is the case when the coherent oscillators perform smaller average phase frequencies, and incoherent oscillators are faster. However, the opposite situation is also possible, when the coherent oscillators perform faster oscillations as the incoherent ones. In the system of nonlocally delay coupled phase oscillators, two types of chimera states were distinguished depending on whether the effective frequencies of the incoherent oscillators are larger or smaller than the frequencies of the coherent ones [35,36]. The regions of stability for these two types of chimera states depend on the time delay and strength of the coupling. Moreover, both types of chimera states can coexist.
The "flipped" phase velocities profile was also observed in systems which do not consider time delay in the coupling, but has not been explained so far. The Kuramoto model with repulsive coupling allows for multi-chimera states for which the phase velocity profiles show larger average frequencies for oscillators that belong to coherent domain [37]. Similar behaviour is also observed for chimera states with one incoherent domain in the complex Ginzburg-Landau equation with nonlocal coupling [7]. In that system, however, chimera states with multiple incoherent domains possess the usually observed mean phase velocity profiles.
Hence, the flip of the phase velocities can not be explained only by the influence of time delay, or strong coupling. This feature was observed in experiments as well, in networks of electrochemical oscillators with nonlocal coupling, the frequencies of the oscillators from the coherent domain of chimera state are larger than the frequencies of the incoherent ones [15].
In our system, we observe direct dependence of the form of the mean phase velocity profile on the parameter b defining the frequency of the local uncoupled unit.

Multistability of patterns: Coexisting chimeras and traveling waves.
The coexistence of different multi-chimeras, traveling waves, and completely synchronized states in the phase space has been observed in many other systems of nonlocal coupled oscillators [20,21,38]. Depending on the initial conditions the stationary state can vary significantly. Such multistable solutions are also possible in system (2) as demonstrated in figure 5. A schematic representation of the identified multi-chimeras in the (b, r)-plane is shown in figure 5(a) . Each region has a different colour associated to a different chimera type as described in the previous section. The black regions correspond to intermittent states, which are mainly desynchronized. Along the white lines, figures 5(b) and (c) display the results of a continuation in b (dashed line) and r (solid line), respectively. The continuation is performed as down sweep in b (or r) and then repeated in the opposite direction. In both cases we find a region where different types of chimera states coexist.
In particular, for intermediate values of the bifurcation parameter b, there is coexistence of a [2]-and {6}-chimera state marked by the shaded (light blue) area of figure 5(b). This area of coexisting chimera states, moreover, marks the transition between "flipped" (∆ω < 0) and "normal" (∆ω > 0) mean phase velocity profile. This transition occurs at a different and, in particular, lower value of b when the continuation is performed in the direction of decreasing b (black dots) than when performed in the opposite direction (red squares), i.e. our system exhibits, apart from multistability, hysteresis phenomena as well.
Coexisting chimera states may also be found by varying parameter r, as shown in figure 5(c): Depending on the choice of initial conditions, one may observe either a [1]-or a [2]-chimera state (shaded, light-blue area) both with ∆ω > 0. In both increasing (red dots) and decreasing r (black dots) directions, there are deviations from the piecewise linear behaviour of ∆ω(r) which correspond to desynchronized states.
The observed multi-chimera states may also coexist with completely synchronized states and traveling waves. One example of such a point in parameter space is marked by the white star in figure 5(a) and the corresponding space-time pattern is shown in figure 5(d). This is a traveling wave solution of wave number 2 coexisting with a {4}chimera state. The time in the vertical axis is scaled by the period T 0 of the uncoupled oscillator. Multistability between traveling waves and breathing states have recently alse been reported for chaotic systems with nonlocal coupling [38].

Role of the coupling strength.
In order to complete our study on the effect of the system parameters on the dynamics of chimera states, we will investigate the role of the coupling strength σ in this section.
Again, we perform a parameter continuation and focus on the behaviour of ∆ω as σ increases for different multi-chimera states. Our findings show that even at large σ the corresponding multi-chimera state is preserved. However, we observe that, for certain values of the bifurcation parameter b and the coupling radius r, the coupling strength may induce a spatial motion of domains of the (in)coherent oscillators. Figure 6 shows the results for the [2]-chimera state associated with the orange regime of figure 5(a). With increasing coupling strength, each oscillator becomes more and more influenced by the dynamics of the remaining oscillators. Therefore, the trajectories of the incoherent oscillators, in particular, begin to deviate significantly from the unit circle as shown in the right plot of figure 6(a) for σ = 1.3.
The corresponding mean phase velocity profiles ω k of the [2]-chimera state can be seen in figure 6(b). For larger σ (right plot) the total number of coherent oscillators increases while the number of incoherent oscillators decreases. Figure 6(c) shows that the difference between the mean phase velocity of the coherent and incoherent oscillators ∆ω linearly increases with the coupling strength, apart for a narrow range of σ ≈ 1 where ∆ω deviates. In this regime, the corresponding space-time plots of the [2]-chimera state reveal that the (in)coherent domains start to move spatially with time (see left inset of figure 6(c)). Beyond this regime of moving patterns, the [2]-chimera state is stationary (see right inset of figure 6(c)). Comparing the two insets of figure 6(c), we observe that for increasing coupling strength the period of the coherent oscillators increases. On the other hand, the period strongly decelerates once the [2]-chimera state introduces spatial motion.
In general, chimera states can be stationary or can perform two types of motion in space, in which the coherent and incoherent domains change their spatial position in time. The first one is a chaotic motion of the position of the chimera observed in nonlocally coupled phase oscillators. Such a motion shows a sensitive dependence on the initial conditions and is a finite-size effect that vanishes in the thermodynamic limit. It can be described as a Brownian motion and depends on the coupling radius, the phase lag parameter, and the shape of the coupling function [34]. The second type is a periodic motion of the coherent and incoherent domains of the chimera state, called "breathing chimera". Breathing chimeras were first observed in the system of two oscillator populations in which each oscillator is coupled equally to all the others in its group, and less strongly to those in the other group [4], and recently in the nonlocal complex Ginzburg-Landau equation with strong coupling limit [7].
The numerical evidence shows, that spatial motion of coherent and incoherent domains in our system is periodic, thus we conclude that for distinct values of parameter b we observe the phenomenon of breathing chimera in our system.

Conclusions
In this work, we have verified the occurrence of clustered chimera states in a generic model for a saddle-node bifurcation on a limit cycle representative for neural excitability type-I. This, along with recent reports on multi-chimera states in nonlocally coupled FitzHugh-Nagumo [20] and Hindmarsh-Rose [21] oscillators provide strong evidence that this kind of symmetry breaking is very relevant for applications in neuroscience.
In particular, we presented a detailed exploration of the parameter space, where chimera states occur, and investigate the dependence on the proximity to the excitability threshold and the range of the nonlocal coupling. We identified chimera states for which the mean phase velocity has a "flipped" profile. A similar result was also reported in a recent study of Kuramoto oscillators with repulsive coupling [37]. Findings of coexisting chimera states and traveling waves in the parameter space establish the existence of multistability in our model. Finally, it was shown that for increasing coupling strength the domains of coherent oscillators become bigger and at the same time spatial motion of the incoherent oscillators is observed.