Emergence of chimera states in a neuronal model of delayed oscillators

Neurons are traditionally grouped in two excitability classes, which correspond to two different responses to external inputs, called phase response curves (PRCs). In this paper we have considered a network of two neural populations with delayed couplings, bound in a negative feedback loop by a positive PRC (type I). Making use of both analytical and numerical techniques, we derived the boundaries of stable incoherence in the continuum limit, studying their dependance on the time delay and the strengths of both interpopulation and intrapopulation couplings. This led us to discover, in a system with stronger delayed external compared to internal couplings, the coexistence of areas of coherence and incoherence, called chimera states, that were robust to noise. On the other hand, in the absence of time delays and with negligible internal couplings, the system portrays a family of neutrally stable periodic orbits, known as “breathing chimeras.”


I. INTRODUCTION
For the past few decades, models of coupled phase oscillators have proved to be particularly successful to describe the emergence of macroscopic rhythmic patterns in a huge variety of natural and artificial contexts [1].In this framework, it is useful to consider the interplay between excitatory and inhibitory (E-I) time-delayed connections, in order to model spatially distributed self-organized systems, such as neuronal networks, that are known to exhibit synchronous behavior.This is associated with many cognitive processes such as memory formation [2], directed attention [3], and the processing of sensory stimuli [4], but it can also be the hallmark of certain disease states such as Parkinson's disease or epileptic seizures [5][6][7].In these cases, synchronization is generally localized to certain cerebral regions only, and coherence and incoherence coexist within the brain.
Synchronization patterns largely depend on the individual oscillators properties.It is known that neurons can be grouped into two excitability classes, that differ in the bifurcation observed while transitioning from the silent mode to the firing mode, in particular, a saddle node on an invariant circle bifurcation or a Hopf bifurcation.Moreover, they are directly linked to two types of phase response curves (PRCs): either always positive (type I) or both positive and negative (type II) [24].The dynamics of networks of type II neurons has been widely explored in the past, e.g., making use of the Kuramoto model and its many generalizations [25][26][27][28][29][30].In this context, chimera states appear commonly in networks of two subpopulations with nonlocal couplings [31], typically with a large ensemble of oscillators (although they have also been observed with as few as two oscillators per group [32]).It has also been shown, in an adaptive Kuramoto model, that asymmetric inter-or intrapopulation couplings enhance the transition from the chimera state to the synchronized state [33].Moreover, the presence of delayed couplings, which enriches the dynamics of the system by making it infinite dimensional, has been also observed to lead to the emergence of chimera states, for instance in Kuramoto-like oscillators [34,35] or in multilayer networks of Fits-Hugh Nagumo oscillators [36].
On the other hand, there have been fewer investigations on the properties of type I neurons, which are generally considered to have low propensity for synchronization when coupled by excitation [37,38].Most of the studies [39] only focused on either purely excitatory or inhibitory couplings, whereas Keane et al. [40] studied a small-world network with a ring of excitatory type I oscillators, with some random inhibitory delayed long interactions.It was shown that, as opposed to the same network of type II neurons, where adding the inhibitory couplings results in a sharp transition to incoherence [41], in this case there were multiple transitions between synchronization and desynchronization.The interplay between E and I was also studied by Ladenbauer et al. [42] in a two-population model of delayed integrate-and-fire neurons, where they observed how adaptation controls synchronization and the formation of cluster states.Moreover, Montbrió and Pazó [43] recently analyzed a network of instantaneously coupled E-I oscillators, showing that synchronization can emerge only if the E population is intrinsically faster than the I one.Finally, clustered chimera states have been observed by Vüllings et al. [44] in a ring of type I neurons, which depended on the distance from the excitability threshold, the range of nonlocal couplings, and the coupling strength.
We build on the basis of these results and apply this to a network of two populations of E-I identical type I oscillators with time delayed couplings [Figs.1(a) and 1(b)].The purpose of this paper is to analyze the dynamics of the system, with a particular focus on determining whether chimera states may arise in this framework and how they depend on the delay and the coupling strengths, both within and between the two populations.
The paper is organized as follows: the stability of the fully incoherent state is analyzed in Sec.II A with identical oscillators, no noise, and time delayed interpopulation connections; in Sec.II B with the addition of delayed or instantaneous internal couplings; and in Sec.II C with heterogeneous natural frequencies, noise, and instantaneous interpopulation couplings.Since it was observed that, starting from random initial conditions, the system approaches a chimera state for periodic values of the external time delay greater than zero, in Sec.II D the basin of attraction of the chimera is discussed making use of the Ott-Antonsen (OA) ansatz [45], in the absence of time delay.Finally, in Sec.II E the optimal parameters to display chimera states and the robustness of the system towards noise are evaluated.

