Introduction

Empirical evidence suggests that brain activity during sleep and rest promotes memory consolidation, knowledge abstraction, and creative insight ^1,2 through the reorganization of neural circuits ^3,4. Such reorganization should rely on mechanisms distinct from those involved in learning during task performance. Indeed, synaptic plasticity during sleep and rest improves task performance ^5,6.

Neuronal synchronization occurs in two distinct states. Collective neuronal firing is desynchronized during active wakefulness, including task engagement, movement, and attentional states, as well as during rapid eye movement (REM) sleep, but more synchronized during quiet wakefulness and non-REM (NREM) sleep. Throughout the paper, we refer to these regimes descriptively as wake-like/desynchronized and sleep/rest-like/synchronized or bistable states, depending on the model dynamics being discussed. The transition between these states is associated with shifts in the excitation–inhibition (EI) balance. Firing rates of inhibitory neurons and synaptic conductance of inhibitory synapses are significantly elevated during wake-like states, while those of excitatory synapses show comparatively smaller changes ^7–10. These observations raise the possibility that, during sleep/rest-like states, EI-balance-dependent changes in collective neuronal synchronization interact with synaptic plasticity to reshape network structure. However, detailed mechanisms by which neuronal synchronization in sleep/rest-like states reshapes network structure through plasticity remain unclear.

Collective synchronization of oscillators, including neurons, has been studied using coupled-oscillator models. Among these, the Kuramoto model is one of the simplest frameworks^11,12, and numerous variants have been proposed. One variant, the EI-Kuramoto model, incorporates the EI balance by dividing oscillators into excitatory and inhibitory groups and assigning four interaction types between them^13,14. We have previously shown that tuning these four strengths toggles the network between synchronized, desynchronized, and bistable states, potentially reflecting a neuronal synchronization mechanism¹⁴. That study, however, kept all coupling weights—defined in this paper as pair-specific interaction strengths—fixed and identical, and therefore could not address how EI-balance-dependent dynamical states shape synaptic coupling through plasticity over time.

In this study, we focused on Hebbian and homeostatic plasticity. Hebbian plasticity, based on the principle “fire together, wire together,” strengthens synapses between neurons that fire in close temporal proximity¹⁵. In contrast, homeostatic plasticity regulates neuronal excitability and synaptic strength to maintain network stability, thereby preventing excessive excitation or quiescence¹⁶. The interplay between these plasticity mechanisms enables the emergence of diverse neuronal functions and firing patterns^17–20. Previous studies have incorporated plasticity into Kuramoto models, enabling coupling weights to adapt based on the oscillator activity to examine their impact on system synchronization^{21– 31}. However, these studies focused on promoting synchronization under a single dynamic state. How distinct collective dynamics, controlled by the EI balance, differentially shape synaptic plasticity remains unexplored.

To address these gaps, we developed a plastic EI (pEI)-Kuramoto model in which excitatory–excitatory couplings undergo Hebbian potentiation and homeostatic regulation. We found that strong inhibition produced a desynchronized regime with relatively stable couplings, whereas weaker inhibition produced a bistable regime in which intermediate-strength couplings fluctuated while strong couplings remained stable. These results suggest that EI-balance-dependent collective dynamics selectively stabilize strong couplings while permitting flexible reorganization of weaker ones, providing a potential mechanism for network reorganization during sleep and rest.

Results

Establishment of the pEI-Kuramoto model

The pEI-Kuramoto model comprised two oscillator groups: excitatory units, referred to as excUnits, and inhibitory units, referred to as inhUnits. These specific terms clarify that the units represented excitatory and inhibitory factors in this model, rather than actual neurons or synapses. Four interaction strengths were assigned between the groups, with only excitatory-excitatory couplings (J_ee) being distributed and plastic. The other interaction strengths, K_ii (repulsive strength between inhUnits), K_ei (attractive strength from excUnits to inhUnits), and K_ie (repulsive strength from inhUnits to excUnits) were constant across time and across unit pairs (Figure 1A; see Materials and Methods). The interactions are modeled as zero-phase-lag sine functions. The interactions from excUnits are purely attractive, while those from inhUnits are purely repulsive. Accordingly, the terms “excitatory” and “inhibitory” are used synonymously with “attractive” and “repulsive,” respectively.

Figure 1.

