# A multi-layer mean-field model of the cerebellum embedding microstructure and population-specific dynamics

Roberta Maria Lorenzi, Alice Geminiani, Yann Zerlaut, Marialaura de Grazia, Alain Destexhe, Claudia A.M. Gandini Wheeler-Kingshott, Fulvia Palesi, Claudia Casellato

<figure><img src="image_1.png" /></figure>

## To cite this version:

Roberta Maria Lorenzi, Alice Geminiani, Yann Zerlaut, Marialaura de Grazia, Alain Destexhe, et al. A multi-layer mean-field model of the cerebellum embedding microstructure and population-specific dynamics. PLoS Computational Biology, 2023, 19 (9), pp.e1011434. 10.1371/journal.pcbi.1011434. hal-04281496

HAL Id: hal-04281496

https://cnrs.hal.science/hal-04281496v1

Submitted on 13 Nov 2023

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
RESEARCH ARTICLE

# A multi-layer mean-field model of the cerebellum embedding microstructure and population-specific dynamics

Roberta Maria Lorenzi<sup>1*</sup>, Alice Geminiani<sup>1</sup>, Yann Zerlaut<sup>2</sup>, Marialaura De Grazia<sup>1</sup>, Alain Destexhe<sup>3</sup>, Claudia A. M. Gandini Wheeler-Kingshott<sup>1,4,5</sup>, Fulvia Palesi<sup>1</sup>, Claudia Casellato<sup>1</sup>, Egidio D'Angelo<sup>1,5</sup>

1 Department of Brain and Behavioural Sciences, University of Pavia, Pavia, Italy, 2 Institut du Cerveau—Paris Brain Institute—ICM, Inserm, CNRS, APHP, Hôpital de la Pitié Salpêtrière, Paris, France, 3 Paris-Saclay University, CNRS, Saclay, France, 4 NMR Research Unit, Queen Square Multiple Sclerosis Centre, Department of Neuroinflammation, UCL Queen Square Institute of Neurology, UCL, London, United Kingdom, 5 Brain Connectivity Center, IRCCS Mondino Foundation, Pavia, Italy

* robertamaria.lorenzi01@universitadipavia.it

<figure><img src="image_3.png" /></figure>

## Abstract

Mean-field (MF) models are computational formalism used to summarize in a few statistical parameters the salient biophysical properties of an inter-wired neuronal network. Their formalism normally incorporates different types of neurons and synapses along with their topological organization. MFs are crucial to efficiently implement the computational modules of large-scale models of brain function, maintaining the specificity of local cortical microcircuits. While MFs have been generated for the isocortex, they are still missing for other parts of the brain. Here we have designed and simulated a multi-layer MF of the cerebellar microcircuit (including Granule Cells, Golgi Cells, Molecular Layer Interneurons, and Purkinje Cells) and validated it against experimental data and the corresponding spiking neural network (SNN) microcircuit model. The cerebellar MF was built using a system of equations, where properties of neuronal populations and topological parameters are embedded in inter-dependent transfer functions. The model time constant was optimised using local field potentials recorded experimentally from acute mouse cerebellar slices as a template. The MF reproduced the average dynamics of different neuronal populations in response to various input patterns and predicted the modulation of the Purkinje Cells firing depending on cortical plasticity, which drives learning in associative tasks, and the level of feedforward inhibition. The cerebellar MF provides a computationally efficient tool for future investigations of the causal relationship between microscopic neuronal properties and ensemble brain activity in virtual brain models addressing both physiological and pathological conditions.

<figure><img src="image_4.png" /></figure>

### OPEN ACCESS

**Citation:** Lorenzi RM, Geminiani A, Zerlaut Y, De Grazia M, Destexhe A, Gandini Wheeler-Kingshott CAM, et al. (2023) A multi-layer mean-field model of the cerebellum embedding microstructure and population-specific dynamics. PLoS Comput Biol 19(9): e1011434. https://doi.org/10.1371/journal.pcbi.1011434

**Editor:** Hermann Cuntz, Ernst-Strungmann-Institut, GERMANY

**Received:** December 22, 2022

**Accepted:** August 15, 2023

**Published:** September 1, 2023

**Peer Review History:** PLOS recognizes the benefits of transparency in the peer review process; therefore, we enable the publication of all of the content of peer review and author responses alongside final, published articles. The editorial history of this article is available here: https://doi.org/10.1371/journal.pcbi.1011434

**Copyright:** © 2023 Lorenzi et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** The codes will be available at https://github.com/RobertaMLo/CRBL_MF.

## Author summary

Whole-brain dynamics can be simulated using cortical and subcortical mean-field models, which provide a population-level description of the underlying neuronal dynamics. While mean-field models of the isocortex have recently been developed, a mean-field model of
**Funding:** This research has received funding from the European Union's Horizon 2020 Framework Program for Research and Innovation under the Specific Grant Agreement No. 945539 (Human Brain Project SGA3) to ED, CGWK, FP and AD, and under the Marie Skłodowska-Curie grant agreement No. 892175 to YZ. CGWK received funding from BRC (#BRC704/CAP/CGW), MRC (#MR/S026088/1), Ataxia UK, MS Society (#77), Wings for Life (#169111). This research has also received funding from Centro Fermi project "Local Neuronal Microcircuits" to ED. Special acknowledgement to EBRAINS and FENIX for informatic support and infrastructure. RL, and AG have been supported by Human Brain Project SGA3. This work was also supported by #NEXTGENERATIONEU (NGEU) and funded by the Ministry of University and Research (MUR), National Recovery and Resilience Plan (NRRP), project MNESYS (PE0000006) - A Multiscale integrated approach to the study of the nervous system in health and disease (DN. 1553 11.10.2022) to ED, CGWK, and CC, and Project EBRAINS-Italy (IR00011) - (M4C2 Line 3.1 of the PNRR, Action 3.1.1 - CUP B51E22000150006) to ED and CC. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

**Competing interests:** I have read the journal's policy and the authors of this manuscript have the following competing interests: CGWK is a shareholder in Queen Square Analytics Ltd.

the cerebellar cortex is still missing but is much needed given its specific structural and functional organization. Thus, we developed the first biologically grounded mean-field model of the cerebellar cortex, which embeds a realistic network architecture with 4 main neuron populations (granule cells, Golgi cells, Purkinje cells, molecular layer interneurons) represented with non-linear neuronal models embedding a set of neuron- and synapse-specific parameters. The model was validated and tuned against experimental data and spiking neural network simulations. The mean-field model can reproduce local neural dynamics elicited by different cortical inputs and accurately predicts population-specific activity patterns. The possibility of tuning multiple neuronal and synaptic parameters allows to capture local neural dynamics both in physiological and pathological conditions. The cerebellar mean-field model is now ready to be integrated into brain dynamic simulators, fostering a deeper understanding of the cerebellar impact on brain dynamics in functional and dysfunctional states.

# 1. Introduction

Realistic modelling of brain function is opening new frontiers for experimental and clinical research towards personalized and precision medicine [1,2]. Brain models can be developed at different scales, ranging from microscopic properties of neurons and microcircuits to the ensemble behavior of the whole brain. Arguably, a model able to reproduce realistic behavior of specific neuronal populations across scales would increase the fidelity of modelling single brain regions and consequently improving the accuracy of whole-brain dynamic simulations [3]. This approach, though, bears conceptual and practical drawbacks. At the microscale level, e.g. at the µm³ scale, Spiking Neural Networks (SNNs) reproduce neural circuits as a set of interconnected neurons [4-6]: the state of each neuron and synapse in the network is updated at each simulation step, allowing the investigation of neural circuits functioning with a high level of granularity and biological plausibility. However, this degree of detail is hard to manage when simulating brain signals, like those derived from electroencephalography (EEG) or functional magnetic resonance imaging (fMRI), representing the activity of large cortical regions, e.g., at the mm³ scale. To manage the high complexity of brain signals, the dynamics of neuronal populations have been condensed into ensemble density models called neural masses. These provide a description of the expected values of neuronal activity states, under the assumption that the equilibrium density has a point mass [7,8]. Neural fields are obtained from neural mass models when considering spatial information: these can be used to model spatial propagation of activity throughout brain volumes [9]. Despite being computationally efficient and easy to fit on brain signal data, neural mass and neural field models lack a direct link to the real behavior of neurons at the microscopic scale; this limits their applicability to investigate the neuronal bases of brain dynamics and the causal relationships between neural mechanisms at different scales. The mean-field (MF) approximation theory describes collective neural activities maintaining a direct link with the corresponding spiking neural network. The MF model formalism approximates high-dimensional random systems by averaging their original properties over degrees of freedom to compute the first two statistical moments of the system variables (e.g., mean and variance of the population firing rates in a neural network). Different mathematical approaches exist to derive such statistical descriptions depending on the complexity of the underlying network. While theoretical derivations in simpler settings allow to capture firing dynamics up to fast synchronized events (e.g., Lorentzian [10,11] or finite-size effect-based formalisms [12,13]), the network
considered here impedes such an approach. Instead, we use a formalism based on the defi-
nition of a transfer function (TF) that models the neuronal input-output relationship to
investigate the properties of network responses to different input [7,14,15]. The MF devel-
oped here is based on the TF approach that summarizes the neuronal and connectivity
properties of an entire SNN through ad-hoc transfer functions [15-17] and captures the sta-
tistical properties of network activity by computing the probabilistic evolution of neuronal
states at subsequent time intervals [14,18,19].

MFs can be used to investigate macroscale phenomena, such as brain rhythms and coherent oscillations [20], and are computationally advantageous, with increased computational speed and low memory requirements compared to SNNs. Among current limitations, this MF does not capture in full the complex properties of specific neuronal populations and are valid only in certain firing regimes, e.g., at low frequency [21].

Moreover, while a diversification of MFs for specific cortical regions has been proposed [22–28], the development of MFs for subcortical regions has seen only a few attempts [29–32], despite the fundamental role such regions exert in controlling brain dynamics and behavior [33–38]. In particular, the cerebellum has a dense connectivity with the thalamus and the cerebral cortex and impacts remarkably on whole-brain dynamics, as shown in resting-state and task-dependent fMRI [36,39,40]. The cerebellum highly specific cortical microcircuit structure and its neuronal population specific dynamics are prompting for the development of specific MFs to be included into whole-brain simulators [41].

The cerebellar cortex, indeed, receives inputs from mossy fibers and climbing fibers and
sends outputs to the deep cerebellar nuclei. Granule Cells (GrC), Golgi Cells (GoC), Molecular
Layer Interneuron (MLI) and Purkinje Cells (PC) constitute the backbone of the cerebellar
cortex, which shows a unique anisotropic geometry supporting a forward architecture with
limited lateral connectivity, recurrent excitation and inhibition, and an inhibitory output pro-
jecting from PC to the Deep Cerebellar Nuclei. These properties, along with the above-men-
tioned neuronal types, differ remarkably from those of the cerebral cortex, which is typically
modelled using MFs including one excitatory and one inhibitory population that are recip-
rocally connected, have recurrent excitation and project an excitatory output. This general
MF design neglects specific properties of the cerebellum such as the multilayer organization
and the lack of recurrent excitation, justifying the need of a specific cerebellar MF model
with appropriate circuitry and temporal dynamics. All the types of neurons of the cerebellar
cortex have been characterized through electrophysiology experiments in *rodents in vitro*
and *in vivo* [42–46] and have been, subsequently, represented by detailed multicompart-
mental models [47–51]; moreover, these detailed models have been simplified into point-
neuron models [52–54], and embedded in network models of the cerebellar microcircuit
[44,52,55–57]. Thus, the cerebellum provides an ideal substrate for generating a MF, with
internal dynamics that can be remapped onto a precise physiological counterpart. Finally,
the MF model of the cerebellum can be validated against the existing rich electrophysiologi-
cal data and detailed cerebellar microcircuit activity simulations. In this work we have
developed and validated a multi-layer MF of the cerebellar cortex, which maintains the
salient properties of the inter-wired cerebellar neuronal populations. Indeed, the MF was
derived from a biology-grounded model of the cerebellar microcircuit, which was used to
define the topology and tune the parameters of the MF, and it was then validated against a
rich set of SNN outputs. The future prospective is to integrate the present mesoscopic cere-
bellar MF as a module of macroscale models, e.g. Dynamic Causal Modeling (DCM)
[58–60] or The Virtual Brain (TVB) [41], to simulate the cerebellar contribution to brain
activity in physiological and pathological conditions.
## 2. Methods

In this section we describe the development, tuning, validation, and application of a multi-layer MF of the cerebellar circuit (Fig 1). The MF formalism provides a statistical summary of a SNN activity through the first two statistical moments (i.e. average and variance) of the population firing rates [14]. Here the SNN bottom-up modelling approach is merged with the standard MF mathematical formalism to obtain a multi-layer MF of the cerebellar cortex. In order, we present the cerebellar SNN model used as the structural and functional reference of the MF (section 2.1), the design of the MF architecture based on cerebellar topology (section 2.2.1 and 2.2.3), the implementation of the MF equations derived with an heuristic approach (section 2.2.2) [19,21,61] the protocols used to optimise the MF time constant (section 2.2.4) and to validate the MF (section 2.3), the applications of the MF to predict the activity modulation induced by different levels of synaptic plasticity and of inhibitory control (section 2.4). For the sake of simplicity, activity patterns in the MF are described using a terminology derived from the SNN, so that a rate-coded waveform with a peak followed by a reduction is called burst-pause to suggest that the generative mechanism is the same as that of a spike-coded pattern.

### 2.1 SNN model

