Point process models for sequence detection in high-dimensional neural spike trains

Point processes are probabilistic models that generate discrete sets ${\{x_{1},x_{2},\ldots\}\triangleq\{x_{s}\}_{s=1}^{S}}$ over some continuous space $\mathcal{X}$. Each member of the set $x_{s}\in\mathcal{X}$ is called an “event.” A Poisson process [58] is a point process which satisfies three properties. First, the number of events falling within any region $V\subset\mathcal{X}$ is Poisson distributed. Second, there exists a function $\lambda(x):\mathcal{X}\mapsto\mathbb{R}_{+}$, called the intensity function, for which $\int_{V}\lambda(x)\,\mathrm{d}x$ equals the expected number of events in $V$. Third, the number of events in non-overlapping regions of $\mathcal{X}$ are independent. A Poisson process is said to be homogeneous when $\lambda(x)=c$ for some constant $c$. Finally, a marked point process extends the usual definition of a point process to incorporate additional information into each event. A marked point process on $\mathcal{X}$ generates random tuples $x_{s}=(\tilde{x}_{s},m_{s})$, where $\tilde{x}_{s}\in\mathcal{X}$ represents the random location and $m_{s}\in\mathcal{M}$ is the additional “mark” specifying metadata associated with event $s$. See Supplement A for further background details and references.

Neyman-Scott Processes [40] use Poisson processes as building blocks to model clustered data. To simulate a Neyman-Scott process, first sample a set of latent events $\{z_{k}\}_{k=1}^{K}\subset\mathcal{Z}$ from an initial Poisson process with intensity function $\lambda(z):\mathcal{Z}\mapsto\mathbb{R}_{+}$. Thus, the number of latent events is a Poisson-distributed random variable, $K\sim\text{Poisson}(\int_{\mathcal{Z}}\lambda(z)\mathrm{d}z)$. Given the latent events, the observed events $\{x_{s}\}_{s=1}^{S}$ are drawn from a second Poisson process with conditional intensity function ${\lambda(x)=\lambda^{\varnothing}+\sum_{k=1}^{{}_{K}}g(x;z_{k})}$. The nonnegative functions $g(x;z_{k})$ can be thought of as impulse responses—each latent event adds to the rate of observed events, and stochastically generates some number of observed offspring events. Finally, the scalar parameter $\lambda^{\varnothing}>0$ specifies a “background” rate; thus, if $K=0$, the observations follow a homogeneous Poisson process with intensity $\lambda^{\varnothing}$.

A simple example of a Neyman-Scott Process in $\mathbb{R}^{2}$ is shown in fig. 1A. Latent events (purple dots) are first drawn from a homogenous Poisson process and specify the cluster centroids. Each latent event induces an isotropic Gaussian impulse response. Observed events (gray dots) are then sampled from a second Poisson process, whose intensity function is found by summing all impulse responses.

Importantly, we do not observe the number of latent events nor their locations. As described below, we use probabilistic inference to characterize this unobserved structure. A key idea is to attribute each observed event to a latent cause (either one of the latent events or the background process); these attributions are valid due to the additive nature of the latent intensity functions and the superposition principle of Poisson processes (see Supplement A). From this perspective, inferring the set of latent events is similar to inferring the number and location of clusters in a nonparametric mixture model.

We model neural sequences as a Neyman-Scott process with marked events in a model we call PP-Seq. Consider a dataset with $N$ neurons, emitting a total of $S$ spikes over a time interval $[0,T]$. This can be encoded as a set of $S$ marked spike times—the observed events are tuples ${x_{s}=(t_{s},n_{s})}$ specifying the time, $t_{s}\in[0,T]$, and neuron, ${n_{s}\in\{1\,,\,\ldots\,,\,N\}}$, of each spike. Sequences correspond to latent events, which are also tuples $z_{k}=(\tau_{k},r_{k},A_{k})$ specifying the time, $\tau_{k}\in[0,T]$, type, ${r_{k}\in\{1\,,\,\ldots\,,\,R\}}$, and amplitude, $A_{k}>0$ of the sequence. The hyperparameter $R$ specifies the number of recurring sequence types, which is analogous to the number of components in convNMF.

