Complex Dynamics in Simplified Neuronal Models: Reproducing Golgi Cell Electroresponsiveness

Brain neurons exhibit complex electroresponsive properties – including intrinsic subthreshold oscillations and pacemaking, resonance and phase-reset – which are thought to play a critical role in controlling neural network dynamics. Although these properties emerge from detailed representations of molecular-level mechanisms in “realistic” models, they cannot usually be generated by simplified neuronal models (although these may show spike-frequency adaptation and bursting). We report here that this whole set of properties can be generated by the extended generalized leaky integrate-and-fire (E-GLIF) neuron model. E-GLIF derives from the GLIF model family and is therefore mono-compartmental, keeps the limited computational load typical of a linear low-dimensional system, admits analytical solutions and can be tuned through gradient-descent algorithms. Importantly, E-GLIF is designed to maintain a correspondence between model parameters and neuronal membrane mechanisms through a minimum set of equations. In order to test its potential, E-GLIF was used to model a specific neuron showing rich and complex electroresponsiveness, the cerebellar Golgi cell, and was validated against experimental electrophysiological data recorded from Golgi cells in acute cerebellar slices. During simulations, E-GLIF was activated by stimulus patterns, including current steps and synaptic inputs, identical to those used for the experiments. The results demonstrate that E-GLIF can reproduce the whole set of complex neuronal dynamics typical of these neurons – including intensity-frequency curves, spike-frequency adaptation, post-inhibitory rebound bursting, spontaneous subthreshold oscillations, resonance, and phase-reset – providing a new effective tool to investigate brain dynamics in large-scale simulations.