This cerebellar cortex model was built using the Brain Scaffold Builder (BSB) (https://bsb.readthedocs.io/en/latest/), a neuroinformatic framework allowing a detailed microcircuit

<figure><img src="image_5.png" /><figcaption>**Fig 1. Pipeline of the multi-layer cerebellar Mean Field (MF) model.** The workflow of the study is represented. MF was designed based on structural and functional parameters extracted from Spiking Neural Network (SNN) simulations. The time constant of the resulting MF was optimized against Local Field Potential (LFP) experimental data. The model was first validated against neural activity of SNN with different stimulation protocols (Constructive & Functional validity) to compare MF simulations with SNN outcomes. MF was also tested to reproduce the effect of synaptic plasticity (Predictive Validity) in Purkinje Cells (PC) by modulating the parallel fibers (pf) and molecular layer interneurons (MLI) synaptic strength to investigate respectively learning mechanisms and feed forward inhibition.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1011434.g001
reconstruction based on neuron morphologies and orientations and the incorporation of active neuronal and synaptic properties [56]. The model was simulated as a SNN with ~3x10⁴ extended-Generalised Leaky Integrate and Fire (E-GLIF) neurons [52,54] and ~1.5x10⁶ alpha-shaped conductance-based synapses [62]. The SNN simulations were performed using NEST version 2.18 (https://zenodo.org/record/2605422) [4,63].

**2.1.1 Neuron model.** The E-GLIF formalism describes the time evolution of membrane potential ($V_m$) depending on two intrinsic currents to generate slow adaptation ($I_{adap}$) and fast depolarisation ($I_{dep}$), using a system of three Ordinary Differential Equations [54]

$$ 
\left\{
\begin{aligned}
\frac{dV_m(t)}{dt} &= \frac{1}{C_m} (C_m(V_m(t) - E_{rev}) - I_{adap}(t) + I_{dep}(t) + I_e + I_{syn}) \\
\frac{dI_{adap}(t)}{dt} &= k_{adap}(V_m(t) - E_{rev}) - k_2 I_{dep}(t) \\
\frac{dI_{dep}(t)}{dt} &= k_1 I_{dep}(t)
\end{aligned}
\right.
\quad (1)
$$

where $I_{syn}$ = synaptic current (it models the synaptic stimulus, see section 2.1.2); $C_m$ = membrane capacitance; $\tau_m$ = membrane time constant; $E_{rev}$ = reversal potential; $I_e$ = endogenous current; $k_{adap}$ and $k_2$ = adaptation constants; $k_1$ = decay rate of $I_{dep}$. When a spike occurs, state variables are updated as follows:

$$ 
\begin{aligned}
V_m(t_{spk}^+) &= V_r \\
I_{adap}(t_{spk}^+) &= I_{adap}t_{spk}^+ + A_2 \\
I_{dep}(t_{spk}^+) &= A_1
\end{aligned}
\quad (2)
$$

where $t_{spk}^+$ = time instant immediately following the spike time $t_{spk}$; $V_r$ = reset potential; $A_2$, $A_1$ = model currents update constants. E-GLIF models were implemented using sets of parameters specific for each neuronal population [52] as shown in S1 Table.

**2.1.2 Synaptic model.** Connections between neural populations were modelled as conductance-based synapses:

$$ I_{syn}(t) = g_{syn}(t)(V_m(t) - E_{rev}) \quad (3) $$

When a spike occurs, the conductance $g_{syn}$ changes according to an alpha function:

$$ g_{syn}(t) = G_{syn} \frac{t - t_{spk}}{\tau_{syn}} e^{\frac{t - t_{spk}}{\tau_{syn}}} \quad (4) $$

where $G_{syn}$ is the maximum conductance change and $\tau_{syn}$ the synaptic time constant. E-GLIF neuron models and conductance-based synaptic models used in SNN simulations provided the functional reference of cerebellar spiking activity for MF development.

## 2.2 Mean Field (MF) design

The design of the cerebellar multi-layer MF was based on the extensive knowledge of cerebellar anatomy and physiology summarized in previous cerebellar cortex network models [52,54,56].

**2.2.1 Architecture.** The cerebellar MF included the main neuronal populations of the cerebellar cortex—GrC, GoC, MLI and PC (Fig 2) and the corresponding excitatory and inhibitory synapses. The granular layer at the cerebellar input stage includes GrC and GoC receiving the external input ($v_{drive}$) from mossy fibers. GrC excite GoC, which, in turn, inhibit themselves and GrC forming recurrent loops. GrC represent the excitatory input for the molecular
<figure><img src="image_6.png" /><figcaption>**Fig 2.** Multi-layer mean-field architecture and parameters. **A:** Spiking Neural Network model (SNN). The cerebellar cortical volume (length x width x height = 300 x = 200 x295 µm3) contained a total of 29230 neurons including 28615 Granule Cells (GrC), 70 Golgi Cells (GoC), 446 Molecular Layer Interneurons (MLI) and 99 Purkinje Cells (PC). **B:** Multi-layer MF architecture with neuronal populations connected according to anatomical knowledge. The main cerebellar neuron types are included: GrC and GoC, receiving input from mossy fibers (mf), MLI, and PC which are the sole output of the cerebellar cortex sending their projections to the Deep Cerebellar Nuclei (DCN). Each population receives excitatory and/or inhibitory input activity) from presynaptic populations, depending on their specific conductance µ and on the synaptic properties of each connection (K = synaptic mean convergence, Q = synaptic conductance, and τ = synaptic time constant).</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1011434.g002

layer constituted by MLI and PC. MLI inhibit PC, which are the sole output of the cerebellar cortex and shape the deep cerebellar nuclei activity through inhibition. Although other neurons have been reported to play a role in the cerebellar microcircuit (e.g. Lugaro cells [64], and unipolar brush cells [65]), for the sake of simplicity we have limited the present model to the canonical architecture that is thought to generate the core network computations.

To connect the nodes of the MF network, synaptic parameters were set according to those of the reference SNN (S2 Table). The connection probability of each connection type (K) was derived from the synaptic convergence on the postsynaptic neuronal population in a cerebellar cortical volume [56]. The quantal synaptic conductance and synaptic time decay (Q, τ) were derived from the weights and time constants of the corresponding synapse models [52].

**2.2.2 Transfer Function (TF) Computation.** The TF is defined for each population as a mathematical construct that takes the activity of the presynaptic population ($v_s$) as input and provides an average population activity signal as output ($v_{out}$) [18].

A purely analytic derivation of the TF (using approximations and stochastic calculus, e.g. [15]) was not possible given the complexity of the neuronal (E-GLIF) and synaptic (alpha-waveform) models considered here. Therefore, TFs, here, relied on a semi-analytical approach that couples an approximate analytical estimate with an optimization step to capture the firing response of analytically intractable models [19], see also [16] for a similar approach). More details can be found in [19,61] but we summarize the approach below.

The analytical template for the TF (indicated with F in the equations for sake of simplicity) of all neuron types is derived from the probability to be above threshold in the fluctuation-
driven regime [18]:

$$v_{out} = F_p(v_s) = \alpha \frac{1}{2\tau_v} \text{erfc}\left(\frac{V_{thre}^{\text{eff}} - \mu_v}{\sqrt{2}\sigma_v}\right) \quad (5)$$

where the complementary erfc is the error function while $\mu_v$, $\sigma_v^2$ and $\tau_v$ are the average, variance and autocorrelation time respectively of the membrane potential fluctuations. Two phenomenological terms were introduced: $V_{thre}^{\text{eff}}$, an effective firing threshold to capture the impact of single cell non-linearities on firing response [19] and $\alpha$, a multiplicative factor to adapt the equations also to high input frequency regimes [21]. Those two terms were optimized for each neuron type (see steps d and e, below) from single neuron simulations of input-output transformation in terms of firing rate (i.e., the numerical TF, see step b). The TF depends on the statistical properties of the subthreshold membrane voltage dynamics (mean = $\mu_v$, standard deviation = $\sigma_v^2$ and autocorrelation time $\tau_v$, calculated in step c). These in turns depend on the average population conductances that are computed with the biologically-grounded functional parameters derived from SNN models at single neuron resolution (step a), bringing the physiological properties into the MF mathematical construct.

a) Equations of Population-specific conductance. For each neuronal population, the average conductance was defined as a function of the presynaptic inputs, according to the topology described in section 2.1 (Fig 3.):

$$\mu_p = \sum_s K_{s-p} \tau_{s-p} Q_{s-p} v_s \quad (6)$$

where, for each population $p$ ($p = \text{GrC}$, $\text{GoC}$, $\text{MLI}$, $\text{PC}$), $K_{s-p}$, $\tau_{s-p}$, $Q_{s-p}$ are the connection probabilities, synaptic decay times and quantal conductances of the connection for each presynaptic population $s$ (e.g., for $p = \text{GrC}$, $s = \text{mf}$ or $\text{GoC}$, hence $s-p$ is either $\text{mf-GrC}$ or $\text{GoC-GrC}$), $v_s$ is the presynaptic population activity in Hz computed as explained in (b).

b) Numerical TF. The reference functional target to reproduce with the MF model was the neuronal spiking activity obtained in SNN simulations (in vivo conditions) as described in section 2.1. The activities of GrC, GoC, MLI and PC embedded in the SNN were simulated for different input amplitudes ($v_{drive}$) with input spikes generated from a Poisson distribution at a rate in the range 0–80 Hz, i.e. within physiological firing rate values of mossy fiber activity at rest and during tasks [66]. For each $v_{drive}$, the simulation lasted 5 seconds with a time resolution of 0.1 ms. The working frequencies of each population were extracted by averaging the spiking neuron firing rates. For each population, the outcome of the numerical TF computation was a template of dimension equal to the number of presynaptic populations, resulting in 2D numerical TFs for GrC, MLI and PC and 3D numerical TF for GoC (Fig 3A).

The 2D numerical TF of each population was computed as the population firing rate when receiving the firing rates of the presynaptic populations, given a certain $v_{drive}$: for example, for the numerical TF template of PC, the average firing rates of MLI and GrC were computed for each $v_{drive}$ in the range 0–80 Hz. Then, these quantities were used as presynaptic signals to stimulate PC and a numerical template was obtained from the resulting PC firing rate, for each combination of presynaptic activities. The 3D TF numerical template of GoC was computed following the same strategy but considering 3 presynaptic signals (GrC, GoC, mossy fibers). GrC excitation and GoC self-inhibition were extracted from SNN simulations and the mossy fibers excitation corresponded to $v_{drive}$.

c) Statistical moments of the MF. The statistical moments included in the MF are $\mu_v$, $\sigma_v$ and $\tau_v$. Starting from the conductances of the presynaptic populations, the average conductance
<figure><img src="image_7.png" /><figcaption>B) Analytical TF</figcaption></figure>

<figure><img src="image_8.png" /></figure>
**Fig 3. Numerical Transfer Functions (TFs) and the corresponding Analytical fitting.** The simulations used to compute the numerical TFs last 5 seconds, with a time step = 0.1 ms. A) A 2D numerical TF template is reported for Granule Cells (GrC), Molecular Layer Interneurons (MLI) and Purkinje Cells (PC), which receive inputs from two presynaptic populations (x axis: excitatory input; y axis: inhibitory input; color bar: numerical TF color code). A 3D numerical TF template is reported for Golgi Cells (GoC), which receive input from 3 presynaptic sources, i.e., mossy fibers, GrC and GoC (showed on the axis, color bar: GoC numerical TF color code). From the 3D domain of frequencies combination, only the physiological working frequencies of GrC and GoC are considered for mossy fibers inputs from 0 Hz to 80 Hz, as obtained in corresponding spiking neural network simulations and 3D numerical TF for GoC which receive three presynaptic inputs. B) Numerical TFs (points) are used to fit the corresponding analytical TFs (lines, v<sub>out</sub> on y axis). Excitatory inputs are reported on the x axis, inhibitory ones are color-coded (color bar). The 2D analytical TF presents a non-linear trend, while the GoC Analytical TF shape is almost linear both for lower v<sub>drive</sub> = 20 Hz and higher v<sub>drive</sub> = 80 Hz.

https://doi.org/10.1371/journal.pcbi.1011434.g003

$\mu_G$ of the target population reads:

$$ \mu_G(v_s) = \sum_s \mu_{Gs} + g_L \qquad (7) $$

Where $\mu_{Gs}$ is the presynaptic population conductance (Eq 4) and $g_L$ is the leak conductance of the target population (Table 1). Then, the effective membrane time constant of the target population, $\tau_m^{eff}$, is computed from $\mu_G$ as:

$$ \tau_m^{eff}(v_s) = \frac{C_m}{\mu_G} \qquad (8) $$

Where $C_m$ is the membrane capacitance (Table 1). The first statistical moment, i.e., the average of the membrane potential fluctuation, $\mu_V$, reads:

$$ \mu_V(v_s) = e^{\frac{\sum_s \mu_{Gs} E_s + g_L E_L}{\mu_G}} \qquad (9) $$

With $E_s$ is the reversal potential of the presynaptic connection (0 mV for presynaptic excitatory populations and -80 mV for presynaptic inhibitory populations), $E_L$ is the rest potential of the target population (Table 1).

This expression is adapted from [61] to model the alpha synapses consistently with the models used in the SNN (Eqs 1, 2, and 3). Consequently, the variance and the autocorrelation

<caption>**Table 1. Fitted coefficients of the Analytical TFs.**</caption>

Phenomenological threshold.

$$ V_{thre}^{eff}(\mu_V, \sigma_V, \tau_V) = P_0 + P_{\mu V} \frac{\mu_V - \mu_V^0}{\mu_V^0} + P_{\sigma V} \frac{\sigma_V - \sigma_V^0}{\sigma_V^0} + P_{\tau V} \frac{\tau_V^N - \tau_V^{N0}}{\tau_V^{N0}} + P_{\mu G} \ln\left(\frac{\mu_G}{g_L}\right) $$

