Bayesian nonparametric (non-)renewal processes for analyzing neural spike train variability
Authors: David Liu
Máté Lengyel
Venue: NeurIPS 2023 poster (2023)
Keywords: Gaussian processes
renewal processes
point processes
neural data analysis
Bayesian machine learning
non-stationary time series
TL;DR: We propose a flexible modulated point process model for characterising variability in non-stationary neural spike train data

Abstract: Neural spiking activity is generally variable, non-stationary, and exhibits complex dependencies on covariates, such as sensory input or behavior. These dependencies have been proposed to be signatures of specific computations, and so characterizing them with quantitative rigor is critical for understanding neural computations. Approaches based on point processes provide a principled statistical framework for modeling neural spiking activity. However, currently, they only allow the instantaneous mean, but not the instantaneous variability, of responses to depend on covariates. To resolve this limitation, we propose a scalable Bayesian approach generalizing modulated renewal processes using sparse variational Gaussian processes. We leverage pathwise conditioning for computing nonparametric priors over conditional interspike interval distributions and rely on automatic relevance determination to detect lagging interspike interval dependencies beyond renewal order.  After systematically validating our method on synthetic data, we apply it to two foundational datasets of animal navigation: head direction cells in freely moving mice and hippocampal place cells in rats running along a linear track. Our model exhibits competitive or better predictive power compared to state-of-the-art baselines, and outperforms them in terms of capturing interspike interval statistics. These results confirm the importance of modeling covariate-dependent spiking variability, and further analyses of our fitted models reveal rich patterns of variability modulation beyond the temporal resolution of flexible count-based approaches.

Bayesian nonparametric (non-)renewal processes for analyzing neural spike train variability
David Liu Department of Engineering University of Cambridge d1543@cam. ac. uk
Máté Lengyel
Department of Engineering University of Cambridge m. m.lengyel@eng. cam ac
Department of Cognitive Science Central European University
Abstract
Neural spiking activity is generally variable, non-stationary, and exhibits complex dependencies on covariates, such as sensory input or behavior. These dependencies have been proposed to be signatures of specific computations, and characteriz- ing them with quantitative rigor is critical for understanding neural computations. Approaches based point processes provide principled statistical framework for modeling neural spiking activity. However, currently, they only allow the instantaneous mean, but not the instantaneous variability, of responses depend on covariates. To resolve this limitation, we propose a scalable Bayesian approach generalizing modulated renewal processes using sparse variational Gaussian pro- cesses. We leverage pathwise conditioning for computing nonparametric priors over conditional interspike interval distributions and rely on automatic relevance determination to detect lagging interspike interval dependencies beyond renewal order. After systematically validating our method on synthetic data, we apply it to two foundational datasets of animal navigation: head direction cells in freely moving mice and hippocampal place cells rats running along a linear track. Our model exhibits competitive or better predictive power compared to state-of-the-ar baselines, and outperforms them terms of capturing interspike interval statistics. These results confirm the importance of modeling covariate-dependen spiking vari- ability, and further analyses of our fitted models reveal rich patterns variability modulation beyond the temporal resolution of flexible count-based approaches.
1 Introduction
Analyses of spike time [49, 76, 86] and count [12, 26, 57] statistics have revealed neural responses in wivo to structured but generally variable or noisy [24, 44, 70, 78]. To model this stochastic aspect of neural spike trains, probabilistic approaches based on temporal point processes have been widely applicd. This in turn has been a major driver point process theory development [42] for capturing spiking variability structure with statistical models [74]. The study of neural computation underlying naturalistic behavior in particular involves non-stationary spike trains, which presents a significant challenge apparent spiking variability is result of both irreducible "intrinsic" neural stochasticity as well as dependencies on behavioral covariates that can themselves vary on multiple time scales.
Different approaches have been proposed for handling non-stationary spike trains, starting with the classical log Cox Gaussian process [14, 55] to allow variations in the local intensity or firing rate while modeling independent spikes. Dependencies on previous spikes can be captured to first order with renewal processes, and these models have been extended to non-stationary cases through modulation of the hazard function with some time-dependent function [43, 79] or through rescaling interspike intervals with a covariate-dependent rate function [2, 6]. Another approach based on Hawkes processes and spike- history filters [45, 83, 91] introduces conditional point processes that go
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

