Abstract
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.
1 Introduction
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 Methods
2.1 EI-Kuramoto Model
Given the rotational symmetry, when the mean natural frequencies of excUnits and inhUnits are equal, , 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 and , 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 s for exploring the parameter (see section 3.1), s for cyclic inhibitory modulation (see section 3.3), time step s, and the of both excUnits and inhUnits populations following a gaussian distribution . The time range for calculating the mean values of the order parameters and was set from to to exclude initial fluctuations.
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. is the attractive interaction strength between excUnits (red arrows). is the attractive strength from excUnit to inhUnit (magenta arrow). is the repulsive strength between inhUnits (blue arrows). is the repulsive strength from inhUnit to excUnit (cyan arrow). (b) Phase of individual units (left) and -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 and in red and blue, respectively. The parameters are as follows: synchronized: , , , ; bistable: , , , ; desynchronized: , , , . (c) Time-averaged derived from the 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.
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. is the attractive interaction strength between excUnits (red arrows). is the attractive strength from excUnit to inhUnit (magenta arrow). is the repulsive strength between inhUnits (blue arrows). is the repulsive strength from inhUnit to excUnit (cyan arrow). (b) Phase of individual units (left) and -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 and in red and blue, respectively. The parameters are as follows: synchronized: , , , ; bistable: , , , ; desynchronized: , , , . (c) Time-averaged derived from the 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.
Detailed model dynamics in the bistable state. (a) Plots of (red) and (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 (-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: (black), (red), and (cyan). The right panels represent the total interaction strength acting toward inhUnits: (black), (magenta), and (blue). In panel g, the lower panels are not present, as in this case.
Detailed model dynamics in the bistable state. (a) Plots of (red) and (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 (-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: (black), (red), and (cyan). The right panels represent the total interaction strength acting toward inhUnits: (black), (magenta), and (blue). In panel g, the lower panels are not present, as in this case.
The equilibrium, where , is attractive if and repulsive if .
2.2 Cyclic Inhibitory Modulation
3 Results
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, ; excUnit-inhUnit attraction, ; inhUnit-excUnit repulsion, ; and inhUnit-inhUnit repulsion, . The 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 parameters (, , , and ) and observed the resulting model dynamics. Collective synchronization was assessed using the order parameter which quantifies the magnitude of the average vector across oscillator units.
Consequently, the dynamics were categorized into three distinct states: synchronized (high ), bistable (alternating high and low ), and desynchronized (low ) (see Figure 1b and supplementary videos 1–3). A systematic exploration of the parameter space for , , , and , 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 (see Figures 2b–2f). 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 ( for excUnits and for inhUnits; see section 2.1) and tracked their changes. If the slope of at (x-intercept) is negative, attracts toward the x-intercept. If the slope is positive, repels from the x-intercept.
The accumulation of excUnits at a specific phase position triggers a repetition of cycles, alternating between high and low 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 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 ( 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 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 is strong, the accumulation of the inhUnits (depicted in Figures 2b–2d) decreases due to mutual repulsion among the inhUnits. When the feedback repulsion to the excUnits through and 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 and hinders the accumulation of the excUnits, resulting in both the excUnits and inhUnits being desynchronized.
3.2 Bifurcation Analysis
The evolution of is negative in the desynchronized state. When has an imaginary number, the and oscillate in time.
The evolution equations for need not be written because they are simply conjugates of equations 3.17 and 3.18. The evolution equation for was zero at .
Note that we focus only on because ; thus, + , and .
Bifurcation of the dynamic states. (a) Bifurcation criteria. The dynamic state converges to a desynchronized state when (yellow rectangle), to a bistable state when and (green rectangle), and to a synchronized state when and (navy blue rectangle). The black lines in each quadrant represent the schematic time course of the order parameter . (b) Bifurcation of dynamic states based on the values of , , and .
Bifurcation of the dynamic states. (a) Bifurcation criteria. The dynamic state converges to a desynchronized state when (yellow rectangle), to a bistable state when and (green rectangle), and to a synchronized state when and (navy blue rectangle). The black lines in each quadrant represent the schematic time course of the order parameter . (b) Bifurcation of dynamic states based on the values of , , and .
3.3 Cyclic Inhibitory Modulation
Cyclic inhibitory modulation. The upper panel shows the time course of , and the middle and lower panels show and , respectively. The blue, green, and yellow background shadows represent synchronized, bistable, and desynchronized states, respectively, of the theoretical prediction based on and . and .
Cyclic inhibitory modulation. The upper panel shows the time course of , and the middle and lower panels show and , respectively. The blue, green, and yellow background shadows represent synchronized, bistable, and desynchronized states, respectively, of the theoretical prediction based on and . and .
4 Discussion
This study demonstrates the various dynamics of the EI-Kuramoto model based on the interaction strength . 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 2b–2f), 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 value, which represents excitatory-excitatory synaptic connections, must be smaller than the other values to induce bistable or desynchronized states. The combined action of and , which provides feedback repulsion to the excUnits, plays a crucial role in driving the bistable and desynchronized states. In particular, facilitates desynchronization. By using and , both derived under the assumption of , we could approximately predict the dynamic states for the finite system (see Figure 4), indicating consistency between the theoretical analysis for the infinite case and the numerical simulation for the finite case. However, slight discrepancies were observed near the transition boundaries (see Figure 4), and it remains unclear whether these differences result from finite- effects or from the time required for the transitions. Overall, the dynamics of the EI-Kuramoto model are primarily influenced by the balance of 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 ( and in this study) without direct attractive interactions ( 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., and 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.
Code Availability Statement
The scripts for this study can be found in the Google Colaboratory: https://colab.research.google.com/drive/1HT7IvHWl1oEiXh6Ka58QvrcZ5kOfQlDE?usp=sharing
Acknowledgments
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.