Evolving Simple Models of Diverse Intrinsic Dynamics in Hippocampal Neuron Types

Introduction

In the last decade, several projects have built large-scale models of brain regions in an attempt to advance our understanding of how the nervous system functions (Izhikevich and Edelman, 2008; Eliasmith et al., 2012; Markram et al., 2015; Hendrickson et al., 2016). The biological realism in these models has been captured in varying levels of detail. One of the characterizing features of biological neural networks is the diversity observed in the intrinsic dynamics of individual neurons. This diversity likely contributes to the emergent properties of neural networks and, consequently, plays a major role in the information processing in the nervous system (Padmanabhan and Urban, 2010; Tripathy et al., 2013; Pozzorini et al., 2015). Therefore, a biologically realistic large-scale network model of a brain region should take into account intrinsic behavioral diversities both within and between neuron types.

Hippocampome.org is a comprehensive knowledgebase of 122 morphologically identified neuron types in the rodent hippocampal formation (Wheeler et al., 2015). One of the motivations behind developing this knowledgebase was to create a real-scale computational model of the entire hippocampus. Toward achieving this goal, we aim to create individual neuronal models using the electrophysiological and spike pattern properties of neuron types available at Hippocampome.org. In deciding which modeling system to use, we considered simulation costs. High simulation costs of biophysically detailed Hodgkin-Huxley-type neuronal models often impose limits on the scale of network models. Conversely, simpler models, such as leaky integrate-and-fire neurons, cannot capture the wide range of dynamics observed in the hippocampus. Models such as Izhikevich (Izhikevich, 2003) and Adaptive Exponential Integrate-and-Fire (AdEx) (Brette and Gerstner, 2005) have been shown to qualitatively reproduce various firing pattern classes observed experimentally in real neurons, while still being computationally efficient. Therefore, these simpler models with lower simulation costs allow large-scale modeling of biological neural networks in a computationally efficient manner. In this work, we create Izhikevich Models (IMs) that reproduce quantitatively comparable features of various hippocampal spike pattern classes through parameters optimization.

The dynamics of Izhikevich models are highly non-linear and error landscapes that are defined over the resulting parameter spaces typically exhibit properties that make them difficult to optimize, such as multiple local optima. As such, several studies have turned to non-convex, derivative-free optimization methods such as evolutionary algorithms (EAs) to fit a neuronal model's responses to experimentally recorded voltage traces. The models used in these studies range from simple spiking models such as AdEx (Rossant et al., 2010, 2011; Lynch and Houghton, 2015) to biophysically detailed Hodgkin-Huxley type models with multiple compartments (Gerken et al., 2005; Keren et al., 2005; Druckmann et al., 2007; Van Geit et al., 2007). Previous studies have also used various techniques such as a feature-based error function (Druckmann et al., 2007) and a phase plane trajectory density method (Van Geit et al., 2008) to create the error landscape for the EA search. Rössert et al. (2016) created an approach to simplify morphologically detailed microcircuit models to their point-neuron counterparts by applying soma-synaptic correction (to account for dendritic attenuation and delay) and constraining Generalized Integrate-and-Fire neurons around an in vivo-like working point. Rounds et al. (2016) used EAs to match firing rates of IMs in a network to experimental recordings in the retrosplenial cortex. However, to our knowledge, optimization techniques have not been successfully used to fit intrinsic IM responses to experimental data. On benchmark optimization tests, the IM showed poor performance compared to other simple models (Rossant et al., 2010, 2011; Lynch and Houghton, 2015). This might be due to the failure to identify an appropriate EA configuration such as the choice of error function and variation operators that are well-suited for the IM parameter space.

Apart from its capability to quantitatively fit IM's responses to experimental voltage traces, the novelty of our automated modeling framework is the integration of spike pattern classification protocols. Previous work (Komendantov et al., in review) identified 23 distinct spike pattern classes overall, among the 89 morphologically distinct hippocampal neuron types in Hippocampome.org for which experimental recordings were available. A behavior for a certain neuron type was defined based on the set of all experimentally recorded spike patterns. If a neuron type exhibited spike patterns of more than one class under different experimental conditions (e.g., bursting and regular spiking for different current stimulation strengths), it was marked as a multi-behavior type. In contrast, a neuron type was marked as a single-behavior type, if all spike patterns recorded from the same neuron under different experimental conditions fell into the same qualitative class. Neuron types with only a single experimentally recorded spike pattern were also marked as single-behavior.