Fitted coefficient P [V]

<table><thead><tr><th></th><th>P<sub>0</sub></th><th>P<sub>μV</sub></th><th>P<sub>σV</sub></th><th>P<sub>τV</sub></th><th>P<sub>μG</sub></th></tr></thead><tbody><tr><th>GrC</th><td>-0.426</td><td>0.007</td><td>0.023</td><td>0.482</td><td>0.216</td></tr><tr><th>GoC</th><td>-0.144</td><td>0.003</td><td>0.011</td><td>0.031</td><td>0.011</td></tr><tr><th>MLI</th><td>-0.128</td><td>-0.001</td><td>0.012</td><td>-0.093</td><td>-0.063</td></tr><tr><th>PC</th><td>-0.080</td><td>0.009</td><td>0.004</td><td>0.006</td><td>0.014</td></tr></tbody></table>

P coefficients computed with the fitting procedure explained in [19] and extended to E-GLIF neurons

https://doi.org/10.1371/journal.pcbi.1011434.t001
time of membrane fluctuations result in:

$$ \sigma_V(v_s) = \sqrt{\sum_s (2\tau_m^{\text{eff}} + \tau_{s-p}) \left( \frac{eU_{s-p}\tau_{s-p}}{2(\tau_m^{\text{eff}} + \tau_{s-p})} \right)^2 K_{s-p}v_{s-p}} \quad (10) $$

$$ \tau_V(v_s) = \frac{1}{2} \frac{\sum_s K_{p-s} v_s (e U_{p-s} \tau_{p-s})^2}{\sum_s \left( (2\tau_m^{\text{eff}} + \tau_{p-s}) \left( \frac{e U_{p-s} \tau_{p-s}}{2(\tau_{p-s} + \tau_m^{\text{eff}})} \right)^2 K_{p-s} v_s \right)} \quad (11) $$

$$ \text{with } U_s = \frac{Q_{s-p}}{\mu_G} (E_s - \mu_V). $$

d) Phenomenological threshold. The ability of the analytical template (Eq 5) to capture different firing behavior is given by the introduction of 5 parameters in the phenomenological threshold term. The phenomenological threshold is expressed as a linear combination of the $V_m$ fluctuations properties whose coefficients are linearly fitted to the numerical TF data [19]:

$$ V_{\text{thre}}^{\text{eff}}(\mu_V, \sigma_V, \tau_V) = P_0 + P_{\mu V} \frac{\mu_V - \mu_V^0}{\partial \mu_V^0} + P_{\sigma V} \frac{\sigma_V - \sigma_V^0}{\partial \sigma_V^0} + P_{\tau V} \frac{\tau_V^N - \tau_V^{N0}}{\partial \tau_V^{N0}} + P_{\mu G} \ln \left( \frac{\mu_G}{g_L} \right) \quad (12) $$

Where $\tau_V^N$ is $\tau_V$ adjusted with the ratio between membrane capacitance and leak conductance ($C_m$), and $\mu_V^0$, $\sigma_V^0$, $\tau_V^{N0}$, $\partial\mu_V^0$, $\partial\sigma_V^0$, $\partial\tau_V^{N0}$ are rescaling constants to normalize the contribution of each term [61]. P are the polynomial coefficients which are the target of the fitting procedure to compute the analytical TF as explained in e) (see Table 1).

e) Analytical TF. The statistical moments in Eqs 9, 10, 11 and the phenomenological threshold in Eq 12 were plugged into Eq 13 and the phenomenological threshold is computed through a fitting procedure described in [19]. Specifically, here the semi-analytical TFs are the output of two minimization algorithm with a tolerance of 10e-5: (i) Sequential Least Squares Programming algorithm and (ii) the Nelder-Mead was used to find out the optimal coefficient of the phenomenological threshold (scipy.optimize.minimize library). The TFs specific for the cerebellar populations are reported in Fig 3B.

The parameter alpha (Eq 5) was set to an optimal value for each population to fit both low and high frequencies [21]. The analytical TF, together with the statistical moments $\mu_V$, $\sigma_V$, and $\tau_V$ defined the cerebellar MF equations.

### 2.2.3 Multi-layer equations.
The multi-layer MF was developed as a set of equations capturing the interdependence of the population-specific TFs, tailoring the isocortical MF described in [14] for excitatory-inhibitory networks to the cerebellar network. This formalism describes the network activity at a time resolution T which is set to ensure a Markovian dynamic of the network: T should be large enough to ensure memoryless activity (e.g., it cannot be much lower than the refractory period, which would introduce memory effects) and small enough so that each neuron fires statistically less than once per time-bin T. The choice of T is quite crucial and here it was tailored to account for cerebellar dynamics as explained in section 2.2.4.

The model describes the dynamics of the first and the second moments of the population activity for each population. The cerebellar network was built up with four interconnected populations (GrC, GoC, MLI, PC) receiving external input from mossy fibers (mf) (Fig 2), thus resulting in twenty differential equations: the four population activities ($v_{GrC}(t)$, $v_{GoC}(t)$, $v_{MLI}(t)$, $v_{PC}(t)$) and the driving input ($v_{mf}(t) = v_{drive}(t)$), the four variances of the population
activities ($c_{GrC-GrC}(t)$, $c_{GoC-GoC}(t)$, $c_{MLI-MLI}(t)$, $c_{PC-PC}(t)$) and the one of the driving input from mossy fibers ($c_{mf-mf}(t)$), the six covariances among population activities ($c_{GrC-GoC}(t)$, $c_{GrC-PC}(t)$, $c_{GrC-MLI}(t)$, $c_{GoC-MLI}(t)$, $c_{GoC-PC}(t)$, $c_{MLI-PC}(t)$) and the four covariances between population activities and the driving input ($c_{GrC-mf}(t)$, $c_{GoC-mf}(t)$, $c_{MLI-mf}(t)$, $c_{PC-mf}(t)$). Einstein summation notation was used to report the differential system in a concise form:

$$ T \frac{dv_{\mu}}{dt} = (F_{\mu} - v_{\mu}) + \frac{1}{2} c_{\lambda\eta} \frac{\partial F_{\mu}}{\partial v_{\lambda} \partial v_{\eta}} \quad (13) \\ T \frac{dc_{\lambda\eta}}{dt} = \delta_{\lambda\eta} \frac{F_{\lambda}(\frac{1}{T} - F_{\eta})}{N_{\lambda}} + (F_{\lambda} - v_{\lambda})(F_{\eta} - v_{\eta}) + \frac{\partial F_{\lambda}}{\partial v_{\mu}} c_{\eta\mu} + \frac{\partial F_{\eta}}{\partial v_{\mu}} c_{\lambda\mu} - 2c_{\lambda\eta} $$

Where $v_{\mu}$ is the activity of population $\mu$; $c_{\lambda\eta}$ is the (co)variance between population $\lambda$ and $\eta$; $N_{\lambda}$ is the number of cells included in population $\lambda$. TF dependencies on the firing rate of presynaptic populations are omitted for simplicity of notation, yielding $F_{\mu}$ instead of $F_{\mu}(v_{s})$ with $\mu = \{GrC, GoC, MLI, PC\}$ and $s$ = presynaptic population (i.e. for example, for $\mu = GoC$, then $F_{GoC} = F_{GoC}(v_{drive}, v_{GrC}, v_{GoC})$. The model equations [13] were numerically solved using forward Euler method with an integration step of 0.1 ms. All the parameters are listed in S3 Table.

### 2.2.4 Timing optimisation.
The MF time-constant T was optimised by comparing the model prediction with experimental data of the cerebellar Granular layer (Fig 4). The simulated average activity was interpolated with the experimental Local Field Potential (LFP) measured with high-density microelectrode arrays (HD-MEA) in the granular layer of acute mouse cerebellar slices [67].

LFP data were recorded at 37°C. The external stimulus consisted in a pulse train of 5 stimuli of 50 Hz amplitude. This stimulation protocol was repeated nine times changing the HD-MEA recording channels across each experiment (Fig 4A). The LFP signals recorded were averaged across the nine experiments resulting in five values that represented the average of each pulse of the input trains. These average records were normalised on the amplitude of the signal recorded after the first stimulus [68].

The cerebellar MF simulation protocol was configured with a $v_{drive} = 50$ Hz for 100 ms, reproducing the experimental protocol (Fig 4B). A range of plausible T values were evaluated according to literature [14,19,21,61]. MF simulations were performed with a systematic change of T values and the granular layer average activity was calculated by a weighted-mean of GrC and GoC activities. The weight of GrC and of GoC was computed according to the ratio of the spiking surfaces (GoC/GrC), which resulted in GoC/GrC = 0.13 [67]. The granular layer average activity was normalised on the maximum peak of the simulated activity, and it was interpolated with LFP recordings (Fig 4C) normalized on the maximum peak of recorded activity. For each simulation the mean absolute error between MF simulation and LFP data was computed to select the T value that minimised the discrepancy between MF and LFP signals. Since the granular layer is the driving layer of the network, the optimal T value obtained for this population was assumed to be the same for the molecular and Purkinje layers.

## 2.3 Constructive and functional validity

For constructive and *functional validity*, the cerebellar MF was tested using stimulation protocols designed to assess its ability in reproducing proper cerebellar dynamics and stimulus-response patterns.

Four different stimulation protocols were defined to reproduce the encoding of stimuli in the mossy fibers, each lasting 500 ms:
<figure><img src="image_9.png" /><figcaption>**Fig 4.** The mean-field time constant. A) Experimental acquisition of LFP. LFP signals were acquired in the cerebellar granular layer of acute mice cerebellar para-sagittal slices using HD-MEA in response to five stimulation pulse trains of 50 Hz. Figure created with Biorender.com. B) MF simulation. The activity of the granular layer was simulated with the cerebellar MF using a stimulation protocol emulating experimental LFP recordings. C) Interpolation of LFP and MF. The weighted average of the predicted Granular Layer activity($v_{GRL}$, violet line) interpolates the LFP data (mean ± SD; dots and bars). The relative weights of GoC and GrC are 13% and 87%, respectively. The optimal T value is 3.5 ms ± 5%, (mean ± mean absolute error between $v_{GRL}$ and LFP).</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1011434.g004

i. $v_{drive}$ = *Step function*. A square wave with steps of amplitude 50 Hz, and lasting 250 ms to reproduce a stimulus like a sound [57].

ii. $v_{drive}$ = *Sinusoidal input*. A sinusoidal input of 20 Hz amplitude and oscillation at 6 Hz to simulate the whisker movements experimental conditions [69].

iii. $v_{drive}$ = *Cortical-like input*. A combination of 3 sinusoidal inputs, with fixed amplitude at 40 Hz and oscillations at 1Hz, 15Hz and 30Hz respectively reproducing an EEG-like pattern.

iv. $v_{drive}$ = *Sensory-like input*. Summation of the step function described in (i) with amplitude of 20 Hz and sinusoidal function with an amplitude of 7.5 Hz and oscillating as described in (iii) to simulate a more complex input like a stimulus overlapped to a realistic basal activity.

Population activities with their standard deviations predicted by the MF were overlayed to the Peristimulus Time Histogram (PSTH) with bin = 15 ms, computed from spiking activities of the SNN in NEST simulations of the same stimulation protocols. Then, the PC activity, i.e., the output of the cerebellar cortex, was quantified as mean ± standard deviation of the firing
rate for both MF and SNN simulations. Boxplots were generated, and the Root Mean Squared Error (RMSE) was computed to quantitatively compare the output activity of MF with SNN. Correlation matrices between population activity were computed both for MF and SNN (See S1 Fig). The computational efficiency of each model was measured as computational time in seconds required for each simulation to be performed.

## 2.4 Predictive validity

For predictive validity, MF parameters were tuned to explore the MF sensitivity to modifications of local mechanisms. These modifications were derived from experimental studies on neural correlates of behavior in functional or dysfunctional conditions, focusing on inhibitory control and long term plasticity on PCs [70,71].

**2.4.1 MLI feed forward inhibition modulation.** Feedforward inhibition from molecular layer interneurons regulates adaptation of PCs. Impact of MLI-PC conductance on PC activity was explored by defining different values of MLI-PC synaptic strength $w_{MLI-PC}$ where $w_{MLI-PC} = 100\%$ represents the standard condition, rates lower than 100% model disinhibited activity, while rates higher than 100% extra- inhibition.

$w_{MLI-PC}$ was added to the Analytical $TF_{PC}$ as a modulatory parameter of the presynaptic input $v_{MLI}$, resulting in a modulation of MLI contribution in PC population conductance [35] defined as follows:

$$ \mu_{PC} = K_{MLI-PC} Q_{MLI-PC} \tau_{MLI-PC} v_{MLI} w_{MLI-PC} + K_{GrC-PC} Q_{GrC-PC} \tau_{GrC-PC} v_{GrC} \quad (14) $$

Each simulation lasted 500 ms with a 50Hz driving input of 250 ms after 125 ms of resting.

The Area Under Curve (AUC) of PC activity, PC peak and the depth of the pause were computed as quantitative scores for each $w_{MLI-PC}$ value AUCs and PC peaks were normalised on the respective values corresponding to the standard condition defined as $w_{MLI-PC} = 100\%$.