II. METHODS AND RESULTS
The governing equations for the phases θ of the ith oscillator in the E and the I populations, with positive PRC, are [24,37,38 where ω where σ = {E , I} and ω σ i = ω σ i + 1 2 (k σ σ + k σ σ ), showing that the excitatory and the inhibitory couplings have the effect to slightly shift the natural frequency of the oscillators.
The order parameter that quantifies the level of synchronization within each population is defined as For later use, taking the long time average • t of the absolute value of the order parameter, it is possible to obtain a single measure for the phase ordering, Z 2 σ , given by Finally, the average order parameter between the two populations is defined as 2) can be expressed in terms of the order parameter, such that where * indicates the complex conjugate.
The analysis of this system can be efficiently carried on in the case N σ → ∞, so that the discrete sets of phases and frequencies turns into the continuum limit {θ σ i , ω σ i } → {θ σ , ω σ }.In other words, for each ω there is a continuum of oscillators distributed along the unit circle; the probability density function is characterized by f σ (θ |t, ω) such that f σ (θ |t, ω)dθ gives the fraction of oscillators with frequency ω which lie between θ and θ + dθ at time t [29].
In general, in this limit the order parameter can be expressed as where g σ (ω) is the frequency distribution.Hence, f σ (θ, ω, t ) satisfies the Fokker-Planck equation It is convenient, for the following analysis, to introduce the Fourier expansion of f σ : Substituting Eq. ( 6) into Eq.( 5) and setting the correspondent modes equal yields To simplify the future notation, we define A linear stability analysis was performed around the fully incoherent state (Z σ = 0), where all the oscillators are randomly distributed along the unit circle ( f σ = 1/2π ).This is always a trivial solution.To do so, a small perturbation was added to the probability density function: with η 1.The analysis, reported in the Supplemental Material [46], reveals that the only possible unstable modes are l = ±1 and considering the case l = 1 results in where • t−τ σ σ means that the quantity within the brackets is evaluated at time t − τ σ σ .

A. Identical oscillators with no internal couplings and positive external delay
First we considered the simplest case with identical oscillators (ω E = ω I ≡ ω 0 ), no noise (D = 0), no internal couplings (k EE = k II = 0), and the same external coupling strength (|k EI | = |k IE | ≡ K).After substituting the ansatz δ f σ 1 (ω, t ) = b σ (ω)e λt , and considering that the frequency distribution is a delta function for both populations (g(ω) = δ(ω − ω 0 )), a brief analysis, reported in the Supplemental Material [46], leads to the characteristic equation for λ: Since the dependence is only on the sum of the two delays, we define τ EI + τ IE ≡ τ for simplicity of notation.
At the bifurcation point between stable and unstable incoherence λ crosses the imaginary axis, λ = iR with R ∈ R, therefore Given that the right-hand side is real, the left-hand side must also be real, so sin(Rτ ) = 0, which leads to the condition Finally, substituting Eq. ( 12) into Eq.( 11) and considering that the natural period is T = 2π/ω 0 , Eq. ( 11) can be solved for K/ω 0 , which yields Therefore, the boundaries of stable incoherence (BSI) are periodic in the delay and only depend on the sum of the delays in the two pathways.This aspect was fully confirmed by the numerical analysis, with the correspondent parameters used listed in the Supplemental Material [46].The simulations were performed in MATLAB2020B, making use of the built-in function dde23 (Fig. 2) and also confirmed via the Euler method (not shown).They were performed first with time delay in the inhibitory connection only [Fig.2(a)], excitatory connection only (shown in the Supplemental Material [46]), and then in both [Fig.2(b)].Moreover, they were repeated over a range of values for K and τ .As initial condition, the phases of both populations followed a normal distribution centered at zero and with standard deviation 2π .Once the system had equilibrated, Z 2 was computed for each population.From the simulations, four possible scenarios emerged, that repeat periodically (denoted I-IV).In particular, the black areas (for instance, area III) were confirmed to be the areas in which both populations are internally not synchronized (Z 2 E ≈ Z 2 I ≈ 0), at the intersection between the curves.In white areas, for instance in area II, both populations were internally synchronized (Z 2 E ≈ Z 2 I ≈ 1).However, the remarkable phenomenon is the emergence of the blue and orange areas such as I and IV, where only one population is coherent while the other is incoherent.These corresponded to chimera states.In Fig. 2(c), a snapshot of the oscillator phases is reported at a fixed time, which confirms the four scenarios depicted above.
Considering only intercouplings between the two populations is a very particular case, that enabled us to highly FIG. 2. Phase diagram obtained by simulating Eq. ( 1) with only external couplings and computing Z 2 E , Z 2 I with (a) τ EI > 0, τ IE = 0 and (b) τ EI = τ IE > 0. The x axis represents the sum of the delays τ in units of the natural period T , while the y axis corresponds to the coupling strength K in terms of the natural frequency ω 0 .The value of the order parameters is indicated by the two-dimensional color gradient on the right.The parameters used are reported in the Supplemental Material [46].The BSI are obtained from Eq. ( 13).(c) Snapshots of the oscillator phases with parameters in regions I-IV which confirm the presence of chimera states in region I and IV.Color scheme as in Fig. 1. simplify the analytical calculations; however, it is not very realistic, therefore we analyzed whether the system phase diagram depicted in Fig. 2 would be robust to the inclusion of small intracouplings.

