Neuroscience data are being synthesized into increasingly complex computational models

Since the early explorations of neurons and their circuits by Golgi [2] and Cajal [3,4], humans have been fascinated by neuron diversity and complexity. In the 1950’s, Hodgkin and Huxley [5] were the first to synthesize electrophysiology data into a mathematical model that accurately predicted the propagation of axonal action potentials. Their approach has been extended to include a wider variety of ion channels and neuronal morphologies [6,7]. More recently, data acquisition has escalated so that large data sets have been collected that describe neurons and the networks the neurons form. In parallel, exponentially increasing computing power has allowed the construction of large, biophysically-realistic network models of connected cells [8–12].

The ability to rapidly find, evaluate, select, and reuse earlier models is becoming more important

Despite the availability of thousands of previously-created neuron models [13,14], it is still relatively difficult and tedious to find, evaluate, and select previously created models or their components like channels and synapses for reuse in a new project. Lack of consistent and easily accessible information about a model’s electrophysiology, morphology, and computational complexity makes it difficult to rapidly evaluate whether a neuron model is fit for a particular modeling purpose, resulting in the problem of model selection. Once a previously published model is found to be suitable, the use of heterogeneous programming and simulator languages makes it difficult to easily reuse models or model components, resulting in the problem of model reuse. These two problems hinder progress within the field. If current trends towards system-level modeling continue to progress toward realistic models of entire brains, the ability to rapidly leverage previously developed models will be paramount.

Online platforms for model sharing such as ModelDB [13,15] and Open Source Brain [16–18] provide rich model search and inspection capabilities to help address the model evaluation problem. Meanwhile, efforts to standardize computational neuroscience model descriptions like PyNN [19] and NeuroML [1,20] help address the problem of model reuse.

ModelDB is an online repository with an extensive database of published computational neuroscience models. While many models are implemented using the NEURON simulator [21,22], other simulators and programming languages are represented as well [14]. ModelDB provides extensive model search and browse capabilities and, for models that are implemented in NEURON, allows users to view some aspects of the model structure like cell morphology and ion channel/synapse files. Similarly, Open Source Brain is an online collaborative environment for the development of multiscale neuroscience models. The platform leverages the structure and relationships of NeuroML models, version control provided by GitHub, and simulation visualization using the web-based Geppetto platform [23]. On Open Source Brain, model files, 3D structure, and simulation results can be inspected, and models or model components can be downloaded. Many of the models available in NeuroML-DB were originally translated to NeuroML at Open Source Brain. Similarly, an emerging online platform like Arkheia [24], where users can examine models, their parameters, and simulation results could also make it easier to evaluate models. A more narrowly focused resource, the ICGenealogy project [25,26] applies consistent and uniform stimulation protocols to ion channel models implemented in the NEURON NMODL language [27,28]. The uniform protocols enable the comparison of channel models to each other, providing insight into channel model taxonomy and publication genealogy. Additionally, a new model or biological channel voltage clamp data can be uploaded to the website [26,29], which identifies similar channel models in the database allowing rapid identification of channel dynamics and a list of seed models which could be used for further fitting. The Allen Brain Atlas Cell Types Database shares models developed based on experimental data from neurons in visual cortex and the corresponding electrophysiology and morphology data [30,31]; however, this resource is limited to data and models from the Allen Institute.

NeuroML model description standardization effort promotes model reuse

The above resources make it easier to find and evaluate previously published models. Meanwhile, the standardization initiatives make it easier to reuse whole or parts of existing models. Through PyNN [19], users can specify network models composed of abstract or single compartment conductance based cell models using an expressive connectivity syntax in Python. Then by adjusting a single line of code, the model can be executed on any compatible software or hardware simulator. On the other hand, NeuroML and associated tools [32] can be used to specify network models composed of multi-compartment cells and biophysically realistic channels and synapses, as well as more abstract model formulations. Though earlier NeuroML models could only be specified using the human-readable XML format, the latest version supports compact storage using the HDF5 format [18], which facilitates the development of large, systems-level models. The modular nature of NeuroML makes it easy to extract subcomponent channel, cell, or synapse models for reuse. Additionally, these extracted components can be converted to a variety of simulator formats using automated tools [32], allowing rapid development of novel models by composition of subcomponents of earlier models. The NeuroML documentation includes examples of how to extend these existing models to create novel cell and network models [33]. The ability to convert NeuroML to a variety of simulator formats also makes it possible to compare simulators using the same models [1].

NeuroML database catalogs over 1,500 NeuroML models and facilitates model evaluation

Ideally, to rapidly develop a novel model, a user would be able to use an online resource to simultaneously evaluate many models and then easily select models or their components for reuse. While the model repositories described above help with locating and evaluating models, and modular and simulator-agnostic languages help with the reuse of model components, no online resource exists that combines rapid search, deep model inspection, and evaluation of features with the modular architecture of NeuroML and exposes the features via an automated interface. To make progress towards this vision, we developed the NeuroML Database and its web-based interface[34,35] (https://neuroml-db.org), cataloged published models translated to NeuroML, added model search and extensive characterization features, and implemented a web accessible web-based Application Programming Interface (API) to make it easier for researchers to evaluate, select, and reuse these models.

Over 1,500 channel, synapse, cell, and network models were added to NeuroML-DB

We identified 45 publications whose models had been translated to NeuroML (see Table 1) and developed a custom, semi-automated procedure to include these models in the NeuroML database. After parsing the individual model components such as the cells, ion channels, and synapses, the database includes a total of 1,222 cell, 183 channel, 141 synapse, 27 ion concentration, and 11 network models (1,584 overall). The translation of some of these models to NeuroML was part of other efforts related to the NeuroML initiative [18,39,40]. Additional models can be incorporated into the database rapidly and easily using this same procedure, and we encourage readers to refer us to any other published models that have been translated to NeuroML.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. The publications with models that have been translated into NeuroML and indexed by NeuroML-DB.

Model counts include the top-level network or cell models as well as any distinct subcomponent synapse and channel models.

https://doi.org/10.1371/journal.pcbi.1010941.t001

Channel electrophysiology was characterized using ICGenealogy project voltage clamp protocols

The ICGenealogy project has developed a uniform set of voltage clamp protocols [26] that can be used to stimulate ion channel models and compare their responses (and compute their similarity indices) in an automated fashion. These protocols can be used to differentiate a variety of ion channel behaviors. While the ion channel equations stored in the NeuroML Database are stored in a structured, cross-platform XML format, it can be difficult to rapidly assess the broad ion channel behavior properties from the model equations and their parameters alone. To facilitate this assessment, we subjected the ion channels stored in the database to a subset of ICGenealogy protocols and made the model response waveforms available online in the form of interactive plots (see Fig 1). These plots allow the user to rapidly gauge the general dynamics of an ion channel model.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. NeuroML-DB displays voltage clamp characterizations of all indexed channel models.

This screenshot from NeuroML-DB depicts the characterization of an example potassium channel model [62]. In addition to the activation protocol, users can view responses to deactivation and inactivation protocols using the drop-down list widget. The slider and plot zoom controls can be used to inspect individual traces. The above plots can be viewed online [63] and also accessed programmatically using the NeuroML-DB API.

https://doi.org/10.1371/journal.pcbi.1010941.g001

For each ion channel model, the plots show the channel model input voltage and output conductance and current levels recorded over the course of simulations. A choice of protocols (Activation, Deactivation, and Inactivation) can be selected from a drop-down list and an interactive slider and plot zoom controls can be used to inspect individual traces. The waveforms are also machine-accessible via the NeuroML database API. The reversal potentials and input voltage waveforms were identical to those used in ICGenealogy.

Cell morphology was characterized using L-Measure and visualized using BlenderNEURON

Cell models vary in their morphological detail. Some are single-compartment “point neurons”, while others include detailed, reconstructed morphologies. Similarly to channel models, while cell models are defined in a machine-readable NeuroML format, it can be difficult to visualize the 3D shape of a cell model just from the coordinates of its neurite compartments. On NeuroML-DB, we made it easy to immediately view an animated 3D rendering of each cell model with more than one compartment (see Fig 2B).

To visualize it, each cell model’s morphology was first aligned along its first principal component and a 360° view of the model is shown from a slightly tilted, horizontal orbit. Additionally, a square current is injected into the model cell’s soma, and the resulting action potential propagation (if any) is shown as illuminated compartments. The animations are in the form of.GIF files, which do not require any additional software or plugins, and can be embedded, in animated form, into presentations and shared on social media platforms. To create the animations, we developed an open-source tool, BlenderNEURON [41], which creates a Python interface between the NEURON simulator and the open-source 3D modeling software Blender.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Screenshots of characterizations of a Descending Axon Cell [8] model available through the NeuroML-DB interface.

A) An electrophysiology tab shows plots of current clamp responses, closest electrophysiology clusters, differentiating property values, and other cell models with similar behaviors. B) A morphology tab shows animated 3D visualizations of multi-compartment cell model geometry and electrical behavior, as well as a list of morphology metrics as computed by the L-Measure software package. C) A computational complexity tab shows model equation counts and cell model simulation speed comparisons relative to the Hodgkin-Huxley model. An example of the above plots can be viewed online [64].

https://doi.org/10.1371/journal.pcbi.1010941.g002

In addition to visual characterization, cell model morphology was characterized by computing the morphology metrics provided by the widely used L-Measure tool [38]. The same metrics that are visible for each reconstructed cell on NeuroMorpho.org [42] are also shown on the database website (Fig 2B bottom). The full set of computed L-Measure metrics is accessible via the NeuroML database API.

Cell model computational requirements were characterized using numerical benchmarks

When evaluating cell models for a particular purpose, one often overlooked aspect is computational complexity–for practical purposes, how “fast” does the model run? The computational complexity of a model will affect the implementation choices needed for the optimization algorithm used for parameter fitting, as well as the overall pace of model development. In general, if all other model aspects are approximately the same, a model with a lower computational complexity likely is preferable to one with higher complexity. Computational efficiency might be especially important in model reuse scenarios where a new candidate model (for example a single component in a larger system) requires a much smaller time step than all the other model components, making it the rate-limiting step in simulating that system.

In modeling publications, the choice of time step is often reported, but rarely rigorously justified. For example, often it is not trivial to know how sensitive a particular model is to deviations from the published time step. To facilitate such assessment and allow for automated evaluation of model simulation efficiency, we developed a simple approach for characterizing computational requirements for variable time step integration methods. Note that variable time step methods are typically used for neuron models, since these methods adjust the size of the time interval depending on how quickly the cell state changes (e.g. smaller step during an action potential, larger step when the cell has reached a steady state).

Cell model complexity using variable time step size integration method can be measured using baseline steps/sec and mean steps/AP at a target firing rate

The popular NEURON simulator enables simulations to use a variable time step integration method (CVODE) [65], which adjusts the size of the time step to maintain a constant local error tolerance. Because the size of integration steps varies during the simulation, we first observed the number of steps the NEURON simulator used during each millisecond to compute the response during current injections that produce an action potential (spike) (Fig 3). We noted that each model tended to have a baseline number of steps, which increased significantly during action potentials. When higher intensity current injections produced more action potentials, the total number of steps to compute 1 s long simulation increased roughly linearly with the number of action potentials, across a variety of cell models. Thus, we modeled the number of steps required to compute 1 s of a cell model’s simulation as a linear function in the form of: . For each model, the slope and intercept terms were fitted using a linear regression of total steps vs. number of action potentials produced in response to a series of square current injections (see Methods). Each model’s absolute 1 s complexity was computed using an assumed 10 Hz (10 APs/sec) rate and scaled to the absolute 1 s complexity of the Hodgkin-Huxley model. See Methods for more details.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Number of integration steps required when using a variable time step integration method (NEURON CVODE) depends on the baseline step rate and additional steps per action potential (AP) for an example model [11].

The plot shows membrane potential (orange) and the number of integration steps per ms computed by the simulator (blue) in response to a square current injection with onset at 50 ms. This “baseline + APs” pattern was similar across other types of input current waveforms.

https://doi.org/10.1371/journal.pcbi.1010941.g003

After assessing each model’s relative computational requirements, we compared the computational complexity of the Izhikevich [48] and Generalized Leaky Integrate and Fire (GLIF) [54] classes of models to the Hodgkin-Huxley model (Table 2). These classes of models are some of the simplest spiking models and are often used in large network simulations. A previous, frequently cited, computational complexity analysis (Fig 2 of [66]) based on floating point operations (FLOPS) suggested a 92-fold speed difference between the Izhikevich and the Hodgkin-Huxley model (13 vs. 1,200 FLOPS). Using the variable time step method, the GLIF models were between 5.5 and 20 times faster than the Hodgkin-Huxley model, while the Izhikevich models were between 1.3 to 3.5 times faster than the Hodgkin-Huxley model.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Computational complexities (“speeds”) of select cell models in the NeuroML database.

See Methods for details. Complexities are relative to the Hodgkin & Huxley model, for which complexity is defined to be 1 HH. Simulations were performed using NEURON CVODE variable step integration method. Top rows: Models with complexities smaller than the Hodgkin & Huxley model (e.g. 0.2 HH means the model was 5 times faster than the Hodgkin & Huxley model). Bottom row: Models more complex than the Hodgkin & Huxley model (e.g. 50 HH means model was 50 times slower than the Hodgkin & Huxley model).

https://doi.org/10.1371/journal.pcbi.1010941.t002

After computing the variable time step computational complexity metrics for each cell model, we made the analyses viewable online (Fig 2C). It should be noted that these performance comparisons were between NeuroML models converted to NEURON using automated conversion tools. These auto-generated models (NEURON mod and hoc files generated from NeuroML) have been robustly tested for correctness against the original model codes using the Open Source Brain Model Validation framework [18]. However, they sometimes are not as computationally efficient as hand written code can be. It is possible that comparing manually optimized model versions would yield different speedup ratios.

Cell model electrophysiology was characterized using Allen Brain Atlas protocols and Human Brain Project electrophysiology properties

When assessing cell model electrophysiology, detailed cell model traces in response to different current clamp stimulation protocols are usually not available in the original publications. One must generally download the full model and perform the current clamp experiments manually. While this is feasible when evaluating a handful of models, it quickly becomes impractical for a large number of models. This limits the number of evaluated models, possibly resulting in a sub-optimal set of initial candidates. To facilitate this task, similarly to channel models, we have made the responses of current clamp responses to a set of standard protocols available online as interactive plots (Fig 2A). We characterized cell model responses using the electrophysiology protocols used by the Allen Cell Type Database [31], which included square, long square, pink noise, ramp, short square, and short square triple protocols (see [67] pages 7 and 15 for protocol details).

Additionally, we computed 38 cell model membrane properties described by Druckmann, et. al. [68], which were used in cell type classification by the Human Brain Project [69]. Broadly, these measures assessed the properties of individual and trains of action potentials. Example action potential properties included amplitude, width, and after hyperpolarization potential, while example action potential train properties included action potential delays, inter-spike interval statistics, and degrees of spike accommodation (for full list of measures and their computation details see Table 1 and Supplementary Methods of [68], respectively). These neuron model characterizations are available as online tables (Fig 2A). Furthermore, to facilitate reuse of these properties, we implemented them as standardized tests within the SciUnit/NeuronUnit framework [70].

