The chimera state on a ring was obtained from the integration of N Kuramoto oscillators with nonlocal coupling. The equations are given by Abrams and Strogatz (2006):

with i = 1, ..., N, where ϕ_i is the oscillator's phase and ρ is the oscillators intrinsic frequency. To obtain a stable chimera we integrated Equation 1 using chimera-like initial conditions and we set A = 0.95, β = 0.2 and N = 500 as in Laing (2009). We chose A and β from the (A, β) parameter plane in which the chimera state exists (see ref. Abrams and Strogatz (2006) for details) and set N large enough (N>50) such that for this type of network the chimera did not collapse (note that the number of oscillators can be reduced when using a different topology as the one used in Equation 1, see ref. Panaggio et al. (2016) for examples). We obtained chimera-like initial conditions by randomly selecting the same phase for half of the network. The phases for the other half were selected from a uniform distribution between [0, 2π]. See ref. Masoliver et al. (2022) for more details. The equations were integrated using the Euler method with an integration step of dt = 10⁻³. Note that all equations are dimensionless, however time can be rescaled so that a single unit of time, which corresponds to 8 cycles of the synchronized population, can be rescaled to 1 second which yields an 8 Hz (theta) oscillation.

The two-population chimera consists of two populations of n Kuramoto oscillators each. The phases of the oscillators for group 1 and group 2 are given by γ={γi}i=1n and ϕ={ϕi}i=1n, which are governed by the following equations:

The coupling within groups is given by μ=1+A2n and between groups by ν=1-A2n, with ν < μ, 0 ≤ A ≤ 1 and n = 3. To obtain a stable chimera, we simulated Equations 2 and 3 with appropriate initial conditions as in ref. Masoliver et al. (2022) and we fixed β = 0.025 and A = 0.1. The temporal component τ=10.012 is used to slow down the chimera dynamics, which is needed in order to successfully train the spiking recurrent network. It is for that value that we get 8 and 9 oscillations, for the synchronized and the unsynchronized populations, respectively, in 1 second. Figures 2E, 3A illustrate the network. The equations were integrated using the Euler method with an integration step of dt = 10⁻³. Note that all equations are dimensionless, but can be rescaled in time as described in the chimera-on-a-ring case.

To investigate the robustness of the chimera state(s) to noise (in Figure 11), an uncorrelated white noise term was added to each oscillator for the chimera on a ring:

where ζ_i(t) has mean 0, and standard deviation σ.

Since the oscillators are periodic and oscillate between 0 and 2π, we can transform the oscillators' time-series (Figures 2B, 3A) into a raster plot. Every time the oscillator's phase ϕ_i = 0 a spike is drawn. The resulting raster plot from the aforementioned time-series is shown in for the chimera on a ring and in Figure 3C for the two populations chimera.

The mean phase velocity (Omelchenko et al., 2013) for a given oscillator with phase θ_i is defined as :

where M_i is the number of complete rotations around the origin performed by the ith oscillator during the time interval Δt = 1000. It acts as a measure of the oscillating frequency for each oscillator. Given a ring of N oscillators, it is denoted as Ω = Ω_i, with i = 1, 2, …N. Having different mean phase velocities for the synchronized and unsynchronized domain is typical for chimera states. In particular, for ring-like topologies it is common to have an arc-like profile of mean phase velocities for the unsynchronized domain, denoted as Ωu=Ωju,…Ωmu (where u stands for unsynchronized and m is the total number of unsynchronized oscillators) and equal mean phase velocities for the synchronized domain (Gu et al., 2013; Kemeth et al., 2016; Sawicki et al., 2017) given by Ωs=Ωjs,…Ωns (where s stands for synchronized and n is the total number of synchronized oscillators). Since the oscillators are synchronized they have the same mean phase velocity Ω^s, therefore we can simplify Ω^s to a unique value given by Ω_s. We will use Ω_s to identify which oscillators synchronize and which do not, given that the synchronized domain oscillates at a slower pace than for the unsynchronized one Ω^s < Ω^u ∀j∈m. In order to identify Ω^s, one can simply compute the minimum of Ω.

For the two-population chimera, both domains (synchronized and unsynchronized) have equal mean phase velocities (different between domains but equal within). Also for that topology, the synchronized population oscillates at a slower pace than for the unsynchronized one: Ω^s < Ω^u.

