# Complex Electroresponsive Dynamics in Olivocerebellar Neurons Represented With Extended-Generalized Leaky Integrate and Fire Models

^{1}NEARLab, Department of Electronics, Information and Bioengineering, Politecnico di Milano, Milan, Italy^{2}Department of Brain and Behavioral Sciences, University of Pavia, Pavia, Italy^{3}IRCCS Mondino Foundation, Pavia, Italy

The neurons of the olivocerebellar circuit exhibit complex electroresponsive dynamics, which are thought to play a fundamental role for network entraining, plasticity induction, signal processing, and noise filtering. In order to reproduce these properties in single-point neuron models, we have optimized the Extended-Generalized Leaky Integrate and Fire (E-GLIF) neuron through a multi-objective gradient-based algorithm targeting the desired input–output relationships. In this way, E-GLIF was tuned toward the unique input–output properties of Golgi cells, granule cells, Purkinje cells, molecular layer interneurons, deep cerebellar nuclei cells, and inferior olivary cells. E-GLIF proved able to simulate the complex cell-specific electroresponsive dynamics of the main olivocerebellar neurons including pacemaking, adaptation, bursting, post-inhibitory rebound excitation, subthreshold oscillations, resonance, and phase reset. The integration of these E-GLIF point-neuron models into olivocerebellar Spiking Neural Networks will allow to evaluate the impact of complex electroresponsive dynamics at the higher scales, up to motor behavior, in closed-loop simulations of sensorimotor tasks.

## Introduction

The variety of neuron types and spiking patterns is thought to play a fundamental role for cerebellar signal processing (Llinás, 1988, 2014) and eventually for motor learning and control. By exploiting pacemaking, bursting, adaptation and more complex properties like oscillation and resonance, cerebellar neurons can precisely encode sensorimotor signals, induce plasticity, filter noise, and efficiently communicate with different cerebellar layers and extra-cerebellar circuits (D’Angelo et al., 2016a).

The electroresponsiveness of cerebellar neurons has been deeply characterized *in vitro* and *in vivo*, allowing to identify, for each neuron type, a set of electrophysiological properties, which can be used as a reference for tuning single neuron models (Table 1). All cerebellar cortical neurons except granule cells show autorhythmic activity that becomes irregular *in vivo* due to synaptic inputs. All cerebellar neurons show an almost linear relationship between input current and firing rate, although with different slopes. In addition, the different cerebellar neurons show specific properties. The Golgi Cells (GoCs) show spike-frequency adaptation (SFA) when depolarized by prolonged currents, post-inhibitory rebound bursts, phase reset, sub-threshold oscillations (STO), and resonance in theta band (Solinas et al., 2007a,b). The granule cells (GRs) exhibit near-threshold oscillations and resonance in theta band (D’Angelo et al., 1998, 2001). The Purkinje Cells (PCs) show a discontinuous *f-I _{stim}* curve, hysteresis following current ramp stimulation and bistability emerging with high stimulus currents (intrinsic bursting) (McKay and Turner, 2005; Masoli et al., 2015; Buchin et al., 2016). Intrinsic bursting is characterized by a sequence of bursts (depolarized spiking states) and pauses (hyperpolarized quiescent states), which correlate with burst-pause responses observed

*in vivo*during behavior (Loewenstein et al., 2005). PC responses consist of simple and complex spikes: simple spikes are high-frequency regular spikes, generated spontaneously or following Parallel Fiber (PF) activation. Complex spikes consist of a burst of action potentials or spikelets, followed by a pause, resulting from Climbing Fiber (CF) excitation (Miall et al., 1998; Rokni et al., 2009). Molecular Layer Interneurons (MLIs) fire spontaneously with an increased firing irregularity

*in vivo*(Lachamp et al., 2009; Jörntell et al., 2010) and have no significant SFA (Galliano et al., 2013). These properties derive from the specific set of ionic channels and from their localization on neuronal dendrites, soma and axons, as well as from the specific nature of synaptic inputs.

The deep cerebellar nuclei cells (DCNs) express SFA and post-inhibitory rebound bursting, which is fundamental *in vivo* to modulate the motor output (Hoebeek et al., 2010; Uusisaari and Knöpfel, 2011; Ten Brinke et al., 2017). Based on the expression of marker proteins, two major types of DCN neurons have been identified, with different morphologies, electrophysiological properties, and connectivity patterns (Uusisaari et al., 2007). Large non-GABAergic DCNs (DCNnL) mainly project to pre-motor areas, adapting motor commands during learning tasks, while small GABAergic DCNs (DCNp) are connected to the Inferior Olive, providing feedback on the learning process (Uusisaari and Knöpfel, 2011).

The olivocerebellar circuit functioning strongly relies on the complex dynamics of Inferior Olive (IO) neurons. They exhibit a stereotyped response with slow STO undergoing phase-reset after impulse currents (Long et al., 2002; Kazantsev et al., 2004; Choi et al., 2010; Lefler et al., 2013). Following hyperpolarization, IO neurons generate rebound spikes (De Zeeuw et al., 2003), while when a depolarizing input is applied, single somatic action potentials are translated into bursts of axonal spikes at instantaneous frequency that can exceed 400 Hz (Maruta et al., 2007; Mathy et al., 2009). IO bursts elicit PC complex spikes and promote plasticity in the cerebellar cortex.

In this scenario, single neuron properties have been described in detailed models based on multi-compartment neurons for the different cerebellar layers (Solinas et al., 2007b; Steuber et al., 2011; De Gruijl et al., 2012; D’Angelo et al., 2013; Masoli and D’Angelo, 2017). However, representing this rich set of electroresponsive patterns through simplified neuron models is fundamental to develop realistic multiscale Spiking Neural Networks (SNNs). To tackle this issue, we here exploited the Extended-Generalized Leaky Integrate and Fire (E-GLIF) point neuron that allows to model single-point neurons while keeping a realistic picture of multiple essential electrophysiological features such as autorhythm, bursting, adaptation, oscillations, and resonance (Geminiani et al., 2018). The E-GLIF, which was originally used to reproduce the GoC electroresponsiveness (Geminiani et al., 2018), was used here to optimize and test the other cerebellar neurons: GRs, PCs, MLIs, DCNs, and IO. The results shown here are fundamental in view of SNNs simulations where the impact of complex single neuron dynamics will be evaluated at the network and, eventually, at the behavioral level (D’Angelo et al., 2016a).