Cluster analysis revealed the structure of model neuron electrophysiology property space

After characterizing the electrophysiology properties of cell models in the NeuroML database, we wanted to explore the structure of the space formed by the electrophysiology measures. We were interested in identifying the features that could be used to summarize cell model behavior and identify any high-density clusters (e.g. cell model types) that existed within this space. These features and cluster memberships could then be displayed online, facilitating model comparison and selection tasks.

In previous work, the nomenclature established during the Petilla convention [71] provided a broad classification scheme of interneuron electrical behavior based on the consensus of the convention attendees. Later, Druckmann and others [68], used a set of 38 action potential and spike train measures of rat cortical interneurons to perform automated classification of interneurons into electrical types, which were then utilized in Human Brain Project cortical column simulations [8]. Similarly, a taxonomy of mouse cortical cells performed at the Allen Institute [72], and available to explore online [73], was constructed using automated clustering of single-cell RNA sequencing data.

We ask whether these classifications are recapitulated in neuron models. For all neuron models in NeuroML-DB, we computed the 38 properties described by Druckmann et al. [68]. We reduced the dimensionality of the features by using principal component analysis [74] and used k-means and HDBSCAN clustering algorithms [75,76] to identify high-density regions of cell models within the reduced space. The analysis resulted in three levels of nested model clusters, with four clusters at the top level, two clusters in the second level, and six clusters in the third level (Fig 4).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. The hierarchy of cell model types identified during clustering analysis of NeuroML-DB neuron models according to electrophysiology features.

The hierarchy was determined by the clustering algorithms as described in the methods section. The first two levels (yellow and blue boxes) were determined using the HDBSCAN algorithm [76]. The last level (orange boxes) was determined using K-means [75]. Numbers to the right of boxes indicate the number of cell models in each cluster.

https://doi.org/10.1371/journal.pcbi.1010941.g004

At the first level of this analysis, we identified four clusters of cell models: Multi-Spikers (MS), Rapidly Adapting (RA), Abstract Spikers (AS), and Non-Positive Rheobase (npRB) (see Fig 5 for examples of models belonging to these clusters). Multi-Spikers are a large group of cell models that produce multiple spikes in response to square current injections. The Rapidly Adapting cluster contains neuron models which produce single spikes in response to square current injections. While some action potential properties can be computed for these cell models (e.g. delay to first spike, amplitude, half-width), spike train properties generally cannot be computed. Abstract Spikers are a group of neuron models for which traditional measures of action potentials do not apply, for example the GLIF cell models [54] have a resting potential at 0 mV and do not define action potential amplitude or width (Fig 5). However, the defining properties of the final six clusters in the third level are predominantly related to spike train patterns (Fig 6). For many applications, spike train patterns are of exclusive interest. While not associated with the final clusters in this analysis, the Abstract Spikers would likely belong to one of the six cluster types. Finally, the Non-positive Rheobase cluster contains neuron models that do not have a positive rheobase current and produce spikes without stimulation. Because the majority of the Druckmann properties assume positive rheobase currents, their values could not be computed (see Discussion).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Responses of neuron models closest to the centers of clusters in the first and second levels of the cluster hierarchy.

awMS: Accommodating Wide Multi-Spikers exhibit spike rate accommodation and long-lasting action potentials (NeuroML DB ID: NMLCL000670). RA: Rapidly Adapting models produce single spikes in response to strong square current injections (ID: NMLCL001126). AS: Abstract Spikers do not have physiologically realistic waveforms (e.g. resting voltage at 0 mV or 0 mV action potential amplitudes. ID: NMLCL001491). npRB: Non-positive Rheobase models produce action potentials at rest (ID: NMLCL001588).

https://doi.org/10.1371/journal.pcbi.1010941.g005

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Responses of neuron models closest to centers of clusters for electrical types located in the third level of the cluster hierarchy.

RS: Regular Spikers exhibit mild spike delay, some spike rate accommodation, and an average steady state spiking rate (NeuroML DB ID: NMLCL000314). dRS: Delayed Regular Spikers are similar to Regular Spikers with longer spike onset delay as indicated by the arrow (ID: NMLCL000468). naRS: Non-accommodating Regular Spikers are similar to Regular Spikers but exhibit little accommodation, as seen in the areas marked by bars (ID: NMLCL000829). FS: Fast Spikers exhibit rapid spike onset (see arrow), little accommodation, and a high spike rate (ID: NMLCL001025). aFS: Accommodating Fast Spikers are similar to Fast Spikers but exhibit some accommodation, as seen in the areas marked by bars (ID: NMLCL000728). B: Bursters exhibit rapid spike onset, high accommodation, and a relatively low spike rate (ID: NMLCL000203).

https://doi.org/10.1371/journal.pcbi.1010941.g006

