Generalized leaky integrate-and-fire models classify multiple neuron types

Abstract

There is a high diversity of neuronal types in the mammalian neocortex. To facilitate construction of system models with multiple cell types, we generate a database of point models associated with the Allen Cell Types Database. We construct a set of generalized leaky integrate-and-fire (GLIF) models of increasing complexity to reproduce the spiking behaviors of 645 recorded neurons from 16 transgenic lines. The more complex models have an increased capacity to predict spiking behavior of hold-out stimuli. We use unsupervised methods to classify cell types, and find that high level GLIF model parameters are able to differentiate transgenic lines comparable to electrophysiological features. The more complex model parameters also have an increased ability to differentiate between transgenic lines. Thus, creating simple models is an effective dimensionality reduction technique that enables the differentiation of cell types from electrophysiological responses without the need for a priori-defined features. This database will provide a set of simplified models of multiple cell types for the community to use in network models.

Fig. 1

Five generalized leaky integrate-and-fire (GLIF) models consisting of different phenomenological mechanisms are fit to electrophysiological data. A schematic describing the mechanisms is shown in a. Example data and models from two neurons of different transgenic lines are shown in b. For all models the input is a current, $I_e\left(t\right)$, injected via a patch electrode illustrated in black at the top of b. Below the current are the voltage traces from four repeats of the same stimulus (here colors represent the different responses to the repeated stimuli and do not adhere to the standard color scheme in the rest of the manuscript). Below the biological data, the GLIF models are plotted. The output of the models is the trans-membrane potential, $V\left(t\right)$, pictured in blue. When $V\left(t\right)$ reaches a threshold, $\Theta ={\Theta }_{\infty }+{\Theta }_S\left(t\right)+{\Theta }_v\left(t\right)$, shown in dashed green, a spike is produced, illustrated by blue dots. Note that the shape of the spike is not plotted as it is not fit by these models. Instead, after a refractory period, $V\left(t\right)$ is reset to a value dependent on the specific model. The GLIF₁ model is equivalent to the traditional LIF model with a refractory period where the model can not spike. This model contains one variable, $V\left(t\right)$, and the threshold is fixed to a value we refer to as Θ_∞. GLIF₂ models include a second variable: a spike-induced threshold ${\Theta }_S\left(t\right)$ which is added to the baseline threshold Θ_∞. When the model spikes, ${\Theta }_S\left(t\right)$, jumps up and then decays. Thus, after a spike, initially the total threshold is higher making it harder for the model to reach threshold. GLIF₃ includes $V\left(t\right)$ and two variables corresponding to two spike initiated after spike currents, $I_j\left(t\right)$, which have different time constants and decay back to zero. The sum of the after-spike currents are illustrated in red. GLIF₄ combines GLIF₂ and GLIF₃ for a total of four variables. GLIF₅ includes an additional threshold component ${\Theta }_v\left(t\right)$. ${\Theta }_V\left(t\right)$ is dependent on the voltage of the model. Scale bars represent all model plots (not the amplitude of the current injection or biological voltage traces)

Fig. 2

Overall, 645 different neurons from 16 transgenic lines containing all the required stimuli on the Allen Cell Types Database are considered in this study. (Left) Illustrated colors correspond to the different transgenic lines in all figures. “n” describes the number of neurons for which the lowest level model (GLIF₁) could be generated. Transgenic lines can identify either or both inhibitory (I) or excitatory (E) cells, which reside in layer 1 (L1) through layer 6 (L6). Note that most Chrna2-Cre-positive neurons are inhibitory and thus are labeled here as inhibitory. (Right) A minimal set of stimuli were required for training and testing different GLIF models. GLIF models were trained using a at least two repeats of pink noise stimuli (3 s each, 1/f distribution of power, 1 – 100 Hz) with amplitudes centered at 75, 100, and 125 percent of action potential threshold, b a short (3 ms) just supra-threshold pulse to fit the instantaneous threshold, Θ_∞, c a long square (1 s) pulse just below threshold to estimate the intrinsic noise present in the voltage traces (used in the post hoc optimization step), and d a series of three peri-threshold short pulse sets for any model with reset rules (GLIF₂, GLIF₄, and GLIF₅). GLIF models were then tested using a hold-out stimulus set e of at least two sweeps of a second pink noise stimuli generated in an identical manner to the training but initialized with a different random seed. Representative data shown from an Htr3a-Cre-positive neuron (resting membrane potential (RMP) −67 mV) and a Ctgf-Cre-positive neuron (RMP = −74 mV). Scale bar in a corresponds to 40 mV (a, e), 50 mV (b, d), or 20 mV (c)

Fig. 3