To draw samples from the PP-Seq model, we first sample sequences (i.e., latent events) from a Poisson process with intensity ${\lambda(z)\triangleq\lambda(\tau,r,A)\!=\!\psi\,\pi_{r}\,\mathrm{Ga}(A;\alpha,\beta)}$. Here, $\psi>0$ sets the rate at which sequences occur within the data, $\pi\in\Delta_{R}$ sets the probability of the $R$ sequence types, and $\alpha,\beta$ parameterize a gamma density which models the sequence amplitude, $A$. Note that the number of sequence events is a Poisson-distributed random variable, $K\sim\text{Poisson}(\psi T)$, where the rate parameter $\psi T$ is found by integrating $\lambda(z)$ over all sequence types, amplitudes, and times.

Conditioned on the set of $K$ sequence events, the firing rate of neuron $n$ is given by a sum of nonnegative impulse responses:

We assume these impulse responses vary across neurons and follow a Gaussian form:

where $\mathcal{N}(t\mid\mu,\sigma^{2})$ denotes a Gaussian density. The parameters ${a_{r}=(a_{1r},\ldots,a_{Nr})\in\Delta_{N}}$, ${b_{nr}\in\mathbb{R}}$, and ${c_{nr}\in\mathbb{R}_{+}}$ correspond to the weight, latency, and width, respectively, of neurons’ firing rates in sequences of type $r$. Since the firing rate is a sum of non-negative impulse responses, the superposition principle of Poisson processes (see Supplement A) implies that we can view the data as a union of “background” spikes and “induced” spikes from each sequence, justifying the connection to clustering. The expected number of spikes induced by sequence $k$ is:

and thus we may view $A_{k}$ as the amplitude of sequence event $k$.

Figure 1B schematizes a simple case containing $K=3$ sequence events and $R=2$ sequence types. A complete description of the model’s generative process is provided in Supplement B, but it can be summarized by the graphical model in fig. 1B, where we have global parameters ${\Theta=(\theta_{\varnothing},\{\theta_{r}\}_{r=1}^{R})}$ with ${\theta_{r}=(a_{r},\{b_{nr}\}_{n=1}^{N},\{c_{nr}\}_{n=1}^{N})}$ for each sequence type, and $\theta_{\varnothing}=\{\lambda^{\varnothing}_{n}\}_{n=1}^{N}$ for the background process. We place weak priors on each parameter: the neural response weights $\{a_{nr}\}_{n=1}^{N}$ follow a Dirichlet prior for each sequence type, and $(b_{nr},c_{nr})$ follows a normal-inverse-gamma prior for every neuron and sequence type. The background rate, $\lambda_{n}^{\varnothing}$, follows a gamma prior. We set the sequence event rate, $\psi$, to be a fixed hyperparameter, though this assumption could be relaxed.

Time-warped sequences PP-Seq can be extended to model more diverse sequence patterns by using higher-dimensional marks on the latent sequences. For example, we can model variability in sequence duration by introducing a time warping factor, $\omega_{k}>0$, to each sequence event and changing eq. 2 to,

This has the effect of linearly compressing or stretching each sequence in time (when $\omega_{k}<1$ or $\omega_{k}>1$, respectively). Such time warping is commonly observed in neural data [43, 44], and indeed, hippocampal sequences unfold $\sim$15-20 times faster during replay than during lived experiences [10]. We characterize this model in Supplement E and demonstrate its utility below.