First, we investigated how the distribution of pair-specific excitatory-excitatory coupling weights J_ee affected the dynamics of the pEI-Kuramoto model in the absence of plasticity (Figure 1B and C). In our previous study, we demonstrated that the Kuramoto model with EI balance (EI-Kuramoto model) exhibits three distinct dynamic states—synchronized, bistable, and desynchronized—depending on the balance between the constant values of the interaction strengths¹⁴. In this study, J_ee followed a Gaussian distribution with a mean of μ_Jee and a standard deviation of σ_Jee, where σ_Jee was set to 0, 1/10, and 1/2 of μ_Jee (Figure 1B). The collective synchronization was quantified using the order parameter R, which represents the magnitude of the average vector between the oscillator units. In this simulation, regardless of the standard deviation value, all systems with a distributed J_ee also exhibited synchronized (high R), bistable (alternating high and low R), and desynchronized (low R) dynamics, as observed with constant (σ_Jee = 0) values (Figure 1C).

Next, we incorporated plasticity into J_ee by applying Hebbian potentiation and homeostatic regularization plasticity rules (see Materials and Methods). Among the various methods available for implementing Hebbian plasticity in the Kuramoto model, we initially selected the simplest method: the cosine function (Figure S1A; see Materials and Methods)^{21–23,26,30,31}. In this process, when a pair of excUnits exhibits closer phases, their J_ee coupling weight is potentiated. Simultaneously, the J_ee distribution was rescaled to maintain the initial mean and standard deviation, thereby preventing overpotentiation through homeostatic regularization plasticity (see Materials and Methods). The system dynamics were evaluated using the parameter set (μ_Jee = 1, σ_Jee = 0.1, K_ei = 3, K_ie = 3, K_ii = 3), at which the system exhibited a desynchronized state in the absence of J_ee plasticity. Under the cosine rule, unit pairs with smaller phase differences acquire a larger J_ee weight. This rule causes the J_ee distribution, which was initially Gaussian, to become bimodal (Figure S1B). Plasticity also influences the distribution of oscillator units, resulting in bimodal distributions of both phase and phase difference (Figure S1C), consistent with previous studies that incorporated the Hebbian rule into the Kuramoto model^21,26. This bimodal distribution of the coupling weight deviates from the long-tailed skewed distributions typically observed in biological synaptic weights^10,32,33.

To address this discrepancy, we introduced an exponential-based Hebbian plasticity rule to adjust the decay rate (sparseness parameter γ) of the gain for the phase difference (Figure 1D). We applied this rule to a range of parameter sets and evaluated the resulting J_ee distribution using skewness, kurtosis, and dip values, which quantify the distribution bimodality³⁴. Although J_ee distributions initially followed a Gaussian pattern, they transitioned into various bimodal or skewed forms (Figure 1E and F). We also quantified the synchronization of the excUnits using the mean and standard deviation of the order parameter R_exc over time, denoted as E[R_exc]_t and S[R_exc]_t, respectively. Plasticity in J_ee affects synchronization (Figure 1G), resulting in uniform, bimodal, or unimodal phase and phase-difference distributions depending on the plasticity parameters (Figure S1D). Moreover, the boundary of the statistical features of the J_ee distribution (Figure 1E) appears to coincide with the synchronization boundary (Figure 1G). For further analysis, we chose a parameter set that resulted in a long-tailed skewed J_ee distribution (Figure 1F, top, magenta circle), which more closely reflects the synaptic weight distributions observed in the brain^10,32,33.

Coupling fluctuations increase in the bistable state

We next investigated how distinct dynamic states associated with different EI balances affect coupling weights. To this end, we altered the EI balance of the pEI-Kuramoto model by manipulating the interaction strength from inhUnits. We unified K_ie and K_ii into a single parameter, K_i, to simplify the simulation, following the approach used in our previous study¹⁴. We then examined how the J_ee fluctuations were associated with the K_i values, quantified as the standard deviation of J_ee over time (S[J_ee]_t; see Materials and Methods) (Figure 2A and B). We determined the system’s dynamic state by analyzing the time course of R_exc (Figure 2C and D), as previously described¹⁴.

Figure 2.

When K_i < 0.4, both the average (E{S[J_ee]_t}_n) and standard deviation (S{S[J_ee]_t}_n) of S[J_ee]_t across unit pairs were almost zero (Figure 2A and B), indicating no plastic changes in J_ee. In this K_i range, the system was in a synchronized state, characterized by high E[R_exc]_t and low S[R_exc]_t (Figure 2C and D). The coupling weights between individual pairs exhibited no fluctuation over time after reaching the synchronized state (Figure 2E and F, K_i: 0.1; Video S1).