The spiking neural network consists of coupled Izhikevich neurons (Izhikevich, 2003), with their dynamics given by the following equations:

The quantity v_i is the voltage variable. Neuron i fires a spike when v_i reaches a voltage peak v_peak and it is instantly reset to a potential v_reset. The adaptation current is given by u_i, which increases an amount d_u every time a spike is fired and which in turn slows down the production of spikes. The current I_i is given by I_i = I_bias+s_i, where I_bias is a fixed value and s_i are the synaptic currents for neuron i, given by

where N is the total number of neurons. The matrix ωij0 controls the magnitude of the postsynaptic currents arriving at neuron i from neuron j. The parameter C represents the membrane capacitance, the parameters v_r and v_t denote the resting and the threshold membrane potential, respectively. The parameter a is an equivalent of the time constant for the adaptation current u_i. The parameter b controls the resonance properties of the model and k controls the half-width of the action potentials. The numeric parameters of the model are listed in Table 1, we used the same parameters as in Nicola and Clopath (2017). The spikes are filtered with a double exponential synapse, given by:

where τ_r = 2 ms is the synaptic rise time, τ_d = 20 ms is the synaptic decay time and t_jk is the time at which the neuron jth fired spike kth. For other synapse types, see Nicola and Clopath (2017).

The output of a spiking neural network is defined as:

where d_j is an m-dimensional vector known as the linear decoder for the firing rate. Here, we want to train the network such that:

where x = (x₁, x₂, …, x_m) are the desired dynamics or the supervisor that the network should mimic. Since the oscillators' phases γ and ϕ are discontinuous and wrapped around the interval [0, 2π), the following supervisor for the chimera was used: x = (cosϕ, sinϕ, cosγ, sinγ). With 2n (n = 3) oscillators, this results in a m = 4n = 12 dimensional supervisor, see Masoliver et al. (2022) for details.

In order to achieve Equation 13 we use the FORCE method (Sussillo and Abbott, 2009), which adds a second set of weights Qηi·djT when defining the synaptic currents. Equation 9 can be rewritten as:

The FORCE method has three phases, the pre-learning, the learning and the post-learning. In the pre-learning phase, the initial synaptic connection matrix wij0 initializes the neurons' dynamics into a well-known high-dimensional chaotic regime (Ostojic, 2014; Harish and Hansel, 2015). The matrix is static and sparse with each element drawn from a normal distribution with mean 0 and variance 1Np2, where p is the sparsity degree (set to 90% sparse or p = 0.1). The variable G controls the network's chaotic behavior and its value depends on the neuronal model, see Nicola and Clopath (2017) for a detailed explanation. Here, we set G = 1.5 × 10³.

The learning phase involves a second set of weights, given by Qηi·djT. Where the parameter Q scales the encoding vector η_i, which has been drawn randomly and uniformly from [−1, 1]^m (where m is the dimensionality of the supervisor). By increasing Q, the feedback applied to the network is strengthened. A value of Q = 1.4 × 10³ was used for all simulations.

In the learning phase, the FORCE method enforces the aforementioned constrain x^≈x by changing d_i online (i.e., as the network is being simulated) with the Recursive Least Squares (RLS) (Sussillo and Abbott, 2009). RLS has an online solution for the optimal d, the one that minimizes the squared error e between the network output x^ and the complex signal or supervisor x. RLS updates to d at each time step n are:

where d₀ = 0 and P₀ = I_n/λ. The parameter λ controls the rate of the error (Sussillo and Abbott, 2009) and we set it to λ = 1. The parameter I_n is a N×N identity matrix.

The third step of the FORCE method is the post-learning phase. RLS is turned-off and the weight matrix Qηi·djT is no longer dynamic but static. The FORCE method is successful if the network is able to reproduce the supervisor for a fixed d.