In principle, it is equally possible to incorporate time warping into discrete time convNMF. However, since convNMF involves a dense temporal factor matrix $\mathbf{H}\in\mathbb{R}_{+}^{R\times B}$, the most straightforward extension would be to introduce a time warping factor for each component $r\in\{1,\ldots,R\}$ and each time bin $b\in\{1,\ldots,B\}$. This results in $O(RB)$ new trainable parameters, which poses non-trivial challenges in terms of computational efficiency, overfitting, and human interpretability. In contrast, PP-Seq represents sequence events as a set of $K$ latent events in continuous time. This ultra-sparse representation of sequence events (since $K\ll RB$) naturally lends itself to modeling additional sequence features since this introduces only $O(K)$ new parameters.

We evaluate model performance by computing the log-likelihood assigned to held-out data. Partitioning the data into training and testing sets must be done somewhat carefully—we cannot withhold time intervals completely (as in fig. 2A, top) or else the model will not accurately predict latent sequences occurring in these intervals; likewise, we cannot withhold individual neurons completely (as in fig. 2A, middle) or else the model will not accurately predict the response parameters of those held out cells. Thus, we adopt a “speckled” holdout strategy [54] as diagrammed at the bottom of fig. 2A. We treat held-out spikes as missing data and sample them as part of the MCMC algorithm. (Their conditional distribution is given by the PP-Seq generative model.) This approach involving a speckled holdout pattern and multiple imputation of missing data may be viewed as a continuous time extension of the methods proposed by [27] for convNMF.

Panels B-E in fig. 2 show the results of this cross-validation scheme on a synthetic dataset with $R=2$ sequence types. The predictions of the model in held-out test regions closely match the ground truth—missing spikes are reliably imputed when they are part of a sequence (fig. 2C). Further, the likelihood of the train and test sets improves over the course of MCMC sampling (fig. 2D), and can be used as a metric for model comparison—in agreement with the ground truth, test performance plateaus for models containing greater than $R=2$ sequence types (fig. 2E).

Finally, to be of practical utility, the algorithm needs to run in a reasonable amount of time. Figure 2G shows that our Julia [55] implementation can fit a recording of 120 hippocamapal neurons with hundreds of thousands of spikes in a matter of minutes, on a 2017 MacBook Pro (3.1 GHz Intel Core i7, 4 cores, 16 GB RAM). Run-time grows linearly with the number of spikes, as expected, but even with a single thread it only takes six minutes to perform 1000 Gibbs sweeps on a 15-minute recording with ${\sim\!1.3\times 10^{5}}$ spikes. With parallel MCMC, this laptop performs the same number of sweeps in under two minutes. Our open-source implementation is available at:

https://github.com/lindermanlab/PPSeq.jl.

We first applied PP-Seq to a recording of HVC premotor neurons in a zebra finch,²²2These data are available at http://github.com/FeeLab/seqNMF; originally published in [4]. which generate sequences that are time-locked to syllables in the bird’s courtship song. Figure 3A qualitatively compares the performance of convNMF and PP-Seq. The raw data (top panel) shows no visible spike patterns; however, clear sequences are revealed by sorting the neurons lexographically by preferred sequence type and the temporal offset parameter inferred by PP-Seq. While both models extract similar sequences, PP-Seq provides a finer scale annotation of the final result, providing, for example, attributions at the level of individual spikes to sequences (bottom left of fig. 3A).

Further, PP-Seq can quantify uncertainty in key parameters by considering the full sequence of MCMC samples. Figure 3B summarizes uncertainty in the total number of sequence events, i.e. $K$, over three independent MCMC chains with different random seeds—all chains converge to similar estimates; the uncertainty is largely due to the rapid sequences (in blue) shown in panel A. Figure 3C displays a symmetric matrix where element $(i,j)$ corresponds to the probability that spike $i$ and spike $j$ are attributed to same sequence. Finally, fig. 3D-E shows the amplitude and offset for each neuron’s sequence-evoked response with 95% posterior credible intervals. These results naturally fall out of the probabilistic construction of the PP-Seq model, but have no obvious analogue in convNMF.