As K_i increased, E{S[J_ee]_t}_n also increased, peaking at K_i = 0.7, then decreased, and eventually stabilized around K_i = 1.5 (Figure 2A). Meanwhile, the S{S[J_ee]_t}_n continued to increase (Figure 2B). For K_i values above 1.5, the system transitioned into a desynchronized state characterized by low E[R_exc]_t and S[R_exc]_t values (Figure 2C and D). Strong or weak couplings remained stable over time, whereas those of intermediate strength fluctuated (Figure 2E, K_i: 1.5; Video S2). The system generated transient clusters of units coupled with strong J_ee weights (Figure 2F, K_i: 1.5, the lower leftmost panel, t = 24). Over time, most units left these clusters, leaving only small, stable subsets (Figure 2F, K_i: 1.5, the lower panels, t = 29–50).

In the intermediate range (0.4 < K_i < 1.5), the system exhibited a bistable state with high S[R_exc]_t (Figure 2D). Notably, the S[J_ee]_t was significantly higher in this bistable state (Figure 2A). Similar to the desynchronized state (K_i > 1.5), the couplings with high or low weights were stable. However, the intermediate-strength couplings fluctuated more (Figure 2E and F, K_i: 0.8 and 1.0). Notably, the number of these intermediate fluctuating couplings was higher in the intermediate K_i range (Figure 2E and F, K_i: 0.8 and 1.0; Videos S3 and S4) than in the low and high K_i ranges. Thus, transitions in the dynamic state were accompanied by changes in the coupling weight fluctuations. Particularly, the bistable state tended to exhibit larger J_ee fluctuations.

Note that Hebbian plasticity affected the dynamic state of the pEI-Kuramoto model (Figure S2A), slightly shifting the bifurcation points between the synchronized and bistable states, and the bistable and desynchronized states toward larger K_i values. Owing to this effect, we observed a small discrepancy from the state predicted by the eigenvalue-based classification for the EI-Kuramoto model¹⁴ (data not shown). However, the overall trend of the state transition from synchronized to bistable and then to desynchronized as K_i values increased remained consistent across different Hebbian learning rates. The enhancement of synchronization induced by Hebbian plasticity was consistent with the findings of previous reports on Kuramoto models incorporating cosine-shaped Hebbian plasticity^26,31. In contrast, σ_Jee did not affect the dynamic states (Figure S2B), consistent with the results in the absence of J_ee plasticity shown in Figure 1C.

Intermediate couplings are selectively destabilized while the strong couplings remain stable

We further investigated coupling weight fluctuations under conditions of slightly weaker inhibition than in the desynchronized state, assuming intermediate brain states between wake-like and sleep/rest-like states. We analyzed the individual J_ee coupling weights at selected points of K_i (K_i = 0.8, 1.0, and 1.5, red crosses in Figure 2A–D, except K_i = 0.1). Histograms of E[J_ee]_t (mean of J_ee over time for each pair of units) and S[J_ee]_t also tended to shift toward higher values as K_i decreased (Figure 3A and B). A scatter plot of E[J_ee]_t and S[J_ee]_t formed a characteristic “bow-shaped” distribution (Figure 3C). In this two-dimensional distribution, data points near the minimum and maximum E[J_ee]_t values had lower S[J_ee]_t values, whereas intermediate E[J_ee]_t showed a higher S[J_ee]_t. This pattern was consistent across all distributions for different K_i values. A similar phenomenon, the stabilization of strong couplings (e.g., synapses with large mushroom-shaped spines), has been observed in brain synapses^35–37.

Figure 3.

Furthermore, the S[J_ee]_t values tended to increase at intermediate E[J_ee]_t values with decreasing K_i (Figure 3C), suggesting that K_i can regulate the fluctuations in individual J_ee plastic couplings. At high K_i values in desynchronized states (K_i = 1.5, green dots), the couplings were relatively stable, whereas at low K_i values in bistable states (K_i = 1.0 and 0.8, orange and blue dots), the couplings of intermediate weights were more unstable. Strong J_ee coupling remained stable with small fluctuations, even when K_i was low. These dynamics may provide a mechanism that preserves the coupling of particularly strong pairs and resets moderate-weight pairs, facilitating the reorganization of the network structure (Figure 2E and 2F).

