Modeling functional cell types in spike train data

Abstract

A major goal of computational neuroscience is to build accurate models of the activity of neurons that can be used to interpret their function in circuits. Here, we explore using functional cell types to refine single-cell models by grouping them into functionally relevant classes. Formally, we define a hierarchical generative model for cell types, single-cell parameters, and neural responses, and then derive an expectation-maximization algorithm with variational inference that maximizes the likelihood of the neural recordings. We apply this "simultaneous" method to estimate cell types and fit single-cell models from simulated data, and find that it accurately recovers the ground truth parameters. We then apply our approach to in vitro neural recordings from neurons in mouse primary visual cortex, and find that it yields improved prediction of single-cell activity. We demonstrate that the discovered cell-type clusters are well separated and generalizable, and thus amenable to interpretation. We then compare discovered cluster memberships with locational, morphological, and transcriptomic data. Our findings reveal the potential to improve models of neural responses by explicitly allowing for shared functional properties across neurons.

Fig 1.

A: The Poisson GLM (2) used to model the spiking response of a single neuron. B: The generative model (11) for the simultaneous method.

Fig 2.. Performance of sequential and simultaneous methods on simulated data.

A: The true stimulus filter ${\mathit{\beta }}_i^{\text{stim}}$ for each simulated neuron. B: The true cluster means ${\mathit{\mu }}_k^{\text{self}}$ used to generate simulated datasets, and those estimated by the sequential and simultaneous methods, ${\widehat{\mathit{\mu }}}_k^{\text{self}}$, fit with the correct K = 5. C-F: Mean ± SEM over 50 simulated datasets of accuracy measures, as a function of σ and shown for both K = 3 and K = 5. When K = 3, the three clusters with leftmost ${\mathit{\mu }}_k^{\text{self}}$ in panel B are used. For each condition (K and σ) and for each measure of accuracy, the simultaneous method’s performance is statistically significantly better than that of the sequential method, except for ARS with K = 3 and σ = 10^−5/6 (evaluated using the Wilcoxon signed rank test: uncorrected P-value<0.002).

Fig 3.. Model selection of K^.

Frequencies of $\widehat{K}$ estimated via Bayesian information criterion over 50 simulated datasets with the same μ_k as in Figure 2, and σ = 10⁻² (A) or 10^−5/6 (B), the maximum value that does not result in degenerate simulations. Black lines indicate the true value of K. The summary below each plot reports the mean ± SEM of the estimated value of $\widehat{K}$ across the 50 datasets, for each value of K and each method.

Fig 4.. IVSCC Dataset: The simultaneous method explains the data with a smaller number of less overlapping clusters.

A: BIC for the simultaneous method over a range of K; BIC determines $\widehat{K}=12$ as optimal. Each dot reports the result of running Algorithm 2 from a different random initialization. B: Cluster centers μ_k of self-interaction filters fit to data using the simultaneous method; shaded region is $\pm \sqrt{diag({\mathrm{\Sigma }}_k)}$. Only the 9 clusters with at least 20 neurons are shown from the model with BIC-selected K = 12. C: By contrast, the standard BIC of a GMM fit to individually fitted self-interaction filters (the sequential method) suggests an optimal $\widehat{K}$ of at least 19. Each dot reports the result of performing the GMM fit, (5), from a different random initialization. D: The clusters of self-interaction filters with at least 20 neurons found by the sequential approach with K = 12 are bunched up closer to the origin, such that the clusters overlap significantly.

Fig 5.. IVSCC generalization performance: the simultaneous method produces single-cell models and clusterings that generalize better, especially when fitted to more neurons.

Additional details about this figure are available in Section S3. A: ANLL (lower is better) for each held-out neuron’s single-cell model, evaluated on responses to the test stimulus (Noise 2), using the MAP ${\widehat{\mathit{\beta }}}_i$ (20) of the simultaneous method with (hyper)parameters $K$, ${\lambda }^{\text{stim}}$, ${\widehat{\mathrm{\Omega }}}_K$ estimated from the training neurons, versus those found by the sequential method. Color encodes number of spikes for each neuron in the evaluation data. B: Median relative difference between methods in ANLL of held-out neurons, evaluated on responses to the test stimulus, as a function of how many neurons and how much data from each were used in training (more negative values indicate that the simultaneous method is better). White asterisks indicate a significant relative difference; white bars indicate adjacent cases of training data subselection where the relative differences were significantly different. Differences pooled across all vertically (horizontally) adjacent conditions showed a significant, p = 3 × 10⁻⁴, (p = 5 × 10⁻²⁶) trend, with more presentations (neurons) yielding greater improvement by the simultaneous method. C: Same analysis as A, but with EV_ratio (see (23); higher values indicate that the simultaneous method is better). D: Same analysis as B, but with EV_ratio. Pooled vertical differences showed no significant trend (p > 0.1); horizontal differences showed a significant (p = 5 × 10⁻³) trend, with more neurons yielding greater improvement by the simultaneous method. E: Relative differences between methods of EV_ratio (shown in C) versus ANLL (shown in A); color encodes number of spikes. Neurons with many (few) spikes show only improved ANLL (EV_ratio) in the simultaneous method. F: Same analysis as B, but for the similarity of cluster assignments ${\widehat{k}}_i$ between model fits with different held-out neurons, measured by ARS (more positive values indicate that the simultaneous method is better). Pooled vertical differences showed a significant (p = 5 × 10⁻²) trend, with fewer presentations yielding greater improvement by the simultaneous method; horizontal differences showed no significant trend (p > 0.1).

Fig 6.. IVSCC metadata is related to discovered cell types.

Z-scored fraction of cells with an attribute in each cluster. Cluster identities in panel A (B) are the same as in Figure 4B (D), obtained using the simultaneous (sequential) method (fitted with BIC-selected K = 12 clusters, showing only clusters with at least 20 neurons). Attributes are spiny or aspiny dendrites, location (hemisphere and cortical layer), and Cre line. Note the ARS between the cluster and metadata labels in each title. Z-scores are calculated as $Z_i^{(a)}=({\widehat{p}}_i^{(a)}-{\widehat{p}}_i)/\sqrt{{\widehat{p}}_i^{(a)}(1-{\widehat{p}}_i^{(a)})/N^{(a)}+{\widehat{p}}_i(1-{\widehat{p}}_i)/N}$, where ${\widehat{p}}_i$ is the empirical probability that a cell is in cluster i and ${\widehat{p}}_i^{(a)}$ is the empirical probability that a cell with attribute a is in cluster i, N is the number of cells, and N^(a) is the number of cells with attribute a.