This article presents the modeling approach and results for single-behavior neuron types. We report at least one example for each of the nine distinct single-behavior types, with the goal of illustrating both the approach and the IM's ability to quantitatively reproduce a variety of neuronal behaviors observed in the hippocampus. The single behaviors reported here include spiking with and without frequency adaptation, delayed spiking, bursting, and intermittent spiking or stuttering. In addition to simple point-neuron (single-compartment) models, multi-compartment IMs were created, where the number of compartments varied from two to four depending on the dendritic invasion of a neuron type across hippocampal layers. For example, the somata of hippocampal pyramidal cells in the principal layer extend basal dendrites in the oriens layer and apical trees in the radiatum layer that reach to the lacunosum-moleculare. Thus, these neurons can be represented as 4 compartments, one for each layer (as illustrated in the Methods below). This stratification is important because it segregates the synaptic inputs: distal lacunosum-moleculare dendrites, for example, are the targets of entorhinal projections, while dendrites in radiatum receive intra-hippocampal connections. Although finer morphological variability observed across various neuron types may also contribute to network dynamics, compartmentalized dendritic integration of distinct laminar inputs is likely to play a crucial computational role in cortical circuits. Furthermore, dendrites located in separate layers typically have different active and passive properties from each other and from the soma. A previous large-scale model of the thalamo-cortical system used multi-compartment IMs (Izhikevich and Edelman, 2008). However, that model did not capture the signal transmission properties between the dendrites and soma in a biologically accurate way. In addition, the dendritic compartments did not reflect the appropriate balance of active and passive properties. Another novelty of our automated modeling approach is its capability to create accurate dendritic representations in the multi-compartment IMs. Our dendritic compartments exhibit generally known active and passive properties of the dendrites of real neurons.

Materials and Methods

In this article, a certain spike pattern class will be used to denote a neuron type's “behavior,” since all the neuron types discussed here were examples of single-behavior types. It is worth mentioning that 14 of the 23 distinct spike pattern classes observed in the hippocampus are part of the multi-behavior types and, hence, not reported in this article. Modeling multi-behavior cases requires a different approach, which we are pursuing but remains beyond the scope of this article.

Spiking Model

We reproduced spike patterns by using the nine-parameter variant of the IM (Izhikevich, 2007) because we found that the EA could reliably find better solutions with this IM than the originally proposed four-parameter formalism. IMs have been shown to reproduce qualitatively many spike patterns observed in biological neurons. The state variables membrane voltage (V) and membrane recovery variable (U) govern this two-dimensional system. The recovery variable U approximates the channel kinetics of Hodgkin-Huxley type models (Hodgkin and Huxley, 1952), making it computationally much cheaper to simulate. Parameter “a” is the time constant for the recovery variable U. Parameter “b” defines the degree of coupling between the state variables V and U. Parameters “b” and “a” collectively determine whether the model is an integrator or resonator (Izhikevich, 2001). Parameters “V_min” and “d” are after-spike reset values for V and U, respectively. Parameter “k” defines the shape of the spike upstroke, and V_peak defines the spike cutoff value. Parameters V_r and V_t are resting and threshold voltages, respectively, and C is cell capacitance.

In addition, we created multi-compartment (MC) models for each neuron type based on the dendritic invasion across the hippocampal layers (Figure 1), whereas each compartment represents the part of the dendritic tree present in a given layer. Compartments were coupled using an asymmetric mechanism. For example, the MC layout for the CA2 pyramidal neuron type depicted in Figure 1 defines the compartment-specific coupling currents as follows:

Here, Stratum Pyramidale (SP) denotes somatic compartment, and Stratum Oriens (SO), Stratum Radiatum (SR) and Stratum Lacunosum-Moleculare (SLM) (Figure 1) are dendritic compartments. I_SP is the total coupling current at compartment SP, which results from the differences in the voltage between SP (V_SP) and SO (V_SO), and between SP and SR (V_SR). G₁and G₂ are coupling constants, and P₁ and P₂ (with values between 0.01 and 0.99) determine the degree of asymmetry in the coupling, where a value of 0.5 denotes symmetric coupling.

Identification of Spike Pattern Classes

To classify the model behaviors, we used the same protocol developed for identifying various transient and steady-state elements of experimentally recorded spike patterns from the hippocampus (Komendantov et al., in review). Transient elements include: Delay (D), if the latency to spike is sufficiently long; Adapting Spiking (ASP), if the spike frequency decreases over time; Transient Stuttering (TSTUT), if a quiescent period follows a cluster of high frequency spikes; and Transient Slow-Wave Bursting (TSWB), if TSTUT is followed by a slow after-hyperpolarizing potential. Steady-state elements include: Silence (SLN), if the quiescence following the last spike is sufficiently long; Non-Adapting Spiking (NASP), if no spike frequency adaptation is identified in non-interrupted firing; Persistent Stuttering (PSTUT), if at least one sufficiently long quiescent period separates two clusters of high frequency spikes; and Persistent Slow-Wave Bursting (PSWB) if slow after-hyperpolarizing potential is present in an otherwise PSTUT pattern.

Given a sufficiently long duration of input current, transients will always be followed by a steady-state in a spike pattern. For example, ASP followed by NASP was identified in the spike pattern experimentally recorded from a CA3 Basket-CCK neuron (Gulyás et al., 2010) and this pattern is an instance of the class ASP.NASP (Figure 2A). The identified class of an experimentally recorded spike pattern represented the target class for model spike pattern. Thus, the criteria that defined a target class were used to validate the model behavior under the given input current. For details on the spike pattern classification criteria, see Komendantov et al. (in review) and Hippocampome.org¹

Evolutionary Optimization of Model Parameters

Many varieties of EAs exist along with numerous ways of implementing their specific components (De Jong, 2006). We employed a non-overlapping generational model of evolution and used elitism to ensure that the best individuals were always preserved in the population. In this section, we describe the specific components of the EA.