## Materials and Methods

### Single Neuron Model

To reproduce the firing patterns described in the Section “Introduction,” single neurons were modeled as E-GLIF point neurons. In previous work, E-GLIF proved able to generate the complete set of GoC spiking responses to different inputs, with a minimum number of equations and free parameters. This makes it the best candidate to be used in SNNs to optimize the compromise between biological plausibility and computational load (Geminiani et al., 2018).

Extended-Generalized Leaky Integrate and Fire couples time-dependent with event-driven algorithmic components and includes three linear Ordinary Differential Equations describing the time evolution of membrane potential (*V _{m}*) and of two intrinsic currents (

*I*and

_{adapt}*I*). These three state variables are updated at spike events, which are generated according to a probabilistic threshold crossing.

_{dep}The model is defined as follows:

Where:

*I _{stim}* = external stimulation current;

*C _{m}* = membrane capacitance;

τ* _{m}* = membrane time constant;

*E _{L}* = resting potential;

*I _{e}* = endogenous current;

*k _{adap}*,

*k*

_{2}= adaptation constants;

*k*_{1} = *I _{dep}* decay rate.

If the neuron is in the refractory period *t _{ref}*, spikes cannot be emitted. Otherwise, a spike is generated stochastically at time

*t*, according to an escape rate noise: the nearer

_{spk}*V*is to the threshold potential

_{m}*V*, the higher the probability to have a spike, depending on an exponential function (Gerstner and Kistler, 2002; Jolivet et al., 2006).

_{th}At each spike event, the state variables are updated according to the rules:

Where:

${t}_{\mathrm{spk}}^{+}$ = time instant immediately following the spike time *t _{spk}*;

*V*_{r} = reset potential;

*A*_{2}, *A*_{1} = update constants of *I _{adap}* and

*I*, respectively.

_{dep}Based on *k*_{2} and *k _{adap}* values, the model exhibits exponential or oscillatory responses (Figure 1A). Elements in the model can be associated to different mechanisms that contribute to the spike patterns. The endogenous current,

*I*, accounts for autorhythm and regulation of the intrinsic steady-state membrane potential; the adaptive current,

_{e}*I*, coupled with

_{adap}*V*accounts for intrinsic sub-threshold oscillations of the membrane potential and represents the slow hyperpolarizing sub-cellular currents, e.g., the K

_{m}^{+}channel currents; the spike-triggered current,

*I*, accounts for fast depolarizing mechanisms, e.g., the Na

_{dep}^{+}and low threshold voltage activated Ca

^{2+}channel currents. For neuron connections within SNNs, conductance-based synapses are used, with spike-triggered change of synaptic conductance,

*g*, according to an alpha function (Cavallari et al., 2014; Geminiani et al., 2018):

_{syn}**Figure 1.** Optimization methods. **(A)** Parameter space for different solution regimes in E-GLIF neuron model, depending on values of parameters *k _{adap}* and

*k*

_{2}. The red line corresponds to not-damped oscillatory solutions, the red region to oscillatory damped solutions, while the green area corresponds to exponential stable solutions. Adapted from Geminiani et al. (2018).

**(B)**Stimulation protocol for evaluation of model analytical solution used for optimization in specific sub-intervals: time to first spike and between first and second spike, at the beginning of zero-/depolarizing current steps and following a hyperpolarizing step (double white arrows); time between two spikes at steady-state (white arrow); time to first spike (pause) at the end of a strong depolarizing current step, exc

_{3}, only for PC E-GLIF optimization, to fit the burst-pause response (black arrow).

where *G _{syn}* is the maximum conductance change and τ

*the synaptic time constant.*

_{syn}### Neuron Model Optimization

Analogously to the GoC E-GLIF optimization, for each cerebellar neuron we derived the parameters related to neurophysiological quantities (i.e., *C _{m}*, τ

*,*

_{m}*E*, Δ

_{L}*t*,

_{ref}*V*,

_{th}*V*) from literature

_{r}*in vitro*experiments (Table 2). For the remaining parameters (i.e.,

*k*,

_{adap}*k*

_{2},

*k*

_{1},

*A*

_{2},

*A*

_{1},

*I*), we used the optimization strategy described in Geminiani et al. (2018), developed in MATLAB, where the cost and constraint functions were adapted to consider the electroresponsive properties of each neuron type as in Table 1.

_{e}**Table 2.** Electrophysiological passive properties chosen from literature for the different cerebellar neurons.

#### Optimization Stimulation Protocol

Exploiting the analytical solution of the model, the optimization algorithm aimed at minimizing the error on spike times during three sub-intervals of a current step stimulation period, where the *V _{m}* solution could be computed: the time to the first spike, the time between first and second spike and the time between two steady-state spikes (Figure 1B). A multi-step stimulation protocol was considered for optimization, including: a zero-current phase, three phases with increasing depolarizing currents (

*exc*

_{1}<

*exc*

_{2}<

*exc*

_{3}), and a zero-current phase following a stimulation interval with a negative current,

*inh*.

#### Cost Function

The cost function evaluated the error on the desired spike times (computed from desired output frequency), in order to fit cell-specific quantitative input–output relationships (Supplementary Table S1): (i) autorhythm frequency, when *I _{stim}* = 0, (ii) response rates (

*freq*

_{1}<

*freq*

_{2}<

*freq*

_{3}), with increasing amplitudes of

*I*(

_{stim}*exc*

_{1}<

*exc*

_{2}<

*exc*

_{3}), and (iii) rebound burst latency and initial frequency, following an inhibitory current step,

*inh*. To take into account SFA during depolarizing current steps, the desired steady-state firing rate was obtained from desired frequencies (

*freq*

_{1}<

*freq*

_{2}<

*freq*

_{3}) multiplied by an attenuating factor (

*factor*

_{1},

*factor*

_{2},

*factor*

_{3}) based on experimental values.

In addition, only PCs exhibit the burst-pause response (Masoli et al., 2015): to account for this specific property, the PC cost function evaluated also the time to the first spike (i.e., the pause), just after the turning off of *I _{stim}* =

*exc*

_{3}(Figure 1B).

