Collective neuronal activity in the brain synchronizes during rest and desynchronizes during active behaviors, influencing cognitive processes such as memory consolidation, knowledge abstraction, and creative thinking. These states involve significant modulation of inhibition, which alters the excitation–inhibition (EI) balance of synaptic inputs. However, the influence of the EI balance on collective neuronal oscillation remains only partially understood. In this study, we introduce the EI-Kuramoto model, a modified version of the Kuramoto model, in which oscillators are categorized into excitatory and inhibitory groups with four distinct interaction types: excitatory–excitatory, excitatory–inhibitory, inhibitory–excitatory, and inhibitory–inhibitory. Numerical simulations identify three dynamic states—synchronized, bistable, and desynchronized—that can be controlled by adjusting the strength of the four interaction types. Theoretical analysis further demonstrates that the balance among these interactions plays a critical role in determining the dynamic states. This study provides valuable insights into the role of EI balance in synchronizing coupled oscillators and neurons.

Collective neuronal oscillations are crucial for processes such as memory consolidation, knowledge reorganization, and creative thinking (Landmann et al., 2014; Lewis et al., 2018). Neuronal activity tends to synchronize during rest and desynchronize during active behavior (Harris & Thiele, 2011). Notably, significant differences in inhibitory neuronal activity and synaptic conductance are observed between synchronized and desynchronized states compared with excitatory activity and conductance (Alfonsa et al., 2022; Haider et al., 2013; Miyawaki & Diba, 2016; Mizuseki & Buzsáki, 2013). These differences contribute to shifts in the excitation–inhibition (EI) balance. However, the influence of EI balance on collective neuronal oscillations remains poorly understood.

Coupled oscillator models are fundamental tools for investigating collective synchronization and have been widely used to simulate various neuronal oscillation phenomena (Bick et al., 2020), such as interregional synchronization (Ódor & Kelling, 2019), traveling waves (Zhang & Jacobs, 2015; Zhang et al., 2018), central pattern generators (Marder & Bucher, 2001), gamma oscillations (Tiesinga & Sejnowski, 2009), and continual attractors composed of grid cells (Campbell et al., 2018). The Kuramoto model, one of the simplest frameworks for studying collective synchronization, incorporates two key components: the natural frequency of each oscillator and the interaction strength between oscillators (Kuramoto, 1975, 1984). Numerous variants of the Kuramoto model with different connection structures, interaction rules, and natural frequency distributions have been developed, yielding complex dynamics (Acebrón et al., 2005). In the Kuramoto model, excitatory and inhibitory interactions are viewed as opposing forces that drive synchronized oscillations through EI-based feedback mechanisms (Montbrió & Pazó, 2018). Nonetheless, the role of EI balance in this context has not been fully explored.

We integrated the EI balance into the Kuramoto model by categorizing oscillators into excitatory and inhibitory groups, each with distinct interaction strengths. This modified framework, termed the EI-Kuramoto model, was used to investigate dynamic states by adjusting the interaction strength parameters. The model demonstrated transitions between synchronized, bistable, and desynchronized states, driven by variations in inhibitory interaction strength, effectively mimicking the dynamics observed in the biological brain.

2.1  EI-Kuramoto Model

The original Kuramoto model (Acebrón et al., 2005; Kuramoto, 1975, 1984) is expressed as:
where θj(t) is the phase of unit j at time t, ωj is the natural frequency of unit j following the arbitrary (gaussian in this study) probability distribution g(ω), K is the interaction strength among the oscillation units, and N is the number of units. Then, the order parameter R and mean phase ψ are defined using the mean vector of the oscillators as
where i denotes an imaginary unit. Given Euler's formula and equation 2.2,
Considering equation 2.3, equation 2.1 can be changed to
The EI-Kuramoto model classifies oscillator units into excitatory (attractive) and inhibitory (repulsive) categories, referred to as excUnits and inhUnits, respectively. These terms are explicitly defined to clarify that the units represent excitatory and inhibitory elements within the model rather than actual neurons or synapses. The model incorporates EI balance by employing the simplest form of opposing sine-wave functions with a phase lag of zero radians, where excUnits engage in purely attractive interactions and inhUnits engage in purely repulsive interactions. Accordingly, the terms excitatory and inhibitory are used interchangeably with attractive and repulsive in this model. This approach, which uses sine functions, differs from the EI-Kuramoto model proposed in a previous study that employed cosine functions (Montbrió & Pazó, 2018). The interaction strength K is categorized into four quantifiable parameters: Kee (attractive strength between excUnits), Kii (repulsive strength between inhUnits), Kei (attractive strength from excUnit to inhUnit), and Kie (repulsive strength from inhUnit to excUnit). Kee, Kii, Kei, and Kie are distinct positive real numbers that are uniformly applicable across all unit pairs. Consequently, equation 2.4 becomes
for excUnits and
for inhUnits. The order parameter and mean phase for each population of excUnits and inhUnits are expressed as, respectively,
and

Given the rotational symmetry, when the mean natural frequencies of excUnits and inhUnits are equal, ω¯=ω¯exc=ω¯inh, the mean natural frequencies can be set to 0 rad/s. This is equivalent to rotating a coordinate by the mean natural frequency ω¯ (Acebrón et al., 2005).

In the numerical simulations, the parameters were set to Nexc=80 and Ninh=20, aligning with approximately 20% of all cortical and hippocampal neurons being inhibitory (Andersen et al., 2006; Braitenberg & Schüz, 1998), with a total simulation duration T=50 s for exploring the K parameter (see section 3.1), T=3000 s for cyclic inhibitory modulation (see section 3.3), time step dt=0.1 s, and the g(ω) of both excUnits and inhUnits populations following a gaussian distribution g(ω)N(0,0.09). The time range for calculating the mean values of the order parameters Rexc and Rinh was set from T/2 to T to exclude initial fluctuations.

