Single-neuron models linking electrophysiology, morphology, and transcriptomics across cortical cell types

Abstract

Which cell types constitute brain circuits is a fundamental question, but establishing the correspondence across cellular data modalities is challenging. Bio-realistic models allow probing cause-and-effect and linking seemingly disparate modalities. Here, we introduce a computational optimization workflow to generate 9,200 single-neuron models with active conductances. These models are based on 230 in vitro electrophysiological experiments followed by morphological reconstruction from the mouse visual cortex. We show that, in contrast to current belief, the generated models are robust representations of individual experiments and cortical cell types as defined via cellular electrophysiology or transcriptomics. Next, we show that differences in specific conductances predicted from the models reflect differences in gene expression supported by single-cell transcriptomics. The differences in model conductances, in turn, explain electrophysiological differences observed between the cortical subclasses. Our computational effort reconciles single-cell modalities that define cell types and enables causal relationships to be examined.

Keywords: CP: Neuroscience; cell types; dimensionality reduction; electrophysiology; high-performance computing; machine learning; modeling; morphology; multimodal cellular data; optimization; transcriptomics.

Figure 1.. Cell types in mouse visual cortex (VIS) and single-cell model generation

(A) Data modalities: single-cell transcriptomics (left) (Tasic et al., 2016), (right) single-cell morpho-electric (ME) data (Gouwens et al., 2019). Sample morphologies arranged according to normalized depth from pia, including putative cortical layer markings and electrophysiological recordings for 3 cells in the same set under subthreshold and suprathreshold current injections (colors: Cre line). (B) Overview of the 230 modeled cells in the ME dataset based on dendrite type and Cre line (cell class— spiny: putative excitatory, sparsely spiny, and aspiny: putative inhibitory). See also Figure S1. (C) Model setup: Morphology (CellID: 483101699) and a predefined set of active conductances/passive properties marked according to their inclusion in each morphology section (apical, basal dendrites, soma, AIS). (D) The 3-stage optimization workflow (stages 0, 1, 2). The parameters added to the variable list at each stage are highlighted in (C) (table, left). (Top) Evolution of the sum of objectives with generation number. The best model (gray) at each generation and the average performance of all of the individual models (red) of that generation (spread: standard deviation of the population). (Bottom) Comparison between experimental traces (black) and the fitted model (gray) at each stage for a representative model (same as in C). (E) Validation comprising novel stimulation protocols (top; noise stimulation) and comparison between model and experiment spiking responses (bottom) under the novel noise stimulus (i.e., standardized colored noise [top]) to evaluate explained variance (see STAR Methods).

Figure 2.. All-active models conserve decision boundaries and cross-class relationships of experimental features that reproduce key differences within and between cell classes