Slices from the parameter space fit from electrophysiological data. Detailed parameter fitting methodology (Fig. 2) is available in the Supplementary Methods. Stars denote example neurons listed in the main text. a Resting potential is measured as the average voltage during rest before training noise (noise 1) current is injected. The threshold relative to rest, ΔV, is measured by subtracting resting voltage from the threshold obtained from the supra-threshold short square pulse. b The spike waveform is removed from the voltage trace by aligning all spikes and fitting a line to the voltage before and after a spike. The best fit line within a window of 10 ms after spike initiation was chosen (Supplementary Fig. 1). This spike cut length is used in all models, and voltage measurements before and after a spike are used to reset voltage in (GLIF₂, GLIF₄, and GLIF₅). c Capacitance and resistance are fit via linear regression to sub-threshold voltage (Supplementary Fig. 2). The membrane time constant, τ = RC is plotted. d Total charges of the fast, Q₁, and slow, Q₂, after-spike currents deposited each time there is a spike. e The amplitude, a_s, and decay, b_s, of the spiking component of the threshold, Δ_s, (used in GLIF₂, GLIF₄, and GLIF₅) is fit to the triple short square data set (Supplementary Fig. 3). f In the GLIF₅, the threshold is influenced by the voltage of the neuron according to (Equation 4). The two parameters of Eq. 4 are plotted here. Colors in all panels correspond to transgenic lines illustrated in Fig. 2

Fig. 4

Rastergrams and explained variance of biological data and all optimized model levels for “hold out” test data. a Data for the two example cells. Injected current shown in black. Black rasters are spikes from recorded neurons to repeated current injections. Colored rasters correspond to the five different, deterministic models. The current injection is 3 s long. As GLIF₁ and GLIF₂ do not have a spike frequency adaptation mechanism (implemented by the summation of after-spike currents over many spikes in GLIF₃), they have trouble reproducing simultaneously the firing patterns at multiple input amplitudes. b Explained variance for the different model levels at different levels of time window resolution Δt. The black lines represent the explained variance of the data (how well the neuron repeats its own spiking behavior). This trace would reach 100% if the spike times of the data were all exactly at the same time in each repeated stimulus: the fact that they are not 100% reflects the intrinsic variation in the spike times within the experimental data. The blue line illustrates the pairwise explained variance of the model with the data. Because the model can not be expected to explain the data better than the data can explain itself, the red line is the ratio of the pairwise explained variance of the model (blue) and the data divided by the explained variance of the data (black). This ratio value at a Δ t = 10 ms time bin is used for the explained variance performance metric in the main text

Fig. 5

Different mechanisms improve model performance for inhibitory and excitatory neurons. The traditional leaky integrate and fire model (GLIF₁) yielded surprisingly high model performance. Overall, inhibitory models were more successful at reproducing spike times than excitatory models. Reset rules implemented on their own (GLIF₂) decreased model performance. After-spike currents (GLIF₃) improved inhibitory model performance, whereas a combination of both after spike currents and reset rules (GLIF₄) were required to gain performance of excitatory models. The voltage-dependent adapting threshold (GLIF₅) improved performance of excitatory models even more, but had only a slight effect on inhibitory models. The thick blue line denotes all excitatory neurons, the thick red line denotes all inhibitory neurons, and the thick black line is for all neurons. Thin lines are different transgenic lines. Accompanying data is available in Table 3. p-values for significant differences between GLIF model levels can be found in Supplementary Figures 8, 9, and 10. Briefly, for the “all”, “excitatory”, and “inhibitory” groupings the p-values are smaller than 0.01 (and often much smaller) for all but between the excitatory GLIF₂ and GLIF₃ (Supplementary Fig. 8). Differences between GLIF levels of different transgenic lines are sometimes statistically significant and sometimes not (Supplementary Figures 9 and 10)

Fig. 6

We identify discrete putative clusters using an iterative binary clustering approach on 645 cells. The top six panels show the summary of clusters obtained by iterative binary clustering using electrophysiological features extracted from the traces and GLIF model parameters. In every panel, each row represents a cluster, and each column a transgenic line. The size of the circle indicates the fraction of cells from a given transgenic line falling into a specific cluster (such that the sum of fractions in a column add up to 1). The dendrogram on the y-axis shows the iterative binary splitting into clusters using the algorithm explained in the text. For each intermediate node, a support vector machine was trained on half the cells at that node and used to classify the remaining cells. The number at each node indicates the minimum percentage of test cells correctly classified over 100 iterations of randomly selected training and test cells. Clustering based on features and using the GLIF₃ and GLIF₄ model parameters shows separation among lines labeling inhibitory and excitatory cells. In addition, transgenic Cre lines marking Pvalb+, Ntsr1+, Nr5a1+, and Ctgf+cells tend to segregate into distinct clusters. The bottom two panels show two measures of overall clustering similarity: the adjusted Rand index (ARI) in red, and the adjusted variation of information (AVOI) in black. The bottom left panel shows similarity between each set of clusters and the transgenic lines. The bottom right panel shows similarity between each set of clusters and the clusters obtained using the features. An ARI of 1 indicates perfect agreement between partitions, whereas 0 or negative values indicate chance levels of agreement. A positive value of the AVOI indicates agreement between partitions that is better than chance (which is indicated by 0). The gray and pink traces in these two panels show the AVOI and ARI values, respectively, for random subsets of the features containing the same number of parameters as each of the four GLIF models

Fig. 7

Iterative binary splitting clustering obtained from using GLIF model parameters plus spike-shape-related feature parameters. The top four panels show cluster versus Cre line composition, similar to Fig. 6. The bottom two panels show the adjusted variation of information metric and the adjusted Rand index for the GLIF model-derived clusters with and without the spike-shape parameters