To enable a direct comparison with Figure 7A and 7B in Yasumatsu et al³⁷, which report the relationship between spine volume and changes in spine volume in hippocampal slices, we quantified the change in coupling weight, ΔJ_ee (= J_ee(t+Δt) – J_ee(t)), as a function of its current value, J_ee, at each time step in our pEI-Kuramoto model (Figure 3D and E). We pooled J_ee values and their changes ΔJ_ee across all excUnit pairs and time points, binned the J_ee values, and computed the mean and the standard deviation of ΔJ_ee for each J_ee bin (E[ΔJ_ee]_t,n and S[ΔJ_ee]_t,n). In the desynchronized state (K_i = 1.5, green lines), E[ΔJ_ee]_t,n was slightly but significantly positive for small J_ee (Wilcoxon signed-rank test, bin=[0.011,0.012), mean=1.1×10^-8, n=6,686,034, W=1.06×10¹³, p<1.0×10^-308), then decreased with increasing J_ee, reaching a local minimum near J_ee≈0.02. Subsequently, E[ΔJ_ee]_t,n increased and approached zero around J_ee≈0.025 (Figure 3D), a region corresponding to the peak of E[J_ee]_t (Figure 3A and C). S[ΔJ_ee]_t,n increased monotonically over the same range of J_ee (Figure 3E). These trends closely mirror observations in hippocampal slices and in vivo cortical imaging^37,38. In comparison, the bistable regimes (K_i=0.8 and 1.0, orange and blue lines) showed large deviations in the mean change (E[ΔJ_ee]_t,n) either in the same direction (K_i=1.0, orange line) or the opposite direction (K_i=0.8, blue line) as the desynchronized state (Figure 3D). Furthermore, the standard deviation (S[ΔJ_ee]_t,n) was generally elevated relative to the desynchronized state across the range of J_ee from 0.01 to 0.025 (Figure 3E).

Factors determining coupling weight properties

To identify the parameters that determine the individual J_ee coupling properties (E[J_ee]_t and S[J_ee]_t) shown in Figure 3, we explored the candidate parameters (natural frequency difference Δω, initial phase difference Δθ, and initial J_ee values of each excUnit pair) that may contribute to J_ee fluctuations (Figure 4). We colored the dots of E[J_ee]_t vs. S[J_ee]_t distribution (Figure 3C) according to each parameter value (Figure 4A–C). For statistical comparison, couplings were roughly categorized based on whether E[J_ee]_t exceeded 0.018 (high-E[J_ee]_t group) or not (low-E[J_ee]_t group), because the correlation between E[J_ee]_t and S[J_ee]_t reversed around this point (Figure 4D–F).

Figure 4.

Regarding differences in the natural frequency Δω, most pairs in the high (≥0.018) E[J_ee]_t group had relatively small Δω, whereas pairs with larger Δω tended to be in the low (<0.018) E[J_ee]_t group (Figure 4A and 4D). In the low-E[J_ee]_t group, S[J_ee]_t seemingly increased with Δω (Figure 4A). The initial phase differences Δθ between the pairs were also significantly different between the high- and low-E [J_ee]_t groups (Figure 4B and E). In the high-E[J_ee]_t group, most pairs exhibited small initial Δθ, whereas in the low-E[J_ee]_t group, the initial Δθ was uniformly distributed. Additionally, we examined the effect of the initial J_ee coupling weight, which showed no noticeable bias (Figure 4C) and followed a normal distribution in both the high and low-E[J_ee]_t groups (Figure 4F). Considering these results, pairs with small differences in natural frequency and initial phase positions maintained more stable and larger J_ee weights (Figure 4G and H). Conversely, the initial J_ee weight was not an important factor in determining the stability and strength of the couplings.

Theoretical analysis of quasi-steady coupling weight distribution

We investigated the theoretical mechanism underlying the fluctuation of coupling weights J_ee by analyzing their distribution in a quasi-steady state, assuming J_ee,kj(t + Δt) ≈ J_ee,kj(t) (Supplemental Text). Although a true static steady weight distribution does not exist in the desynchronized or bistable states due to the non-constant phases of each oscillator unit, the desynchronized condition allows for a quasi-steady state where the phase and weight distributions evolve on a sufficiently slow timescale, resulting in minimal fluctuation of the J_ee distribution. Under this assumption, the coupling weights do not change as a result of the plastic modulation in Equations M11 and M14, leading to the following relationship (Supplemental Text, Equation S3).

