Edited by: Gaute T. Einevoll, Norwegian University of Life Sciences, Norway
Reviewed by: Astrid A. Prinz, Emory University, United States; Pablo Varona, Universidad Autónoma de Madrid, Spain
†Colast authorship
This is an openaccess 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.
Brain neurons exhibit complex electroresponsive properties – including intrinsic subthreshold oscillations and pacemaking, resonance and phasereset – which are thought to play a critical role in controlling neural network dynamics. Although these properties emerge from detailed representations of molecularlevel mechanisms in “realistic” models, they cannot usually be generated by simplified neuronal models (although these may show spikefrequency adaptation and bursting). We report here that this whole set of properties can be generated by the
The causal relationship between components of the nervous system at different spatiotemporal 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, bottomup modeling is an advanced strategy that allows to propagate lowlevel cellular phenomena into largescale brain networks (Markram,
The first effort toward dimensionality reduction has been made through singlecompartment neurons based on the HodgkinHuxley (HH) model (Golomb et al.,
We propose here an
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 wholecell patchclamp 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 this work, taking the move from previous GLIF neurons (Mihalaş and Niebur,
EGLIF couples timedependent and eventdriven algorithmic components to generate a rich set of electrophysiological behaviors, while keeping the advantages of LIF neuron models in terms of simplicity and analytical solvability. The EGLIF neuron includes 3 linear Ordinary Differential Equations (ODEs) describing the time evolution of membrane potential (
The model is defined as follows:
State variables evolution (secondorder 3dimenstional ODE system):
where:
τ_{m} = membrane time constant;
Spike generation: if the neuron is not in the refractory interval Δ
where:
λ_{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.,
Update: following a spike, the state variables are modified according to the rules:
where:
The parameters in the model include those directly related to neurophysiological quantities (
In addition to the leaky current term,
By computing the analytical solution of the model (see Appendix
Mathematical properties of the model.
To generate a neuronspecific model, we first considered the parameters that are directly measured as neurophysiological quantities (



Inputs and expected outputs for the optimization algorithm with corresponding target electrophysiological property.
Firing at 
Autorhythm  
Firing at 

Firing at 

Firing at 

Silence period during hyperpolarization and return to spiking with at least 2spike burst (faster than 
Postinhibitory rebound bursting 
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
At the beginning of each depolarizing phase (i.e., during Δ
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 steadystate condition when the adaptive current,
Analogously, to evaluate rebound bursting, we computed the solutions with
Starting from all the solutions computed in terms of
The parameters
Each error term compared the actual area
for each input current
where:
computing only the area from the relative initial time t_{0} = 0 to the desired spiking time, before reaching the threshold, since spikereset 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
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.,
EGLIF 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 singlecell recordings, GoCs show spontaneous firing around 8 Hz, a nearlylinear inputoutput relationship (about 0.25 Hz/pA), inputdependent spikefrequency adaptation when depolarization is maintained, rebound bursting after hyperpolarization, phaseresetting, subthreshold selfsustained oscillations and resonance in the theta band (around 3–6 Hz) (Forti et al.,
As reported in section Model parameters, the values of the electrophysiological parameters were taken from literature (Forti et al.,
Electrophysiological parameters of the cerebellar GoC taken from literature.
145 pF  Solinas et al., 

44 ms  Solinas et al., 

−62 mV  Solinas et al., 

2 ms  Tripathy et al., 

−75 mV  Tripathy et al., 

−55 mV  Forti et al., 
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.,
Inputoutput relationship for GoC model optimization with corresponding electroresponsive properties (Forti et al.,
Autorhythm  
adaptation 

adaptation 

adaptation 

Latency of 1st spike lower than 0.5·(1/ 
Postinhibitory rebound bursting  
Δ 
In order to account for the wholeset of GoC electrophysiological properties, the cost function included all the terms reported in
The optimal parameter set was chosen as the median of the final values over the 5 optimization runs with different random parameter initialization.
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 multistep protocol with the same input currents used for optimization, taken from literature (see Table
Secondly, we added stimulation patterns required to evaluate the emergence of model features (autorhythmicity, rebound bursting,
an initial zerocurrent 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 1s zerocurrent phases, to test intrinsic excitability;
two zeroinput 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 phasereset mechanism;
two 1s hyperpolarizing phases with inhibitory current −200 and −100 pA, followed by a 1s zerostimulation phase, where evaluating rebound bursting properties;
a sequence of 5 steps, each of amplitude 600 pA, lasting 30 ms, at increasing frequencies: 0.523.556.37.7101215 Hz, to evaluate resonance.
The resulting total duration of the whole protocol was 48.93 s (Figure
Finally, we assessed the capability of the model to intrinsically generate selfsustained
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 conductancebased, with a spiketriggered change of the synaptic conductance
This change of synapse conductance caused a change in the
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 EGLIF with currentbased synapses, which are less realistic but also require a lower computational load and can be suitable for largescale SNN simulations (Cavallari et al.,
The outcome of simulations using the long multiphase protocol was compared to real data acquired
The experiments have been conducted on 16to21daysold (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.,
Wholecell patchclamp recordings were performed with Multiclamp 700B (3dB; cutoff frequency 10 kHz), sampled with Digidata 1550 interface, and analyzed offline 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, ATPNa_{2} 4, GTPNa_{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 lowpass filtered at 10 kHz and acquired at 50 kHz. The stability of wholecell recordings can be influenced by modification of series resistance (R_{s}). To ensure that R_{s} remained stable during recordings, passive electrodecell parameters were monitored throughout the experiments.
In each recording, once in the wholecell configuration, the current transients elicited by 10 mV hyperpolarizing pulses from the holding potential of −70 mV in voltageclamp 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 steadystate current flowing after termination of the transient (195.9 ± 97.8 GΩ). The 3 dB cutoff frequency of the electrodecell 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 1s steps of current (from −200 to 600 pA in a 100pA increment). Resonance was investigated by applying sequences of 30ms 600 pA current steps repeated 5 times at different frequencies (range from 2 to 10 Hz) or injecting a single 0.5ms step of 4 nA. Finally, the cells were maintained at their resting membrane potential for evaluating subthreshold
Using an automatic spike detection algorithm, spike times were extracted from the experimental recordings, and used for electrophysiological feature extraction (see section Feature Extraction).
To evaluate GoC electroresponsive behavior, multiple parameters were computed according to (Solinas et al.,
the tonic firing rate,
the coefficient of variation of interspike intervals (CV_{ISI}) was measured during the 10s zeroinput step to quantify the irregularity of firing;
the firing rate,
the
the
the parameter
latency and frequency of 1^{st} spike were measured in the rebound burst after hyperpolarization (
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 (Gandolfi et al.,
Parameter values for simulations and recordings are reported as mean ±
In addition, for quantifying experimental subthreshold oscillations, the power spectral density of 1 s
The EGLIF neuron is a linear monocompartment neuron model using three differential equations for membrane potential (
In the present development, an optimization strategy using a gradientbased 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 EGLIF 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 multiphase 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.
In order to reproduce nondamped membrane potential oscillations,
Parameter optimization.
The optimized model was tested by 10 simulations in PyNEST (Eppler et al.,
Using the parameters resulting from the optimization process, the model was able to generate a linear relationship between
PyNEST simulations.
The model response at increasing
After a hyperpolarizing current step (200 pA), the model showed a postinhibitory 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
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 depolarizationinduced excitation and spikefrequency adaptation in response to positive current steps, rebound doublet bursts after negative current steps (Figures
Responses to current steps.
Excitability and adaptation.
Validation and prediction (e.g., phase reset with higher pulse) of rebound bursting and phase reset phenomena.
after 
20.5 ± 8.5  21.9 ± 5.6  
after 
24.1 ± 8.3  26.3 ± 10.9  
after 
36.8 ± 17.2  30.4 ± 3  
after 
44.2 ± 18.7  45.6 ± 7.7  
after pulse = 4000 nA  1.65 ± 0.38  1.84 ± 0.09  
after pulse = 6800 nA  not acquired  1.9 ± 0.16 
A linear fitting on the experimental data resulted in the same
An important phenomenon evident in GoCs (as well as in some other brain neurons) is
Phase reset. Phase reset in a series of 10 experimental acquisitions
Finally, another property evident in GoCs and often observed in central nervous system neurons is the presence of endogenous
Subthreshold oscillations and resonance.
As shown in Figure
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 conductancebased, with an alphashaped conductance change. As shown in
This work reports the development, optimization and testing of the
By exploiting its analytical solution, the EGLIF GoC model was optimized using gradient descent minimization methods. This allowed to fasttune a unique set of parameters generating the appropriate spiking responses to various input patterns. The cost function was designed to evaluate subthreshold membrane potential dynamics using the integral of membrane potential all over the interspike interval: in this way, the cost function was still differentiable, while taking into account the full history of the signal preceding the spike. Suprathreshold dynamics may not be relevant since EGLIF used a spikeupdatereset 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
Compared to other optimization strategies, which are based on semiautomatic fitting on a training set of neural traces or spiking patterns (Pozzorini et al.,
EGLIF 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,
First, EGLIF reproduces subthreshold responses and spike timing rather than the shape of action potentials. 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.,
Secondly, the EGLIF analytical tractability can be exploited to implement versions designed for timedriven simulations, where the exact spike time can be computed through iterative methods like bisection (Hanuschkin et al.,
Thirdly, EGLIF exploits an adaptive current coupled with membrane potential to model spiking responses driven by adaptation mechanisms, like spikefrequency adaptation and postinhibitory rebound bursting, rather than a sliding threshold depending on the actual membrane voltage (Mihalaş and Niebur,
Finally, for testing the neuron with synaptic inputs, conductancebased synaptic receptors were embedded into the model. It should be noted that currentbased synapses may be used to further improve computational efficiency at the expense of biophysical realism (Cavallari et al.,
While appropriately reproducing the physiological spiking patterns in response to various input stimuli, the EGLIF 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 secondorder dynamics, the subthreshold EGLIF
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 highamplitude After Hyperpolarization potential at the end of depolarizing steps does not overlap the experimental behavior. The reason is the high steadystate value of
Given the specific design of EGLIF (see above), it is possible to trace the relationship between EGLIF 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.,
Correspondence between subcellular mechanisms that generated specific electrophysiological properties and how they were simplified in EGLIF model elements.
I_{Na−t} ↑ / I_{K−V} ↓ balance  Action potential  Digital 
I_{h} ↑I_{Na−p} ↑ / I_{K−slow} ↓ balanceI_{Ca−HVA} ↑/ I_{K−AHP} ↓ balance  Autorhythmicity  
I_{Na−r} ↑ and I_{K−A} ↓  Depolarizationinduced burst  
I_{Ca−HVA} ↑/ I_{K−AHP} ↓ balanceI_{K−slow} ↓  Spikefrequency adaptation  
I_{Ca−HVA} ↑/ I_{K−AHP} ↓ balance  Phasereset  
I_{h} ↑I_{Ca−LVA} ↑  Postinhibitory rebound burst  
I_{K−slow} ↓ and I_{Na−p} ↑  Subthreshold oscillations  
I_{K−slow} ↓ and I_{Na−p} ↑  Resonance 
While, in order to simultaneously generate 8 different electrophysiological properties (action potential, autorhythmicity, depolarizationinduced burst, postinhibitory rebound burst, spikefrequency adaptation, phasereset, 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}), EGLIF has just 3 currents (
This work shows that it is possible to represent complex neuronal firing dynamics through a monocompartment neuron model by updating the GLIF into EGLIF at the modest computational expense of 3 ODEs. Yet there is a remarkable efficiency gain with respect to realistic multicompartmental models. For example, compared to the realistic version of the GoC model (Solinas et al.,
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.
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.
This project has been developed within the
The Supplementary Material for this article can be found online at: