Olivocerebellar learning is highly adaptable, unfolding over minutes to weeks depending on the task. However, the stabilizing mechanisms of the synaptic dynamics necessary for ongoing learning remain unclear. We constructed a model to examine plasticity dynamics under stochastic input and investigate the impact of inferior olive (IO) reverberations on Purkinje cell (PCs) activity and synaptic plasticity. We explored Upbound and Downbound cerebellar micromodules, which are organized loops of IO neurons, cerebellar nuclei neurons and microzones of PCs characterized by their unique molecular profiles and different levels of baseline firing. Our findings show synaptic weight convergence followed by stability of synaptic weights. In line with their relatively low and high intrinsic firing, we observed that Upbound and Downbound PCs have a propensity for potentiation and depression, respectively, with both PC types reaching stability at differential levels of overall strength of their parallel-fiber (PF) inputs. The oscillations and coupling of IO neurons participating in the Upbound and Downbound modules determine at which frequency band PFs can be stabilized optimally. Our results indicate that specific frequency components drive IO resonance and synchronicity, which, in turn, regulate temporal patterning across Upbound and Downbound zones, orchestrating their plasticity dynamics.

The olivocerebellar system is a part of the brain that facilitates learning and controlling movements. To coordinate movements it integrates sensorimotor information with motor command signals. The resulting behaviour needs to be continuously adjusted during motor learning. The leading hypothesis is that changes in strength of the synaptic connections between neurons underlie the learning process. The challenge is to elucidate the factors that determine the exact timing and precision of the learned movements. To answer this question we developed a computational model of the two main types of modules of the olivocerebellar system that can control different types of movements in a bidirectional fashion. We found that the rhythm and coupling of the olivary neurons play an important role in controlling and stabilizing plasticity in the cerebellar cortex of both types of modules, together shaping learning-dependent timing of motor behaviour.

Designing an olivocerebellar loop model to examine plasticity dynamics

Our hypothesis is that PF-PC synaptic weights are stabilized during procedural learning by the resonant dynamics of the olivocerebellar loops. To test this, we developed a biologically plausible model of the olivocerebellar system (Fig 1A ), which is decoupled from the rest of the brain to focus on the olivocerebellar resonant network dynamics that drive plasticity. The questions we tested with this model are centered on exploration of synaptic stability during ongoing plasticity. We emulated ongoing arbitrary cortical and subcortical inputs to the cerebellum via 5 uncorrelated PF bundles, which is conveyed by 100 PCs via 40 CN cells to 40 cells in the IO [53–60]. These IO cells have the ability to respond by advancing or delaying their phase as a function of when inhibition arrives in the oscillation cycle, as observed in in vivo experiments [16,61–65].

As input to our model network we used bundles of PFs as a noisy current input rather than individual synapses. A bundle of PFs represents the combination of numerous synaptic sources, as expected during the awake and behaving state [66,67]. To avoid biasing the input, we modelled background activity as the sum of several independently generated Ornstein-Uhlenbeck (OU) processes, taken to represent the integrated, continuous current from large numbers of PF EPSCs [68]. This use of stochastic but temporally correlated input allowed us to assess whether synaptic weights stabilize under a fixed input structure, enabling us to isolate the contribution of internal loop dynamics, such as IO resonance and feedback, to plasticity. This choice was motivated by theoretical considerations and in vivo experimental data of granule cell and PF activity during learning, at rest and during walking [66,69,70], which shows complex, high dimensional firing patterns transmitted to PCs in awake animals. We refer to this temporally varied OU input as "frozen input" because it was generated once and then held constant across all simulations. The input was not repeated or manipulated across trials; it remained fixed throughout all simulations, serving as a controlled condition. This design allows us to compare network behavior across simulation conditions while ensuring that any observed differences could be attributed exclusively to the effects of plasticity mechanisms, rather than to variations in the external input. This is a crucial aspect of our model as we can test the hypothesis that under the stimulus representing sensorimotor contingencies there would be synaptic adjustment and weight stabilization in the PF bundles.Our PF bundles can be seen as bundles of correlated input originating via MF rosettes [10,71]. Given that around 90% of PF synapses are silent [72], we modelled 5 bundles with approximately 200 active PFs each (Fig 1A–C). The PF-PC synapse acts as a "meta-synapse", representing the average direction of change of the individual synaptic weights in the PF bundle. To represent plasticity, we added a weight value for each PF bundle, which could change independently. Model parameters were tuned to match spontaneous activity statistics reported experimentally — such as PC SSpk and CSpks firing rates, CN output, and IO STOs. No tuning was performed to enforce the hypothesized plasticity directions. Rather, the observed phenomena (e.g., synaptic weight dynamics, CF pause modulation, and frequency-dependent IO entrainment) emerged from the interactions of BCM and LTD plasticity, IO coupling, and intrinsic PC excitability under fixed input conditions.

We tuned the PCs firing frequency to represent a mixture of their intrinsic firing frequency and their response to the PF input. As shown experimentally [9], when the PF input was completely absent, PCs had an intrinsic firing frequency that was about 20 Hz lower than when the PF input was present. For simplicity we focussed on the primary drivers of the olivocerebellar loop dynamics - the GABAergic CN projection cells, and excluded cells involved in other loops. CN AdEx model cells were primarily designed to have robust transmission of rate coding [73]. Our CN cells received divergent input from the PC layer, which modulated CN responses in an anti-correlated fashion [73,74].

It is known that the IO contains a mixture of cells: some of which oscillate robustly, some display dampened oscillations, and others that do not seem to display any resonant behavior. Among the oscillating cells, many exhibit a low-frequency rhythm in the 4–10 Hz range [52]. Our modelled IO cells oscillated robustly around 5.8-6 Hz and were connected with electrotonic Connexin-36 synapses [52,75]. IO spiking synchrony was disrupted in absence of gap junction coupling (S1 Fig). Though anatomical clusterization of cells is highly plausible [55,76], here we examined a small set of cells without explicit clustering.

Homeostatic plasticity and CSpk-triggered plasticity mechanisms

To capture the essential interactions between LTP and LTD we modelled the PF-PC synapse with a mixture of a BCM learning rule and a PF input dependent CSpk-triggered LTD mechanism. This CF-induced LTD operates inversely to the classical BCM mechanism and changes the probability of BCM-LTP [12]. In our model, the synaptic weights were potentiated when PC activity was low to maintain excitability. This complex combination of synaptic plasticity mechanisms enabled us to inspect their interactions and dynamics.

The BCM mechanism (Fig 2A ) in our model encapsulates the effect of changes in the PF activity by either potentiating or depressing the synapse independently of CF input. This mechanism depends on 3 factors: the average of the recent activity of the PC, ; the BCM threshold — which also depends on — and the instantaneous activity of the PF, (Fig 2B ). After an IO spike (Fig 2C ), a pause in SSpks leads to a CF pause, which causes a decrease in . The LTD component of plasticity due to the CSpk gives a negative change in synaptic weight also proportional to the (Fig 2D ). The plasticity threshold distinguishing depression from potentiation follows , that is lower than results in potentiation. Conversely, higher than leads to depression. Both these changes are proportional to activity for each synapse (see Fig 2E ). The total change in synaptic weight is computed as the sum of two plasticity components: a BCM-based continuous update and a CSpk-triggered term (Fig 2F ). The BCM plasticity mechanism has only one hard constraint limiting the potentiation of the PC to frequencies of 250 Hz, as firing rate beyond this rate is considered pathological.

After an initial transient, synaptic weights converged to stable distributions. However, we observed small residual stabilizing oscillations around this stable mean, which we refer to as weight fluctuations. These reflect the continuous balancing process of BCM potentiation and CSpk-triggered LTD under stochastic input. Importantly, these fluctuations were minute relative to the absolute weight values and did not show cumulative drift. While shaped by the structure of the input —such as frequency content and temporal correlations— the fluctuations were not dependent on the input duration, indicating a steady-state regime rather than input-length dependence.

While the BCM rule used is not fitted to in vitro data, it captures key features of activity-dependent potentiation and depression shaped by CSpk activity. This allowed us to examine how loop-level dynamics guide synaptic trajectories. Different BCM parameter choices may alter the speed of convergence but not the main results: weights stabilized, CSpk frequency remained frequency-dependent, and Up/Downbound modules retained distinct profiles. For example, we selected the threshold time constant () within a range of 10–20 ms, based on the typical duration of CF pauses. This parameter governs the low-pass filtering of recent PC activity used to compute the threshold . We did not tune to produce specific outcomes; rather, this range was chosen a priori based on known physiology, and only a posteriori did we find that it yielded expected plasticity directions —net depression in Downbound zones and potentiation in Upbound zones.