The first term, , represents the value of J_ee,kj, which is rescaled with the homeostatic regulation plasticity rule so that the mean and variance become and , respectively. The second term is the Hebbian potentiation component, scaled by the ratio between the Hebbian and homeostatic plasticity learning rates c. This equation indicates that the weight J_ee,kj depends on the distribution of and the distribution of phase difference between unit k and j Δθ_kj (c exp(−γ|Δθ_kj|)). Assuming that this rule remains applicable during the approach to a quasi-steady state, successive iterations of the update rule progressively integrate the second Δθ-dependent Hebbian term into the J_ee distribution. Furthermore, the distribution of Δθ is dependent on that of θ (Supplemental Text). Therefore, when the weight dynamics are sufficiently slow, the J_ee distribution is determined by the phase distribution. Indeed, the characteristics of the J_ee distribution appear to be correlated with the synchronization state of the system, which reflects the shape of the phase distribution (Figure 1E, 1F, and S1). Specifically, a uniform phase (and phase difference) distribution (desynchronized state) corresponds to a right-skewed, long-tailed J_ee distribution (magenta circle in Figure 1E, F and S1D). A bimodal phase distribution corresponds to a bimodal J_ee distribution (green circle in Figure 1E, F and S1D and the cos rule in S1B and S1C). A unimodal phase distribution (synchronized state) corresponds to a left-skewed J_ee distribution (yellow circle in Figure 1E, F and S1D).

Next, we considered the change in the phase difference for a pair of excUnits, Δθ_kj. Assuming that J_ee,kj and Δθ_kj are constant, the following relationship can be derived from Equations 1 and M9 (Supplemental Text, Equation S18):

Δω is the natural frequency difference between units k and j, F_inh is the interaction strength from the inhUnits, and is the interaction strength from the excUnits, except for units k and j. is the R_exc- and R_inh-dependent term; its direction can be either attractive or repulsive, depending on the phase configuration of the system (Supplemental Text). The right-hand side of this equation represents the interaction strength between units k and j, whereas the left-hand side represents the net forces of the other factors. A stable fixed point for Δθ_kj exists, where the functions for the left- and right-hand sides intersect (Figure S3).

On the right-hand side, the attractive force is driven by both the homeostatic set-point, , and the Hebbian potentiation term. Near Δθ≈0, the Hebbian effect is strong, yielding a large attractive force. As the phase difference ∣Δθ∣ increases, this Hebbian potentiation weakens, and the coupling weight relative position within the J_ee distribution falls to a lower value, causing an overall attractive force reduction.

On the left-hand side, while the influence of F_inh and is limited in a desynchronized state where R_exc and R_inh are near zero, their impact becomes non-negligible in the bistable state, where these order parameters undergo periodic modulation. Sufficiently large fluctuations in these forces cause the stable fixed point to disappear or be substantially displaced, thereby leading to large fluctuations in the phase difference between the unit pair.

Cyclic modulation of inhibition

Finally, we varied K_i to mimic the sleep–wake cycle^39,40. J_ee fluctuated (Figure 5A, Video S5) in response to K_i modulation (Figure 5B, top panel). R_exc alternated between the desynchronized and bistable states (Figure 5B, middle), with corresponding changes in E{S[J_ee]_t}_n (Figure 5B, bottom). When K_i was high, the system was in a desynchronized state, and E{S[J_ee]_t}_n was small. Conversely, when K_i was low, the system exhibited a bistable state, and E{S[J_ee]_t}_n increased. These dynamics can be considered as modeling how the circuit structure is stabilized during the wake-like state and destabilized during the sleep/rest-like state, thereby promoting circuit reorganization.

Figure 5.

Discussion

The main conceptual advance of this study is to show that EI balance can determine not only the synchronization state of a network but also whether plastic couplings are stabilized or rendered labile. In the model, strong inhibition favored a desynchronized regime with relatively stable couplings, whereas weaker inhibition produced a bistable regime that selectively destabilized intermediate couplings while preserving the strongest ones. This provides a candidate systems-level mechanism by which sleep/rest-like network states can promote reorganization without erasing established strong connections.

Although the sleep/rest-like state is often associated with synchrony, it is not uniformly synchronized. Near the transition between synchronized and desynchronized states, the dynamics tend to become irregular or weakly stable, as captured by two-population coupled-oscillator models ^41,42. Sleep onset can be formalized as a noisy saddle-node bifurcation ⁴³, and macroscopic brain-state transitions are often argued to occur near the bifurcation point ^44,45. Beyond the transitions between wakefulness and sleep, even within sleep, the microarchitecture exhibits intrinsic instability, such as the NREM-REM sleep cycle^3,4,39,40, the cyclic alternating pattern in NREM sleep⁴⁶, UP–DOWN state alternations^47,48, and microarousals⁴⁹, which may reflect switching among synchronized and desynchronized states. These fluctuation-driven transitions align with the bistable regime of our model under reduced inhibition.