We then performed the same PCA and clustering analysis on the property values of the neuron models that belong to the Multi-Spiker cluster. At this second level of analysis, we identified a large cluster of Regular Multi-Spikers (rMS) and a smaller cluster of Accommodating Wide Multi-Spikers (awMS) (Fig 5). The Accommodating Wide Multi-Spikers are different from the Regular Multi-Spikers in that they display much longer peak-to-trough action potential widths and exhibit spike accommodation.

Finally, we performed an additional cluster identification procedure on the models that belong to the Regular Multi-Spiker cluster. At this third level, we identified six clusters of models: Regular Spikers (RS), Delayed Regular Spikers (dRS), Non-accommodating Regular-Spikers (naRS), Fast Spikers (FS), Accommodating Fast Spikers (aFS), and Bursters (B) (see Fig 6 for examples of each). To provide a general overview of each neuron model’s electrical behavior, we identify which cluster center is closest to the model and display the information for that cluster on the NeuroML database interface (see Fig 2A for an example).

Component loadings identified features responsible for clustering of Regular Multi-Spiker models

The space formed by the first few principal components of the PCA procedure is somewhat abstract and lacks an intuitive mapping to the individual features that are analyzed. To gain a more intuitive understanding of what the components represent for Regular Multi-Spikers, we first selected the first three principal components and identified single features with the highest absolute principal component weights. For the first component, the delay to first action potential had the highest weight, the second component was most heavily weighted by the median inter-spike interval, and the third component was most heavily weighted by the mean steady state accommodation percent.

When the values of these three features are plotted for the models in the cluster of Regular Multi-Spikers, we observe a roughly tetrahedral shape (Fig 7). The Non-Accommodating Regular Spikers and Regular Spikers form the sides of the tetrahedron. The sides of the tetrahedron suggest that the three properties of models belonging to the Regular Spiker sub-clusters exhibit linear relationships. We found that the Delay to First AP was strongly negatively correlated to the Mean Accommodation at Steady State for models within the Delayed Regular Spikers (Fig 7B and 7C left) and Regular Spikers (Fig 7B) clusters (dRS: Pearson r = 0.94, p<0.001; RS: Pearson r = 0.77, p<0.001). Similarly, the Non-accommodating Regular Spikers (Fig 7B and 7C right) exhibit a strong linear relationship between the Delay to First AP and Median Inter-spike Interval (naRS: Pearson r = 0.96, p<0.001). Interestingly, these relationships are weak or not significant within the models belonging to the Fast Spiker, Accommodating Fast Spiker, and Burster clusters.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Tetrahedron of the models assigned to the high-density cluster cores of Regular Multi-Spikers.

A) The level three cell model clusters plotted within the space of the three most informative properties (Delay, Accommodation, Median ISI) form a tetrahedral structure (gray lines). B) Left, Top-Down, and Right views of the tetrahedral structure. Note the strong linear relationships within the dRS and RS (delay vs. accommodation), and naRS clusters (delay vs. median ISI). C) Left: the responses of two extreme models of the dRS cluster (circled in B Left View). Note the relationship between spike delay (indicated by arrows) and accommodation (compare interspike intervals indicated by bars). Right: Extreme models of the naRS cluster (circled in B Right View). Note the relationship between median ISI and delay. For clarity, several scattered cell models that were not in the vicinity of the cluster cores are not shown (they were included in analysis, see Methods. Interactive plots that include all models can be viewed at https://tabsoft.co/32nHnNd).

https://doi.org/10.1371/journal.pcbi.1010941.g007

The NeuroML-DB interface facilitates rapid integration of model characterizations and metadata

As illustrated above, various voltage clamp responses of ion channel models are openly available from NeuroML-DB. Here we demonstrate the utility of NeuroML-DB for facilitating analyses of neuron models across multiple scales. We ask whether the transitions along the tetrahedral structure of the discovered model property space may be explained in terms of model ion channel densities–thus bridging channel mechanisms to neuron electrophysiology. To that end, we first characterized channel model voltage-responses using the ICGenealogy analysis protocols [26] as described in detail in Methods. Our characterization protocol resulted in compact, quantitative PCA-based representations of responses for each individual channel model organized by ionic channel families (Kv, Nav, Cav, Ih, and KCA). Finally, these PCA-based representations were used to cluster channel models into channel model sub-types using agglomerative hierarchical clustering. The resulting 22 channel model sub-types were then used to probe differences between cell model clusters for the six multiple regular spiking clusters (last level of Fig 4).

The 22 resultant channel model sub-types were not easily interpretable due to a lack of correlations to recognizable computed features (e.g., mean delay or interspike intervals). However, channel model metadata is programmatically available via the NeuroML-DB API. Channel model metadata were downloaded, and channel model names were grouped by their associated sub-types. Using these names, channel models could be characterized by common channel model names. Such names reflect the intended purpose of the models, the expertise of the original modelers, and common naming conventions used in the literature.

Consistency in channel model descriptions is observable for the larger channel model sub-type clusters. The three largest clusters belonged to the Kv (N = 22 and N = 19) and Nav (N = 32) ion channel families, while the rest of the cluster sizes ranged from one to seven models in size (Fig 8A). The largest clusters were populated with 1) delayed rectifiers and slow non-inactivating potassium channels (Kv Cluster A, Fig 8B), 2) A-type and slow inactivating potassium channels (Kv Cluster B, Fig 8B), and 3) fast transient inactivating sodium channels.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8. Analysis of properties of Regular Multi-Spikers based on discovered channel model sub-types.