#### Optimization Constraints

The cell-specific constraints (Supplementary Table S2) were customized to obtain:

• Neurophysiological ranges of currents in the model;

• Neurophysiological steady-state value of the membrane potential during inhibition (*V _{m_inh}*);

• Oscillatory damped or not (red area in Figure 1A) or exponential (green area in Figure 1A) *V _{m}* dynamics (Geminiani et al., 2018), based on

*k*

_{2}and

*k*ranges as in Figure 1A;

_{adap}• Neurophysiological values of oscillation frequency, in case of oscillatory neurons, i.e., GRs and IOs;

• Sub-threshold value of the steady-state membrane potential (*V _{m_ss_tonic}*) and limited amplitude of oscillations (

*A*) to prevent spontaneous firing in oscillatory neurons without autorhythm – GRs and IOs, in case of zero external input.

_{osc_tonic}The mathematical expression of the cost function, the fitted input-output quantitative patterns and the values of the constraints are reported with proper details in Supplementary Material.

#### Optimization Implementation

For each neuron type, we ran five optimizations with different random initializations of parameters within their ranges, to test the robustness of results with respect to initialization. We chose the optimal parameter set as the median of the final parameters in each optimization run.

### Neuron Model Validation

To validate the outcome of optimization and test the effective proper functioning of the model based on literature data, we simulated the E-GLIF responses during a continuous stimulation protocol with current steps in PyNEST (Diesmann and Gewaltig, 2002). This validation was fundamental to assess the result of optimization that was based on the evaluation of the neuron response only in sub-sampled intervals of a continuous simulation. In order to evaluate all the electroresponsive properties in Table 1, the stimulation protocol included a first phase with zero external current, where to measure autorhythm and irregular firing, followed by three depolarizing phases lasting 1 s and interleaved with 1-s zero-current intervals, to measure intrinsic excitability and adaptation. Afterward, a 1-s inhibitory current was applied and turned off in the subsequent step, to test rebound bursting (Figure 2A, left panel). The amplitudes of current steps in each phase were the same used during optimization, but the whole continuous response was here assessed, and not just the sub-intervals included in the optimization. The stimulation protocol was then customized with additional or modified phases for neurons with specific electroresponsive patterns:

**Figure 2.** Stimulation protocol for E-GLIF model validation in PyNEST simulations. **(A)** General *in vitro* protocol with the three depolarizing current steps (exc_{1,2,3}) and the inhibitory step (inh) used for PyNEST simulations of MLI and DCN E-GLIF (left panel); a shorter exc_{3} current step is used for PCs to test the burst-pause response (right panel). The current amplitude values are the same used in the optimization process, where only sub-intervals of each stimulation phase were considered. **(B)** Customized protocol for GR E-GLIF to test resonance through a stimulation phase with periodic spike trains at increasing frequencies. **(C)** Customized protocol for IO E-GLIF with one shorter depolarizing step and an impulse stimulus to evaluate phase reset of membrane potential oscillations.

• For PCs, we reduced the third depolarizing interval from 1 to 0.01 s (Figure 2A, right panel) to test the burst-pause response with high input currents (McKay and Turner, 2005) and evaluate the effect of current pulses (analogous to CF bursts);

• For GRs, we included an additional phase with input current step trains at increasing frequencies (0.3-3-6-9-12-15 Hz), to evaluate resonance (Figure 2B);

• For IOs, we considered only one depolarizing phase lasting 0.05 s, to adapt to literature reference protocols for *in vitro* experiments. Then, we tested the effect of different current amplitudes on burst response properties and we evaluated phase reset of STO, following a current impulse (amplitude = 1 nA, duration = 5 ms), during a zero-current interval lasting 1.5 s (Figure 2C).

We ran 10 simulations for each neuron and computed the mean ± Standard Deviation (SD) of activity parameters (see section “Validation Data Analysis”).

### Validation Data Analysis

Significant parameters were extracted from spiking time instants to evaluate single neuron firing patterns in validation protocols:

• The tonic firing rate, *f _{tonic}*, as the inverse of the mean inter-spike interval (ISI), and the coefficient of variation of inter-spike intervals (CV

_{ISI}) to quantify the irregularity of firing, during the zero-current phase;

• The firing rate, *f*, as the inverse of the mean ISI, during the first three spikes of each depolarizing phase;

• The steady-state firing rate, *f _{ss}*, as the inverse of the mean ISI, during the last six spikes of each 1-s depolarizing phase;

• The *f-I _{stim}* slope derived from initial responses to the excitatory step currents;

• The SFA gain, computed as the ratio between *f* and *f _{ss}*;

• Latency and initial frequency (i.e., inverse of the first burst ISI), measured in the rebound burst after hyperpolarization (*lat_rebound* and *rebound_freq*, respectively). Post-inhibitory activity was considered a rebound burst if *lat_rebound* and *rebound_freq* were lower than the autorhythm ISI and higher than the autorhythm frequency, respectively.

To quantify resonance in GRs, we also computed the response speed as the inverse of the mean spike latency in each resonance step; the values from multiple simulation tests and frequencies were fitted through a smoothing spline in order to obtain the resonance curve (Gandolfi et al., 2013).

## Results

The single-point models of cerebellar neurons were generated using E-GLIF protocol (Geminiani et al., 2018) and were tuned toward their specific neurophysiological response patterns. For GoCs, we used the same optimal parameters reported in Geminiani et al. (2018). For the other neurons, after fixing the passive properties from literature data (Table 2), the optimization algorithm was used to tune the remaining model parameters toward specific electrophysiological features. In most cases, the algorithm converged to the same region of the parameter space over the five optimization runs (Supplementary Figures S1, S2). The resulting parameter sets achieved the optimal compromise between minimum cost function and constraint violation (below 1.0 and 0.1, respectively), best reproducing the electroresponsiveness of each neuron type (Table 3).

Tuned E-GLIF neurons were then tested in PyNEST simulations with the stimulation protocol described in the Section “Neuron Model Validation.” The model was able to capture the intrinsic excitability of all neurons, generating linearly increasing firing rates with depolarizing current steps. As shown in Figure 3, frequencies values and *f-I _{stim}* slope were close to the target values for all neurons or within acceptable ranges. For GRs the