EA Configuration

Each individual in the evolutionary population consisted of a complete configuration of the IM we are seeking to tune. We represented these configurations as vectors of floating-points. During the search of the parameter space, we bounded each value within an allowed range. Choosing a biophysically reasonable range for each parameter up front has a significant effect on the efficiency of the optimization procedure, and through preliminary EA runs, we found that different behavior classes required slightly different ranges for some parameters.

When tuning a single-compartment model, each parameter vector contained 9 + n_I elements (genes), representing the 9 parameters of the IM (k, a, b, d, C, V_r, V_t, V_peak, V_min) and n_I input currents. n_I equaled the number of voltage traces (which were recorded for different input currents) a model was fit to. A small range encompassing each experimentally injected input current (I_exp.), [I_exp. – 10, I_exp. + 10] pA was included in the EA search. By allowing the EA to search within a small range, we boosted the exploration and identified more optimal points (across multiple EA trials) that are very similar in the phenotype. This design also helped to achieve more reliable fitting in cases where a single model was fit to multiple traces (see section Quantitative Comparison of Spike Pattern Features). In rare cases, where the experimental input current was unknown, an unbounded range of [50, 800] pA was included for the search. The multi-compartment models had a larger number of parameters: if c is the number of compartments, we require 8c + 1 parameters representing c compartments (the parameter V_r is the same for all c compartments), plus 2(c−1) parameters representing coupling parameters for consecutive compartment pairs, and n_I input currents.

We first initialized a population of these vectors by sampling uniformly from within each parameter's allowed range. We used a fixed population size of 120 individuals for single-compartment IMs and of 400 for multi-compartment IMs. An exception was the 4-compartment ASP.SLN fast-spiking model, which we found was easier to optimize with a larger population size (800). After initializing the population, and at each generational cycle thereafter, we immediately selected the 10% of the population with the lowest error to survive to the next generation.

We filled the remainder of the child population by selecting parents via binary tournament and by applying two-point crossover and a mutation operator. Each gene had a probability of mutation between 0.1 and 0.3. For the parameters d, C, G, and I, we applied an integer random-walk mutation: when selecting one of those genes for mutation, an integer increment or decrement was applied with equal probability. All remaining parameters were mutated by reset: a new floating-point value was randomly chosen out of that gene's allowed range. The EA was run until a maximum number of generations was reached. This number varied between 500 and 5,000 depending on the number of compartments in the model and the class of behavior the model was fit to.

Error Function

We employed a feature-based error function to quantitatively reproduce the spike pattern. Features for more than 250 experimentally recorded spike patterns are available at Hippocampome.org. These features include first-spike latency (fsl), post-spike silence (pss), spike frequency adaptation parameters (sfa), burst width (bw), post-burst interval (pbi), and rebound voltage amplitude (rbv) (Figure 2). Spike frequency adaptation (sfa) was quantified as previously detailed (Komendantov et al., in review) with a piecewise linear regression on the inter-spike intervals (ISIs) by extracting the parameters of linear fits such as slopes (m) and Y-intercepts (c) (Figure 2B). The error in the model sfa was calculated as follows: the two parameters of linear fit (for NASP and ASP. class) or three parameters of piecewise linear fits (for ASP.NASP class) obtained by plotting ISI's against their latencies to the second spike were compared between experimental and model spike patterns. In addition, the number of ISIs (nisi) corresponding to a linear fit was also compared. For a bursting/stuttering class, the number of bursts (nbs) and the number of spikes (nspikes) within each burst were included.

Spike pattern classification protocols were also incorporated into the error function by dynamically assigning different weight factors to different features. This reduced the number of generations required for the EA to find an acceptable solution. In addition, for certain spike pattern classes, this approach more reliably found solutions across multiple stochastic trials. The error function was defined as:

S: {fsl, pss, m, c, nisi} for continuous spiking, and S: {fsl, pss, nbs, bw, pbi, nspikes} for interrupted spiking (Figure 2).

Using the spike pattern classification protocols, the qualitative class of a candidate model's spike pattern was first identified. Then, the weight factor W_f was calculated for each feature by comparing the target class with the model spike pattern class. During the EA search, each feature weight changed based on that feature's distance from the target class boundary. These class boundaries are given by the set of criteria that define that class (see Supplementary Material section 1 for pseudocode description of feature weight calculation). This accelerated the search during earlier generations of an EA, when many candidate solutions were outside the target class boundary (fast-exploitation toward narrow regions of interest). Once the population began converging within the target class boundary, this dynamic weight assignment scheme allowed slower exploration and ensured heterogeneity among the best models from within a class (see section Variabilities in the Intrinsic Properties within a Neuron Type).

In order to identify several possibilities, we ran a total of one thousand EA instances for each neuron type, yielding several best models due to the stochastic nature of the EA (different initial populations, stochasticity in variation operators, and selection pressure). At the end of an EA search, the best model was accepted only if its spike pattern class under the given input current matched the target class. If a single-compartment IM failed to reproduce a firing pattern class, two identical compartments were symmetrically coupled. We noticed that coupling effects enriched IM dynamics and were useful to reproduce certain subclasses of stuttering behavior (see Results section for examples).

