Realistic modeling allows a faithful reconstruction of neuronal excitable properties based on the principles of neuronal biophysics (Koch, 1999; De Schutter, 2001). This approach requires a precise representation of the electrotonic structure of neurons and of their ionic membrane mechanisms through a variety of ionic channels. Thus, realistic modeling is, in essence, a modern expansion of the approach developed by Hodgkin and Huxley (1952) for the action potential (AP) in the squid giant axon. Despite the amount of parameters populating realistic models is huge, most of them are constrained by experimental measurements and the parameters that remain free are basically the maximum ionic conductances (G_i-max). Experimentally, G_i-max can rarely be measured reliably due to space-clamp problems. Moreover, the ionic current identified by electrophysiological and pharmacological tools often reflects activation of a blend of different channel molecules rather than correspond to a single type of genetically identified channel. Thus, once the ionic conductances in a neuron have been identified and represented in a Hodgkin-Huxley-like style (HH), what is usually done is to empirically adjust their G_i-max until matching the neuronal firing pattern. This “iterative multiparametric matching” with large experimental datasets can lead to precise models (e.g., see D'Angelo et al., 2001; Solinas et al., 2007a,b; Diwakar et al., 2009; Subramaniyam et al., 2014; Masoli et al., 2015), but it is slow and laborious.

A recent technique that allows rapid and automatic parameter estimation is based on multi-objective evolutionary algorithms (MOEA), such as the “Non-dominated Sorting Genetic Algorithms-II” (NSGA-II; Deb et al., 2002) and “Indicator-Based Evolutionary Algorithm” (IBEA) (Zitzler and Künzli, 2004). These are based on a genetic approach, in which the unknown parameters are treated like the genes forming a chromosome (Deb et al., 2002; Zitzler and Künzli, 2004; Druckmann et al., 2007, 2008). To find the solutions, these algorithms require specific target parameters, for example “features” extracted from experimental traces. The features are used to define the basic properties of APs (such as amplitude and hyperpolarization depth), of neuron discharge (e.g., frequency and first-spike delay) and of subthreshold responses, in order to assess the fitness function. From the features one or multiple “objective” functions can be computed, that have to be minimized simultaneously and represent the “fitness” of the individuals. Each individual carrying a specific combination of parameters is part of a population. In each generation, the algorithm performs a ranking of the individuals of the population and removes a predefined number of the worst individuals. This makes room for new individuals derived from the retained individuals using genetic principles (e.g., cross-over, mutation, elitism). Retained and new individuals make up the next generation. While evolutionary algorithms can lead to fast approximation of neuronal firing patterns, some G_i-max combinations could be non-physiological (e.g., non-unique solutions). Therefore, a stringent test is required to assess whether, among the solutions provided by the optimization procedure, there are (at least) some that match the biological G_i-max distribution.

In order to face the issue, we need a neuron in which the G_i-max values have been precisely estimated in electrophysiological recordings providing stringent constraints to the mechanisms of AP generation. The GrC offers this unique opportunity. The GrC is one of the smallest neurons of the brain and this has allowed to achieve exceptionally good voltage-clamp conditions leading to a precise determination of ionic current gating kinetics and G_i-max values. These include the high-voltage activated Ca²⁺ current (Ca-HVA; Rossi et al., 1994), the Na⁺ current (Na, Nap, Nar) (Magistretti et al., 2006; Goldfarb et al., 2007; Dover et al., 2010), the inward rectifier K⁺ current (Kir) (Rossi et al., 1998, 2006), the A-type K⁺ current (KA), the voltage-dependent outward-rectifier K⁺ current (KV), and the K⁺ calcium dependent (KCa) (Bardoni and Belluzzi, 1994), the M-type slow K⁺ current (Kslow; D'Angelo et al., 2001). In addition, some models that have been previously developed using iterative multiparametric matching can be used for comparison (Gabbiani et al., 1994; D'Angelo et al., 2001; Nieus et al., 2006; Diwakar et al., 2009). Therefore, the question is whether advanced optimization procedures can capture the whole set of GrC properties through a set of maximum ionic G_i-max values compatible with those measured experimentally.

Here we show that an automatic parameter estimation procedure, the Optimizer Framework (OF) which is based on Druckmann et al. (2007, 2008, 2011) and improved to run with IBEA (Zitzler and Künzli, 2004), can indeed provide GrC models with a biologically plausible set of G_i-max values. These models can predict electroresponsive properties like inward rectification, near-threshold oscillations, theta-frequency resonance and AP conduction velocity that were not set as features. These results indicate that the OF generates biophysically accurate models endowed with appropriate ionic mechanism, providing the basis for reconstructing large-scale neuronal networks operating with arbitrary firing patterns.

In this paper, a pipeline was developed to generate families of GrC mono-compartmental and multi-compartmental models, to optimize their G_i-max complement, and to validate the models through the simulation of electroresponsive properties not considered for model construction. The features used as templates were extracted from GrC spike discharges under the assumption that these contain all the information required to optimize G_i-max values. The present models can be defined “realistic” as far as they reflect a modeling strategy that implements neuronal membranes with biophysically-detailed mechanisms (see discussion in De Schutter, 2001; Santamaria et al., 2007; D'Angelo et al., 2016).

In this work we used OF (Druckmann et al., 2007, 2008; Van Geit et al., 2016) to generate optimized GrC models that accounted for the ionic currents of real GrCs and to compare them to previous mono-compartmental [D'Angelo et al. (2001); updated in Nieus et al. (2006); Solinas et al. (2010)] and multi-compartmental GrC models [Diwakar et al. (2009); updated in Dover et al., 2016] tuned by iterative multiparametric matching. The features used here for optimization addressed GrC resting membrane potential, spike shape and firing pattern (Table 1). The matching of optimized model responses with the template is illustrated in Figure 1A. The correspondence of G_i-max in the optimized models with those measured experimentally in GrCs is shown in Figure 1B, along with a comparison of parameters in the optimized models with respect to previous models. Strikingly, the experimental and modeled G_i-max values (both for the optimized and previous models) showed an almost linear correspondence across 4 orders of magnitude, ranging from about 10⁻² to 10² mS/cm². This precise correspondence implied that in most cases model parameters differed from the experimental ones by <20%, and the same was true for the difference between model parameters obtained using OF or procedures of iterative multiparametric matching. The largest variations were observed for the multicompartmental models, probably reflecting an imperfect knowledge about the localization of some ionic channels. The solutions converged rapidly and were almost at steady-state already after about 20 iterations for either models (Figure 1C). It should be noted that the non-0 steady-state level demonstrated a continuous remixing of genes over the entire population and the maintenance of genetic variability.

This paper shows that feature extraction from firing pattern templates combined with IBEA optimization/selection methods (OF; Druckmann et al., 2007, 2008) allow to generate neuronal models with ionic mechanisms that closely reflect those of real cerebellar GrCs (cf. Figure 1). These models are capable of uncovering properties that were not evident in the templates and derive uniquely from the appropriateness of the whole set of G_i-max values. To our knowledge, this is amongst the first OF applications to neurons other than those of the neocortex and hippocampus (for a recent application see Eyal et al., 2016).

The use of OF to parameterize cerebellar GrC models, for which a previous experimental determination of G_max−i is available (D'Angelo et al., 2001; Diwakar et al., 2009), illustrates some general aspects of the procedure. First, firing patterns recorded from the soma were sufficient to provide the information needed to fully parameterize the whole set of G_i-max values. The average G_i-max fluctuation was 12.3% for the mono-compartment models and 19.7% for the multi-compartment models. In both cases, most channels G_i-max fluctuations are between the variation limits suggested as acceptable for realistic modeling (≤ ± 20%, Bower and Beeman, 2007). Secondly, the OF neuron models proved correct against a set of properties used as features, including the spike shape and the dependence of firing frequency and first-spike delay on the injected current. This is akin with the empirical practice of matching features by adjusting model parameters by trial-and-error iterations and parametric fitting (e.g., Solinas et al., 2007a,b; Subramaniyam et al., 2014; Masoli et al., 2015). Thirdly, the model predicted properties not represented in the firing pattern and not used as features, including inward rectification, resonance and near-threshold oscillations. This is a proof of the plausibility of the underlying biophysical mechanisms, that allowed the whole set of neuronal properties to emerge. Consistently, G_i-max and membrane currents were similar to those recorded experimentally. Of particular relevance are the size and temporal dynamics of small ionic currents, that do not contribute remarkably to the total current (I_Nap and I_Kslow are 3 orders of magnitude smaller than the I_Na and I_KV currents) but were non-etheless matched by OF. Finally, with the GrC multi-compartmental model, indications about ionic channel localization had to be anticipated in order to allow OF to appropriately reconstruct the spike generation mechanism. This was an essential constraint, otherwise the models may still be able to generate a seemingly good firing pattern but the underlying mechanism would be incorrect (for example Na⁺ channels may well-generate spikes while staying in the soma, data not shown).

Beside the effectiveness of this optimization approach, there are some aspects that deserve attention. The identification of models with plausible biophysical properties required a validation going beyond the matching of features characterizing the firing pattern. For example, for the GrC models to be valid, the presence of resonance and near-threshold oscillation needs to be ascertained even if it is not included into the features. Therefore, following a selection process based on quality indicators (Sutskever et al., 2015), the final validation can further limit realistic models to a subset of all optimized models. Validation proved critical for both the multi-compartment and the mono-compartment models, with the percentage of neurons showing a plausible firing and spike shape being 13 and 6.5%, respectively. The rendering and requirements for validation in very complex neuron models (e.g., in the pyramidal neuron; Traub et al., 1991; Druckmann et al., 2011) or in the Purkinje cell neuron models (De Schutter and Bower, 1994a,b; Achard and De Schutter, 2006; Masoli et al., 2015) needs to be evaluated.

Given that all granule cells normally show stereotyped firing patterns under a unified mechanism, the fact that not all models in the final population do the same raises a relevant issue. Is it possible that, in biology, specific mechanisms constrain the solution toward the ionic channel asset needed to reach a specific firing pattern? These mechanisms may reside in yet undiscovered feed-back biochemical processes regulating G_i-max values through channel expression or modulation. Moreover, the fact that we have accepted only those solutions conforming to the most typical or “canonical” description of GrCs might have restricted the acceptance criteria. The OF may actually be able to predict model variants, whose correspondence with existing biological states is currently unclear and requires experimental assessment.

Special consideration deserves the AP generating process. The careful investigation of GrC electrophysiology has revealed that Na⁺ channels are almost absent from soma and are maximally concentrated in the AIS (Magistretti et al., 2006; Goldfarb et al., 2007). In the present model, AP passive back propagation to soma and dendrites occurred with proper delay and the active propagation in the axon occurred at the proper speed. These can again be considered as emerging properties depending on model structure and ionic channel distribution and density. A very recent study showed that the GrC axon contains Na⁺ channels that have different gating kinetics from those of the AIS and that the specific axon leak resistance is several times higher than in AIS, soma and dendrites (Dover et al., 2016). This allows the axon to lower by several time the energy consumption compared to the classical HH mechanism. These special properties of the axon, that have not been considered here, may be used as further constraints in the optimization process.

One may speculate on the way the optimization algorithm predicted the emerging phenomena. Probably, as far as inward rectification is concerned, the Kir conductance was set by extracting information from the current needed to move from rest to AP threshold and from the anomalous slow-down in first spike delay as the current injection was increased. And since the same feature is also controlled by KA, this could have generated some uncertainty in the conjoint estimation of these two parameters. Likewise, the I_Nap/I_Kslow balance was probably determined by setting the firing threshold and the subsequent f/I relationship. Therefore, given representative templates of the firing patterns, the OF could predict the cell response in the subthreshold and near-threshold functional regimes, that were not considered as features. The same consideration applies to resonance, which depends on I_Kslow and is amplified by I_Nap.

In conclusion, the OF using IBEA provided an objective automatic strategy for reconstructing the biophysical properties of neurons through a realistic set of ionic currents, reducing optimization times by orders of magnitude compared to traditional operator-guided procedures (days vs. months or years). The OF-derived models required a validation accounting for response patterns not used for model construction. This validation, by being itself based on biological constraints, does not introduce any arbitrary selections. The unsupervised optimization through OF confirmed the precision of multi-parametric matching procedures used previously for the same neurons. The ability of the OF models to respond to low-frequency oscillations and bursts of various frequency and duration makes them suitable for reconstructing large-scale models of neuronal microcircuits.

SM wrote the model codes and the performed most of the simulations with the contribution of MR. MS performed the experiments. WV and FS supported model implementation and the discussion of principles. ED coordinated the work and wrote the paper with the contribution of all other Authors.