An important question in neuroscience research is understanding the functional connectivity between neurons in different parts of the brain. Spike trains based on simultaneously recorded neurons provide information about population coding and neuronal interaction. Both model-free and modelbased spike-train analysis tools have been developed to answer this question. While model-free methods (Perkel et al., 1967;Ventura et al., 2005;Fujisawa et al., 2008;Humphries, 2011;Lopes-dos-Santos et al., 2013) are typically more efficient and convenient to implement, they often fail to uncover more complex underlying neuronal relationships beyond correlation at the level of observed data. In contrast, latent factor models (LFMs) are able to discover patterns which modelfree algorithms cannot, thanks to their ability to specify different structures in the latent layers. The remarkable success of spike-train LFMs in predictive tasks (Yu et al., 2008;Gao et al., 2016;Wu et al., 2017;Pandarinath et al., 2018) and neuronal clustering (Buesing et al., 2014;Wei et al., 2022) motivates our work in this paper to study the functional connectivity between neurons. As all the aforementioned LFMs are specified in discrete time, they are commonly applied upon binning the experimental time to obtain spike counts, which are then modelled using, e.g., a Poisson likelihood. However, the bin size is often chosen arbitrarily-despite this having a significant impact on parameter estimation (Nelson, 2002;Kass & Ventura, 2006;Ramezan et al., 2014). Although there are Poisson process models set in continuous time without binning (e.g., Duncker & Sahani, 2018;Williams et al., 2020), these models cannot be readily applied to the functional connectivity analysis of multiple neurons.
In this work we present a continuous-time multivariate point process LFM to study neuronal interactions based on simultaneously recorded spike trains in a neural population. To the best of our knowledge, this is the first continuous-time LFM proposed for the analysis of neural spike trains. In our model, the activities of a neuronal population are described by correlated Wiener processes with resetting. Each of these processes is viewed as a proxy of the evolving membrane potential of a neuron which resets after reaching a threshold. Crucially, we assume that the high-dimensional multivariate latent process can be summarized by a small number of dynamic factors. Not only does this factor analysis framework provide an interpretable low-dimensional representation of neural activities, it also serves as a means of dimension reduction for studying large neuronal populations. Our model generalizes the limiting case of the multivariate Skellam point process with resetting of Ramezan et al. (2022). However, by passing to the Brownian limit, we are able to develop an efficient algorithm for parameter inference, which both circumvents the choice of bin size (it can be made so small as to approximate the continuous time process arbitrarily well) and scales linearly in the number of neurons in the analysis. The applicability of the proposed model and inference procedure is demonstrated in simulated and experimental data analyses.
A graphical summary of the proposed model is presented in Figure 1. Let Y t = (y 1,t , . . . , y q,t ) ∈ {0, 1} q denote the observed binary spike trains of a population of q neurons at time t, and X t = (x 1,t , . . . , x q,t ) ∈ R q be the unobserved latent dynamic processes of the neuronal population. We assume that the neuronal activities are governed by X t , which are modelled as correlated Wiener processes with resetting. Each process i is reset to its initial value whenever it crosses some neuron-specific threshold, and then a spike occurs, i.e., y i,t = 1 is observed. We let θ θ θ denote the threshold parameters and the Wiener process drift vector. The threshold-crossing and resetting processes mimic the spike-generating mechanism and the refractoriness of the neuron, respectively. We further assume that the q dynamic processes can be represented as X t = Λ Λ ΛF t + ε ε ε t , where F t = ( f 1,t , . . . , f d,t ) are d (d ≪ q) dynamic factors, and ε ε ε t = (ε 1,t , . . . , ε q,t ) are neuron-specific idiosyncrasies. The loading matrix Λ Λ Λ can be interpreted within the traditional factor analysis framework-that is, it identifies a small number of factors driving the neuronal dynamic processes. It is also used to model the correlation matrix between latent neuron processes, i.e., cor(X t ) = Σ Σ Σ = Λ Λ ΛΛ Λ Λ ′ + Ψ Ψ Ψ, where Ψ Ψ Ψ is a diagonal matrix with elements determined by Λ Λ Λ and the fact that Σ Σ Σ has unit diagonal.
Model inference is carried out in two steps. First, the drift of the Wiener process and the threshold parameters, sum- marized by θ θ θ, can be estimated analytically since the first passage time of a Wiener process follows an Inverse Gaussian distribution (Brown, 2005). However, closed-form estimation of Λ Λ Λ is not available. One common solution is to employ MCMC sampling on p(X 1:T , Λ Λ Λ | Y 1:T ), where time has been discretized to an arbitrarily fine grid. However, this can be prohibitively slow in high dimensions. Instead, the latent processes X 1:T are integrated out via the Laplace approximation (Tierney & Kadane, 1986;Skaug & Fournier, 2006;Koyama et al., 2010)
so that an estimate of Λ Λ Λ can be obtained by maximizing the Laplace-approximated marginal likelihood p LA (Y 1:T | Λ Λ Λ).
An efficient gradient-based nested optimization algorithm for this is implemented in our R/C++ library fastr (Chen et al., 2023). The algorithm involves repeatedly solving linear systems of the form Σ Σ Σx = b, which due to the factor structure Σ Σ Σ = Λ Λ ΛΛ Λ Λ + Ψ Ψ Ψ, scales linearly in the number of neurons q for fixed number of factors d.
To assess the performance of the proposed model, we first applied it to a biologically plausible simulated dataset generated in NEURON (Hines & Carnevale, 1997). We found that the K-means clustering results based on Λ Λ Λ are far more ac- curate than those based on (convolved) multiple neuron spike trains Y 1:T (see Figure 2(a)-(c)). Therefore, the estimate Λ Λ Λ
uncovers the underlying associations between neurons beyond what is observed at the data level. Next, we applied our model on neuronal ensembles recorded from a rat's medial prefrontal cortex (mPFC) and orbitofrontal cortex (OFC) during a classical conditioning experiment. The rat was trained to recognize a cue predictive of an appetitive outcome (sucrose).
We modelled neuronal population dynamics in relation to encoding appetitive outcomes and functional connectivity, and found that coordinated activities are modulated by brain areas
We have proposed a novel continuous-time multivariate point process latent factor model for simultaneously recorded spike
• Figure 2: Clustering performance. The (i, j)-th entry of the matrix is black if neuron i and j are assigned to the same cluster, and white otherwise.
• • Figure 3: Estimated correlation between neurons. "O" stands for OFC neurons and "M" stands for mPFC neurons. The numbers are the neuron indices. For visual presentation, nonsignificant values and the diagonal elements are set to zero. trains. Downstream analyses using the proposed model can reveal neuronal clustering and estimate correlations between neurons. Computational challenges in model inference are addressed by carefully designing and implementing an efficient model-fitting procedure based on the Laplace approximation. We have confirmed, via simulation studies, that our algorithm achieves satisfactory accuracy and precision in parameter estimation. Finally, compared to black-box deep learning models for neural spike trains, our proposed model is able to provide more reliable statistical inference results with uncertainty quantification, which is vital for making scientifically sound conclusions. One immediate future direction is to apply our model to appetitive and aversive classical conditioning outcomes (available within the same experimental data) to investigate reward-value coding, and to identify value signals for the most relevant contextual features.