The parameters of MC models were optimized such that the somatic compartment reproduced features of the experimentally recorded voltage traces both qualitatively and quantitatively using the same techniques discussed previously for single-compartment models. Furthermore, four additional constraints were enforced in order to capture the known general active and passive properties of dendrites in the additional compartments. Unlike the somatic compartment constraints, all MC models shared the same dendritic constraints because of the lack of sufficient experimental dendritic voltage recordings. These general constraints include excitability and input resistance of dendrites relative to the soma as well as forward propagation of spikes and subthreshold signals.

Firstly, the dendritic compartments in MC models were constrained to be less excitable than the somatic compartment when they were decoupled. The minimum depolarizing current (I^rheo) required to elicit a spike at a compartment was used as the measure of its excitability (Aou et al., 1992). During the EA search, a ramp current rather than step currents was used to measure compartment excitabilities. This avoided the need for a local search for the minimum step current magnitude required to elicit a spike in each compartment. To reduce capacitive effects in measuring excitability, the ramp current had minimal slopes and high resolution of discrete increments (+0.1 pA/1 ms). Secondly, the decoupled dendritic compartments in MC models were constrained to have higher input resistances than the somatic compartment. The amplitudes of steady state voltage deflections from resting voltage (V^def) during a strong hyperpolarizing current input were compared between compartments to measure their relative input resistances. The spike propagation rate (R) was defined as the ratio between the number of spikes observed at the destination compartment and the number of spikes initiated at the source compartment. A few hundred excitatory synapses were stimulated at a dendritic compartment for spike initiation. On the other hand, a single AMPA synapse was stimulated at a dendritic compartment and the amplitude of the excitatory post-synaptic potential (EPSP) was measured at the somatic compartment. A range of (0.1, 0.9) mV was enforced for the EPSP amplitude. All synapses used a value of 10 for the weight, and this value is based on the range used for the multicompartment models by Izhikevich and Edelman (2008). The following errors were calculated for each dendritic compartment and added to the somatic spike pattern error described earlier:

Model and Algorithm Implementations

We used the open-source Java-based evolutionary computation system ECJ (Luke et al., 2015) to tune IM parameters. Single compartment models were simulated using the Apache Commons Mathematics Library². The MC models with up to 39 open parameters were tuned using the parameter tuning interface of CARLsim, an open-source high performance GPU-based spiking neural network simulator (Beyeler et al., 2015). The EA and the single compartment model simulations were run on distributed CPU nodes, and the MC models were run on the GPU nodes available at the Office of Research Computing at George Mason University. All scripts necessary to reproduce the results reported in this article are publicly available³.

Results

Models of Distinct Single Behavior Types

A total of 33 of 122 neuron types in Hippocampome.org version 1.0 (Wheeler et al., 2015) exhibit single behavior. Nine distinct single-behavior classes exist among these neuron types, and, in this article, we present at least one model for each of those classes. It is worth mentioning that different neuron types that exhibit the same qualitative behavior class typically exhibit different quantitative features and excitability levels. Figure 3 illustrates an exemplar neuron type for each of the nine distinct single-behavior classes and the corresponding best model from all EA trials. Our simple models were able to reproduce quantitatively comparable spike pattern features for all these classes (see section Quantitative Comparison of Spike Pattern Features). While earlier models reproduced seven qualitatively different classes of spike patterns (Izhikevich, 2003), our models capture the spike pattern features of all observed single-behavior spike patterns in hippocampal neuron types both qualitatively and quantitatively. Importantly, our systematic and more detailed spike pattern classification identifies distinct hippocampal spike pattern classes within general firing behaviors (Komendantov et al., in review). For instance, among the adapting spike patterns, our approach distinguishes between the patterns that reach a specific steady state such as non-adapting or silence (ASP.NASP and ASP.SLN classes, respectively) and those with experimental recordings that only allow determination of the transient state (ASP. class). Our models effectively reproduced the features of these classes (Figure 3: ASP., ASP.SLN, and ASP.NASP).

All models shown in Figure 3 are single compartment IMs except for PSTUT and TSTUT.NASP, which were reproduced by coupling two homogenous compartments. Although stuttering behavior can be modeled in a single compartment IM, multiple compartments (coupling effects) were required to accurately capture various subclasses of stuttering behavior such as TSTUT.NASP and TSTUT.ASP. However, the number of compartments for a multi-compartment IM is determined based on the neuronal morphology (see section Spiking Model) and CA1 O-LMR and CA1 Oriens-Bistratified neurons have both their soma and dendrites in the oriens layer. Thus, MC models were created by symmetrically coupling two identical compartments, unlike the MC IMs with morphologically defined layouts (section Constrained Multi-Compartment Models). These two-compartment IMs were able to capture the classes PSTUT and TSTUT.NASP by integrating coupling effects into the IM dynamics. For the EA search, this simply means inclusion of an additional parameter (coupling constant). In many cases, the EA population converged in less than 500 generations, but certain classes required more generations (Figure 4). As mentioned in section Error Function, we reject the best solution found from a single EA run, if its spike pattern features do not meet the target class criteria (see section Diversity in the Intrinsic Properties across Neuron Types for discussion on the number of accepted models for different classes). The IM parameters of the nine models from Figure 3 are given in Table 1.