In section 3.1, we derive the total interaction strengths for excUnit and inhUnit, Fexc(θ,ψ) and Finh(θ,ψ), from the summation of the second and third terms in equations 2.5 and 2.6. In this case, θ represents the phase of a single representative unit. We assume ψψexcψinh, because excUnits attract both themselves and inhUnits to a nearly identical phase (see Figures 1b and 2b2f). Fexc(θ,ψ) and Finh(θ,ψ) can be written as:
Figure 1:

Dynamics of the EI-Kuramoto model. (a) Schematic representation of the model. Oscillator units are categorized into excUnits (warm-colored circles) and inhUnits (cool-colored circles). These units, each characterized by distributed natural frequencies (indicated by warm and cool color gradients), oscillate and interact. Kee is the attractive interaction strength between excUnits (red arrows). Kei is the attractive strength from excUnit to inhUnit (magenta arrow). Kii is the repulsive strength between inhUnits (blue arrows). Kie is the repulsive strength from inhUnit to excUnit (cyan arrow). (b) Phase of individual units (left) and R-values (right) as a function of time under three conditions (synchronized, bistable, desynchronized). The left panels show phase color plots of individual units over time, with excUnits indexed from 1 to 80 and inhUnits from 81 to 100. The right panels depict Rexc and Rinh in red and blue, respectively. The K parameters are as follows: synchronized: Kee=3, Kei=1, Kie=1, Kii=1; bistable: Kee=1, Kei=3, Kie=3, Kii=1; desynchronized: Kee=1, Kei=3, Kie=3, Kii=3. (c) Time-averaged Rexc derived from the K parameter validation. Magenta, light green, and yellow circles correspond to parameters that approximately match those in the synchronized, bistable, and desynchronized states shown in panel b, respectively. The left panel provides a zoomed-in view of the red-framed region in the right panel.

Figure 1:

Dynamics of the EI-Kuramoto model. (a) Schematic representation of the model. Oscillator units are categorized into excUnits (warm-colored circles) and inhUnits (cool-colored circles). These units, each characterized by distributed natural frequencies (indicated by warm and cool color gradients), oscillate and interact. Kee is the attractive interaction strength between excUnits (red arrows). Kei is the attractive strength from excUnit to inhUnit (magenta arrow). Kii is the repulsive strength between inhUnits (blue arrows). Kie is the repulsive strength from inhUnit to excUnit (cyan arrow). (b) Phase of individual units (left) and R-values (right) as a function of time under three conditions (synchronized, bistable, desynchronized). The left panels show phase color plots of individual units over time, with excUnits indexed from 1 to 80 and inhUnits from 81 to 100. The right panels depict Rexc and Rinh in red and blue, respectively. The K parameters are as follows: synchronized: Kee=3, Kei=1, Kie=1, Kii=1; bistable: Kee=1, Kei=3, Kie=3, Kii=1; desynchronized: Kee=1, Kei=3, Kie=3, Kii=3. (c) Time-averaged Rexc derived from the K parameter validation. Magenta, light green, and yellow circles correspond to parameters that approximately match those in the synchronized, bistable, and desynchronized states shown in panel b, respectively. The left panel provides a zoomed-in view of the red-framed region in the right panel.

Close modal
Figure 2.

Detailed model dynamics in the bistable state. (a) Plots of Rexc (red) and Rinh (blue) in the bistable state. This plot shows an expanded period range from 26.5 to 35.5 s, as indicated in the middle-right panel of Figure 1b. The dashed vertical lines correspond to the time points referenced in panels b–g. (b–g) The upper panels display the phases of each unit at specific time points. The warm- and cool-colored circles represent the excUnits and inhUnits, respectively. The central red and blue lines represent the mean vector angles and amplitudes (R-values) of the excUnits and inhUnits, respectively. The lower panels display the interaction strengths at each time point. The left panels show the total interaction strength acting toward excUnits: Fexc (black), -KeeRexcsin(θ-ψ) (red), and KieRinhsin(θ-ψ) (cyan). The right panels represent the total interaction strength acting toward inhUnits: Finh (black), -KeiRexcsin(θ-ψ) (magenta), and KiiRinhsin(θ-ψ) (blue). In panel g, the lower panels are not present, as ψexcψinh in this case.

Figure 2.

Detailed model dynamics in the bistable state. (a) Plots of Rexc (red) and Rinh (blue) in the bistable state. This plot shows an expanded period range from 26.5 to 35.5 s, as indicated in the middle-right panel of Figure 1b. The dashed vertical lines correspond to the time points referenced in panels b–g. (b–g) The upper panels display the phases of each unit at specific time points. The warm- and cool-colored circles represent the excUnits and inhUnits, respectively. The central red and blue lines represent the mean vector angles and amplitudes (R-values) of the excUnits and inhUnits, respectively. The lower panels display the interaction strengths at each time point. The left panels show the total interaction strength acting toward excUnits: Fexc (black), -KeeRexcsin(θ-ψ) (red), and KieRinhsin(θ-ψ) (cyan). The right panels represent the total interaction strength acting toward inhUnits: Finh (black), -KeiRexcsin(θ-ψ) (magenta), and KiiRinhsin(θ-ψ) (blue). In panel g, the lower panels are not present, as ψexcψinh in this case.

Close modal

The equilibrium, where F=0, is attractive if F'lt;0 and repulsive if F'gt;0.

2.2  Cyclic Inhibitory Modulation

In section 3.3, we present a simulation scenario where the repulsive interaction strength from inhUnits fluctuates sinusoidally over time. Kii and Kie are implemented with the same parameter Kito simplify the simulation.
where Ki(t) is modulated as a sine function:
where As (= 1.485), ωs (= 0.004π rad/s), and Cs (= 2.505) are the amplitude, frequency, and baseline constant for the sine modulation, respectively. The other interaction strengths are held constant to simplify the simulation; Kee=1 and Kei=3.