*f-I*slope was lower than in the reference study (D’Angelo et al., 1998) but still consistent with experimental ranges (Spanne et al., 2014; Masoli et al., 2017). In DCNnL, depolarization frequencies were higher than target values, but linearly increasing with an acceptable

_{stim}*f-I*slope (Table 4). SFA was present for PCs and DCNnL with average SFA gain of 1.1 at all

_{stim}*I*values, close to the target values of adaptation gain from electrophysiological recordings (1.1 and 1.2, respectively) (Uusisaari et al., 2007; Kim et al., 2013). In DCNp, SFA was more pronounced, with an average gain of 1.3 for I

_{stim}_{stim}=

*exc*

_{2,3}(Uusisaari et al., 2007). In absence of external stimuli, PC, MLI and both DCN E-GLIF produced irregular autorhythm at physiological frequencies, while GRs and IOs generated STO at 6 and 7 Hz, respectively (Figure 4). At the end of a hyperpolarizing current step, PCs and DCNs exhibited rebound excitation (doublets/bursts), which is fundamental for efficient signal transmission (Figure 4). In IOs, post-inhibitory rebound spikes were generated with 50% probability, as in experiments (De Zeeuw et al., 2003; Mathy et al., 2009). When stimulating PC with current pulses of 2.4 nA, the typical intrinsic bursting (burst-pause response) was generated. This was achieved thanks to the balance of model currents,

*I*and

_{dep}*I*that accounted for subcellular mechanisms leading to PC complex spikes (De Zeeuw et al., 2011). A 10-ms pulse caused a burst at 254.58 ± 18.26 Hz followed by a pause of 23.47 ± 2.38 ms, longer than the tonic ISI (Figure 5A); with a 50-ms current step the neuron was silent for 32.46 ± 1.22 ms after a burst at 234.87 ± 2.70 Hz (Figure 5B; Grasselli et al., 2016). This spiking pattern well fits with the PC response to dendritic current injection; however, the typical PC bistable regime caused by a continuous high-amplitude stimulation could not be reproduced in the model without losing other electroresponsive properties (Masoli et al., 2015). Intrinsic STO in GRs lead to resonance at 6 Hz, when stimulating the GR neuron model with periodic spike trains at increasing frequencies (Figure 6A). Finally, the optimized E-GLIF model was able to generate also the typical IO bursting response (193.91 ± 24.58 Hz) in case of current step input, thanks to the rapid effect of

_{adap}*I*at the beginning of stimulation and the slower accumulation of

_{dep}*I*that blocked the firing (Figure 5B). Increased amplitudes of the input current caused a non-linear increase of the burst frequency, within physiological ranges; instead, lower currents (i.e., 200 pA) were not sufficient to activate bursts, but they only produced single spikes followed by a pause. Current pulses in the IO E-GLIF induced a spike and a subsequent phase reset of STO, independent from the phase of the stimulus (Figure 6B). Consistently with experimental results, post-impulse STO phase in the model was (0.87 ± 0.02)⋅T for pre-stimulus phases ranging from 0.06⋅T to 0.92⋅T, being T the period of oscillations (Kazantsev et al., 2004; Lefler et al., 2013).

_{adap}**Figure 3.** Plots of *f-I _{stim}* relationships for GRs, PCs, MLIs, and DCNnL neurons, comparing outcome of PyNEST simulations (black markers) with literature target values used for optimization (blue markers), at the beginning (circles), and after 1-s (squares) current step stimulation. Experimental data taken from D’Angelo et al. (2001), McKay and Turner (2005), Uusisaari et al. (2007), and Galliano et al. (2013).

**Figure 4.** E-GLIF responses to zero-input current (left column) and following a hyperpolarizing current step (right column) for the main olivo-cerebellar neurons. Zero-current inputs cause STO in GR and IO neurons, and autorhythm in the others. Post-inhibitory rebound excitation (burst or spike) is highlighted in the blue circle, where present.

**Figure 5.** Bursting responses in E-GLIF simulations. **(A)** Burst-pause in PC E-GLIF with a 10-ms input current step (left panel) and 50-ms input current step (right panel). *V _{m}* and input current traces are reported in top panels, showing the burst during the stimulation phase and the subsequent pause (blue segment) when the current goes back to 0 nA. Model current traces are reported in bottom panels, with respect to their steady-state value (ΔI). The I

_{adap}current is reported in negative values as it has a hyperpolarizing effect in the neuron model. At the end of the stimulation, the accumulated inhibitory effect of I

_{adap}causes the pause, until it decays, and the tonic balance of currents is restored.

**(B)**Bursting response in IO E-GLIF during a 50-ms current step stimulus, showing a first doublet (zoom in the inset) followed by a pause (blue segment); even in this case, the intrinsic model currents drive the

*V*response (bottom panel).

_{m}**Figure 6.** Oscillation-driven properties in E-GLIF simulations. **(A)** Resonance in GR E-GLIF following spike train stimulation at increasing frequencies. The resonance curve was obtained fitting the data points from 10 simulations (black dots) with a smoothing spline. **(B)** Phase-reset of IO E-GLIF *V _{m}* during 10 simulations with the same current pulse, causing a spike and a subsequent reset of oscillation phase, independent from the phase before the stimulus.

Therefore, the whole set of olivo-cerebellar cells could be modeled with E-GLIF neurons, generating realistic spiking patterns and capturing crucial electroresponsive properties for cerebellar functioning.

## Discussion

In this paper, the E-GLIF model (Geminiani et al., 2018), that was previously developed and validated for Golgi cells, was tuned toward the unique electroresponsive properties of granule cells, Purkinje cells, molecular layer interneurons, deep cerebellar nuclei cells and inferior olivary cells. In these neurons, E-GLIF effectively reproduced pacemaking, adaptation, bursting, post-inhibitory rebound excitation, subthreshold oscillations, resonance, and phase reset. Therefore, for the first time, a whole set of single point neurons is made available to investigate the functional dynamics of the olivocerebellar circuit (Voogd and Glickstein, 1998; Ruigrok, 2011; D’Angelo et al., 2013; Witter et al., 2013; Zhou et al., 2014). These include oscillations and resonance, which are thought to play a critical role for network entraining into large-scale brain oscillations (De Zeeuw et al., 2011; Courtemanche et al., 2013; Llinás, 2014), and long-term synaptic plasticity, which is considered the main mechanism underlying the cerebellar role in motor control and learning (Ito et al., 2014; D’Angelo et al., 2016b).