A) Dimensionality reduction and clustering reveal multiple channel model sub-types across different channel model families. B) NCDs of large potassium channel model sub-types (Z-Score<3) plotted against the electrophysiological feature space reveals global patterns for cell model parameterization. (Selected channel family clusters are indicated by gray arrows.) Note the increasing absence along the delay to 1^st AP axis for Kv Cluster A (Left) and the mixed patterns for Kv Cluster B (Right). C) Reference schematic of Regular Multi-Spiker high-density cluster cores in the tetrahedral property space. D) Average NCDs for each Regular Multi-Spiker cluster using all models within the base (Top), middle (Bottom-Left and Bottom-Middle), and apex (Bottom-Right) property regions (with bootstrapped 95% confidence intervals).

https://doi.org/10.1371/journal.pcbi.1010941.g008

Joint large-scale analysis integrated model characterizations across multiple scales

Channel models are used to model active conductances, which are spread across the channel membrane. The conductance densities, which may be non-uniform across the cell, can be extracted from the NeuroML model’s description. Understanding the emergence of the global structure of cell model electrophysiology requires mapping the contributions of model mechanisms at these lower scales to the computed electrophysiological features at the scale of the cell model. To gain such an understanding, we integrated the channel densities for co-clustered channel models into a single normalized channel density (NCD) representing the summed presence of the channel model cluster within the cell model as described in the Methods section. This was done for each channel model cluster and enabled joint visualization of the parameter space of channel models over the property space of cell model behaviors (Fig 8B).

The tetrahedral structure of the electrophysiological property space suggests a straight-forward way of relating the underlying mechanisms shared within a cell model cluster to this global structure. By examining how changes in average NCD correlate to movements within different property regions of the tetrahedron (Fig 7C), we observed distinct trends for the two largest potassium channel clusters. These channel clusters are known to influence long timescale dynamics. Kv Cluster A consists of voltage-gated channels with non-activating outward currents with slow activation. Kv Cluster B consists of voltage-gated channels with activating outward currents on similar timescales to Kv Cluster A. While Kv Cluster A models only have one non-activating subprocess, Kv Cluster B models have two independent subprocesses that become present at different levels of depolarization. Additionally, the inactivation subprocess of Kv Cluster B models can have time constants upward of 40–50 ms providing a slowly decreasing outward current. The differences between these two channel model clusters results in competing resonance effects that can both decrease and increase cell model spike frequency, respectively.

To what extent are transitions along the tetrahedral structure explainable by the relative contribution of the two potassium channel model sub-types? First, transitions along the base corners indicate trade-offs between accommodation and interspike intervals (Fig 8C). These transitions of the average model in the dominant base clusters (FS, aFS, B) revealed first that increasing accommodation was associated with decreasing mean NCD for Kv Cluster B shifting the NCD ratio of Kv Cluster A to Kv Cluster B to greater than 1—illustrated in moving from the FS corner to the aFS corner. This ratio also appeared to be conserved when examining the shift from the aFS corner to the B corner. The transition was also associated with increasing mean NCD for both Kv channel model sub-types (Fig 8D, top) and increasing median ISI.

Secondly, the transition into the middle property region (increasing delay to first action potential) toward the RS cluster also conserved the mean NCD ratio of Kv Cluster A to Kv Cluster B (greater than 1) but was associated with decreasing mean NCD for both Kv channel model sub-types (Fig 8D, bottom-left). However, the transition into the middle property region between the FS corner to the naRS cluster was associated with a decrease in the mean NCD of Kv Cluster A and an increase in the mean NCD of Kv Cluster B (Fig 8D, bottom-middle). This shifted the NCD ratio to be less than 1. Interestingly, the dRS cluster in the apex of the tetrahedron also exhibits a decreasing mean NCD of Kv Cluster A, an NCD ratio less than 1, and spans the range of accommodation between the two middle clusters. These results suggest that the ratio of channel densities between the two largest potassium channel sub-types influences the degree of accommodation. Finally, a decrease in the overall channel density of delayed rectifiers and slow non-inactivating potassium channels (Kv Cluster A) increases the first spike latency evident in both sides of the tetrahedral property space.

NeuroML models on NeuroML-DB can be easily extended and reused

Models in the NeuroML database are stored and made available for download in the NeuroML format. Software tools that are able to take this format as input can then be used to further process these models. For example, the jNeuroML [78] and pyNeuroML [79] libraries can be used to convert the NeuroML models to simulator formats such as NEURON [65], NetPyNE [80], XPP [81], and MOOSE [82]. Additionally, it is possible to convert spiking neuron network models to PyNN scripts [19], allowing for simulations across NEST [83], NEURON, and Brian [84].

NeuroML-DB provides tested, pre-converted versions of channels (‘.mod’) and cells (‘.hoc’) in formats compatible with NEURON; these can then be used in any other software tools that can take NEURON files as inputs. Furthermore, pyNeuroML has additional features which allow automated analysis of channel model dynamics beyond what is provided by NeuroML-DB. Additionally, model construction tools like neuroConstruct [85] and NetPyNE [80] allow the composition of larger models from NeuroML components. NeuroML files can be visualized using the Open Source Brain Model Explorer [16–18] powered by the Geppetto [23] platform.