Preserving robust couplings while randomizing weaker ones simplifies the network, rendering it sparse. The plasticity parameters used here, which yield a long-tailed coupling-weight distribution (Figure 1E and F) similar to that observed in the brain^10,32,33 and mirror the EI balance during state transitions and in sleep/rest-like states^7–10, appear to facilitate network sparsification. The selective preservation of strong synaptic connections is also frequently observed in the brain^35–37 and is understood through molecular mechanisms, including the dynamics of NMDA receptors^5,50 and actin filaments⁵¹. Our results introduce a complementary, system-level mechanism for sustaining strong coupling, viewed through the dynamics of collective oscillators. Additionally, our results may be mechanistically linked to the systems-level observation that only a subset of neuronal populations remain stable while others drift^52–54. This process resembles pruning strategies used in artificial neural networks that facilitate knowledge generalization^55,56. Such pruning-like dynamics are particularly relevant during sleep/rest-like states when synaptic downregulation is dominant^5,50, as fluctuations in coupling strength can weaken intermediate connections and eventually reset them toward baseline.

Fluctuations in network connectivity can promote network reorganization. In our model, bistable states under low inhibition conditions, mirroring the EI balance during state transitions and sleep/rest-like states^7–10, promoted fluctuating coupling weights between excUnits. It has also been proposed that the brain enhances neuronal circuit reorganization during the sleep/rest-like state^3,4. Combined with mechanisms that evaluate and stabilize reorganized circuits⁵⁷—potentially including dopaminergic signaling during NREM sleep⁵⁸, which may provide a reinforcement-like signal and gate plastic changes—the sleep/rest-like state could, in principle, contribute to the generation of novel associations, creative ideas, or innovative problem-solving strategies.

The characteristic “bow-shaped” relationship between coupling weight mean (E[J_ee]_t) and standard deviation (S[J_ee]_t) (Figure 3C) arises from the stability of phase differences between excUnits (Equation 2). Pairs with strong coupling maintain small and stable phase differences (Δθ≈0), ensuring that the Hebbian potentiation term (c exp(−γ|Δθ|)) and their relative position in the J_ee distribution remain consistently high. In contrast, pairs with weak couplings evolve independently with fluctuations driven mainly by natural frequency differences Δω. The largest S[J_ee]_t was found in pairs with intermediate weights. In this regime, the intrinsic synchronizing force between the two units is closely counterbalanced by external driving influences, chiefly their natural frequency difference, Δω. Furthermore, Hebbian potentiation and homeostatic depression are also balanced. Small fluctuations in Δθ can frequently tip these balances and reverse the direction of the weight change, leading to large fluctuations in the coupling weight J_ee.

From a dynamics perspective, these mechanisms explain not only the characteristic bow-shaped distribution in Figure 3C but also a key feature of biological synapses: larger synapses are more stable than smaller ones^35–37. Moreover, particularly in the desynchronized state, our model reproduces the relationship between synaptic size and changes in synaptic size reported in experimental studies^37,38 (Figure 3D and E, green lines; see Results), further supporting its biological relevance.

The quasi-steady J_ee distribution is determined by the phase-difference distribution (Figure 1E, 1F, S1C and S1D), as it reflects a transformation of the Δθ distribution via the Hebbian function (Equation 1). In this study, the exponential rule c exp(−γ|Δθ_kj|) transforms the uniform phase-difference distribution into a hyperbolic weight distribution (Supplemental Text, Equation S8–12). Modifying the Hebbian function would allow the generation of alternative J_ee distributions, such as log-normal distributions commonly observed in the brain^10,32,33.

Equation 2 further predicts that unit pairs with small natural frequency differences (Δω) and small initial phase differences (Δθ) are more robust and more likely to return to a stable, synchronized state even after large perturbations, consistent with the results in Figure 4. In our simulation (Figure 1F, magenta circle), the Hebbian plasticity-induced enhancement of interaction strength near Δθ=0 is limited, even though the plasticity slightly strengthens the restoring force near the stable point (Figure S3). Thus, the tendency for pairs with small Δω and small initial Δθ to be more robust and stable arises primarily from the underlying sinusoidal form of the baseline interaction function rather than plasticity-induced modifications. This finding resonates with experimental observations that neurons exhibiting synchronized firing prior to learning (analogous to a small Δθ) are preferentially incorporated into the same assembly during learning^59–61.