Quantitative Comparison of Spike Pattern Features

Our approach can reliably fit a model's responses to multiple experimental voltage traces. As an illustration, the model of a CA1 OR-LM neuron type (a variant of the O-LM interneuron superfamily with dendrites in oriens and axons in both radiatum and lacunosum-moleculare) was created by fitting its responses to four distinct experimental voltage traces recorded for different current stimulation strengths (Figure 5A). The model reproduces features of spike pattern and subthreshold voltage traces that are quantitatively comparable to the experimental traces (Table 2). The model spike pattern features are reported for the input currents that were selected by the EA (see Materials and Methods). In addition, only the minimum set of features required to fully capture the temporal properties of spike patterns were included in the error function. For instance, single spike traces do not require pss as an objective feature, when fsl and nspikes are included. By allowing a narrow range for input current, the EA was able to reliably fit the model responses to multiple voltage traces. Although the voltage sag is not as clearly visible as in the hyperpolarized experimental trace, the corresponding model response nevertheless has a post-inhibitory rebound potential with a 7 mV amplitude. It should be noted that for multiple voltage trace fitting, we only considered traces that were recorded under the same experimental conditions (except for the strength of current stimulation), such as animal species (rat vs. mouse), electrode type (patch vs. sharp), and temperature (room vs. body). As mentioned in section Models of Distinct Single Behavior Types, our approach does not only differentiate between different classes of frequency adapting spike patterns, but also reproduces quantitatively comparable parameters of sfa (Figure 5B). See supplementary material (section Materials and Methods) for quantitative comparison between experimental and model traces for all nine classes from Figure 3.

Furthermore, our approach does not simply identify a single optimal point in the IM parameter space, but instead identifies several possibilities that correspond to the known behaviors of a certain neuron type. The size of such region of possibilities in the parameter space depends on the target behavior class to which the model is constrained as well as the amount of experimental data available for each neuron type. For instance, NASP behavior roughly correspond to the range (0.01, 0.1) for the IM parameter “a,” whereas ASP., especially a strongly adapting behavior, significantly reduces the possibilities to the range (0, 0.005) for “a.” Similarly, if multiple experimental voltage traces were recorded for different input current strengths, the possibilities in the IM parameter space are reduced.

The variability in the quantitative features among all accepted models is given in Figure 5C. The experimentally observed feature typically lies within the range of features observed in the corresponding models, with few exceptions. For instance, while the IMs for the CA3 Giant neuron type exhibited a range of sfa slopes that encompassed the experimentally observed sfa slope (Figure 5C bottom), its fsl lies outside the model range (Figure 5C top). These models were nonetheless accepted because they all satisfy the criteria for the target class (ASP.).

The best model for CA3 Basket (ASP.SLN) showed the highest error in the sfa slope among all the adapting classes (Figure 5C bottom). Yet, the accepted models for this neuron type not only exhibited the desired class (Figure 3, ASP.SLN), but also captured fast-spiking behavior, which plays an important role in network synchronization (Traub et al., 1996; Cardin et al., 2009). It should be emphasized that there is no guarantee that a model fit to a single experimentally recorded fast-spiking trace is indeed a fast-spiking model. For example, if the model of a CA3 Basket fast-spiking neuron type was created by simply fitting to the only available trace (Figure 3, ASP.SLN), it might still exhibit non-fast spiking behavior for a lower current input. To avoid this discrepancy, in fast-spiking models we enforced a minimum instantaneous frequency of 25 Hz for a step current close to the rheobase (Figure 6). In the end, a neuron type is represented by a set of heterogeneous models with similar behavioral features.

Constrained Multi-Compartment Models

In addition to the simple point-neuron models described in previous sections, we create MC Izhikevich models with heterogeneous compartments for all neuron types with dendrites spanning multiple hippocampal layers. These models capture the differences in the active and passive properties between soma and dendrites as well as coupling mechanisms that allow biologically realistic signal transmission between compartments. However, our MC models do not have branching dendritic arbors, and only consists of up to four compartments. This is because each additional compartment adds 10 new parameters for optimization, and tuning hundreds of compartments for each neuron type is an unrealistic goal. We assume that layer-level segregation of synaptic inputs is sufficient to significantly increase the computational power of the models in a network.

As an illustration, we present a four compartment model of CA2 pyramidal neuron type (Figure 7). The somatic compartment (SP) reproduced features of experimentally recorded voltage trace [see Figure 3(2B)] both qualitatively and quantitatively (Figure 7A). Furthermore, we enforced four additional constraints for MC models as detailed in section Error Function. Decoupled dendritic compartments are less excitable than the somatic compartment (Figure 7B) and have higher input resistances (Figure 7C). It should be noted that the dendritic input resistances and excitabilities were only enforced qualitatively (the quality of being higher or lower) relative to the somatic compartment. Neither the absolute values nor the magnitudes of the differences were enforced in the error function (see Equations 6 and 7 in section Error Function).