(A) Left: The parameters across the 40 hall of fame (HOF) models corresponding to each cell creates a cluster in ñ-dimensional euclidean space (ñ : common parameters across all cells). The centroid of this cluster (color: euclidean distance between HOF index and its centroid). The sum of all distances measures cluster dispersion (bar plot top of the distance matrix). (Right) t-Distributed stochastic neighbor embedding (t-SNE) plot highlights 4 cases of maximally (red: excitatory cell; green: inhibitory cell) and minimally (blue: excitatory cell; magenta: inhibitory cell) dispersed parameter vectors in 2 dimensions (2D). (B) Left: Heatmap showing euclidean distances between HOF models for the 3 broad inhibitory classes; each block diagonal represents 40 HOF models for a single cell. Right: Interpretation for different regions of the heatmap, namely, intracell intraclass (blue), intercell intraclass (orange), and intercell interclass (green). Darker block diagonals and the distribution of the intra/inter distances for the 3 broad classes indicate an ordered structure in the parameter dispersion. Degenerate parameters for a single cell are tightly clustered compared to parameters of the same broad class followed by model parameters between different classes (p < 0.01; Mann-Whitney U test). (C) For cells belonging to the 4 broad subclasses, HOF model E-features at the maximal amplitude stimulus protocol and projected onto UMAP embedding of the corresponding features at the experiment level (optimal number of clusters: n_clusters = 3 via gap statistic (Tibshirani et al., 2001). The corresponding k-means decision boundary is drawn on the embedded space for the experiment, the best model (HOF index = 0), and all HOF models (40 models per cell). The number of detected clusters, their composition, and associated decision boundaries remain unaltered. (D) UMAP embedding of HOF parameters colored according to the broad subclasses. (E) Similar to (D), the separation in the electrophysiology and model parameter space is also preserved between two putative excitatory types, namely L2/3 and L5 PCs. (F) The separation in the HOF model parameter space between two excitatory L5 excitatory types, IT and PT.

Figure 3.. Model-based prediction in Kv3.1 differences between GABAergic cell classes are supported by single-cell RNA-seq and explain divergent electrophysiology properties

(A) Comparison between single-cell model, somatic ion channel conductance profiles, and single-cell expression of associated genes via single-cell RNA-seq for Kv3.1 (associated gene: Kcnc1); boxplot: median, 1.5 interquartile range; cpm, counts per million. (B) Expression profiles of marker genes (Pvalb, Sst, Vip) and specific ion channel (Kv3.1, KP – persistent K and KT – transient K) genes (Kcnc1, Kcna1–3,6 Kcnd1–3, respectively) for a set of GABAergic lines in mouse V1 via single-cell RNA-seq. (C) Ion channels predicted by the models and confirmed by sequencing to differentially express between inhibitory cell types affect E-features. Sensitivity analysis of single-cell models (line: mean; error bar: 95% confidence interval). (D) Top: Analysis of E-features from in vitro experiments. Bottom: Class-specific characteristics are preserved in the corresponding all-active models, providing a validation step for the ion channel predictions (boxplot, circles: single-cell data). Statistical testing: Mann-Whitney U test; statistical significance: *p < 0.05, **p < 10⁻², ***p < 10⁻³, adjusted for a false discovery rate (FDR) of 5%.

Figure 4.. Model-based predictions in Ih and NaT differences between 3 excitatory subclasses link transcriptomic and electrophysiological properties

(A) The 3 excitatory subclasses L4 IT, L5 IT, and L5 PT highlighted in the transcriptomic tree derived from single-cell transcriptomics (Tasic et al., 2016). (B) Gene expression of the top differentially expressed (marker) genes between the 3 excitatory subclasses (Rspo1, Fezf2, Fam19a1) and ion channel genes implicated by the model-based predictions via pairwise comparisons: Ih (HCN1–3) and transient Na (Scn8a) (heatmap). (C) Left: h-Channel conductance density comparison between Nr5a1 (L4 IT) and L5 Rbp4 (L5 PT) models (p < 0.05, 1-sided Mann-Whitney U test; Cliff’s delta effect size analysis of the medians (Ho et al., 2019). Right: Same comparison for gene expression from the single-cell transcriptomics between Nr5a1 (L4 IT) and L5 Rbp4 (L5 PT). (D) Analysis of sag in the E-responses. Sag_ratio (Table S2) comparison at −80 pA between (left) all cells of the two Cre lines, (center) the subset of experiments for which models were generated, and (right) the corresponding model responses. (E) Pairwise comparison between L5 IT and PT models reveals elevated somatic, apical h channel, and transient axonal Na conductance density for L5 PT cells (with the same Hcn1, Scn8a expression pattern between the 2 classes; line, mean; error bar, 95% confidence interval) (F) Experimental E-features comparison between L5 IT and PT cells (line, mean; error bar, 95% confidence interval). Statistical testing: Mann-Whitney U test; statistical significance: *p < 0.05, **p < 10⁻², ***p < 10⁻³, adjusted for a false discovery rate (FDR) of 5%.

Figure 5.. Increased apical Ih leads to preservation of synaptic patterns in deep versus superficial pyramidal (Pyr) neurons

(A) Left: Pairwise comparison between L2/3 and L5 excitatory models across conductances. Right: Single-cell transcriptomics comparison. (B) Simulation of synaptic activation along the apical dendrites of L23 and L5 Pyr neurons. Individual synapses are distributed and activated along the apical dendrite, and their postsynaptic effect is measured at the soma. Incoming presynaptic spikes along the apical dendrite (left column) and the resultant post-synaptic somatic depolarization (right column) for a single cell (CellID: 488679042). (C) Postsynaptic somatic depolarization for superficial and deep Pyr cell types normalized by the peak postsynaptic potential (PSP) amplitude when the synapse is at soma (line: mean; shaded area: standard deviation; green: L2/3 Pyr models; purple: L5 Pyr models; 15 L2/3 models: 33 L5 models; mean peak PSP amplitude of 0.75 mV and 0.56 mV for L2/3 and L5 models, respectively; synapse location: at soma and ~90 μm from soma on the apical dendrites). (D) Comparing against the nominal case (e.g., postsynaptic event if presynaptic spike was directly injected into the soma), the percentage change in somatic postsynaptic attributes (left to right: somatic PSP response amplitude, width, and latency) as function of synapse location along the apical dendrite (line, mean; error bar, standard deviation).

Figure 6.. All-active models offer increased bio-realism

(A) Performance of 112 spiny, 109 aspiny, and 5 sparsely spiny all-active models. Left: E-features for the best models ranked according to performance (Z score) during training. Right: The histogram of the training error (Z score) for each E-feature and the associated lognormal fit. (B) Reconstructed neuron morphology (CellID: 483101699). (C) Both all-active (blue) and perisomatic (red) models capture the experimental (black) somatic voltage response, but for all-active models, spikes are initiated from the AIS; in perisomatic models, spikes start at the soma. Arrows denote the propagation direction of the spike. Right: Simulated recordings at AIS and soma depict the direction of spike propagation. ${\text{Δ}}_{\text{t}}:={\text{t}}_{\text{soma}}^{\text{i}}-{\text{t}}_{\text{axon}}^{\text{i}}$, for the ith spike, is positive for all-active models and negative for perisomatic. (D) Performance of all models on E-features (x axis: experimentally measured feature; y axis: model-produced feature after training). Each point corresponds to the feature from a single stimulus protocol. Data along the diagonal of each panel reflect 1-to-1 agreement with experimental values. (E) Parameter combinations for all-active models for spiny (left) and aspiny (right) are different from their perisomatic counterpart (bar, mean; error bar, standard deviation). (F) The relative contribution of AIS and somatic ion channel conductances for each intracellular somatic feature derived from the sensitivity analysis (parameter perturbation: 10% about optimized value; relative contribution of each conductance expressed in terms of Sobol indices). Results shown for all-active models of 10 excitatory cells (bar: mean; error bar: standard deviation).