### Modeled Single Neuron Dynamics

Extended-Generalized Leaky Integrate and Fire (Geminiani et al., 2018) is a simplified point-neuron based on a system of three linear ordinary differential equations and its analytical tractability allows to define different solution regimes and to tune model parameters through a generalizable optimization algorithm. In the current work, E-GLIF was able to simulate complex input-output relationships of cerebellar and IO neurons, generating cell-specific intrinsic excitability and non-linear firing properties that would not be possible using previous GLIF models (Mihalaş and Niebur, 2009).

For neurons with oscillatory *V _{m}*, the second order dynamics of the model allowed to simulate intrinsic self-sustained STO. Second order dynamics allowed to reproduce also other non-linear electroresponsive behaviors like resonance in GRs and phase reset of STO in IO neurons. These properties have been measured in single-neuron experiments and are probably amplified at network level (D’Angelo et al., 2001). Specifically, the feedback inhibitory loop from GoCs to GRs is supposed to contribute to resonance and oscillations in the Granular layer network, enhancing theta-band signals coming from extra-cerebellar regions (D’Angelo and Casali, 2013; Gandolfi et al., 2013). Future simulations of the granular layer network with E-GLIF neurons will help to elucidate the different contribution of single cell and circuit properties on network oscillations and resonance. This would extend the results of previous studies where detailed microcircuit models and SNNs with Leaky Integrate-and-Fire units were exploited (D’Angelo et al., 2013; Casali et al., 2019). In the IO circuit, phase reset of STO has been measured in single neurons (Kazantsev et al., 2004), but synchronous stimulation of an olivary area was shown to amplify this response (Lefler et al., 2013). The IO E-GLIF could reproduce the first response during simulation of

*in vitro*protocols. In principle, adding gap junctions to the neuron model would account also for the phase-reset amplification at network level, thanks to the intrinsic communication within IO nuclei.

To simulate IO neurons, E-GLIF was optimized taking the axonal bursting regime as the target behavior (Maruta et al., 2007; Mathy et al., 2009). This aspect challenges the traditional view of CFs as a low-frequency all-or-none signaling pathway: indeed, bursting and rebound activity in IO is fundamental for information encoding, as rebound excitation amplifies the feedback from DCNp cells and olivary bursts elicit complex spikes at PC level. PC E-GLIF successfully reproduced regular firing and the burst-pause pattern following dendritic current stimulation *in vitro*, which can be associated to simple and complex spikes *in vivo* (Masoli et al., 2015). However, bistability and spiking patterns with longer bursts and pauses could not be obtained in the E-GLIF model without losing intrinsic excitability properties. For simulations in SNNs, this is a sufficient approximation since it allows to generate the typical PC network spiking patterns, as shown in the Section “Results.” However, for a more detailed representation even of axonal responses, a multi-compartment version of the PC E-GLIF could be implemented, where multiple E-GLIF neurons are optimized to reproduce the electroresponsiveness of the main PC compartments.

In cerebellar nuclei neurons, rebound excitation has been widely proven *in vitro* but long debated *in vivo* (Alviña et al., 2008). However, recent experimental findings demonstrate that rebound bursting correlates with motor responses and is fundamental for integrating synaptic inputs from PCs, MFs, and IO neurons that all converge in the cerebellar nuclei (Hoebeek et al., 2010; Manto and Oulad Ben Taib, 2010; Witter et al., 2013; Sarnaik and Raman, 2018). Rebound excitation also contributes to cerebellum-driven learning, as demonstrated for associative learning (Ten Brinke et al., 2017). Single-neuron rebound properties are thus crucial in SNNs aimed at multiscale simulations of sensorimotor tasks.

This scenario shows the capability of the E-GLIF point neuron to reproduce the variety of olivo-cerebellar spiking responses following different input stimuli, through a single optimal set of model parameters. Conversely, the traditional approach for single neuron modeling aims at identifying different regions of the parameter space corresponding to different spiking behaviors (Izhikevich, 2003). This makes E-GLIF a best candidate for simulations of SNNs, where neuron response needs to depend on the received input, rather than on the parameter values, achieving higher neurophysiological realism without increasing computational load.

## Conclusion

The E-GLIF single-point neuron models were able to capture the complex non-linear dynamics of olivocerebellar neurons including spontaneous firing, subthreshold oscillations, bursting, phase-reset, and resonance. These ingredients, coupled to algorithms accounting for synaptic integration over dendrites (e.g., Marasco et al., 2012; Rössert et al., 2016), will provide the fundamental ingredients to reconstruct non-linear dynamics in extended spiking cerebellar networks. Future work will include embedding these neuron models into cerebellar SNNs to simulate cerebellum-driven motor paradigms and evaluate the impact of single neuron electroresponsiveness on network dynamics, plasticity and, eventually, motor behavior.

## Data Availability

All datasets generated for this study are included in the manuscript and/or the Supplementary Files.

## Author Contributions

AG and CC elaborated the mathematical model and optimization, designed and carried out the simulations for each neuron, performed the data analysis, and wrote the manuscript. ED and AP coordinated the whole work and substantially contributed to the writing of the final manuscript.

## Funding

This project has been developed within the *CerebNEST* HBP Partnering Project and has received funding from the European Union’s Horizon 2020 Framework Programme for Research and Innovation under Grant Agreement No. 785907 (Human Brain Project SGA2).

## Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fncom.2019.00035/full#supplementary-material

## References

Aizenman, C. D., and Linden, D. J. (1999). Regulation of the rebound depolarization and spontaneous firing patterns of deep nuclear neurons in slices of rat cerebellum. *J. Neurophysiol.* 82, 1697–1709. doi: 10.1152/jn.1999.82.4.1697

Alviña, K., Walter, J. T., Kohn, A., Ellis-Davies, G., and Khodakhah, K. (2008). Questioning the role of rebound firing in the cerebellum. *Nat. Neurosci.* 11, 1256–1258. doi: 10.1038/nn.2195