3.1  Numerical Simulations

In the EI-Kuramoto model, oscillator units were categorized as excUnits and inhUnits. EI balance was implemented using the simplest form of opposing forces, with excUnits exhibiting purely attractive sine-shaped interactions and inhUnits displaying purely repulsive sine-shaped interactions. Four distinct types of interaction strengths were defined: excUnit-excUnit attraction, Kee; excUnit-inhUnit attraction, Kei; inhUnit-excUnit repulsion, Kie; and inhUnit-inhUnit repulsion, Kii. The K values were determined independently (unless otherwise noted) and applied uniformly to all unit pairs to preserve the simplicity of the model (see Figure 1a and section 2.1).

In this study, we did not explicitly define each oscillator unit as representing specific neurons or neuronal populations, as our focus was on investigating the dynamics of simple oscillators in modulating EI balance. Nor were the specific frequencies of neuronal oscillations a primary focus. However, each oscillator in the model can be conceptually associated with the burst activity of a neuron or a group of neurons within the slow-wave (0.5–4 Hz) neuronal oscillation, which undergoes substantial changes during the sleep–wake cycle.

First, we explored the parameter space of the four K parameters (Kee, Kii, Kei, and Kie) and observed the resulting model dynamics. Collective synchronization was assessed using the order parameter R, which quantifies the magnitude of the average vector across oscillator units.

Consequently, the dynamics were categorized into three distinct states: synchronized (high R), bistable (alternating high and low R), and desynchronized (low R) (see Figure 1b and supplementary videos 1–3). A systematic exploration of the parameter space for Kee, Kii, Kei, and Kie, revealed substantial differences in dynamic states depending on the specific combination of parameter values (see Figure 1c).

To further analyze the system, we qualitatively examined the dynamics. Figure 2 highlights the dynamics of the bistable state (as shown in Figure 1b, bistable). At most time points, the phases of the mean vectors for excUnits and inhUnits are closely aligned ψψexcψinh (see Figures 2b2f). This alignment arises from the attractive force exerted by the excUnits on both themselves and the inhUnits, pulling them toward the same phase position. We extracted the interaction terms F (Fexc for excUnits and Finh for inhUnits; see section 2.1) and tracked their changes. If the slope of F at F=0 (x-intercept) is negative, F attracts θ-ψ toward the x-intercept. If the slope is positive, F repels θ-ψ from the x-intercept.

The accumulation of excUnits at a specific phase position triggers a repetition of cycles, alternating between high and low R values (see Figure 2b). During this stage, the excUnits attract both themselves and the inhUnits toward the mean phase ψ. Concurrently, the inhUnits, which are gradually attracted toward ψ, amplify the repulsive strength over time (see Figure 2c). At a certain point, the repulsive force from the inhUnits to the excUnits becomes stronger than the attractive force between the excUnits, causing a shift in the total interaction strength from attraction to repulsion for the excUnits (see Figure 2d). As a result, the excUnits disperse, diminishing Rexc and the attractive force (see Figure 2e). As the attractive force weakens and only the repulsive force remains, the inhUnits also begin to disperse (see Figure 2f). Eventually the excUnits accumulate at another phase position (ψexcψinh in this case; see Figure 2g) and start attracting both the excUnits and inhUnits toward the same mean phase ψ, thus returning to a state similar to that shown in Figure 2b.

When Kee is sufficiently strong, the repulsive force on the excUnits does not surpass the attractive force, preventing the model from reaching the state shown in Figure 2d, where the repulsive forces on the excUnits exceed the attractive forces. As a result, the model continues to accumulate oscillator units, resembling the states in Figures 2b and 2c, ultimately reaching a synchronized state.

When Kii is strong, the accumulation of the inhUnits (depicted in Figures 2b2d) decreases due to mutual repulsion among the inhUnits. When the feedback repulsion to the excUnits through Kei and Kie is weak, the repulsive force on the excUnits is reduced, leading to the synchronization of the excUnits, while the inhUnits remain desynchronized. Conversely, strong feedback repulsion with large Kei and Kie hinders the accumulation of the excUnits, resulting in both the excUnits and inhUnits being desynchronized.

3.2  Bifurcation Analysis

We aimed to identify the parameter conditions for bifurcation points between dynamic states through theoretical analysis. Introducing the density distribution ρ(θ,t,ω) is useful for assessing the dynamic states. This distribution defines the distribution of phase θ at each time t and natural frequency ω, assuming an infinite number of oscillator units (N) and treating them continuously (Bick et al., 2020; Montbrió & Pazó, 2018; Strogatz & Mirollo, 1991). The sum of the density distribution at each time t and each natural frequency ω is 1, that is,
The evolution of density distribution can be expressed using the transport equation (Bick et al., 2020):
Based on equations 2.5 and 2.6, the velocity v, which represents the phase change over time, for the excUnits and inhUnits is given as
The order parameter R and mean angle ψ can be redefined as
To simplify the analysis, g(ω) is approximated as a Dirac delta function, corresponding to the limit of the infinitely narrow peak:
We can ignore ω0 and rewrite equation 3.5 as
When the phase is uniformly distributed, ρ0=1/2π in the N limit, and the desynchronized state maintains a uniform distribution of phases. This implies that when a small perturbation εη(ε1) from a uniform distribution exists,

The evolution of η is negative in the desynchronized state. When η has an imaginary number, the η and ρ oscillate in time.