INTRODUCTION
The causal relationship between components of the nervous system at different spatio-temporal scales, from subcellular mechanisms to behavior, still needs to be disclosed and this represents one of the main challenges of modern neuroscience. To this aim, bottom-up modeling is an advanced strategy that allows to propagate low-level cellular phenomena into large-scale brain networks (Markram, 2013;Markram et al., 2015;D'Angelo and Gandini Wheeler-Kingshott, 2017).
Precise biophysical representations can be generated by "realistic" neuron and network models, but these need then to be simplified to achieve computational efficiency (Gerstner and Naud, 2009;D'Angelo et al., 2016a). Simplified neuron models are fundamental for studying the emergent properties of neural circuits in large-scale simulations and for summarizing in a principled way the electrophysiological intrinsic neural properties that drive network dynamics and high-level behaviors (Gerstner et al., 2014). The specific electroresponsive properties of single neurons are crucial for efficient signal processing, e.g., contributing to noise filtering, signal coding, and synaptic plasticity. The expression of detailed neuron dynamics in simplified models would allow to analyze physiological and pathological phenomena of spiking networks during simulations of sensorimotor tasks in closed-loop control systems, allowing the inference of causal relationships across scales. A critical issue is therefore to obtain simplified neuronal models, that should be at the same time biologically meaningful and computationally efficient.
The first effort toward dimensionality reduction has been made through single-compartment neurons based on the Hodgkin-Huxley (HH) model (Golomb et al., 1993;Guckenheimer et al., 1993;Hill et al., 2001;Doloc-Mihu and Calabrese, 2011). While being able to simulate realistic membrane potential shape and complex patterns like tonic bursting and network oscillations, these models still include more than 5 differential equations, resulting in a significant computational load when simulating large-scale circuits. Other models abstract the biophysical properties of neurons, simplifying the geometry-depended propagation of action potential (Rall, 1962). A compromise between accuracy and efficiency has been reached with simplified mono-compartmental neuron models (or point neuron models), which neglect morphology and loose some functionalities compared to detailed multi-compartmental models but gain computational efficiency. A mono-compartmental neuron model that has been widely used as the basic element of Spiking Neural Networks (SNNs) in different brain areas is the Leaky Integrate and Fire (LIF) (Lennon et al., 2014). LIFs represent neurons as first-order capacitive circuits and embed a threshold-based reset mechanism to reproduce spiking activity (Burkitt, 2006). LIFs are able to generate simple subthreshold dynamics and spike patterns but, in their original formulation, cannot reproduce smooth spike initiation zone, firing adaptation and bursting properties. Non-linear adaptive LIFs have been developed to enhance electrophysiological realism. In the Izhikevich model of cortical neurons, the dynamics of membrane potential, V m , depends on both V 2 m and a membrane recovery variable (Izhikevich, 2003). By introducing an exponential term in the differential equation of membrane potential and an adaptive current with slow dynamics, the action potential shape was well-fitted without the need of a threshold-reset mechanism (Brette and Gerstner, 2005). However, the nonlinearity entailed more difficulties in optimizing model parameters and in computational efficiency. Therefore, recently, new linear adaptive point models have been developed (Generalized LIF, GLIF), with spike-triggered currents and moving threshold as the source of adaptation (Mihalaş and Niebur, 2009) and with stochastic processes in firing emission (Pozzorini et al., 2015;Rössert et al., 2016). The possibility to use a linear and analytically solvable neuron model is fundamental when simulating large-scale SNNs, since computational efficiency can be enhanced without severe loss in spike time accuracy and realism (Hanuschkin et al., 2010). However, unless an oscillating input current is applied (Brunel et al., 2003), these GLIF with first-order dynamics can hardly generate phenomena like intrinsic subthreshold oscillations, resonance and phase-reset, which are critical for large-scale network entrainment and communication (Buzsáki, 2004;Buzsáki and Draghun, 2004).
We propose here an extended GLIF (E-GLIF) model, which achieves a sound compromise between model complexity, biological plausibility and computational efficiency (Figure 1) (Geminiani et al., 2018b). The model has only three state variables, the membrane potential and two currents, which can be associated to main biophysical subcellular mechanisms. Thanks to its mathematical structure, which is similar to GLIF and analytically solvable, E-GLIF can be optimized by traditional optimization algorithms (Pozzorini et al., 2015;Teeter et al., 2018) avoiding metaheuristic methods, like Genetic Algorithms, used for multi-compartment realistic neurons with highdimensional parameter search space (Masoli et al., 2015) or nonlinear LIF models (Venkadesh et al., 2018). E-GLIF can reproduce a rich variety of electroresponsive properties with a single set of optimal parameters: autorhythmicity, depolarization-induced excitation and post-inhibitory rebound bursting, specific inputoutput (f-I stim ) relationships, spike-frequency adaptation, phasereset, sub-threshold oscillations and resonance. A comprehensive example of this entire set of excitable properties was given by the E-GLIF of a cerebellar Golgi cell (GoC), whose electrophysiological properties have been carefully investigated and modeled previously using realistic multi-compartmental approach (Forti et al., 2006;Solinas et al., 2007a,b). And since the GoC expresses among the most common and varied electroresponsive properties, the E-GLIF model should be easily applied to various central neurons that can be represented with mono-compartmental models and promote the investigation of complex brain dynamics in large-scale simulations (Jordan et al., 2018).

METHODS
In this work, taking the move from previous GLIF neurons (Mihalaş and Niebur, 2009;Hertäg et al., 2012;Pozzorini et al., 2015), we have developed and tested the E-GLIF neuron. The model was implemented in the Neural Simulation Tool (NEST) (Diesmann and Gewaltig, 2002), using NESTML (Plotnikov et al., 2016) and the C++ core of NEST. Experimental recordings were performed from cerebellar GoCs using patch-clamp recording techniques for validation.

The Model
E-GLIF couples time-dependent and event-driven algorithmic components to generate a rich set of electrophysiological behaviors, while keeping the advantages of LIF neuron models FIGURE 1 | From complex to simplified neuron model. Neurons exhibit complex properties that affect synaptic plasticity, circuit dynamics and behavior. The point neuron described here shows that it is possible to capture this complexity through a simplified model. The top left panel shows a Golgi cell of mouse cerebellum stained with red tomato and reconstructed using a fluorescence confocal microscope (courtesy of Prof. Javier De Felipe, Department of Neuroanatomy and Cell Biology, Instituto Cajal (CSIC), Madrid, Spain). The Golgi cell shows a complex set of dendritic and axonal ramifications. In whole-cell patch-clamp electrophysiological recordings from the soma, the Golgi cell behaves as a pacemaker generating adapting spike trains in response to step current injection (bottom left) and other more complex properties that are explained in this paper. As shown in the picture, membrane potential firing and subthreshold properties recorded in vitro (Left) can be observed in simulated spiking patterns and subthreshold dynamics obtained through the E-GLIF point neuron (Right).
in terms of simplicity and analytical solvability. The E-GLIF neuron includes 3 linear Ordinary Differential Equations (ODEs) describing the time evolution of membrane potential (V m ) and of two intrinsic currents (I adap and I dep ). Each of these three state variables is modified by an update rule at spike events, which are generated according to a probabilistic threshold crossing controlled by an escape noise.
The model is defined as follows: • State variables evolution (second-order 3-dimenstional ODE system): 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.
• Spike generation: if the neuron is not in the refractory interval t ref , a spike is generated stochastically at t spk , based on a point process that depends on the escape rate function λ(t) as explained more in detail in Appendix I (Gerstner and Kistler, 2002;Jolivet et al., 2006): where: V th = threshold potential; λ 0 , τ V = escape rate parameters. Spike times do not correspond strictly to when the potential reaches the threshold: they occur with higher probability if the membrane potential is near the threshold, depending on the parameters τ V and λ 0 that define the minimum distance from threshold corresponding to the maximum probability of having a spike (Gerstner et al., 2014).
• Update: following a spike, the state variables are modified according to the rules: where: t + spk = time instant immediately following the spike time t spk V r = reset potential; A 2 , A 1 = model currents update constants.
The parameters in the model include those directly related to neurophysiological quantities (C m , τ m , E L , t ref , V th , V r in blue), that are fixed for each specific cell type, and the more abstract ones related to neuron-specific functional mechanisms, that need to be optimized (k adap , k 2 , k 1 , A 2 , A 1 , I e , in red). In addition to the leaky current term, C m τ m (V m − E L ), each one of the membrane currents defined in the model (I e , I adapt , I dep ) accounts for a different mechanism that can be properly parameterized: • I e is an endogenous current modeling the net contribution of depolarizing ionic currents generating autorhythmicity (Mihalaş and Niebur, 2009). • I adap is an adaptive current, usually hyperpolarizing, which is characterized by a small spike-triggered increment (A 2 ) that decays thereafter according to k adap and k 2 . I adap models the activation of potassium channels generating a slow hyperpolarizing current. Since I adap activates slowly while I dep is already decaying, the balance between the two currents generates spike-frequency adaptation and after Hyperpolarization potential. Moreover, by being coupled with V m by k adap , I adap endows the model with the capability of generating post-inhibitory rebound bursting, intrinsic subthreshold oscillations and resonance (Brette and Gerstner, 2005;Hertäg et al., 2012). • I dep is a depolarizing spike-triggered current, which has a larger spike-triggered increment (A 1 ) and faster decay (k 1 ) compared to I adap . I dep mimics the fast (almost instantaneous) activation and deactivation of sodium channels. I dep can generate depolarization-induced excitation and sustain post-inhibitory rebound bursts (Mihalaş and Niebur, 2009).
By computing the analytical solution of the model (see Appendix I), we were able to associate different regions in the parameter space to different system responses (i.e., exponential or oscillatory and stable or unstable). Specifically, as reported in Figure 2A, the k 2 -k adap plane includes an area corresponding to exponential and stable solutions (green in the figure) and an area with oscillatory stable solutions (red in the figure). Within the latter, the red line corresponding to k 2 = 1 τ m defines the condition for not-damped oscillatory solutions, which allow to reproduce self-sustained oscillations of the membrane potential.

Model Parameters
To generate a neuron-specific model, we first considered the parameters that are directly measured as neurophysiological quantities (fixed parameters, highlighted in blue in section The Model): C m , τ m , E L , t ref , V th , V r were fixed to biological values taken from literature or available from animal experiments or databases (Tripathy et al., 2014). For the other neuronspecific functional parameters (tunable parameters, highlighted in red in section The Model), k adap , k 2 , k 1 , A 2 , A 1 , I e , we developed an optimization strategy based on a desired input-output relationship, considering a current step I stim as the input and spike times as the output. Specifically, we supposed to evaluate the neuron response to the inputs listed in Table 1: -I stim = 0 pA: zero current (zero_stim) generating spikes at frequency tonic_freq to evaluate autorhythm; -I stim at three increasing excitatory current steps (exc1 < exc2 < exc3) producing firing with increasing frequency (freq1 < freq2 < freq3) to reproduce the f-I stim relationship and spikefrequency adaption (i.e., steady-state decreased frequency with gain1 > gain2 > gain3); -I stim = inh, an inhibitory input current to evaluate the occurrence of an inhibition-induced silence followed by rebound burst, made of at least 2-spikes.
Silence period during hyperpolarization and return to spiking with at least 2-spike burst (faster than tonic_freq) when hyperpolarization stops.

Post-inhibitory rebound bursting
Cost Function, Constraints, and Algorithm To evaluate different parameter sets by computing the corresponding cost function, we exploited the analytical tractability of the model and we evaluated the model solution V m (t), within the most significant time windows (initial, transitory, and steady state) during each stimulation current step. For each I stim = (i) = zero_stim, exc1, exc2, exc3, the three time windows taken into account were: start to first spike t (i) 1 , t (i) 2 between t (i) 1 and second spike time t (i) 2 , and t (i) ss between two spikes at Steady-State (ss).
At the beginning of each depolarizing phase (i.e., during t (i) 1 ), we supposed the state variables to be initialized as follows: Then, for the following time intervals, the initial conditions were derived from the update rules of the model, supposing that the system had reached the steady-state condition when the adaptive current, I adap , decayed during the interspike interval of an amount equal to A 2 , i.e., its update constant. Analogously, to evaluate rebound bursting, we computed the solutions with I stim = 0 pA, just after a hyperpolarization (I stim = inh < 0 pA), within 2 consecutive time windows: t (inh) lat_reb , from the end of the inhibitory current stimulus to the first rebound spike time t lat_reb , and during the 1 st Inter-Spike Interval (ISI) of rebound burst (t (inh) 1st_reb ). Again, the initial conditions were derived using the update rules.
Starting from all the solutions computed in terms of V m (t), for each input I stim in Table 1, we derived the desired spike times (t (i) 1_des , t (i) 2_des , t (i) ss_des , or t (inh) lat_reb_des and t (inh) 1st_reb_des ), by imposing The error is computed as the difference of blue and green areas, neglecting trend above V th . (C) The short continuous stimulation protocol is represented. This is used for assessing model optimization, applying the same input current values used during optimization, i.e., a zero-current phase, 3 depolarizing steps with input exc1, exc2, exc3 = 200, 400, 600 pA, and a hyperpolarizing step with input inh = −200 pA, all interleaved with zero-current phases. (D) The full multi-phase in vitro stimulation protocol is reported. This has been used for the continuous long simulations (48.93 s) to validate the model and consequently for the input definition in the experimental protocol. After 10 s at zero current (spontaneous activity), 6 increasing current steps are imposed (from 100 to 600 pA, with 100 pA of increase), each lasting 1 s; they are interleaved by 1 s of zero-current. Then, two impulses are provided (4,000 and 6,800 pA, each lasting 0.5 ms). Then, 2 negative current steps are imposed (−100 and −200 pA), each lasting 1 s. Finally, 13 sequences of pulse trains are provided. Each sequence is made up of 5 pulses (each 600 pA lasting 30 ms); the time interval between pulses is constant within each sequence (2, 0.5, 0.28, 0.2, 0.17, 0.13, 0.1, 0.08, 0.06 ms).
a spike event at the time when V m (t) = V th . Therefore, the spike generation during optimization was assumed deterministic.
To take into account the variability in spike generation due to the stochasticity into the model, for each stimulation input during training, we used a distribution of 10 desired firing frequencies with specific mean and Standard Deviation (SD), and thus a distribution instead of a single target spike time.
The parameters k adap , k 2 , A 2 , k 1 , A 1 , I e were optimized through the Sequential Quadratic Programming (SQP) algorithm from the MATLAB R2015b Optimization Toolbox, with normalized parameter values and constraints. SQP optimization aims at the simultaneous minimization of a cost and constraint function, using a gradient-based minimization method. We set the stopping criteria to the iteration when the variations in parameter search, cost or constraint functions were below 10 −3 or the cost function reached a value lower than 0.1. We chose the cost function for gradient minimization, summing all the errors relative to the analyzed properties along the stimulation protocol as in Table 1 (zero input, 3 excitatory currents, first and second spikes of the post-inhibitory rebound burst).
Each error term compared the actual area A act , i.e. the area between the actual V m (t) curve and the V m (t) = V th horizontal line, during each desired ISI, to the target area A target ( Figure 2B): Frontiers in Neuroinformatics | www.frontiersin.org for each input current (i) where: computing only the area from the relative initial time t 0 = 0 to the desired spiking time, before reaching the threshold, since spike-reset mechanisms were enabled in the model simulations.
that can be considered an approximation of the ideal subthreshold membrane potential for a linear LIF neuron not firing at very low spiking frequencies.
considering two aspects of the rebound burst following inhibition, the time of the first spike after hyperpolarization and the distance of the second spike after it.
Based on the mathematical considerations in section The Model and Appendix I about the dynamics of the solutions, the parameter space needed to be constrained to obtain the desired membrane voltage evolution (Figure 2A). Further limits in the parameter space could be included to take into account neuron-specific information from neurophysiology.
In order to evaluate the convergence and stability of the optimization process, for each optimization run, parameters were initialized to random values within their ranges and 5 optimizations were performed with different random initializations. The median of all the resulting parameters was then considered as the final optimal set and used to run complete simulations using the simulator PyNEST (Eppler et al., 2009) for assessing the optimization results and for further validation.

E-GLIF Model of the Cerebellar Golgi Cell
E-GLIF model and optimization were applied to reproduce the complex electroresponsiveness of cerebellar Golgi cells (GoC). GoCs are the main inhibitory neurons in the granular layer of cerebellum and are responsible for reshaping the input signals coming from mossy fibers. In single-cell recordings, GoCs show spontaneous firing around 8 Hz, a nearly-linear input-output relationship (about 0.25 Hz/pA), input-dependent spike-frequency adaptation when depolarization is maintained, rebound bursting after hyperpolarization, phase-resetting, subthreshold self-sustained oscillations and resonance in the theta band (around 3-6 Hz) (Forti et al., 2006;Solinas et al., 2007a,b;D'Angelo et al., 2016a). A multi-compartmental realistic model (Solinas et al., 2007a,b) assumed that dendrites were passive and used them to redistribute the passive electrotonic load while placing all the ionic channels in the soma, suggesting that an appropriate single point model could have been effective as well. In the present E-GLIF model, all electrical properties are collapsed in a point and gating kinetics of ionic channels are substituted by lumped and simplified membrane mechanisms.

Model Construction and Optimization: Physiological Parameters, Cost Function and Constraints
As reported in section Model parameters, the values of the electrophysiological parameters were taken from literature (Forti et al., 2006;Solinas et al., 2007a,b) and databases (Tripathy et al., 2014) describing experimental GoC properties in vitro; their values, used in optimization and simulations, are listed in Table 2.
The optimization of the remaining tunable parameters was achieved setting proper values to the target behavior of the optimization algorithm, derived from literature (Forti et al., 2006;Solinas et al., 2007a,b), as shown in Table 3.
In order to account for the whole-set of GoC electrophysiological properties, the cost function included all the terms reported in cost_function in section Cost function, constraints, and algorithm. In addition, the parameter space was limited to fulfill mathematical and neurophysiological constraints. First, once we had set the neurophysiological parameters to literature values, we constrained k 2 and k adap to obtain not-damped oscillations of the membrane potential at a preferred frequency. We also constrained V m oscillation amplitude during zero and hyperpolarizing input. Further boundaries were applied on A 1 , A 2 , I e, and k 1 to obtain neurophysiological values of currents into the neuron. Additional details on the parameters constraints are described in Appendix II.
The optimal parameter set was chosen as the median of the final values over the 5 optimization runs with different random parameter initialization.

Model Simulations
Following parameter optimization in MATLAB, we derived irregular firing parameters (λ 0 and τ V ), with the aim to obtain a physiological variability of spike events during autorhythmicity. In fact, in the PyNEST model simulations, we activated the firing stochasticity that was disabled during optimization. Then, the GoC simulations in NEST proceeded in two phases.
First, we implemented a multi-step protocol with the same input currents used for optimization, taken from literature (see Table 3 and section Model Construction and Optimization:   (Forti et al., 2006;Solinas et al., 2007a,b).

Input current step (I stim ) Expected output (mean ± SD)
Desired spike times Corresponding properties Post-inhibitory rebound bursting rebound_freq > 2·tonic_freq t 1st_reb_des < 62.5 ms (uniform distribution) Input: current step values (zero_stim, exc1, exc2, exc3, inh). Output: the distribution of desired firing rates during the autorhythmic phase (tonic_freq), at the beginning and end of depolarizing phases (freq1, freq2, freq3), and during the rebound burst (latency and rebound_freq compared to tonic_freq). For each output, the desired spike intervals in the significant time windows ( t 1st_reb_des ) are computed for each input current (superscript in brackets).
Physiological Parameters, Cost Function and Constraints). This phase was fundamental to assess the effectiveness of parameter tuning in continuous simulations with escape noise during application of dynamic input patterns, and not only in sample intervals (i.e., t (i) 1 , t (i) 2 , t (i) ss , or t (inh) lat_reb and t (inh) 1st_reb ). In these sequences, after a 10-s zero-current phase used to evaluate firing irregularity, we delivered three 1 s steps of increasing amplitude (exc1, exc2, exc3 = 200, 400, 600 pA) to monitor the f-I stim slope and spike-frequency adaptation, interleaved with 1-s zero-current to let the neuron recover to its spontaneous activity. We then stimulated the neuron with a 1-s inhibitory step current (−200 pA) for evaluating rebound bursting during a subsequent zero-current step ( Figure 2C).
Secondly, we added stimulation patterns required to evaluate the emergence of model features (autorhythmicity, rebound bursting, f-I stim , adaptation) at higher resolution and to check for further emergent properties that were not considered during the optimization process (e.g., resonance, phase-reset), including: • an initial zero-current phase of 10 s to evaluate frequency and irregularity of intrinsic firing; • six depolarizing steps lasting 1 s, with input currents ranging from 100 to 600 pA (increments of 100 pA), interleaved with 1-s zero-current phases, to test intrinsic excitability; • two zero-input phases lasting 2.5 s, where we provided a short pulse excitation (amplitude 4 nA and 6.8 nA, for 0.5 ms), to measure phase-reset mechanism; • two 1-s hyperpolarizing phases with inhibitory current −200 and −100 pA, followed by a 1-s zero-stimulation phase, where evaluating rebound bursting properties; • a sequence of 5 steps, each of amplitude 600 pA, lasting 30 ms, at increasing frequencies: 0.5-2-3.5-5-6.3-7.7-10-12-15 Hz, to evaluate resonance.
The resulting total duration of the whole protocol was 48.93 s ( Figure 2D). Each simulation with the same optimized set of Golgi neuron parameters was run 10 times with different seeds of the random number generator used to produce the escape noise, and thus spike stochasticity. Finally, we assessed the capability of the model to intrinsically generate self-sustained V m oscillations, running a 1-s simulation with I stim = 0 pA, where V th was increased to −5 mV to avoid the spiking mechanisms, which partially hide the underlying oscillations.

Model Synaptic Activation
Synaptic mechanisms were added to the model in order to allow its connection with different input neural populations and simulate the GoC response to network activity. The synapses were conductance-based, with a spike-triggered change of the synaptic conductance g syn according to an alpha function (Roth and van Rossum, 2013): This change of synapse conductance caused a change in the V m through the input synaptic current: This kind of synaptic model was chosen to maximize the realism in synapses behavior, despite losing the neuron model analytical tractability when connecting to other neurons and increasing the computational load. This solution can be considered acceptable when using a small-/medium-scale SNN. However, the NEST platform flexibility guarantees the possibility to use the E-GLIF with current-based synapses, which are less realistic but also require a lower computational load and can be suitable for large-scale SNN simulations (Cavallari et al., 2014). In simulations of synaptic responsiveness, synaptic strength G syn and delay τ syn were set to different values for excitatory and inhibitory synapses (excitatory synapses: G syn = 40 nS and τ syn = 0.1 ms, inhibitory synapses: G syn = 10 nS and τ syn = 0.1 ms). GoC E-GLIF was connected to two different spike generators, excitatory and inhibitory, respectively. We provided a 50 Hz spike train on the excitatory synapse lasting for 800 ms, while the inhibitory synapse received a short spiking burst, for evaluating the capability of the neuron model to produce a rebound burst after an inhibitory spiking input.

Experimental Data Acquisition and Analysis
The outcome of simulations using the long multi-phase protocol was compared to real data acquired ad hoc from mice GoCs in acute cerebellar slices. These data are not identical to those derived from literature, allowing therefore to test the generalization capability of GoC E-GLIF beyond the specific dataset used for model construction. As for the model, also for the data we evaluated the same electrophysiological properties and used similar inputs to the neurons described in section Model Simulations (see the protocol in Figure 2D).
The experiments have been conducted on 16-to-21-days-old (P0 = day of birth) male and female mice heterozygous for the bacterial artificial chromosome insertion of EGFP under the control of the glycine transporter type 2 gene (Zeilhofer et al., 2005) (GlyT2-GFP mice). All procedures were conducted in accordance with European guidelines for the care and use of laboratory animals (Council Directive 2010/63/EU), and approved by the ethical committee of Italian Ministry of Health (628/2017-PR). The mice were anesthetized with halothane (1 ml in 2 L administered for 1-2 min) and killed by decapitation in order to remove the cerebellum for acute slice preparation according to a well-established technique (Forti et al., 2006;Cesana et al., 2013). The cerebellum was gently removed, and the vermis was isolated and fixed on the stage of a vibroslicer (Leica VT1200S) with cyanoacrylic glue. Acute 220-µm-thick slices were cut in the parasagittal plane in cold Kreb's solution and maintained at 32 • C before being transferred to a 1.5 ml recording chamber mounted on the stage of un upright epifluorescence microscope (Axioskop 2 FS; Carl Zeiss, Oberkochen, Germany) equipped with a 63, 0.9 NA water-immersion objective. Slices were perfused with Kreb's solution and maintained at 32 • C with a Peltier feedback device (TC-324B, Warner Instrument Corp., Hamden, CT). The Kreb's solution contained the following (in mM): 120 NaCl, 2 KCl, 1.2 MgSO 4 , 26 NaHCO 3 , 1.2 KH 2 PO 4 , 2 CaCl 2 , and 11 glucose, and was equilibrated with 95% O 2 and 5% CO 2 , for pH 7.4.
Whole-cell patch-clamp recordings were performed with Multiclamp 700B (-3dB; cutoff frequency 10 kHz), sampled with Digidata 1550 interface, and analyzed off-line with pClamp10 software (Molecular Devices). Patch pipettes were pulled from borosilicate glass capillaries (Sutter Instruments, Novato, CA) and filled with internal solution containing (in mM): potassium gluconate 145, KCl 5, HEPES 10, EGTA 0.2, MgCl 2 4.6, ATP-Na 2 4, GTP-Na 2 0.4, adjusted at pH 7.3 with KOH. Pipettes had a resistance of 3-5 M when immersed in the bath. Signals were low-pass filtered at 10 kHz and acquired at 50 kHz. The stability of whole-cell recordings can be influenced by modification of series resistance (R s ). To ensure that R s remained stable during recordings, passive electrode-cell parameters were monitored throughout the experiments.
In each recording, once in the whole-cell configuration, the current transients elicited by 10 mV hyperpolarizing pulses from the holding potential of −70 mV in voltage-clamp mode showed a biexponential relaxation. The recording properties for 5 cells from 3 mice were measured as follows (values reported as mean ± Standard Error of Mean): membrane capacitance was evaluated from the capacitive charge (37.1 ± 8.5 pF), while the membrane resistance was computed from the steady-state current flowing after termination of the transient (195.9 ± 97.8 G ). The 3 dB cutoff frequency of the electrode-cell system, f VC , was calculated as f VC = (2π • 2τ VC ) −1 , with τ VC = 215.7 ± 7.4 µs, resulting in 0.7 ± 0.02 kHz. These properties were constantly monitored during recordings. Cells showing variation of R s higher than 20% were discarded from analysis.
After switching to current clamp, GoCs were maintained at resting membrane potential by setting the holding current at 0 pA. Intrinsic excitability was investigated by injecting 1-s steps of current (from −200 to 600 pA in a 100-pA increment). Resonance was investigated by applying sequences of 30-ms 600 pA current steps repeated 5 times at different frequencies (range from 2 to 10 Hz) or injecting a single 0.5-ms step of 4 nA. Finally, the cells were maintained at their resting membrane potential for evaluating subthreshold V m oscillations.
Using an automatic spike detection algorithm, spike times were extracted from the experimental recordings, and used for electrophysiological feature extraction (see section Feature Extraction).

Feature Extraction
To evaluate GoC electroresponsive behavior, multiple parameters were computed according to (Solinas et al., 2007a,b), from both simulation and recording spike trains: • the tonic firing rate, f tonic , was the inverse of the mean ISI, during the initial zero-current phase; • the coefficient of variation of inter-spike intervals (CV ISI ) was measured during the 10-s zero-input step to quantify the irregularity of firing; • the firing rate, f, was the inverse of the mean ISI, during the first 2 spikes of each depolarizing phase; • the f-I stim slope was derived from initial responses to the excitatory step currents; • the f ss /f ratio was used to evaluate spike-frequency adaptation, where f ss was computed at the end (last 5 spikes) of the first 1 s interval of depolarizing stimulation; • the parameter phase_param = ISI post−pre ISI tonic was computed to evaluate phase-reset; ISI post−pre was the interval between the 2 spikes preceding and following the impulse current, while ISI tonic was the average ISI during the 2.5-s zero-input intervals of phase-reset testing ( Figure 2D); • latency and frequency of 1 st spike were measured in the rebound burst after hyperpolarization (lat_rebound and rebound_freq, respectively); • the response speed was computed to evaluate resonance, as the inverse of the mean spike latency in each resonance step; then the values from multiple simulation tests and frequencies were fitted through a smoothing spline in order to obtain the resonance curve .
Parameter values for simulations and recordings are reported as mean ± SD. In addition, for quantifying experimental subthreshold oscillations, the power spectral density of 1 s V m traces was computed and, for each recording, the main oscillation frequency was associated to the spectrum peak.

RESULTS
The E-GLIF neuron is a linear mono-compartment neuron model using three differential equations for membrane potential (V m ), a depolarizing current (I dep ), and an adaptation current (I adap ). These state variables are updated at each spike event, which occurs according to a probabilistic threshold crossing controlled by escape noise. E-GLIF contains 6 fixed parameters (C m , τ m , E L , t ref , V th , V r ) and 6 tunable parameters (k adap , k 2 , A 2 , k 1 , A 1 , I e ) and is designed to obtain a flexible representation of complex firing dynamics. The equation system is analytically solvable. The interdependency across state variables makes the model able to show oscillatory behaviors. The update mechanisms occurring at spike events, that are specific for each state variable, generate adaptive behaviors on multiple time scales.
In the present development, an optimization strategy using a gradient-based minimization method allowed full control over the search of optimal parameters. This operated with explicit solutions related to specific time frames making unnecessary to compute the state variables at each time instant throughout a simulation. Desired spike times were imposed depending on input stimulation.
To challenge E-GLIF in a complex and meaningful case, the cerebellar GoC was chosen as the target neuron to be modeled. Once extracting the optimal parameters, a complete simulation of the neuron could be carried out by applying a multi-phase stimulation protocol able to reveal the emergence of the various spiking response patterns. An identical protocol was used to stimulate real GoCs in acute cerebellar slices, so as to obtain experimental data for robust model validation.

Optimization
In order to reproduce non-damped membrane potential oscillations, k 2 was set to 1 τ m = 0.02 ms −1 (section Model Construction and Optimization: Physiological Parameters, Cost Function and Constraints). For the remaining parameters, multiple optimization runs with random initialization converged toward similar values for all the 5 runs within 200 iterations, without achieving values on the boundary of the search space ( Figures 3A-E). Optimization, exploring the 5-D space of tunable parameters, stably met minimization of the cost function and compliance with all the mutual parameter constraints: all the optimization runs ended with a cost function value lower than 0.5 and all the constraints verified ( Figure 3F). Considering the median of the final parameter sets, the resulting parameters were: k adap , A 2 , k 1 , A 1 , I e = [0.22 MH −1 , 178.01 pA, 0.03 ms-1, 259.99 pA, 16.21 pA]. These values were used for all simulations, thus making the different neuron responses (e.g., autorhythm, adaptation, or rebound firing) depending only on the input stimulus rather than on property-specific parameter sets.

Model Responsiveness During Current
Step Protocols The optimized model was tested by 10 simulations in PyNEST (Eppler et al., 2009), letting V m and the two synthetic currents to evolve in time and to be updated at each spike, thus generating corresponding spike patterns, following the stochastic rule (see section The Model).
Using the parameters resulting from the optimization process, the model was able to generate a linear relationship between I stim and response frequency f, with a constant slope of 0.2 Hz/pA ( Figure 4A). The autorhythm was at frequency f tonic = 12.8 ± 0.02 Hz (Figures 4A,B). The escape noise process caused a slight ISI variability (CV ISI = 3.4 ± 1.4% during the zero-current phase) with the parameters λ 0 and τ V set to 1 ms −1 and 0.4 mV, respectively.
The model response at increasing I stim lasting 1 s was (mean ± SD): I stim = 200 pA, firing rate 49 ± 6 Hz at the beginning of the current step and 36 ± 0.2 Hz at the end; I stim = 400 pA, firing rate 90 ± 10 Hz at the beginning of the current step and 53 ± 0.2 Hz at the end; I stim = 600 pA, firing rate 134 ± 8 Hz at the beginning of the current step and 68 ± 0.2 Hz at the end ( Figure 4A). Figure 4C, the higher initial frequencies reflected the depolarization-induced excitation driven by I dep , while the decrease of the firing rate along the stimulation step was mainly caused by the slightly slower I adap increase and corresponds to spike-frequency adaptation, which becomes more pronounced at higher stimuli ( Figure 4A) (Solinas et al., 2007a). At the end of the stimulus current step, when I stim goes back to zero, the balance among model currents result in a hyperpolarization of the neuron followed by a return to autorhythm.

As shown in
After a hyperpolarizing current step (-200 pA), the model showed a post-inhibitory rebound doublet (i.e., a couple of spikes at a higher frequency followed by a short quiescent period) before recovering to spontaneous firing rate, with latency 30 ± 13 ms and frequency 47 ± 5 Hz ( Figure 4D) (Izhikevich, 2006). This effect reflected the different dynamics of the two currents, I adap and I dep , affecting V m : the current I adap , being coupled with V m , reached negative values during inhibitory stimulation and thus contributed to depolarize the neuron when the stimulation stopped (section The Model), affecting the latency of rebound burst. After the first spike, the fast current I dep sustained burst persistence; after the first 2 spikes I adap attained a steady state balance with I dep , bringing back the activity to the autorhythm (Figure 4D). It should be noted that during the FIGURE 3 | Parameter optimization. (A-E) Tunable parameters evolution over 5 optimization runs. Normalized parameter values (k adap , A 2 , k 1 , A 1 , I e ) are represented along iterations. In each run (different colors), the parameter starts from random values within permitted ranges (shown normalized on y-axis). In each panel, the red point represents the optimal final parameter value (median across optimization runs). (F) Cost and constraints functions over 5 optimization runs. The cost and the constraint functions are reported for the iterations in 5 optimization runs (different colors). Especially in the starting iterations, the optimization algorithm tries to minimize the cost function while respecting the constraints. Parameter search can cause evident changes in the cost function value, especially in the first explorative iterations.
hyperpolarizing phase, the subthreshold V m is not constant but exhibits oscillations due to the second-order dynamics of the system. Nevertheless, these oscillations do not generate spikes (see Appendix II) and the realistic spiking behavior of the neuron is preserved.

Model Validation Against Experimental Data
Experimental data from 5 GoCs showed physiological inter and intra variability in recorded voltage traces under different stimulation protocols. They exhibited autorhythm at a rate of 11.5 ± 8 Hz, increasing depolarization-induced excitation and spike-frequency adaptation in response to positive current steps, rebound doublet bursts after negative current steps (Figures 5A-C). These features were tested in the GoC model in a validation test using more numerous and different current inputs (see Figure 2D). The autorhythm was stably reproduced (Figures 5D, 6A); the linearity between I stim and response f was maintained over multiple I stim levels (Figures 5E, 6A), confirming the results of the first simulations with a shorter stimulation protocol (Figure 4). With increasing input current steps, the balance of model currents at the end of stimulation resulted in a pronounced hyperpolarization before returning to autorhythm. This effect, that was not observed in physiological recordings, was due to the high value of I adap at the end of current steps required to achieve spike-frequency adaptation.
The rebound burst systematically occurred as a doublet after a hyperpolarization and its internal speed and latency increased, with higher absolute values of the preceding negative current steps, consistent with experimental results (Figure 5F, Table 4).
A linear fitting on the experimental data resulted in the same f-I stim slope (0.2 Hz/pA) as when fitting simulation data ( Figure 6A). Experimental recordings and model behaviors evidently differed in the steady-state response rate (f ss ) that may be related to experimental mechanisms not modeled in simulations. Indeed, the experimental f ss values were lower than in the model (Figure 6A), thereby reducing adaptation in the model compared to experiments ( Figure 6B).
An important phenomenon evident in GoCs (as well as in some other brain neurons) is phase-reset, which allows desynchronization within sub-circuits triggered by strong impulses (Figure 7A; Buzsáki, 2004;Buzsáki and Draghun, 2004). The GoC E-GLIF was able to reproduce this feature (Figure 7B), thanks to the coupling between I adap and V m that caused a rapid increase of I adap when V m value raised following a huge external pulse; this blocked spike generation for the same time interval independent from the phase of autorhythm before the pulse, resetting the cells phase of the autorhythm. Impulse amplitude was an element affecting the after hyperpolarization duration following the pulse-triggered spike, slightly increasing the phase-reset parameter as the pulse charge increased ( Table 4). The E-GLIF simulations matched the experimental results with the impulse at 4,000 pA and predicted the slightly increased phase-reset in case of higher pulse current (6,800 pA), confirming literature results (Solinas et al., 2007b).
Finally, another property evident in GoCs and often observed in central nervous system neurons is the presence of endogenous subthreshold oscillations, which provide a fundamental mechanism for efficient network intercommunication and plasticity. Oscillations are correlated to resonance, and both have been shown to depend on intrinsic membrane properties (Hutcheon and Yarom, 2000). The cerebellar GoCs exhibit intrinsic theta-frequency subthreshold oscillations and resonance (Solinas et al., 2007b) that are thought to be instrumental to propagate similar properties throughout the entire cerebellum granular layer . GoCs work as a band-pass filters by amplifying the input in the theta band. The present experimental data showed V m subthreshold oscillations at 4.2 ± 1.4 Hz ( Figure 8A): the normalized power spectral density from 9 experimental V m traces had a peak between 2 and 6 Hz (theta band). As a result, during 5 recordings with the periodic stimulation protocol, the average response speed curve (evaluating resonance) exhibited a maximum at 3.5 Hz ( Figure 8B) and all the curves in each individual recording also had the highest value at 3.5 Hz. Consistently, the simulated GoC model behaved as a band-pass filter with a prominent peak at 3.5 Hz (Figure 8D), while exhibiting pure V m self-sustained oscillations at 5.5 Hz (Figure 8C) (section Model Construction and Optimization: Physiological Parameters, Cost Function and Constraints). In order to verify the low-and high-pass properties of the f stim -response rate curves, the model was simulated with additional stimulation frequencies and generated points falling within the predicted band-pass filter region ( Figure 8D). It should be noted that the slight difference among the frequencies of autorhythm, subthreshold oscillations, and resonance (Hutcheon and Yarom, 2000) observed in the experimental data [for a similar effect see (Solinas et al., 2007a,b)] was also observed in model simulations.

Model Responsiveness to Synaptic Stimulation
As shown in Figure 9B, synaptic stimulation of E-GLIF caused an alpha-shaped change in synaptic conductance. The excitatory spike train increased irregularity in neuron firing (CV ISI = 39%), consistent with the irregular spiking of GoCs in vivo (Cerminara and Rawson, 2004;D'Angelo, 2009). The inhibitory spike train generated a rebound burst in E-GLIF after the end of the input stimulus, due to the intrinsic changes of the model currents ( Figure 9A). This result confirmed the ability of GoC E-GLIF to reproduce the rich variety of electrophysiological properties of GoCs.

DISCUSSION
This work reports the development, optimization and testing of the extended generalized leaky integrate-and-fire neuron model, E-GLIF. E-GLIF is a simplified point-neuron based on a system of 3 linear ordinary differential equations that is able to represent multiple complex electrophysiological mechanisms at different levels of abstraction. Like GLIF, E-GLIF maintains analytical tractability, allows to define different solution regimes and to optimize model parameters through minimization methods. The main improvement provided by E-GLIF is to generate a richer set of neuronal dynamics beyond depolarizationinduced excitation and adaptation, which includes also rebound bursting, phase-reset, intrinsic sub-threshold oscillations, and resonance. Thus, E-GLIF covers almost the whole set of neuronal discharge properties relevant for microcircuit functioning and network entrainment. Moreover, E-GLIF is designed to maintain a traceable correspondence between lumped parameters and ionic conductances of the neuronal membrane. In this way, E-GLIF allows to combine GLIF computational efficiency with the ability of reproducing the salient dynamic properties of neuronal discharge, while keeping insight into the underlying cellular mechanisms. E-GLIF appears therefore suitable to bridge the gap between biophysically detailed realistic models and computationally efficient simplified models, and could be used to investigate the impact of neuronal dynamics in large-scale networks. As a prove of its validity, E-GLIF was shown to reproduce the main spiking discharge properties of the cerebellar GoC, a prototype of a neuron with complex electroresponsiveness.

Model Parameterization and Optimization Strategy
By exploiting its analytical solution, the E-GLIF GoC model was optimized using gradient descent minimization methods. This allowed to fast-tune a unique set of parameters generating the appropriate spiking responses to various input patterns. The cost function was designed to evaluate sub-threshold membrane potential dynamics using the integral of membrane potential all over the inter-spike interval: in this way, the cost function was still differentiable, while taking into account the full history of the signal preceding the spike. Supra-threshold dynamics may not be relevant since E-GLIF used a spike-update-reset mechanism.
In addition, the general cost function could be customized by disabling the terms corresponding to the properties that are not exhibited by the specific neuron model. For example, for a neuron without autorhythm, the error zero_stim term would not be included in the cost function. For the remaining terms, the desired output parameters listed in Table 1 were derived from electrophysiological studies reported in literature and generally applicable to many different neuron types. It should also be noted that, although stochasticity in spike generation was not expressed in the explicit model solution used for optimization, it was accounted for by optimizing the neuron on a distribution of desired spike times. Compared to other optimization strategies, which are based on semi-automatic fitting on a training set of neural traces or spiking patterns (Pozzorini et al., 2015;Rössert et al., 2016) or realistic model mapping (Marasco et al., 2012), E-GLIF tuning was based on features related to neurophysiological activity (primarily membrane capacitance and time constant, resting potential, refractory period, spike threshold) and on inputoutput relationships. Then, the model was validated against real experimental data. This feature-based approach guarantees high generalization capability and would be particularly suitable when the data needed for SNN models reconstructions have been obtained in different laboratories and therefore do not share the uniform database required for automatic fitting and validation.

Comparison to Other Simplified Models
E-GLIF integrates elements taken from the theory of point neuron models in order to reach a good compromise between computational efficiency, number of tunable parameters and biological plausibility (Izhikevich, 2004;Pozzorini et al., 2015;Teeter et al., 2018).
First, E-GLIF reproduces sub-threshold responses and spike timing rather than the shape of action potentials. All plots show a peak in theta band (between 2 and 6 Hz). (B) Response speed values for 5 experimental recordings with the resonance protocol (using f stim from 2 to 10 Hz), and fitted resonance curve showing a maximum at 3.5 Hz, despite data variability. (C) Model membrane potential self-sustained sinusoidal subthreshold oscillations in theta band, when the spiking mechanism is blocked. To avoid firing, we increased the V th value from −55 mV (dotted green line), to −5 mV. (D) Response speed during 10 simulations with the periodic input pattern fitted with a smoothing spline, showing resonance in theta band (higher peak at 3.5 Hz). Simulation results confirm the low-and high-pass effect in the frequency ranges not acquired during experiments (gray areas in panel (B)).
Nonetheless, some spike properties, such as spike width and after hyperpolarization, are accounted for by the reset rules, which map the state after the spike to that before it (Teeter et al., 2018). It should be noted that monocompartment and simplified HH models or non-linear point neuron models would be useful to reproduce supra-threshold dynamics or accurate action potential and bursting initiation (Fitzhugh, 1961;Hindmarsh and Rose, 1984;Guckenheimer et al., 1993;Brette and Gerstner, 2005); however, it is not needed in E-GLIF, which is conceived for embedding into SNN and therefore privileges spike timing and population spiking patterns.
Secondly, the E-GLIF analytical tractability can be exploited to implement versions designed for time-driven simulations, where the exact spike time can be computed through iterative methods like bisection (Hanuschkin et al., 2010). This would be useful to decrease simulation time resolution without losing spike timing accuracy in large-scale SNNs.
Thirdly, E-GLIF exploits an adaptive current coupled with membrane potential to model spiking responses driven by adaptation mechanisms, like spike-frequency adaptation and post-inhibitory rebound bursting, rather than a sliding threshold depending on the actual membrane voltage (Mihalaş and Niebur, 2009;Pozzorini et al., 2015). This E-GLIF feature, in addition to decreasing the number of ODEs in the model, correlates effectively to the membrane mechanisms of firing, as exemplified in the case of cerebellar GoCs (cf. Solinas et al., 2007a,b). In GoCs, both the f-I stim curve shift during adaptation and post-inhibitory rebound bursting are driven by the adapting current I adap , in agreement with ionic mechanisms revealed by data-driven realistic modeling (Solinas et al., 2007a) and with mechanisms adopted in other simplified models (Benda et al., 2010;Naud and Gerstner, 2012). Moreover, the adaptive current coupled with membrane potential allowed to simulate intrinsic self-sustained subthreshold oscillations at a preferred frequency and thus to generate intrinsic resonance (Richardson, 2003). Importantly, this is not possible either in the traditional GLIF neurons (Pozzorini et al., 2015) with firstorder dynamics or in other non-linear adaptive point neurons (Touboul, 2012). The spike-triggered updates of I adap and I dep can be associated to mechanisms of ion-channel activation/deactivation, in particular concerning K + (slow) and Na + (fast) currents (Mihalaş and Niebur, 2009), as further considered below.
Finally, for testing the neuron with synaptic inputs, conductance-based synaptic receptors were embedded into the model. It should be noted that current-based synapses may FIGURE 9 | Responses to synaptic inputs. Simulation of the neuron model response to synaptic inputs on two synaptic receptor mechanism (R exc and R inh ). The synaptic models are conductance-based, with an alpha-shaped conductance change. As shown in (A), the excitatory spike train (bottom panel-red) increases irregularity in neuron firing (top panel), while the inhibitory input burst (bottom panel-blue) produces the rebound burst. Synaptic conductances are modified according to an alpha function, with delay 0.1 ms and maximum change depending on the weight of the synapse (middle panel). A zoom on the conductance change is shown in (B). be used to further improve computational efficiency at the expense of biophysical realism (Cavallari et al., 2014).

Limitations of the E-GLIF Model
While appropriately reproducing the physiological spiking patterns in response to various input stimuli, the E-GLIF neuron presents some limitations. As an Integrate and Fire model, the action potential shape is not reproduced but is approximated to a spike event followed by a refractory period that takes into account the duration of the action potential. In addition, thanks to second-order dynamics, the subthreshold E-GLIF V m intrinsically oscillates at a preferred frequency, also in response to a step current. This allows obtaining spontaneous subthreshold oscillations (Figure 8), but, unavoidably, causes non-physiological oscillations in membrane potential with constant hyperpolarizing currents (Figures 4, 5). In addition, oscillations are highly sensitive to k adap and k 2 : to obtain notdamped oscillations, k 2 needs to be fixed to 1 τ m (Figure 2A) and consequently the frequency of oscillations depends only on k adap (see Appendix II). This is far from the cellular behavior where a balance across multiple mechanisms is responsible for generating the neuron response (Solinas et al., 2007b). These limitations can be considered minor though, since this simple second-order model is able to produce the neurophysiological output activity pattern.
We chose to optimize the neuron model with a multiobjective approach, since we aimed at generating different spiking patterns with a single set of parameters based only on the input stimulus. As a consequence, some quantitative properties do not match exactly the experimental values. For example, the amplitude of subthreshold oscillations is higher in the model than in recordings, but this mismatch does not compromise the realistic spiking behavior. Also, the high-amplitude After Hyperpolarization potential at the end of depolarizing steps does not overlap the experimental behavior. The reason is the high steady-state value of I adap : this is necessary for spikefrequency adaptation, but it unavoidably strongly hyperpolarizes the neuron when the external excitatory current stops. This hyperpolarization causes also the lower values of resonance parameters in simulations (Figure 8), because it slows down the membrane potential dynamics during repetitive current step stimulation. However, this does not affect the final resonance frequency value. Therefore, we preferred to find the best combination of parameters to generate multiple spiking patterns that could be crucial for SNN simulations, instead of fitting a specific single electroresponsive property.

Correspondence of Model Parameters to Subcellular Mechanisms
Given the specific design of E-GLIF (see above), it is possible to trace the relationship between E-GLIF parameters and the ionic and membrane mechanisms generating the spiking response in real neurons and realistic models. This comparison is wellexemplified by considering the ionic mechanisms of cerebellar GoCs, which have been previously determined in great detail (Solinas et al., 2007a,b). Table 5 shows a list of membrane currents generating various electrophysiological properties in the cerebellar GoC (D' Angelo et al., 2013) and their remapping onto the lumped parameters of E-GLIF, at a different level of abstraction.
While, in order to simultaneously generate 8 different electrophysiological properties (action potential, autorhythmicity, depolarization-induced burst, post-inhibitory rebound burst, spike-frequency adaptation, phase-reset, subthreshold oscillations, resonance), the realistic GoC model engages 10 ionic currents (I h , I Na−t , I Na−p , I Na−r , I Ca−HVA , I Ca−LVA , I K−V , I K−A , I K−AHP , I K−slow ), E-GLIF has just 3 currents (I e , I adap , I dep ) controlling a digital spike generator.
Actually, there are correlations between ionic mechanisms (e.g., I Ca−HVA and I K−AHP are coupled), the same ionic mechanisms can influence multiple electrophysiological properties (e.g., the I Ca−HVA /I K−AHP balance influences both autorhythmicity, adaptation, and phase-reset), and some electroresponsive properties are at least partly bound one to each other (e.g. subthreshold oscillations and resonance) reducing the dimensionality of the real neuron parameter space (Solinas et al., 2007a,b). De facto, with its 3 currents and spike-reset mechanisms, E-GLIF can effectively abstract the high-dimensional response pattern of the GoC, supporting the concept that appropriate models can provide a mathematical interpretation of complex neuronal properties (Gerstner and Naud, 2009). It should be noted that the association between model elements and firing responses is more difficult in other simplified non-linear models (e.g., Izhikevich, 2003), which appear to be less interpretable in mechanistic terms. Nevertheless, those models can more accurately reproduce the membrane potential shape and phenomena like tonic bursting (Izhikevich and Hoppensteadt, 2004). In addition, they associate regions of the parameter space to different spiking regimes, instead of generating different firing responses based only on the input stimulus: this makes it difficult to exploit them in large-scale cerebellar simulations where neurons are supposed to generate multiple spiking patterns based on the received external stimuli.

CONCLUSIONS
This work shows that it is possible to represent complex neuronal firing dynamics through a mono-compartment neuron model by updating the GLIF into E-GLIF at the modest computational expense of 3 ODEs. Yet there is a remarkable efficiency gain with respect to realistic multi-compartmental models. For example, compared to the realistic version of the GoC model (Solinas et al., 2007a,b) with 23 ODEs, there is an 80-fold computational time reduction to simulate the same stimulation protocol with the E-GLIF. Together with the computational efficiency, the E-GLIF was able to reproduce multiple electrophysiological features with a single set of model parameters, moving forward the traditional approach of neuron modeling and resulting in a higher biological plausibility (Izhikevich, 2004). Specific advantages of E-GLIF are the second-order dynamics and the linearity: the model admits an oscillatory response to constant inputs and an analytical solution that allows to extend the theoretical analysis of complex firing dynamics. Moreover, E-GLIF keeps a correspondence between lumped model parameters and electrophysiological mechanisms, limiting black-box fitting and supporting the interpretation of neuronal physiological properties and their changes by neuromodulation and plasticity. There is a large category of neurons that could be represented as point processes and could indeed be modeled with E-GLIF. For example, good candidates are the granule cells and the stellate and basket cells in the cerebellum as well as various interneurons in the neocortex and hippocampus. Nonetheless, when dendritic or axonal computations become critical, E-GLIF modifications would be needed, e.g., by connecting multiple E-GLIF style compartments one to each other and adopting dendritic compression and synaptic efficacy scaling procedures (Marasco et al., 2012). This could be the case of Purkinje cells (Masoli et al., 2015;Masoli and D'Angelo, 2017) in the cerebellum and of pyramidal cells in the neocortex and hippocampus (Migliore et al., 2008). Together with the other GLIF neurons, E-GLIF could contribute to generate a standardized database of computationally efficient models capable of generating a rich repertoire of non-linear firing patterns applicable to diverse brain regions and scientific issues by the community (Pozzorini et al., 2015;Tiesinga et al., 2015;Rössert et al., 2016;Teeter et al., 2018). Along this line, the implementation of E-GLIF on the NEST platform (Diesmann and Gewaltig, 2002) is expected to bring salient aspects of neuronal firing dynamics into large-scale network simulations, where different point neuron models can be combined based on the complexity of the neural population to be represented (Destexhe et al., 1996;Geminiani et al., 2018a). The resulting SNN models could help understanding how brain networks generate their computations. Indeed, E-GLIF has been designed to investigate cerebellar network dynamics during closed-loop behavioral testing in neurorobots (Casellato et al., 2014(Casellato et al., , 2015Antonietti et al., 2016;D'Angelo et al., 2016b). In conclusion, E-GLIF could effectively bridge bottom-up and top-down modeling approaches (Herz et al., 2006) paving the way to the establishment of a set of simplified yet biologically meaningful spiking neuron models as the fundamental elements of multi-scale brain modeling.

AUTHOR CONTRIBUTIONS
AG and CC elaborated the mathematical model and optimization, carried out the simulations analyzing them together with the experimental dataset, and substantially contributed to text writing. FL and FP carried out experimental recordings and analysis. AP coordinated mathematical modeling and simulations and substantially contributed to the writing of the final manuscript. ED coordinated the whole work and wrote the final version of the manuscript.