B. Identical oscillators with both internal and external couplings
First, we considered the case of equal external coupling strengths k EI = k IE ≡ k ext , as well as equal internal ones k EE = k II ≡ k int .Then we analyzed the robustness of the phase diagram depicted in Fig. 2, computing the number of overall synchronous, asynchronous, and chimera states while increasing the values of the internal couplings.As shown in Fig. 3(a), the number of asynchronous and chimera states are robust to the inclusion of internal couplings until k int ≈ 10 −1 k ext , when the system starts to transition to an overall synchronous state.
Secondly, we considered the case of internal and external couplings of equal strength (|k EI,IE,EE,II | ≡ K) while assuming instantaneous intracouplings (τ EE,II = 0) and the same delayed cross-couplings (τ EI,IE > 0), which could represent the case of two neuronal populations that are spatially separated by some distance.In this case, Eq. ( 10) becomes (shown  1) with external delayed and internal instantaneous couplings and computing Z 2 E , Z 2 I .The x axis represents the sum of the delays τ in units of the natural period T , while the y axis corresponds to the coupling strength K in terms of the natural frequency ω 0 .The value of the order parameters is indicated by the two-dimensional color gradient on the right.The parameters used are reported in the Supplemental Material [46].The BSI are obtained from Eq. ( 15).Color scheme as in Fig. 1.
in the Supplemental Material [46]) where again τ = τ EI + τ IE .Following the same reasoning as before, we found the BSI as Therefore, the equations representing the BSI are hyperbolas and vertical lines in the (τ /T, K/ω 0 ) space.This is in accordance with the simulations, as depicted in Fig. 3(b).In this case, compared to Fig. 2, the regions of incoherence shrink around integer values of τ /T , while the area corresponding to the chimera states decreases for increasing values of the external delay.Because the inclusion of small delayed internal couplings did not alter significantly the plot depicted in Fig. 2, while it highly complicated the analytical calculations, we decided, for the rest of the paper, to consider negligible internal couplings.
From the simulations it seemed that the system would never display chimeras without external time delay, hence we analyzed whether the delay was the only necessary condition or whether other parameters could play a role in it as well.In order to do so, we relaxed some of the constraints, in particular that of identical oscillators and identical couplings, and included also some external noise in the system.

C. Heterogeneous oscillators without time delay
We considered the case with no time delay in the connections (τ EI,IE = 0), some external Gaussian noise (D > 0), and heterogeneity in the frequencies of the oscillators.Considering symmetric frequencies over the origin is equivalent to considering a system of reference that is rotating with frequency = (ω E + ω I )/2 and does not modify the analysis of the stability.Moreover, we define |k EI | = |k IE | = k both for simplicity of notation and because we will be interested in observing the ratio of the two coupling strengths, .
The following analysis builds on the work done by Montbrió and Pazó [43].A brief linear stability analysis (see the Supplemental Material [46] for more details), again close to the incoherent state, leads to the condition where ≡ λ + D + γ for simplicity of notation.Solving for λ results in Finally, imposing the condition that Re(λ + ) = 0 gives the BSI: (18) This result confirms, as already shown by Montbrió and Pazó [43], that without time delay, coherence is only possible for positive values of ω, namely, when the E population is faster than the I population.Moreover, increasing the disparity between the E-I couplings causes the area of coherence to shrink.
The simulations were performed with D = 0, spanning the parameter space given by (k/γ , ω/γ ), and they were repeated first with = 1 [Fig.4(a)], namely, with the same coupling strength in the two pathways, and then with = 10 [Fig.4(b)]; again, the precise parameters used are reported in the Supplemental Material [46].In both cases the simulations were in good agreement with the analytical curves, reported as continuous lines.In this case, though, no chimera states were observed.Indeed, the two populations either both synchronized (in the dark green areas where Z 2 ≈ 1) or both did not (in light areas where Z 2 ≈ 0).