In addition, the dendritic compartments allow forward propagation of spikes to the adjacent compartments. A few hundred excitatory synapses were simultaneously stimulated in order to initiate a spike at a dendritic compartment, and forward spike propagation (in the direction toward soma) was verified at the adjacent compartment (Figure 7D). Interestingly, the SLM compartment required an additional depolarizing current of 1,200 pA in order to initiate a spike, consistent with experimental observations (Jarsky et al., 2005). Although we enforced a spike propagation rate of 1 for isolated spikes initiated at a dendritic compartment (see section Error Function), we noticed that the rate was less than 1 for high frequency dendritic spike trains. Finally, the amplitude of unitary EPSP measured at the somatic compartment was constrained to be in the biologically realistic range of (0.1, 0.9) mV (Figure 7E).

Even though not directly enforced in the error function, our MC models qualitatively exhibited the known directional voltage attenuation properties of hippocampal neurons: voltage attenuation from a dendritic location to the soma is much higher than in the opposite direction (Mainen et al., 1996; Carnevale et al., 1997; Chitwood et al., 1999; Golding et al., 2005). This behavior was observed in the models because of the higher input resistances of the dendritic compartments and the asymmetric coupling between the compartments. The EA always selected weaker coupling toward the soma than away from it. Optimization of a 4-compartment model required ~20 h of total execution time on the GPU. This is roughly a 15X speedup from CPU execution. The 39 parameters of the model from Figure 7 are given in Table 3.

Variabilities in the Intrinsic Properties Within a Neuron Type

Our models of hippocampal neuron types were constrained using voltage traces digitized from figures in the published literature. It is thus natural to ask: how faithfully does a representative recording from a single neuron, which the authors chose to include in an article, reveal the real intrinsic property of that neuron type? Under the same experimental conditions, a different neuron of the same type might behave slightly differently due to experimental noise and biological variability. A truly accurate model of a neuron type should take into account (and even represent) such intra-neuron type behavioral variabilities. For instance, if the experimental data consist of spike pattern traces recorded from several neurons of the same type under the same experimental conditions, the error function could use statistical measures such as z-score to capture the variability in the models (Druckmann et al., 2007; Markram et al., 2015). However, we ultimately strive to create spiking models for over a hundred hippocampal neuron types based on the available experimental data as gleaned in Hippocampome.org. Except for a few well-studied neuron types such as CA1 pyramidal neurons, DG granule neurons, and a handful of GABAergic interneurons, the vast majority of neuron types identified in the hippocampal formation to date lack adequate data to represent their behavioral variabilities. While such a paucity of empirical evidence might pose the risk of overfitting the model to experimental noise, two key aspects of our approach synergistically reduce that risk: (i) inclusion of qualitative class criteria in the error function and (ii) parameter space exploration using the EA.

Firstly, we dynamically weigh the feature errors during the EA search with weights determined by comparing the model's spike pattern class to the experimental target (see section Materials and Methods). This ensures that several near-optimal points in the error landscape represent the appropriate class, even though the exact feature errors might be higher than the globally optimal point (Figure 8A). Without such a weight-assignment scheme, a near-optimal point might not necessarily represent the target class, because both the feature that defines a boundary between classes (e.g., fsl between ASP. and D.ASP.) and the feature that does not (e.g., pss between ASP. and D.ASP.) would equally contribute to the error. Thus, explicitly integrating qualitative definitions in the error function sharply distinguishes the near-optimal points that satisfy class criteria from the ones that do not (Figure 8B). This increases the EA's chances of finding the models that reproduce the target class.

Secondly, rather than just exploiting the search space to identify a single optimal point that precisely reproduces the spike times, our approach explores the search space and identifies numerous points that elicit a similar behavior. The similarity is governed by the qualitative class definitions, which are inherent to the error landscape as described before. The EA exploration was boosted by a high-rate reset mutation along with a two-point crossover. The downside of such a configuration is the reduced EA reliability in finding acceptable models in certain cases. For instance, only 651 out of 1,000 trials found best models that satisfy the target class criteria (ASP.) for the CA3 Giant neuron type (Figure 8B). It is possible to increase this EA reliability by using a step mutation with a lower rate; however, this will be achieved at the cost of global exploration, ultimately resulting in reduced heterogeneity in the accepted models. In the end, a subset of all best models exhibiting quantitative features with a certain degree of variability is chosen to represent a neuron type (section Quantitative Comparison of Spike Pattern Features). Those features strictly adhere to the criteria for the qualitative class of the spike pattern recorded from that neuron type.

There is no guarantee to avoid over-fitting the best model from a single EA trial to the experimental noise. However, by reducing the acceleration of evolution within the bounds of target class (dynamic feature weight assignment), and identifying several near-optimal points within these bounds (parameter space exploration), we reduce the risk of the best models from all EA trials converging to a single globally optimal point, which might or might not represent a noisy feature. It is worth remembering that a feature threshold of a spike pattern class was statistically inferred from the distribution of that feature from all neuron types (Komendantov et al., in review).