Fluctuations in J_ee at intermediate K_i levels (Figs. 2–4, K_i=0.8, 1.0) arise from bistability between synchronized and desynchronized states. This alternation perturbs phase distributions, which directly drive changes in J_ee. Concurrent periodic fluctuations of the order parameters R_exc and R_inh constantly modulate the interaction term, (Equation 2), further destabilizing phase differences Δθ of moderately coupled pairs.

While the Kuramoto model provides conceptual insight, it is substantially simpler than biological neural circuits. The model reduces neuronal interactions to attractive and repulsive forces, neglecting pulse-based communication, conduction delays, and complex network topology⁶². Additionally, our model simplifies the complexity of neuroplasticity—excluding neurotransmitters⁶³, inhibitory plasticity^64,65, and adaptation⁶⁶—and does not distinguish between specific sleep stages (e.g., NREM vs. REM), which are thought to play distinct roles in memory processing and circuit reorganization^3,4. Addressing these limitations is essential for a comprehensive understanding of neural dynamics.

Future work can extend the framework by introducing explicit evaluation and stabilization signals that select among reorganized circuits^57,67, potentially linking fluctuation-driven reorganization to selective retention during sleep and rest. Another important direction is to examine how synchrony transitions influence interregional plasticity between functionally distinct networks, such as hippocampus and cortex⁶⁸, thereby linking circuit-level reorganization to systems-level consolidation.

Materials and Methods

pEI-Kuramoto model

The original Kuramoto model^11,12,69 is expressed as:

where θ_j(t) is the phase of unit j at time t, ω_j is the natural frequency of unit j following a probability distribution g(ω), N is the number of units, and K is the interaction strength among the oscillation units. In the original Kuramoto model, K is applied constantly and uniformly to all unit pairs. The order parameter R and the mean phase ψ are then calculated from the mean vector of the oscillators as:

where i denotes an imaginary unit. Given Euler’s formula and Equation M2,

Considering Equation M3, Equation M1 becomes

In our previous EI-Kuramoto model¹⁴, the oscillation units were categorized into two groups: excitatory and inhibitory units, referred to as excUnits and inhUnits, respectively. Employing these terms clarified that the units in the model represented excitatory and inhibitory factors, distinguishing them from actual neurons or synapses. Zero-phase-lag sine interactions are purely attractive for excUnits and purely repulsive for inhUnits. Hence, the terms “excitatory” and “inhibitory” are used interchangeably with “attractive” and “repulsive,” respectively. The interaction strength K is divided into four quantifiable parameters: K_ee (attractive strength between excUnits), K_ii (repulsive strength between inhUnits), K_ei (attractive strength from excUnits to inhUnits), and K_ie (repulsive strength from inhUnits to excUnits). K_ee, K_ii, K_ei, and K_ie are distinct positive real numbers that are uniformly applicable across all the unit pairs. Consequently, Equation M4 becomes:

for excUnits, and

for inhUnits, where exc and inh represent the parameters for excUnits and inhUnits, respectively. The order parameter and mean phase of each excUnit and inhUnit are

and

In the present study, we developed a pEI-Kuramoto model (Figure 1A). A key modification of our previous EI-Kuramoto model is that the constant and uniform interaction strength value K_ee is replaced with the heterogeneous and plastic coupling weight (interaction strength of an individual unit pair) matrix J_ee. The J_ee distribution is defined by its mean μ_Jee and standard deviation σ_Jee, which are scaled with the inverse of N_exc to adjust for the system size and to facilitate comparisons with other Kuramoto-model variants. Accordingly, the initial J_ee distribution follows the Gaussian distribution , where . The other interaction strengths, K_ii, K_ei, and K_ie, were uniform and consistently applicable across all units. Equation M5 becomes

J_ee,kj is the coupling weight between the excUnits k (pre) and j (post).

Given the rotational symmetry, when the mean natural frequencies of excUnits and inhUnits are equal, , they can be set to 0 rad/s, which is equivalent to rotating the coordinate frame by the mean natural frequency . For all simulations in this study: N_exc = 80, N_inh = 20; total simulation duration T = 50 s; time step dt = 0.01 s. The g(ω) followed the Gaussian distribution g(ω)∼𝒩(0, 0.01) for both excUnits and inhUnits. In Figs. 2–5, K_ii and K_ie are unified with the same parameter K_i to simplify the simulation.

