A universal probabilistic spike count model reveals ongoing modulation of neural variability
Authors: David Liu
Máté Lengyel
Venue: NeurIPS 2021 Poster (2021)
Keywords: neural coding
neural variability
neural correlations
Gaussian processes
discrete distributions
latent variable modelling

Abstract: Neural responses are variable: even under identical experimental conditions, single neuron and population responses typically differ from trial to trial and across time. Recent work has demonstrated that this variability has predictable structure, can be modulated by sensory input and behaviour, and bears critical signatures of the underlying network dynamics and computations. However, current methods for characterising neural variability are primarily geared towards sensory coding in the laboratory: they require trials with repeatable experimental stimuli and behavioural covariates. In addition, they make strong assumptions about the parametric form of variability, rely on assumption-free but data-inefficient histogram-based approaches, or are altogether ill-suited for capturing variability modulation by covariates. Here we present a universal probabilistic spike count model that eliminates these shortcomings. Our method builds on sparse Gaussian processes and can model arbitrary spike count distributions (SCDs) with flexible dependence on observed as well as latent covariates, using scalable variational inference to jointly infer the covariate-to-SCD mappings and latent trajectories in a data efficient way. Without requiring repeatable trials, it can flexibly capture covariate-dependent joint SCDs, and provide interpretable latent causes underlying the statistical dependencies between neurons. We apply the model to recordings from a canonical non-sensory neural population: head direction cells in the mouse. We find that variability in these cells defies a simple parametric relationship with mean spike count as assumed in standard models, its modulation by external covariates can be comparably strong to that of the mean firing rate, and slow low-dimensional latent factors explain away neural correlations. Our approach paves the way to understanding the mechanisms and computations underlying neural variability under naturalistic conditions, beyond the realm of sensory coding with repeatable stimuli.