The accepted models for a single behavior showed notable variation in their parameters, except for TSTUT.NASP (Figure 9A). Such a variation was most prominent for parameters “a,” “b,” and “d,” but only the dimensions “a” and “b” are shown in Figure 9A for the nine single-behavior types. Thus, a wide range of parameters yielded similar behaviors, demonstrating the robustness of our EA in exploring the parameter space. This is also consistent with the notion that a given neuron behavior may result from multiple distinct combinations of ion-channel conductance densities (Marder and Prinz, 2002).

Diversity in the Intrinsic Properties Across Neuron Types

In section Variabilities in the Intrinsic Properties within a Neuron Type, we discussed within-neuron type variabilities, where several IMs for a single neuron type reveal slightly different quantitative features for similar input currents. In addition to this, feature diversities across different neuron types, both qualitative and quantitative, also likely play a major role in the emergent properties in a network.

Even different neuron types that exhibit similar qualitative behavior might reveal substantial diversity in their quantitative features. For instance, the neuron types that include transient ASP. in their behavior are CA3 Giant (ASP.), CA3 Basket CCK (ASP.NASP), CA3 Basket (ASP.SLN), and CA2 Pyramidal (D.ASP.). The magnitudes of sfa experimentally observed from these neuron types were ~0.1 for 100 pA, ~0.06 for 400 pA, ~0.02 for 400 pA, and ~0.2 for 400 pA, respectively (Figure 5C). The IMs constrained using these features reveal considerable diversity in the magnitudes of sfa among these four models when the input current is gradually increased, as evidenced by plotting sfa against a range of input currents “I” (Figure 9B). For all cases, sfa decreased exponentially with linear increases of “I.” However, these models showed notable differences in their excitabilities and their sfa ranges. Most of the variance in sfa slopes for the fast spiking CA3 Basket PV+ model could be explained by a narrow range of inputs (325–425 pA). This window is much larger for the regular spiking CA2 pyramidal model. Figure 9B illustrates the diversity of input-dependent sfa ranges among these four models. Although experimental data are too sparse to validate such ranges of sfas in the models, the above results demonstrate that our approach can create models with remarkable quantitative diversities, even with limited amounts of data.

The diverse single behavior classes were most separated along the dimensions “a” and “b” of the parameter space (Figure 9A). The four behaviors that include ASP. were restricted to very small values of “a” (<0.01), whereas the NASP models converged to a broad range (0.02, 0.1). As mentioned in section Spiking Model, “a” is the time constant for recovery variable “U,” and lower values for “a” results in stronger sfa. The region “b” > “a” correspond to Andronov-Hopf bifurcation (Izhikevich, 2003) and all the fast-spiking ASP.SLN models were identified in the range (9, 90) for “b”. Although only 25 best models from 1,000 EA trials satisfied the criteria for ASP.SLN and fast-spiking, these models encompassed a broad range for “b” (Figure 9A top).

The optimal region for each class is shown in Figure 10. There is a significant overlap between the regions for the classes ASP. and ASP.NASP (see also Figure 9A). This is because the difference between these two classes often depends on the input conditions rather than the nine parameters of the model. It is worth mentioning again that in the ASP. class only the transient element is present in the spike pattern. Given a longer duration of input, this pattern will most likely show a steady-state of NASP. The classes NASP and D.NASP encompassed larger regions in the parameter space (Figures 9A, 10). In the case of D.NASP, this is likely due to the fact that the experimentally injected input current was unknown (see Supplementary Table A1), and the EA identified several possibilities for similar behavior under a broad range of input currents.

Discussion

A major motivation behind the current work is the intent to create large-scale network models using IMs with both biologically realistic within-neuron type behavioral variabilities and experimentally validated between-neuron type diversity. Our compact model representations of diverse neuron behaviors allow the implementation of hippocampal circuit simulations in a computationally efficient manner. More importantly, our results offer a sampling range for a neuron group in a network model (Figures 9A, 10). Several studies have shown that neurons have intrinsic plasticity and undergo homeostatic regulatory mechanisms, which modify their non-synaptic ion-channels such as sodium and delayed-rectifier potassium channels, in order to maintain a certain target activity level in the network (Desai et al., 1999; Aizenman and Linden, 2000; Desai, 2003). This implies similar intrinsic properties or behaviors can arise from various combinations of ion channel conductance densities (Foster and Ungar, 1993; Marder and Goaillard, 2006; Schulz et al., 2006). In the mathematically abstracted IM, this is equivalent to various combinations of parameter interactions. Although it might be difficult to describe such interactions mathematically, a robust EA can identify several optimal points in the multi-dimensional search space that correspond to the known behaviors of a neuron type (Figure 8B). Thus, our method represents a neuron type as possibilities in the model parameter space (Figure 9A). Such a representation is crucial for a thorough and systematic investigation of the contributions of neuronal intrinsic properties to network behavior and function. Our multi-compartment models extend this platform to investigate the effects of dendritic filtering on the emergent network properties, while still being reasonably compact.