At each time step Δt (= dt = 0.01 s in this study), the J_ee values were first modulated by Hebbian potentiation and then rescaled using homeostatic regularization rules. The Hebbian potentiation plasticity is implemented as an exponential function:

α is the learning rate for Hebbian plasticity. Note that α should scale proportionally with an inverse of N_exc as α ∝ 1/N_exc to maintain consistency across different system sizes. γ determines the sharpness of the exponential decay from the peak. Δθ_kj is the phase difference between excUnit k and j

The implementation of Hebbian plasticity with a cosine function (Figure S1A) is defined as

In any Hebbian learning rule in this paper, weight changes of J_ee,kj and J_ee,kj are the same due to the symmetry with respect to the sign of Δθ_kj.

After the Hebbian modulations, J_Hebb was rescaled with the homeostatic regularization plasticity rule to prevent overpotentiation. Plasticity was implemented as follows:

represents the normalized value (homeostatic set-point) of J_Hebb under the initial setting of the mean and standard deviation . β is the learning rate for homeostatic plasticity. This procedure scales J_ee,kj toward the initial targets, and .

The change in coupling weight at each time step ΔJ_ee (Figure 3D and E) was defined as

For plotting, pairs of ΔJ_ee and J_ee(t) values were pooled across all excUnit pairs and time points, and the J_ee(t) values were binned (bin size = 0.001).

Cyclic inhibitory modulation

In Figure 5, K_i is modulated as a sine function

where A_s (= 0.32), ω_s (= 0.004π), and C_s (= 1.28) are the amplitude, frequency, and baseline offset for the sine modulation, respectively. The other parameters for interaction strengths are held constant to simplify the simulation: μ_Jee = 1, σ_Jee = 0.1, K_ei = 3, and T = 3,000 s.

Statistics

E[·]_t (also E{·}_t) and S[·]_t are the mean and the standard deviation over time, respectively, where t indicates the number of time points. E[·]_n and S[·]_n are the mean and the standard deviation among unit pairs, respectively, where n indicates the number of unit pairs. E[·]_t,n and S[·]_t,n represent mean and standard deviation, respectively, pooling all values across time and unit pairs. Here, · represents a certain variable. The range of time for the time series analyses (E[·]_t, S[·]_t, E[·]_t,n, and S[·]_t,n) is [T/2, T] to exclude the initial fluctuations. For Figure 5B lower panel, value of E{S[J_ee]_t}_n was estimated in the time window [-15s, +15 s] for each time point.

Descriptive statistics (mean, standard deviation, kurtosis, and skewness) and statistical tests were performed using the Python NumPy and SciPy libraries. The Dip test³⁴ was performed using the Python Diptest Library (https://github.com/RUrlus/diptest). The dip statistic quantifies the maximum distance between the observed distribution and the closest unimodal distribution. If the data distribution is close to unimodal distribution, the dip value is low. If the distribution is bimodal, the dip value is high. Hierarchical clustering was performed using the linkage function of the Python SciPy library.

Data availability

The current manuscript is a computational study, so no data have been generated for this manuscript. Modelling code is uploaded at https://colab.research.google.com/drive/1v8JF97G9YdwUOh_Wuh9i_TvuGa42fLAu?usp=sharing.

Acknowledgements

We would like to thank Susumu Setogawa for his insightful comments on the manuscript. This study was supported by Japan Society for the Promotion of Science KAKENHI [grant 24K18243 to S.K.; grants 23K27479 and 25H02508 to K.M.], Japan Agency for Medical Research and Development (grant JP24wm0625216 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 paper for publication.

Additional information

Author contributions

Conceptualization: SK; Methodology: SK; Investigation: SK; Visualization: SK; Supervision: KM; Writing—original draft: SK; Writing—review & editing: SK, KM; Funding acquisition: SK, KM

Funding

MEXT | Japan Society for the Promotion of Science (JSPS) (24K18243)

• Satoshi Kuroki

MEXT | Japan Society for the Promotion of Science (JSPS) (23K27479)

• Kenji Mizuseki

MEXT | Japan Society for the Promotion of Science (JSPS) (25H02508)

• Kenji Mizuseki

Japan Agency for Medical Research and Development (AMED) (JP24wm0625216)

• Kenji Mizuseki

Takeda Science Foundation (TSF)

• Kenji Mizuseki

Additional files

Supplemental video 1.

Supplemental video 2.

Supplemental video 3.

Supplemental video 4.

Supplemental video 5.

Supplemental information.