Finally, Dale's law can also be enforced in trained spiking neuronal networks. In Dale's law, a neuron can only be either inhibitory or excitatory, not both. Dale's Law was enforced by constraining ω to the inhibitory/excitatory nature of each individual neuron. If neuron i is inhibitory (excitatory), all of its outgoing connections will be negative (positive): ω1i1,⋯,ωNi1<0 (ω1i1,⋯,ωNi1>0). We first define ωij0 such that ωij0≥0 ∀j ∈[0, N] (the first half of the population of neurons only projects positive weights, i.e., excitatory neurons) and ωij0≤0 ∀j ∈[0, N] (the second half of the population of neurons only projects negative weights, i.e., inhibitory neurons). Second, the trained matrix Qηi·djT is limited to project either positive or negative weights. We obtain that by defining η as η = η₋+η₊, where η₋ and η₊ are unequivocally defined as negative and positive matrices, respectively. And finally, d^T is defined such that d_ij≥0 ∀i∈[0,N2] and d_ij ≤ 0 ∀i∈[N2,N]. For the exact implementation refer to Additional Information where the link to the code is available and for more details see Nicola and Clopath (2019, 2017).

To investigate if chimera dynamics are a potential mechanism for the neuronal dynamics associated with the hippocampal theta oscillation, chimera dynamics were first simulated in pre-existing models to test the hypothesis that parameter ranges that exhibit hippocampal-like dynamics (i.e., phase precession and sequential content) could readily be determined.

Two standard model versions, each generating a different chimera state—a chimera on a ring and a two-population chimera, respectively—were considered. The chimera on a ring arises for N = 500 non-locally coupled identical Kuramoto oscillators (see Figure 2A and Methods, Equation 1), whereas the two-population chimera arises for two weakly coupled populations, each one formed by 3 globally coupled Kuramoto oscillators (see Figures 2E, 3E Methods, Equations 2, 3). Depending on the parameter values and on the initial conditions, both network topologies can display different dynamics: Either a fully synchronized state where all oscillators are in phase or chimera states where one sub-population of neurons is synchronized while the other sub-population oscillates asynchronously (see Supplementary Video S1).

First, the model parameters of the dynamical equations (Equation 1 and Equations 2, 3, respectively) were set to well known or classical parameter regimes where chimera dynamics readily emerge (Laing, 2009; Panaggio et al., 2016). The parameters A and β affect the coupling strength and the phase difference, respectively, and take different values for the two different systems, see Table 2 for details. For the chimera on a ring, the synchronous subpopulation of oscillators is non-static, and drifts slowly around the ring. Oscillators drift in and out of the synchronous sub-population, while the others oscillate asynchronously (see Supplementary Video S2, Figures 2B, C) with a narrow distribution of frequencies (Figure 2D). In contrast, for the two-population chimera, the chimera state is static: one population fully synchronizes (triangles in Figure 2F) while the other one does not (diamonds in Figure 2F). Unless the system is perturbed, the synchronized and unsynchronized populations remain fixed, each with a fixed oscillation frequency (Figure 2H). The identity of the synchronized or unsynchronized population depends on the initial conditions (Figures 2F, G). The synchrony profile between the two populations can be exchanged by externally perturbing the system, where the synchronized and unsynchronized populations swap. For example, in Figure 4, the triangle population is synchronized before a perturbation, and after a perturbation, the oscillators move to an asynchronous regime (vice versa for the diamond population).

For the chimera on a ring, the synchronized and unsynchronized populations drift (Supplementary Videos S1, S2 and Figures 2B, C). An external perturbation, in this case, is not necessary to change the oscillators' synchrony profile. The identity of the neurons that constitute the synchronized population slowly drifts around the ring as a slowly moving traveling wave. As the drift's period is much larger than the oscillations' period, we can study the differences between the two domains, synchronized and unsynchronized (see ref. Abrams and Strogatz (2006) for details on the drift).

With the classical chimera dynamics reproduced, we investigated how to explicitly draw a mapping between the Kuramoto networks, specifically the ring network (Figure 5A), and hippocampal dynamics. Each neuron has more complex dynamics than those of a Kuramoto oscillator which is a simple oscillator where the frequency is integrated to arrive at the oscillator phase (Figure 5B). Specifically, neurons emit spikes when their inputs are sufficient to reach a threshold. Thus, each oscillator's continuous time-series was converted into spike trains via a Poincare Map. Each time any Kuramoto oscillator's phase reaches 2π, a “spike” is generated at the time that this occurred (θj(t*)=2π) as depicted in Figure 5C (see Methods for details). With a spike-generating Poincare map, the “spikes” generated by the chimera on a ring (Figure 5D) and for the two-population chimera (Figure 3) can be analyzed.