**2.4.2 PC Learning.** Long term potentiation and depression (LTP and LTD) are forms of synaptic plasticity at the basis of brain learning processes [72]. In the cerebellum, motor learning is driven by PC activity modulation, regulated by the plasticity of the synapses between parallel fibers (projecting from GrC) and PC, resulting in a reduction of PC activity due to LTD and in an increase of PC activity due to LTP [73,74]. To simulate the effects of GrC-PC plasticity on PC dynamics, different synaptic strengths ($w_{GrC-PC}$) were explored, including $w_{GrC-PC} = 65\%$ that corresponds to the decrease of GrC-PC synaptic strength during motor learning according to animal experiment [70]. The values smaller than 100% means LTD occurred, while the values higher than 100% represent LTP occurrence. The strategy applied was analogous to equation (14).

PC AUC and PC peaks were computed for each $w_{GrC-PC}$ value as quantitative score to correlate the amount of LTD and LTP with the output activity of the cerebellar cortex.

## 2.5 Hardware and software

The SNN was built with the BSB release 3.0 (https://bsb.readthedocs.io/en/v3-last) and the numerical simulations were performed with NEST version 2.18 (https://zenodo.org/record/2605422).

The MF design, the timing optimization, the MF validation, and the MF predictive simulations were implemented in Python 3.8. Functions packages written for the present work are available on https://github.com/RobertaMLO/CRBL_MF and a DEMO is released on EBRAINS platform at https://wiki.ebrains.eu/bin/view/Collabs/mean-field-cerebellar/.

All optimisation procedures and simulations were run on a Desktop PC provided with AMD Ryzen 7 2700X CPU @ 2.16GHz with 32 GB RAM in Ubuntu 16.04.7 LTS (OS).
# 3 Results

## 3.1 The cerebellar MF

The workflow for reconstructing the cerebellar MF is shown in Fig 1 leading to a condensed representation with 4 neuronal populations for GrC, GoC, PC, MLI neurons (Fig 2). The MF was designed based on structural and functional parameters extracted from SNN simulations and the time constant was optimized with LFP experimental data. The MF working frequencies were extracted from NEST simulations of the cerebellar SNN, exploring multiple $v_{\text{drive}}$ from 4 to 80 Hz. Then these frequency ranges were used to set different plausible presynaptic signals in defining the Numerical TFs of each population (Fig 3A). The ranges were [0.42, 24.17] Hz for GrCs, [3.63, 183.15] Hz for GoCs, and [3.27, 41.66] Hz for MLIs. Note that PC working frequencies were not computed since PC activity is only projected forward to the cerebral cortex, therefore PCs never play the role of presynaptic population in this cerebellar cortical microcircuit. The $\alpha$ parameters that maximised the fitting performance for each population were $\alpha_{\text{GrC}} = 2$, $\alpha_{\text{GoC}} = 1.3$, $\alpha_{\text{MLI}} = 5$ and $\alpha_{\text{PC}} = 5$. The fitted coefficients P are reported in Table 1

2D Analytical TFs show a sigmoidal trend in relation with excitatory inputs (Fig 3B). GoC inhibition strongly affects the GrC Analytical TF; for $v_{\text{GoC}}$ higher than 100 Hz, GrC Analytical TF is almost 0 Hz. For low inhibition, e.g., $v_{\text{GoC}} = 13$ Hz, GrC Analytical TF is almost linear. MLI Analytical TF presents a well-defined sigmoidal trend depending on $v_{\text{GrC}}$ and modulated by the auto-inhibition, with resulting activity frequency spanning from 0 Hz up to 200 Hz. PC Analytical TF presents an increasing trend ranging from 0 to 100 Hz in relation to $v_{\text{GrC}}$ from 0 to 25 Hz, with the modulation due to the inhibitory control from MLIs. 3D Analytical TF of GoC shows a linear trend both for low and high $v_{\text{drive}}$.

The cerebellar MF resulted in a set of 20 second order differential equations including the specific population TFs, where 4 equations described the time variation of population activity, and the remaining 16 equations modelled the covariances of the interconnected populations. Fig 4 shows the result of T optimisation: for $T = 3.5$ ms the average granular layer activity (purple line) interpolates the experimental LFPs (red dots) with a mean absolute error of 3%.

Additionally, the overall input/output relationship of the system was analyzed by considering the input drive as a step function (see Section 2.3, protocol i) at different rates, and the output system response as the resulting PC activity. First, when input drive maintained an intensity constantly above 70 Hz, the system became unstable, and the outcome is not compatible with the physiological ranges for PC activity anymore (Fig 5A). While, within the cerebellar MF validity regime (drive in [0,70] Hz), the gain was found to be linearly dependent on the input intensity but three different sub-regimes emerged with evident discontinuities at 12 and 35 Hz (Fig 5B).

## 3.2 Constructive and functional validity

The validation of the cerebellar MF was obtained generating neuronal population dynamics with different stimulation protocols and comparing them with the corresponding SNN simulations (Fig 6). For all the protocols, the simulation lasted 500 ms and was performed with the hardware and software specified in section 2.5.

*Step function (i).* The MF was tested with a 250 ms@50 Hz step function on the mossy fibers, simulating a stimulus like a sound (Fig 6A) [75]. GrC activity rapidly raised at the beginning of the step input, then strongly decreased due to inhibitory GoC activity. GoCs, after an initial small peak due to both the direct incoming input and the GrC excitation, maintained a steady-state activity throughout the step duration. The dynamics of GrC-GoC interplay faithfully
<figure><img src="image_10.png" /><figcaption>**Fig 5.** Cerebellar mean-field overall input/output relationship. A) The average response of output PC population to different levels of the driving input ($v_{drive}$) was computed, using steps at different rates in the range [0, 100] Hz as input on mossy fiber. The shaded areas show the standard deviation of PC activity computed when the $v_{drive}$ is active. As expected from the used MF formalism, the cerebellar MF works only for low input rates and specifically the critical point above which the system is not in a physiological regime anymore is at 70 Hz [48]. B) Within the MF validity regime, cerebellar MF I/O relationship resulted always linear, but we can identify 3 different sub-regimes depending on the input range: low firing rate in [0,12] Hz, intermediate firing in [13,35] Hz, high firing in [36,70] Hz. Each sub-regime is characterized by a different gain quantified in terms of fitting parameters obtained from a simple linear fitting procedure ($f = ax+b$) performed separately on each sub-regime (dashed lines).</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1011434.g005

reproduced the feedback loop between GoCs and GrCs. GrC was the excitatory input for the molecular layer, and both MLI and PC activity arose in correspondence of the GrC initial peak. Thus, exploiting network di-synaptic delays, MLIs reduced PC activity soon after its maximum, generating the typical burst-pause pattern of these neurons. The PC pause is due to both the single neuron parameters and to inhibitory local connectivity in the microcircuit. After this rapid transient, the activity of MLIs and PCs reached a steady-state. In the MF, fast dynamics at the input step onset and at the steady-state matched SNN simulations for all neuronal populations.

**Sinusoidal input (ii).** Simulated dynamics of all cerebellar populations showed oscillations driven by the input reproducing whisker movements (Fig 6B) [38]. GrC activity projected to the molecular layer a sinusoidal-shaped signal at 0.05-5Hz, contributing to an oscillatory behaviour in GoCs, and causing an oscillation in MLI between at 23-41 Hz, and in PC activity at 42-69 Hz. Oscillations had comparable amplitude in MF and SNN simulations and occurred in the same frequency ranges (except for MLI activity that was slightly higher in MF that SNN).

**Cortical-like input (iii).** Combination of sinusoids (with EEG-like frequency [76] of 1 Hz, 15 Hz, 30 Hz, respectively) caused an irregular oscillation in the input carried by mossy fibers (2-47 Hz range) (Fig 6C) Oscillations had comparable amplitude in MF and SNN simulations and occurred in the same frequency ranges.

**Sensory-like input (iv).** The summation of step function (i) with cortical-like input (iii) resulted in an irregular input (Fig 6D) depicting in-vivo noisy baseline activity with a sound
<figure><img src="image_11.png" /><figcaption>**Fig 6. Constructive and functional validity.** Comparison of Spiking Neural Network (SNN) and mean-field (MF) activity with standard deviation in cerebellar cortical populations (GrC = Granule Cells, GoC = Golgi Cells, MLI = Molecular Layer Interneurons, PC = Purkinje Cells) in response to different driving input ($v_{drive}$) patterns, lasting 500 ms with time resolution = 0.1 ms. **A)** step function (sound-like stimulus); **B)** sinusoidal input (whisker-like); **C)** Cortical-like input **D)** combined stimulus summing (A) step function and (C) cortical-like sinusoidal. The trace of MF activity (line = average, shade = standard deviation) is overlayed to the spiking activity, which is represented as a Peristimulus Time Histogram (PSTH, time bins of 15ms). In all cases, MF activity is within physiological ranges, capturing also fast changes of activity due to instantaneous input changes in step-function input protocols. The boxplots of PC simulated activity with SNN and MF (inset in panels), shows that the MF is able to respond to the different stimulation patterns within the same frequency ranges of SNN.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1011434.g006

superimposed. GrC activity faithfully transmitted the driving input, with peaks at ~21 Hz. In correspondence with the GrC excitatory peak, MLIs peaked at ~130 Hz and PCs at ~100 Hz.

For each stimulation protocol, the PC activity predicted by the MF laid within the variability range of the corresponding SNN (the RMSE between MF and SNN was ~ 30%).
### 3.3 Predictive validity

The PC response to a 50Hz step-stimulus was described by a peak at 97 Hz followed by a pause down to 68 Hz; then a steady-state of 78 Hz was attained. MLI-PC and GrC-PC modulation (Fig 7) perturbed this reference condition.

#### 3.3.1 MLI-PC feed forward inhibition.
Inhibitory interneurons control the generation of burst-pause patterns in PC, which is fundamental for shaping the cerebellar output during motor learning [77,78]. For instance, knock-out of MLI inhibition on PC impacts on vestibulo-ocular reflex adaptation [71]. Here in MF simulations, when the MLI-PC conductance was reduced to 5% of the reference condition, the burst-pause dynamics of PC was lost, so that the PC firing settled directly back to baseline (which was elevated due to lack of inhibitory control). Conversely, when PCs were over-inhibited by the MLIs (MLI-PC conductance increased to 250% of reference condition), the pause was deeper. The PC overall activity (AUC), and PC Peak reveal an exponential trend that decays for higher MLI-PC conductances. The PC pause shows a decreasing sigmoidal trend for higher MLI-PC conductances (Fig 7A).

#### 3.3.2 PC plasticity.
Long Term Depression (LTD) and Long Term Potentiation (LTP) at GrC-PC synapses are thought to drive cerebellar adaptation and learning. The overall activity and the peak of PCs showed a linear and a sigmoidal trend, respectively, as the GrC-PC weight increased. With decreased GrC-PC strength, the peak was reduced or disappeared, and the steady-state activity reached lower levels. With increased GrC-PC strength, both the peak and

<figure><img src="image_12.png" /><figcaption>**Fig 7. Predictive validity.** Mean-field simulations using different strengths at PC connections. Simulations last 500 ms with a time resolution of 0.1 ms and vdrive is a step at 50 Hz lasting 250 ms starting at 125 ms. Quantitative score normalized to baseline (w = 100% in blue). **A) MLI-PC feed forward inhibition.** The change of total activity (measured as Area Under Curve (AUC)) and the initial peak amplitude decrease exponentially with the MLI-PC strength. The amplitude of the pause after the peak response decreases with MLI-PC strength following an inverse sigmoidal function. **B) PC Learning** with different GrC-PC plasticity conditions (Long Term Potentiation (LTP) and Long Term Depression (LTD)). AUC linearly increases with the GrC-PC strength, matching the experimental values in experimental learning protocols. The initial peak caused by the step stimulus onset increases with GrC-PC until saturation, following a sigmoidal function.</figcaption></figure>
steady-state values were increased (Fig 7B). During a typical cerebellum-driven behaviour, the
eyeblink classical conditioning (EBCC), a level of suppression of about 15% has been reported
and correlated with a stable generation of associative blink responses at the end of the learning
process [70]. SNN simulations obtained the same result by setting AMPA-mediated parallel fie-
bers-PC synapses = 35% [56]. In MF simulations, for $w_{GrC-PC}$ = 65%, which corresponded to a
reduction of GrC-PC conductance of 35%, PC activity presented a 22% reduction of the peak
and a 10% reduction of the AUC, falling into the experimental range of PC suppression [70].

**4. Discussion**

This work shows, for the first time, a MF of cerebellar cortex. According to its bottom-up
nature, the MF transfers the microscopic properties of neurons (including GrC, GoC, MLI,
and PC) and synapses of the cerebellar cortex into a condensed representation of neural activ-
ity through its two main statistical moments, mean and variance. The construction and valida-
tions strategies adopted here make the present MF an effective representation of the main
physiological properties of a canonical cerebellar module [79].

**4.1 MF design and validation**

**The cerebellar network and TF formalism.** The MF of the cerebellar cortex was based on the same general formalism previously developed for the MFs of the isocortex [21,80–82]. However, the cerebellar cortex MF benefitted from a previously validated cerebellar SNN to enabling to remap the multi-layer topology and the synaptic biophysical properties, avoiding collapsing these properties into an excitatory/inhibitory mono-layer network. Moreover, the input/output relationship of neurons was represented using non-linear E-GLIF models and the synapses with alpha-based conductance functions tuned on cerebellar parameters. This results in a cerebellar MF that has three key advantages over using a generic MF.

First, the parameters of population specific TFs were validated against biophysically detailed
models of neurons and the connectome was derived from precise scaffold model reconstruc-
tions providing a direct link to the biological microcircuit [56].

Secondly, the equations of $\mu_V$, $\sigma_V$ and $\tau_V$ included in the TF formalism were adapted to model the alpha-shaped synapses and to maintain rise-times in synaptic dynamics. For comparison, previous MFs [19,61] used exponential synapses, which provide a less realistic approximation due to their instantaneous rise time [83].