Novel methods quantify computational requirements and reveal speed differences among model classes

We developed a set of measures of model computational requirements, where complexity is approximated relative to the well-known Hodgkin-Huxley model. Although this approach is not machine invariant, it provides a useful estimate that it is intuitive for the user. For the variable step integration method, we developed a protocol to identify the baseline number of simulation steps required to compute a unit of simulation time and the number of additional steps required to compute each additional action potential.

While it is possible to estimate model complexity by examining a model’s time constants or performing an analysis of required FLOPs [66], the measures we developed here are strictly empirical and amenable to automation because they do not require an examination of each model’s equations and parameters. Using this automated method, we were able to assess the computational complexity of 1,000+ neuron models and identify the Izhikevich and GLIF models to be the fastest classes of models in the database.

The relative complexity measure we developed has some limitations. Perhaps different clusters of cell models should have different, more tailored input waveforms. If shown to be the case, we could add different protocols in future versions of the database. Similarly, the variable step measure assumes a 10 Hz spiking rate to compute the relative model complexity. It’s possible that this spiking rate is not appropriate for some cells or provides an unfair assessment. Using an expected firing rate for each cell could be a more accurate way to assess the model’s complexity. Finally, the variable-step complexity measure ignores any computations and delays due to synaptic activity and inputs, which, depending on network design, may dominate the total network computational complexity.

Automated assessment identifies key model electrophysiology properties and behavior clusters

Using raw traces of cell model responses to a standardized set of protocols, it is possible to identify a small number of features that capture a large portion of variability in model electrical behavior. If high density regions are present within the space defined by these leading features, knowing to which cluster a given model is closest also provides information about the model’s behavior. To identify such features and clusters, we performed a nested PCA and clustering analysis of the cell model electrophysiology features. This analysis was consistent with previous findings that suggested that cell electrical behavior is not uniformly distributed [68], that it exhibited clusters of models with familiar spiking behaviors such as Regular and Fast Spiking, and reflected some of the important features (e.g. accommodation, delay) identified during the Petilla convention [71]. Importantly, it identified the delay to first spike, median inter-spike interval, and steady state spike accommodation as three individual measures which are highly informative of cell model behavior. By making the values of these three features and the cluster membership of each model available on NeuroML-DB, we have provided an efficient public summary of each cell model’s electrical behavior. The cell model pages also provide a list of cell models with similar spiking properties, which can be used to identify cell models that have similar spiking behavior but have other properties (e.g. different channel types, morphologies, complexity) that make the models more useful for a particular modeling purpose.

Top level neuron model clusters reveal opportunities to improve electrophysiology protocols

The electrophysiology protocols and computed properties used here effectively assessed the electrical behavior of multi-spiker models–neuron models with negative resting potentials that produced multiple spikes in response to 1.5x and 3.0x rheobase square current injections with no spontaneous spiking. However, this protocol and its set of electrophysiology properties did not effectively characterize cell models which did not meet these criteria.

For example, cell models in the Rapidly Adapting cluster do not produce spikes at 1.5x rheobase (1.5x rheobase is sub-threshold stimulation, e.g. RA in Fig 5). For such models, properties that rely on the production of more than one spike (e.g. second spike amplitude, or median inter-spike interval) cannot be computed. Because some of these properties are highly informative within the large group of multi-spikers, the inability to compute them makes it difficult to place the rapidly adapting cell models within the multi-spiker space. Similarly, the Hodgkin-Huxley model does not produce multiple spikes at 1.5x rheobase but does at 3.0x. It is known that electrical behavior of cells can differ under “normal” vs. “strong” stimulation conditions [71]. However, the definition of “normal” vs. “strong” is not clearly defined. Druckmann and others, as used here, defined it as 1.5x vs. 3.0x rheobase. But is there anything intrinsically special about these rheobase multiples? For example, could the 1.5x value be already too strong for some models, and the behavior at 3.0x will not be qualitatively different? A similar issue would occur if 3.0x were not strong enough to produce behavior different from the 1.5x stimulation. Ideally, these stimulation values would be based on the dynamical bifurcation structure of the cell membrane or cell model, sampling the regions of qualitatively different behaviors. An automated method which could identify such regions in a black-box manner (e.g. without knowing the governing equations) would greatly facilitate the electrophysiology assessment of both cell models and large numbers of cells.

Another example of models which do not easily lend themselves to analysis under the current protocol is the cluster of intrinsically spiking neuron models. Because these models spontaneously produce spikes, their rheobase currents are negative (e.g. some amount of hyperpolarizing current must be injected to reduce the number of spikes from non-zero to zero). Because the protocols used here use positive multiples of rheobase as stimulation, little information about spikes or spike trains of such models or biological cells can be gained from using the protocol. An understanding of the diversity of intrinsically spiking or bursting cells could be used to develop an automated experimental protocol for stimulation and feature extraction to assess such cells and their models.