By substituting equation 3.8 into equation 3.2, the evolution equation of η becomes
By substituting equation 3.8 into equation 3.7, the order parameter can be expressed as
Upon the substitution of equation 3.11 into equations 3.3 and 3.4, followed by the substitution of the outcomes into equation 3.10, the evolution equations on the order of ε(O(ε)) change to
Using the Fourier method, η can be expressed as an infinite summation of different oscillation terms because η is real and 2π-periodic in θ. We assume that ψexc=ψinh because the excUnits attract both types of units, leading ψexc and ψinh to similar phases in most cases (see Figure 2). η changes to
where c is the coefficient of the first-order oscillation, c* is the complex conjugate of c, and η is the second and higher harmonic of η.
When equation 3.14 is substituted into equation 3.11, it becomes
Given equation 3.15 and R1cos(ψ-θ)= Re [R1eiψe-iθ], we obtain
where c.c. denotes the complex conjugate of the first term.
We then derive the evolution equation for each c by substituting equations 3.14 and 3.16 into equations 3.12 and 3.13, and equating the coefficients of eiθ on both sides of the resulting equation. Consequently, we obtain

The evolution equations for c* need not be written because they are simply conjugates of equations 3.17 and 3.18. The evolution equation for η was zero at O(ε).

We rewrote equations 3.17 and 3.18 using a Jacobian matrix (Izhikevich, 2010):
The eigenvalues λ of equation 3.19 can be expressed as
To obtain a convergent solution (corresponding to the desynchronized state), the two eigenvalues λ± must be negative. Based on equation 3.20, it can be deduced that

Note that we focus only on λ+={(Kee-Kii)+2D}/4 because Re[D]0; thus, + D-D, and λ+λ-.

Furthermore, from equation 3.22, we have (Kee-Kii)2>(Kee+Kii)2-4KieKei, leading to
When at least one λ (in this case λ+) is nonnegative, the condition Dlt;0 must be met to obtain an oscillating solution (corresponding to the bistable state). Based on equation 3.21 and the condition that all K values are nonnegative, we obtain
Under this assumption, it can be observed that the bifurcation criteria depend solely on the balance between the interaction strengths. These results are summarized in Figure 3.
Figure 3:

Bifurcation of the dynamic states. (a) Bifurcation criteria. The dynamic state converges to a desynchronized state when λ+lt;0 (yellow rectangle), to a bistable state when λ+gt;0 and Dlt;0 (green rectangle), and to a synchronized state when λ+gt;0 and Dgt;0 (navy blue rectangle). The black lines in each quadrant represent the schematic time course of the order parameter R. (b) Bifurcation of dynamic states based on the values of Kee, Kii, and KeiKie.

Figure 3:

Bifurcation of the dynamic states. (a) Bifurcation criteria. The dynamic state converges to a desynchronized state when λ+lt;0 (yellow rectangle), to a bistable state when λ+gt;0 and Dlt;0 (green rectangle), and to a synchronized state when λ+gt;0 and Dgt;0 (navy blue rectangle). The black lines in each quadrant represent the schematic time course of the order parameter R. (b) Bifurcation of dynamic states based on the values of Kee, Kii, and KeiKie.

Close modal

3.3  Cyclic Inhibitory Modulation

Inhibition differs significantly across neuronal collective synchronization states (Alfonsa et al., 2022; Haider et al., 2013; Iwase et al., 2024; Miyawaki & Diba, 2016; Mizuseki & Buzsáki, 2013). In this section, we demonstrate the ability of the EI-Kuramoto model to transition between dynamic states (synchronized, bistable, and desynchronized) by modulating the strength of the inhUnits (see Figure 4 and supplementary video 4). To simplify the simulation, we consolidate the repulsive interaction strengths Kii and Kie into a single parameter. Ki(Ki=Kie=Kii) and gradually vary the Ki value sinusoidally (see Figure 4, upper panel). Rexc and Rinh consistently increase or decrease according to the magnitude of Ki. The decrease in Ki is associated with an increase in both Rexc and Rinh. The dynamic states were accurately predicted using D and λ+, despite a minor timing discrepancy between the theoretical predictions (backgrounds) and numerical simulations (red and blue lines). These results demonstrate that despite its simplicity, the EI-Kuramoto model can effectively capture dynamic shifts in response to changes in repulsive interaction strength, thereby reflecting the neuronal dynamics associated with inhibitory function in the biological brain.
Figure 4:

Cyclic inhibitory Ki modulation. The upper panel shows the time course of Ki, and the middle and lower panels show Rexc and Rinh, respectively. The blue, green, and yellow background shadows represent synchronized, bistable, and desynchronized states, respectively, of the theoretical prediction based on λ+and D. Kee=1 and Kei=3.

Figure 4:

Cyclic inhibitory Ki modulation. The upper panel shows the time course of Ki, and the middle and lower panels show Rexc and Rinh, respectively. The blue, green, and yellow background shadows represent synchronized, bistable, and desynchronized states, respectively, of the theoretical prediction based on λ+and D. Kee=1 and Kei=3.

Close modal

This study demonstrates the various dynamics of the EI-Kuramoto model based on the interaction strength K. Furthermore, by modulating the inhibitory interaction strength, the shift in dynamic states can be modeled, potentially mirroring the biological regulatory mechanisms observed in neural collective oscillations.

The spectrum of dynamics governing collective neuronal oscillations in the brain likely mirrors the continuum between desynchronized and synchronized states observed in our model. A desynchronized state, especially under strong inhibition, is advantageous for the precise processing of external information in the brain (Harris & Thiele, 2011; Nobre & Kastner, 2014). As inhibition gradually decreases toward synchronization, resembling a bistable state in our model, it may create conditions that enable the simultaneous activation of related information, facilitating associative memory and inference. Further reduction of inhibition leads to increased synchronization, similar to conditions observed in sleep, epilepsy (Shao et al., 2019), or pathological symptoms (Jahangir et al., 2021), potentially hindering the processing of external information. The behavior exhibited by our model, where excUnits attract both themselves and inhUnits to the same phase (see Figures 2b2f), may reflect the correlated excitatory and inhibitory neuronal activities observed in the biological brain (Atallah & Scanziani, 2009). Furthermore, numerous neuronal computational models highlight the role of inhibition in synchronizing neurons (Börgers & Kopell, 2003; Brunel & Wang, 2003; Chalk et al., 2016; Hansel & Mato, 2003; Liang et al., 2020; Tetzlaff et al., 2012; Wilson & Cowan, 1972). These parallels between our model and brain activities provide a promising avenue for deepening our understanding of neuronal oscillations and their cognitive implications.