Finally, the cerebellar MF included 4 different species of neurons that were modelled using either 2D or 3D TFs to account for the multiplicity of their inputs (Fig 3). It is worth noting that the analysis of both 2D TFs of GrC, MLI and PC, and 3D TF of GoC were fitted considering only physiological input combinations computed from single-neuron computational models. The fitting procedure allows to fine-tune the single neuron firing responses (via the TF coefficients, Eq 12) to the set of neuronal and synaptic parameters of the network. The 3D dimension of the GoC TF allowed us to keep the excitatory input from GrC and mossy fibers separate ($v_{GrC}$ and $v_{drive}$, respectively), enabling us to investigate distinct excitatory input contributions to granular layer dynamics and to the whole cerebellar MF. By fixing the excitatory mossy fibers driving input, we assessed the power of the Analytical TF in simulating spiking network activity for inputs at both low and high frequency (e.g., see Fig 3 with $v_{drive}$ = 20 Hz and $v_{drive}$ = 80 Hz).

A technical issue incurred while fitting the numerical TF. The TF formalism models the difference between the phenomenological threshold ($V_{thre}^{eff}$) and population average responses ($\mu_V$) through the complementary error function (erfc in Eq (5)), providing an immediate interpretation of how a single neuron activity was related with statistics of population dynamics
($\mu_V$, $\sigma_V$, $\tau_V$). Since the erfc is a sigmoidal-shaped function, it did not accurately follow the numerical TF distribution at the boundaries, loosing precision at high frequencies. This problem was circumvented by tuning a factor, alpha (Eq 5), similarly to what was done for isocortical models [21], which does not represent a physiological quantity, but is able to extend the TF reliably over the whole range of physiological input frequencies.

**MF tuning.** The inclusion of specific structural and functional parameters in the design of cerebellar MF (see Figs 2 and 3) generated a biology grounded model that could be validated at a mesoscopic scale with the prediction of cerebellar dynamics (multi-layer Eq 13). The dynamics of the cerebellar cortex are several times faster compared to those of the cerebral cortex [84], so that the MF time constant, T, must be optimized accordingly. The MF time constant was optimized using experimental LFP recordings from the cerebellar granular layer acquired on the same spatial scale of the MF. The best fitting was obtained by accounting for the smaller contribution of GoC than GrC activity (13% vs. 87% [85, 86]) to LFPs [67], revealing that the MF time constant of the cerebellar cortex is T = 3.5 ms with mean absolute error of 3% (Fig 3). T is definitely smaller than in cerebral cortex MFs, which can assume values up to 20 ms [21, 61, 81], and captures the specific high speed dynamics of the cerebellum [84]. This result further confirms the need of a MF specifically tailored on the cerebellum functional and topological parameters.

The optimal T value was plugged into Eq 13 resulting in a second order differential equation system of interdependent TFs capturing the dynamics of multiple cerebellar populations and their covariances. This rich pool of equations allowed our cerebellar MF to reproduce a variety of cerebellar dynamics in response to different inputs (see **section 2.3** and shown in Fig 5) which, by comparison with the equivalent SNN output, provided the benchmark for constructive and functional validity. A rapidly changing input like a step function reproduced a stimulus [75] carried by the peripherical mossy fibers, causing rich dynamics in the cerebellar cortex including the typical PC burst-pause responses [87]. Adding a multi-sinusoidal input to the stimulus replicated a more physiological condition accounting for background activity. This resulted in rich PC dynamics, which still maintained burst-pause responses. A sinusoidal input was meant to emulate more complex experimental conditions, like those determined by whisker movement [38, 50, 88–91]. In particular, a multi-sinusoidal waveform was used to emulate composite inputs from the cerebral cortex [76, 92]. This makes the MF model of the cerebellum an appropriate tool to investigate different pathological conditions, that could eventually be tested in the context of a cerebro-cerebellar loop.

The cerebellar MF reproduced the known aspects of circuit physiology revealed experimentally. These include: the ability of GrCs both to respond to impulsive inputs with bursts curtailed by GoC inhibition and to faithfully follow slow input fluctuations [93]; the ability of PCs to generate burst-pause responses accentuated by MLIs [48, 51]. The cerebellar MF behaved as a smoothing filter reproducing the trend of SNN simulations without capturing the intra-population variability (e.g., GrC fluctuations in response to a constant input—Fig 6A). Since the intra-population variability was not accounted for, the maximal GoC responses were underestimated by the MF, which approximated GoG to a homogeneous population neglecting their biological heterogeneity [87]. The differences that could be observed between the MF and the corresponding SNN simulations could thus be explained as a drawback of the homogeneity assumption. Despite these aspects (further discussed in section 4.4) the MF predicted with good approximation the SNN cerebellar output, i.e., the PC activity (Boxplot in Fig 6).

## 4.2 MF predictions

A critical step in model validation is to demonstrate its ability to predict functional states not used for model construction. The cerebellum is well known for the ability to change its
network dynamics because of synaptic and non-synaptic plasticity. However, differently from SNN, the MF did not include plasticity mechanisms. Thus, we directly tested the MF ability to predict the effects of plasticity expression by mapping a set of precomputed synaptic changes on the MF itself. MF reproduced the impact of changes in MLI-PC conductance confirming that it can reproduce the complex burst-pause behaviour of PC, tuned through the MLI-PC connectivity.

The cerebellar MF reproduced the experimental recordings in EBCC experiments on behaving mice [70] pointed out a PC LTD of 10% in terms of overall activity and 22% for the peak (Fig 7 $w_{GrC-PC}$ = 65%). This protocol corresponds to a reduction of 35% of AMPA-mediated parallel fibers-PC in the corresponding SNN simulation, therefore our MF capability of capturing synaptic mechanism is further validated against in-vivo recordings.

In aggregate, the cerebellar MF was able to reproduce complex physiological mechanisms and predict the activity changes caused by synaptic modulation [56,70]. Although the cerebellar MF was able to replicate changes in activity patterns characterizing learning states, it was not possible yet to simulate the long-term plasticity rules that would be needed to autonomously generate the PC learning curves. Therefore, this MF is a flexible tool that can be used to investigate physiological or pathological properties by tuning the input of the TF on a target population (TF<sub>PC</sub> in these simulations), without complicating neither the TF fitting procedure nor the model equations. In contrast to conventional MFs, which are often a black box in terms of physiological parametrization, the cerebellar MF internal model parameters inform about the average properties and variance of fundamental mechanisms in the circuit (e.g., intrinsic and synaptic excitation) and can therefore be used to remap biological properties.