Songbird HVC is a specialized circuit that generates unusually clean and easy-to-detect sequences. To compare the robustness of PP-Seq and convNMF under more challenging circumstances, we created a simple synthetic dataset with $R=1$ sequence type and $N=100$ neurons. We varied four parameters to manipulate the difficulty of sequence extraction: the rate of background spikes, $\lambda^{\varnothing}_{n}$ (“additive noise,” fig. 4A), the expected value of sequence amplitudes, $A_{k}$ (“participation noise”, fig. 4B), the expected variance of the Gaussian impulse responses, $c_{nr}$ (“jitter noise”, fig. 4C), and, finally, the maximal time warping coefficient (see eq. 4; fig. 4D). All simulated datasets involved sequences with low spike counts ($\mathbb{E}[A_{k}]<100$ spikes). In this regime, the Poisson likelihood criterion used by PP-Seq is better matched to the statistics of the spike train. Since convNMF optimizes an alternative loss function (squared error instead of Poisson likelihood) we compared the models by their ability to extract the ground truth sequence event times. Using area under reciever operating characteriztic (ROC) curves as a performance metric (see Supplement F.1), we see favorable results for PP-Seq as noise levels are increased.

We demonstrate the abilities of time-warped PP-Seq further in fig. 4E-F. Here, we show a synthetic dataset containing sequences with 9-fold variability in their duration, which is similar to levels observed in some experimental systems [10]. While convNMF fails to reconstruct many of these warped sequences, the PP-Seq model identifies sequences that are closely matched to ground truth.

Finally, we tested PP-Seq and its time warping variant on a hippocampal recording in a rat making repeated runs down a linear track.³³3These data are available at http://crcns.org/data-sets/hc/hc-11; originally published in [56]. This dataset is larger ($T\approx 16$ minutes, $S=\,$137,482) and contains less stereotyped sequences than the songbird data. From prior work [56], we expect to see two sequences with overlapping populations of neurons, corresponding to the two running directions on the track. PP-Seq reveals these expected sequences in an unsupervised manner—i.e. without reference to the rat’s position—as shown in fig. 5A-C.

We performed a large cross-validation sweep over 2,000 random hyperparameter settings for this dataset (see Supplement F.2). This confirmed that models with $R=2$ sequence performed well in terms of heldout performance (fig. 5D). Interestingly, despite variability in running speeds, this same analysis did not show a consistent benefit to including larger time warping factors into the model (fig. 5E). Higher performing models were characterized by larger sequence amplitudes, i.e. larger values of $\mathbb{E}[A_{k}]=\alpha/\beta$, and smaller background rates, i.e. smaller values of $\lambda^{\varnothing}$. Other parameters had less pronounced effects on performance. Overall, these results demonstrate that PP-Seq can be fruitfully applied to large-scale and “messy” neural datasets, that hyperparameters can be tuned by cross-validation, and that the unsupervised learning of neural sequences conforms to existing scientific understanding gained via supervised methods.

A.H.W. received funding support from the National Institutes of Health BRAIN initiative (1F32MH122998-01), and the Wu Tsai Stanford Neurosciences Institute Interdisciplinary Scholar Program. S.W.L. was supported by grants from the Simons Collaboration on the Global Brain (SCGB 697092) and the NIH BRAIN Initiative (U19NS113201 and R01NS113119). We thank the Stanford Research Computing Center for providing computational resources and support that contributed to these research results.

Understanding neural computations in biological systems and ultimately the human brain is a grand and long-term challenge with broad implications for human health and society. The field of neuroscience is still taking early and incremental steps towards this goal. Our work develops a general-purpose, unsupervised method for identifying an important structure—neural sequences—which have been observed in a variety of experimental datasets and have been studied extensively by theorists. This work will serve to advance this growing understanding by providing new analytical tools for neuroscientists. We foresee no immediate impacts, positive or negative, concerning the general public.