Furthermore, the diversity captured in our models may help experimentalists identify and distinguish real neurons in finer electrophysiological terms. Our models of sfa suggest the existence of different critical input windows for different neuron types as explained in section Diversity in the Intrinsic Properties Across Neuron Types (Figure 9B). Thus, a neuron exhibiting sfa could be characterized by the range of sfa slopes and its critical input window. The sfa makes a neuron act as high-pass filter (Benda and Herz, 2003) and plays a role in emergent network synchronization (Ermentrout et al., 2001). Two neurons with different input-dependent sfa ranges will likely have different filtering properties, and consequently, may contribute to network synchronization in different ways.

An advantage of using a Pareto-optimal front approach for model optimization (Druckmann et al., 2007) is that it avoids the need to weigh different feature errors. However, the performance of such multi-objective optimization techniques is affected by the number of objectives (Khare et al., 2003) and exponentially increasing population sizes are required to represent high-dimensional Pareto-optimal fronts (Deb, 2014). A single interrupted spike pattern trace (e.g., TSTUT.NASP in Figure 3) presents at least eight objectives for optimization (see Supplementary Table A2, CA1 O-LMR). Moreover, some of our simple models are constrained using several spike pattern traces (Figure 5A). With population sizes as small as 120, our approach can efficiently optimize model parameters for several objectives.

On the other hand, the approach created by Rössert et al. (2016) requires a data-driven microcircuit model constructed from morphologically detailed neuron models (such as the one in Markram et al., 2015) as a reference, and such reference models are computationally very expensive. Compared to this approach, our simple models might be less constrained in some cases, but they significantly reduce the open parameter space size to create biologically accurate circuit models. Furthermore, our simplified multi-compartment models intrinsically capture the dendritic voltage attenuation properties without a need for synaptic correction (section Constrained Multi-Compartment Models).

The precise shape of the spike was not captured in some of our models (for example, D.NASP in Figure 3). This could be attributed to the quadratic voltage dependence in the IM voltage equation. The AdEx model (Brette and Gerstner, 2005), which replaces the quadratic term in the IM with an exponential term for the voltage dependence, has been shown to reproduce more realistic spike shapes (Badel et al., 2008). However, our selection criteria for the models were entirely based on the temporal features of the overall spike pattern, and do not include characteristics of individual spikes. More important for the information processing in a neural network are the excitability of neurons, the precise timing of spikes, and the properties of connections. The shape of the spike is unlikely to play an equally prominent role in network dynamics. In fact, in the nine-parameter IM formalism “k” and “Vt” collectively determine the shape of the spike. It is thus possible to obtain realistic spike shape by restricting the ranges for these parameters (e.g., Figure 3, ASP.SLN). However, we did not explore these parameter ranges and interactions for all the cases for the reasons mentioned above.

Although only IMs have been presented in this article, our framework can be easily enhanced to include any phenomenologically rich model of spiking behavior. The only part of this framework that is specific to the IM is the EA configuration presented in section EA Configuration. This configuration was identified partly based on the topographical features of error landscape created by the IM parameters. Once an appropriate EA configuration is identified, our error function and spike pattern classification procedures are readily applicable to any alternative model.

In the future, we will enhance our framework to model multi-behavior neuron types. At least 15 morphological neuron types in the Hippocampome.org exhibit sharply distinguishable qualitative features under different experimental conditions. One of the commonly occurring multi-behavior types in the hippocampus is stuttering and spiking observed in a single neuron for different current stimulation strengths. For example, a CA1 Bistratified neuron exhibited stuttering and regular spiking behaviors for 400 and 600 pA, respectively (Pawelzik et al., 2002). Similarly, a CA1 Neurogliaform projecting neuron exhibited this multi-behavior for 500 and 700 pA (Price et al., 2005). Our preliminary work with multi-behavior types revealed vast possibilities for modeling such behaviors using IM, which could also provide insights into the existence of electrophysiological subtypes for a given morphological type.

The categorization of neuron type behaviors as either single-behavior or multi-behavior is solely based on the currently available experimental data. Consequently, it is possible that additional qualitative behaviors will be observed in future experiments from neuron type currently considered to display a single-behavior based on available data. An advantage of our modeling approach is that it identifies many possibilities for the known behaviors of a neuron type in the IM parameter space. Furthermore, the flexibility of our framework allows easier addition of newly observed behaviors from a neuron type to improve the accuracy of its model representations. Eventually, we plan to create models for over a hundred hippocampal neuron types and to make them freely available at Hippocampome.org. Nevertheless, although the modeling framework and the results presented in this article pertain to the hippocampus, our approach could be easily adapted to other brain regions.

Author Contributions

SV, AK, and GA created the overall approach and all the models. SL and JK extended the CARLsim framework to include simulation of multicompartment Izhikevich models on the GPU. ES and KD developed the interface between ECJ and CARLsim.

Funding

This work is supported by grants R01NS39600 (NIH) and IIS-1302256 (NSF). Publication of this article was funded in part by the George Mason University Libraries Open Access Publishing Fund.

Conflict of Interest Statement

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

Acknowledgments

We thank Charise White, Diek Wheeler, Ernest Barreto, and Sumit Nanda for their contribution and valuable comments on improving this article.

Supplementary Material

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

Footnotes

References