The procedure of parameters tuning might pave the way for further manipulation to remap physiological and/or pathological features onto the MF. As an example, identifying and extracting biophysical meaningful features from atlases, such as the receptor maps (e.g., extracted from Allen Brain Atlas, http://ww.brain-map.org), might allow combining mesoscopic MF simulation with characteristics of molecular functions.

## 4.3 Performance vs. realism

The MF approximated a complex SNN of ~3x10⁴ neurons and ~1x10⁶ synapses (section 2.1) with 20 equations reducing the computational time by 60% with lower memory requirements. Nevertheless, TFs fitting could be improved replacing the procedure in section 2.2.2 with a lookup table-based algorithm, which might yield to a gain in computational time up to an order of magnitude. This will represent a definite advantage when performing long-lasting simulations reflecting the acquisition time of *in vivo* recordings like EEG and fMRI or when simulating learning processes in closed-loop controllers. On the other hand, thanks to its bottom-up nature it was possible to maintain the biological realism in responses to various stimuli, including simulations of learning-induced firing rate changes and pathological conditions at the neuronal population level, obtaining a good balance between computational load and biological plausibility. This will allow to make predictions on the underlying neural bases of ensemble brain signals and to identify the elementary causes of signal alterations in pathology.

Theoretically, the cerebellar MF is not just a replacement for SNN, but rather a specific instrument for addressing questions at a mesoscale. SNNs and MFs indeed cover complementary aspects of brain function and dysfunction. For example, in Parkinson disease, a SNN was successfully used to simulate the thalamic and basal ganglia microcircuits in order to control deep brain stimulation electrodes [94,95]. The SNN provided a description of microcircuit activity and of the involvement of dopaminergic neurons in tremor [96,97]. On the other
hand, MFs of the basal ganglia and thalamocortical circuit were used to embed neural mecha-
nisms into the motor network involved in Parkinson disease and investigate large-scale oscilla-
tions [98]. In dystonia, a cerebellar SNN embedding Dystonic alterations was successfully used
to simulate alterations in PC firing rate, olivocerebellar connectivity, and parallel-fiber/PC
spike-timing dependent plasticity (STDP) [99]. While these microscopic properties may not
be reproduced by the cerebellar MF, it may however reproduce the overall PC output alter-
ations allowing to simulate the large-scale propagation of cerebellar abnormalities to other
brain regions.

## 4.4 Limits and further developments

The MF presented here simplifies the representation of neural activity and enables computa-
tionally efficient analysis of network dynamics. The tradeoff of simplification is the loss of
some properties that the real system has. In general, MFs present a coarser granularity than
SNNs, limiting studies in which single neurons, single spikes, or neural heterogeneity matters.
Furthermore, MFs and SNNs perform a different information encoding. SNNs encode infor-
mation with spike-time precision, while MFs are rate-coded machines modeling the average
activity of entire neuronal populations. For example, STDP may not be simulated with a MF,
where single spike times are not represented.

The chosen mathematical formalism can also impact on the performance of the MF. MFs based on stationary TFs are not intended to process frequency oscillations unlike MFs resulting from analytical derivations of the time-varying rate (e.g. available for two populations model) [100,101]. Specifically, translating the cerebellar SNN into the cerebellar MF involves simplifications (e.g., lack of plasticity rules) and assumptions (e.g., homogeneity inside the populations). In practice, this results in a cerebellar MF with the predictive capability of reproducing the effects of complex dynamics like burst-pause of PCs at a population level, but without modelling the neuronal coding over an heterogeneous neuronal population [102,103]. For example, a cerebellar SNN was used to simulate the control of saccades movements exerted by two subpopulation of PCs [104]. Conversely, the cerebellar MF assumes intra-population homogeneity and can't be used to study such phenomena.

Besides the limitations intrinsic to a population-level formalism, the cerebellar MF could be further developed by adding climbing fibers (cfs) as an additional input to PCs. The PC TFs including cfs input could catch further complexity and MF predictions might generate burst-pause responses. Another improvement of the cerebellar MF could include synaptic plasticity rules like in the corresponding SNN [105]. Such mechanisms at a mesoscopic scale could be modeled by implementing changes of synaptic strength based on the firing rates.

Considering the regime of applicability, the formalism of our cerebellar MF is intended to work at an asynchronous state with irregular firing rate. Therefore, input oscillations are not processed by this formalism [106,107]. Oscillations in Fig 6 are merely input-driven. To include frequency elaborations as intrinsic processes of the cerebellar MF itself, further improvements in the MF formalism should be introduced, e.g., embedding mechanisms to account for adaptation and neuronal heterogeneity [108]. The proposed improvements of the cerebellar MF presented here are all possible in principle and their implementation will depend on the specific use of the MF in future applications.

## 4.5 Conclusions and perspectives

In aggregate, the cerebellar cortex MF enforces a bottom-up approach tailored to the specific structural and functional interactions of the local neuronal populations and has a robust constructive and functional validity. By accounting for a variety of representative patterns of
discharge in cerebellar cortical neurons, the present MF can be considered a proxy of the real biological network. The cerebellar MF model developed here was validated against mice data and our promising results prompt for its translation to humans. In principle, the fundamental network architecture will remain the same, while parameters of neuron and synaptic models will require specific retuning. Further model validations should be implemented on histological and electrophysiological data from human cerebellar slices. The internal model parameters inform about the average properties and variance of fundamental mechanisms in the circuit, namely intrinsic and synaptic excitation, and can therefore be used to remap biological properties onto the MF [28]. In future applications, this will allow to tune the MF towards specific functional or dysfunctional states that affect the cerebellum. Among these it is worth mentioning ataxias [109,110], paroxysmal dyskinesia [111,112], dystonia [113,114], autistic spectrum disorders [115,116] as well as other neurological pathologies like multiple sclerosis [117–119], dementia [39] and Parkinson disease [120,121], in which a cerebellum involvement has been reported. The cerebellar MF could be applied to whole-brain simulators using TVB, similarly to what has been done already for the isocortical MF [14,122–125]. Different modules of cerebellar MFs could be connected one to another and with other cortical and subcortical MFs following subject-specific connectome rules that could be extracted from imaging techniques such as connectivity matrices computed from diffusion weighted imaging [126]. Considering the specificity of signal processing in different brain regions, this approach represents a leap ahead towards adopting specific neural mass models for different brain regions.

On the theoretical side, TVB simulations using classical neural masses [127] for all brain nodes [39,40] can now be compared to those using the cerebellum MF. At the other extreme of the spectrum, TVB with embedded cerebellar MF can be compared to TVB-NEST co-simulations [95], in which spiking neurons are represented explicitly [52,54,56,57].

In conclusion, the cerebellar MF represents the first step toward a new generation of models capable of bearing biological properties into virtual brains that will allow to simulate the healthy and pathological brain towards the overarching aim of a personalized brain representation and the technology of brain digital twins [128].

**Supporting information**

**S1 Table. Neuron parameters.** Parameters specific of the type of neurons included in the multi-layer MF populations. The parameters in the top part are chosen according to literature, while the parameters at the bottom were extracted from spiking neural network simulating the cerebellar cortex spiking activity. mf = mossy fibers, GrC = Granule Cells, GoC = Golgi Cells, MLI = Molecular Layer Interneurons (Basket cells and Stellate cells) (DOCX)

**S2 Table. Presynaptic parameters.** Parameters used to set up the inter-population connectivity of the multi-layer MF cerebellar network. Parameters were extracted from the spiking neural network simulating the cerebellar cortex spiking activity. mf = mossy fibers, GrC = Granule Cells, GoC = Golgi Cells, MLI = Molecular Layer Interneurons (Basket cells and Stellate cells). K = pre-synaptic connectivity resulting by weighting the mean synaptic convergence with the number of synapses; Q = quantal synaptic conductance in nS; τ = synaptic decay time constant; E<sub>rev</sub> = reversal potential that is 0 V for excitatory synaptic connections and -80 V for inhibitory synaptic connections. (DOCX)

**S3 Table. Mean-field symbols.**
(DOCX)
**S1 Fig. Correlation matrices between population activities, in MF and SNN networks.** Pop-
ulation activities predicted by cerebellar MF were highly cross-correlated, resulting in Pearson
Correlation Coefficients always > 0.7 (Panel A). Correlation between population activities
simulated by SNN is reported in Panel B; variability intrinsic in SNN led to lower correlations
between pairs of population activities than in MF. These correlations matrices are reported for
each input pattern. The input frequencies might be replaced with a probabilistic kernel to
introduce variability in the MF formalism and reduce the inter-population correlations.
(TIFF)

Acknowledgments

We thank Robin De Schepper for useful discussion on Brain Scaffold Builder framework and
Alessio Marta for an insight into the mathematical formalism of mean field.

**Author Contributions**

**Conceptualization:** Claudia A. M. Gandini Wheeler-Kingshott, Fulvia Palesi, Claudia Casellato, Egidio D'Angelo.

**Data curation:** Roberta Maria Lorenzi.

**Formal analysis:** Roberta Maria Lorenzi, Marialaura De Grazia.

**Funding acquisition:** Yann Zerlaut, Alain Destexhe, Claudia A. M. Gandini
Wheeler-Kingshott, Fulvia Palesi, Egidio D'Angelo.

**Methodology:** Roberta Maria Lorenzi, Alice Geminiani.

**Project administration:** Claudia Casellato, Egidio D'Angelo.

**Supervision:** Yann Zerlaut, Alain Destexhe, Fulvia Palesi, Claudia Casellato, Egidio D'Angelo.

**Validation:** Roberta Maria Lorenzi.

**Visualization:** Roberta Maria Lorenzi.

**Writing – original draft:** Roberta Maria Lorenzi, Alice Geminiani, Egidio D'Angelo.

**Writing – review & editing:** Roberta Maria Lorenzi, Alice Geminiani, Yann Zerlaut, Marialaura De Grazia, Alain Destexhe, Claudia A. M. Gandini Wheeler-Kingshott, Fulvia Palesi, Claudia Casellato, Egidio D'Angelo.

**References**

1. Schirner M, Rothmeier S, Jirsa VK, McIntosh AR, Ritter P. An automated pipeline for constructing personalized virtual brains from multimodal neuroimaging data. Neuroimage [Internet]. 2015; 117:343-57. Available from: https://doi.org/10.1016/j.neuroimage.2015.03.055 PMID: 25837600

2. Amunts K, Lepage C, Borgeat L, Mohlberg H, Dickscheid T, Rousseau MÈ, et al. BigBrain: An ultra-high-resolution 3D human brain model. Science (80-). 2013; 340(6139):1472-5. https://doi.org/10.1126/science.1235381 PMID: 23788795

3. D'Angelo E, Jirsa V. Neurosciences The quest for multiscale brain modeling. Trends Neurosci [Internet]. 2022; xx(xx):1-14. Available from: https://doi.org/10.1016/j.tins.2022.06.007

4. Plesser HE, Diesmann M, Gewaltig M-O, Morrison A. NEST: the Neural Simulation Tool BT—Encyclopedia of Computational Neuroscience. In: Jaeger D, Jung R, editors. New York, NY: Springer New York; 2015. p. 1849–52. Available from: https://doi.org/10.1007/978-1-4614-6675-8_258

5. Yavuz E, Turner J, Nowotny T. GeNN: a code generation framework for accelerated brain simulations. Sci Rep. 2016; 6(1):1–14.

6. Knight JC, Komissarov A, Nowotny T. PyGeNN: a Python library for GPU-enhanced neural networks. Front Neuroinform. 2021; 15:659005. https://doi.org/10.3389/fninf.2021.659005 PMID: 33967731
7. Wilson HR, Cowan JD. Excitatory and Inhibitory Interactions in Localized Populations of Model Neurons. Biophys J [Internet]. 1972; 12(1):1-24. Available from: https://www.sciencedirect.com/science/article/pii/S0006349572860685 https://doi.org/10.1016/S0006-3495(72)86068-5 PMID: 4332108

8. Jansen BH, Rit VG. Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns. Biol Cybern [Internet]. 1995; 73(4):357-66. Available from: https://doi.org/10.1007/BF00199471 PMID: 7578475

9. Deco G, Jirsa VK, Robinson PA, Breakspear M, Friston K. The dynamic brain: From spiking neurons to neural masses and cortical fields. PLoS Comput Biol. 2008; 4(8). https://doi.org/10.1371/journal.pcbi.1000092 PMID: 18769680

10. Montbrió E, Pazó D, Roxin A. Macroscopic description for networks of spiking neurons. Phys Rev X. 2015; 5(2):1-15.

11. Coombes S, Byrne Á. Next Generation Neural Mass Models BT—Nonlinear Dynamics in Computational Neuroscience. In: Corinto F, Torcini A, editors. Cham: Springer International Publishing; 2019. p. 1-16. Available from: https://doi.org/10.1007/978-3-319-71048-8_1

12. Schwalger T, Deger M, Gerstner W. Towards a theory of cortical columns: From spiking neurons to interacting neural populations of finite size. Vol. 13, PLoS Computational Biology. 2017. 1-63 p. https://doi.org/10.1371/journal.pcbi.1005507 PMID: 28422957

13. Fellner A, Heshmat A, Werginz P, Rattay F. A finite element method framework to model extracellular neural stimulation. J Neural Eng. 2022; 19(2). https://doi.org/10.1088/1741-2552/ac6060 PMID: 35320783

14. El Boustani S, Destexhe A. A master equation formalism for macroscopic modeling of asynchronous irregular activity states. Neural Comput. 2009; 21(1):46-100. https://doi.org/10.1162/neco.2009.02-08-710 PMID: 19210171

15. Amit DJ, Brunel N. Model of global spontaneous activity and local structured activity during delay periods in the cerebral cortex. Cereb Cortex. 1997; 7(3):237-52. https://doi.org/10.1093/cercor/7.3.237 PMID: 9143444

16. Brunel N, Sergi S. Firing frequency of integrate-and-fire neurons with finite synaptic time constants. J Theor Biol. 1998; 195(1):87-95.

17. Kumar A, Schrader S, Aertsen A, Rotter S. The High-Conductance State of Cortical Networks. Neural Comput. 2008; 20(1):1-43. https://doi.org/10.1162/neco.2008.20.1.1 PMID: 18044999

18. Kuhn A, Aertsen A, Rotter S. Neuronal Integration of Synaptic Input in the Fluctuation-Driven Regime. J Neurosci. 2004; 24(10):2345-56. https://doi.org/10.1523/JNEUROSCI.3349-03.2004 PMID: 15014109

19. Zerlaut Y, Teleńczuk B, Deleuze C, Bal T, Ouanounou G, Destexhe A. Heterogeneous firing rate response of mouse layer V pyramidal neurons in the fluctuation-driven regime. J Physiol. 2016; 594 (13):3791-808. https://doi.org/10.1113/JP272317 PMID: 27146816

20. Cakan C, Obermayer K. Biophysically grounded mean-field models of neural populations under electrical stimulation. PLoS Comput Biol [Internet]. 2020; 16(4):1-30. Available from: https://doi.org/10.1371/journal.pcbi.1007822 PMID: 32324734

21. Carlu M, Chehab O, Dalla Porta XL, Depannemaecker D, Héricé XC, Jedynak M, et al. A mean-field approach to the dynamics of networks of complex neurons, from nonlinear Integrate-and-Fire to Hodgkin-Huxley models. J Neurophysiol. 2020; 123(3):1042-51. https://doi.org/10.1152/jn.00399.2019 PMID: 31851573

22. Marreiros AC, Kiebel SJ, Friston KJ. Dynamic causal modelling for fMRI: A two-state model. Neuroimage. 2008; 39(1):269-78. https://doi.org/10.1016/j.neuroimage.2007.08.019 PMID: 17936017

23. Bastos AM, Usrey WM, Adams RA, Mangun GR, Fries P, Friston KJ. Canonical Microcircuits for Predictive Coding. Neuron [Internet]. 2012; 76(4):695-711. Available from: https://doi.org/10.1016/j.neuron.2012.10.038 PMID: 23177956

24. Deco G, Ponce-Alvarez A, Hagmann P, Romani GL, Mantini D, Corbetta M. How local excitation-inhibition ratio impacts the whole brain dynamics. J Neurosci. 2014; 34(23):7886-98. https://doi.org/10.1523/JNEUROSCI.5068-13.2014 PMID: 24899711

25. Auksztulewicz R, Friston K. Attentional enhancement of auditory mismatch responses: A DCM/MEG study. Cereb Cortex. 2015; 25(11):4273-83. https://doi.org/10.1093/cercor/bhu323 PMID: 25596591

26. Glomb K, Ponce-Alvarez A, Gilson M, Ritter P, Deco G. Resting state networks in empirical and simulated dynamic functional connectivity. Neuroimage. 2017 Oct; 159:388-402. https://doi.org/10.1016/j.neuroimage.2017.07.065 PMID: 28782678

27. El Houssaini K, Bernard C, Jirsa VK. The Epileptor Model: A Systematic Mathematical Analysis Linked to the Dynamics of Seizures, Refractory Status Epilepticus, and Depolarization Block. eNeuro. 2020; 7(2).
28. Naskar A, Vattikonda A, Deco G, Roy D, Banerjee A. Multiscale dynamic mean field (MDMF) model relates resting-state brain dynamics with local cortical excitatory-inhibitory neurotransmitter homeostasis. Netw Neurosci (Cambridge, Mass). 2021; 5(3):757-82. https://doi.org/10.1162/netn_a_00197 PMID: 34746626

29. Levenstein D, Buzsáki G, Rinzel J. NREM sleep in the rodent neocortex and hippocampus reflects excitable dynamics. Nat Commun [Internet]. 2019; 10(1):1-12. Available from: http://dx.doi.org/10.1038/s41467-019-10327-5

30. van Wijk BCM, Cagnan H, Litvak V, Kühn AA, Friston KJ. Generic dynamic causal modelling: An illustrative application to Parkinson's disease. Neuroimage [Internet]. 2018; 181(January):818-30. Available from: https://doi.org/10.1016/j.neuroimage.2018.08.039 PMID: 30130648

31. Moran RJ, Mallet N, Litvak V, Dolan RJ, Magill PJ, Friston KJ, et al. Alterations in brain connectivity underlying beta oscillations in parkinsonism. PLoS Comput Biol. 2011; 7(8). https://doi.org/10.1371/journal.pcbi.1002124 PMID: 21852943

32. Saggar M, Zanesco AP, King BG, Bridwell DA, MacLean KA, Aichele SR, et al. Mean-field thalamocortical modeling of longitudinal EEG acquired during intensive meditation training. Neuroimage. 2015 Jul; 114:88-104. https://doi.org/10.1016/j.neuroimage.2015.03.073 PMID: 25862265

33. Schutter DJLG, van Honk J. An electrophysiological link between the cerebellum, cognition and emotion: Frontal theta EEG activity to single-pulse cerebellar TMS. Neuroimage [Internet]. 2006; 33 (4):1227-31. Available from: https://www.sciencedirect.com/science/article/pii/S1053811906006768 https://doi.org/10.1016/j.neuroimage.2006.06.055 PMID: 17023183

34. Castellazzi G, Palesi F, Casali S, Vitali P, Wheeler-Kingshott CAM, Sinforiani E, et al. A comprehensive assessment of resting state networks: Bidirectional modification of functional integrity in cerebro-cerebellar networks in dementia. Front Neurosci. 2014; 8(8 JUL):1-18. https://doi.org/10.3389/fnins.2014.00223 PMID: 25126054

35. Castellazzi G, Bruno SD, Toosy AT, Casiraghi L, Palesi F, Savini G, et al. Prominent Changes in Cerebro-Cerebellar Functional Connectivity During Continuous Cognitive Processing. Vol. 12, Frontiers in Cellular Neuroscience. 2018. p. 331. https://doi.org/10.3389/fncel.2018.00331 PMID: 30327590

36. Casiraghi L, Alahmadi AAS, Monteverdi A, Palesi F, Castellazzi G, Savini G, et al. I See Your Effort: Force-Related BOLD Effects in an Extended Action Execution-Observation Network Involving the Cerebellum. Cereb Cortex. 2019; 29(3):1351-68. https://doi.org/10.1093/cercor/bhy322 PMID: 30615116

37. Andersen LM, Jerbi K, Dalal SS. Can EEG and MEG detect signals from the human cerebellum? Neuroimage. 2020 Jul; 215:116817. https://doi.org/10.1016/j.neuroimage.2020.116817 PMID: 32278092

38. Popa D, Spolidoro M, Proville RD, Guyon N, Belliveau L, Léna C. Functional role of the cerebellum in gamma-band synchronization of the sensory and motor cortices. J Neurosci. 2013; 33(15):6552-6. https://doi.org/10.1523/JNEUROSCI.5521-12.2013 PMID: 23575852

39. Monteverdi A, Palesi F, Costa A, Vitali P, Pichiecchio A, Ramusino MC, et al. Subject-specific features of excitation / inhibition profiles in neurodegenerative diseases. 2022;(August):1-17. https://doi.org/10.3389/fnagi.2022.868342 PMID: 35992607

40. Palesi F, Lorenzi RM, Casellato C, Ritter P, Jirsa V, Gandini Wheeler-Kingshott CAM, et al. The Importance of Cerebellar Connectivity on Simulated Brain Dynamics. Front Cell Neurosci. 2020; 14(July):1-11. https://doi.org/10.3389/fncel.2020.00240 PMID: 32848628

41. Sanzleon P, Knock SA, Woodman MM, Domide L, Mersmann J, McIntosh AR, et al. The virtual brain: A simulator of primate brain network dynamics. Front Neuroinform. 2013; 7(MAY).

42. Lachamp PM, Liu Y, Liu SJ. Glutamatergic Modulation of Cerebellar Interneuron Activity Is Mediated by an Enhancement of GABA Release and Requires Protein Kinase A/RIM1α Signaling. J Neurosci [Internet]. 2009 Jan 14; 29(2):381 LP-392. Available from: http://www.jneurosci.org/content/29/2/381 abstract

43. D'Angelo E, Nieus T, Maffei A, Armano S, Rossi P, Taglietti V, et al. Theta-frequency bursting and resonance in cerebellar granule cells: experimental evidence and modeling of a slow k+-dependent mechanism. J Neurosci Off J Soc Neurosci. 2001 Feb; 21(3):759-70. https://doi.org/10.1523/JNEUROSCI.21-03-00759.2001 PMID: 11157062

44. Solinas S, Forti L, Cesana E, Mapelli J, De Schutter E, D'Angelo E. Computational reconstruction of pacemaking and intrinsic electroresponsiveness in cerebellar Golgi cells. Front Cell Neurosci. 2007; 1:2. https://doi.org/10.3389/neuro.03.002.2007 PMID: 18946520

45. McKay BE, Turner RW. Physiological and morphological development of the rat cerebellar Purkinje cell. J Physiol. 2005 Sep; 567(Pt 3):829-50. https://doi.org/10.1113/jphysiol.2005.089383 PMID: 16002452

46. Molineux ML, McRory JE, McKay BE, Hamid J, Mehaffey WH, Rehak R, et al. Specific T-type calcium channel isoforms are associated with distinct burst phenotypes in deep cerebellar nuclear neurons.
Proc Natl Acad Sci U S A. 2006; 103(14):5555-60. https://doi.org/10.1073/pnas.0601261103 PMID: 16567615

47. Masoli S, Solinas S, D'Angelo E. Action potential processing in a detailed Purkinje cell model reveals a critical role for axonal compartmentalization. Front Cell Neurosci. 2015; 9:47. https://doi.org/10.3389/fncel.2015.00047 PMID: 25759640

48. Masoli S, D'Angelo E. Synaptic activation of a detailed Purkinje cell model predicts voltage-dependent control of burst-pause responses in active dendrites. Front Cell Neurosci. 2017; 11(September):1-18. https://doi.org/10.3389/fncel.2017.00278 PMID: 28955206

49. Masoli S, Tognolina M, Laforenza U, Moccia F, D'Angelo E. Parameter tuning differentiates granule cell subtypes enriching transmission properties at the cerebellum input stage. Commun Biol [Internet]. 2020; 3(1):222. Available from: https://doi.org/10.1038/s42003-020-0953-x PMID: 32385389

50. Masoli S, Ottaviani A, Casali S, D'Angelo E. Cerebellar Golgi cell models predict dendritic processing and mechanisms of synaptic plasticity. PLoS Comput Biol. 2020 Dec; 16(12):e1007937. https://doi.org/10.1371/journal.pcbi.1007937 PMID: 33378395

51. Rizza MF, Locatelli F, Masoli S, Sánchez-Ponce D, Muñoz A, Prestori F, et al. Stellate cell computational modeling predicts signal filtering in the molecular layer circuit of cerebellum. Sci Rep [Internet]. 2021; 11(1):1-18. Available from: https://doi.org/10.1038/s41598-021-83209-w

52. Geminiani A, Casellato C, D'Angelo E, Pedrocchi A. Complex Electroresponsive Dynamics in Olivocerebellar Neurons Represented With Extended-Generalized Leaky Integrate and Fire Models. Front Comput Neurosci [Internet]. 2019 Jun 6; 13(June):1-12. Available from: https://www.frontiersin.org/article/10.3389/fncom.2019.00035/full

53. Lennon W, Hecht-Nielsen R, Yamazaki T. A spiking network model of cerebellar Purkinje cells and molecular layer interneurons exhibiting irregular firing. Front Comput Neurosci. 2014; 8:157. https://doi.org/10.3389/fncom.2014.00157 PMID: 25520646

54. Geminiani A, Casellato C, Locatelli F, Prestori F, Pedrocchi A, D'Angelo E. Complex dynamics in simplified neuronal models: Reproducing golgi cell electroresponsiveness. Front Neuroinform. 2018; 12 (December):1-19. https://doi.org/10.3389/fninf.2018.00088 PMID: 30559658

55. Casali S, Tognolina M, Gandolfi D, Mapelli J, D'Angelo E. Cellular-resolution mapping uncovers spatial adaptive filtering at the rat cerebellum input stage. Commun Biol [Internet]. 2020; 3(1):635. Available from: https://doi.org/10.1038/s42003-020-01360-y PMID: 33128000

56. De Schepper R, Geminiani A, Masoli S, Rizza MF, Antonietti A, Casellato C, et al. Model simulations unveil the structure-function-dynamics relationship of the cerebellar cortical microcircuit. Commun Biol. 2022; https://doi.org/10.1038/s42003-022-04213-y PMID: 36376444

57. Geminiani A, Pedrocchi A, D'Angelo E, Casellato C. Response Dynamics in an Olivocerebellar Spiking Neural Network With Non-linear Neuron Properties. Front Comput Neurosci. 2019; 13(October):1-15.

58. Friston KJ, Harrison L, Penny W. Dynamic causal modelling. Neuroimage [Internet]. 2003 Aug; 19 (4):1273-302. Available from: https://linkinghub.elsevier.com/retrieve/pii/S1053811903002027 https://doi.org/10.1016/S1053-8119(03)00202-7 PMID: 12948688

59. Friston KJ, Preller KH, Mathys C, Cagnan H, Heinzle J, Razi A, et al. Dynamic causal modelling revisited. Neuroimage [Internet]. 2019; 199(February 2017):730-44. Available from: https://doi.org/10.1016/j.neuroimage.2017.02.045 PMID: 28219774

60. Parr T, Sajid N, Friston KJ. Modules or mean-fields? Entropy. 2020; 22(5):1-25.

61. Zerlaut Y, Chemla S, Chavane F, Destexhe A. Modeling mesoscopic cortical dynamics using a mean-field model of conductance-based networks of adaptive exponential integrate-and-fire neurons. J Comput Neurosci. 2018; 44(1):45-61. https://doi.org/10.1007/s10827-017-0668-2 PMID: 29139050

62. Roth A, van Rossum MCW. Modeling Synapses. Comput Model Methods Neurosci. 2013;(November 2014):139-60.

63. Jordan J, Mørk H, Vennemo SB, Terhorst D, Peyser A, Ippen T, et al. NEST 2.18.0. 2019 Jun 27 [cited 2022 Oct 17]; Available from: https://zenodo.org/record/2605422

64. Melik-Musyan AB, Fanardzhyan V V. Morphological Characteristics of Lugaro Cells in the Cerebellar Cortex. Neurosci Behav Physiol [Internet]. 2004; 34(6):633-8. Available from: https://doi.org/10.1023/b:neab.0000028297.30474.f9 PMID: 15368913

65. Mugnaini E, Sekerková G, Martina M. The unipolar brush cell: A remarkable neuron finally receiving deserved attention. Brain Res Rev [Internet]. 2011; 66(1):220-45. Available from: https://www.sciencedirect.com/science/article/pii/S0165017310001128 https://doi.org/10.1016/j.brainresrev.2010.10.001 PMID: 20937306

66. Powell K, Mathy A, Duguid I, Häusser M. Synaptic representation of locomotion in single cerebellar granule cells. Raman IM, editor. Elife [Internet]. 2015; 4:e07290. Available from: https://doi.org/10.7554/eLife.07290 PMID: 26083712
67. Mapelli J, D'Angelo E. The spatial organization of long-term synaptic plasticity at the input stage of cerebellum. J Neurosci. 2007; 27(6):1285-96. https://doi.org/10.1523/JNEUROSCI.4873-06.2007 PMID: 17287503

68. Monteverdi A, Di Domenico D, D'Angelo E, Mapelli L. Anisotropy and Frequency Dependence of Signal Propagation in the Cerebellar Circuit Revealed by High-Density Multielectrode Array Recordings. Vol. 11, Biomedicines. 2023. https://doi.org/10.3390/biomedicines11051475 PMID: 37239146

69. Romano V, de Propris L, Bosman LWJ, Warnaar P, Ten Brinke MM, Lindeman S, et al. Potentiation of cerebellar purkinje cells facilitates whisker reflex adaptation through increased simple spike activity. Elife. 2018; 7:1-33. https://doi.org/10.7554/eLife.38852 PMID: 30561331

70. ten Brinke MM, Boele HJ, Spanke JK, Potters JW, Kornysheva K, Wulff P, et al. Evolving Models of Pavlovian Conditioning: Cerebellar Cortical Dynamics in Awake Behaving Mice. Cell Rep [Internet]. 2015; 13(9):1977-88. Available from: https://doi.org/10.1016/j.celrep.2015.10.057 PMID: 26655909

71. Wulff P, Schonewille M, Renzi M, Viltono L, Sassoè-Pognetto M, Badura A, et al. Synaptic inhibition of Purkinje cells mediates consolidation of vestibulo-cerebellar motor learning. Nat Neurosci. 2009 Aug; 12(8):1042-9. https://doi.org/10.1038/nn.2348 PMID: 19578381

72. Bliss TVP, Cooke SF. Long-term potentiation and long-term depression: A clinical perspective. Clinics. 2011; 66(SUPPL.1):3-17. https://doi.org/10.1590/s1807-59322011001300002 PMID: 21779718

73. Mittmann W, Häusser M. Linking synaptic plasticity and spike output at excitatory and inhibitory synapses onto cerebellar Purkinje cells. J Neurosci. 2007; 27(21):5559-70. https://doi.org/10.1523/JNEUROSCI.5117-06.2007 PMID: 17522301

74. Prestori F, Bonardi C, Mapelli L, Lombardo P, Goselink R, De Stefano ME, et al. Gating of Long-Term Potentiation by Nicotinic Acetylcholine Receptors at the Cerebellum Input Stage. PLoS One. 2013; 8 (5). https://doi.org/10.1371/journal.pone.0064828 PMID: 23741401

75. Jirenhed DA, Hesslow G. Time course of classically conditioned purkinje cell response is determined by initial part of conditioned stimulus. J Neurosci. 2011; 31(25):9070-4. https://doi.org/10.1523/JNEUROSCI.1653-11.2011 PMID: 21697357

76. Del Percio C, Drinkenburg W, Lopez S, Infarinato F, Bastlund JF, Laursen B, et al. On-going electroencephalographic rhythms related to cortical arousal in wild-type mice: the effect of aging. Neurobiol Aging [Internet]. 2017; 49:20-30. Available from: https://www.sciencedirect.com/science/article/pii/S0197458016302172 https://doi.org/10.1016/j.neurobiolaging.2016.09.004 PMID: 27728831

77. Kim J, Augustine GJ. Molecular Layer Interneurons: Key Elements of Cerebellar Network Computation and Behavior. Neuroscience [Internet]. 2021; 462:22-35. Available from: https://www.sciencedirect.com/science/article/pii/S0306452220306618 https://doi.org/10.1016/j.neuroscience.2020.10.008 PMID: 33075461

78. Casali S, Marenzi E, Medini C, Casellato C, D'Angelo E. Reconstruction and simulation of a scaffold model of the cerebellar network. Front Neuroinform. 2019; 13(May):1-19.

79. D'Angelo E, Casali S. Seeking a unified framework for cerebellar function and dysfunction: from circuit operations to cognition [Internet]. Vol. 6, Frontiers in Neural Circuits. 2013. Available from: https://www.frontiersin.org/articles/10.3389/fncir.2012.00116

80. Huang CH, Lin CCK. A novel density-based neural mass model for simulating neuronal network dynamics with conductance-based synapses and membrane current adaptation. Neural Networks [Internet]. 2021; 143:183-97. Available from: https://doi.org/10.1016/j.neunet.2021.06.009 PMID: 34157643

81. Di Volo M, Romagnoni A, Capone C, Destexhe A. Biologically Realistic Mena-Field Models of Conductance-Based Networks of Spiking Neurons with Adaptation. Neural Comput. 2019; 31:653-80.

82. Moran R, Pinotsis DA, Friston K. Neural masses and fields in dynamic causal modelling. Front Comput Neurosci. 2013; 7(APR 2013):1-12.

83. Brette R, Gerstner W. Adaptive exponential integrate-and-fire model as an effective description of neuronal activity. J Neurophysiol. 2005; 94(5):3637-42. https://doi.org/10.1152/jn.00686.2005 PMID: 16014787

84. D'Angelo E. Neural circuits of the cerebellum: hypothesis for function. J Integr Neurosci. 2011 Sep; 10 (3):317-52. https://doi.org/10.1142/S0219635211002762 PMID: 21960306

85. Dieudonné S. Submillisecond kinetics and low efficacy of parallel fibre-Golgi cell synaptic currents in the rat cerebellum. J Physiol. 1998; 510(3):845-66. https://doi.org/10.1111/j.1469-7793.1998.845bj.x PMID: 9660898

86. D'Angelo E, Rossi P, Armano S, Taglietti V. Evidence for NMDA and mGLU receptor-dependent long-term potentiation of mossy fiber-granule cell transmission in rat cerebellum. J Neurophysiol. 1999; 81 (1):277-87. https://doi.org/10.1152/jn.1999.81.1.277 PMID: 9914288
87. Herzfeld DJ, Kojima Y, Soetedjo R, Shadmehr R. Encoding of action by the Purkinje cells of the cerebellum. Nature. 2015 Oct; 526(7573):439–42. https://doi.org/10.1038/nature15693 PMID: 26469054

88. Yamazaki T, Igarashi J. Realtime cerebellum: a large-scale spiking network model of the cerebellum that runs in realtime using a graphics processing unit. Neural networks Off J Int Neural Netw Soc. 2013 Nov; 47:103–11. https://doi.org/10.1016/j.neunet.2013.01.019 PMID: 23434303

89. Antonietti A, Casellato C, D'Angelo E, Pedrocchi A. Model-driven analysis of eyeblink classical conditioning reveals the underlying structure of cerebellar plasticity and neuronal activity. IEEE Trans Neural Networks Learn Syst. 2017; 28(11):2748–62. https://doi.org/10.1109/TNNLS.2016.2598190 PMID: 27608482

90. Gagliano G, Monteverdi A, Casali S, Laforenza U, Wheeler-Kingshott CAMG, D'Angelo E, et al. Non-linear frequency-dependence of neurovascular coupling in the cerebellar cortex implies vasodilation-vasoconstriction competition. bioRxiv [Internet]. 2021; 2021.07.30.454430. Available from: https://www.biorxiv.org/content/10.1101/2021.07.30.454430v1%0Ahttps://www.biorxiv.org/content/10.1101/2021.07.30.454430v1 abstract

91. Gagliano G, Monteverdi A, Casali S, Laforenza U, Gandini Wheeler-Kingshott CAM, D'Angelo E, et al. Non-Linear Frequency Dependence of Neurovascular Coupling in the Cerebellar Cortex Implies Vasodilation&ndash;Vasoconstriction Competition. Vol. 11, Cells. 2022.

92. Tzvi E, Gajiyeva L, Bindel L, Hartwigsen G, Classen J. Coherent theta oscillations in the cerebellum and supplementary motor area mediate visuomotor adaptation. Neuroimage [Internet]. 2022; 251 (February):118985. Available from: https://doi.org/10.1016/j.neuroimage.2022.118985 PMID: 35149228

93. D'Angelo E, De Zeeuw CI. Timing and plasticity in the cerebellum: focus on the granular layer. Trends Neurosci [Internet]. 2009; 32(1):30–40. Available from: https://www.sciencedirect.com/science/article/pii/S0166223608002452 https://doi.org/10.1016/j.tins.2008.09.007 PMID: 18977038

94. Sarma S V, Cheng M, Hu R, Williams Z, Brown EN, Eskandar E. Modeling Neural Spiking Activity in the Sub-Thalamic Nucleus of Parkinson's Patients and Healthy Primates. IFAC Proc Vol [Internet]. 2008; 41(2):9081–6. Available from: https://www.sciencedirect.com/science/article/pii/S147466701640412X

95. Meier JM, Perdikis D, Blickensdörfer A, Stefanovski L, Liu Q, Maith O, et al. Virtual deep brain stimulation: Multiscale co-simulation of a spiking basal ganglia model and a whole-brain mean-field model with The Virtual Brain. Exp Neurol. 2022; 354(April). https://doi.org/10.1016/j.expneurol.2022.114111 PMID: 35569510

96. Halliday GM. Thalamic changes in Parkinson's disease. Parkinsonism Relat Disord [Internet]. 2009; 15:S152–5. Available from: https://www.sciencedirect.com/science/article/pii/S1353802009708041 https://doi.org/10.1016/S1353-8020(09)70804-1 PMID: 20082979

97. McGregor MM, Nelson AB. Circuit Mechanisms of Parkinson's Disease. Neuron [Internet]. 2019; 101 (6):1042–56. Available from: https://www.sciencedirect.com/science/article/pii/S0896627319302119 https://doi.org/10.1016/j.neuron.2019.03.004 PMID: 30897356

98. van Albada SJ, Robinson PA. Mean-field modeling of the basal ganglia-thalamocortical system. I. Firing rates in healthy and parkinsonian states. J Theor Biol. 2009; 257(4):642–63.

99. Geminiani A, Mockevičius A, D'Angelo E, Casellato C. Cerebellum Involvement in Dystonia During Associative Motor Learning: Insights From a Data-Driven Spiking Network Model. Front Syst Neurosci. 2022; 16(June):1–15. https://doi.org/10.3389/fnsys.2022.919761 PMID: 35782305

100. Bandyopadhyay A, Bernard C, Jirsa VK, Petkoski S. Mean-field approximation of network of biophysical neurons driven by conductance-based ion exchange. bioRxiv [Internet]. 2022; 2021.10.29.466427. Available from: https://www.biorxiv.org/content/10.1101/2021.10.29.466427v3%0Ahttps://www.biorxiv.org/content/10.1101/2021.10.29.466427v3 abstract

101. Bick C, Goodfellow M, Laing CR, Martens EA. Understanding the dynamics of biological and neural oscillator networks through exact mean-field reductions: a review. J Math Neurosci [Internet]. 2020; 10 (1). Available from: https://doi.org/10.1186/s13408-020-00086-9 PMID: 32462281

102. Jelitai M, Puggioni P, Ishikawa T, Rinaldi A, Duguid I. Dendritic excitation-inhibition balance shapes cerebellar output during motor behaviour. Nat Commun [Internet]. 2016; 7:1–13. Available from: https://doi.org/10.1038/ncomms13722 PMID: 27976716

103. Sedaghat-Nejad E, Pi JS, Hage P, Fakharian MA, Shadmehr R. Synchronous spiking of cerebellar Purkinje cells during control of movements. Proc Natl Acad Sci [Internet]. 2022 Apr 5; 119(14): e2118954119. Available from: https://doi.org/10.1073/pnas.2118954119 PMID: 35349338

104. Fruzzetti L, KalindিHT, Antonietti A, Alessandro C, Geminiani A, Casellato C, et al. Dual STDP processes at Purkinje cells contribute to distinct improvements in accuracy and speed of saccadic eye movements. PLoS Comput Biol. 2022 Oct; 18(10):e1010564. https://doi.org/10.1371/journal.pcbi.1010564 PMID: 36194625
105. Casellato C, Antonietti A, Garrido JA, Ferrigno G, D'Angelo E, Pedrocchi A. Distributed cerebellar plasticity implements generalized multiple-scale memory components in real-robot sensorimotor tasks [Internet]. Vol. 9, Frontiers in Computational Neuroscience. 2015. Available from: https://www.frontiersin.org/articles/10.3389/fncom.2015.00024

106. Ros H, Sachdev RNS, Yu Y, Sestan N, McCormick DA. Neocortical networks entrain neuronal circuits in cerebellar cortex. J Neurosci Off J Soc Neurosci. 2009 Aug; 29(33):10309–20. https://doi.org/10.1523/JNEUROSCI.2327-09.2009 PMID: 19692605

107. Barmack NH, Yakhnitsa V. Vestibulocerebellar Connections BT—Handbook of the Cerebellum and Cerebellar Disorders. In: Manto M, Schmahmann JD, Rossi F, Gruol DL, Koibuchi N, editors. Dordrecht: Springer Netherlands; 2013. p. 357–75. Available from: https://doi.org/10.1007/978-94-007-1333-8_18

108. Di Volo M, Destexhe A. Optimal responsiveness and information flow in networks of heterogeneous neurons. Sci Rep [Internet]. 2021; 11(1):1–11. Available from: https://doi.org/10.1038/s41598-021-96745-2

109. Rosenthal LS. Neurodegenerative Cerebellar Ataxia. Continuum (Minneap Minn). 2022 Oct; 28 (5):1409–34. https://doi.org/10.1212/CON.0000000000001180 PMID: 36222772

110. Pedroso JL, Vale TC, Braga-Neto P, Dutra LA, França MCJ, Teive HAG, et al. Acute cerebellar ataxia: differential diagnosis and clinical approach. Arq Neuropsiquiatr. 2019 Mar; 77(3):184–93. https://doi.org/10.1590/0004-282X20190020 PMID: 30970132

111. Ekmen A, Meneret A, Valabregue R, Beranger B, Worbe Y, Lamy J-C, et al. Cerebellum Dysfunction in Patients With PRRT2-Related Paroxysmal Dyskinesia. Neurology. 2022 Mar; 98(10):e1077–89. https://doi.org/10.1212/WNL.00000000000200060 PMID: 35058336

112. Mendonça MD, Alves da Silva J. A Cerebellar Wave for Paroxysmal Dyskinesia. Mov Disord. 2021 Dec; 36(12):2767. https://doi.org/10.1002/mds.28848 PMID: 34755914

113. Morigaki R, Miyamoto R, Matsuda T, Miyake K, Yamamoto N, Takagi Y. Dystonia and cerebellum: From bench to bedside. Life. 2021; 11(8). https://doi.org/10.3390/life11080776 PMID: 34440520

114. Mahajan A, Gupta P, Jacobs J, Marsili L, Sturchio A, Jinnah HA, et al. Impaired Saccade Adaptation in Tremor-Dominant Cervical Dystonia-Evidence for Maladaptive Cerebellum. Cerebellum. 2021 Oct; 20 (5):678–86. https://doi.org/10.1007/s12311-020-01104-y PMID: 31965455

115. Bruchhage MMK, Bucci M-P, Becker EBE. Cerebellar involvement in autism and ADHD. Handb Clin Neurol. 2018; 155:61–72. https://doi.org/10.1016/B978-0-444-64189-2.00004-4 PMID: 29891077

116. Kelly E, Meng F, Fujita H, Morgado F, Kazemi Y, Rice LC, et al. Regulation of autism-relevant behaviors by cerebellar-prefrontal cortical circuits. Nat Neurosci. 2020 Sep; 23(9):1102–10. https://doi.org/10.1038/s41593-020-0665-z PMID: 32661395

117. Tornes L, Conway B, Sheremata W. Multiple sclerosis and the cerebellum. Neurol Clin. 2014 Nov; 32 (4):957–77. https://doi.org/10.1016/j.ncl.2014.08.001 PMID: 25439291

118. Schreck L, Ryan S, Monaghan P. Cerebellum and cognition in multiple sclerosis. J Neurophysiol. 2018 Dec; 120(6):2707–9.

119. Savini G, Pardini M, Castellazzi G, Lascialfari A, Chard D, D'Angelo E, et al. Default mode network structural integrity and cerebellar connectivity predict information processing speed deficit in multiple sclerosis. Front Cell Neurosci. 2019; 13(February):1–15. https://doi.org/10.3389/fncel.2019.00021 PMID: 30853896

120. Wu T, Hallett M. The cerebellum in Parkinson's disease. Brain. 2013 Mar; 136(Pt 3):696–709. https://doi.org/10.1093/brain/aws360 PMID: 23404337

121. Shen B, Pan Y, Jiang X, Wu Z, Zhu J, Dong J, et al. Altered putamen and cerebellum connectivity among different subtypes of Parkinson's disease. CNS Neurosci Ther. 2020 Feb; 26(2):207–14. https://doi.org/10.1111/cns.13259 PMID: 31730272

122. Goldman JS, Tort-Colet N, di Volo M, Susin E, Bouté J, Dali M, et al. Bridging Single Neuron Dynamics to Global Brain States. Front Syst Neurosci. 2019; 13(December):1–9. https://doi.org/10.3389/fnsys.2019.00075 PMID: 31866837

123. Ruffini G, Deco G. The 2D Ising model, criticality and AIT. 2021;

124. Sadeghi S, Mier D, Gerchen MF, Schmidt SNL, Hass J. Dynamic Causal Modeling for fMRI With Wilson-Cowan-Based Neuronal Equations. Front Neurosci. 2020; 14(November). https://doi.org/10.3389/fnins.2020.593867 PMID: 33328865

125. Pinotsis DA, Moran RJ, Friston KJ. Dynamic causal modeling with neural fields. Neuroimage [Internet]. 2012; 59(2):1261–74. Available from: https://doi.org/10.1016/j.neuroimage.2011.08.020 PMID: 21924363

126. Palesi F, Tournier J-D, Calamante F, Muhler N, Castellazzi G, Chard D, et al. Contralateral cerebello-thalamo-cortical pathways with prominent involvement of associative areas in humans in vivo. Brain
Struct Funct. 2015 Nov; 220(6):3369–84. https://doi.org/10.1007/s00429-014-0861-2 PMID: 25134682

127. Wong KF, Wang XJ. A recurrent network mechanism of time integration in perceptual decisions. J Neurosci. 2006; 26(4):1314–28. https://doi.org/10.1523/JNEUROSCI.3733-05.2006 PMID: 16436619

128. Amunts K, DeFelipe J, Pennartz C, Destexhe A, Migliore M, Ryvlin P, et al. Linking Brain Structure, Activity, and Cognitive Function through Computation. eNeuro. 2022; 9(2). https://doi.org/10.1523/ENEURO.0316-21.2022 PMID: 35217544