The Kee value, which represents excitatory-excitatory synaptic connections, must be smaller than the other K values to induce bistable or desynchronized states. The combined action of Kie and Kei, which provides feedback repulsion to the excUnits, plays a crucial role in driving the bistable and desynchronized states. In particular, Kii facilitates desynchronization. By using λ+ and D, both derived under the assumption of N, we could approximately predict the dynamic states for the finite N system (see Figure 4), indicating consistency between the theoretical analysis for the infinite N case and the numerical simulation for the finite N case. However, slight discrepancies were observed near the transition boundaries (see Figure 4), and it remains unclear whether these differences result from finite-N effects or from the time required for the transitions. Overall, the dynamics of the EI-Kuramoto model are primarily influenced by the balance of K values in both finite and infinite systems.

While our model shares the fundamental premise with the one proposed by Montbrió and Pazó (2018), in considering a heterogeneous population of excitatory and inhibitory oscillators, it diverges in several important ways. Montbrió and Pazó focused on synchronization driven solely by feedback mechanisms (Kie and Kei in this study) without direct attractive interactions (Kee in this study). In contrast, we conduct a comprehensive parameter exploration, independently varying all four types of interaction strengths to examine how their balance influences transitions between synchronization patterns.

Additionally, our parameter exploration differs in scope. While Montbrió and Pazó combined certain interaction strengths into composite parameters (e.g., K and εK in Montbrió & Pazó, 2018) and allowed for a natural frequency difference (Δω) between excitatory and inhibitory populations, we set the mean natural frequencies fixed at zero and systematically modulate each interaction strength independently.