beyond the first order Markov assumption. Approaches based on recurrent networks [51, 90] and neural ODEs [10, 41] can in theory capture arbitrarily long dependencies on past spikes and input covariates, but provide more limited descriptive interpretabilityy
However, not only the rate but also the variability of spiking encodes task-relevant information [34, 60], and bears signatures of the underlying computations [LI]. Importantly this variability has stimulus- and state-dependent structure [12, 21, 48, 67]. In statistical modeling language, this corresponds to heteroscedastic or input-dependent observation noise. Such structure reflects computations performed in the underlying neural circuit, and thus characterizing it from data in a flexible and robust manner is critical for advancing theories of neural computation. The classical approaches reviewed above do not attempt to characterize such covariate-depp -eeiaatteeeeeddeet changes in variability. Flexible count models have been introduced to more faithfully capture variability at the count level [28], and recent work has extended this to the general case of input-dependent variability [48]. Count approaches however are limited in the resolution of the analysis set by the time bin size. In addition, the resulting count statisties strongly depend on the chosen time bin size [48, 73].
While firing rates are routinely modeled as input-dependent, extending point process models with input-dependent variability has not been widely explored in the literature. Rate-rescaled and mod- ulated renewal processes rely on fixed base renewal densities. Allowing the shape parameters of the renewal density to vary with covariates in the corresponding hazard function is one potential approach, but this still relies on a commitment to particular parametric family of renewal densities. -pike-history filters in conditional point processes are conventionally fixed and thus do not directly model input-dependent spiking variability, though dependence on observed and unobserved covariates [91] and switching filters based on discrete states [23] have been considered. Recent work has moved away from parametric filters to nonparametric Gaussian processes [18], which can be extended to flexibly model dynamic filters as functions of external covariates with a spatio-temporal Gaussian process. However, any modulation of the filter will no longer permit fast convolutions, and such models will be computationally expensive as the filter needs to be recomputed every time step. The direct nonparametric estimation of conditional intensity functions based on maximum likelihood has been explored [13, 82], but scalable Bayesian approaches have remained absent.
Contribution To enable flexible modeling as well as modulation of instantancous point process statistics for analyzing neural spike train variability, we introduce the Bayesian nonparametric non- renewal process (NPNR). NPNR builds on sparse variational Gaussian processes and defines a nonparametric prior over conditional interspike interval distributions, generalizing modulated renewal processes with nonparametric renewal densities and spike-history dependencies beyond renewal order. In particular, our point process can flexibly model modulations of not only spiking intensity but also variability. We validate our model using parametric inhomogeneous renewal processes, recovering conditional interspike interval distributions and identifying renewal order in spike- dependence. On neural data from mouse thalamus and rat hippocampus, our method has competitive predictive power while being superior in capturing interspike interval statisties from the non-stationary data. In particular, our method provides instantaneous measures of spike train variability that are modulated by covariates, and shows rich variability patterns in both datasets consistent with previous studies at coarser timescales. We provide JAX [4] implementation of our method as well as established baseline models within a scalable general variational inference scheme."
2 Background
We start with a brief overview of the theoretical foundations and related point process models, as well as their combination with Gaussian processes introduce non-stationarity.
2.1 Temporal point processes
Statistical modeling of events that occur stochastically in time is handled by the general framework of temporal point processes [59, 72]. Denoting the number of events that occurred until time t by N (f). a temporal point process model is completely characterized by its conditional intensity function (CIF)
H() = lim [N + öt) - [H.
5t-+0
8t
'Code available https::/ http://hhhhh...ceoearrradaddd//s hww.:parrrrratttieooooneeeeaaala process

---

where A(t)ot is the probability to emit a spike event in the infinitesimal interval [t. t 6t) conditioned on H4, the spiking history before We can write the point process likelihood for a single neuron spike train consisting of an ordered sequence of S spike events at times 1 as [5]
(slx()) = II A(4 ((AAFH.)) E A(e
In neuroscience applications, one often wants to describe modulation of the point process statistics with some time-varying covariates æ=(), such animal head direction body position, which leads to generalized CIF A(t)H,, ægg)). Several classes of models have been proposed that are defined by particular restrictions on the functional form of A(t)HL, z(..).
2.1.1 Inhomogeneous renewal processes
The statistical model in Eq. 2 describes dependencies between all spikes. One common simplification the renewal assumption: interspike intervals (ISIs) A = 1++1 t; drawn i.i.d. from an interval distribution called the renewal density g( g(A;6) with parameters A. This induces a Markov structure
Renewal assumption
((44)4-1) in the spike train likelihood
S-1
p(1 p(::...,(s;;)
1s; 0) = p(2))) ....:tstssss)s)) = p(t) (u;;
Common renewal densities used for neural data are the exponential (equivalent to a Poisson process), gamma, and inverse Gaussian distributions [2, 6].
Hazard function modulation Non-stationary point processes need to model changes in statistics with time, and combined with the renewal assumption obtains inhomogeneous renewal processes. classical approach that dates back to Cox [14] is to modulate the hazard function (Appendix A)
((t h(F) p(t)
with time since last spike 5 t t, and the modulation factor p(t). In context, this can be replaced by some function covariates p(x;) [43]. A multiplicative interaction between o and h as above is typically considered, though this framework allows general parametric forms [71, 79].
Rate-rescaling Another approach that has been widely applied the neuroscience community is rate-rescaling [2, 6]. closely related to time-rescaling (see Appendix B.3). Here, modulation is achieved with a rate function r(æ,) 2 0 that transforms time t into rescaled time
f r r(zv)) dt'
E(t) f
(5)
This maps all spike times 4; to rescaled times E, and will be one-to-one as long as r(t) > 0. By modeling the rescaled ISIs (+) = 4++1 - 1, as drawn from stationary renewal density g(.), we obtain inhomogeneous renewal process from a homogeneous one. The CIF becomes dependent on the covariate path since last spike Pe (æ(n)|u E (f). see Appendix B.3.
2.1.2 Conditional point processes
Conditional Poisson processes The renewal assumption ignores correlations between ISIs, which generally are observed in both the peripheral and central nervous system [1, 25] and can be com- putationally relevant for signal detection and encoding [8. 9, 22, 69]. Going beyond Markovian dependencies, a tractable approach, similar to Hawkes processes [51], is to introduce causal linear filter (t) that is convolved with spikes and added to the log CIF, giving conditional Poisson processes
log (t)),, xx) h y(t) + f(xx), with y(t)
Σ 8(t t:)
(6)
where denotes temporal convolution. These models are closely linked to mechanistic integrate-and- fire models [52, 85) and have a long history the neuroscience literature [18, 31, 35, 47, 62, 66, 83), appearing generalized linear models (GLMs) and spike response models (SRMs).

---

Conditional renewal processes An even more expressive model can be obtained by replacing the Poisson spiking process with rese-lescaled renewal process. This results conditional renewal process, where the rate function has history dependence r(4 (H,, 2tx) as in Eq. 6.
2.2 Gaussian process modulated point processes
Gaussian processes (GPs) represent data-efficient alternative neural networks, which have been widely used to model the CIF [59, 72, 90]. When combining GPs with point process likelihoods, the resulting generative model leads to doubly stochastic processes for event data. Placing a Gaussian process prior [87] over the log intensity function leads to the classic log Cox Gaussian processes [14, 55], and in the same spirit one can modulate renewal hazard functions [79] or perform rate- rescaling [16] with GPs. Such constructions form the basis of many widely used Bayesian neural encoding models for spike trains, both for modeling single neuron responses [15, 16, 68] as well as population activity [19, 93]. Combining the flexibility offered by renewal and conditional point processes with GP rate or modulation functions within variational framework has been impeded by the fact that the original papers were built on a GLM framework with parametric covariate mappings [6. 7, 18, 65]. To provide fair comparison of NPNR to these baselines, we implement a general scalable variational inference framework for the construction and application of such models (see Appendix B for details on the baseline models).
Method
We now introduce the nonparametric non-renewal (NPNR) process and present the approximate Bayesian inference scheme used for model fitting, noting connections to related works in the literature. Our NPNR model provides nonparametric generalization of modulated renewal processes heyond renewal order, and adds suitable inductive biases for neural spike train data. It implicitly defines a flexible prior over conditional ISI distributions that can be computed using pathwise conditioning, which enables one to analyze spiking variability modulation with minimal parametric constraints. Furthermore, the Bayesian framework provides an elegant data-driven approach inferring the lagging ISI order of the spike-histon dependence.
3.1 Generative model
Conditional intensity surface priors directly model the CIP, or more precisely its logarithm, of the form
To obtain flexible modulated point process models,
Xt) zg:) = A(t)t,, f...........
(7)
where t, is the most recent spike at current time t. First considering the renewal case, note the spatio-temporal structure in the log CIF using time since last spike T = t - t4 is
log A(t)He, a..) log A(t|44, 2x) = f 2x))
(8)
which suggests placing spatio-tempont GP prior the log CIF to describe a log intensity surface f(r, ~ gP (m(r,)), ku(T, (x,
This generalizes the parametric forns of modulation considered in previous approaches [43, 79], in particular allowing modulation of the effective instantaneous renewal density by covariates Ty. We can introduce lagging ISIs covariates (t) with lag ls as depicted in Fig. 1A to extend the model to non-renewal process
log A(t|te, 44-1, 44-1,44222211.. f(r, (T, A1 (t), At(),....)
(10)
and for maximum ISI lag we denote the lagging ISIs A, (t)(),..... AK(t) to obtain
log X(1 = f(T, A...,)
(11)
Inductive biases for neural data for building more realistic prior. a spike, though in practice neural recordings may not respect this due to contamination in spike sorting [38]. Another potentially useful inductive bias is that changes in the spiking intensity of
For neural spiking data, there are biological properties to consider Firstly, neurons have refractory periods immediately following

---

neurons fluctuate mostly at shorter ISI timescales [32], whereas at longer delays they tend to be temporally smoother. The latter suggests non-stationary GP kernels to be more suitable for modeling the spike-history dependencies [18]. However, non-stationary kernels do not allow straight-forward use of random Fourier features for evaluating GP posterior function samples at many locations with pathwise conditioning [88, 89]. This in particular is needed to compute the conditional ISI distributions _(.]...) in Eq. 15, see Appendix B,5 for details. To achieve the desired non-stationarity for modeling Eq. 11 while maintaining the ability to draw samples using pathwise conditioning, we apply time warping on t from [0, oo) + [0, 1] with some warping timescale Tw
-=
- log (1 7)
(12)
and place stationary Gaussian process prior over the warped temporal dimension f(f....) GP
with a temporal kernel k(F, F') = k(f- This transformation is monotonic (see Fig. IB), and
hence we can easily compute the transformation on the CIF
di
= ef (-1(1)/ dr
- x(7 -
(13)
Tw
Similarly, we apply time warping to the A& dimensions which we also place stationary kernels A(,, Á') = k(A Á'). We note that unlike spike-history filters in conditional point processes (Eq. 6) which do not change with inputs, the resulting coupling to past activity Eq. 11 dependent on covariates æt. This allows one to capture spiking variability modulation via the conditional ISI distribution perspective discussed below. The refractory nature of real neurons can be addressed by the mean function
m(f, a) m(() am - f//mm + bm
(14)
with parameters am, Tm,y bm. where refractory periods can be modeled with large negative am.
Conditional ISI distributions Instead of looking the CIF, can view the model as prior over conditional ISI distributions as depicted in Fig. IC using the relation (see Appendix A)
g A. F(r,)) x A(t)PLe. 2x) 11 X((IH,
(15)
where one drops the dependence on Hi in for the modulated renewal case. If one fixes the lagging ISIs and picks a constant covariate path A.,... this can be interpreted as an instantaneous ISI distribution of neuron at the conditioned inputs A. and æ.. Moments of the conditional ISI distribution are computed using Leguss- quadratures in warped time (Appendix B.5)
g(T) TThh dr = L.'' dr di
E [T(+)(T"]
((T(F)) di
(16)
and this can be used to compute tuning curves of spike train statistics, such as the mean ISI E|r] and coefficient of variation CV / Var [r]] [58], as a function of x.. This approach generalizes the homogencous case considered in the literature, and in particular allows one compute instantaneous measures of non-stationary spike train variability that are otherwise non-trivial to estimate [58, 75].
3.2 Inference
Temporal discretization The generative model is formulated as continuous time model. In practice, neural and behavioural data are typically recorded with finite temporal resolution at small regular intervals At. The cumulative intensity integral has to be approximated by a sum, though note that directly modeling the cumulative hazard function [59] elegantly avoids this for purely temporal point processes. Spike times are now discretized as a binary vector y = .......99 .r] where y/ = 1 if I has spike event, zero otherwise. Overall, this discretizes the point process likelihood Eq. 2 as
p(y II th ΣT ΣT, A.
p(t . tts.....ssl(((t))
(17)
where we have T time steps in total. Note that the discretization scheme implies we do not have observations at 0, since the time step immediately after observed spike has = At.
5

---

A
B
C
non-refractory prior D m 0)
refractory prior (3m -6)
0.3
Snseea time
log CIF
density
4
= (191)
(a.u.)
distribution
f(au.
X3
WWW
time
time
density
(8)
Figure I: Schematic of our proposed model. (A) Time since last spike 7 and lagging ISIs for an observed spike train (top row) alongside covariates æ. (B) Illustration of the time warping procedure. We fix the warping parameter Tw to the empirical mean ISI, which leads to more uniform distributions i E [0, 1] suitable for stationary GP kernel. (C) Prior samples from the generative model for two values am characterized by the lack and presence of a refractory period. The transformation Eq. 15 links the log CIF (top rows) with conditional ISI distributions (bottom rows).
Variational lower bound We use stochastic variational inference [39] with batches obtained from consecutive temporal segments and sparse variational GPs [37], giving the loss objective
L=
By Bg(salum)) (u)) -Wneenn + At
Σ for + DKl(g(un)) p( p(un))
(18)
1=1
where n indexes neurons (of which there are N)², points, p(u) the GP prior at inducing locations, ((us the variational posterior, and q( ()) the conditional posterior (see Appendix B. for details). Combined with temporal mini-batching to fit batch segments of length T, we can fit to very long time series given the O(N M + N M3) computational complexity. We also no longer rely on computing hazard functions of parametric renewal densities to obtain the CIF, which can be numerically unstable. Modulated renewal processes instead rely on a specialized thinning procedure [79], but we take a more scalable and general variational approach. Overall, we optimize the kernel hyperparameters, variational posterior mean and covariance, inducing point locations, and mean function parameters am, bm and Tm using gradient descent with Adam [46] (see Appendix C for details). The time warping parameter Tw is fixed in our experiments to the empirical mean ISI, and the hyperparameter K is fixed and chosen in advance (see also subsection automatic relevance determination below).
[Ui...., UM]] denotes the set of M inducing
Automatic relevance determination The Bayesian framework with GPs enables us to perform automatic relevance determination (ARD) over the input dimensions to automatically select relevant input [37, 80]. Applied to lagging ISI dimensions in our NPNR model, this provides an elegant approach making a data-driven renewal assumption and generally determining the spike-history dependence of the CIF. We choose to fix Tu to the empirical mean ISI as shown in Fig. 1B (rather than learning to achieve interpretability of kernel timescales for ARD (Fig. 6A) at small cost of performance (Fig. 12). For a chosen maximum lag K, there is no need for manual selection of the history interaction window for GLM spike-history filters, though recent work on nonparametric GLM filters provides a related window size selection procedure [18]. In the spirit of Bayesian models, we choose K give sufficiently high capacity model [40, 80] to be able to flexibly capture history dependence, as seen in panel D of Fig. 3 and Fig. 4.
4 Results
All datasets discretize spike trains and input time series at regular intervals of A/ = ms. We use a product kernel for k( x²) with periodic kernels for angular dimensions, and squared exponential kernels in other cases. For k(f, (") and ke A Á'), pick product kernel with Matérn-3/2 (see Fig. 12 for different kernel choices) and set the maximum ISI lag K = 3. For illustration, conditional
Note the simple summation over n in Eq. 18, as our model does not capture neural correlations without introducing latent covariates [48]. In other words, in its current form, our model treats neurons independent.

---

A
maps
B
inferred
Poisson
conditional Poisson
kemal timescales
G
frequency
onngn omegince de ollo fun ursion
10
10"
1234
neuron index
order
rescaled gamma
NPNR (ours)
true renewal shapes
log
density
probabiltt
log
ISI (a.u.)
(s)
quantle
Figure 2:
Validation on synthetic data. (A) True rate maps (brighter is higher) defined over square environment (top) and base renewal densities (bottom) each column given by gamma (left), log normal (middle) and inverse Gaussian (right) distributions with various shape parameters. Each color corresponds
separate neuron. (B) Posterior mean rate maps (top left) and conditional ISI distribution samples in gray overlaid on true renewal densities at various locations (bottom) for the NPNR process fit to synthetic data. The relevance boundary (dotted line) for kernel timescales (top right) is placed at i = 3 (dimensionless). (C) QQ-plots of fitted models (each curve is a neuron).
ISI distributions and corresponding tuning curves computed by fixing A to be the mean ISI per neuron. Firing rates are defined as 1 1/E[r]] since this corresponds to the number of spikes fired per unit time in infinitely large time bins for renewal process (Eq. 26). GP inducing points were randomly initialized, and for fair comparison, all models used inducing points for each covariate dimension (including temporal dimensions and A in the NPNR process). For each experiment, we repeat model fitting with 3 different random seeds and pick the model with the best training likelihood. Further details experiments presented in Appendix C.
4.1 Validation on synthetic data
For validating our approach, we generate 1000 of data using rate-rescaling [65] mimicking place cell population of 9 neurons for an animal moving in a 2D square arena, each with a unique rate map and renewal density (Fig. 2A, details in Appendix C). The models applied are baseline Poisson, raised cosine filter conditional Poisson and rate-rescaled gamma processes (Appendix C), and our NPNR process. Note that the rescaled gamma process is within-model class for 3 of the synthetic neurons. Inferred conditional ISI distributions and rate maps of our NPNR process in Fig. 2B are close to ground truth, showing the ability of our model to capture modulated spiking statistics drawn from various parametric families. To assess how well ISI statistics are captured, we apply time-rescaling using the GP posterior mean functions (Appendix A) which visualize with coantile-quantile (QQ) plots [6] in Fig. 2C. Again, we see an excellent fit our model compared to baseline models, indicating that only the NPNR is capable of satisfactorily capturing the empirical ISI statistics. Learned temporal kernel timescales of the NPNR process in Fig. 2B show a clear separation between the time since last spike dimension (lag 0) and lagging ISI A dimensions (lag 2 1) with the dotted relevance boundary at = 3 (dimensionless), as expected for renewal processes.
4.2 Neural data
Now we apply our method to head direction cells freely moving mice [63, 64] and place cells in rats rumning along linear track [54]. We select 33 units from the mouse and 35 units from the rat data, which leads to around and 68 million data points to the training set, respectively (see Appendix C for preprocessing details and Appendix B.2 on data scaling). Experiments involve fitting to the first half of dataset (~~18 min. for mouse, 32 min. for rat), and testing on the second half split into 5 consecutive segments. The split into 5 test folds is used to guantify dataset variability. For
7

---

A
(ours)
EDF
RC cond.
Qquantie
vials leg CIF O
cond.
Gauss.
Polason
10th
10-1
300
400
KS p-values
test ELL (nats/s)
D
10
nouron
kernel timescales
n0"
lag
CV rale (H22
frequency
average CV
R
probabiitty
leg
P.TEOOOOEEEEEE neuron Index
density
1234
average rate (Hiz)
head direction
2n
ISI
Figure 3: Application to mouse head direction cell data. (A) Violin plot of KS p-values per neuron (left, lines marking quartiles) and test expected log likelihoods with errorbars showing s.c.m. across test folds (right). Larger values in both metrics indicate better model fit to data. (B) QQ-plots for various models (each curve a neuron) identified by color (panel A left). (C) Predicted log CIF (middle) for observed spike train (top) and posterior spike train samples (bottom) conditioned on the same covariates 2; for various models identified by color. (D) Left: temporal kernel timescales for r (lag 0) and A₁ dimensions with the relevance boundary at / 3 (dimensionless)) Right: histogram of "[SI-order" (1 + largest lag & for which k is below the boundary) across neurons. (E) Time average of estimated instantaneous rates and CVs from the training data (left) and R2 values of CV-rate regression with a linear and GP model (right). (F) Posterior median and 95% intervals of tuning curves over head direction for the rate and CV, with posterior ISI distribution samples (right) at dashed locations.
prediction, we evaluate the expected log likelihood summed over all neurons
ELL Bog
(19)
using Gauss- Herlmie guadrature with 50 points (Monte Carlo for renewal processes, see Appendix B). To assess goodness-of-fit to ISI statistics, we compute QQ plots as before and apply the Kolmogorov- Smirnov (KS) test, giving p-value per neuron that indicates how likely the data came from the model (Appendix A). Baselines are the inhomogeneous Poisson (P), rate-rescaled gamma (G) and inverse Gaussian (IG) renewal [6, 83], raised cosine (RC) filter conditional Poisson [66, 85] and renewal [65], and nonparametric (NP) filter conditional Poisson processes [18] (details in Appendix C).
4.2.1 Mouse head direction cell data
We choose the animal head direction as our ID input covariate x. From the KS test p-value distribution, we see that our model outperforms all basclines capturing ISI statistics (Fig. 3A left, higher is better; Fig. 3B, QQ plots closer to diagonal). It performs competitively to conditional Poisson processes in terms of predictive performance (Fig. 3A right), but those models fail capture ISI statisties. In addition, we note the spiking saturation in some samples of the nonparametric conditional Poisson model (Fig. 3C, purple) due known instability [3, 29]. Samples from our model (Fig. 3C, gray) exhibit visually similar spike patterns to the real spike train segment. Furthermore, kernel timescales Fig. 3D show a sizable fraction of the population is characterized by a non-renewal spiking process.
Neural dispersion regimes From Fig. 3E and we observe both under- and overdispersion (CV smaller and bigger than one) consistent with previous study based on spike counts [48]. Estimated instantaneous rates and CVs Fig. 3E are computed using the conditional ISI distribution evaluated along the time series of covariates in the training data. One can regress instantaneous CV against rate, and the CV-rate R² shows some cells with near linear relations and some with nonlinear trends that can be captured by GP. Despite that, many cells still show low overall R², implying there is
8

---

A
(ours)
EDF
RC
quantie
vials leg CIF
Gauss.
Nww Mww
Polason
10Th
10-1
-20
-100
KS p-values
tast ELL (nata/s)
D
neuron
16
kernel timescales
log
37.11
fraquenny
average CV
R
LFP L-P--phase
leg
2.2 0.4
4.6 0.2
max
crmmermonnnoyyo nouron Index
dersity
average rate (Hiz)
X posiioon
ISI (8)
Figure 4:
Application to rat hippocampal place cell data. (A)-(E) Similar to Fig. 3A-E. (F)
Posterior mean tuning maps over z- c--position and 9-phase for the rate and CV (left) for left-to-right runs (based on head direction) with posterior ISI distribution samples (right) at marked locations.
parametric relation. Most neurons increase CV with rate, while some show slight linear decrease reminiscent of refractory Poisson processes (Eq. 31, see Fig. 7 for CV-rate patterns).
4.2.2 Rat place cell data
In this case, æ is 3D consisting of the body position along the track, head direction, and local field potential (LFP) 8-phase. Our model significantly outperforms all baselines, having both better KS test p-value distribution (Fig. 4A left, higher is better; Fig. 4B, QQ plots closer to diagonal) and predictive performance (Fig. 4A. right). Note that the ELLs for this dataset differ from those shown in Fig. 3A due to the fundamentally different (less predictable) spiking statistics of place cells compared head direction cells (e.g. due theta oscillations). As the nonparametric conditional Poisson model introduces nonparametric but covariate independent variability patterns, these results highlights the importance of modeling coveriate-depennent spiking variability. Note that the rate-rescaled renewal processes struggle to fit this data, with test ELLs of inverse Gaussian models below -200 nats/s (Fig. 4A right). Samples from our model (Fig. 4C, gray) show it captures the characteristic bursting nature of the real spike train segment. Kernel timescales in Fig. 4D show most cells are described well by renewal process, different to mouse data Fig. 4D.
covariate-independen
Capturing overdispersion We CV values in Fig. 4E higher than the mouse thalamus dataset, consistent with overdispersion of place cell discharge in 2D open field navigation [26]. Similar to the mouse data, the CV-rate R² again shows there is generally observe larger increases in CV with firing rate compared to mouse data (Fig. 9).
parametric relation. We also tend to
0-modulation and phase precession Spiking activity modulation during 0-cycles [53] prominent in rat hippocampus, and is visible here in the log CIF (Fig. 4C). We also see phase precession [76] in Fig. 4F (top), classical example where spike timing relative some rhythm has coding significance [33]. Our model enables one to extract not only spiking intensity but also variability, and shows that variability generally inherits the phase precession pattern (Fig. 4F bottom).
Discussion
5.1 Limitations and further work
ISI statistics Apart from the coefficient of variation, there are other ISI statistics that characterize spiking dynamics aspects such as bursting regularity. Of particular interest the local coefficient of variation [75], which involves joint statistics of consecutive ISIs (A (9), A (I-11) that can be computed

---

from our model (see Appendix D). The same applies to serial correlations [25], which may provide insights into biophysical details [77]. Furthermore, quantifying the shape of ISI distributions is of interest as it is associated with various properties of the underlying neural circuit dynamics [61].
Neural correlations To capture correlations in multivariate spike train data, direct spike couplings in GLMs are less suitable for current neural recordings compared to latent variable models due to the sparse sampling of populations by electrodes [50]. Combining the latter alongside observed covariates [48] with our point process provides a powerful framework for capturing correlations [81], which can have significant impact on neural coding [56]. To perform goodness-of-He tests, the Kolmogorov-Sming Smirnov test with time-rescaling can be extended the multivariate case [30, 92].
5.2 Conclusion and impact
We introduced the Bayesian nonparametric non-renewal (NPNR) process for flexible modeling of variability in neural spike train data. On synthetic renewal process data, NPNR successfully captures spiking statistics and their modulation by covariates, and finds renewal order in the spike-history dependence. When applied to mouse head direction cells and rat hippocampal place cells, NPNR has competitive or improved predictive performance to established baseline models, and is superior in terms of capturing ISI statistics, establishing the importance of capturing saariateedeeenden variability. NPNR-based analyses recover known behavioral tuning, while also revealing novel patterns of spiking variability at millisecond timescales that are compatible with count-based studies.
Neural firing rates traditionally characterize most computational functions and information encoded by neurons [16, 17, 33], but recent work on VI [20, 27, 36, 60) and hippocampal place cells [84] have started to assign computationally well-defined roles to variability in the context of representing uncertainty. Our method introduced in this paper is a principled tool for empirically characterizing neural spiking variability and its modulation at the timescales of individual spikes, and we hope model will be useful for revealing new aspects of neural coding. Such findings are foundational to advances in computational and theoretical neuroscience, and may have downstream practical applications in designing and improving algorithms for brain-machine interfaces.
Acknowledgments and Disclosure of Funding
This work was supported by the Cambridge European and Wolfson College Scholarship by the Cambridge Trust (D.L.) and by the Wellcome Trust (Investigator Award in Science 21 212262/7/18/7 to M.L.). We are grateful to Kristopher Jensen, Marine Schimel and Valentina Njaradi for helpful comments on the manuscript. We would also like to thank Alexander Terenin and Jonathan So for helpful discussions.
References
[1] Avila-Akerberg, O. and Chacron, M. J. (2011). Nonrenewal spike train statistics: causes and functional consequences on neural coding. Experimental brain research, 210(3))353-371.
[2] Barbieri, R., Quirk, M. C., Frank, L. M., Wilson, M. A., and Brown, E. N. (2001). Construction and analysis of non- non-poisson stimulus-response models neural spiking activity. Journal of neuroscience methods, 105(5):2
[3] Bellec, G., Wang, S., Modirshanechi, A., Brea, J., and Gerstner, W. (2021). Fitting summary statisties of neural data with a differentiable spiking network simulator. Advances in Neural Information Processing Systems. 14:18855--188
[4] Bradbury, J., Frostig, R., Hawkins, P., Johnson, M. J., Leary, C., Maclaurin, D., Necula, G., Paszke, A., VanderPlas. J., Wanderman-Miin S., and Zhang, Q. (2018). JAX: composable transformations of Python+NumPP programs.
[5] Brown, E. N., Barbieri, R., Eden, U. T., and Frank, L. M. (2003). Likelihood methods for neural spike train data analysis. Computational neuroscience: comprehensive approach, pages 253-286.
10

---

[6] Brown, E. N., Barbieri, R.. Ventura, V., Kass, R. E., and Frank, L. M. (2002). The time-rescaling theorem and its application to neural spike train data analysis. Neural computation, 14(2):325-346. -366.
[7] Brown, E. N., Frank, L. M., Tang, D., Quirk, M. C., and Wilson, M. A. (1998). A statistical paradigm for neural spike train decoding applied to position prediction from ensemble firing patterns of rat hippocampal place cells. Journal of Neuroscience, 18(18):7411-745
[8] Chacron, M. J., Lindner, B., and Longtin, A. (2004). Noise shaping by interval correlations increases information transfer. Physical review letters, 92(8):080601.
[9] Chacron, M. J., Longtin, A., and Maler, L. (2001). Negative interspike interval correlations increase the neuronal capacity for encoding time-dependent stimuli. Journal of Neuroscience, 21(14)5328-5344
[10] Chen, T. Q., Rubanova, Y., Bettencourt, J., and Duvenaud, D. K. (2018). Neural ordinary differential equations. In Advances neural information processing systems, pages 6571 6551-6
[11] Churchland, A. K., Kiani, R., Chaudhuri, R., Wang, X.-J,, Pouget, A., and Shadlen, M. N. (2011). Variance signature of neural computations during decision making. Neuron, 69(4))818-83)
[12) Churchland, M. M., Byron, M. Y., Cunningham, J. P., Sugrue, L. P., Cohen, M. R., Corrado, G. S., Newsome, W. T., Clark, A. M., Hosseini, P., Scott, B. B., et al. (2010). Stimulus onset quenches neural variability: widespread cortical phenomenon. Nature neuroscience, 13(3):369.
[13] Coleman, T. P. and Sarma, S. S. (2010). A computationally efficient method for nonparametric modeling neural spiking activity with point processes. Neural Computation, 22(8):0002-003)
[14] Cox, D. R. (1972). The statistical analysis of dependencies in point processes. Stochastic Point Processes. Wiley: New York, pages 55-66.
[15] Cunningham, J. P., Shenoy, K. V., and Sahani, M. (2008). Fast gaussian process methods for point process intensity estimation. In Proceedings of the 25th international conference on Machine learning, pages 192-199.
[16] Cunningham, J. P., Yu, B. M., Shenoy, K. V., and Sahani, M. (2007). Inferring neural firing rates from spike trains using gaussian processes. Advances in neural information processing systems, 20.
[17] DePasquale, B., Sussillo, D., Abbott, L., and Churchland, M. M. (2023). The centrality of population-level factors to network computation is demonstrated by a versatile approach for training spiking networks. Neuron.
[18] Dowling, M., Zhao, Y., and Park, I. M. (2020). Non-parametric generalized linear model. arXiv preprint arXiv:2009.00 016.
[19] Duncker, L. and Sahani, M. (2018). Temporal alignment and latent gaussian process factor inference in population spike trains. Advances in neural information processing systems, 31.
[20] Echeveste, R., Aitchison, I.., Hennequin, G., and Lengyel, M. (2020). Cortical-like dynamics in recurrent circuits optimized for sampling-bases probabilistie inference. bioRxiv, page 696088.
[21] Ecker, A. S., Berens, P., Cotton, R. J., Subramaniyan, M., Denfield, G. H., Cadwell, C. R., Smirnakis, S. M., Bethge, M., and Tolias, A. S. (2014). State dependence of noise correlations in macaque primary visual cortex. Neuron, (2)))2335-
[22] Engel, T. A., Helbig, B., Russell, D. F., Schimansky- -eier, L., and Neiman, A. B. (2009). Coherent stochastic oscillations enhance signal detection in spiking neurons. Physical Review E, 80(2):22111..
[23] Escola, S., Fontanini, A., Katz, D., and Paninski, (2011). Hidden markov models for the seimuuus relationships of multistate neural systems. Neural computation, 1375) 1132.
[24] Faisal, A. A., Selen, I.. P., and Wolpert, D. M. (2008). Noise in the nervous system. Nature reviews neuroscience, 9(4):292-303.

---

[25] Farkhooi, F., Strube-Bloss, M. F., and Nawrot, M. P. (2009). Serial correlation in neural spike trains: Experimental evidence, stochastic modeling, and single neuron variability. Physical Review E, 79(2):221905
[26] Fenton, A. A. and Muller, R. U. (1998). Place cell discharge is extremely variable during individual passes of the rat through the firing field. Proceedings of the National Academy of Sciences,
[27] Fiser, J., Berkes, P., Orbán, G., and Lengyel, M. (2010). Statistically optimal perception and learning: from behavior to neural representations. Trends in cognitive sciences, 14(3):119-133
[28] Gao, Y., Busing, I.., Shenoy, K. V., and Cunningham, P. (2015). High-dimensional neural spike train analysis with generalized count linear dynamical systems. In Advances in neural information processing systems, pages 2044-2052.
[29] Gerhard, F., Deger, M., and Truccolo, W. (2017). On the stability and dynamics stochastic spiking neuron models: Nonlinear hawkes process and point process glms. PLoS computational biology, 1305)3e0
[30] Gerhard, F., Haslinger, R., and Pipa. G. (2011). Applying the multivariate time-rescaling theorem to neural population models. Neural computation, 23(6):1452-445
[31] Gerstner, W. (2001). A framework for spiking neuron models: The spike response model. In Handbook of biological physies, volume 4, pages 469-516. Elsevier.
[32] Gerstner, W., Kistler, W. M., Naud, R., and Paninski, L. (2014). Neuronal dynamics: From single neurons to networks and models of cognition. Cambridge University Press.
[33] Gerstner, W., Kreiter, A. K., Markram, H., and Herz, A. V. (1997). Neural codes: firing rates and beyond. Proceedings of the National Academy of Sciences, 17404111740-122411
[34] Ghanbari, A., Lee, C. M., Read, H. L., and Stevenson, H. (2019). Modeling stimulus- dependent variability improves decoding of population neural responses. Journal of Neural Engineering, 16(6):066018.
[35] Harris, K. D., Henze, D. A., Hirase, H., Leinekugel, X., Dragoi, G., Czurkó, A., and Buzsáki, G. (2002). Spike train dynamics predicts theta-related phase precession in hippocampal pyramidal cells. Nature, 41766890):78--74.
[36] Hennequin, G., Ahmadian, Y., Rubin, D. B., Lengyel, M., and Miller, K. D. (2018). The dynamical regime of sensory cortex: stable dynamics around single stimulus-tuned attractor account for patterns of noise variability. Neuron, 98(4):846-866 860.
[37] Hensman, J., Fusi, N., and Lawrence, N. D. (2013). Gaussian processes for big data. In Nichol- son, A. E. and Smyth, P., editors, Proceedings of the Twenty-Ninth Conference on Uncertainty in Artificial Intelligence, UAI 2013, Bellevue, WA, USA, August l 11-15, 2013. AUAI Press.
[38] Hill, D. N., Mehta, S. B., and Kleinfeld, D. (2011). Quality metrics to accompany spike sorting of extracellular signals. Journal of Neuroscience, 31924):8699-8755
[39] Hoffman, M. D., Blei, D. M., Wang, C., and Paisley, J. (2013). Stochastic variational inference. The Journal of Machine Learning Research, 14(1::130311447
[40] Jensen, K., Kao, T.-C., Stone, J., and Hennequin, G. (2021). Scalable bayesian gpfa with automatic relevance determination and discrete noise models. Advances in Neural Information Processing Systems, 34: 10613-10666.
[41] Jia, J. and Benson, A. R. (2019). Neural jump stochastic differential equations. Advances in Neural Information Processing Systems, 32.
[42] Kass, R. E., Kelly, R. C., and Loh, W.-L. (2011). Assessment of synchrony in multiple neural spike trains using loglinear point process models. The annals of applied statistics, 5(2B):1262.
[43] Kass, R. F. and Ventura, V. (2001). A spike-train probability model. Neural computation, 1713-177.
12

---

[44] Kayser, C., Logothetis, N. K., and Panzeri, S. (2010). Millisecond encoding precision of auditory cortex neurons. Proceedings of the National Academy of Sciences, 107(39):16976-1699
[45] Kim, S., Putrino, D., Ghosh, S., and Brown, E. N. (2011), A granger causality measure for point process models of ensemble neural spiking activity. PLoS computational biology, 7(3)=1001110 1100
[46] Kingma, D. and Ba, J. (2014), Adam: A method for stochastic optimization. arXiv preprint arXiv: 11Xiv:141..6980.
[47] Linderman, S. W. Adams, R. P., and Pillow, J. W. (2015). Inferring structured connectivity from spike trains under negative-binomia generalized linear models. Cosyne Abstracts.
[48] Liu, D. and Lengyel, M. (2021). A universal probabilistic spike count model reveals ongoing modulation of neural variability. Advances in Neural Information Processing Systems, 34:13392- 13405.
[49] London, M., Roth, A., Beeren, L., Häusser, M., and Latham, P. E. (2010). Sensitivity to perturbations in vivo implies high noise and suggests rate coding in cortex. Nature, 466(7302):123- 127.
[50] Macke, J. H., Buesing, L.,, Cunningham, P., Yu, B. M., Shenoy, K. V., and Sahani, M. (2011). Empirical models of spiking in neural populations. Advances in neural information processing systems, 24.
[51] Mei, H. and Eisner, J. M. (2017). The neural hawkes process: A neurally self-modulating multivariate point process. Advances in neural information processing systems. 30.
[52] Mensi, S., Naud, R., and Gerstner, W. (2011). From stochastic nonlinear integrate-and-fire to generalized linear models. Advances in Neural Information Processing Systems, 24.
[53] Mizuseki, K., Sirota, A., Pastalkova, E., and Buzsáki, G. (2009). Theta oscillations provide temporal windows for local circuit computation in the entorhinal-hippocampp loop. Neuron, 64(2):267-280.
[54] Mizuseki, K., Sirota, A., Pastalkova, E., Diba, K., and Buzsáki, G. (2013). Multiple single unit recordings from different rat hippocampal and entorhinal regions while the animals were performing multiple behavioral tasks. CRCNS org.
[55] Møller, J.. Syversveen, A. R., and Waagepetersen, R. P. (1998). Log gaussian processes. Scandinavian journal of statistics, 25(3):451-482.
[56] Moreno-Bote, R., Beck, J., Kanitscheider, I., Pitkow, X., Latham, P., and Pouget, A. (2014). Information- limiting correlations. Nature neuroscience, 17(10):1410.
[57] Nagele, J., Herz, A. V., and Stemmler, M. B. (2020). Untethered firing fields and intermittent silences: Why grid-cell discharge is variable. Hippocampus.
[58] Nawrot, M. P., Boucsein, C., Molina, V. R., Riehle, A., Aertsen, A., and Rotter, S. (2008). Measurement of variability dynamics in cortical spike trains. Journal of neuroscience methods, 169(2):77--39.
[59] Omi, T., Ueda, N., and Aihara, K. (2019). Fully neural network based model for general temporal point processes. In Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché-Bue Buc. F., Fox, E., and Garnett, R., editors, Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc.
[60] Orbán, G., Berkes, P., Fiser, J., and Lengyel, M. (2016). Neural variability and sampling- probabilistic representations in the visual cortex. Neuron, 92(2):530-543
[61] Ostojic, S. (2011). Interspike interval distributions of spiking neurons driven by fluctuating inputs. Journal of neurophysiology 106(1):361-33)
(62] Paninski, I.. (2004). Maximum likelihood estimation of cascade point-process neural encoding models. Network: Computation in Neural Systems, 15(4):243-26..
13

---

[63] Peyrache, A. and Buzsáki, G. (2015). Extracellular recordings from multi-site silicon probes in the anterior thalamus and subicular formation of freely moving mice. CRCNS.
[64] Peyrache, A., Lacroix, M. M., Petersen, P. C., and Buzsáki, G. (2015). Internally organized mechanisms of the head direction sense. Nature neuroscience, 18(4):569-575.
[65] Pillow, J. W. (2009). Time-rescaling methods for the estimation and assessment of non-poisson neural encoding models. Advances in neural information processing systems, pages 1473-1481.
[66] Pillow, J. W., Shiens, J., Paninski, L.,, Sher, A., Litke, A. M., Chichilnisky, E., and Simoncelli, P. (2008). Spatio-temporal correlations and visual signalling in a complete neuronal population. Nature, 454(7207)
[67] Ponce-Alvarez. A., Thiele, A., Albright, T. D.. Stoner, G. R., and Deco, G. (2013). Stimulus- dependent variability and noise correlations in cortical mt neurons. Proceedings of the National Academy of Sciences, 110632):1312--13167
[68] Rad, K. R. and Paninski, L.. (2010). Efficient, adaptive estimation of two-dimensional firing rate surfaces via gaussian process methods. Network: Computation in Neural Systems, 21(3-4):14- 168.
[69] Ratnam, R. and Nelson, M. E. (2000). Nonrenewal statistics of electrosensory afferent spike trains: implications for the detection of weak sensory signals. Journal of Neuroscience, 20(17)6672-6683.
[70] Reich, D. S., Victor, D., Knight, B. W., Ozaki, T., and Kaplan, E. (1997), Response variability and timing precision of neuronal spike trains in vivo. Journal of neurophysiology. 77(5):2836 2841.
[71] Sahin, I. (1993). A generalization of renewal processes. Operations research letters, 13(4):259- 263.
[72] Shchur, O., Gao, N., Biloš, M., and Günnemann, S. (2020). Fast and flexible temporal point processes with triangular maps. Advances in Neural Information Processing Systems, 33:73-84.
[73] Shimazaki, H. and Shinomoto, S. (2007). A method for selecting the bin size of a time histogram. Neural computation, 19(6); 1503-1577
[74] Shimokawa, T., Koyama, S., and Shinomoto, S. (2010). A characterization of the time-rescaled gamma process as a model for spike trains. Journal of computational neuroscience, 29(1-2):88- 191.
[75] Shinomoto, S., Miura, K., and Koyama, S. (2005). A measure of local variation of inter-spike intervals. Biosystems, 79(1-3):67-72
[76] Skaggs, W. E., McNaughton, B. L., Wilson, M. A., and Barnes, C. A. (1996). Theta phase precession in hippocampal neuronal populations and the compression of temporal sequences. Hippocampus, 149-172.
[77] Song, S., Lee, J. A., Kiselev, I., Iyengar, V., Trapani, J. G., and Tania, N. (2018). Mathematical modeling and analyses of interspike-interaal of spontaneous activity in afferent neurons of the zebrafish lateral line. Scientific Reports, 8(1):1-14.
[78] Stein, R. B., Gossen, E. R., and Jones, K. E. (2005). Neuronal variability: noise or part of the signal? Nature Reviews Neuroscience, 6(5):389-397.
[79] Teh, Y. W. and Rao, V. (2011). Gaussian process modulated renewal processes. In Advances in Neural Information Processing Systems, pages 2474-2482
[80] Titsias, M. and Lawrence, N. D. (2010). Bayesian gaussian process latent variable model. In Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, pages 844-851.
14

---

[81] Truccolo, W. (2016). From point process observations to collective neural dynamics: Nonlinear hawkes process glms, low-dimensional dynamics and coarse graining. Journal of Physiology-Par Paris, 11 11044):33--37.
[82] Truccolo, W. and Donoghue, J. P. (2007). Nonparametric modeling of neural point processes via stochastic gradient boosting regression. Neural computation, 19(3):672-705.
[83] Truccolo, W., Eden, U. T., Fellows, M. R., Donoghue, J. P., and Brown, E. N. (2005). A point process framework for relating neural spiking activity to spiking history, neural ensemble, and extrinsic covariate effects. Journal of neurophysiology. 1874::074--189.
[84] Uifalussy, B. B. and Orbán, G. (2022). Sampling motion trajectories during hippocampal theta sequences. Elife. 11:77458..
[85] Weber, A. and Pillow, J. W. (2017). Capturing the dynamical repertoire of single neurons with generalized linear models. Neural computation, (9(1))3320-3328)
[86] Williams, A. H., Poole, B., Maheswaranathan. N., Dhawale, A. K., Fisher, T., Wilson, C. D., Brann, D. H., Trautmann, E. M., Ryu, S., Shusterman, R., et al. (2020). Discovering precise temporal patterns in large-scale neural recordings through robust and interpretable time warping. Neuron, 105(2):246-255))
Williams, C. K. and Rasmussen, C. E. (2006). Gaussian processes for machine learning, [87]
volume 2. MIT press Cambridge, MA.
[88] Wilson, J., Borovitskiy, V. Terenin, A., Mostowsky, P., and Deisenroth, M. (2020). Efficiently sampling functions from gaussian process posteriors. In International Conference on Machine Learning, pages 10292-10308 PMLR.
[89] Wilson, J. T., Borovitskiy, V., Terenin, A., Mostowsky, P., and Deisenroth, M. (2021). Pathwise conditioning of gaussian processes. The Journal of Machine Learning Research, 22(1):444-44787.
[90] Xiao, S., Yan, J., Yang, X., Zha, H., and Chu, S. (2017). Modeling the intensity function of point process via recurrent neural networks. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 31.
[91] Yang, R., Gupta, G., and Bogdan, (2019). Data-driven perception of neuron point process with unknown unknowns. In Proceedings of the 10th ACM/IEEE International Conference on Cyber-Physice Systems, pages 259-269.
[92] Yousefi, A., Amidi, Y., Nazari, B., and Eden, U. T. (2020). Assessing goodness-of-fit in marked point process models of neural population coding via time and rate rescaling. Neural Computation, 32(1 (0(1))2145-2186
[93] Yu, B. M., Cunningham, J. P., Santhanam, G., Ryu, S., Shenoy, K. V., and Sahani, M. (2008). coussian-process factor analysis for www-dimessional single-trial analysis of neural population activity. Advances in neural information processing systems, 21:1881-1888
15


Reviews:
[{"review_id": "QecPjTQB7b", "reviewer": "Reviewer_3xxY", "summary": "This paper proposes a Bayesian nonparametric approach using modulated renewal process to model neural spike train, and capable of modeling the covariability. The method includes a nonparametric priors on conditional interspike interval distribution and automatic relevance determination for lagging interspike interval based on renewal order. The method is evaluated on one synthetic data and two real datasets on animal navigation. It demonstrates better performance than the current SOTA baselines in its capability of capturing the interspike interval statistics.  ", "questions": "1. Give a detailed introduction about the parameters that would be optimized under eqn 18. and list other hyper-parameters for reproducibility. \n\n2. Evaluate the method on larger scale neural datasets. ", "limitations": "1. There is no potential negative societal impact of their work.\n2. The limitations about scalability of the proposed approach should be carefully addressed.  ", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. Motivation: the paper is well-motivated, modeling the spike train statistics is an important question, and the interspike interval statistics is an important property of neural dynamics, and potentially leads to identification of cell types, and functionality, etc. Designing a model that captures this property well is important.\n\n2. Method: it uses Bayesian nonparameteric approach that could fit complex data structures and patterns well, and it could infer the spike-history dependence using a data-driven approach.\n\n3. Results: the method demonstrate better accuracy than the current SOTA methods in multiple tasks and datasets.", "weaknesses": "1. Method: the model requires hyperparameter tuning of critical components includes $\\tau_w$  $K$, which might be hard to optimize, given its variability across neurons and datasets.\n\n2. Evaluation: the scalability of the method is a major concerns, as the datasets evaluated in this paper only has 9 neurons, ~30 units in each datasets. It's important to show how well the model performs in larger neural datasets.\n\n3. Complexity and computational cost: add evaluations based on the cost and speed of the proposed model and other baselines.\n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the Reviewer for their time and helpful feedback.\n\n**Q&A**\n\n*1. Give a detailed introduction about the parameters that would be optimized under eqn 18. and list other hyper-parameters for reproducibility.*\n\nThe NPNR model parameters consist of:\n- Gaussian process hyperparameters (kernel hyperparameters $\\theta_{GP}$, variational posterior mean $\\mu$ and covariance $\\Sigma$, inducing point locations $Z$, Gaussian process mean function parameters $\\tau_m$ and $b_m$). All these parameters are learned with gradient descent\n- Time warping timescale $\\tau_w$. This parameter is fixed in our experiments to the (empirical) mean ISI, but we also compared the cases where we jointly learn this parameter with the rest (see appendix figure 12)\n- Maximum lagging ISI order $K$ (integer). This hyperparameter is chosen in advance to a number we expect to be sufficiently large (in general, it is unlikely to have many neurons with significant history dependence beyond the past 2 ISIs based on studies of single neuron modeling). Indeed, figure 3D and 4D show this applies to most neurons in both datasets\n\nFor the camera-ready version, we will provide a short description at the end of the inference Section summarizing all the parameters and hyperparameters, stating which are learned and which are not.\n\n*2. Evaluate the method on larger scale neural datasets.*\n\nThe datasets in this work each contain around 30 selected neurons that are fitted simultaneously and have 1 to 2 million time steps in the training set (1 ms resolution), giving around 30-60 million data points in total (see appendix C.3). Most recent probabilistic models applied to neural data take in 100s of neurons but often with binned spike counts, reducing the temporal length to O(10^5) to O(10^6) in typical studies (e.g. “Scalable Bayesian GPFA with automatic relevance determination and discrete noise models”, Jensen et al. 2021; “Linear Time GPs for Inferring Latent Trajectories from Neural Spike Trains”, Dowling et al. 2023). Hence overall, in terms of number of data points this qualifies as a scalable method by recent standards. \n\nWe fit did not fully utilize the 11 GB GPU RAM that we ran the experiment on. With the same temporal batch size and other settings, we generally observe OOM issues around ~80 neurons, which can be increased even further by reducing the temporal batch size. We encourage the reviewer to try out the JAX implementation that is provided, which also provides fully reproducible runs of the experiments presented in this study.\n\nFor the camera-ready version, we will add these numbers (50-60 millions data points and total time points from appendix C.3) to Section 4.2 to explicitly demonstrate the scale of our datasets, which can be misleading when just looking at the number of neurons (~30 neurons).\n\n**Weaknesses and limitations**\n\n*1. Method: the model requires hyperparameter tuning of critical components includes tau_w, K, which might be hard to optimize, given its variability across neurons and datasets.*\n\nThe time warping process is chosen based on inductive biases from the neuroscience literature (see Section 3.1). We fix $\\tau_w$ instead of learning it as this leads to automatic relevance determination in the temporal kernel timescales (see figure 6A). In appendix figure 12, we show that the loss of performance due to fixing this parameter is minimal.\n\nChoosing $K$ is in a sense similar to choosing the number of layers in a deep neural network, or the dimensionality of latent spaces in latent variable models like VAEs. A rigorous approach would involve a grid search, but in the spirit of Bayesian models one can use a single high capacity model and perform automatic relevance determination (“Scalable Bayesian GPFA with automatic relevance determination and discrete noise models”, Jensen et al. 2021; “Bayesian Gaussian Process Latent Variable Model”, Titsias and Lawrence 2010) to get rid of redundant degrees of freedom as shown in figure 2B, 3D and 4D.\n\nIn the camera-ready paper, we will elaborate on the issue of hyperparameter selection and tuning as discussed above in Section 3.2, and include relevant references to similar procedures done in the literature. We will explicitly mention and discuss all our hyperparameter tuning results (appendix figure 12 and automatic relevance determination figure 2B, 3D and 4D) in the light of hyperparameter tuning.\n\n*2. Evaluation: the scalability of the method is a major concerns, as the datasets evaluated in this paper only has 9 neurons, ~30 units in each datasets. It's important to show how well the model performs in larger neural datasets.*\n\nSee response to question 2 in Q&A.\n\n*3. Complexity and computational cost: add evaluations based on the cost and speed of the proposed model and other baselines.*\n\nThe stochastic sparse variational Gaussian processes (Hensman et al. 2013) has computational complexity $O(N T M^2 + N M^3)$, where $M$ is the number of inducing points, $N$ the number of output dimensions (neurons) and $T$ is the temporal batch size. As we can shorten the batch size $T$, the predominant factor is $O(M^3)$ as we mention in Section 3.2 in the main text. Our method has the same computational scaling as stochastic sparse variational Gaussian processes (“Gaussian processes for big data”, Hensman et al. 2013), which can be applied to very large datasets due to mini-batching. Baseline models from the literature are adapted to use the same underlying SVGP machinery to model the input-output mapping or spike-history dependence, and hence inherits the similar complexity scaling up to multiplicative factors. A detailed discussion of model implementation and inference is presented in appendix B.2-5, where we present a large amount of details not suitable for the main text.\n\nFor the camera-ready paper, we will extend the discussion on computational complexity in Section 3.2 and appendix B.2-5, with the full $O(N T M^2 + N M^3)$ expression.\n\n"}, {"review_id": "aiqpA47M9L", "reviewer": "Reviewer_yMLW", "summary": "The variability of neural data is widely observed in many neuroscience experiments. Using statistical model to capture the variability structure plays an essential role in understanding neural computations. Generally, the variability of neural data is a result of non-stationary activities and dependencies on behavioral covariates. To tackle these challenges, the authors proposes a scalable Bayesian approach generalizing modulated renewal processes (NPNR) to analyze neural data variability. They develops a nonparametric generalization of modulated renewal processes beyond renewal order, which makes their method flexibly model irreducible “intrinsic” neural stochasticity and input-dependent variability. Furthermore, the authors apply stochastic variational inference and can fit the model to long time neural data given cubic time complexity in the number of inducing points. The performace of NPNR is evaluated on both synthetic and real neural datasets.", "questions": "* The proposed method captures ISI statistics using a spatio-temporal GP prior over CIF. Could the Neural Temporal Point Process (NTPP) perform similarly to your method in capturing ISI statistics? The CIF in NTPP is usually modeled by neural networks, and could this be more powerful to represent ISI distributions and capture ISI statistics?", "limitations": "N/A", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "* The proposed method can model two types of neural variability: (1) capturing spiking statistics from non-stationary data; (2) capturing modulation by behavioral covariates.\n* To achieve the desired non-stationarity, the proposed method uses time warping on $\\tau$, which avoids the use of non-stationary kernls and maintains the ability to draw samples by pathwise conditioning.\n* The proposed inference method provides an elegant approach to determine the spike-history dependence in ISI statistics.\n* The authors' exposition of their motivation, contribution, and conclusions from the experiments are comprehensive and clear.\n* The proposed method would provide an important set of contributions to the field of neural coding.", "weaknesses": "The proposed method is clearly written and well supported by experiments. I have nothing further to add here.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the Reviewer for their time and helpful feedback. \n\n*1. The proposed method captures ISI statistics using a spatio-temporal GP prior over CIF. Could the Neural Temporal Point Process (NTPP) perform similarly to your method in capturing ISI statistics? The CIF in NTPP is usually modeled by neural networks, and could this be more powerful to represent ISI distributions and capture ISI statistics?*\n\nIndeed, we are also familiar with the NTPP works and have considered them as candidates. However, the Gaussian process approach offers a few unique properties that are difficult to achieve with neural networks:\n- Automatic relevance determination of the lagged ISI dimensions by learning temporal kernel lengthscales (see figures 2B, 3D, 4D and 6A)\n- Principled regularization using the evidence lower bound and smoothness constraints on the CIF dependence on covariates $\\mathbf{x}_t$ and past ISIs $\\mathbf{\\Delta}_t$ (e.g. figure 12 comparing against different kernels that impose different levels of differentiability of the CIF (i.e. choosing different order Matern kernels)\n\nIn general, any method that involves flexible modeling of the CIF directly is closely related to the main idea/contribution of this work. Our work introduces nonparametric autoregressive dependence on past ISIs, which adds $K$ more input dimension compared to baselines which handle past dependencies using parametric spike history filters or renewal assumptions. For this reason, kernel methods such as Gaussian processes provide an elegant framework for handling this without overfitting as much as neural networks.\n\nA detailed comparison study of neural networks replacing the Gaussian process for modeling the CIF is out of scope for the current work and is more related to the general argument of Gaussian processes versus deep neural networks, which is an entire subfield on its own. Small examples of neural networks overfitting in cases when Gaussian processes do not are for instance given in Ref. 43 (“A universal probabilistic spike count model reveals ongoing modulation of neural variability”, Liu and Lengyel 2021). Additionally, there are many influential works investigating the relationship between these two classes (“Approximate Inference Turns Deep Networks into Gaussian Processes”, Khan et al. 2019 ; “Wide Feedforward or Recurrent Neural Networks of Any Architecture are Gaussian Processes”, Yang 2019).\n\nAn interesting case mentioned in our work is Ref.54 (“Fully Neural Network based Model for General Temporal Point Processes”, Omi et al. 2019), where the authors chose to model the cumulative hazard function (integral of the CIF) for temporal point processes. This requires a monotonic function approximator, which is most easily done with constrained neural networks and difficult to achieve with Gaussian process priors.\n\nIn the camera-ready paper, we will expand on the comparison to neural network approaches similar to the discussion above, and add relevant references on these topics.\n"}, {"review_id": "9ecU4LMgLG", "reviewer": "Reviewer_hah4", "summary": "The authors proposed a scalable Bayesian approach which generalizes modulated renewal processes using sparse variational Gaussian processes. They applied the proposed method to simulated and two real neural datasets and showed that the proposed method is effective on these datasets and outperforms other baseline methods.", "questions": "- Method section: all the formulas are described for 1 neuron. It might be clearer to write likelihood and loss function in terms of multiple neurons.\n\n- I wonder if the authors have done any comparisons to the parametric methods? e.g. Gao Y*, Archer E*, Paninski L, Cunningham JP (2016) Linear dynamical neural population models through nonlinear embeddings. NIPS 2016.\n\n- The proposed method doesn't incorporate any latent variables. It might be worth adding the latents to discover useful representations from the data and fit the data variability better. I wonder if the model would fit data worse if miss one covariate in the inputs? (e.g. only includes location and direction in covariate not theta phase in the hippocampus data.) I feel that adding latents might help with this as well.", "limitations": "The authors have discussed the limitations and future work of their proposed method.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The paper is well written. The authors have done extensive experiments to show the effectiveness of the proposed method. ", "weaknesses": "The proposed method doesn't incorporate any latent variables in the model. (see more in the questions section below.)", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the Reviewer for their time and helpful feedback. \n\n*1. Method Section: all the formulas are described for 1 neuron. It might be clearer to write likelihood and loss function in terms of multiple neurons.*\n\nThe current work inherently models each neuron separately while fitting the neural population simultaneously, hence the likelihood for all neurons is simply the sum of the log likelihoods in equation 18 \n\n$\\mathcal{L} = \\sum_n \\mathbb{E}\\_{ q( f_n \\| u_n) } \\left[ - y_{nt} f_{nt} + \\Delta t \\sum_{t=1}^T e^{f_{nt}} \\right] + D_{\\text{KL}}( q(\\mathbf{u}_n) |\\ p(\\mathbf{u}_n) ) )$\n\nwhere in $f_{nt}$ etc. we add the neuron index to the quantities. We will modify this in the camera-ready version to make this more explicit. As you suggest below, introducing latent variables would provide the ability to capture neural (noise) correlations. See our response to question 3 below for a brief discussion on this topic.\n\n*2. I wonder if the authors have done any comparisons to the parametric methods? e.g. Gao Y\\*, Archer E\\*, Paninski L, Cunningham JP (2016) Linear dynamical neural population models through nonlinear embeddings. NIPS 2016.*\n\n(Gao et al. 2016) is a model for spike counts with parametric noise models (Poisson and generalized count distributions). Our work is in the realm of spike trains where we have essentially binary time series after time discretization (see Section 3.2). We discussed related works in both realms, see the introductory paragraphs starting at line 40 to line 61, i.e. last paragraphs before the “contribution” paragraph. In this case of working with point processes/spike trains, spike count models are not appropriate as the only “spike count distribution” in a time bin is a Bernoulli distribution. Note that the main contribution of this paper is the flexible statistical modeling of the spiking process of a neural, or in machine learning terms flexible modeling of the noise distribution, per neuron. Works like (Gao et al. 2016) are latent variable models that model neural correlations despite having the same separable structure across neurons in the noise model.\n\nIn terms of parametric models, we conducted an extensive comparison of current parametric models as GLMs/SRMs under the name conditional processes (Ref 80, “Capturing the dynamical repertoire of single neurons with generalized linear models”, Weber and Pillow 2017; Ref. 17, “Non-parametric generalized linear model”, Dowling et al. 2020), renewal processes and other variants described in Section 2 (details in appendix B.2-3 and C). Results are shown for various models in figures 3 and 4, with the abbreviations described in the last sentence of the first paragraph in Section 4.2 (which include the simple GLMs as “RC cond. P” or radial cosine conditional Poisson processes).\n\nIn the camera-ready version, we will include (Gao et al. 2016) in the introductory paragraphs starting at line 40 to line 61 and more explicitly mention the lack of neural correlations in Section 3.1 (similar to the discussion in section C.1, “Vectorization over neurons”).\n\n*3. The proposed method doesn't incorporate any latent variables. It might be worth adding the latents to discover useful representations from the data and fit the data variability better. I wonder if the model would fit data worse if miss one covariate in the inputs? (e.g. only includes location and direction in covariate not theta phase in the hippocampus data.) I feel that adding latents might help with this as well.*\n\nIndeed, we agree that adding latent variables and thus making the model suitable for analyzing neural noise correlations is a major avenue for further research (see Section 5.1). Our current contribution is a model that is able to capture single neuron statistics with a sufficiently flexible noise model, which shows up as explaining all the data variance as measured by usual distribution comparison methods for point processes (see the KS plots figures 2C, 3B and 4B). However, such measures do not take into account neural correlations, and extensions to the KS framework have been proposed as in Ref. 28 (“Applying the multivariate time-rescaling theorem to neural population models”, Haslinger and Pipa 2011). Importantly, our current implementation is suitable for latent variable augmentations as we fit neurons in parallel like a coherent neural population (see appendix C.1 and the JAX implementation provided).\n\nAdding more covariates generally will improve the test ELL/predictive performance from the following perspective: our model can capture arbitrary spike train noise models (in theory), and the signal that is modeled is the modulation of the flexible noise distribution by those covariates. Adding more covariates allows the model to discover more fluctuating signals since we have more freedom if the covariates are uncorrelated. When the signal can account for more of the data variance, the noise model inferred can be sharper, i.e. more precise, and hence this will mathematically lead to higher log likelihoods. The extreme case here would be a noise model leading to conditional ISI densities that are delta distributions at the average ISI given some covariates (such models are approximated by Gamma renewal processes with shape parameters well above 1).\n\nIn the Reviewer’s example, including theta allows one to capture oscillations in the firing output of neurons with the theta rhythm (one can see the wiggled curve in figure 4C, note the contrast to smoother rates in figure 3C where no fast oscillating covariates as theta phase are included). If we fit the same model leaving out the theta phase covariate, we obtain similar goodness-of-fit measures (as the flexible noise model/ISI densities adapt to account for the reduction in signal) but lower test ELLs:\n\n| model | test ELL (nats / s) | \n|:--------:|:--------:|\n| $x$ position + head direction  | $13.92 \\pm 3.90$  | \n| $x$ position + head direction + theta phase | $14.39 \\pm 4.11$ |\n"}, {"review_id": "nydR6rYsti", "reviewer": "Reviewer_Mj1g", "summary": "This paper proposes the Bayesian nonparametric non-renewal process (NPNR) for inferencing both neural spiking intensity and variability. The tuning curve is based on a sparse variational Gaussian process (GP) prior, considering both spatial and temporal factors. They compare NPNR with other competitors on a synthetic dataset showing the capability of NPNR in inferencing the rate map and renewal density. On the two real-world neural datasets, they show that NPNR outperforms lots of competitors in terms of event prediction and interspike interval (ISI) recovery by different statistics combined with visualizations.", "questions": "* I think Eq. 15 should be $=$ rather than $\\propto$.\n* I'm wondering if this model can report a predictive log-likelihood using the mean as the estimation for the intensity in each time bin. In such a case, we can compare this model with other models (especially the simple GLM) to show that the proposed NPNR outperforms GLM? I'm expecting that this model will be slow but if the firing rates (tuning curve) recovery is better than GLM, the predictive log-likelihood (which is actually a golden criterion in neural latent variable models) should be better.\n* Can this model get information on the causal relationships between neurons? Is this model only dependent on time $t$ and external input $\\boldsymbol x$, but does not consider influences between neurons? From my understanding, this is not a latent variable model doing information extraction from coupled neurons (like dimensionality reduction), but getting the firing rate for each neuron, and the firing rate is mainly affected by the neuron itself and the external input $\\boldsymbol x$.", "limitations": "/", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 4, "strengths": "* Clear logical flow and presentation. Key maths are derived elegantly with lots of necessary details in the Appendix.\n* Both synthetic and real-world experiments are good and solid, and also supported by the code.\n* The literature review by the authors are exhaustive so that the comparison between different kind of models are clear and in detail.\n* The idea is new and intuitive, both preserve the interpretability and are not too simple.", "weaknesses": "* I'm very excited when seeing the model part. But when I get to Section 3.2, I feel a bit pity that we still need to do time discretization (convert the continuous timestamps TPP data to spike counts in time bins).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their time and helpful feedback on the work.\n\n**Q&A**\n\n*1. I think Eq. 15 should be = rather than ∝.*\n\nThis would indeed be true if the denominator (normalization constant) in equation 20 (appendix A) was 1. For a Poisson process with $t_i = 0$ (i.e. the most recent spike happened at time 0), this is true and the equality holds. In general, this does not hold and we need to compute the normalization constant which leads to valid ISI densities as plotted in the paper.\n\n*2. I'm wondering if this model can report a predictive log-likelihood using the mean as the estimation for the intensity in each time bin. In such a case, we can compare this model with other models (especially the simple GLM) to show that the proposed NPNR outperforms GLM? I'm expecting that this model will be slow but if the firing rates (tuning curve) recovery is better than GLM, the predictive log-likelihood (which is actually a golden criterion in neural latent variable models) should be better.*\n\nIn fact, the expected log likelihood (ELL) on the test sets is exactly the predictive log-likelihood you mention (see Section 4.2), defined as \n\n$ \\text{ELL} = \\sum_n \\mathbb{E}_{q(\\mathbf{f}_n)} [ \\log p(\\mathbf{y}_n | \\mathbf{f}_n) ] $ \n\nInstead of the predictive log-likelihood, we use the term test ELL as we are working in a variational framework rather than with a maximum likelihood approach. This metric also incorporates the variational posterior uncertainty in the intensity estimation in each bin, similar to metrics used in the Gaussian process literature (“Scalable Exact Inference in Multi-Output Gaussian Processes”, Bruinsma et al. 2020; “Dual Parameterization of Sparse Variational Gaussian Processes”, Adam et al. 2021).\n\nThe simple GLM (and more advanced variants of it) are actually already presented in figures 2, 3 and 4. These appear under the name “conditional Poisson” or “cond. P” for short (see Section 4.2 first paragraph last sentence). We apologize for the slightly unorthodox nomenclature, but this is necessary as we are comparing many variants and extensions of the GLM/SRM family (see Section 2.1.2). The Poisson GLM model has been extended with different kinds of spike-history filters, hence we also add the type of filter to the name. For example, “RC cond. P” is the radial cosine conditional Poisson in figures 3A and 4A, the classical GLM as referred to in the literature (Ref. 80, “Capturing the dynamical repertoire of single neurons with generalized linear models”, Weber and Pillow 2017).\n\nSince we run all baseline models using Gaussian process mappings in the same variational framework (see appendix B.2, B.3 and B.4), we also compute the test ELL in a similar way to our method (conventional GLMs use a maximum likelihood linear feature model for the input-output mapping). Since this metric is computed over the whole population, by adding latent input dimensions we would arrive at the predictive likelihood that you mentioned for latent variable models (although this is not done in this study due to the focus of this work).\n\n*3. Can this model get information on the causal relationships between neurons? Is this model only dependent on time $t$ and external input $x$, but does not consider influences between neurons? From my understanding, this is not a latent variable model doing information extraction from coupled neurons (like dimensionality reduction), but getting the firing rate for each neuron, and the firing rate is mainly affected by the neuron itself and the external input $x$.*\n\nGenerally, causal relationships are very hard or impossible to estimate from statistics alone (see Ref. 42, “Inferring structured connectivity from spike trains under negative-binomial generalized linear models”, Linderman et al. 2015). However, we could capture correlations between neurons through latent variables that would co-modulate the CIFs of individual neurons in the population. Alternatively, we could include the recent spiking history of other neurons in the input covariates in addition to the self-history of a neuron, similar to the GLM framework (“Spatio-temporal correlations and visual signaling in a complete neuronal population”, Pillow et al. 2008; “Non-parametric generalized linear model”, Dowling et al. 2020). In the NPNR framework the latter would naively lead to a significant increase in the input dimensions, and one interesting further work direction is to encode population spike-history in an efficient way.\n\nIndeed, the firing properties of each neuron in this study are only affected by its own history and the external covariates/inputs. For demonstrating our contribution, which is a flexible spike train noise model, this suffices and indeed improves over SOTA methods applied in the same manner.\n\n**Weaknesses and limitations**\n\n*I'm very excited when seeing the model part. But when I get to Section 3.2, I feel a bit pity that we still need to do time discretization (convert the continuous timestamps TPP data to spike counts in time bins).*\n\nWe agree with the reviewer and indeed drafts of this work were considering continuous time techniques that are faithful to the continuous-time formulation. However, exact methods for continuous time processes often rely on parametric assumptions that we try to alleviate (“Variational Inference for Gaussian Process Modulated Poisson Processes”, Lloyd et al. 2015; Ref. 18, “Temporal alignment and latent Gaussian process factor inference in population spike trains”, Duncker and Sahani 2018), or use MCMC, generalized thinning and other techniques (Ref. 74, “Gaussian process modulated renewal processes”, Teh and Rao 2011) that are not as straightforward to scale up to millions of time steps as in our work using stochastic variational inference.\n\nDo note that the discretization is done at 1 ms, which for neural data that has been sorted properly (no significant spike contamination) should not contain more than one spike (a binary time series).\n"}, {"review_id": "EbRdBVsK5x", "reviewer": "Reviewer_C2R4", "summary": "In this paper, the authors introduce the Bayesian nonparametric non-renewal (NPNR) process to model variability in neural spike trains with covariate dependence. The method generalizes modulated renewal processes using sparse variational Gaussian processes. Tested on synthetic data, as well as mouse head direction cell data and rat hippocampal place cell data, NPNR shows its superiority in terms of capturing interspike interval statistics and predictive power.", "questions": "1. The authors mention that the variational inference framework is scalable. In the synthetic and real-world experiments, the number of neurons/units are relatively small (less than 50). I wonder if there's any limitation of the method to scale up to larger ensemble systems with much more neurons?\n2. How does the proposed model perform in terms of predictive power when applied to synthetic data?\n3. The experiments on mouse head direction cell data and rat hippocampal place cell data shows the predictive performance of NPNR and baseline models by showing the expected log-likelihood. For the two datasets, the range of ELL for the models vary a lot and is hard to interpret. I wonder what value of ELL can be considered as a good prediction? For example, if NPNR is used for k-step ahead forecasts, how many steps can it achieve with at least 50% accuracy?\n4. The manuscript states that “Extending point process models with input-dependent variability has not been widely explored…” Multivariate auto-regressive frameworks and multiple covariates based models have been considered in \"A Granger causality measure for point process models of ensemble neural spiking activity.\" PLoS computational biology 7, no. 3 (2011): e1001110. \"Data-driven perception of neuron point process with unknown unknowns.\" In Proceedings of the 10th ACM/IEEE International Conference on Cyber-Physical Systems, pp. 259-269. 2019. \"Variance as a signature of neural computations during decision making.\" Neuron 69, no. 4 (2011): 818-831. In general the prior work needs to be more exhaustively checked and discussed, as of now it is biased and solely based on one group while there are similar and related works from other groups.\n5. Within the context of multiple neuronal recordings there is always the issue of interference and the problem that we cannot with certainty measure exactly N number of neurons. The activity of N neurons may be influenced by another P neurons so the question is how we can subtract the effect or perturbations in order to accurately model the N neurons and their covariates, etc. This again has been tackled in the neuroscience literature and the authors should check this related problem of understanding neural computations with unknown influences.\n6. In the experiments, 1-D, 2-D and 3-D x_t are considered for the NPNR modeling. I wonder if and how the inference can be affected by the dimension of x_t, number of samples and number of units/channels? I.e., how much data (in space and time) is needed for an efficient inference?\n", "limitations": "Not applicable in my opinion, this is a mathematical modeling paper with applications in neuroscience.", "rating": 5, "confidence": 5, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "+ Bayesian nonparametric non-renewal (NPNR) process to model variability in neural spike trains;\n+ Validation of the approach on synthetic data and mouse head direction cell data and rat hippocampal place cell data;\n", "weaknesses": "- There are some aspects that while they may be hard to analyze theoretically at least in the experimental part can be investigated through simulations. For example, how is the inference affected by the dimension of x_t, number of unmeasured or unknown neurons  that can act as perturbations, the number of samples and number of units/channels? I.e., how much data (in space and time) is needed for an efficient inference?\n- Within the experimental investigation, it is unclear what ELL values can be considered as good predictions. For example, if NPNR is used for k-step ahead forecasts, how many steps can it achieve with at least 50% accuracy?\n- Although this is a minor issue, the reader has some difficulties at seeing all the plots well from Figure 2 in particular the ISI probability plots. I assume they fit well gamma or exponential distribution.\n- A major issue is to check the literature and in unbiased way provide an accurate description even if the problems considered by prior work have a different or more advanced setup. \n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the Reviewer for their time and helpful feedback. We provide brief comments on the mentioned weaknesses in the relevant question.\n\n*1. The authors mention that the variational inference framework is scalable...*\n\nThe datasets in this work each contain around 30 selected neurons that are fitted simultaneously and have 1 to 2 million time steps in the training set, giving around 30-60 million data points in total (see appendix C.3). Most recent probabilistic models applied to neural data take in 100s of neurons but with binned spike counts, reducing the temporal length to $O(10^5)$ typically (“Scalable Bayesian GPFA with automatic relevance determination and discrete noise models”, Jensen et al. 2021; “Linear Time GPs for Inferring Latent Trajectories from Neural Spike Trains”, Dowling et al. 2023). Hence overall, in terms of number of data points this qualifies as a scalable method by recent standards.\n\nFor the camera-ready version, we will add these numbers to Section 4.2 to explicitly demonstrate the scale of our datasets.\n\n*2. How does the proposed model perform in terms of predictive power when applied to synthetic data?*\n\nWhen extending the synthetic experiment on 5 different held-out synthetic datasets:\n\n| model | test ELL (nats / s) | \n|:--------:|:--------:|\n|  Poisson   | $40.55 \\pm 2.09$  | \n|  Gamma renewal   | $47.07 \\pm 2.10$ | \n|  conditional Poisson   | $40.82 \\pm 2.22$ | \n|  NPNR (ours)    | $44.93 \\pm 2.30$ | \n\nNote the Gamma renewal process is within model class for 3 of the synthetic neurons. We will include these numerical results in the camera-ready version in the appendix.\n\n*3. The experiments on mouse head direction cell data and rat hippocampal place cell data...*\n\nThe two datasets are very different in terms of the spiking statistics when comparing the figures 3 and 4, hence the differences in the ELL values across the two datasets is consistent. Such numbers nevertheless can be compared in a relative sense within datasets. The test ELLs are also useful for fair comparison of predictive powers for qualitatively different baseline models (Ref. 2, “Construction and analysis of non-Poisson stimulus-response models of neural spiking activity”, Barbieri et al. 2001; Ref. 60, “Time-rescaling methods for the estimation and assessment of non-Poisson neural encoding models”, Pillow 2009; Ref. 17, “Non-parametric generalized linear model”, Dowling et al. 2020).\n\nOur model requires autoregressive sampling of posterior spike trains as seen in figures 3C and 4C, with each draw having a different CIF. The k-step ahead forecast you proposed is therefore not clear to us; prediction accuracy of the held-out spikes after k steps will vary across different posterior spike train samples.\n\nWe will point out the difference in ELL values between the two datasets more explicitly in the camera-ready paper with a brief discussion similar to the above.\n\n*4. The manuscript states that...*\n\nWe would like to thank the Reviewer for pointing out these works, some of which were not known to us.\n\n- **Kim et al. 2011**: a parametric CIF with a fixed interval history dependence, with a focus on Granger causality structure inference\n\n- **Yang et al. 2019**: a parametric GLM framework with unobserved or latent inputs, while considering the dependence of the spike-history on the external/unobserved covariates\n\n- **Churchland et al. 2011**: great example of the relevance of studying neural variability for neural coding\n\nFor the camera-ready paper, we will add these papers to the literature review in the relevant introduction paragraphs with appropriate discussion. Our contribution remains unique in the scalable Bayesian nonparametric modeling of spike-history dependence with emphasis on covariate-dependent modulation of instantaneous spiking variability, the autoregressive dependence in terms of interspike intervals rather than spike-train windows, and finally our perspective on the connection to conditional ISI distributions when conceptualizing spiking variability extracted with our model. We will qualify the incriminated sentence in our Introduction accordingly.\n\n*5. Within the context of multiple neuronal recordings...*\n\nWe agree strongly with the Reviewer on the importance of investigating these issues.\n\nThe effect of ignoring subsets of neurons generally leads to apparent noise correlations in the neurons we record, which we can capture by augmenting our model with latent variables shared across the population. Overall, capturing neural correlations is orthogonal to our main contribution of modeling the single neuron responses in a more flexible manner, and hence beyond the scope of this paper. Another related tangential topic is the downstream effect of spike sorting contamination (e.g. Ref. 86, “Assessing goodness-of-fit in marked point process models of neural population coding via time and rate rescaling”, Yousefi, A. et al. 2020).\n\n*6. In the experiments, 1-D, 2-D and 3-D x_t are considered for the NPNR modeling...*\n\nIn this study, our synthetic data has 2D inputs $x_t$ and has a million time points, and we were able to accurately recover the ground truth (figure 2). Real data had either 1D or 3D $x_t$ and ~2 million time points. Based on the synthetic experiment, we can be quite confident that this was in a sufficient data regime. Our method builds on sparse variational Gaussian processes, and hence inherits the convergence properties and data scaling from such models (Ref. 63, “Efficient, adaptive estimation of two-dimensional firing rate surfaces via Gaussian process methods”, Rad and Paninski 2010; “Rates of Convergence for Sparse Variational Gaussian Process Regression”, Burt et al. 2019; and the references in appendix B.1 on technical papers introducing sparse variational GPs).\n\nWe will add and expand the relevant references in appendix B.1 to address these questions for the camera-ready paper, as well as add this discussion topic to Section 3.2 with the relevant references."}]

Author rebuttal:

Meta-review: This is a solid piece of work using time-warped Gaussian processes to model the conditional intensity function / inter-spike interval distribution. While the approach draws on standard machinery, it shows promising performance on a variety of multi-neuronal spike train recordings. The experimental validation is thorough and the presentation is clear. This paper is a valuable contribution to the statistical neuroscience literature.
Decision: Accept (poster)