Finally, neuron models in the Abstract Spiker cluster have unusual membrane potential properties such as positive resting potential or 0-width or undefined-amplitude action potentials (Fig 5). Given the variety of different possible abstract models, it’s not clear if a single protocol could be developed that would allow placing all such cells within the same space as the more physiologically realistic neuron models. One partial solution would be to identify separate sub-spaces for evaluating model spike train properties and action potential shape properties.

Identifying the function of Regular Multi-Spiker model clusters requires further comparisons to experimental recordings

Model clusters from the lowest level of the hierarchy form a tetrahedron in 3D space, where the leading features associated with the axes are: 1) delay to first AP, 2) steady state accommodation, and 3) median inter-spike interval of regular multi-spikers. This shape appears due to the strong linear relationships between delay and accommodation and between delay and median inter-spike interval within Regular Spiking models (RS, dRS, and naRS, Fig 6). The overall significance of this structure is not completely clear. If confirmed with larger and more diverse datasets of recordings from cells or simulation results from associated models, it may indicate fundamental dynamical system constraints among these three properties. AP onset, accommodation, and spike rate play important roles in neural representations of the magnitude and variation of an input signal. Since many of the models were constructed to exhibit predefined types of firing patterns (e.g., continuous accommodation or burst accommodation), the identified linear subspaces that capture the diverse types of firing patterns suggest possible compatibility between different modes of neuronal output.

The clusters that are missing are also interesting. For example, there are no high-density clusters of what could be called “non-accommodating slow spikers”, or “delayed fast spikers”, or “delayed bursters” (Fig 7). It’s possible that the tetrahedron is a product of an undersampling of potential models by the database. As additional, novel models are added, this structure might disappear. Similarly, because stochastic channels are currently not supported by NeuroML, models with stochastic firing patterns could affect the results. As the capabilities of NeuroML are further developed, the effects of such stochasticity could be tested. Finally, the tetrahedron might only reflect a structure within neuron model space, and it might not exist within the space occupied by experimental data. This could be tested with access to a diverse database of raw electrophysiology recordings, similarly to the Allen Cell Type Database but with a larger variety of represented brain regions. More speculatively, the tetrahedron might form a volume that real biological neurons occupy due to constraints on proper neuronal function. Future investigations of how neurons with electrical behavior that falls outside the tetrahedron affect network level dynamics or how neurons of various stages of disease migrate within this space could help test this hypothesis.

Large-scale analysis of models elucidates mechanisms across multiple scales

The joint visualization of cell model electrophysiological properties and their associated aggregated channel densities in this study is unique in that it incorporated additional information shared across models, e.g., the statistical structure of model channel densities across the database. While other high-dimensional visualizations of channel model mechanisms for neuron model databases exist [86,87], this study aimed to find common channel model mechanisms across a variety of cell and channel models spanning brain regions. Our approach demonstrated the feasibility of a functional and mechanistic mapping between the (in)activation properties of channel family model subtypes and the electrophysiological properties of cell model subtypes within the NeuroML database.

The primary goal of the combined analysis is to relate cell model electrophysiology to underlying channel mechanisms. We find that the two largest potassium channel model sub-types are associated with: 1) delayed rectifiers and slow non-inactivating potassium channels, and 2) A-type and slow inactivating potassium channels.

Intuitively, we can see that the ratio of these two discovered channel model sub-types indeed can affect the degree of accommodation (relative deceleration of spiking) of their associated cell models (Fig 8B). Further, the difference in activation thresholds provides useful information for interpreting the delay to first spike under “strong” current injection conditions (3x the rheobase). Kv Cluster B models are activated under weak stimulation while Kv Cluster A are activated only under strong stimulation of the cell model. This additional activation of outward currents from Kv Cluster A is one of the dominant causes of the delay to first spike (Fig 8B).

Future directions

The current release of the NeuroML database focuses on the characterization of ion channel and neuron model electrophysiology, neuron model morphology, and model complexity. Based on the requests of researchers, we could add additional features to the database to help make the selection and reuse of previously published models more efficient. In future releases, based on user demand, we will provide characterizations of synapse models, and similar dimensionality reduction and clustering analysis could be performed with synapse and cell morphology data. For network models, we could provide visualizations of their structure and connectivity and display traces of network outputs in response to specific inputs. Another useful feature might be to identify cell models in the database that have similar responses when provided with a set of user-uploaded electrophysiology recordings from experiments.

Finally, we’ve made the analysis code available online [88], which allows other researchers to compute the electrophysiology, morphology, and complexity properties from responses of other models or even experimental data.

NeuroML tools and simulator

jNeuroML v0.8.3 was used to parse and convert all NeuroML models to NEURON format (.mod and.hoc). All simulations were performed using NEURON simulator v7.5 on an Intel(R) Xeon(R) CPU E3-1240 V2 @ 3.40GHz, 32 GB RAM machine running Ubuntu Linux 16.04 LTS. For variable step simulations, the NEURON cvode_active variable was set to 1 before starting simulations.

Nested cell electrophysiology dimensionality reduction and clustering

Channel model electrophysiology dimensionality reduction, clustering, and parameter analysis