Furthermore, while Montbrió and Pazó (2018) employed cosine functions for their interaction terms, we chose sine functions. The choice of interaction form is crucial and can significantly alter the resulting synchronization patterns (Omel'chenko & Wolfrum, 2012; Sakaguchi & Kuramoto, 1986). The cosine form, introducing a phase lag of −π/2 relative to the sine form, does not promote synchronization in the Kuramoto model (Clusella et al., 2022; Montbrió & Pazó, 2018; Sakaguchi & Kuramoto, 1986). Interactions between actual neurons can be assessed using phase-response curves (PRCs), which characterize how neuronal phase shifts depend on perturbations such as synaptic input at different points in the oscillatory cycle. PRCs are typically classified into two major types (Breakspear et al., 2010; Smeal et al., 2010). Cosine-wave interactions can approximate type 1 PRCs under specific conditions (Clusella et al., 2022; Montbrió & Pazó, 2018), while sine wave interactions more closely resemble type 2 PRCs (Breakspear et al., 2010; Smeal et al., 2010). Overall, Montbrió and Pazó (2018) derived their cosine interaction from a neural model, while we selected sine functions because they offer a more intuitive representation of basic attractive or repulsive forces, providing a clearer interpretation of the underlying mechanisms and helping to elucidate the relationship between the balance of interaction strengths and the synchronization phenomena in the system. Through these differences from Montbrió and Pazó (2018), our findings and theirs complement each other, enhancing our understanding of the EI-Kuramoto model.

Several models of coupled oscillators, including our own, have demonstrated bistable or desynchronized states. Features such as subgrouping (Abrams et al., 2008; Montbrió et al., 2004; Okuda & Kuramoto, 1991; Petkoski et al., 2016), second-order interactions (Ji et al., 2014; Kori & Kuramoto, 2001; Xu & Skardal, 2021), or interaction delays (Abrams et al., 2008; Montbrió et al., 2004; Petkoski et al., 2016; Yeung & Strogatz, 1999) can induce fluctuations in the order parameter of the Kuramoto model. Additionally, coupled oscillators with repulsive interactions have been shown to undergo desynchronization (Ryu et al., 2021; Salova & D'Souza, 2020). However, many of these studies focus on populations subject to same-directional interactions (excitatory or inhibitory interactions in our study) under various parameter settings, with relatively few exploring how opposite-directional interactions (excitatory and inhibitory interactions in our study) influence dynamic states. In this regard, Montbrió and Pazó (2018) and our work are notable exceptions. Investigating whether our findings are qualitatively similar to or differ from these earlier efforts will help clarify how collective oscillator dynamics relate to neuronal synchronization.

The EI-Kuramoto model in this study provides valuable insights but only partially captures the complexity of neuronal activity. In particular, synaptic connections in the brain are sparse, and their strengths follow highly skewed distributions (Buzsáki & Mizuseki, 2014; Iwase et al., 2024; Mizuseki & Buzsáki, 2013), unlike the uniform application of each interaction strength in our study. The brain network exhibits scale-free dynamics, which play a significant role in synchronization processes (Ódor & Kelling, 2019). Neurons, particularly inhibitory ones, have diverse connections and firing patterns (Zeng & Sanes, 2017). Excitatory and inhibitory inputs from neuronal spikes do not function as purely attractive or repulsive forces, as their effects depend on the phase relationships between presynaptic and postsynaptic neurons during spike events (Breakspear et al., 2010; Izhikevich, 2010; Smeal et al., 2010). In the future, we should refine the model to incorporate these biological features and validate it by examining brain dynamics.

Even during sleep, synaptic connections undergo plastic changes that affect learning and memory consolidation (Goto et al., 2021; Norimoto et al., 2018). Different resting states, such as resting awake, rapid eye movement (REM) sleep, or non-REM sleep, may influence the plasticity of neural circuits in distinct ways (Landmann et al., 2014; Lewis et al., 2018). Investigating modified EI-Kuramoto models that incorporate plasticity (Aoki & Aoyagi, 2009; Bronski et al., 2017; Fialkowski et al., 2023; Ha et al., 2016, 2018; Niyogi & English, 2009; Park et al., 2021; Picallo & Riecke, 2011; Poyato, 2019; Ren & Zhao, 2007; Seliger et al., 2002) could enhance our understanding of how synchronized neural activity influences learning-related processes.

The scripts for this study can be found in the Google Colaboratory: https://colab.research.google.com/drive/1HT7IvHWl1oEiXh6Ka58QvrcZ5kOfQlDE?usp=sharing

We thank Takuya Isomura and Hayato Chiba for their valuable advice on the theoretical analysis. Our gratitude also goes to Hiroyuki Miyawaki and Susumu Setogawa for their insightful contributions to the discussion. This study was supported by JSPS KAKENHI (grant 24K18243 to S.K.; grant 23H02788 to K.M.) and the Takeda Science Foundation (to K.M.). The funding sources had no involvement in the study design, data collection, interpretation, or the decision to submit this letter for publication.

Abrams
,
D. M.
,
Mirollo
,
R.
,
Strogatz
,
S. H.
, &
Wiley
,
D. A.
(
2008
).
Solvable model for chimera states of coupled oscillators
.
Physical Review Letters
,
101
(
8
), 084103.
Acebrón
,
J. A.
,
Bonilla
,
L. L.
,
Pérez Vicente
,
C. J.
,
Ritort
,
F.
, &
Spigler
,
R.
(
2005
).
The Kuramoto model: A simple paradigm for synchronization phenomena
.
Reviews of Modern Physics
,
77
(
1
),
137
185
.
Alfonsa
,
H.
,
Burman
,
R. J.
,
Brodersen
,
P. J. N.
,
Newey
,
S. E.
,
Mahfooz
,
K.
,
Yamagata
,
T.
, . . .
Akerman
,
C. J.
(
2022
).
Intracellular chloride regulation mediates local sleep pressure in the cortex
.
Nature Neuroscience
,
26
(
1
),
64
78
.
Andersen
,
P.
,
Morris
,
R.
,
Amaral
,
D.
,
Bliss
,
T.
, &
O'Keefe
,
J.
(Eds.)
. (
2006
).
The hippocampus book
.
Oxford University Press
.
Aoki
,
T.
, &
Aoyagi
,
T.
(
2009
).
Co-evolution of phases and connection strengths in a network of phase oscillators
.
Physical Review Letters
,
102
(
3
), 034101.
Atallah
,
B. V.
, &
Scanziani
,
M.
(
2009
).
Instantaneous modulation of gamma oscillation frequency by balancing excitation with inhibition
.
Neuron
,
62
(
4
),
566
577
.
Bick
,
C.
,
Goodfellow
,
M.
,
Laing
,
C. R.
, &
Martens
,
E. A.
(
2020
).
Understanding the dynamics of biological and neural oscillator networks through exact mean-field reductions: A review
.
Journal of Mathematical Neuroscience
,
10
(
1
), 9.
Börgers
,
C.
, &
Kopell
,
N.
(
2003
).
Synchronization in networks of excitatory and inhibitory neurons with sparse, random connectivity
.
Neural Computation
,
15
(
3
),
509
538
.
Braitenberg
,
V.
, &
Schüz
,
A.
(
1998
).
Cortex: Statistics and geometry of neuronal connectivity
.
Springer
.
Breakspear
,
M.
,
Heitmann
,
S.
, &
Daffertshofer
,
A.
(
2010
).
Generative models of cortical oscillations: Neurobiological implications of the Kuramoto model
.
Frontiers in Human Neuroscience
,
4
, 190.
Bronski
,
J. C.
,
He
,
Y.
,
Li
,
X.
,
Liu
,
Y.
,
Sponseller
,
D. R.
, &
Wolbert
,
S.
(
2017
).
The stability of fixed points for a Kuramoto model with Hebbian interactions
.
Chaos
,
27
(
5
), 053110.
Brunel
,
N.
, &
Wang
,
X. J.
(
2003
).
What determines the frequency of fast network oscillations with irregular neural discharges? I. Synaptic dynamics and excitation-inhibition balance
.
Journal of Neurophysiology
,
90
(
1
),
415
430
.
Buzsáki
,
G.
, &
Mizuseki
,
K.
(
2014
).
The log-dynamic brain: How skewed distributions affect network operations
.
Nature Reviews Neuroscience
,
15
(
4
),
264
278
.
Campbell
,
M. G.
,
Ocko
,
S. A.
,
Mallory
,
C. S.
,
Low
,
I. I. C.
,
Ganguli
,
S.
, &
Giocomo
,
L. M.
(
2018
).
Principles governing the integration of landmark and self-motion cues in entorhinal cortical codes for navigation
.
Nature Neuroscience
,
21
(
8
),
1096
1106
.
Chalk
,
M.
,
Gutkin
,
B.
, &
Denève
,
S.
(
2016
).
Neural oscillations as a signature of efficient coding in the presence of synaptic delays
.
eLife
,
5
, e13824.
Clusella
,
P.
,
Pietras
,
B.
, &
Montbrió
,
E.
(
2022
).
Kuramoto model for populations of quadratic integrate-and-fire neurons with chemical and electrical coupling
.
Chaos
,
32
(
1
), 013105.
Fialkowski
,
J.
,
Yanchuk
,
S.
,
Sokolov
,
I. M.
,
Schöll
,
E.
,
Gottwald
,
G. A.
, &
Berner
,
R.
(
2023
).
Heterogeneous nucleation in finite-size adaptive dynamical networks
.
Physical Review Letters
,
130
(
6
), 067402.
Goto
,
A.
,
Bota
,
A.
,
Miya
,
K.
,
Wang
,
J.
,
Tsukamoto
,
S.
,
Jiang
,
X.
, . . . &
Hayashi
,
Y.
(
2021
).
Stepwise synaptic plasticity events drive the early phase of memory consolidation
.
Science
,
374
(
6569
),
857
863
.
Ha
,
S. Y.
,
Lee
,
J.
,
Li
,
Z.
, &
Park
,
J.
(
2018
).
Emergent dynamics of Kuramoto oscillators with adaptive couplings: Conservation law and fast learning
.
SIAM Journal on Applied Dynamical Systems
,
17
(
2
),
1560
1588
.
Ha
,
S. Y.
,
Noh
,
S. E.
, &
Park
,
J.
(
2016
).
Synchronization of Kuramoto oscillators with adaptive couplings
.
SIAM Journal on Applied Dynamical Systems
,
15
(
1
),
162
194
.
Haider
,
B.
,
Häusser
,
M.
, &
Carandini
,
M.
(
2013
).
Inhibition dominates sensory responses in the awake cortex
.
Nature
,
493
(
7430
),
97
100
.
Hansel
,
D.
, &
Mato
,
G.
(
2003
).
Asynchronous states and the emergence of synchrony in large networks of interacting excitatory and inhibitory neurons
.
Neural Computation
,
15
(
1
),
1
56
.
Harris
,
K. D.
, &
Thiele
,
A.
(
2011
).
Cortical state and attention
.
Nature Reviews Neuroscience
,
12
(
9
),
509
523
.
Iwase
,
M.
,
Diba
,
K.
,
Pastalkova
,
E.
, &
Mizuseki
,
K.
(
2024
).
Dynamics of spike transmission and suppression between principal cells and interneurons in the hippocampus and entorhinal cortex
.
Hippocampus
,
34
(
8
),
393
421
.
Izhikevich
,
E. M.
(
2010
).
Dynamical systems in neuroscience: The geometry of excitability and bursting
.
MIT Press
.
Jahangir
,
M.
,
Zhou
,
J.-S.
,
Lang
,
B.
, &
Wang
,
X. P.
(
2021
).
GABAergic system dysfunction and challenges in schizophrenia research
.
Frontiers in Cell and Developmental Biology
,
9
, 663854.
Ji
,
P.
,
Peron
,
T. K. D. M.
,
Rodrigues
,
F. A.
, &
Kurths
,
J.
(
2014
).
Low-dimensional behavior of Kuramoto model with inertia in complex networks
.
Scientific Reports
,
4
, 4783.
Kori
,
H.
, &
Kuramoto
,
Y.
(
2001
).
Slow switching in globally coupled oscillators: Robustness and occurrence through delayed coupling
.
Physical Review. E, Statistical, Nonlinear, and Soft Matter Physics
,
63
(
4 Pt. 2
), 046214.
Kuramoto
,
Y.
(
1984
).
Chemical oscillations, waves, and turbulence
.
Springer
.
Kuramoto
,
Y.
(
1975
).
Self-entrainment of a population of coupled non-linear oscillators
. In
International Symposium on Mathematical Problems in Theoretical Physics
(pp.
420
422
).
Landmann
,
N.
,
Kuhn
,
M.
,
Piosczyk
,
H.
,
Feige
,
B.
,
Baglioni
,
C.
,
Spiegelhalder
,
K.
, . . .
Nissen
,
C.
(
2014
).
The reorganisation of memory during sleep
.
Sleep Medicine Reviews
,
18
(
6
),
531
541
.
Lewis
,
P. A.
,
Knoblich
,
G.
, &
Poe
,
G.
(
2018
).
How memory replay in sleep boosts creative problem-solving
.
Trends in Cognitive Sciences
,
22
(
6
),
491
503
.
Liang
,
J.
,
Zhou
,
T.
, &
Zhou
,
C.
(
2020
).
Hopf bifurcation in mean field explains critical avalanches in excitation-inhibition balanced neuronal networks: A mechanism for multiscale variability
.
Frontiers in Systems Neuroscience
,
14
, 580011.
Marder
,
E.
, &
Bucher
,
D.
(
2001
).
Central pattern generators and the control of rhythmic movements
.
Current Biology
11
(
23
),
R986
R996
.
Miyawaki
,
H.
, &
Diba
,
K.
(
2016
).
Regulation of hippocampal firing by network oscillations during sleep
.
Current Biology
,
26
(
7
),
893
902
.
Mizuseki
,
K.
, &
Buzsáki
,
G.
(
2013
).
Preconfigured, skewed distribution of firing rates in the hippocampus and entorhinal cortex
.
Cell Reports
,
4
(
5
),
1010
1021
.
Montbrió
,
E.
,
Kurths
,
J.
, &
Blasius
,
B.
(
2004
).
Synchronization of two interacting populations of oscillators
.
Physical Review. E, Statistical, Nonlinear, and Soft Matter Physics
,
70
(
5 Pt. 2
), 056125.
Montbrió
,
E.
, &
Pazó
,
D.
(
2018
).
Kuramoto model for excitation-inhibition-based oscillations
.
Physical Review Letters
,
120
(
24
), 244101.
Niyogi
,
R. K.
, &
English
,
L. Q.
(
2009
).
Learning-rate-dependent clustering and self-development in a network of coupled phase oscillators
.
Physical Review. E, Statistical, Nonlinear, and Soft Matter Physics
,
80
(
6 Pt. 2
), 066213.
Nobre
,
A. C.
, &
Kastner
,
S.
(Eds.)
. (
2014
).
The Oxford handbook of attention
.
Oxford University Press
.
Norimoto
,
H.
,
Makino
,
K.
,
Gao
,
M.
,
Shikano
,
Y.
,
Okamoto
,
K.
,
Ishikawa
,
T.
, . . .
Ikegaya
,
Y.
(
2018
).
Hippocampal ripples down-regulate synapses
.
Science
,
359
(
6383
),
1524
1527
.
Ódor
,
G.
, &
Kelling
,
J.
(
2019
).
Critical synchronization dynamics of the Kuramoto model on connectome and small world graphs
.
Scientific Reports
,
9
(
1
), 19621.
Okuda
,
K.
, &
Kuramoto
,
Y.
(
1991
).
Mutual entrainment between populations of coupled oscillators
.
Progress of Theoretical Physics
,
86
(
6
),
1159
1176
.
Omel'chenko
,
O. E.
, &
Wolfrum
,
M.
(
2012
).
Nonuniversal transitions to synchrony in the Sakaguchi-Kuramoto model
.
Physical Review Letters
,
109
(
16
), 164101.
Park
,
J.
,
Poyato
,
D.
, &
Soler
,
J.
(
2021
).
Filippov trajectories and clustering in the Kuramoto model with singular couplings
.
Journal of the European Mathematical Society
,
23
(
10
),
3193
3278
.
Petkoski
,
S.
,
Spiegler
,
A.
,
Proix
,
T.
,
Aram
,
P.
,
Temprado
,
J. J.
, &
Jirsa
,
V. K.
(
2016
).
Heterogeneity of time delays determines synchronization of coupled oscillators
.
Physical Review. E
,
94
(
1–1
), 012209.
Picallo
,
C. B.
, &
Riecke
,
H.
(
2011
).
Adaptive oscillator networks with conserved overall coupling: Sequential firing and near-synchronized states
.
Physical Review. E, Statistical, Nonlinear, and Soft Matter Physics
,
83
(
3 Pt. 2
), 036206.
Poyato
,
D.
(
2019
).
Filippov flows and mean-field limits in the kinetic singular Kuramoto model
.
arXiv [math.AP]. arXiv
.
Ren
,
Q.
, &
Zhao
,
J.
(
2007
).
Adaptive coupling and enhanced synchronization in coupled phase oscillators
.
Physical Review. E, Statistical, Nonlinear, and Soft Matter Physics
,
76
(
1 Pt. 2
), 016207.
Ryu
,
H.
,
Miller
,
J.
,
Teymuroglu
,
Z.
,
Wang
,
X.
,
Booth
,
V.
, &
Campbell
,
S. A.
(
2021
).
Spatially localized cluster solutions in inhibitory neural networks
.
Mathematical Biosciences
,
336
, 108591.
Sakaguchi
,
H.
, &
Kuramoto
,
Y.
(
1986
).
A soluble active rotater model showing phase transitions via mutual entertainment
.
Progress of Theoretical Physics
,
76
(
3
),
576
581
.
Salova
,
A.
, &
D'Souza
,
R. M
. (
2020
).
Decoupled synchronized states in networks of linearly coupled limit cycle oscillators
.
Physical Review Research
,
2
(
4
), 043261.
Seliger
,
P.
,
Young
,
S. C.
, &
Tsimring
,
L. S.
(
2002
).
Plasticity and learning in a network of coupled phase oscillators
.
Physical Review. E, Statistical, Nonlinear, and Soft Matter Physics
,
65
(
4 Pt. 1
), 041906.
Shao
,
L. R.
,
Habela
,
C. W.
, &
Stafstrom
,
C. E.
(
2019
).
Pediatric epilepsy mechanisms: Expanding the paradigm of excitation/inhibition imbalance
.
Children
,
6
(
2
), 23.
Smeal
,
R. M.
,
Ermentrout
,
G. B.
, &
White
,
J. A.
(
2010
).
Phase-response curves and synchronized neural networks
.
Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences
,
365
(
1551
),
2407
2422
.
Strogatz
,
S. H.
, &
Mirollo
,
R. E.
(
1991
).
Stability of incoherence in a population of coupled oscillators
.
Journal of Statistical Physics
,
63
(
3
),
613
635
.
Tetzlaff
,
T.
,
Helias
,
M.
,
Einevoll
,
G. T.
, &
Diesmann
,
M.
(
2012
).
Decorrelation of neural-network activity by inhibitory feedback
.
PLOS Computational Biology
,
8
(
8
), e1002596.
Tiesinga
,
P.
, &
Sejnowski
,
T. J.
(
2009
).
Cortical enlightenment: Are attentional gamma oscillations driven by ING or PING?
Neuron
,
63
(
6
),
727
732
.
Wilson
,
H. R.
, &
Cowan
,
J. D.
(
1972
).
Excitatory and inhibitory interactions in localized populations of model neurons
.
Biophysical Journal
,
12
(
1
),
1
24
.
Xu
,
C.
, &
Skardal
,
P. S.
(
2021
).
Spectrum of extensive multiclusters in the Kuramoto model with higher-order interactions
.
Physical Review Research
,
3
(
1
), 013013.
Yeung
,
M. K. S.
, &
Strogatz
,
S. H.
(
1999
).
Time delay in the Kuramoto model of coupled oscillators
.
Physical Review Letters
,
82
(
3
),
648
651
.
Zeng
,
H.
, &
Sanes
,
J. R.
(
2017
).
Neuronal cell-type classification: Challenges, opportunities and the path forward
.
Nature Reviews. Neuroscience
,
18
(
9
),
530
546
.
Zhang
,
H.
, &
Jacobs
,
J.
(
2015
).
Traveling theta waves in the human hippocampus
.
Journal of Neuroscience
,
35
(
36
),
12477
12487
.
Zhang
,
H.
,
Watrous
,
A. J.
,
Patel
,
A.
, &
Jacobs
,
J.
(
2018
).
Theta and alpha oscillations are traveling waves in the human neocortex
.
Neuron
,
98
(
6
),
1269
1281
.
This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. For a full description of the license, please visit https://creativecommons.org/licenses/by/4.0/legalcode

Supplementary data