A universal probabilistic spike count model reveals ongoing modulation of neural 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 uc
Department of Cognitive Science Central European University
Abstract
Neural responses are variable: even under identical experimental conditions, single neuron and population responses typically differ from trial trial and across time. Recent work has demonstrated that this variability has predictable structure, can be modulated by sensory input and behaviour, and bears critical signatures of the underlying network dynamics and computations. However, current methods for characterising neural variability are primarily geared towards sensory coding in the laboratory: they require trials with repeatable experimental stimuli and behavioural covariates. In addition, they make strong assumptions about the parametric form of variability, rely on assumption- or are altogether ill-suited for capturing variability modulation by covariates. Here we present universal probabilistic spike count model that eliminates these short- comings. Our method builds on sparse Gaussian processes and can model arbitrary spike count distributions (SCDs) with flexible dependence observed well latent covariates, using scalable variational inference to jointly infer the covariate- to-SCD mappings and latent trajectories in data efficient way. Without requiring repeatable trials, it can flexibly capture -epeedeed joint SCDs, and pro- vide interpretable latent causes underlying the statistical dependencies between neurons. We apply the model to recordings from a canonical non-sensory neural population: head direction cells in the mouse. We find that variability in these cells defies a simple parametric relationship with mean spike count as assumed standard models, its modulation by external covariates to that of the mean firing rate, and slow low-dimensional latent factors explain away neural correlations. Our approach paves the way understanding the mecha- nisms and computations underlying neural variability under naturalistic conditions, beyond the realm of sensory coding with repeatable stimuli.
but data-inefficient histogram-based approaches,
be comparably strong
Introduction
Classical analyses of neural coding are based mean spike counts neural firing rates. Indeed, some of the most paradigmatic examples of the neural code were discovered by regressing neural to characterize their tuning properties. However, neural spiking is generally not regular. Recordings from many cortical areas show significantly different activity patterns within and across identical trials 7 despite fixing experimentally controlled variables. This irregularity is also seen in continual neural recordings without trial structure [8]. The resulting variability has classically been characterised as 'Poisson' with Fano factor (variance mean ratio) of one [9] but experimental data also often exhibits significantly more [10 S [2] and sometimes less [13 [4] variability, respectively referred to as over- underdispersion. Moreover, experimental studies have revealed that neural variability generally depends on stimulus input and behaviour [5]
firing rates to particular sensory stimuli 21 or behavioural covariates 3 5
[8], and exhibits structured shared
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

variability ('noise correlations') across neurons even after conditioning on such covariates. Such cor- relations can have important consequences for decoding information from neural population activity
119 20 21 and reveal key properties of the underlying circuit dynamics [22]. Moreover, theories of neural representations of uncertainty have assigned computational significance to variability as
signature of Bayesian inference [23 24 25 26]. Thus, just as classical tuning curves for firing
rates have been crucial for understanding some of the fundamental properties of the neural code, principled statistical characterisation of neural variability, and its dependence on stimulus and behavioral covariates, is a key step towards understanding the dynamics of neural circuits and the computations they subserve.
The traditional approach to characterising neural variability has been pioneered in sensory areas, and relies on repeatable trial structure with sufficiently large number of trials using identical Variability in this case can be quantified by simple summary statistics of spike counts across trials of the same condition. However, this approach does not readily generalise to more naturalistic conditions where covariates cannot be precisely controlled and repeated in experiment. This more general setting requires statistical methods that take into account temporal variation of covariates for predicting neural count activity. Generalised Linear Models are popular choice [29], but they only model the dependence of firing rates on covariates with changes in variability directly coupled to changes in the rate inherent to Poisson spiking. More complex methods for inferring neural tuning 33 B1 and latent structure [32 33 34 35
stimulus and behavioral correlates 27 28
similarly use restrictive parametric families for spike count distributions, and thus also cannot model changes in variability that are not 'just' consequence of changes in mean counts firing rates. Conversely, statistical models capable of capturing arbitrary single neuron count statisties, such as histogram-based approaches or copulas 36), do not incorporate dependencies on covariates.
Here we unify these separate approaches. resulting in a single framework for jointly inferring neural tuning, single neuron count statistics, neural correlations, and latent structure. Our semi-parametric approach leads to the universal count model (UCM) for counts ranging from 0 to K, in the sense that we can model arbitrary distributions over the joint count space of size (K + 1) of N neurons. The trade-off between computational overhead and model expressivity is controlled by hyperparameters, with expressivity upper bounded by the true universal model. Our approach extends the idea of a universal binary count model BV to a finite range of integer counts, while allowing flexible dependence observed and latent covariates to model non-stationary neural activity and correlations.
The flexibility reduces biases from restrictive assumptions in any of the model components. Scalability is maintained by leveraging sparse Gaussian processes (38 with mini-batching (39 40 to handle the
size of modern neural recordings.
We first define the UCM, and then describe how to interpret as well as evaluate model fits. As our model is able to capture arbitrary single neuron statisties, we build on the Kolmogorov-Sming test to construct more absolute goodness-of-fi measures. After validating our method on synthetic data that cannot be captured by currently used methods, we apply the model to electrophyyiological recordings from two distinct brain regions in mice that show significant tuning the head direction of the animal [11 [22. We find that (1) neural activity tends to be less dispersed than common Poisson-like models at higher firing rates, and more dispersed low rates; (2) mean and variance of counts defy a simple parametric relationship imposed by parametric count distribution families; (3) variability modulation by behaviour can be comparable or even exceed that of the mean count or firing rate; (4) a two-dimensional latent trajectory varying on timescales of ~ 1 is sufficient to explain away neural correlations but not the non-Poisson nature of single neuron variability. Finally, we discuss related work, limitations and proposed extensions of our model.
2 Universal count model
Spike count activity of N neurons recorded into T' time bins is formally represented as an N-dimensional time series of non-negative integers. Due to biological constraints, the possible spike counts have some finite upper bound K, taken as the highest observed count. We denote probabilities of spike count distribution (SCD) by a vector of length K + 1, and use II to denote the collection of vectors The for neurons n and time steps l. Additionally, we denote the count activity by a matrix Y E [0, K] NxT XI with elements Ynt. Input covariates are observed X E Tx (e.g. animal speed)
Notation
or latent Z E xxD: (to capture e.g. attention), with range depending on topology [3]. We denote their clements ztd and 2iq with observed and latent dimension d and q. respectively.

---

X
8
time
firing
Fant
2n
TTI (rate)
5.0 Hz
Tl (rate)
Tit
W
b
factor
prob.
ITL (FF)
max
T
spikes K
min
X
Figure 1: Schematic of the UCM and the workflow. Left: graphical model corresponding to Equation with shaded circles observed, open circles latent, and squares deterministic variables. Filled dots represent fixed quantities. Middle: example inference of model posterior Equation with inferred latent trajectories (green, top) and covariate-dependent SCDs (blue, bottom) that depend on both observed and latent covariates. Note we only show the posterior over single SCD evaluated on (x,2) grid, whereas the full posterior defines SCDs over all neurons. Right: obtaining interpretable spike count statistics from the SCDs (see (sessectcn 2.3). Examples show firing rate and Fano factor tuning curves over observed x and latent 2 covariates, either jointly (heatmaps) or marginalized (grey curves). The depth of modulation in marginalized tuning curves is used extract a tuning index (TT) for the chosen subsets of covariates, see Equation
2.1 Generative model
The big picture is to model counts Y with dependence on X. For each neuron, our model consists of C' Gaussian process (GP) priors, basis expansion RRC C and a linear-softmax mapping
Ztq p(Z; 0Pt)
hen()) GP((, kien GPP)
fent = hen (Iz, 2,)
(1)
Wnt softmax (,, (ffrt) bn))
Discrete (T..)
where W (K++))x b IR K+1 and the GP covariance functions ken have hyperparameters 6GP The use of non-parametric GP mappings with point estimates for W and leads to semi remi-parametric model with parameters 0, see details in Appenddi E The overall generative model (Y is depicted schematically in Figure conceptually similar to Dirichlet priors (37 but allowing non non-paraeetric dependence
Note the model specifies prior p(TT[X) over joint SCDs,
X. With
latent input Z, our model can flexibly describe multivariate dependencies in joint SCDs as conditional addition, p(Z)
independence across neurons longer holds when marginalizing over 44
models temporal correlations in the latent states. We use Markovian priors (details in subsecctiion E.2
T
pe (Z) po(2))
II po p=(2)|z111 2=2
with absorbed into model parameters 9 for compactness. This allows the model to flexibly capture both neural and temporal correlations in Y. To attain scalability,
sparse GPs (38).
Depending on C and basis functions Φ(-), we obtain an approximation the universal prior on joint SCDs, with 'universal' referring to the ability to capture any joint SCD over all neurons. Arbitrary single neuron statistics can be captured when
with Φ(f) f. but fitting is computationally expensive when N xxG. > 1. For capturing all correlations, the model also requires a sufficiently large latent space. One controls the trade-off between model expressiveness and computational overhead through C' and Φb Larger expansions & allow one to model count distributions more expressively with small C, e.g. the element-win Emeaa-exponeenial Φ(f) = (f1,eh ef,..... Sc. ele) covers a range of distributions including the truncated Poisson with only C = 1 (see subsection A.3

---

2.2 Stochastic variational inference and learning
For the joint model distribution Po(Y, II, Z|X) P(Y IT) po (I) 2) Po(Z), with count distribu-
tions PII we approximate the posterior by x0.x.p (II, x that factorizes in the form
.N
(II ([(
40.x (T|x, (z)) = II
90.x (II, X, 2)
(3)
with tp and the variational parameters for latent states and the sparse Gaussian process posterior (Appendix E), respectively. Note that we use factorized normal q(Z) for Euclidean Z, and wrapped normal for circular Z based on the framework of reparameterized Lie groups 45 43 The posterior over count probabilities _0.x (II] Z) is defined as mapping the sparse Gaussian process posterior 20.x( (F) X, ?) through IT(F) (Equation 1). deterministic mapping. This is analytically
intractable, so in practice it is represented by Monte Carlo samples. Using stochastic variational inference [46], we minimize an upper bound on the negative log marginal likelihood
P[I][I pa (II Z) pa (Z)
Fo Fa (Y -Ez-g (2) Bi (
log
ga.x (1Ix Z) gz(()
known as the variational free energy. This objective leads to tractable terms subsection E.I). allowing
us to infer the approximate posterior well lower bound of the log marginal likelihood 47 40
We use Adam [88 for optimization, see details of implementation and model fitting in Appendix E
2.3 Obtaining interpretable spike count statisties from the model
Characterizing spike count distributions From the posterior q(II D we can compute samples of the posterior of any statistic of spike counts function of covariates. Single neuron statistics in particular can be characterized by tuning curves for both mean firing rates and Fano factors (FF)
1 A EPIrI [Y]
Var [Y
p(X)
FF(X)
(5)
EP(rmn [Y]
with time bin length A. The model also quantifies private neuron variability that cannot be explained away by regressing to shared input (both observed and latent) through (y (.)).
To quantify the sensitivity of a some aspect of neuron activity to a set of covariates I., we define a tuning index (TT) with respect to count statistic T2 (x.)
TI maxæ, Ta(x:)) ming. (x.))
(6)
maxe. (x,) + ming. Ty (x.)
with (x.) evaluated under the mean posterior SCD marginalized over all other covariate dimensions complementary to I.. These marginalized distributions are computed using observed input x (subsection D.5). Resulting marginalized tuning curves for TIs are depicted conceptually in (igure
Generalized Z-scores and noise correlations The deviation of activity from the predicted statis- ties is commonly quantified through Z-scores X 49
which are computed as (y (v))/ / (v) with (y) being the mean count in some time bin. If neural activity follows a Poisson distribution, the distribution of Z-scores asymptotically tends unit normal when average counts (3) > 1 (Appendix Bp. To generalize the normality of the 7-score for arbitrary counts and SCDs, we use
pyt
8-1
< = <= with u(y))
p(7) dū
Σ P(k) +eP(v), 100
(7)
k=0
with (-) the unit Gaussian edf., and dequantization noise to get continuous quantiles u U(0,1) and generalized Z-scores & ~ NO0 1) from the probability integral transform.
With ç. one can completely describe single neuron statistics with respect to the model. Correlations in the neural activity however will cause & to be correlated. We define generalized lagged correlations riy (A) E [-1, 1] and Fisher Z E R that is more convenient for statistical testing
1 +
ri (A) (E:(1) + A)), Zrisher log
r
which describes spatio-temporal correlations not captured by the model. Noise correlations [50] refer the case of A when rij becomes symmetric.
For notational convenience, X denotes both observed and latent covariates here.

---

2.4 Assessing model fit
Our model depends on a hyperparameter C Σ K that trades off flexibility with computa- tional burden. In practice one likely captures the neural activity accurately with C' well be- low K and simple basis expansion as the iinearexponential above or quadratic-exponentia ((f) = JR, e! ..... Ji ....... This can be quantified by the statistical measures provided below and allows us to select appropriate hyperparameters to capture the data sufficiently well.
To assess the model fit to neural spike count data, a conventional machine learning approach is to evaluate the expected lig-liielihood of the posterior predictive distribution on held-out data Y. leading to the cross-validated log-likelihood
evLL =
Dlog P(Y[T))
where cross-validate over the neuron dimension by using the majority of neurons to infer the latent states (Z) in the held-out segment of the data, and then evaluate Equation over the remaining neurons. Without latent variables, we simply take the expectation with respect x)) However, the cvLL does not reveal how well the data is described by the model in an absolute sense. Likelihood bootstrap methods are possible (28), but become cumbersome for large datasets. To assess whether the neural data is statistically distinguishable from the single neuron statistics predicted by the model, we use the Kolmogorov-Smin framework
(5l) with from Equation across time steps
Tks 1 (ur) Uz
(10)
with empirical distribution function (u) (for details see Appendix B. This scalar number is positive and does not indicate whether the data is under- or overdispersed relative to the model. For this, useful measure of dispersion is the logarithm of the variance of & with correction
( 372
1
Tos log (5t):
(11)
which provides a real number with positive and negative sign indicating over- and underdispersion, respectively. Its sampling distribution under N 1) asymptotically normal, centered around 0 (due to the additive parenthetical term) with a variance depending on the number of timesteps T (subsection B.3) This extends the notion of over- and underdispersion beyond the usual definition of dispersion relative to Poisson models [52]. To quantify whether the model has captured noise correlations in the data, we compute & with respect to the mean posterior predictive distribution
N (x = / II II
Eg( (2) [P(yn TTe)]
(x(2)) dz,
(12)
which whitens E. hence reducing noise correlations in Equation 8 if correlations are explained away by co-modulation of neurons due to shared low-dimensional factors [44). These are captured through latent states Z using the posterior q,p (Z(Z). inferred from the same data used to compute ç.
Overall, this can be interpreted as treating Z as if it was part of the observed input to the model.
3 Results
In the following, we use C 3 3 with an wlement-wise linear-exponential basis expansion (subsec [ion 2.. This empirically provided sufficient model flexibility as seen in goodness-of-fit metrics. We use an RBF kernel with cosine distances in case of angular input dimensions (subsection E.3).
3.1 Synthetic data
Animals maintain an internal estimate of their head direction statistical models of head direction cell populations [54] for validating the ability of the UCM to capture complex count statistics, as well as neural correlations through latent structure. The task is to jointly recover the ground truth count likelihoods, their tuning to covariates, and latent trajectories relevant from activity generated using two synthetic populations. The first population was generated with a parametric heteroscedastic Conway-Maxwell-Possse mean and variance modulation as well as simultaneously over- and underdispersed activity (Fano
42 53
Here, we extend simple
(CMP) model (55), which has decoupled
5

---

A
Universal (GP)
posterior q₂(2)
neuron neuron 17
2m
truth
truth
pood
fit
time
U (GP)
(ANN)
Tus ADvLL (10%)
count
5
10
ADYLL (10")
2.5
mean
mean
mouron 7
0.05
hNB @a
RMSE
11.1
2
t
1.0
F
0.00
0.5
Poleson
U (GP) (ANN))
2rT
Polsson
hNB
(GP)
(
2nz
head direction
head direction (Iruth)
Poisson
Universal (X)
Universal (X,Z)
DovLL (10º)
Poisson (X)
(X.))
Yo
N
N
Universal (X,Z)
0.1
inferred
0.0
2m
-0.1
truth
time (50
Fiehar
neuron
0.1
mean
0-
0
-5
neuron
Figure 2: Model validation with two synthetic head direction cell populations. (A) Applying Poisson, heteroscedastic negative binomial (hNB) and universal (U, with either GP or ANN mappings) models to synthetic data from the heteroscedastic CMP population. Error bars indicate s.e.m. over cross-validation runs. Shaded regions for tuning curves and Tks indicate 95% CI. AcvLL the difference w.r.t. Poisson baseline. Left: regression models, visualizing SCDs (top right), tuning curves (bottom right), and fitting scores (bottom left) for two representative cells. Right: models with angular latent variable, visualizing inferred latent states (top right), comparing the GP to ANN model (bottom right), and plotting fitting scores (bottom left). between the inferred latent and ground truth uses the geodesic distance the ring (subsection D.9). (B) Applying regression (Poisson, Universal (X)) and joint latent-observed (Universal (X,Z)) models the modulated Poisson population data. Left: three columns showing progressively how single neurons variability, with & and Tps (middle), and noise correlations, with rij and Fisher Z (bottom), captured (see subsection 2.4). Right: AcvLL for all models (top) and visualization of the joint latent-observed model input spaces (middle) and tuning curves (bottom).
-ee---mann-uuueee errors (RMSE)
factors above and below 1). The second population consists of Poisson neurons tuned to head direction and additional hidden signal, which gives rise apparent overdispersion (28) well as noise correlations when only regressing to observed covariates. For mathematical details the count distributions and synthetic populations, see Appendix and Appendix D respectively.
compare our UCM to the Poisson GP model 33 and the heteroscedastic negative binomial GP model, pon-maraettiic extension of [55]. To show data-efficiency and regularization benefits of

---

GPs, we also compare to UCM with an artificial neural network (ANN) mapping replacing the GP mapping. For details of the baseline models, see Appenndid D. For cross-validation we split the data into 10 roughly equal non-overlapping segments, and validated on 3 chosen segments that were evenly spread out across the data. When a latent space was present, we used 90% of the neurons to infer the latent signal while validating on the remaining neurons, and repeated this for non-overlapping subsets. We rescale the log-likelshoos by the ratio of total neurons to neurons in subset and then take the average over all subsets to obtain comparable cross-validation runs to regression.
Figure 2/A shows that only the UCM successfully captures the heteroscedastic CMP data, indicated by Tks. Baseline models cannot capture frequent cases where the Fano factor drops below 1. In addition, we observe that using a Bayesian GP over an ANN mapping in the model leads to a reduction in overfitting, especially in the latent setting (Figure 2A)). Therefore, all other analyses with the UCM reported here used the GP mapping. Figure 2B shows that the modulated Poisson population activity is seen by a Poisson regression model as overdispersed, indicated by TDs- Our model flexibly captures the overdispersed single neuron statistics, independent from noise correlations raj that are captured when we introduce Euclidean latent dimension. As expected, the scatter plots shows whitening under the posterior predictive distribution when the correlations are captured.
3.2 Mouse head direction cells
We apply the UCM to a recording of 33 head direction cells in the anterodorsal thalamic nucleus (ADn) and the postsubiculum (PoS) of freely moving mice [41 see see|subsection D.2 Neural data was binned into 40 ms intervals, giving K- 11. Note that observed count statistics differ with bin size, see subsection C.2 for a discussion. Regression was performed against head direction (HD), angular head velocity (AHV), animal speed, two-dimensional position in arena, and absolute time, which collectively form X in this model. We used 64 inducing points for regression, and added 8 for every latent dimension added (Appendix E}. valis-valiittion was performed as in synthetic experiments, but with 6 validation segments and subsets of 85% of neurons to infer the latent signal.
Figure 3/4 shows that for regression hNB overfits and performs worst, despite containing Poisson as a special case. However, this limit is not reached in practice due to the numerical implementation, sce|Appendix A] Only the UCM captures the training data satisfactorily with respect to confidence bounds for Iss and Lis, although the data remained slightly underdispersed to the model with Tbs values slightly skewed to negative. Compared to the Poisson model, the evLL only slightly higher for the UCM as the data deviates from Poisson statistics in subtle ways. We see both FF above and below 1 (over- and underdispersed) across the neural firing range in FFigure 3B, with quite some neurons crossing 1. Correspondingly. FF-mean correlations coefficients are often negative. Their spread away from 111 indicates firing rate and FF do not generally satisfy a simple relationship, especially for examples such as cell 27. Furthermore, ADn neurons seem to deviate less from Poisson statistics. From Figure 3C, we note in particular that FFs tend to decrease at the preferred head direction, but rise transiently that tuning to speed and time primarily modulates variability rather than firing rates. All of this is impossible to pick up with baseline models, which constrain FF > 1 as well as FF increasing with firing rate AAppendix A). Finally, we see more tuning of the firing rate to position in PoS cells.
the head direction approaches the preferred value. We also see
When adding latent dimensions, Figure 3D shows peak in the cvLL at two dimensions. where correspondingly the Fisher Z distribution starts to match the unit normal well. Kernel length scales however did not indicate redundant latent subspaces for higher dimensions as expected for automatic relevance determination, likely due to mixing of latent dimensions. Notice the noise correlation patterns in FFiuree 3 tend to show positive correlations for similarly tuned neurons roughly around the diagonal of blocks, as expected from ring attractor models (22). Intrinsic neuron variability, roughly quantified by the average FF, further decreased and thus become even more underdispersed when considering additional tuning to latents, in particular for ADn. In addition, latent signals primarily modulate firing rate as seen from TIs in Figure 3F. When looking at time scales of covariates in Figure 3G (computed as the decay time constant of the autocorrelogram (Appendix D), the latent processes seem to vary on time scales right in the gap of behavioural time scales.
Tangential to our main contribution of characterising the detailed structure of neural variability, our results have another novel element that does not specifically rely on the UCM. Using GP-based non-parametric methods, we successfully estimated the tuning of cells to as many as 8 different covariates (6 observed + 2 latent, see Figure 3G) in a statistically sound fashion, while previous
7

---

A
B
Tis (10 (10-))
cell
AoUL (103)
Fano factor
0.0
ADn
Po5
Poisson hN Universal
-0.1
hNB U
firing rate (Hz)
50
average
1.0
FF-mean
C
cell 12
cell 27
rate (HZ)
20
position
25.5
94.0
25
100
25
100
25
100
0
0
0-
heed direction
x position
AHV (radis)
speed (om/s)
time (min)
0.6
0.2
0.2
T(
0.0
0.0
0.0
0.0
0.6
0.0
0.2
0.0
0.2
(rate)
D
E
F
G
D
D==2
Fisher2 dovLL (10P)
2
min drift
0.05
2.5
0
pos.
-0.05
HD
0.0
neuron j
loghe
Z?
average FF
Fisher
([F
0-
speed
0.8
D2
2
3
0.8
1.2
D2
0
-50 5
time (8)
10
TI (rate)
time bin
Figure 3: Application to mouse head direction cells in the anterodorsal thalamic nucleus (ADn) and the postsubiculum (PoS). (A) Goodness-of-fi for Poisson (P), heteroscedastic negative binomial (hNB) and universal
(U) regression with covariates described in the text. Our method (U) outperforms baselines (p 7 10 for AcvLL w.r.t. Poisson, one-sample t-test). (B) Left: Fano factor (FF) versus mean count of the predictive distribution for time steps (scatter dots) for two representative cells. Red lines are linear regression fits to scatter points. Right: average FF across time steps and Pearson r correlation between FF and mean count for all cells, with colour indicating region ADn (orange) or PoS (green). (C) Top: conditional tuning curves (subsection D.5) of the UCM with complementary x at preferred HD, positioned at centre of arena, with zero speed, AHV and at time t = 0. Bottom: TIs w.r.t. corresponding covariates for the FF and firing rates. (D) Adding Euclidean latent dimensions D.. Top: AcvLL w.r.f. regression only. Bottom: Fisher histogram, with red curves its sampling distribution. (E) Comparison between models with latent space dimensions D= = 2 and D2 0. Top: off-diagonal noise correlations rig. with neurons ordered by region first (PoS and ADn), and then by preferred head direction within region. Bottom: average FF (same as in (B)) and Fisher Z under the posterior predictive distribution Equation 12 (F) Left: inferred latent trajectories for D: 2. Right: corresponding TIs for FF and firing rates. (G) Time scales for covariates computed from autocorrelograms (Appendix D). Note pos. refers to both and position (near-identical values). The top horizontal line is the minimum kernel length scale over absolute time across neurons. The bottom line depicts the 40 ms time bin. Error bars in A and D show s.e.m. over cross-validation runs. Shaded regions for Tks, TDs, tuning curves and latent states show the 95% CL.
Specifically, one of our covariates was absolute experimental time to capture stattiitttiiitess in neural tuning. As result, our model captured several experimental phenomena that are studied separately in the literature:
GP-based approaches typically only consider around 2 to 3 input dimensions 56 35 43
drifting neural representations [57 58 59 anticipatory time intervals (54] and conjunctive
tuning to behaviour [60]. We also applied the model in a purely latent setting similar to the example
Figure 2A, with the UCM uncovering latent signal more closely correlated to the head direction compared to baseline models. These additional results are presented in Appendixx C
8

---

4 Discussion
Related work Neural encoding model provide a statistical description of neural count activity, and typically rely parametric count likelihood, such as Poisson [33] or negative binomial [30 31 61],
that is often mismatched to empirical count statistics. Heteroscedastic versions additionally regress the dispersion parameter of count distributions to covariates [62] 63 This has shown improvements
in stimulus decoding and more calibrated posterior uncertainties [55]. Copula-based models [36
64] separate marginal distributions of single neurons from the multivariate dependency structure parameterized by the copula family, and thus remove parametric constraints on single neuron count statistics. Generally, the idea of universal model that can capture arbitrary joint distributions has been explored for binary spike trains (37] However, neither approach naturally incorporates modulation of spike count distributions by input covariates.
Our model deals with discrete spike counts ranging from 0 to in a manner similar to categorical output variables in classification, where existing GP-based approaches pass function points directly
through softmax nonlinearity [65 66). Our approach instead relies on a basis expansion and
linear-softmax mapping. At small time bins where K = 1, our model becomes identical to Bernoulli
models [67] and comparable to point process models [5T 68 69 70 In these cases, modulation of count variability becomes inseparable from firing rate modulation, making it difficult to generalize for heteroscedasticity in an interpretable manner and thus inconvenict for studying response variability. Introducing unobserved input variables incorporates aspects of Gaussian process latent variable models 771 (2). Such models have been applied to neural data to perform dimensionality reduction 33 with extensions non- oon-Euclidean latent spaces and non-reversible temporal priors 4 73
Limitations and further work The empirical choice of hyperparameters C' and basis functions based on achieving sufficient model flexibility, as confirmed with the Kolmogorov-S -miiroy approach. Recently, a multivariate extension of the Kolmogorov-Smirng test has been proposed to directly test multivariate samples against the model 74., instead of looking at single neuron statistics. A Alternatively, one could perform ARD 77 76 61 by placing a Gaussian prior on W, allowing automatic selection of relevant dimensions once basis expansion is chosen. Another avenue for future work could consider going completely non-parametric by adding count dimension to the input space, which is evaluated at counts 0 to K for every time step. This however increases the number of evaluation points by factor K + 1. For high-dimensional input stimuli common in sensory areas, deep kernels 777 provide a scalable modification of our framework. In addition, extending our model with more powerful priors for latent covariates, such Gaussian process priors [33 73 can improve latent variable analysis, especially at smaller time bins where the temporal prior influence becomes more important. Regularization methods may help decorrelate inferred trajectories 778 (9).
Conclusion and impact We introduced a universal probabilistic encoding model for neural spike count data. Our model flexibly captures both single neuron count statisties and their modulation by covariates. By adding latent variables, one can additionally capture neural correlations with potentially interpretable unobserved signals underlying the neural activity. We applied our model to mouse head direction cells and found count statistics that cannot be captured with current methods. Neural activity tends to be less variable at higher firing rates, with many cells showing both over- and underdispersion. Fano factors and mean counts generally do not show a simple relation and can even be decoupled, with Fano factor modulation comparable or in some cases even exceeding that of the rate. Finally, we found that a two-dimensional latent trajectory with timescale of around a second explained away noise correlations in these cells.
Neural variability is usually not considered on the same footing as mean firing rates, with models assigning most computational relevance to rates 80 8I). However, recent work on V1 has started to explore variability as playing a computationally well-defined useful role in the representation of uncertainty 24 25 22 26). The framework introduced in this paper provides a principled tool for empirically characterising neural variability and its modulations without the biases inherent in traditional approaches, which would likely miss potentially meaningful patterns in neural activities beyond mean rates. Our model has the potential to reveal new aspects of neural coding, and may find practical applications in designing and improving algorithms for brain-machine interfaces. As progress is made in scaling and applying such technology beyond research environments (82), it becomes increasingly more important to maintain transparency, e.g. through open source code, and to raise awareness of potential ethical issues (

---

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 212262/2/18/Z to M.L). We are grateful K.T. Jensen and Melkonyan for helpful feedback on the manuscript.
References
[1] David H Hubel and Torsten N Wiesel. Receptive fields of single neurones in the cat's striate cortex. The Journal of physiology, 148(3):574, 1959.
[2] Frédéric E Theunissen, Kamal Sen, and Allison J Doupe. Spectral-temporal receptive fields of nonlinear auditory neurons obtained using natural sounds. Journal of Neuroscience, 20(6):2315- 2331, 2000.
[3] John O' Keefe and Jonathan Dostrovsky. The hippocampus spatial map: preliminary evidence from unit activity in the freely-moving rat. Brain research, 1971.
[4] Jeffrey S Taube, Robert U Muller, and James B Ranck. dead-directions cells recorded from the postsubiculum in freely moving rats. description and quantitative analysis. Journal of Neuroscience, 10(2):420-4355 1990.
[5] Torkel Hafting, Marianne Fyhn, Sturla Molden, May-Britt Moser, and Edvard Moser. Mi- crostructure of a spatial map in the entorhinal cortex. Nature, 436(7052):801- 806, 2005.
[6] Colin Lever, Stephen Burton, Ali Jeewajee, John O'Keefe, and Neil Burgess. Boundary vector cells in the subiculum of the hippocampal formation. Journal of Neuroscience, 29(31)9771197777 2009.
[7] Michael N Shadlen and William T Newsome. The variable discharge of cortical neurons: implications for connectivity, computation, and information coding. Journal of neuroscience, 18(10):3870-3880 1998.
[8] André A Fenton and Robert U Muller, Place cell discharge is extremely variable during individual passes of the rat through the firing field. Proceedings of the National Academy of Sciences, (5(6)33122-1187 1998.
[9] Alexandre Pouget, Peter Dayan, and Richard Zemel. Information processing with population codes. Nature Reviews Neuroscience, 1(2):125-132, 2000.
[10] George J Tomko and Donald R Crapper. Neuronal variability: non-stationary responses to identical visual stimuli. Brain research, 79(3):405-418, 1974.
[11] A Aldo Faisal, Lue PJ Selen, and Daniel M Wolpert. Noise in the nervous system. Nature reviews neuroscience, 9(4):292-30 303, 2008.
[12] Johannes Nagele, Andreas VM Herz, and Martin B Stemmler. Untethered firing fields and intermittent silences: Why grid-cell discharge is so variable. Hippocampus, 2020.
[13] Justin Keat, Pamela Reinagel, R Clay Reid, and Markus Meister. Predicting every spike: model for the responses of visual neurons. Neuron, 30(3):803-811 2001.
[14] Gaby Maimon and John A Assad. Beyond poisson: increased spike-tim regularity across primate parictal cortex. Neuron, 62(3):426-440, 2009.
[15] Mark M Churchland, M Yu Byron, John Cunningham, Leo P Sugrue, Marlene R Cohen, Greg S Corrado, William T Newsome, Andrew M Clark, Paymon Hosseini, Benjamin Scott, et al. Stimulus onset quenches neural variability: a widespread cortical phenomenon. Nature neuroscience, 13(3):369, 2010.
[16] Adrián Ponce-Alvarez, Alexander Thiele, Thomas D Albright, Gene R Stoner, and Gustavo Deco. Stimulus-dependent variability and noise correlations in cortical mt neurons. Proceedings of the National Academy of Sciences, 11003)): 2013.
10

---

[17] Alexander S Ecker, Philipp Berens, R James Cotton, Manivannan Subramaniyan, George H Denfield, Cathryn R Cadwell, Stelios M Smirnakis, Matthias Bethge. and Andreas S Tolias. State dependence of noise correlations in macaque primary visual cortex. Neuron, (2(1)2335,,48 2014.
[18] Neil C Rabinowitz, Robbe Goris, Marlene Cohen, and Eero P Simoncelli. Attention stabilizes the shared gain of v4 populations. Elife, 4:c08998, 2015.
[19] Larry F Abbott and Peter Dayan. The effect of correlated variability on the accuracy of population code. Neural computation, 1 1111:91-101, 1999.
[20] Bruno B Averbeck, Peter E Latham, and Alexandre Pouget. Neural correlations, population coding and computation. Nature reviews neuroscience, 7(5):358-335 -666, 2006.
[21] Rubén Moreno-Bote, Jeffrey Beck, Ingmar Kanitscheider, Xaq Pitkow, Peter Latham, and Alexandre Pouget. Information-limiting correlations. Nature neuroscience, 17(10):1410, 2014.
[22] Guillaume Hennequin, Yashar Ahmadian, Daniel B Rubin, Máté Lengyel, and Kenneth D Miller. The dynamical regime of sensory cortex: stable dynamics around a single stimulus-tuna attractor account for patterns of noise variability. Newron, 9844))846-860, 2018.
[23] Wei Ji Ma, Jeffrey M Beck, Peter E Latham, and Alexandre Pouget. Bayesian inference with probabilistic population codes. Nature neuroscience, 911):1432-1438 2006.
[24] József Fiser, Pietro Berkes, Gergö Orbán, and Máté Lengyel. Statistically optimal perception and learning: from behavior neural representations. Trends in cognitive sciences, 14(3):11-130, 2010.
[25] Gergö Orbán, Pietro Berkes, József Fiser, and Máté Lengyel. Neural variability and sampling- based probabilistic representations in the visual cortex. Neuron, 92(2):530-54) 2016.
[26] Rodrigo Echeveste, Laurence Aitchison, Guillaume Hennequin, and Máté Lengyel. Cortical-like dynamics in recurrent circuits optimized for sampling-based probabilistic inference. Nature Neuroscience, 23(9): 138-1111111 149, 2020.
[27] David Tolhurst, Anthony Movshon, and Andrew F Dean. The statistical reliability of signals in single neurons in cat and monkey visual cortex. Vision research, 2378)7775-785, 1983.
[28] Robbe LT Goris, J Anthony Movshon, and Eero Simoncelli. Partitioning neuronal variability. Nature neuroscience, 17(6):858, 2014.
[29] Jonathan W Pillow, Jonathon Shlens, Liam Paninski, Alexander Sher, Alan M Litke, EJ Chichilnisky, and Eero P Simoncelli. Spatio-temporal correlations and visual signalling in a complete neuronal population. Nature, 454(7207):995-999 2008
[30] Jonathan Pillow and James Scott. Fully bayesian inference for neural models with negative- binomial spiking. Advances in neural information processing systems, 25:1898-1999 2012.
[31] Yuanjun Gao, Lars Busing, Krishna V Shenoy, and John P Cunningham. High-dimensional neural spike train analysis with generalized count linear dynamical systems. In Advances in neural information processing systems, pages 2044-2052, 2015.
[32] Byron M Yu. John P Cunningham, Gopal Santhanam, Stephen Ryu, Krishna V Shenoy, and Maneesh Sahani. Gaussian-process factor analysis for low-dimensional single-trial analysis of neural population activity. Advances in neural information processing systems, 21:1881- 1111881--888 2008.
[33] Anqi Wu, Nicholas A Roy. Stephen Keeley, and Jonathan W Pillow. Gaussian process based nonlinear latent structure discovery multivariate spike train data. In Advances neural information processing systems, pages 3496-3505, 2017.
[34] Ryan J Low, Sam Lewallen, Dmitriy Aronov, Rhino Nevers, and David W Tank. Probing variability in a cognitive map using manifold inference from neural dynamics. bioRxiv, page 418939, 2018

---

[35] Chethan Pandarinath, Daniel J O' Shea, Jasmine Collins, Rafal Jozefowicz, Sergey D Stavisky, Jonathan C Kao, Eric M Trautmann, Matthew Kaufman, Stephen Ryu, Leigh R Hochberg. et al. Inferring single-trial neural population dynamics using sequential auto-encoders. Nature methods, 15110)8005881, 2018.
[36] Pietro Berkes, Frank Wood, and Jonathan Pillow. Characterizing neural dependencies with copula models. Advances in neural information processing systems, 21:129-1360 2008
[37] Il Memming Park, Evan W Archer, Kenneth Latimer, and Jonathan W Pillow. Universal models for binary spike patterns using centered dirichlet processes. In Advances in neural information processing systems, pages 2463-2471, 2013.
[38] Michalis Titsias. Variational learning of inducing variables in sparse gaussian processes. In Artificial Intelligence and Statistics, pages 567-574, 2009.
[39] James Hensman, Nicolo Fusi, and Neil Lawrence. Gaussian processes for big data, arXiv preprint ar.iv::309.0.... 6835, 2013.
[40] James Hensman, Alexander Matthews, and Zoubin Ghahramani. Scalable variational gaussian process classification. In Artificial Intelligence and Statistics, pages 351-360. PMLR, 2015.
[41] Adrien Peyrache and György Buzsáki. Extracellular recordings from multi-site silicon probes in the anterior thalamus and subicular formation of freely moving mice. CRCNS, 2015.
[42] Adrien Peyrache, Marie M Lacroix, Peter C Petersen, and György Buzsáki. Internally organized mechanisms of the head direction sense. Nature neuroscience, 18(())569-57, 2015.
[43] Kristopher Jensen, Ta-Chu Kao, Mareo Tripodi, and Guillaume Hennequin. Manifold gplvms for discovering non-euclidean latent structure in neural data. Advances in Neural Information Processing Systems, 33, 2020.
[44] Jakob H Macke, Lars Buesing, John P Cunningham, Byron M Yu, Krishna V Shenoy, and Maneesh Sahani. Empirical models of spiking in neural populations. In Advances in Neural Information Processing Systems 24: 25th conference on Neural Information Processing Systems (NIPS 2011), pages 1350-1358, 2012.
[45] Luca Falorsi, Pim de Haan, Tim R Davidson, and Patrick Forré Reparameterizing distributions on lie groups. arXiv preprint arXiv:1903.00 02958, 2019.
[46] Matthew D Hoffman, David M Blei, Chong Wang, and John Paisley. Stochastic variational inference. The Journal of Machine Learning Research, 14(1):33331134, 2013.
[47] Diederik P. Kingma and Max Welling. Auto-encoding variational bayes. In 2nd International Conference on Learning Representations, 2014.
[48] Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv: 1412,690,, 2014.
[49] Adam Kohn and Matthew A Smith. Stimulus dependence of neuronal correlation in primary visual cortex of the macaque. Journal of Neuroscience, (1114)6661-36773 2005.
[50] I-Chun Lin, Michael Okun, Matteo Carandini, and Kenneth D Harris, The nature of shared cortical variability. Neuron, 87(3):644-556, 2015.
[51] Jonathan W Pillow. Time-rescaling methods for the estimation and assessment of non-poisson neural encoding models. In Advances neural information processing systems, pages 1473-1481, 2009,
[52] Adam S Charles, Mijung Park, J Patrick Weller, Gregory D Horwitz, and Jonathan W Pillow. Dethroning the fano factor: a flexible, model-based approach to partitioning neural variability. Neural computation, 3014):0012-1045, 2018.
[53] Rishidev Chaudhuri, Berk Gerçek, Biraj Pandey, Adrien Peyrache, and lla Fiete. The intrinsic attractor manifold and population dynamics of a canonical cognitive eircuit across waking and sleep. Nature neuroscience, 1512::152-15200 2019.
12

---

[54] Johannes Zirkelbach, Martin Stemmler, and Andreas VM Herz. Anticipatory neural activity improves the decoding accuracy for dynamic -radddirettion signals. Journal of Neuroscience, 39(15):2847-2859 2019.
[55] Abed Ghanbari, Christopher M Lee, Heather L Read, and Ian H Stevenson. Modeling stimulus- dependent variability improves decoding of population neural responses. Journal of Neural Engineering, 16(6):066018, 2019.
[56] Kamiar Rahnama Rad and Liam Paninski. Efficient, adaptive estimation of two-dimensional firing rate surfaces via gaussian process methods. Network: Computation in Neural Systems, 21 (1334))4422-68 2010.
[57] Yaniv Ziv, Laurie D Burns, Eric D Cocker, Elizabeth O Hamel, Kunal K Ghosh, Lacey Kitch, Abbas El Gamal, and Mark J Schnitzer. Long-term dynamics of cal hippocampal place codes. Nature neuroscience, 16(3):264, 2013.
[58] Michael E Rule, Adrianna R Loback, Dhruva V Raman, Laura N Driscoll, Christopher D Harvey, and Timothy O'Leary. Stable task information from unstable neural population. Elife, 9:e51 121, 2020.
[59] Daniel Deitch, Alon Rubin, and Yaniv Ziv. Representational drift in the mouse visual cortex. bioRxiv, 2020.
[60] Francesca Sargolini, Marianne Fyhn, Torkel Hafting, Bruce L McNaughton, Menno P Witter, May-Britt Moser, and Edvard Moser. Conjunctive representation of position, direction, and velocity in entorhinal cortex. Science, 312(5774) 2006.
[61] Kristopher T Jensen, Ta-Chu Kao, Jasmine Talia Stone, and Guillaume Hennequin. Scalable bayesian gpfa with automatic relevance determination and discrete noise models. bioRxiv, 2021.
[62] Seth D Guikema and Jeremy Goffelt. A flexible count data regression model for risk analysis. Risk Analysis: An International Journal, 28(1):213-223, 2008
[63] Kimberly F Sellers and Galit Shmueli. A flexible regression model for count data. The Annals of Applied Statistics, pages 943-961, 2010.
[64] Meng Hu, Kelsey Clark, Xiajing Gong, Behrad Noudoost, Mingyao Li, Tirin Moore, and Hualou Liang. Copula regression analysis of simultaneously recorded frontal eye field and inferotemporal spiking activity during object-based working memory. Journal of Neuroscience. (5123))8745-875) 2015.
[65] Kian Ming A Chai. Variational multinomial logit gaussian process. The Journal of Machine Learning Research, 13:1745-1808, 2012.
[66] Yarin Gal, Yutian Chen, and Zoubin Ghahramani. Latent gaussian processes for distribution estimation of multivariate categorical data. In International Conference on Machine Learning, pages 645-654, PMLR, 2015.
[67] Robert Haslinger, Gordon Pipa, and Emery Brown. Discrete time rescaling theorem: determin- ing goodness of fit for discrete time statistical models of neural spiking. Neural computation, 22(10):4477-2566 2010.
[68] Shuai Xiao, Junchi Yan, Xiaokang Yang. Hongyuan Zha, and Stephen Chu. Modeling the intensity function of point process via recurrent neural networks. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 31, 2017.
[69] Takahiro Omi, Naonori Ueda, and Kazuyuki Aihara. Fully neural network based model for general temporal point processes. arXiv preprint 1rXiv:190S.09691 2019.
[70] Oleksandr Shchur, Nicholas Gao, Marin Bilos, and Stephan Günnemann. Fast and flexible temporal point processes with triangular maps. Advances in Neural Information Processing Systems, 33, 2020.
[71] Neil D Lawrence. Gaussian process latent variable models for visualisation of high dimensional data. In Nips, volume 2, page 5, Citeseer, 2003.
13

---

[72] Michalis Titsias and Neil D Lawrence. Bayesian gaussian process latent variable model. In Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, pages 844 851, 2010.
[73] Virginia Rutten, Alberto Bernacchia, Maneesh Sahani, and Guillaume Hennequin. Non- reversible gaussian processes for identifying latent dynamical structure in neural data. Advances in Neural Information Processing Systems, 33, 2020.
[74] Michael Naaman. On the tight constant in the multivariate dvoretzky ity. Statistics & Probability Letters, 173:10908$ 2021.
inequal-
[75] Christopher M Bishop. Bayesian pea. Advances in neural information processing systems, pages 382-388, 1999.
[76] Andreas Damianou, Carl Henrik Ek, Michalis K. Titsias, and Neil D. Lawrence. Manifold relevance determination. In Proceedings of the 29th International Conference on Machine Learning, ICML 2012, Edinburgh, Scotland, UK, June 26 July 1, 2012. ieml.ee Omnipress, 2012.
[77] Andrew Gordon Wilson, Zhiting Hu, Ruslan Salakhutdinov, and Eric P Xing. Deep kernel learning. In Artificial intelligence and statistics, pages 370-378. PMLR, 2016.
[78] Irina Higgins, Loic Matthey, Arka Pal, Christopher Burgess, Xavier Glorot, Matthew Botvinick, Shakir Mohamed, and Alexander Lerchner. beta-vae: Learning basic visual concepts with constrained variational framework. In 5th International Conference on Learning Representations, ICLR 2017, Toulon, France, April 24-26, 2017, Conference Track Proceedings. OpenReview.net, 2017.
[79] Ricky TQ Chen, Xuechen Li, Roger Grosse, and David Duvenaud. Isolating sources of disentanglement in vaes. In Proceedings of the 32nd International Conference on Neural Information Processing Systems, pages 2615-2625, 2018
[80] Guillaume Hennequin, Tim P Vogels, and Wulfram Gerstner. Optimal control of transient dynamics in balanced networks supports generation of complex movements. Neuron, 82(6))1394 1406, 2014.
[81] Klaus Wimmer, Duane Q Nykamp, Christos Constantinidis, and Albert Compte. Bump attractor dynamics in prefrontal cortex explains behavioral precision in spatial working memory. Nature neuroscience, 17(3):331-439 2014.
[82] Elon Musk et al. An integrated train-machine interface platform with thousands of channels. Journal of medical Internet research, 21(10)xe1619, 2019.
[83] Jens Clausen. Man, machine and in between. Nature, (08--223)))00000000000000 2009.
14


Reviews:
[{"review_id": "yBlHOjyQ7Ud", "reviewer": "Reviewer_oo8R", "summary": "This study proposes a statistical model to jointly estimate the mean firing rate and covariance structure based on the spike count of neural populations, and test the method on mouse head-direction cells. ", "questions": "", "limitations": "", "rating": 6, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper proposes a new descriptive statistical model to fit neural data. By using Gaussian process, the model is able to jointly fit mean firing rate (1st order statistics) and covariance (2nd order statistics) simultaneously, which is  the main novelty of this paper. And the author used this model to fit synthetic data and real neural data, and compared the model performance with other alternative models used in neural data analysis. The study is solid, and figures look clean and great. However, since some key math notations are not defined yet, I am not sure the math used in this study is correct or not although I hope so.\n\n### Writing:\nThe writing and the definition of math notations can be improved a lot in a revised version. For example, some underlying motivations of proposing a math equation are lacking, and some math notations are not defined at all when using.\n\n- I don’t understand what are observed X and latent Z, and what is the difference between the two variables. From the generative model shown in Fig. 1, the observed X seems to be a latent variable generating recorded spike count Y. Or X is regarded as the presented stimulus to an animal? Please give a concrete example of what X and Z can be referred to in neurophysiology experiments and otherwise it is not easy to understand.\nis the motivation to include those two variables.\n\n- Line 100: It is not unclear what $\\Pi$ is. It is not defined. The $\\Pi$ seems a lumped variable including all observations but I am not sure.\n- Eq.8: what is \\phi^{-1}? it is not defined.\n- Eq.11: it is not clear what F(q) is. What is the difference between $F(q)$ and $q$? Is $F(q)$ the empirical distribution while $q$ the parametric distribution defined in Eq. 2?\n- Eq. 12: please add the motivation of including $1/T+ 1/3T^2$ in the equation. \n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes", "review_text": "", "extra_scores": {"ethical_concerns": "No", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2.5"}, "extra_text": {}}, {"review_id": "okg2MtC1ey", "reviewer": "Reviewer_Saxk", "summary": "In this paper authors provide a method for modelling spike counts as a function of input covariates and latent covariates. A Markovin prior is assumed over latent parameters, and the effect of latent parameters and input covariates on spike counts are modelled using a Gaussian process. The authors also suggest a method for calculating Z-scores and tuning index based on their method. The results include analyses of a synthetic dataset and an experimental data, which generally shows the proposed method outperforms considered baselines. ", "questions": "", "limitations": "", "rating": 6, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper is generally clear (but see below). The statistical model is interesting and suitable for modelling the joint effect of latent and observed covariates on the neural activities. (the main text needs to be updated based on the typo mentioned in the supplementary materials). The results are also encouraging and show that the model is improving on baseline methods.\n\nI have the following concerns about the study that require further improvements:\n\nSome of the results are presented in a rather confusing way. For example cvLL results are presented relative to the Poisson model in some panels (e.g,, Fig 2A) but not in other places (e.g., Fig 3 in supple). In Fig 2A for example it's hard to see whether differences with the performance of the Poisson model is significant or not. Also RMSE is reported for synthetic experiments, but not for the experimental data.\n\nThe baselines are also different between the experiments; the performance of ANN is not reported for the experimental data. Also, I couldn’t find how the architecture of the neural networks was selected; how many neurons? How were hyper-parameters tuned? This is particularly important since recent approaches report promising results using recurrent neural networks, e.g. below, and the reflective performance of the suggested method w.r.t.to these methods is unclear.\n\nPandarinath, Chethan, et al. \"Inferring single-trial neural population dynamics using sequential auto-encoders.\" Nature methods 15.10 (2018): 805-815.\n\nSections 2.4 and 2.5 are rather hard to follow. For example it would be good to expand why the definition of Z in eq 8 is a sensible one and how it is related to the common definition of z score. Below eq 11, for example it's unclear what does data being more/less regular than the data mean.\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "They are adequately discussed.", "review_text": "", "extra_scores": {"ethical_concerns": "Please include relevant ethics information regarding the experimental data.", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}, {"review_id": "gWjs5b9fJzr", "reviewer": "Reviewer_S8Uc", "summary": "The authors present a very broad framework to model multi-neuronal spike trains (multivariate time series over the nonnegative integers). The model class is sufficiently broad that it can express any arbitrary distribution over N-dimensional vectors containing the integers {0, ..., K}, and thus model the activity of N neurons within a single time bin (assuming an upper bound on the number of spikes, K).", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "This paper is well-written, the figures look polished, and I believe the core claims of the paper are correct. Overall, I would like to see this paper published if the authors can address my concerns below, which mostly deal with the interpretation of the model, relation to existing methods, and framing of the paper. In short, I think the model could be useful in cases where one would like to specify large time bin sizes, but the utility of the method is less clear in the common case where one takes the time bins to be very short.\n\nMost existing methods, the authors claim, have a \"fundamental underlying assumption of Poisson spiking statistics.\" The model they propose does not have this assumption if time bins are large and K > 1. However, if we consider the limit of very small time bins (on the order of 1 millisecond) then it is reasonable to set K=1, since the probability of a neuron firing more than one spike is negligible. The multinomial distribution in thier model becomes a binomial distribution. In the limiting case of of time bin duration going to zero, you recover Poisson statistics (see, e.g., sec 5.2 of Kass et al. \"Analysis of Neural Data\"). In this limit, the model proposed by the authors looks very similar to existing models, like GPFA or GLMs.\n\nIf the above analysis is correct, then I think the authors should acknowledge it in the paper and explain that their method is a generalization that retains \"universal\" modeling capabilities when time bins are not very small and when K > 1. Then, they should explain what the relative advantages and disadvantages are to choosing larger time bins.\n\nThere are other parts of this paper that I really like, such as the combination of supervised learning and latent variable modeling to capture representational drift and head direction tuning. However, these modeling innovations don't seem to be sold as the key results of the paper and it doesn't seem like the \"universal\" spike count model is really necessary to enable this analysis. Again, it seems like one could choose small time bins and use GPFA (w/ poisson noise) + GLMs to arrive at a similar result.\n\nOverall, I thought this paper was very well-written and presented some elegant analysis of real experimental data, so I have scored it as a marginal accept. The modeling innovations that are most heavily emphasized in the manuscript seem to be dependent on having large time bin sizes. If my interpretation here is wrong I am happy to raise my score higher.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The limitations of the model are clear and I do not foresee any negative societal impacts.", "review_text": "", "extra_scores": {"ethical_concerns": "n/a", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}, {"review_id": "YODdA_iMqr", "reviewer": "Reviewer_XYHL", "summary": "The paper presents a model of spike count distributions in a population of neurons. The model’s spike count distribution depends both on stimulus information and latent variables through a Gaussian process to capture a large range of joint spike count distributions (as the number of factors in the model increases). Moreover, the authors present a set of model fit statistics for evaluating the model fit, and how correlations factor into the model’s predictions.", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The model improves upon previous non-parametric approaches (such as the universal binary count model) to modeling spike count distributions in neural population recordings: importantly, the proposed model includes stimulus information and temporal correlations. The example applications are very thorough in examining the model fit and comparisons to the Poisson and negative binomial models.\n\nThe presentation was generally clear. However, I am not convinced I understand the procedure for quantifying correlations given on lines 183-190. I’m certain I’m overthinking things here, but I do think this one portion of the methods could be presented more clearly.\n\nMinor point: one limiting factor not mentioned in the discussion is that the stimulus term in the examples here are low-dimensional compared to stimuli in sensory systems.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2.5"}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper presents a universal spike-count model for flexibly describing both over-dispersion and under- dispersion in neural spiking data. The model is presented with an efficient variational fitting method, and can be applied to data without known stimuli, as is often the case for overdispersion metrics. The reviewers recognized the novelty of the model, and the topic fits nicely in the probabilistic modeling and neuroscience focuses in NeurIPS. Moreover spike-count models are an ongoing discussion in the computational neuroscience community (heavily represented at NeurIPS) and this work presents a new and interesting addition to the literature. I therefore recommend this work be accepted in NeurIPS.
Decision: Accept (Poster)