In order to measure phase-precession, an equivalent component to the hippocampal local field potential in the Kuramoto network is required. The hippocampal LFP is a macroscopic observable that is a complex synthesis of propagating action potentials, and synaptic activity. While there is some debate as to whether or not the LFP is reflective of underlying oscillations, or indeed organizes the timing of spikes, it is convenient to measure other oscillation frequencies (i.e., the oscillations of individual units) relative to the LFP (Buzsáki et al., 2012). During in vivo recordings, the hippocampal LFP is typically converted into a phase (for example with a Hilbert transform). Interneurons and sometimes pyramidal neurons lock to phases of the hippocampal LFP, while other pyramidal neurons fire at a slightly faster rate.

Given the locking of synchronized sub-populations to the hippocampal LFP, a phenomenological LFP can be computed as follows: the cosine of the phase of each oscillator in the synchronized population is obtained (cosϕ_j) and globally averaged over the synchronized population. The LFP can also be computed as the mean over all oscillators (both synchronized and unsynchronized): 1N∑i=0Ncosϕi. Both methods of computing the LFP product qualitatively similar results (Figures 6A, B). We note that there are more direct, biophysically based models of LFPs considered in the literature (Mazzoni et al., 2015).

Interestingly, we observed phase advancement from the unsynchronized oscillators when compared to the synchronized ones (Figures 5D, C). While this is similar in principle to phase precession, where the unsynchronized pyramidal neurons fire slightly faster than the local-field-potential, the frequency ratio between the synchronized oscillators and the unsynchronized is different from those observed experimentally. For example, in Figure 5D, for every synchronized spike we get approximately three unsynchronized ones, which roughly gives us a ratio of ≈0.33. In the hippocampus, pyramidal neurons fire at approximately 9 Hz, while the theta oscillation observed in the LFP is approximately 8 Hz, which yields a a ratio of ≈ 0.88.

However, chimera states are solutions to coupled oscillator networks that are parameter dependent. Indeed, this is similar to limit cycles, chaotic solutions, or fixed points. The precise characteristics of all of these solutions depend on the chosen parameters for the underlying network. For example, in Panaggio et al. (2016), the chosen system parameters yield three unsynchronized spikes to one synchronized spike ratio as mentioned above. This ratio is also approximately the mean-phase velocity ratio between synchronized and unsynchronized populations. As another example, in Gu et al. (2013), it was found that the difference in mean-phase velocities can be very small, with only a 2% difference in the frequencies between the synchronized and unsynchronized populations. Finally, in Sawicki et al. (2017), the mean-phase velocity ratio is more intermediate in range, between 50%–100%. In some cases, the synchronized population can also oscillate faster than the unsynchronized population. All these differences in chimera dynamics arise from differences in the underlying models and model parameters. In the next section, we show that two well established chimera-capable models (Panaggio et al., 2016; Abrams and Strogatz, 2006) can yield phase-precession like spiking dynamics as in the hippocampus.

Next, we investigated if the parameters in both models could be varied to both preserve the chimera state, and obtain a frequency ratio closer to that of hippocampal phase precession (≈0.88). Accomplishing this in both models would indicate that one can generically obtain hippocampal-like dynamics in chimera systems. To start, the intrinsic frequency parameter (ρ) was varied in Equations 1–3. This acts as the fundamental driving force for an oscillator and causes the oscillator to intrinsically oscillate when no coupling is present. Thus, it is directly comparable to the applied current I typically considered in neuron models as higher applied currents lead to faster neuronal oscillations.

As ρ was varied, the oscillating frequency for each oscillator was quantified as follows: the mean phase velocity Ω_i was computed for oscillator i to determine its frequency. As the driving frequency ρ interacts with the coupling in a non-trivial way, the frequencies must be computed numerically. For a given oscillator i and a given amount of time Δt, the number of rotations around the origin (or equivalently, the number of spikes fired) was summed and multiplied by 2π (see methods and Equation 5 for details). This was then divided by Δt to yield the rotations.

To see if the chimera dynamics could mimic hippocampal observations, we focused on the mean phase velocity ratio 〈Ω〉_ratio. This ratio was computed as the average of the mean phase velocity ratio between the synchronized and unsynchronized populations as a function of ρ, as shown below (see Methods for details):