Synaptic weights are stable across plasticity epochs

We hypothesized that synaptic weights would reach a new level of stability and stay at the level for subsequent plasticity events. To investigate these changes in synaptic weights across plasticity events we did the following simulations (see Methods), where for each epoch we gave the same frozen input as PF input. We started with the first epoch ("no plasticity") with the plasticity turned off and with static synaptic weights which were constant throughout this 120 s experiment. After this "no plasticity" epoch finished, BCM and LTD plasticity mechanisms were enabled and the same frozen input was replayed. At the end of this plasticity epoch, we took the last value of the filtered synaptic weights and turned it into the new static weights for the final epoch (after plasticity). Again the same input was given one last time with the new static weights. This configuration enabled us to compare synaptic weights at the PF-PC synapse before and after plasticity. To evaluate the impact of plasticity on the weights, each experiment was performed 4 times with the same frozen input (Fig 2G-L).

The same frozen input was used in each epoch. It provides consistent, varying input statistics and temporal correlations across epochs ensuring reproducibility. When plasticity is switched off, the synaptic weights remain fixed, allowing us to compare resonant activity at stable weight configurations before and after plasticity. In contrast, during plasticity-on epochs, the BCM and CSpk-triggered LTD mechanisms allowed the weights to settle into a new homeostatic equilibrium.

For both Up/Downbound zones the direction of change of the 500 synaptic weights for the 100 PC population (5 synapses per PC) was defined during the first "plasticity" epoch. The synaptic weights at this epoch were reduced in the Downbound zone and increased in the Upbound zone (Fig 2I and 2J ), suggesting that the intrinsic PC currents determined the plasticity direction. From the second epoch after plasticity we found that the weights were mostly stable, though a few changed in rank (Fig 2G and 2H ), indicating that the interaction between CSpk and PF inputs may drive these rank changes. In Downbound zones weights fluctuated around the mean, indicating a faster rate of weight change, due to the higher intrinsic PC current (Fig 2J ). The Upbound zone, with a lower intrinsic drive, exhibited more synaptic weight stability after relaxation in epochs 3 and 4 (Fig 2G and 2H ).

To test the robustness of these weight dynamics, we repeated the simulation protocol using 10 different random seeds, each with independently generated OU input and randomly sampled network parameters. For each seed, we computed the median weight change per plasticity epoch (Fig 2K and 2L ). The same bidirectional trend was observed across all simulations: Upbound modules exhibited early potentiation followed by stabilization, while Downbound modules showed early depression followed by stabilization. These results support the generalizability of our findings and indicate that the observed stabilization arises from the intrinsic dynamics of the model rather than from input-specific contingencies.

Gap junction coupling creates a global phase and leads to synchronized firing

We expected gap junctions (i.e., "coupling") to have an impact on IO STO synchrony with strong coupling leading to higher synchrony. To quantify this effect, we looked at the STO oscillation phases over time. These were decomposed into a network level (global) phase, defined as the mean Hilbert phase across all IO neurons, and into an individual-level (delta) phase, defined as the difference between the Hilbert phase of each IO neuron and the global phase at each time point. This was done to compare the difference between the individual and global phases (S2A Fig). The global phase of STOs, representing average phase progression, was very stable in both coupled and uncoupled conditions (S2B, S3B Fig). In the condition with coupled IO cells, oscillators phase-locked onto the global phase creating synchrony (S2C Fig) not seen in the uncoupled case (S3C Fig), where oscillators behaved independently. To investigate whether the coupling also led to dynamically organized clusters in the IO, we calculated the pairwise Kuramoto phase (S2D Fig). At the scale of one second (S2E Fig), transient clustering could be observed, however clusters were not stable (S2F–G Fig) and reflected the random phase fluctuations driven by GABAergic CN input. At the scale of 30 s (quarter-time of the full simulation) no clear clustering was observed. Network synchrony peaked around IO spikes, but only for the coupled scenario (S3A Fig).

Upbound PCs have a propensity for potentiation and Downbound PCs for depression