D. Existence of chimera states without time delay
In the end, we investigated the possibility that chimera states theoretically exist without time delay in the couplings, even if they have not been observed in the previous simulations.In this context, inspired by Abrams et al. [47], we considered a special class of density functions f σ , that have the form of a Poisson kernel, such that they satisfy the so-FIG.4. Phase diagram obtained by simulating Eq. ( 1) and computing Z 2 with τ EI = τ IE = 0, Lorentzian distributed frequencies D = 0, and (a) same coupling strengths ( = 1) and (b) inhibitory coupling strengths greater than the excitatory one ( = 10).On the x and y axis there are the coupling strength k and the difference in the natural frequencies ω, respectively, both in units of γ , the spread of the frequency distribution.The value of the order parameter is indicated by the color gradient.The parameters used are reported in the Supplemental Material [46].The BSI are obtained from Eq. ( 18).called OA ansatz [45]: These functions are characterized by having the same coefficients in all the Fourier harmonics, except that they are raised to the lth power for the lth harmonic, (19) in Eq. ( 7) with D = 0, and given that Therefore, if this condition is fulfilled, such kernels satisfy the governing equations exactly.The OA ansatz is useful because it reduces the dynamics from infinite to finite dimensional.It is convenient to introduce polar coordinates such as (ρ, φ), defined by Therefore, φ σ can be considered the "center" of the density f σ and ρ σ measures how sharply peaked it is [47].Hence, in terms of its real and imaginary part, Eq. ( 20) can be written as For simplicity, we considered the case with |k EI | = |k IE | ≡ K and we defined ψ ≡ φ E − φ I .Therefore, Eq. ( 21), written explicitly in terms of E and I, reduces to where ω ≡ ω E − ω I .We then considered the case in which the E population is perfectly synchronized (ρ E = 1) and the other is not (ρ I ≡ r), as shown in Fig. 5(a).This further (e) Simulation of the system in the space (ρ E , ρ I , ψ ) following Eq.( 22) represented here in cylindrical coordinates.Color scheme as in Fig. 1.
reduces the dynamics to a two-dimensional system: In this way, the chimera states correspond to the fixed points, namely, r(t ) = const and r = 1 (since r = 1 corresponds to the perfectly synchronized case, that is always a trivial solution) and ψ (t ) = const.Calculating the values of the fixed points yields As expected, considering the opposite case with ρ I = 1 led to the same results [46], given that the system is symmetric.
When ω = 0, that is, with identical oscillators, (r * , ψ * ) = (1/3, 2mπ ), independently of K, which means that chimera states could in principle exist even without time delay, because fixed points different from the perfectly synchronized case exist.To find out the stability of these fixed points, we linearized around them, by computing the Jacobian J. Given the system of differential equations of the form { ṙ = f (r, ψ ) ψ = g(r, ψ ) , the Jacobian J is defined as As the trace is null and the determinant is = K 2 /3 > 0, the linearization predicts a linear center.Moreover, since the system is invariant under the change of variables t → −t, ψ → −ψ, the fixed point is also a nonlinear center [48].Hence, a family of periodic orbits surrounds the chimera, which can be defined as neutrally stable "breathing chimeras" [47].
Both the presence and the nature of the chimeras were fully confirmed by the numerical simulation, which was performed in the case of ω = 0, K = 0.5 [Fig.5(b)].At the same time, the numerical analysis also showed the presence of the perfectly synchronized state, in which both populations synchronize to the same phase, that corresponds to r = 1, ψ = 0, which is a saddle on an invariant circle.
Consequently, we aimed to test whether this descriptionin the continuum limit-would agree with the simulations performed with finite N.For this purpose, it could be exploited that the distribution function f , in the OA ansatz, has the shape of a Poisson kernel, such as which, in this case, becomes Hence, the simulations were performed with the I population initially distributed according to f I for some chosen (ρ I , φ I ), whereas the E population was perfectly synchronized (ρ E = 1), as shown in Fig. 5(a).We obtained, as expected, that while Z E (t ) = const = 1, Z I (t ) oscillates over time, corresponding to the periodic orbit in the (r, ψ ) plane [Fig.5(c)].To explain this concept further, in Fig. 5(d) we report the raster plot of the two populations, with the corresponding order parameter Z 2 (t ) superimposed.Every dot represents the moment each neuron is firing, i.e., when the phase of each neuron completes one full rotation: therefore, the blue dots are perfectly aligned, since the E neurons fire in synchrony, whereas the I neurons oscillate between moments of higher and lower synchrony, corresponding to peaks and troughs of Z 2 I (t ).As expected, the situation was totally symmetric when starting from ρ I = 1 [46].
Moreover, the whole system was simulated in the coordinates (ρ E , ρ I , ψ) following Eq.( 22), shown in the figure as cylindrical coordinates [Fig.5(e)].This revealed that the manifold ρ E = 1, in which the system was restricted in Fig. 5(b), is not attracting.In other words, the system should be "prepared" in this initial condition in order to observe it, which is coherent with our previous simulations.

E. Robustness of chimera states
First, we considered what are the most likely configurations in the parameter space to display the coexistence of coherence and incoherence, by calculating the fraction of chimera states per time delay and per coupling strength observed in Fig. 2(b).This showed that a small (nonzero) time delay, for any value of K, would result in a chimera state; moreover, increasing either τ or K progressively leads to the loss of these states [Fig.6(a)].
Secondly, the robustness of the chimeras towards external perturbations was also evaluated, simulating the system with and without time delay in the cross-couplings with increasing values of D, which quantifies the level of stochastic noise [Fig.6(b)].As expected, it appears that the system without time delay is not very robust, as the manifold that corresponds to the chimera is not attracting, whereas time delays help to stabilize the system until higher values of noise.

III. CONCLUSIONS
In this paper we have sought to shed light on the properties of type I oscillators, that are traditionally believed to have low propensity to synchronize.To do so, we modeled a feedback loop between an excitatory and an inhibitory population, which is a known neural mechanism that produces oscillations.The analysis revealed that including time delayed couplings highly enriches the dynamics of the network.In case of stronger external than internal couplings, we observed the emergence of stable periodic chimera states, that are robust to noise.On the other hand, with negligible internal connections and instantaneous external ones, a family of "breathing chimeras" was observed, that is neutrally stable and less robust to external perturbations.
Future work may investigate whether the chimera states are preserved also in more realistic scenarios, such as when considering local couplings in a spatially distributed system.Moreover, this paper offers a theoretical framework to investigate more biologically inspired neuronal models and testable predictions that can be potentially verified experimentally.

FIG. 1 .
FIG. 1.(a) Schematic representation of the system with two subpopulations [blue, excitatory (E); orange, inhibitory (I)] having N E ,I number of oscillators with intrinsic frequency ω E ,I .They are bound in a negative feedback loop with delayed connections (with time delay τ EI,IE,EE,II , interpopulation coupling strengths k EI,IE , and intrapopulation coupling strengths k EE,II ).(b) Type I PRC, with the form (1 − cos θ )/2.

FIG. 3 .
FIG. 3. (a) Percentage of "synchronized" (squares), "asynchronized" (crosses), and "chimera" states (circles) among all the measured states at increasing values of the internal couplings (k int ) compared to the external ones (k ext ).(b) Phase diagram obtained by simulating Eq. (1) with external delayed and internal instantaneous couplings and computing Z 2 E , Z 2 I .The x axis represents the sum of the delays τ in units of the natural period T , while the y axis corresponds to the coupling strength K in terms of the natural frequency ω 0 .The value of the order parameters is indicated by the two-dimensional color gradient on the right.The parameters used are reported in the Supplemental Material[46].The BSI are obtained from Eq. (15).Color scheme as in Fig.1.

FIG. 5 .
FIG. 5. (a) Distributions of E-I: I follows a Poisson kernel [Eq.(24)] centered at φ I , with the spread dependent on ρ I = r; E is a δ function centered at φ E , with ρ E = 1, as E is perfectly synchronized.The distance between the centers of the two distributions is ψ = φ E − φ I .(b) Trajectories in the (r, ψ ) plane, in polar coordinates, obtained by simulating Eq. (23) with = chimera state and • = perfectly synchronized case.( , , ) correspond to the parameters used in (c).(c) Order parameter plotted as a function of time for the E, I populations with N E ,I = 250.(d) Raster plot of the E-I neurons at the breathing chimera state for t ∈ [100-150 s].(e) Simulation of the system in the space (ρ E , ρ I , ψ ) following Eq.(22) represented here in cylindrical coordinates.Color scheme as in Fig.1.

FIG. 6 .
FIG. 6.(a) Fraction of chimera states registered in Fig. 2(b) calculated per coupling strength, which corresponds to the light-blue curve, and per time delay, shown by the red curve.(b) Noise robustness comparison of the chimera states between delayed (dashed and dotted lines) and instantaneous couplings (solid line), obtained computing the order parameter of the synchronized population (Z 2 E ) for increasing values of D.