Buchin, A., Rieubland, S., Häusser, M., Gutkin, B. S., and Roth, A. (2016). Inverse stochastic resonance in cerebellar purkinje cells. *PLoS Comput. Biol.* 12:e1005000. doi: 10.1371/journal.pcbi.1005000

Casali, S., Marenzi, E., Medini, C., Casellato, C., and D’Angelo, E. (2019). Reconstruction and simulation of a scaffold model of the cerebellar network. *Front. Neuroinform.* 13:37. doi: 10.3389/fninf.2019.00037

Cavallari, S., Panzeri, S., and Mazzoni, A. (2014). Comparison of the dynamics of neural interactions between current-based and conductance-based integrate-and-fire recurrent networks. *Front. Neural Circ.* 8:12. doi: 10.3389/fncir.2014.00012

Choi, S., Yu, E., Kim, D., Urbano, F. J., Makarenko, V., Shin, H. S., et al. (2010). Subthreshold membrane potential oscillations in inferior olive neurons are dynamically regulated by P/Q- and T-type calcium channels: a study in mutant mice. *J. Physiol.* 588, 3031–3043. doi: 10.1113/jphysiol.2009.184705

Courtemanche, R., Robinson, J. C., and Aponte, D. I. (2013). Linking oscillations in cerebellar circuits. *Front. Neural Circuits* 7:125. doi: 10.3389/fncir.2013.00125

D’Angelo, E., Antonietti, A., Casali, S., Casellato, C., Garrido, J. A., Luque, N. R., et al. (2016a). Modeling the cerebellar microcircuit: new strategies for a long-standing issue. *Front. Cell. Neurosci.* 10:176. doi: 10.3389/fncel.2016.00176

D’Angelo, E., Mapelli, L., Casellato, C., Garrido, J. A., Luque, N., Monaco, J., et al. (2016b). Distributed circuit plasticity: new clues for the cerebellar mechanisms of learning. *Cereb.* 15, 139–151. doi: 10.1007/s12311-015-0711-7

D’Angelo, E., and Casali, S. (2013). Seeking a unified framework for cerebellar function and dysfunction: from circuit operations to cognition. *Front. Neural Circuits* 6:116. doi: 10.3389/fncir.2012.00116

D’Angelo, E., De Filippi, G., Rossi, P., and Taglietti, V. (1998). Ionic mechanism of electroresponsiveness in cerebellar granule cells implicates the action of a persistent sodium current. *J. Neurophysiol.* 80, 493–503. doi: 10.1152/jn.1998.80.2.493

D’Angelo, E., Nieus, T., Maffei, A., Armano, S., Rossi, P., Taglietti, V., et al. (2001). Theta-frequency bursting and resonance in cerebellar granule cells: experimental evidence and modeling of a slow K-dependent mechanism. *J. Neurosci.* 21, 759–770. doi: 10.1523/jneurosci.21-03-00759.2001

D’Angelo, E., Solinas, S., Mapelli, J., Gandolfi, D., Mapelli, L., and Prestori, F. (2013). The cerebellar Golgi cell and spatiotemporal organization of granular layer activity. *Front. Neural Circuits* 7:93. doi: 10.3389/fncir.2013.00093

De Gruijl, J. R., Bazzigaluppi, P., de Jeu, M. T. G., and De Zeeuw, C. I. (2012). Climbing fiber burst size and olivary sub-threshold oscillations in a network setting. *PLoS Comput. Biol.* 8:e1002814. doi: 10.1371/journal.pcbi.1002814

De Schutter, E., and Steuber, V. (2009). Patterns and pauses in Purkinje cell simple spike trains: experiments, modeling and theory. *Neuroscience* 162, 816–826. doi: 10.1016/j.neuroscience.2009.02.040

De Zeeuw, C. I., Chorev, E., Devor, A., Manor, Y., Van Der Giessen, R. S., De Jeu, M. T., et al. (2003). Deformation of network connectivity in the inferior olive of connexin 36-deficient mice is compensated by morphological and electrophysiological changes at the single neuron level. *J. Neurosci.* 23, 4700–4711. doi: 10.1523/JNEUROSCI.23-11-04700.2003

De Zeeuw, C. I., Hoebeek, F. E., Bosman, L. W. J., Schonewille, M., Witter, L., and Koekkoek, S. K. (2011). Spatiotemporal firing patterns in the cerebellum. *Nat. Rev. Neurosci.* 12, 327–344. doi: 10.1038/nrn3011

Diesmann, M., and Gewaltig, M.-O. (2002). “NEST: An environment for neural systems simulations,” in *Forsch. und Wisschenschaftliches Rechn. Beitrage Zum Heinz-billing-pr. 2001*, Vol. 58, eds T. Plesser and V. Macho (Göttingen: GWDG-Bericht), 43–70.

Forti, L., Cesana, E., Mapelli, J., and Angelo, E. D. (2006). Ionic mechanisms of autorhythmic firing in rat cerebellar Golgi cells. *J. Physiol.* 3, 711–729. doi: 10.1113/jphysiol.2006.110858

Galliano, E., Gao, Z., Schonewille, M., Todorov, B., Simons, E., Pop, A. S., et al. (2013). Silencing the majority of cerebellar granule cells uncovers their essential role in motor learning and consolidation. *Cell Rep.* 3, 1239–1251. doi: 10.1016/j.celrep.2013.03.023

Gandolfi, D., Lombardo, P., Mapelli, J., Solinas, S., and D’Angelo, E. (2013). Theta-frequency resonance at the cerebellum input stage improves spike timing on the millisecond time-scale. *Front. Neural Circuits* 7:64. doi: 10.3389/fncir.2013.00064

Geminiani, A., Casellato, C., Locatelli, F., Prestori, F., Pedrocchi, A., and D’Angelo, E. (2018). Complex dynamics in simplified neuronal models: reproducing Golgi cell electroresponsiveness. *Front. Neuroinform.* 12:88. doi: 10.3389/fninf.2018.00088

Grasselli, G., He, Q., Wan, V., Adelman, J. P., Ohtsuki, G., and Hansel, C. (2016). Activity-dependent plasticity of spike pauses in cerebellar purkinje cells. *Cell Rep.* 14, 2546–2553. doi: 10.1016/j.celrep.2016.02.054

