Combined statistical-biophysical modeling links ion channel genes to physiology of cortical neuron types

Abstract

Neural cell types have classically been characterized by their anatomy and electrophysiology. More recently, single-cell transcriptomics has enabled an increasingly fine genetically defined taxonomy of cortical cell types, but the link between the gene expression of individual cell types and their physiological and anatomical properties remains poorly understood. Here, we develop a hybrid modeling approach to bridge this gap. Our approach combines statistical and mechanistic models to predict cells' electrophysiological activity from their gene expression pattern. To this end, we fit biophysical Hodgkin-Huxley-based models for a wide variety of cortical cell types using simulation-based inference, while overcoming the challenge posed by the mismatch between the mathematical model and the data. Using multimodal Patch-seq data, we link the estimated model parameters to gene expression using an interpretable sparse linear regression model. Our approach recovers specific ion channel gene expressions as predictive of biophysical model parameters including ion channel densities, directly implicating their mechanistic role in determining neural firing.

Figure 1. Sketch of the statistical-biophysical hybrid model.

Neuronal gene expression levels (left) as well as electrophysiological patterns (right) are obtained experimentally with Patch-seq. The electrophysiology is fitted with a conductance-based biophysical model (middle right). The estimated model parameters are then predicted with sparse reduced-rank regression (middle left), completing the statistical-biophysical bridge from neuronal genotype to phenotype.

Figure 2. Exemplary experimental observations and their closest simulations from the prior.

a Middle: t-SNE embedding of $n=955$ MOp cells based on their transcriptome. Surround: one example neuron for each of the six transcriptomic families. For each neuron, we show the experimental observation (below) and the biophysical model simulation from the prior with the smallest Euclidean distance in standardized electrophysiological feature space (above). b Comparison of nine electrophysiological feature values between experimental observations and best prior simulations shown in (a).

Figure 3. Neural Posterior Estimation vs Neural Posterior Estimation with Noise

a The MAP parameter set simulation derived with NPE-N is closer to the experimental reference (in blue below, L4/5 IT_1 pyramidal cell) than derived with NPE. Residual distance of the MAP parameter set simulation to the experimental observation (model fit error) shown for each electrophysiological feature (0 corresponds to a perfect fit). b 1- and 2-dimensional marginals together with 3 simulations generated from parameter combinations with highest probability under the posterior (out of 10 000 samples); NPE (left) vs NPE-N (right) setting. 7 out of 13 model parameters have been selected for illustration.

Figure 4. Neural posterior estimation of conductance-based model parameters in the presence of model misspecification.

a Sketch illustrating model misspecification: in electrophysiological feature space, not enough simulations cover the space of experimental observations. NPE-N introduces isotropic noise to the summary statistics of simulations (dotted arrows). b Simulations are further away from experimental observations (blue) than from other simulations (orange). Qualitatively, simulations increasingly further away from an experimental observation look more dissimilar than from another simulation. Numbers 1 and 2 refer to the 2nd and 501st closest simulations respectively.

Figure 5. Two-dimensional embedding reveals difference in HH-based parameters between neural families.

a T-SNE embedding of $n=955$ MOp neurons based on transcriptomic data. Colors like in Fig. 2a. Cells in the middle of the embedding had lower quality transcriptomic data and therefore grouped together. b Marker genes expression levels overlayed and interpolated on embedding confirm known families (dark purple: low expression, yellow: high expression). c Uncertainty of MAP parameters for each cell overlayed on the embedding. The uncertainty was calculated as the posterior entropy $-{\sum }_{k=1}^{1000}\text{log}{q}_{\varphi }\left({\theta }_k\mid {\mathbf{x}}_o\right)$, where we sampled ${\theta }_k~q_{\varphi }\left(\theta \mid {\mathbf{x}}_o\right)$. d Selection of summary statistics derived from simulations corresponding to MAP estimates, overlayed on the embedding. e Selection of summary statistics describing observed electrophysiology, overlayed on the embedding. f Selection of MAP parameters, overlayed on the embedding.

Figure 6. Prediction of MAP parameter estimates from gene expression with sparse reduced-rank regression.

a sRRR schematic. A linear combination of selected genes is used to predict fitted HH-based model parameters. b Cross-validation performance for rank-2 and full-rank sRRR models with elastic net penalty. The dashed vertical line shows the performance with 25 genes. c Rank-2 sRRR model predictive performance for each model parameter, using the entire data set. d Middle: rank-2 sRRR model latent space visualization. All 955 MOp neurons are shown. Left: Selected ion channel and marker gene overlays. Right: Predicted model parameter overlays.

Figure 7. MAP parameter estimates and sRRR predictions for each family and cell type.

a MAP parameter estimates averaged over cells belonging to a family and belonging to a transcriptomic cell type (left and right respectively). We Z-scored all values by subtracting and scaling with the mean and standard deviation of $n=955$ MAP estimates, respectively. b Analogous to a, but with rank-2 sRRR predicted model parameter values.

Figure 8. Family representation of MAP estimates together with sRRR predictions.

a Analogous to Fig. 2a except that the simulation on top is derived from the family-average MAP estimate calculated as in Fig. 7a, left. Simulation on the bottom is derived from the family-average sRRR prediction calculated as in Fig. 7b, left. b Comparison of 9 electrophysiological feature values derived with the MAP estimate versus sRRR-based estimate.