As it has been shown that PCs adaptively adjust their SSpk firing rate via CSpks to elicit blinks during eyeblink conditioning [45,81,82], we hypothesize that the direction of change of the PC firing rate depends on the intrinsic frequency of each PC microzone. Eyeblink conditioning is a classical associative learning paradigm where a neutral stimulus, such as light, is paired with an air puff that evokes reflexive eyelid closure so that the initially neutral stimulus can eventually trigger blinking. As argued previously, in this associative task, bidirectional changes in the SSpk activity are expressed in the Downbound zone during initial stages of training, whereas the Upbound zone controls acceleration of movement in paradigms such as the vestibulo-ocular learning (Fig 3A) [44]. In fact, several studies support the distinction between these zones, proposing that they constitute complementary modules that may act in concert to encode prediction and feedback signals during sensorimotor learning [5,9,43,47–49]. Notably, while eyeblink conditioning involves learning over many trials with discrete, stimulus-locked inputs and explicit error signals, our model instead uses fixed, temporally varied, OU input, which was generated once and reused across simulations to explore whether intrinsic firing properties and internal dynamics alone can explain the directionality of plasticity. Fig 3A provides a schematic summary of previously published experimental findings [44], illustrating the typical direction of SSpk modulation in Up/Downbound PCs during eyeblink conditioning and vestibulo-ocular learning. In contrast, Fig 3B shows how the same directionality emerges in our model when plasticity mechanisms (BCM, LTD, LTP) are activated, but over much longer timescales (seconds rather than hundreds of milliseconds. For display, we plot the raw 1-s binned population mean together with a Savitzky-Golay-smoothed trace. The smoothing applied across the full concatenated simulation produces small boundary effects at the epoch transition. Nevertheless, the comparison highlights that, despite the different inputs and timescales, the baseline level of PC activity is sufficient for the BCM mechanism to establish the direction of synaptic weight change. PCs in our model were tuned to reproduce the experimental data (Fig 3C and 3D), whereby Downbound PCs showed relatively high SSpks (~90 Hz) and CSpks (~1.3 Hz) firing frequencies, while Upbound PCs had lower average SSpk (~60 Hz) and CSpk [(~0.9 Hz), [9,44]]. Tuning PCs intrinsic currents resulted in the CF pause reproducing experimental observations (Fig 3E).

The direction of plasticity was the same irrespectively of the absence or presence of IO gap junction coupling (S4A-E Fig). For the Downbound micromodule, there was a small increase in the mean of the CSpk firing frequency in the uncoupled case with respect to the coupled case (Figs 3D and S4B), which led to a wider distribution of the CSpk frequency squared coefficient of variation [(CV2); (S4F Fig)].

Inhibition followed by disinhibition of CN elicits CSpk

We hypothesized that specific characteristics of the input would elicit CSpks. We started by analysing the CSpk-triggered PSC (Fig 4A) and we observed that the CSpks were promoted by specific oscillations present in the input. In our model, the PSC refers to the current that enters the PC soma,specifically, the sum of the 5 PF currents, each scaled by their corresponding synaptic weight. CSpks were evoked when the input carried specific phases and frequency (around 5–6 Hz), with a preference for reduced input just before the CSpk. At the level of the CN this translated to inhibition followed by disinhibition, which had previously been empirically observed by [83,84]. The decrease (following an increase) in PSC happened 60 ms before the bottom of the subthreshold oscillation (STO) of the IO membrane response, which was closely related to the synaptic delay in the loop (10 ms between PC and CN and 50 ms between CN and IO). This suggests that CSpks had a preference for particular frequencies and phases of the input. A clear contrast is observed after plasticity, as in Upbound zones the sensitivity to the PSC oscillations was increased and reduced for Downbound zones.

Every CN spike contributed to a reset of the IO STOs (Fig 4B), in line with experimental results [16]. We found that in the Upbound zone there was a lower CN spike frequency and lower amplitude of IO oscillations before the CN spike. As there was more CN inhibition in the Downbound zone, a lower CN spike frequency led to reduced dampening of the IO oscillation. In addition, we found that increasing CN inhibition lowered the frequency of IO STOs while increasing IO spike (CSpk) output (S5E Fig). In the Upbound zone, we observed a lower CN spike frequency and lower amplitude of IO oscillations before the CN spike, indicating weaker modulation of IO dynamics. In contrast, in the Downbound zone, stronger CN inhibition slowed the IO STOs and reduced their amplitude, creating conditions that increased the likelihood of IO spiking. This may be due to reduced oscillation amplitude, resulting in milder hyperpolarization and thus a greater chance of crossing spike threshold. Thus, CN activity modulates IO excitability indirectly, through effects on STO phase and amplitude, rather than through direct excitation. These nonlinear interactions were further shaped by PF-PC plasticity: IO responses to CN spikes varied depending on prior synaptic changes, showing that the same input could produce variable IO responses depending on network history. This is consistent with experimental observations of trial-to-trial variability in CSpk timing [7,85], and supports the idea that IO STO phase and amplitude influence CSpk generation [52]. After plasticity, CN spikes no longer significantly affected IO STO frequency in the Downbound zone. This result suggests that the IO STOs are affected but not driven by CN inhibition.

After the first "plasticity" epoch, the PC SSpk activity was increased for the Upbound zone and decreased for Downbound zone, resulting in both PC populations firing at about 60 Hz after plasticity (Fig 4C). As a result, the CN activity was decreased for the Upbound and increased for the Downbound zone (Fig 4D). Furthermore, the CSpk frequency was decreased for Upbound PCs and increased for Downbound PCs (Fig 4E). The length of the CF pause decreased for both Up/Downbound zones (Fig 4F).

Investigating the effect of gap junction coupling, we observed that the uncoupled Downbound IO had on average a higher variance and higher CSpk firing frequency (S5B Fig), leading to a higher variance of SSpk firing frequency for Downbound PCs (S5A Fig). Nevertheless, the CN (S5C Fig) and CF pause (S5D Fig) were the same for both coupled and uncoupled conditions.

Specific input frequency bands differentially drive plasticity

Our previous hypothesis was that specific characteristics of the input would directly influence CSpk generation. Here, we expand our focus to include the full olivocerebellar loop, proposing that distinct input frequency bands engage different resonant dynamics of the circuit, thereby shaping both CSpk timing and the resulting plasticity. As a continuation of the results shown in Fig 4, we examined how different frequency bands shape loop-level responses under plasticity. Signals from different sensorimotor sources have different spectral properties, and PC population SSpks activity has been reported to contain bands from low delta [86] to high gamma [50] frequencies.

To investigate how these spectral differences influence the loop dynamics under plasticity we filtered the original PF signal into distinct frequency bands. The filtered signals were normalized to keep the same spectral power (Fig 5A), and were used in exact copies of the same network to compare the before and after plasticity activity of the cells across the different frequency bands (Fig 5B-D). We found that CN dynamics were hardly modulated by the filter classes, indicating high level of signal divergence and mixing (Fig 5C).

As judged from both SSpk (Fig 5B) and CSpk (Fig 5D) distributions, Up/Downbound PC zones showed marked differences in how their firing rates responded to different PF input frequency bands. These differences were most pronounced for SSpks in the upbound zones.

While the direction of plasticity of CSpks activity qualitatively resembled physiological observations [9,43,44], frequency distributions were sensitive to plasticity under different bands, with lower input frequencies driving most of the variability. This was true for low frequencies (5–10 Hz band), where the CSpks activity was higher for Up/Downbound zones. Furthermore, frequency bands larger than 10 Hz evoked similar SSpks and CSspks distributions after plasticity, indicating that the PC network response was robust to frequency variations.

Differential CSpk response to PF input changes drives plasticity in 5–10 Hz band

The largest differential plasticity responses were observed for the 5–10 Hz filtered input in the Downbound zone. Here we observed a large decrease in SSpk frequency after plasticity (Fig 5B), with a twofold increase of CSpks (S6A Fig). Before plasticity, the SSpk frequency was approximately the same for the 5–10 Hz filtered input and the "no filter" case, as the current input power spectrum and weights were the same. We assumed that the decrease in SSpk frequency after plasticity originated from CSpk-triggered LTD. However, this assumption did not explain the increase in CSpk frequency for the 5–10 Hz filtered PF input.

We looked at the relationship between the STO frequency and CSpk frequency (S6B Fig). We found that STO frequency and frequency-distribution width was narrower in the Downbound than in the Upbound zone. For both Up/Downbound the STO frequency and IO spike frequency were inversely related. This may seem unintuitive but is explained by channel dynamics of IO cells, as more CN inhibition both increases the driving force and reduces Ca-T type inactivation, leading to more prompt spiking upon disinhibition from DCN (S5E Fig) [16,52,87]. For both the Up/Downbound zone, 5–10 Hz input lowered the STO frequency further, while at the same time increased the IO spike frequency. The lowering of the STO frequency was more pronounced in the Upbound zone. We examined the IO spike-triggered PSC to understand further how the input characteristics affected CSpk occurrences (S6C Fig). In general, for both IO unfiltered and 5–10 Hz input, IO spikes were preceded by oscillations below or at the lower range of the STO frequency (5.6 Hz for Upbound, 5.4 Hz for Downbound). For unfiltered noise, IO spikes were triggered by a small decrease in PSC amplitude, whereas for 5–10 Hz input, IO spikes were triggered by a shortening of the PSC oscillation period by 10 ms (S6C Fig). Our results suggest that, for the unfiltered input (OU input), the smaller frequencies, which are similar to the STO frequency, have a bigger impact on the system.

Gap junction coupling drives frequency selectivity

Following our previous hypothesis that gap junctions affect IO STO synchrony, we hypothesized that this change in STO synchrony will be reflected in the loop’s response to different frequencies. We observed that the frequency selectivity of the IO was lost when the gap junction was uncoupled (Fig 6A). The differences in the mean and distribution of CSpk and SSpk frequency after plasticity were much less pronounced within the zones in the uncoupled condition than in the coupled case. Furthermore, the higher IO spike frequency and lower STO frequency was completely lost for the Upbound zone, but was increased for the Downbound micromodule when the IO neurons were uncoupled (Fig 6B). The loss of the inverse STO - IO spike frequency relationship in the Upbound zone co-occurred with a loss of the oscillatory response in the IO-spike triggered PSC (Fig 6C). A higher CSpk frequency for the Downbound zone (S4B Fig) led to a less pronounced loss of the inverse STO - IO spike frequency relationship. This suggests that the gap junction coupling has an effect upstream to the PC-CN and back to the IO as there is less variation in the CN inhibitory input onto the IO.

Baseline Up/Downbound states support task-related plasticity during conditioning

Our model reproduced bidirectional plasticity in PC activity in response to arbitrary inputs (Fig 3B ). While frozen noise may approximate inputs in the awake and behaving state, it does not reflect specific motor behaviors. To test whether the stable Up/Downbound baseline states generated by the model could support learning-related plasticity, we reproduced a classical eyeblink conditioning paradigm [88]. A conditioned stimulus (CS, 250 ms) was delivered to one PF bundle and co-terminated with an unconditioned stimulus (US, 30 ms) delivered to another PF bundle. A copy of the US was sent to the IO to elicit CF input, as in standard protocols (Fig 7A). Each block consisted of an initial US-only trial, ten CS + US pairings, and a final CS-only trial, separated by 8-12s intervals (Fig 7B).

The model captured key features of learning-related plasticity. The US elicited CSpks across the IO population (Fig 7C), and STO synchrony increased in the Upbound zone after plasticity (Fig 7D). Before learning, SSpks increased in response to the CS; this increase disappeared after learning (Fig 7E). These dynamics were consistent with the baseline states in Fig 4, as seen in before/after plasticity comparisons of SSpks, CSpks and CF pause (Fig 7F-H).

Some discrepancies with experimental data remained. We did not observe a post-learning decrease in SSpks, likely due to the absence of MLIs, which can suppress PC activity in response to increased PF input [89,90]. Nor did we observe the CS-aligned CSpk increase reported experimentally, which may depend on modulatory input from the striatum or amygdala [91], both absent from our model. However, our model did replicate the experimentally observed decrease in CSpk frequency in response to the US in the Downbound zone after learning, suggesting that this effect arises from intrinsic olivocerebellar dynamics.

STO frequency decreases after plasticity for CS-US experiment

To better understand the learned response during CS-US trials, we examined PF weight changes over time. PF1–3 and PF5 (US) exhibited similar weight dynamics (S7A Fig), while the PF4 (CS) weight showed a pronounced decrease during plasticity, likely explaining the lack of CSpks at CS onset (Fig 7E). Trial-averaged weights (S7B Fig) confirmed this pattern: the PF4 weight declined steadily across CS trials, whereas the other PF weights remained stable. During the US, PF1–3 weights briefly decreased, then rapidly returned to baseline, while PF5 remained unaffected. Before plasticity, increased PC firing transiently accelerated STO frequency and global oscillator phase (S7C Fig), an effect abolished after plasticity.

When gap junctions were uncoupled, IO STO synchrony decreased overall but still peaked near the US (S8A Fig). SSpk responses were similar to the coupled case, but the CSpk pattern reversed: the Downbound peak increased while the Upbound decreased post-learning (S8B Fig). Uncoupling also led to higher IO firing rates (S8C Fig).

Steady state distributions of PC activity

Firing rates measured under resting conditions in awake animals likely reflect the stable physiological operating point of the olivocerebellar loop, shaped by continuously adapting plasticity mechanisms. These baseline configurations should not be interpreted as ‘pre-plasticity’ or naïve states, but steady states that emerge from continuous input, intrinsic excitability, and ongoing synaptic plasticity. This interpretation is supported by experimental evidence: spontaneous firing rates in Up/Downbound zones are stable over tens of seconds in adult animals [9], and zonal differences in firing rate and expression of plasticity-related markers are already present during early postnatal development [92]. These data suggest that zone-specific firing characteristics are not simply acquired through experience-dependent learning but are likely developmentally programmed and maintained through homoeostatic plasticity.

In our model, PF-PC synapses were subject to two types of plasticity, BCM and CSpk-dependent LTD. The BCM rule is a Hebbian activity dependent learning rule, and as such is expected to reflect statistical regularities in the input [93,94]. However, in a biophysically plausible setting with a large set of parameters and feedback there is no guarantee that IO and PC responses remain bounded over time. To examine weight convergence, we verified the weight values after multiple epochs of repeated, stochastic input, modeled as OU noise. We found that synaptic weight distributions and firing rates stabilized across epochs for both Up/Downbound micromodules. After an initial transient, Upbound zones displayed more weight stability across epochs, while Downbound modules showed small, bounded fluctuations around a stable mean, and thus were more sensitive to loop dynamics. In effect, ongoing plasticity tuned PC activities in the same direction as observed in physiological experiments [9,43,44] (Fig 3); thus our model showed how module-specific intrinsic PC firing rates shaped PC-PF plasticity direction. Plasticity did not simply preserve the initial intrinsic firing asymmetry, but instead reconfigured PC activity through closed-loop dynamics. Preferred activity levels emerged as a result of interactions between intrinsic excitability, ongoing input, and network feedback. In this sense, the loop does not enforce a fixed firing rate signature per module, but rather shapes zone-specific attractors under homeostatic plasticity. The intrinsic current established a rate preference that acted as a corrective "elastic force" on the weights. In the case of the Downbound module, this corrective drive was less efficient at damping weight fluctuations, leading to slower convergence and greater residual variability, whereas in the Upbound zone, convergence was faster and more stable. This leads us to predict a higher changeability of synaptic weights in Downbound zones, as the PC intrinsic rate is kept farther from circuit equilibrium in this case. This model prediction is in line with the biological observation that proteins implicated in synaptic plasticity are more prominently expressed in the Downbound zones, than those of the Upbound zones [43,44,95]. Notably, recent transcriptomics studies indicate that the subdiversity of different PC microzones is even greater than the two main categories of Up/Downbound modules, particularly in primates [42,96,97]. Hence, future models will need to consider reflecting the full heterogeneity of all PC populations.

Using our model’s simulations we can conclude that while zone-specific PC characteristics are known to be intrinsic [9,42–44], the olivocerebellar loop mechanisms stabilize these differences through the interaction of BCM and CSpk-dependent long-term plasticity. However, the precise direction of baseline firing rate divergence between Up/Downbound modules may also depend on additional biological mechanisms not captured in our model, such as MLI inhibition or structured MF input.

Input frequency dependent CSpk responses

Changes of weight of specific PCs in our model were strongly determined by their recent input, both in terms of PF and CF input. Recent firing rate acted as a low pass proxy to history of inputs, while the CSpk timing was strongly determined by the phase state of the IO. We observed a strong preference for specific frequencies and phases of the stochastic PF input and showed that CSpks probabilities were modulated by specific temporal features of this input (Fig 4A and 4B). As the state of the IO was subject to CN population drive, we showed that the plasticity of the overall system was frequency-dependent.

Our choice of using OU PF noise [68] was motivated by emulating arbitrary input that may exhibit any frequency and any phase combination. In the long interval of our plasticity epochs (120 s), this OU input contained very different frequency compositions. This modeling strategy was chosen to reflect in vivo experimental data of granule cell and PF activity during learning, at rest and during walking [66,67,69,70], which show complex, high dimensional firing patterns transmitted into PCs in awake animals. The variability of the IO feedback and the resulting jitter in the latency of the CSpk response across trials, which we show in our model, has also been observed experimentally [7,85] and can be explained by phase shifts of IO STOs [52]. We analyzed the extent of the input frequency dependence by filtering inputs to include only specific frequency bands. We discovered significant differences in SSpk distributions to filtered inputs, which was particularly apparent in Downbound zones. Specifically, we found that our model responded preferentially to the 5–10 Hz frequency band, which overlapped with the IO STO frequency (5–6 Hz), indicating that specific ranges of PF frequency bands led this model closer to physiological responses.

This frequency preference suggests that IO output is maximized when PF input contains components near the intrinsic STO frequency of the IO (4–6 Hz). When STOs are aligned with this input band, IO neurons are more likely to spike, thereby increasing the probability of IO spiking. When we further investigated the relation between IO STO frequency and IO spike output, we found them to be inversely related in both Up/Downbound modules. This arose from differences in inhibitory CN feedback: stronger inhibition in the Upbound zone slowed the STOs, effectively shifting them away from the preferred frequency range and thus reducing the likelihood of IO spiking. In contrast, the Downbound zone received weaker CN inhibition, resulting in faster STOs that remained better aligned with the preferred frequency range, facilitating more frequent IO spikes. This mechanism demonstrates how CN inhibition indirectly gates plasticity through STO tuning, and aligns with in vivo recordings showing that CSpk rates differ across modules [7,85], with Upbound zones typically displaying lower spontaneous CSpk frequencies than Downbound zones.

Stability of PF-PC synaptic weights

The overarching question that motivates our model is the apparent dilemma between stability of acquired sensorimotor patterns and ever changing synaptic strengths. However, we did not specifically test the loop in a behaviorally-relevant setting, rather we examined the general case of any arbitrary input, reflecting ongoing sensorimotor behavior. Future research with our model could study input signals that more closely represent the contingencies of sensorimotor learning, with their specific frequencies and phases, related to the interaction to the body and other brain regions. Upcoming effort could focus on adding the excitatory input to the IO from the mesodiencephalic junction (MDJ) [16] and even the CN [13]; oscillatory gating of the granular layer via Golgi cells [10,98]; and/or inhibition from the MLIs [89,90]. While we deliberately excluded the excitatory CN-MDJ-IO pathway and MF collaterals to the CN to focus on the core feedback loop, we acknowledge that these inputs are likely to modulate IO excitability and contribute to zonal differences. Experimental evidence shows that excitatory input via the MDJ can evoke CSpks with much shorter latencies (4–8 ms) [15] than the slower GABAergic inhibition of over 50 ms [79,99,100], and that MF input to the CN may shift baseline excitability.

Our model was not tuned for specific outcomes but to test whether stable synaptic states could emerge from internal loop dynamics. Resonance between PF and CF activity drove frequency-dependent PF-PC synapses stabilization without task-specific tuning. We interpret this as a baseline regime of the olivocerebellar loop, where Up/Downbound modules settle into distinct steady states under stationary input and ongoing plasticity.

We extended this framework by incorporating a classical eyeblink conditioning paradigm. The baseline states established by ongoing plasticity persisted when task-related CF inputs were introduced. During CS-US pairings, Up/Downbound modules underwent distinct learning-related adaptations, demonstrating that homeostatic stabilization provides a functional substrate for supervised, CF-driven plasticity. This reconciles experimental observations of stable baseline states with the emergence of biased plasticity during learning.

Impact of gap junction coupling on stabilization

Gap junctions can cause synchronized IO spikes, in both Up/Downbound zones, leading to a global phase in which many individual IO neurons are coupled [16,101]. When the gap junctions were abolished in our model, the frequency selectivity of the IO response was disrupted. In the Upbound module, this disruption was accompanied by a reduction in the variability of IO spike rates and a loss of oscillatory structure in IO-triggered PSC responses. However, in the Downbound module, IO spike variability persisted to a larger degree, likely due to its intrinsically higher IO and PC firing rates and more robust CSpk activity. This differential effect underscores the zone-specific contributions of IO electrotonic coupling to timing and frequency selectivity. Moreover, this result suggested that the gap junctions have a critical influence on the frequency selectivity of the system. Indeed, when tested, we found that the frequency selectivity due to the difference in strength of the CN input, highlighted above, was lost in the absence of IO gap junction coupling. As CN inhibition turned out to not completely shape the IO STOs, this could indicate that the interplay between CN inhibition and IO gap junction coupling is at the core of this STO - IO frequency relationship. It should be noted that the interaction between GABAergic CN input to the IO and regulation of the coupling strength as well as the relevance of coupling in the IO for the prominence of its STO’s has been the subject of many experimental studies. Both relations have been confirmed at the morphological, electrophysiological and pharmaceutical level, in vivo as well as in vitro [15,16,52,77–80,86,102–106]. The relationship between the level of STO’s and IO frequency band is underscored by the finding that Up/Downbound modules showed different correlation strengths in this respect. Thus, as gap junction coupling facilitates IO resonance and as there are module-specific interactions, our data imply that the temporal pattern selectivity processing across Up/Downbound micromodules is enhanced by dynamic regulation of electrotonic coupling in the IO. This was further confirmed in the eyeblink conditioning paradigm: uncoupling reduced overall IO STO synchrony but preserved peak synchrony around the US. Interestingly, CSpk responses reversed post-learning Downbound increase and Upbound decrease, highlighting how coupling biases plasticity directionality. IO firing rates also increased without coupling, suggesting that electrotonic interactions constrain IO excitability and shape zone-specific plasticity. These modeling data align well with experiments in which olivary coupling is affected either genetically or pharmacologically; impairing dendrodendritic IO coupling hampers cerebellar sensorimotor learning [75].

Limitations and biological generalization

Our model captures key internal dynamics of the olivocerebellar loop, however due to epistemic considerations several biological features were simplified or explicitly excluded from the design. The PF-PC synaptic plasticity was modeled using a BCM-inspired homeostatic rule in combination with a CSpk-triggered LTD component, while leading to LTP in the absence of CSpks. While this captures bidirectional plasticity and weight stabilization, it remains a phenomenological approximation. Experimental data on PF-PC plasticity, particularly under in vivo or awake conditions, suggest a large set of mechanisms that our simplified rule does not fully represent [23–27]. For instance, the CF CSpk itself was not modeled with detailed waveform dynamics, and thus the influence of CSpk duration, calcium influx or burst structure on plasticity was not explored [107,108]. In terms of input, we employed stochastic noise (OU process) to test how arbitrary fluctuations at different frequencies shape loop dynamics and plasticity. This does not capture structured MF activity, such as bursts or sequences associated with sensory events or behavior. Nor does it reflect modulations by Golgi-Granule cell-mediated gating oscillations in the 8–12 Hz range [10,98]. However, we observed that our model responded most strongly to frequencies near 10 Hz — consistent with the filtering properties of the granular layer [98]. This suggests that Golgi cell-induced frequency gating may reinforce IO resonance at these frequencies and further constrain learning windows.

To assess whether the stable Up/Downbound configurations could support learning, we modeled a classical eyeblink conditioning paradigm typically associated with Downbound zones [5]. We applied the paradigm to both modules and found that CS-induced SSpk increases were suppressed after learning, and US-triggered CSpks decreased in the Downbound zone. Notably, these effects reversed when IO gap junctions were removed. While the conditioning paradigm captures key elements of learning-related CF modulation, it does not implement full IO-gated sensorimotor feedback loops that modulate CF activity. A recent model, incorporating MLIs, Golgi and granule cells replicated CS-induced SSpk suppression and CSpk increase during conditioning [45], but did not examine the role of IO reverberations on PC homeostatic plasticity. Our model shows that the reduced synchrony from removing gap junctions increased Downbound CSpks and reduced Upbound ones, suggesting IO coupling is required for the timing of acquired eyeblink responses [75].

The simplifications in our model were designed to probe ongoing learning properties of the resonant olivocerebellar loop, from a reverse-engineering perspective, under minimal assumptions about what is to be learned. The model captures homeostatic stabilization under persistent input, not task-driven synaptic reorganization. Future versions could test how structured task-dependent sensorimotor contingencies shape weight homeostasis, clarifying stabilization versus learning-related plasticity. Overall, our work highlights the dynamic nature of this loop, essential for understanding how the cerebellum "learns to learn".

Network architecture

To balance anatomical accuracy with computational tractability, we modeled each micromodule with 100 PCs, 40 CN neurons and 40 IO neurons. These ratios are grounded in known anatomical estimates: within a module, approximately 20–30 PCs converge onto a single CN neuron, while each PC projects to ~30–50 CN neurons [58]. Half of IO neurons connect to the PCs. Each of those IO neurons project to ~7–10 PCs, and each PC receives CF input from a single IO neuron [59,60]. We modeled only the inhibitory CN neurons, as these are the ones that project to the IO [53–57]. For the CN-IO synapse, we assumed a 25% reciprocal connectivity probability, resulting in each CN neuron connecting ~10 IO neurons and each IO neuron receiving from ~10 CN neurons. This structured yet simplified architecture ensures that core convergence and divergence patterns of the circuit are preserved while enabling efficient simulation (see Fig 1).

Neuron models and justification

PCs and CN neurons were modeled using simplified adaptive exponential integrate-and-fire neurons. This abstraction was sufficient to capture the key dynamics relevant for our study: the modulation of PC firing rate by PF and CF input, and the inhibitory feedback from PCs to CN. Furthermore, IO neurons are endowed with rather unique properties. They are one of the most densely coupled via dendro-dendritic gap junctions than any other type of neuron in the mammalian brain [102,103,109] and they express high levels of voltage-gated and calcium-dependent ion-channels that allow them to oscillate at frequencies of up to 10 Hz [52,77–80]. The interplay between the excitatory and inhibitory inputs to the IO as well as the temporal relation with the subthreshold oscillations (STOs) of the coupled olivary cells determines whether and when an olivary action potential will be generated [16,64,105,110]. Because our primary focus was on IO-driven resonances and loop-level plasticity, we prioritized a more detailed multi-compartment IO model to accurately reproduce STOs and phase-resetting phenomena. Modeling PC morphology or PC calcium dynamics was not necessary for our hypothesis (see Introduction), which depended primarily on CF timing and population-level PC output.

Parallel fiber bundles

We simulate the PF input, representing the activity of PF bundles, as an Ornstein-Uhlenbeck (OU) process:

where is a Gaussian white noise with mean and standard deviation that scales with units of . The mean is and the variance —the integral of the power spectral density— is equal to . The noise intensity is defined as . Parameter values in the model are shown in Table 1.

PC and CN neuronal models

The PC and CN populations are modeled as adaptive exponential integrate-and-fire (AdEX) models [111]. The general cell equations for the PC and CN are described as:

Here, Eq. 2 describes the membrane potential dynamics, while Eq. 3 shows the dynamics of the adaptation value. The refers to the sum of the currents specific to cell type X (where and indexes the current types affecting that cell); is the slope factor; is the threshold potential; is the adaptation variable; is the adaptation coupling factor, and is the adaptation time constant. The slope factor is a quantification of the sharpness of the spike, which can be seen as the sharpness of the sodium activation curve if the activation time constant is ignored [112]. The parameters of all AdEx cells (both CN and PC) were sampled using evenly spaced values across the listed range in Table 2. For each parameter we generated a list of linearly spaced values equal to the population size (e.g., 100 for PCs), and assigned them to individual cells without replacement. This approach ensures a uniform spread across the specified range and introduces controlled variability in intrinsic properties across the population. Hence, each cell in the population will display individual differences in firing properties. When a spike happens both the membrane potential and the adaptive variables are updated. The actual threshold potential is used for the neuron models. When the potential reaches (see Table 2), a spike occurs and both the membrane potential is reset to , and the adaptation variable is updated to the reset adaptation current and then decays until either reaching 0 or the next spike time [113].

Our goal for using the randomization of parameters is to have variability across the cells in each population. Since the CN is inhibited by the PC we will already have enough variability in the inhibitory current as each PC is different. Also, the CN cells have their own variability across the parameters. Hence, we do not randomize the intrinsic current of the CN as we can then see the direct effect of one PC onto the CN.

Purkinje cell Currents

For the PC, there are 2 currents denoted and . The former is the activation (intrinsic) input current given to the cell to tune the desired intrinsic firing frequency distribution. The latter is the postsynaptic current (PSC) input entering the PC soma from the PFs (Fig 1A).

Cerebellar nuclei Synaptic Currents

The CN currents are shown in Eq. 6, where is the intrinsic current of the CN and is the inhibitory current from the PC to the CN. It is a filtered version of the PC spike train that decays exponentially in the absence of spikes. Upon each PC spike, is updated by a constant (see Table 6). This means that the total current will be smaller and the CN will spike less (inhibition from PC).

Inferior olive model

The IO model is based on the De Gruijl model [108] —a modification of the Schweighofer model [87] with an added axonic current. The general cell equations for the model are shown below:

Where, the membrane potential specific to each compartment, somatic s, dendritic d and axonic a is represented with (where ). The equations for the ionic conductances of the general cell are shown in Table 3 and their parameters in Table 4. All parameters are randomized in the same way as for the AdEx parameters. The IO receives an a noise input modeled as an OU process (similarly to the PF input, see Table 1). With this input we model the baseline firing rate of the IO as well as modeling inputs from other brain regions. Furthermore, the inhibitory input from the CN onto the IO is modeled by the current , as described in the Synapses subsection. Finally, two applied currents (where ) are in the dendrites and somatic compartments. These are used when adding external current into the cell. The dynamics of each activation or inactivation variable (see Table 5) is given by the differential equation below. In this equation indicates the variables and . There are 2 exceptions, and which have a small enough that they are approximately equal to their steady-states ( and , respectively).

Three morphological parameters determine the electrotonic properties of the three-compartmental model: the ratio of the somatic area to the total somatic and dendritic surface area ; the ratio of the somatic area to the total somatic and axonic surface area , and the electrotonic coupling conductance between the 3 compartments. The current flowing from the somatic into the dendritic compartment and from the dendritic into the somatic compartments , as well as the leakage currents for the soma and for the dendrites are represented in this model [87]. We also have the current flowing out from the soma into the axon and from the axon into the soma as well as the leakage for the axonic compartment [108]. Furthermore, the dynamics of the calcium concentration is defined as follows:

Two distinct types of calcium currents exist in the IO neuron: a low-threshold current located in the soma and a high-threshold current located in the dendrites [77,114]. As the window of conductance of is around the resting membrane potential, the IO cell is excited in response to hyperpolarizing currents [77,114]. On the other hand, is non-inactivating. Thus, a depolarizing dendritic input results in a prolonged plateau potential. As calcium enters the dendrites, the calcium-activated potassium current is activated. This current abruptly terminates the plateau potential after about 30 milliseconds [77,114]. Due to the long time constant (several hundred milliseconds) of , there is a long afterhyperpolarization which, in turn, de-inactivates resulting in a postinhibitory rebound spike. This current contributes to the subthreshold oscillations as it contributes to amplitude and frequency [87,115]. Furthermore, somatic sodium spikes are generated with the sodium current . These are terminated by an outward delayed rectifier potassium current [77,114]. Finally, the axon hillock includes a sodium current (whose inactivation function was altered to allow for the generation of bursts of spikes) and potassium current [108].

Synapses

Table 6 shows the parameters and synaptic delays for each synapse in the model. The synaptic delay and time constants were chosen from literature for the PC to CN and the CN to IO delays [100,108,116–119]. Other parameters were selected to reproduce the electrophysiological results shown in Fig 3. Synaptic parameters in Table 6 were held constant across simulations to isolate the effects of plasticity mechanisms and network architecture. While neuronal parameters in Tables 2 and 4 were varied to reflect known biological heterogeneity, fixed synaptic parameters ensured consistent scaling of synaptic inputs, allowing clearer interpretation of plasticity-driven dynamics. When applicable, these values are overridden by those in Table 7 for zone-specific simulations. While the synaptic values chosen fall within biologically plausible ranges, we manually verified that the key model outcomes —including the frequency preference for CSpk generation and PF-PC weight stabilization— were robust to variation in these delay parameters.

PF-PC connectivity

Each PC receives input from 5 PF bundles with different (meta)synaptic weights. These weights are sampled from the uniform distribution and normalized to a sum of 5. This resembles a draw from a scaled symmetric Dirichlet distribution with concentration parameter , i.e., the values are favored to the center of the simplex. The PSC (see PC model) is calculated as:

PF-PC Plasticity

The plasticity mechanism between the PF bundles and PCs consists of an intrinsic homeostatic mechanism and a long-term depression (LTD) mechanism . The homeostatic mechanism is based on the Bienenstock-Cooper-Munro-type (BCM) plasticity mechanism [120]. LTD happens at the IO spike, whereas LTP happens during the absence of IO spikes. The final weight is as shown in Eq. 14. Where, is the PC (where ) and is the synapse of that PC (where ).

BCM

The homeostatic component of synaptic plasticity is governed by a modified BCM rule. This formulation builds on the original BCM model [120] but introduces additional damping to match cerebellar PC firing ranges.

Let denote the effective firing rate (Hz) of the j-th PF bundle, computed from the amplitude of the OU-driven synaptic current, and let be the moving average firing rate (Hz) of the i-th PC over a 200 ms window. The sliding plasticity threshold (Hz) for the i-th PC evolves as shown in Eq. 16.

where, is the threshold time constant. The effective PF firing rate, , is computed by transforming the instantaneous synaptic current (, in nA) into units of Hz by scaling with the amplitude of the OU process and assuming a linear mapping between current magnitude and presynaptic rate. This transformation is valid under the model assumption that PF input represents the activity of a large bundle of asynchronously active fibers whose summed conductance mimics a continuous input current. Moreover, represents the general change in plasticity specific to the i-th PC, as it depends on the firing rate , and the sliding threshold , both of which vary across PCs. This term is then modulated by for each PF bundle:

This classical parabolic term (Eq. 17) has zeros at and a minimum of at . However, left unbounded, this can lead to unrealistically large potentiation at high firing rates. To limit this, we introduced a damping factors:

where is a shape constant. In implementation (Brian2), we ensure the argument is dimensionless by multiplying the rate term by a time unit (second). This hyperbolic term flattens the curve at high firing rates, ensuring that potentiation and depression remain within physiologically plausible bounds.

This formulation yields bidirectional weight changes depending on the position of relative to , with potentiation when and depression when . The tanh function ensures saturation and boundedness of weight change. The only hard constraint is imposed on , limited to 250 Hz to prevent pathological firing, but no explicit bounds are placed on .

To assess the robustness of the BCM mechanism, we varied the threshold time constant, a key parameter controlling the speed at which the PC activity influences synaptic change. The overall behavior of the model remained robust as the PF-PC weights stabilized but the Up/Downbound modules diverged in their plasticity outcomes only for a narrow range of time constant values. Namely, net depression in Downbound modules and net potentiation in Upbound modules. Outside this range, homeostasis was still achieved, but the direction of synaptic change was not aligned with the expected biological behavior. These results suggest that, while the system operates in a resilient homeostatic regime, certain features of learning dynamics are parameter-sensitive and emerge only within a constrained space of biologically plausible plasticity kinetics.

LTD vs LTP

LTD is modeled as a negative change of the synaptic weight proportional to the activity of the PF and with a time decay of 350 ms (Table 7) simulating the calcium release time (Fig 2D). Upon each IO spike, the LTD term was incremented by a constant by a constant is scaled by the absolute deviation from the instantaneous PF input (Hz) from its mean level (converted to Hz by explicitly multiplying the appropriate unit scaling factors). As a consequence, LTP occurs in periods without IO spikes.

As the CSpk-triggered LTD component was incremented by a negative amount and decays exponentially back toward zero between IO spikes, it always remains negative. This decay can allow the BCM-driven component to dominate when its positive contribution exceeds . The total synaptic weight change is the sum of these two contributions.

PC-CN Synapse

The inhibitory synapse from PC to the CN is modeled by a discrete increase in at the time of a presynaptic PC spike. After each increase, decays exponentially back toward zero, resulting in a temporary suppression of CN firing.

Each PC projects to 16 CN cells and each CN cell receives on average from 40 PCs, ranging from 30 to 52 (Fig 1A). Each PC source indices is repeated for each PC and CN target indices are sampled without replacement from the CN population.

CN-IO Synapse

Analogously to the PC-CN synapse, the new synaptic current increases with CN spikes. The decay time of is small (30 ms), and enters the IO at the soma (). This inhibitory input can be seen as the asynchronous release of GABA onto the IO cell [99]. The update value of is seen in Table 6 and the synaptic update is done as follows:

Each CN projects to 10 IO cells and each IO cell receives on average from 10 CN cells, ranging from 6 to 16 (Fig 1A). The connectivity is generated ensuring these ratios. CN cell indices are repeated for each CN cell, while IO indices are sampled without replacement from the IO population, avoiding duplicate connections.

IO-IO Gap Junction

The electrotonic coupling between the cells is represented by a current (see Table 6) in the dendrites and the transjunctional voltage dependence of the gap junction conductance [87] is defined below. In the coupled case, each IO neuron connects to all other neurons in the population (Fig 1A). The transjunctional voltage dependence of the gap junction conductance (Eq. 26) is defined as , where both and are dendritic membrane potentials.

IO-PC Synapse

The modeling of the IO-PC synapse focuses on the pause mechanism following a CSpk. Each PC receives a CF input from one IO as seen in Fig 1A. As the PC AdEX model already has an adaptation variable w that increases the interspike interval (ISI), the synaptic modeling only requires the increase of w following a presynaptic IO spike. The update of the w adaptation variable is increased as shown in Eq. 27. As w is increased to a higher value, the PC cannot spike until w decays to a lower value.

Only half of the IOs synapse with PCs. Each of these IO CF connects to on average 5 PCs (ranging 2–9) and each PC receives from only 1 IO cell (Fig 1A). IO cell indices are sampled without replacement from the IO population. Each PC is either assigned a unique IO source, or —when convergence exceeds the IO population size— assigned an IO neuron randomly selected from a pre-sampled pool, ensuring the intended level of convergence while avoiding duplicate connections.

Upbound and Downbound zones

To reproduce the physiological results of the different zones we did a grid search over 2 synaptic currents, 1 IO conductance and the IO OU current (see Table 7). This grid search was done in steps of 0.1 for all parameters except for the PC to CN synaptic current which was done in steps of 0.001. Then we tried to find a set of parameters such that we could reproduce physiological firing ranges with minimal parameter changes between cases. Table 7 shows the different parameters that need to be changed to reproduce electrophysiological results as shown in Fig 3.

The values of Table 7 override the corresponding baseline parameters listed in Table 6 for the indicated zone and coupling condition. These variations reflect biologically plausible differences in inhibitory strength, dendritic conductance, and IO noise that give rise to module-specific dynamics. Parameters not listed here retain their default values from Table 6.

Filtered inputs

Filtered inputs were obtained via bandpass filtering. For each noise input and frequency range ( to ), a forward 6th butterworth filter was applied on the mean-offset removed original noise signal, obtaining . Filtering was done as cascaded second-order sections, taking into account the rate at which the noise was sampled as with being the time step of the simulation. To account for power loss due to filtering, a correction was applied by multiplying the frequency-domain energy:

CS-US Experiment

We used step currents of different amplitude and duration (Table 8) to simulate the CS and US. The US was copied onto the IO with a 4ms delay and multiplied by 2.8 to have the desired effect on the IO spike response. The ITI was randomly chosen between 8–12 seconds.

Simulations

For all simulations the Euler-Maruyama method was used to simulate the frozen noise. For the rest of the neuron groups, the forward Euler method was used to integrate the differential equations. Two integration timesteps were used, one for the simulations (0.025 ms) and one for the recording of the neuron output (1 ms). All experiments had a simulation length of 120 s. We first started by generating a seed which included generating a network with connectivities, cell parameters and frozen PF input. When modeling the different micromodules we used the parameters shown in Table 7.

The "before plasticity" condition was used to initialize the network and is meant to represent a biologically plausible but non-stationary starting configuration. It does not reflect a steady-state physiological condition. Physiological comparisons should instead be made using the "after plasticity" condition, which represents the stable operating point of the model under continuous adaptation.

The frozen PF input was used for the "no plasticity" experiment. For the subsequent "plasticity" epoch, we enabled the BCM and LTD mechanisms and applied the same frozen PF input. At the end of each plasticity epoch, the evolving PF-PC synaptic weights were smoothed using a Savitzky–Golay filter (window size = 100 samples, polynomial order = 2) to reduce variability caused by the instantaneous state of the LTD re LTP mechanism. Without smoothing, the final weight could be biased if the simulation ended immediately after a CSpk, when the LTD term is transiently active. The final smoothed weight value was then used as the static weight for the subsequent "after plasticity" simulation with the same frozen input. This procedure was repeated for the second, third and fourth plasticity epochs (see Fig 2G-L).

The experiments with filtered signals were done in the same way. The only difference was that the new PF input was a filtered version of the same frozen input. We still had the same network seed. We redid the experiments for the following frequency bands: 5–10 Hz, 10–15 Hz, 15–20 Hz, 25–30 Hz, 50–75 Hz, 100–150 Hz and 1–800 Hz (see Fig 5A).

The CS - US experiment (light/puff) was done by adding a step current of 250 ms to the PF 4 input which represented the CS. The same was done for PF 5 but with a 30 ms step current input that represented the US. This input was also copied (with a 4ms delay) to the IO to elicit CSpks. The two signals co-terminated. A block consisted of an initial US followed by 10 CS-US stimuli with varying ITIs and it finished with a final CS.

For each of the simulations we also did a copy simulation with the uncoupled scenario. This is the same seed, meaning that we had the same network, same connectivity, same cell parameters, and same frozen PF input. The only difference was that the coupling conductance (see Table 3) of the coupling current is lowered to 0 ().

Quantification of firing rates

All reported firing rates (SSpk, CSpk, IO, CN) were computed using 1-second time bins for the entire experiment time. Spike trains were binned, converted to firing rates (Hz), and averaged across all recorded cells in the given population. For statistical comparisons between plasticity conditions, outlier exclusions and variance thresholding were performed using interquartile range filtering (5–95th percentiles), as described in the code and consistent across figures.

Code

All code used for simulations and analysis is available at: https://github.com/eliasmateo95/CerebellarLoop.

Statistical Analysis

To assess differences across simulation conditions, we used two types of parametric statistical tests depending on the data distribution. Two-sample Student’s t-tests were applied for comparisons of means across groups. F-tests were used to compare variances across groups. All tests were two-sided and group sized (n = 40 for CN and IO cells, n = 100 for PCs) are reported in figure captions.

S1 Fig. Gap junctions increase spiking and subthreshold synchrony (A) Gap junctions are disabled by setting the gap junction current = 0.

For the Downbound micromodule, we compare the IO membrane response (B), raster (C), and correlation matrix of membrane potentials (D) in coupled (top) and uncoupled networks (bottom). Synchrony directly impacts the correlations in weight changes expected at the PF-PC synapse.

https://doi.org/10.1371/journal.pcbi.1013609.s001

(TIF)

S2 Fig. Synchrony and Clustering.

In a single micromodule we observe a few clusters in the short time scale whose stability fades exponentially. (A) Kuramoto synchrony metric for STOs is derived from IO Vm in 5 transformation steps. (B) Stability of global oscillator (i.e., usability as a low dimensional global clock): global phase (mean Hilbert phase of all IO cells) from start to end of an 80 second period. The deviation from horizontal indicates phase drifts across phase in radians. Both Upbound and Downbound IO have small drifts and maintain a robust global phase, equating with a single synchronized cluster. (C) Delta phase plot. Difference of phase of each IO neuron to the global phase over time. Plotted modulo 32pi radians to preserve space. Phase-locked oscillators are clearly visible as overlapping lines, while lag/leap events are shown as vertical jumps. Comparing long (D, 30s) and short term (E, 1s) clusterization with the pairwise kuramoto phase reveals a few clusters but a single global oscillator in the 30s frame. (F) Voltage trace of the clusters found in E. The black bars indicate the 1s used for clustering. (G) Cluster stability of 1s clusters. In intervals of 10ms, 1s clusters were found with HDBSCAN [121] and compared with each other via NMI (normalized mutual information). At a distance of 1s, clusters are barely related.

https://doi.org/10.1371/journal.pcbi.1013609.s002

(TIF)

S3 Fig. Uncoupled Synchrony.

When uncoupled, the global oscillator still exhibits stable clock behavior. In Upbound zones there is less phase dispersion, but no phase locking appears in uncoupled scenarios. Downbound zone oscillators act completely independently. (A) IO Spike triggered Kuramoto synchrony with average 25/75 percentiles in Downbound, after plasticity. Average STO synchrony is higher for the coupled case. STO synchrony peaks before and after IO spikes in the coupled case, but not for the uncoupled case. This indicates that gap junctions not only set a base level of synchrony, but they also induce synchronous spikes. Note the increase of synchrony before the IO spike in the coupled case. (B) Stability of the global oscillator (see S2.1B), for the uncoupled case. Mean oscillator phase has smaller variation. (C) Delta phase plot (see S2.1C) for the uncoupled case (without modulo because phase locking does not occur as in the coupled case (S2 Fig). IO neurons in the Upbound zone still keep some synchrony while for Downbound zones, oscillators behave independently.

https://doi.org/10.1371/journal.pcbi.1013609.s003

(TIF)

S4 Fig. Difference between coupled and uncoupled scenarios.

(A) SSpks frequency is significantly different across zones (coupled: t = 13.9, p < 0.001, uncoupled: t = 13.61, p < 0.001; two-sample Student’s t-test; n = 100). (B) CSpks frequency (coupled: t = 3.98, p < 0.001, uncoupled: t = 2.5, p = 0.014; two-sample Student’s t-test; n = 40). The variance of CSpks is significantly different in the uncoupled case across zones(F = 8.65, p < 0.001; F-test using CCDF of F-distribution; n = 40). (C) The CN firing frequency (coupled: t = -3.26, p = 0.002, uncoupled: t = -2.98, p = 0.004; two-sample Student’s t-test; n = 40). (D) The length of the CF Pause is significantly different across zones for both coupling cases (coupled: t = -12.52, p < 0.001, uncoupled: t = -14.39, p = 0.014; two-sample Student’s t-test; n = 40). (E) SSpk CV is significantly different for the Downbound zone across coupling cases (t = 5.94, p < 0.001; two-sample Student’s t-test; n = 100). For the coupled case, across both zones, the SSpk CV is significantly different (t = -4.27, p < 0.001; two-sample Student’s t-test; n = 100). (F) SSpk CV2 is significantly different for both coupled and uncoupled cas, across both zones, (Coupled: t = -13.06, p < 0.001, uncoupled: t = -12.69, p < 0.001; two-sample Student’s t-test; n = 100). (G) CSpk CV is significantly different for the Downbound zone across coupling cases (t = -3.69, p < 0.001; two-sample Student’s t-test; n = 100. For the coupled case, across both zones, the CSpk CV is significantly different (t = 7.16, p < 0.001; two-sample Student’s t-test; n = 100). (H) CSpk CV2 is significantly different for the coupled case across both zones, (t = -3.57, p < 0.001; two-sample Student’s t-test; n = 100). For the uncoupled case, the variance is significantly different across both zones (F = 6.86, p < 0.001; F-test using CCDF of F-distribution; n = 100).

https://doi.org/10.1371/journal.pcbi.1013609.s004

(TIF)

S5 Fig. Increased CN inhibition leads to reduced STO frequency and to more CSpk frequency.

(A) After Plasticity difference for SSpks is significant across coupled cases for the Upbound zone (t = 4.51, p < 0.001; two-sample Student’s t-test; n = 100) and for coupled case across zones (t = -5.36, p < 0.001; two-sample Student’s t-test; n = 100). The variance of SSpk frequency is significantly different for the coupled case across both zones (F = 4.85, p < 0.001; F-test using CCDF of F-distribution; n = 100). (B) CSpk frequency is significantly different across coupling cases (coupled: t = 9.44, p < 0.001, uncoupled: t = 5.84, p < 0.001; two-sample Student’s t-test; n = 40) and the variance is also significantly different (coupled: F = 2.91, p < 0.001, uncoupled: F = 12.19, p < 0.001; F-test using CCDF of F-distribution; n = 40). (C) CN firing frequency is significantly different across coupling cases (coupled: t = 9.15, p < 0.001, uncoupled: t = 8.87, p < 0.001; two-sample Student’s t-test; n = 40). (D) The length of the CF pause is significantly different across coupling cases for Upbound and Downbound zones (Upbound: t = -6.25, p < 0.001, Downbound: t = 2.5, p = 0.013; two-sample Student’s t-test; n = 40). It is also significantly different across zones for the uncoupled case (t = -7.76, p < 0.001; two-sample Student’s t-test; n = 40). (E) Increase in CN-to-IO inhibitory synapses, and thus inhibition strength, led to a decrease in IO STO frequency and a concurrent increase in CSpk frequency.

https://doi.org/10.1371/journal.pcbi.1013609.s005

(TIF)

S6 Fig. Input Filter 5–10Hz.

(A) Firing frequency distributions for each population before and after plasticity for different frequency bands. (B) STO frequency vs IO spike frequency for Upbound (left) and Downbound (right) zones for before (top) and after (bottom) plasticity. There seems to be an inverse relationship between STO frequency and IO spike frequency. To reduce the effect of the spike refractory period on STO frequency, STO frequency is here defined as the mean of the 25–75 percentile of the inverse of time differences between zero crossings in the Hilbert phase. For the 5–10 Hz filter in the Upbound zone, STO frequency seems to decrease while STO spikes increase. For the Downbound zone, STO frequency is lowered after plasticity, correlated to a higher baseline CSpk frequency and thus a less pronounced difference to the 5–10 Hz filter. (C) IO spike triggered PSC for NF (top), 5–10 (middle) and 10–15 Hz (bottom) input filters. In general, IO spikes for The NF case seem to derive from entrainment followed by a final larger decrease in PSC. For 5–10 Hz filter, IO spikes seem to derive purely from entrainment to the oscillatory signal, with also a much clearer (larger amplitude) template before the spike.

https://doi.org/10.1371/journal.pcbi.1013609.s006

(TIF)

S7 Fig. CS-US experiment.

(A) Average weight shift after plasticity for each parallel fiber. The CS pf sees a large decrease in weight. (B) Weight changes during plasticity, averages for each CS/US trial. PF4 (CS) weight decreases continuously during CS. PFs 1,2,3 (noise) decrease during US and return quickly to baseline. PF5(US) responds the least to US. (C) Mean instantaneous Hilbert frequency of STOs, smoothed using 51ms 3rd order Savitzky–Golay filter.

https://doi.org/10.1371/journal.pcbi.1013609.s007

(TIF)

S8 Fig. Uncoupled CS-US experiment.

(A) Kuramoto synchrony for the uncoupled scenario (B) CS-triggered frequency density responses of the populations for SSpk (top) and CSpk (bottom). (C) SSpk frequency response after plasticity is significantly different across zones for both coupling cases (coupled: t = -8.72, p < 0.001, uncoupled: t = -3.02, p = 0.003; two-sample Student’s t-test), as well as the variance (coupled: F = 5.87, p < 0.001, uncoupled: F = 1.94, p < 0.001; F-test using CCDF of F-distribution). (D) CSpk frequency after plasticity is significantly different across coupling cases in the Downbound zone (t = -2.28, p = 0.026; two-sample Student’s t-test). It is also significantly different across zones for both coupling cases (coupled: t = 7.39, p < 0.001, uncoupled: t = -5.51, p < 0.001; two-sample Student’s t-test), as well as the variance (coupled: F = 1.83, p < 0.031, uncoupled: F = 10.4, p < 0.001; F-test using CCDF of F-distribution).

https://doi.org/10.1371/journal.pcbi.1013609.s008

(TIF)

S1 Table. PC SSpk frequency statistics (n = 100).

https://doi.org/10.1371/journal.pcbi.1013609.s009

(DOCX)

S2 Table. CN and IO spike frequency statistics across.

Upbound and Downbound modules

https://doi.org/10.1371/journal.pcbi.1013609.s010

(DOCX)