Hoebeek, F. E., Witter, L., Ruigrok, T. J. H., and De Zeeuw, C. I. (2010). Differential olivo-cerebellar cortical control of rebound activity in the cerebellar nuclei. *Proc. Natl. Acad. Sci. U.S.A.* 107, 8410–8415. doi: 10.1073/PNAS.0907118107

Hourez, R., Servais, L., Orduz, D., Gall, D., Millard, I., de Kerchove d’Exaerde, A., et al. (2011). Aminopyridines correct early dysfunction and delay neurodegeneration in a mouse model of spinocerebellar ataxia type 1. *J. Neurosci.* 31, 11795–11807. doi: 10.1523/JNEUROSCI.0905-11.2011

Houston, C. M., Diamanti, E., Diamantaki, M., Kutsarova, E., Cook, A., Sultan, F., et al. (2017). Exploring the significance of morphological diversity for cerebellar granule cell excitability. *Sci. Rep.* 7:46147. doi: 10.1038/srep46147

Hoxha, E., Boda, E., Montarolo, F., Parolisi, R., and Tempia, F. (2012). Excitability and synaptic alterations in the cerebellum of APP/PS1 mice. *PLoS One* 7:726. doi: 10.1371/journal.pone.0034726

Ito, M., Yamaguchi, K., Nagao, S., and Yamazaki, T. (2014). Long-term depression as a model of cerebellar plasticity. *Prog. Brain Res.* 210, 1–30. doi: 10.1016/B978-0-444-63356-9.00001-7

Izhikevich, E. M. (2003). Simple model of spiking neurons. *IEEE Trans. Neural Netw.* 14, 1569–1572. doi: 10.1109/TNN.2003.820440

Jolivet, R., Rauch, A., Lüscher, H. R., and Gerstner, W. (2006). Predicting spike timing of neocortical pyramidal neurons by simple threshold models. *J. Comput. Neurosci.* 21, 35–49. doi: 10.1007/s10827-006-7074-5

Jörntell, H., Bengtsson, F., Schonewille, M., and De Zeeuw, C. I. (2010). Cerebellar molecular layer interneurons - computational properties and roles in learning. *Trends Neurosci.* 33, 524–535. doi: 10.1016/j.tins.2010.08.004

Kazantsev, V. B., Nekorkin, V. I., Makarenko, V. I., and Llinás, R. (2004). Self-referential phase reset based on inferior olive oscillator dynamics. *Proc. Natl. Acad. Sci. U.S.A.* 101, 18183–18188. doi: 10.1073/pnas.0407900101

Kim, Y. S., Kang, E., Makino, Y., Park, S., Shin, J. H., Song, H., et al. (2013). Characterizing the conductance underlying depolarization-induced slow current in cerebellar Purkinje cells. *J. Neurophysiol.* 109, 1174–1181. doi: 10.1152/jn.01168.2011

Lachamp, P. M., Liu, Y., and Liu, S. J. (2009). Glutamatergic modulation of cerebellar interneuron activity is mediated by an enhancement of GABA release and requires protein kinase A/RIM1 signaling. *J. Neurosci.* 29, 381–392. doi: 10.1523/JNEUROSCI.2354-08.2009

Lefler, Y., Torben-nielsen, B., and Yarom, Y. (2013). Oscillatory activity, phase differences, and phase resetting in the inferior olivary nucleus. *Front. Syst. Neurosci.* 7:22. doi: 10.3389/fnsys.2013.00022

Lennon, W., Hecht-nielsen, R., and Yamazaki, T. (2014). A spiking network model of cerebellar Purkinje cells and molecular layer interneurons exhibiting irregular firing. *Front. Comput. Neurosci.* 8:157. doi: 10.3389/fncom.2014.00157

Llinás, R. (1988). The intrinsic electrophysiological properties of mammalian neurons: insights into central nervous system function. *Science* 242, 1654–1664. doi: 10.1126/science.3059497

Llinás, R. R. (2014). The olivo-cerebellar system: a key to understanding the functional significance of intrinsic oscillatory brain properties. *Front. Neural Circuits* 7:96. doi: 10.3389/fncir.2013.00096

Loewenstein, Y., Mahon, S., Chadderton, P., Kitamura, K., Sompolinsky, H., Yarom, Y., et al. (2005). Bistability of cerebellar Purkinje cells modulated by sensory stimulation. *Nat. Neurosci.* 8, 202–211. doi: 10.1038/nn1393

Long, M. A., Deans, M. R., Paul, D. L., and Connors, B. W. (2002). Rhythmicity without synchrony in the electrically uncoupled inferior olive. *J. Neurosci.* 22, 10898–10905. doi: 10.1523/jneurosci.22-24-10898.2002

Manto, M., and Oulad Ben Taib, N. (2010). Cerebellar Nuclei: key roles for strategically located structures. *Cerebellum* 9, 17–21. doi: 10.1007/s12311-010-0159-8

Marasco, A., Limongiello, A., and Migliore, M. (2012). Fast and accurate low-dimensional reduction of biophysically detailed neuron models. *Sci. Rep.* 2, 1–7. doi: 10.1038/srep00928

Maruta, J., Hensbroek, R. A., and Simpson, J. I. (2007). Intraburst and interburst signaling by climbing fibers. *J. Neurosci.* 27, 11263–11270. doi: 10.1523/JNEUROSCI.2559-07.2007

Masoli, S., and D’Angelo, E. (2017). Synaptic activation of a detailed purkinje cell model predicts voltage-dependent control of burst-pause responses in active dendrites. *Front. Cell. Neurosci.* 11:278. doi: 10.3389/fncel.2017.00278

Masoli, S., Rizza, M. F., Sgritta, M., Van Geit, W., Schürmann, F., and D’Angelo, E. (2017). Single neuron optimization as a basis for accurate biophysical modeling: the case of cerebellar granule cells. *Front. Cell. Neurosci.* 11:14. doi: 10.3389/fncel.2017.00071

Masoli, S., Solinas, S., and Angelo, E. D. (2015). Action potential processing in a detailed Purkinje cell model reveals a critical role for axonal compartmentalization. *Front. Cell. Neurosci.* 9:47. doi: 10.3389/fncel.2015.00047