As the intrinsic oscillation frequency (ρ) increases, the oscillation frequency of both the synchronized and unsynchronized oscillators in the coupled network increases, but the frequency difference between the synchronized and unsynchronized domains decreases. This was quantified for ρ = 1.8 (Figures 7A, C) and ρ = 2.8 (Figures 7B, D) and, more generally, for the mean phase velocity ratio as a function of ρ (Figure 7E). As ρ was varied, Ω^u varied over a range which was bounded by a minimum Ω_min and a maximum value Ω_max. For ρ = 1.8, (Ω_min, Ω_max) = (1.056, 1.565) while for ρ = 2.8, they increase to (Ω_min, Ω_max) = (2.055, 2.545) and we achieve 〈Ω〉_ratio≈0.88 for that value (Figure 7E). As ρ is increased further past this value, the ratio slowly increases until the chimera state collapses and all oscillators synchronize (Figure 8).

Next, we tested if this was a generic response by considering the two-population chimera model (Figure 9). Once again, we found that the 〈Ω〉_ratio≈0.88 can occur for a specific ρ due to the slow gradual increase in 〈Ω〉_ratio as a function of ρ. Thus, the phase precession regime of classical chimera models is seemingly robust and generic.

Finally, we investigated what the net impact of the coupling was. That is, we considered how the mean phase velocity for both domains (synchronized and unsynchronized) and for both network topologies compares to the mean phase velocity of an uncoupled oscillator. In the latter case, the mean-phase velocity is given by

Interestingly, regardless of the network topology, the net effect of the coupling was always inhibitory: The oscillators fire at a faster frequency when uncoupled, rather than when coupled into a chimera state in both network topologies (Figure 10).

Finally, we investigated the impacts of noise on the Kuramoto system (Figure 11) in the phase precessing regime. We found that injecting white noise into each oscillator for the chimera on a ring, with mean 0 and standard deviation σ did not substantially impact the results, with similar phase precession dynamics and mean-phase velocities when σ was small.

Chimera dynamics in networks of Kuramoto oscillators with different network topologies can be altered by changing one parameter, the intrinsic driving frequency, to mimic a hippocampal-like phase precession regime. Despite the general nature of these results, the Kuramoto-oscillator network is phenomenologically different from the neurons and synaptic connections in the hippocampus in addition to having the property that all of the oscillators are homogeneous. Thus, we sought to determine if embedding a chimera state into a spiking-neural-network would still yield hippocampal phase precession, and a global theta-oscillation.

A chimera state can be “embedded” in a recurrent neural network by training the network to output a chimera, as seen in (Masoliver et al. 2022). To test if such an embedding was applicable in a spiking network, we trained a spiking neural network using the FORCE method (Sussillo and Abbott, 2009; Nicola and Clopath, 2017) to output the two-population chimera, described by Equations 2, 3. This network was constrained with Dale's law, with a proportion of the neurons being excitatory, and the rest inhibitory. Initially, the individual neurons (modeled using the Izhikevich model, see Methods for details) are sparsely connected [to support the learning process (Sussillo and Abbott, 2009)] with a set of static weights Gω₀ which initiate the neurons' rate r(t) into a high-dimensional chaotic regime. During the training period a second set of weights Qηd^T is added to Gω₀ and changes the connections between neurons such that the network's output (defined as d^Tr) equals the desired dynamics. The desired dynamics or supervisor are cosines of the phases of a two-population Kuramoto oscillator network in the chimera regime. At each time step, d is updated using the Recursive Least Squares (RLS), which minimizes the sum-squared difference between the network output and the two-population chimera. The network has learned when for a fixed value of d it is able to mimic the desired chimera dynamics (Figure 12). A specific example is shown in Figure 13. We remark that while the chimera state supervisors have homogeneous oscillators, the trained neurons whether in a rate or spiking network are heterogeneous, as they receive a combination of randomly generated, and trained weights which alters their activity levels and how they encode the chimera dynamics.

In order to assess if the individual neurons of the spiking network show phase precession, the voltage traces were transformed into phases (Figure 13C). The spike times were transformed into phases by using a linear interpolation to approximate the phase at each spike time with the phase of one of the synchronized components of the network output ϕ^syn(t). We found that phase precession occurred generically for many of the neurons sampled, where neurons displayed decreasing burst phases based on subsequent cycles (Figure 13D). However, some of the neurons were primarily phase locked (Figure 13E).