Mathy, A., Ho, S. S. N., Davie, J. T., Duguid, I. C., Beverley, A. C., and Häusser, M. (2009). Encoding of oscillations by axonal bursts in inferior olive neurons. *Neuron* 62, 388–399. doi: 10.1016/j.neuron.2009.03.023

McKay, B. E., and Turner, R. W. (2005). Physiological and morphological development of the rat cerebellar Purkinje cell. *J. Physiol.* 567, 829–850. doi: 10.1113/jphysiol.2005.089383

Miall, R. C., Keating, J. G., Malkmus, M., and Thach, W. T. (1998). Simple spike activity predicts occurrence of complex spikes in cerebellar Purkinje cells. *Nat. Neurosci.* 1, 13–15. doi: 10.1038/212

Mihalaş, Ş., and Niebur, E. (2009). A generalized linear integrate-and-fire neural model produces diverse spiking behaviors. *Neural Comput.* 21, 704–718. doi: 10.1162/neco.2008.12-07-680

Molineux, M. L., McRory, J. E., McKay, B. E., Hamid, J., Mehaffey, W. H., Rehak, R., et al. (2006). Specific T-type calcium channel isoforms are associated with distinct burst phenotypes in deep cerebellar nuclear neurons. *Proc. Natl. Acad. Sci.* 103, 5555–5560. doi: 10.1073/pnas.0601261103

Rokni, D., Tal, Z., Byk, H., and Yarom, Y. (2009). Regularity, variabilty and bi-stability in the activity of cerebellar Purkinje cells. *Front. Cell. Neurosci.* 3:12. doi: 10.3389/neuro.03.012.2009

Rössert, C., Pozzorini, C., Chindemi, G., Davison, A. P., Eroe, C., King, J., et al. (2016). Automated point-neuron simplification of data-driven microcircuit models. *arxiv*Google Scholar

Ruigrok, T. J. H. (2011). Ins and outs of cerebellar modules. *Cerebellum* 10, 464–474. doi: 10.1007/s12311-010-0164-y

Sarnaik, R., and Raman, I. M. (2018). Control of voluntary and optogenetically perturbed locomotion by spike rate and timing of neurons of the mouse cerebellar nuclei. *eLife* 7:e29546. doi: 10.7554/eLife.29546

Solinas, S., Forti, L., Cesana, E., Mapelli, J., De Schutter, E., and D’Angelo, E. (2007a). Computational reconstruction of pacemaking and intrinsic electroresponsiveness in cerebellar golgi cells. *Front. Cell. Neurosci.* 1:1. doi: 10.3389/ne

Solinas, S., Forti, L., Cesana, E., Mapelli, J., De Schutter, E., and D’Angelo, E. (2007b). Fast-reset of pacemaking and theta-frequency resonance patterns in cerebellar golgi cells: simulations of their impact in vivo. *Front. Cell. Neurosci.* 1:4. doi: 10.3389/neuro.03.004.2007

Spanne, A., Geborek, P., Bengtsson, F., and Jörntell, H. (2014). Simulating spinal border cells and cerebellar granule cells under locomotion – a case study of spinocerebellar information processing. *PLoS One* 9:e107793. doi: 10.1371/journal.pone.0107793

Steuber, V., Schultheiss, N. W., Silver, R. A., De Schutter, E., and Jaeger, D. (2011). Determinants of synaptic integration and heterogeneity in rebound firing explored with data-driven models of deep cerebellar nucleus cells. *J. Comput. Neurosci.* 30, 633–658. doi: 10.1007/s10827-010-0282-z

Ten Brinke, M. M., Heiney, S. A., Wang, X., Proietti-Onori, M., Boele, H. J., Bakermans, J., et al. (2017). Dynamic modulation of activity in cerebellar nuclei neurons during pavlovian eyeblink conditioning in mice. *eLife* 6, 1–27. doi: 10.7554/eLife.28132

Tripathy, S. J., Savitskaya, J., Burton, S. D., Urban, N. N., and Gerkin, R. C. (2014). NeuroElectro: a window to the world’s neuron electrophysiology data. *Front. Neuroinform.* 8:40. doi: 10.3389/fninf.2014.00040

Uusisaari, M., and Knöpfel, T. (2011). Functional classification of neurons in the mouse lateral cerebellar nuclei. *Cerebellum* 10, 637–646. doi: 10.1007/s12311-010-0240-3

Uusisaari, M., Obata, K., Kno, T., and Morphological, T. (2007). Morphological and electrophysiological properties of GABAergic and Non-GABAergic cells in the deep cerebellar nuclei. *J. Neurophysiol.* 97, 901–911. doi: 10.1152/jn.00974.2006

Van Der Giessen, R. S., Koekkoek, S. K., van Dorp, S., De Gruijl, J. R., Cupido, A., Khosrovani, S., et al. (2008). Role of olivary electrical coupling in cerebellar motor learning. *Neuron* 58, 599–612. doi: 10.1016/j.neuron.2008.03.016

Voogd, J., and Glickstein, M. (1998). The anatomy of the cerebellum. *Trends Cogn. Sci.* 2, 307–313. doi: 10.1016/0166-2236(95)93903-B

Witter, L., Canto, C. B., Hoogland, T. M., de Gruijl, J. R., and De Zeeuw, C. I. (2013). Strength and timing of motor responses mediated by rebound firing in the cerebellar nuclei after Purkinje cell activation. *Front. Neural Circuits* 7:133. doi: 10.3389/fncir.2013.00133

Keywords: neuronal modeling, point neuron, neuron model simplification, neuronal electroresponsiveness, olivocerebellar neurons

Citation: Geminiani A, Casellato C, D’Angelo E and Pedrocchi A (2019) Complex Electroresponsive Dynamics in Olivocerebellar Neurons Represented With Extended-Generalized Leaky Integrate and Fire Models. *Front. Comput. Neurosci.* 13:35. doi: 10.3389/fncom.2019.00035

Received: 15 April 2019; Accepted: 20 May 2019;

Published: 06 June 2019.

Edited by:

Mario Senden, Maastricht University, NetherlandsReviewed by:

Christian Hansel, The University of Chicago, United StatesPaolo Bazzigaluppi, Toronto Western Hospital, Canada

Copyright © 2019 Geminiani, Casellato, D’Angelo and Pedrocchi. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Alice Geminiani, alice.geminiani@polimi.it