Bayesian Clustering of Neural Spiking Activity Using a Mixture of Dynamic Poisson Factor Analyzers
Authors: Ganchao Wei
Ian Stevenson
Xiaojing Wang
Venue: NeurIPS 2022 Accept (2022)
Keywords: neural spike data
clustering
dynamic Poisson factor model (DPFA)
mixture of finite mixtures model
Markov chain Monte Carlo(MCMC)
Laplace Approximation

Abstract: Modern neural recording techniques allow neuroscientists to observe the spiking activity of many neurons simultaneously. Although previous work has illustrated how activity within and between known populations of neurons can be summarized by low-dimensional latent vectors, in many cases what determines a unique population may be unclear. Neurons differ in their anatomical location, but also, in their cell types and response properties. Moreover, multiple distinct populations may not be well described by a single low-dimensional, linear representation.	To tackle these challenges, we develop a clustering method based on a mixture of dynamic Poisson factor analyzers (DPFA) model, with the number of clusters treated as an unknown parameter. To do the analysis of DPFA model, we propose a novel Markov chain Monte Carlo (MCMC) algorithm to efficiently sample its posterior distribution. Validating our proposed MCMC algorithm with simulations, we find that it can accurately recover the true clustering and latent states and is insensitive to the initial cluster assignments. We then apply the proposed mixture of DPFA model to multi-region experimental recordings, where we find that the proposed method can identify novel, reliable clusters of neurons based on their activity, and may, thus, be a useful tool for neural data analysis.

Bayesian Clustering of Neural Spiking Activity Using a Mixture of Dynamic Poisson Factor Analyzers
Ganchao Wei Department of Statistics University of Connecticut ganchao.. wei @uconn. edu
Ian. H. Stevenson
Department of Psychological Sciences
University of Connecticut
ian. stevenson@uconn.cmm edu
Xiaojing Wang
Department of Statistics University of Connecticut xiaojing wang@uconn. edu
Abstract
Modern neural recording techniques allow neuroscientists to observe the spiking activity of many neurons simultaneously. Although previous work has illustrated how activity within and between known populations of neurons can be summarized low-dimensional latent vectors, in many cases what determines unique popu- lation may be unclear. Neurons differ in their anatomical location, but also, in their cell types and response properties. Moreover, multiple distinct populations may not be well described by single fow-dimensional, linear representation. To tackle these challenges, we develop a clustering method based mixture of dynamic Poisson factor analyzers (mix DPFA) model, with the number of clusters treated as unknown parameter. To do the analysis of DPFA model, we propose a novel Markov chain Monte Carlo (MCMC) algorithm to efficiently sample its posterior distribution. Validating our proposed MCMC algorithm with simulations, we find that it can accurately recover the true clustering and latent states and is insensitive to the initial cluster assignments. We then apply the proposed mixDPFA model to multi-region experimental recordings, where we find that the proposed method can identify novel, reliable clusters of neurons based on their activity, and may, thus, be useful tool for neural data analysis.
Introduction
With modern high-density probes (Jun et al, 2017), neuroscientists can observe the spiking activity of many neurons from many different anatomical regions simultaneously. With these expanding capabilities, new methods to analyze neural data at the population-level and at the level of multiple populations become necessary. Several recent models have been developed to extract shared latent structures from simultaneous neural recordings, assuming that neural activity can be described through low-dimensional latent states. Many existing approaches are extensions of two basic models: the linear dynamical system (LDS) model (Macke et al. 2011) and a Gaussian process factor analysis (GPFA) model (Yu et al. 2009). The LDS model is built on the state-space model and assumes latent factors evolve with linear dynamics. On the other hand, GPFA models the latent vectors by non non-paramen parammrrcttic Gaussian processes. However, in both cases, the observation model is generalized linear. Several variants of these models have been implemented to analyze multiple neural populations and their interactions (Semedo al.) 2019 Glaser et al. 2020). However, in many cases, the total number of distinct populations and which neurons belong population unclear.
36th Conference on Neural Information Processing Systems (NcurIPS 2022).

---

Neurons in different anatomical locations may interact with each other or receive common input from unobserved brain areas, sharing the same latent structure. On the other hand, neurons of different cell-types within the same brain area may be better described by distinct latent structures. From a functional point of view, neither the anatomical location nor cell type (Fig. [TA indicates which neurons should be grouped into the same populations. The incorrect population assignments can lead to biased and inconsistent inference on the latent structure (Ventura 2009). If we instead ignore multi-population structure and treat all neurons single population, then using linear model based methods may not describe their activity well, especially when the input is nomooooggsssusss Besides, nonlinear models such deep learning (Pandarinath et al. 2018 W Witeway et al., (0199) and Gaussian processes (Wu et al. 2017) have been developed, but these models do not explicitly distinguish among distinct populations of neurons.
Motivated by the mixture of (Gaussian) factor analyzers (MFA Arminger et al. 1999 Ghahramani and
Hinton 1996. Fokoué and Titterington 2003). which describes globally nonlinear data by combining a number of local factor analyzers, here we group neurons based on the latent factors (Fig. BB). A
similar idea was previously implemented using mixture of Poisson linear dynamical system (PLDS)
model (mixPLDS, Buesing et al. 2014). The mixPLDS model infers the subpopulations and latent factors using deterministic variational inference Wainwright and Jordan (008)) Jordan et al. (1999); ;
Emtiyaz Khan et al. (2013) and the model parameters are estimated by Expectation Maximization (EM). Unlike MFA, the mixPLDS can capture temporal dependencies of neural activity as well as interactions between clusters over time. However, there are several limitations for mixPLDS: 1) requires we predetermine the number of clusters, and 2) the clustering results often sensitive to the initial cluster assignment.
Here we cluster the neurons by a mixture of dynamic Poisson factor analyzers (mixDPFA). The DPFA model takes the advantages of both Poisson factor analysis (FA) and PLDS and includes both population baseline and baselines for individual neurons. The number of clusters is treated as an unknown parameter in the mixDPFA, and the posteriors are sampled using Markov Chain Monte Carlo (MCMC). To sample high dimensional latent factors, we approximate the full conditional distribution of the latent state by a Gaussian, which is similar to results by sampling from exact full conditional distribution. To improve mixing in the cluster assignments, we marginalize the loading out for clustering by Poisson-Gamm for successful sampling of the proposed models. After validating the proposed model with simulated data, we apply it to analyze multi-region experimental recordings from behaving mice: the Visual Coding Neuropixels Dataset from the Allen Institute for Brain Science. Overall, the proposed method provides way to efficiently cluster neurons into populations based on their activity.
conjugacy. We also discuss the constraints necessity
2 Methods
Here we introduce a mixture of dynamic Poisson factor analyzers (mixDPFA) to cluster neurons based on multi-population latent structure. The number of mixture components is treated as an unknown parameter and the posteriors are sampled by MCMC. In this section, we first provide the single population DPFA for a given cluster. Then, we introduce prior on the number of clusters and describe how we use the mixture of finite mixture model (MFM) to efficiently sample the posterior of the miXDPFA.
2.1 Dynamic Poisson Factor Analyzer
Denote the observed spike count of neuron E [I,...,N] NY at time bin 1 E (l,..., T} as you (a non-negative integer), and let y: (yar,..., (er)). Further, let be the cluster indicator of neuron i. Motivated by the nature of neural activity and the former PLDS model (Macke et al., 2011). propose new Poisson FA model by adding individual baselines ő,. The proposed model is
a combination of PLDS and Poisson FA, which includes both population baseline and individual baseline. Assume neuron belongs to the j-th cluster (i.c., 2 = j), and its spiking activity is
independently Poisson distributed, conditional on the low-dimensional latent state n..") E RPs and population baseline p,(") as follows:
Di ~ Poi(Aa))
log Are = 6, + ub + caa

---

A
anatomy
cell type
C
fx
B
pop. baseline
latent X
1
spikes
1 x
11
Yi Duu
6,
Pr
I....
Figure 1: Model overview. A. There are multiple potential ways to define neural populations. For instance, populations could be defined by anatomical regions (left) or by cell types (right). Since the same latent structure could be shared across anatomical sites and cell types, useful alternative may be define populations based on neural activity directly. B. The main goal for the proposed method is to cluster neurons according their activity and extract functional grouping structure, based on spike train observations. The activity of each neuron is determined by a low dimensional latent state, specific that neuron's 's cluster assignment (e.g. yellow, red, blue). C. Graphical representation of the mixture of finite mixtures (MFM) of dynamic Poisson factor analyzers (DPFA) generative model.
Here the cluster number is treated random variable. The population baseline (4.2)) and the latent factor (x)))) for each cluster is generated by linear dynamics, with a Gaussian noise.
with c; ~ N (0, I, The neuron-specific baseline & is constant across time for the ith neuron and unrelated the cluster assignment. For simplicity, we assume the dimension of latent factors (states) is the same for all clusters, s.t. Pj = however, our method easily extends to the situation when Pys differs across clusters (see Discussion). Further, we assume the population baseline (t)
and the latent state x. (5) evolve lincarly time with Gaussian noise following
P 90 h1) + G..
(1) B(5) +A A0)x((0) + ni (
x+
where e ~ N'0, 26)) and m
(i) N0,00Q(1)) Q())..
If denote X; (An, A proposed model can be rewritten as
()) (v( :....,I 1 and X(i (x
ay Llly the
Y: Poi(>:),
(1)
log A 6,17 + p (i) + x (.)c.
Generally, factor model is consistent only when T/N- 0 (Johnstone and Lu 2009). but this is often not the case for most neural spike data. However, when we assume linear dynamics on (,()) and X (i), it resolves the consistency issue. As known in a FA model, when > 1, the model is only identifiable up to orthogonal rotation on X() with e; ~ N(0, Ip)- With including an individual baseline & 1r in proposed DPFA model (T) it further makes the model invariant translation of (i) and ( That means if ( and ( is set solution, then ( + aly and X X & m' also satisfy the model, for any a, m and orthogonal matrix U. Thus, to make the model identifiable, need to add several constraints. Although the clustering is invariant to orthogonal rotation, how put constraints translation will influence the cluster assignments. Here assume A and

---

Q are diagonal (Peña and Poncela 2004 Lopes et al. 2008). and, to encourage clustering based on the trajectories of latent factors, we set E Mi 0 and Σ x 0. See Section 5 for more discussions about the choice of these constraints.
Given the parameters of the j-th cluster ((i) {ul) . 1 ggD ,0 A B0) Q }. the spike counts of neuron i are generated by the dynamic Poisson factor analyzer (DPFA) model as We 2; j DPF A C,, ).))). To faciliate the Bayesian computation, we have to impose priors on 0 see more details of prior settings in Appendix.
2.2 Clustering by Mixture of Finite Mixtures Model
When the population labels 2,5 are unknown, we cluster the neurons by a mixture of DPFA (mixDPFA). Since the number of neural populations finite but unknown, we need to put priors on To make the Bayesian computation more efficient, we utilize the idea from the mixture of finite mixtures (MFM, Miller and Harrison 2018) model, by assigning the priors for the clusters in the following way:
k ~ fe. Dire(s.....
fa is p.m.f. on given given n.
(T1,
TTp)
i.i.d. i.i.d
π 21 2N
(2)
(())
0 (kk)
given k,
(Bi1...., V DPFA(8, c,, (())))))) given 8, C, (())) vi = I....,N,
where p. denotes the probability mass function. By using the MFM, we can integrate the field knowledge about the number of neural populations into our analysis. In the analysis of this paper, we assume & follows a geometric distribution, i.e., k ~ Geometric(a) with its density defined as (kla)) (1 a a for & 1, and let y 1. The complete generative model is summarized in a graphical form shown in Fig. FIg.
2.3 Inference
Here the posteriors of the proposed mixDPFA model sampled by MCMC algorithm (see Appendix I.I). In each iteration, we sample the model parameters assuming the known cluster indices first, and then sample the cluster indices given the model parameters. When sampling the (labeled) model parameters, the latent state X)) and population baseline (كبر have closed-for full conditional distributions. Although we could sample the posterior by particle MCMC directly, convergence may be too slow for clustering. Here, we approximate the full conditional distribution for X and (ن) by Gaussian distribution (a Laplace approximation) and generate samples according to this approximation. This Laplace approximation is widely used for EM (Macke et al., (011)) and variational inference (Glaser al. 2020) with PLDS models and their variants, and Gaussian approximation of the intractable full conditional distribution of latent effects has also been used in Bayesian mixed effects binomial regression, where Berman et al. 2022 found that the approximation provided reasonable estimation accuracy with substantial computational speedups. We, thus, use a global I Laplace approximation that be efficiently computed in ((T) (Paninski et al. 2010).. To help convergence, sampling on (labeled) model parameters is repeated several times before updating the cluster indices. Approximating the intractable full conditional with a Laplace approximation also makes computation of the proposed mixDPFA more efficient. However, to assess the accuracy of the approximation we also compare our approach to directly sampling from the exact posterior of the model. We develop a Polya- (PG) data augmentation approach (Windle et al. 2013
Linderman al. 2017 2016 Polson al. 2013) with an additional Metropolis-
(MH) step (Metropolis et al. 1953 Hastings 1970 sample exactly from the full conditional of X ()) and fe(o).
We find that the proposed method using a Laplace approximation is faster but performs similarly as sampling from the exact joint posterior (Fig. Fig. [ Appendix A.2, and Table (D).
Once we update the latent state (i) and population baseline fe the cluster index is then sampled by the analogy of partition- algorithm Dirichlet process mixtures (DPM, Neal 2000). See details in Miller and Harrison 2018 and the Appendix [). When doing the clustering, we need to evaluate the likelihood for neurons under each cluster. Although we can sample e; directly and evaluate the full likelihood as in MCMC for Gaussian MFA (data-augmentattion/ imputation-poster

---

algorithm, Fokoué and Titterington 2003). the chain has poor mixing and stops after a few iterations, because of the high dimensionality. The heavy dependency on the starting point when fitting the mixture of PLDS (mixPI.SS, Buesing et al. (2014) model may suggest similar prohlem. To resolve this, we evaluate the marginal likelihood by integrating out the neuron-specifie C;, i.e., the marginal likelihood of neuron i in cluster j is computed by
Mau (3:) = P(y ,6.) - P(y 6,, c.)
(3)
However, this marginal likelihood has closed form. Though may evaluate by Laplace approximation, but iterating over all potential clusters for each neuron is computationally intensive. To make faster clustering, we approximate the marginal likelihood by utilizing Poisson-Gamma conjugacy. This approach has been previously utilized to approximate posteriors (EI-Sayyad 1973)
and predictive distributions (Chan and Vasconcelos 2009). In our situation, since e, N (0, 1p).
lognormall 6; + p4 (x("))) x(1)) xx")), and then we can approximate this lognormal distribution by gamma distribution, i.e., assume X follows Gamma (a., bae) with
have X = cxp(6; +4 + tc,x,
axt (ax (6) ) and bye ((3) (1)
=()))(()) .eitrl" Then, by the conjugate property with Poisson and
Gamma random variables, we have
P(uul)a))P(i))) dha z NB(yalaa Psc).
with V/2 = au and Pa = 1/(1 + bee). Further, noticing that we have the conditional independence assumption for P (()), 8,), that is (0), (1)6),66)))))))) 6,) = I P(ys (0()) 52). then have a closed-
form for Equation 3
Another possible idea is to approximate the log- le-likeliioodd by second-order polynomials, with coefficients determined by Chebyshev polynomial approximation (Keeley et al. 2020). However, we find that this approximation doesn't work well in practice when spike counts have a wide range. The model is implemented MATLAB and the code available /https://gt tht.g.thhuu...... com/we wwi/ww/gddsb////// MFM DPFA clean
3 Simulations
To validate and illustrate the proposed clustering method, we simulate neural data directly from the generative model (1) The labels for each neuron are assumed known and fixed at first check convergence and model identifiability. We then infer the labels to evaluate clustering performance. All experiments in this paper were performed using 3.40 GHz processor with 16 GB of RAM.
3.1 Labeled data
We first simulate clusters with neurons in each, with recording length T = 1000 and p = dimensional latent factors for each cluster. Individual baselines generated by 5 N(0, 0.5²). and the loading for the latent states generated by c; N N(0, I2). The population baseline pu(i) and latent vector (i) are generated by the spline interpolation on 10 to 30 evenly spaced knots. The simulations are conducted 50 times with different seeds. Here we show results for one simulation, and the performance for the rest is similar.
Since the labels are known, each whole simulation is equivalent to 10 independent simulations, with 5 neurons in each. Running MCMC for 10,000 iterations, find that the log- likeliiod per spike converges rapidly for individual clusters and overall (Fig2A) Trace plots (Fig 2 of the (Frobenius) norms for linear dynamics samples (h:)), (3), 200) A (7) b()) (()) show rapid mixing and convergence for each DPFA. The fitted mean firing rate (mean response) (Fig. 2C), (5) and x() (Fig. 2D) match the ground truth well. The convergence is fast, especially terms of mean response and population baseline (6) Together, these results demonstrate the identifiability of the DPFA with appropriate constraints.
We then compare the 10-cluster model with the simplified model ignoring the clustering structure (1-cluster model). The p for 1-cluster model 14, chosen by 5-fold speckled -ras--valiaation described in (Williams al. 2020). To evaluate the fitting performance of these two models,
5

---

hold out 1/2 of the data and compare the held-out log- logllikelihood per spike. The distribution training log- likelihood is evaluated by averaging over the samples in short chains (e.g. iteration 50 to 100), which is justified by the observed fast convergence. The same procedure is replicated for 50 times. In this case, ignoring the clustering structure leads to worse performance FF) and the single population analysis cannot describe the data as accurately, since the input is non-homogeneous and the data is global nonlinear.
A
Cluster 3,
overall
B
[dllggt. All
ll.. 2/11,
[|diaglo", all
p=2
p-14
duserr duster
0.02
0.02
log-hhu Froood splk
1.8
hed ouu thei spae
0.01
-1 1.000
1.7
0.025
-1.01
-1.1-
duater
0.01
0.02
-1.00
-1.2%
5000
10000
5000
10003
502D
10000
c0-sketerr
1-dlrrrr
D
true
Aug. (9K 10K)
Avg. (34 12X0
duser
duaae2
0.77
0.24
057
10
0.35
dUso9
Oaters
0.00
0.34
200
100
1000
200
620
1003
1000
500
10001 500
Figure 2: Bayesian inference with labeled data Here we simulate 10 clusters and assess convergence and mixing of for the DPFA. A. Traceplot of the log-likelihood per spike for all neurons and the first 3 clusters. B. The traceplot of (Frobenius) norms of linear dynamics of (2(t) and X for each cluster (showing the first 3). C. The true and the fitted mean firing rate, showing the averages over samples from iteration 5000 to 10000. D. The true (black) and the fitted (colored) population baseline and latent factor. The X.: denotes the I-th latent factor (i.e. the (-th column of X). The dashed lines show the 95% highest posterior density (HPD) interval. The cosine (the "overlap") between true values and posterior means shown besides. E. Comparison of the held-out likelihood per spike when fitting to 1/2 of the data: 1) 10-cluster model where each cluster has p = (true value), and 2) 1-cluster model where single DPFA describes all neurons, with 14 selected by 5-fold speckled cross volidataii Dots denote results from individual short (iterations 50- 60-100), independent chains.
3.2 Clustering
Using the same simulation, we now infer the cluster labels. The latent factor dimension was first optimized with p selected by 5-fold speckled CV on small chains (100 iterations). We then compare three chains fitting with all data: two unique chains initialized using a single cluster and one chain initialized with N 50 clusters (i.e. all clusters are singletons). Traceplots (Fig. 3A) of training log-likelihood and number of clusters show that all chains converge. When fit to the full data or only half, the number of clusters converges to 10, although the prior over the number of clusters is K- ~ Geometric(0.22 When the recording length is sufficiently long, the likelihood will dominate, and the number of clusters will not be much affected by the prior setting. The true mean firing rate for each neuron (Fig. (3B) can be well recovered, even with half data held out. To evaluate cluster membership, here we show a similarity matrix where the entry (i. is the posterior probability that data points i and / belong to the same cluster. The clustering results for all 4 chains recover the true clusters, no matter what the the starting assignment is (Fig. 3C), which suggest the convergence of MCMC. The overall performance of the full model, where cluster membership is inferred alongside the latent states, similar to the case when cluster labels are known and substantially higher than the 1 1-cluster model (Fig. 3D and E).

---

A
full data
1/2 data
true
full data, chain
1/2 data
custart
0.045
chaar sien
custer2
14
4.93
12
11
cluster10
800
600
ECO
1000
C
Full data, chain
full data, chain 2
E
36
10 dusters, tru (p=22
ddantencen (p=2)
dlsster (pm14)
Inldooul logo hotholl spake
-1
Frequenc)
20
full data, chain
1/2 data
-1.05
500
1000
-1.00
Eoration
ing-linlPhosk
Figure 3: Bayesian clustering The same simulation setting in Fig. 2 but inferring cluster labels from spike observations alone. A. The trace plots of training log-likelihood per spike and the number of clusters. The model is fitted by using all and half data as training with three chains shown for full data (initialized with single cluster, Cl or N = 50 clusters, CN). B. The true and fitted mean firing rate, averaging from iteration 500 to 1000. C. The posterior similarity matrix for all chains (rows and columns ordered according to the ground truth). D. Traceplot of the held-out (1/2) log- oigllikelhhood for 1) the 10-cluster model with labels known, 2) the full model estimating labels and optimizing p. and 3) the 1-cluster, single population model with 14 chosen by CV. E. Held-out log- log-likelihoo per spike for each model (samples for iteration 500 to 1000).
4 Multi-region neural spike recordings
We then apply the proposed clustering method to the Allen Institute Visual Coding Neuropixels dataset. The dataset contains spiking activity from hundreds of neurons from multiple brain regions of an awake mouse. See detailed data description in (Siegle et al.) 2021). Here we investigate the clustering structure of neurons from four anatomical sites (83 neurons): 1) hippocampal CAI (24 neurons), 2) dorsal part of the lateral geniculate complex (LGd, 36 neurons). 3) lateral posterior nucleus of the thalamus (LP,, 12 neurons) and 4) primary visual cortex (VISp. neurons). And we analyze responses to 20s epochs during three visual stimuli: drifting gratings, spontaneous activity. and natural movies. Only neurons with rates > 1Hz within the selected epochs are included (72% of 115 neurons) and we analyze data with 40ms bins. We use Geometric(O.33) prior over the number of clusters, such that p(k S 4) 0.8.
In responses to drifting gratings. the 5-fold speckled CV log-likelihood optimized with p = 2, and, as in the simulations, the lokelikeld and number of clusters show rapid convergence and mixing 4A and B). Low firing rates and short recording lengths tend to cause confusions in clustering, reflecting uncertainty
cluster membership for neurons with little information. Here the average number of clusters is 16. To summarize the clustering results stored posterior samples in MCMC, we give the single estimate for cluster indices & by maximizing the posterior expected adjusted Rand index (maxPEAR, Fritsch and Ickstadt 2009). The maxPEAR--ortee neural activity and posterior similarity matrix are shown Fig. 4C and Results sorted by Maximum a posteriori (MAP) estimate are similar and are shown in the Appendix (Fig. 5). To examine the relationship between the clustering results and anatomy, we additionally sort neurons according anatomical labels (upper left panel in Fig. 4E). Although many identified clusters neurons from the same anatomical area,
7

---

clusters also include neurons from different regions and neurons within a region are often clustered into separate populations (P(neuron i, in the same region| 2; 21, {y: ) 0.57). Together, these results suggest that simple assignment of populations based on anatomy many not accurately represent the latent structure.
We then evaluate the clustering patterns for different visual stimuli. We run 2 independent chains for each epoch (results from the second chain in Fig. 2D). The similarity matrices show that the pattern consistent for the same epoch. but will change along the time even under the same experimental settings (D1 vs. D2 and S1 vs. S2). The changes in the clustering patterns may suggest long-term drift for neuron interactions. To quantify the observations, we evaluate the adjusted Rand index (ARI) of maxPEAR estimates (bottom right panel in Fig. 4). Between-epoch comparisons tend to have lower similarity (average ART from comparing 2 chains for each epoch) than within-epoch comparisons (different chains) for both maxPEAR and MAP (Fig. BC).
The MAP number of clusters is largest (18) for the natural movie, suggesting this epoch has the most severe global non-linearity issue. Here we compare three models: 1) clustering model with p 2, 2) single cluster model, with p = 8 selected by 5-fold speckled CV and 3) anatomical cluster model. The anatomical cluster model fits single DPFA for each region, using the anatomical labels to define the clusters explicitly. Using cross-vali,ation, we find that the optimized dimension for each region is p (1, 11, 5, 3), respectively. We find that the single cluster model tends to underfit the data, while the anatomical cluster model tends to overfit the data (Fig.
A
B
E
TORD 2y- BEEM
TAMM
-1.
10-
500
C
TR EAAA
sorted by RASPER
Soor
120
300
530
Figure 4: Application in Neuropixels data. A. and B. The trace plots of logelikelihood per spike and number of clusters for drifting grating responses. All results are averages from iteration 500 to 1000. C. The observed spikes and fitted mean firing rate, sorted by the maxPEAR label. D. The posterior similarity matrix, sorted by the maxPEAR label. E. The posterior similarity matrices for 4 adjacent epochs and I further epoch with different visual stimuli, sorted sorted according to maxPEAR estimate and anatomical label in the first drifting grating epoch. The last panel shows the adjusted Rand index of the maxPEAR estimates. The diagonal the ARI between two chains for the same data, while off- dif-oiagonal values show the mean ARI of maxPEAR for the four comparisons between two chains from two different epochs. F. For the natural epoch, we hold out 1/2 data as the training, and show the histograms of training and held-out log-likelihooe per spike, from iteration 500 to 1000 for three models: 1) clustering model, 2) single cluster model, and 3) anatomical cluster model.
5 Discussion
Here we introduce a Bayesian approach to cluster neural spike trains by MCMC. Previous approaches to multi-population latent variable modeling have used anatomical information to label distinct groups of neurons, but this choice is somewhat arbitrary. Brain region and cell-type, for instance, can give
8

---

contradictory population labels. The proposed method groups neurons by common latent factors, which may be useful for identifying "functional populations" of neurons. Here we use a mixDPFA model and infer the number of clusters by MFM with a partition- algorithm similar to DPM. MFM may be more conceptually appropriate than DPM, since the number of neural "populations" is unknown but finite. A Additionally, MFM produces more concentrated, evenly dispersed clusters (see Mileer and Harrison 201 8 for detailed discussion). The mixture modeling approach may also be appropriate in cases where neurons share non-homogeneous inputs, since it can approximate global nonlinearity with mixture of locally linear models. Here we find that the mixture model outperforms globally linear (1 (1-cluster) models in simulations and with experimental data.
Although the proposed method can describe data and cluster neural spiking activity successfully, there are some potential improvements. Firstly, as mentioned above, the unconstrained model does not have unique solutions. To ensure model identifiable, we put diagonal constraints on A and Q and constrain (.( and ( to have mean zero The assumption that A and QG) are diagonal does not allow interaction between latent factors. However, these interactions could be allowed by instead constraining ( X to be diagonal (Krzanowski and Marriott 1994b a Fokoué and Titterington 2003). Such a constraint could allow unique solutions for the (P)LDS and GPFA. A second potential improvement would be to automatically infer the dimension of the latent factors (states). In this paper, we assume p; is the same for all clusters, for convenience. p is a pre-selected value or can be selected by velidaaaiot (CV). This may limit the accuracy of the model, since populations of neurons in experimental data are likely to have different latent dimensionalities. In future work, would also be possible treat pj as a parameter and sample the posterior by reversible-jump MRJM (Lopes and West 2004). birth-death (BD)MCMC (Stephens) 2000) Fokoué and Titterington 2003). or
adaptive Gibbs sampling with shrinkage prior on X (Bhattacharya and Dunson 2011). Although RJMCMC and BDMCMC can be easily implemented, they perform poorly for high dimensional data and may be sensitive to priors. Adaptive Gibbs sampling with shrinkage, on the other hand, has been implemented with the infinite mixture of infinite factor analyzers (IMIFA, MMurphy et al. (020). The
same idea may be useful here with an additional prior on linear dynamics (AO), 6(3) and Q()) to
encourage shrinkage in X Finally, a deterministic approximation of MCMC, such as variational inference may more computationally efficient. Standard methods for fitting the PLDD could be used directly in the VI updates, and if we further use a stick- tick-breaking representation for the MFM model, it would be straightforward to use VI for clustering as well, similar to (Blei and Jordan 2006).
As the number of neurons and brain regions that neuroscientists are able to record simultancously continues grow, understanding the latent structure of multiple populations will be a major statistical challenge. The Bayesian approach to clustering neural spike trains introduced here converges fast and is insensitive to the initial cluster assignments, and may, thus, be a useful tool for identifying "functional populations" of neurons.
Acknowledgments
This material is based upon work supported by the National Science Foundation under Grant No. 193 1249.

---

References
Arminger, G., Stein, P., and Wittenberg, J. (1999). Mixtures of conditional mean-and covariance-ssrrutuue models. Psychometrika, 64(4):475-49) 4944.
Berman, B., Johnson, W. O., and Shen, W. (2022). Normal Approximation for Bayesian Mixed Effects Binomial Regression Models. Bayesian Analysis, pages 1-21.
Bhattacharya, A. and Dunson, D. B. (2011). Sparse Bayesian infinite factor models. Biometrika, 98(2):291.
Blei, D. M. and Jordan, M. I. (2006). Variational Inference for Dirichlet Process Mixtures. Bayesian Analysis, 1(1):121-144.
Buesing, L., Machado, T. A., Cunningham, J. P., and Paninski, L. (2014). Clustered factor analysis of multineuronal spike data. In Advances in Neural Information Processing Systems, volume 27.
Chan, A. B. and Vasconcelos, N. (2009). Bayesian poisson regression for crowd counting. Proceedings of the IEEE International Conference on Computer Vision, pages 545-551.
El-Sayyad, G. M. (1973). Bayesian and Classical Analysis of Poisson Regression. Journal of the Royal Statistical Sociery: Series B (Methodological), 35(3):445
Emtiyaz Khan, M., Aravkin, A., Friedlander, M., and Seeger, M. (2013). Fast dual variational inference for non-conjugate latent gaussian models. In Proceedings of the 30th International Conference on Machine Learning, volume 28, pages 951-959.
Fokoué, E. and Titterington, D. M. (2003). Mixtures of factor analysers. Bayesian estimation and inference by stochastic simulation. Machine Learning, (011-2)733
Fritsch, A. and Ickstadt, K. (2009). Improved criteria for clustering based on the posterior similarity matrix. Bayesian Analyvis, 4(2):367-391.
Ghahramani, Z. and Hinton, G. E. (1996). The EM Algorithm for Mixtures of Factor Analyzers.
Glaser, J.. Whiteway, M.. Cunningham, P., Paninski. L., and Linderman, S. (2020). Recurrent Switching Dynamical Systems Models for Multiple Interacting Neural Populations. In Advances in Neural Information Processing Systems, volume 33, pages 14867-14878.
Hastings, W. K. (1970). Monte Carlo sampling methods using Markov chains and their applications. Biometrika, 57(1):97-109.
Johnstone, I. M. and Lu, Y. (2009). On Consistency and Sparsity for Principal Components Analysis in High Dimensions. Journal of the American Statistical Association, 104(486):682.
Jordan, M. I., Ghahramani, Z., Jaakkola, T. S., and Saul, K. (1999). An Introduction to Variational Methods for Graphical Models. Machine Learning 1999 37:2, 37(2):183-233.
Jun, J. J., Steinmetz, N. A., Siegle, H., Denman, D. J., Bauza, M., Barbarits, B., Lee, A. K., Anastassiou, C. A., Andrei, A., Aydin, Ç.. et al. (2017). Fully integrated silicon probes for high-density recording of neural activity. Nature, 55 551(7679):232225
Keeley, S.. Zoltowski, D., Yu, Y., Smith, S., and Pillow, J. (2020). Efficient non- non-conjugate Gaussian process factor models for spike count data using polynomial approximations. In Proceedings of the 37th International Conference on Machine Learning, volume 119, pages 5177-5186.
Krzanowski, W. J. and Marriott, F. H. C. (1994a). Multivariate analysis. Part 2: Classification, covariance structures and repeated measurements.
Krzanowski, W. J. and Marriott, F. H. C. (1994b). Multivariate analysis. Part I. Distributions, ordination and inference.
Linderman, S., Adams, R. P., and Pillow, J. W. (2016). Bayesian latent structure discovery from multi-neuron recordings. In Advances in Neural Information Processing Syxtems, volume 29.
Linderman, S.. Johnson, M.. Miller, A., Adams, R., Blei, D., and Paninski, L. (2017). Bayesian Learning and Inference in Recurrent Switching Linear Dynamical Systems. In Proceedings of the 20th International Conference on Artificial Intelligence and Statistics, volume 54, pages 914 9222.
Lopes, H. F., Salazar, E., and Gamerman, D. (2008). Spatial dynamic factor analysis. Bayesian Analysis. 3(4):759-792.
10

---

Lopes. H. F. and West, M. (2004). Bayesian Model Assessment in Factor Analysis. Statistica Sinica, 14:41-67.
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.
Mctropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., and Teller, E. (1953). Equation of state calculations by fast computing machines. The Journal of Chemical Physics, 21(6):1087-1092
Miller, W. and Harrison, M. (2018). Mixture models with a prior on the number of components. Journal of the American Statistical Association, 113(52)):340.
Murphy, K., Viroli, C., and Gormley, (2020). Infinite Mixtures of Infinite Factor Analysers. Bayesian Analysis, 15(3):937-96
Neal, R. M. (2000). Markov Chain Sampling Methods for Dirichlet Process Mixture Models. Journal of Computational and Graphical Statistics, 9(2):299-265.
Pandarinath, C., O'Saa,, D. J., Collins, J., Jozefowicz, R., Stavisky, S. D., Kao, J. C., Trautmann, E. M., Kaufman, T., Ryu, S. I., Hochberg, L. R., et (2018), Inferring single-trial neural population dynamics using sequential auto-encoders. Nature Methods, 15(10)8005 815.
Paninski, L.,, Ahmadian, Y., Ferreira, D. G., Koyama, S., Rahnama Rad, K., Vidne, M., Vogelstein, J., and Wu, W. (2010). new look at state-space models for neural data. Journal of Computational Neuroscience, 29(1-2):107-126.
Peñaa, and Poncela, (2004). Forecasting with nonstationary dynamic factor models. Journal of Econometrics. 119(2):291-321.
Polson, N. G., Scott, J. G., and Windle, J. (2013). Bayesian inference for logistic models using pólya-gamma latent variables. Journal of the American Statistical Association, 108(504):1339-1349)
Semedo, J. D., Zandvakili, A., Machens, C. K., Yu, B. M., and Kohn, A. (2019). Cortical areas interact through communication subspace. Neuron, 102(():249-255))
Siegle, J. H., Jia, X., Durand, S., et al. (2021). Survey of spiking in the mouse visual system reveals functional hierarchy. Nature, 592(7852):86-92.
Stephens, M. (2000). Bayesian analysis of mixture models with an unknown number of components an alternative reversible jump methods. The Annals of Statistics, 28( 28)4440-7
Ventura, V. (2009). Traditional waveform based spike sorting yiclds biased rate code estimates. Proceedings of the National Academy of Sciences of the United States of America, (06(17):6921-692))
Wainwright, M. J. and Jordan, M. I. (2008), Graphical models, exponential familics, and variational inference. Foundations and Trends in Machine Learning, 1(1-2):1-305.
Whiteway, M. R., Socha, K., Bonin, V., and Butts, D. A. (2019). Characterizing the nonlinear structure of shared variability in cortical neuron populations using latent variable models. Neurons, behavior, data analysis and theory, 3(1).
Williams, A., Degleris, A., Wang. Y., and Linderman, S. (2020). Point process models for sequence detection in high-dimensional neural spike trains. In Advances in Neural Information Processing Systems, volume 33, pages 14350-14361.
Windle, J., Carvalho, C. M., Scott, J. G., and Sun, L. (2013). Efficient Data Augmentation in Dynamic Models for Binary and Count Data.
Wu, A., Roy, N. A., Keeley, S., and Pillow, J. W. (2017). Gaussian process based nonlinear latent structure discovery in multivariate spike train data. In Advances in Neural Information Processing Systems, volume 30.
Yu, B. M., Cunningham, J. P., Santhanam, G., Ryu, S. L., Shenoy, K. V., and Sahani, M. (2009). Gaussian- process factor analysis for din-ennenii single-trial analysis of neural population activity. Journal of Neurophysiology, ()2(1))614-633

---

Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's
contributions and scope? [Yes]
(b) Did you describe the limitations of your work? [Yes]
(c) Did you discuss any potential negative societal impacts of your work? [No]
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes]
2. If you are including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [N/A]
(b) Did you include complete proofs of all theoretical results? [N/A]
3. If you experiments..
(a) Did you include the code, data, and instructions needed to reproduce the main experi- mental results (either the supplemental material or as URL)? [Yes]
(b) Did you specify all the training details (e.g., data splits. hyperparameters, how they were chosen)? [Yes]
(c) Did you report error bars (e.g.. with respect the random seed after running experi- ments multiple times)? [Yes]
(d) Did you include the total amount of compute and the type of resources used (c.g., type of GPUs, internal cluster, or cloud provider)? [Yes]
4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...
(a) If your work uses existing assets, did you cite the creators? [Yes]
(b) Did you mention the license of the assets? [Yes]
(c) Did you include any new assets either the supplemental material or as URL? [No]
(d) Did you discuss whether and how consent was obtained from people whose data you're using/curating? [Yes]
(e) Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content? [N/A]
5. If you used crowdsourcing or conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, applicable? [N/A]
(b) Did you describe any potential participant risks, with links Institutional Review Board (IRB) approvals, if applicable? [N/A]
(c) Did you include the estimated hourly wage paid to participants and the total amount spent on participant compensation? [N/A]
12


Reviews:
[{"review_id": "oyBrKEbhwIE", "reviewer": "Reviewer_QfU5", "summary": "The authors propose a model to identify distinct functional clusters of neurons in neural population activity without an an a priori determined number of clusters. Each cluster is characterized by a latent LDS and mapped through an exponential nonlinearity to poisson rates. The prior over number of clusters and cluster assignment is computed efficiently borrowing from recent work on a mixture of finite mixture (MFM) model. The authors use a bespoke procedure for using MCMC to evaluate posteriors for their model, which involves a laplace approximation to estimate the population states, and a poisson-gamma conjugacy to marginalize over of the latent loadings $c_i$ to evaluate the probability of a neuron belonging to a given cluster. The authors show results on simulated data where the true number of clusters is known, and then evaluate the model on neuropixel recordings. They show on the neuropixel data that the clusters identified using their method outperform anatomical clustering or a single cluster on held-out data. This work is an important step in identifying functional groups of neurons in large multi-region neural population recordings.", "questions": " Some more detailed clarification on the inference would be helpful. Could you expand on how the laplace approximation is being used in 122-124? More organization in the introductory part of the inference section would I think go a long way. An equation showing the P(x|..) being approximated by Laplace would be helpful. Also it appears there is a typo in the third line of the laplace approximation in the appendix (floating closed parentheses).\n\nCan you also clarify the difference between how the marginal likelihood in (3), evaluated using a conjugacy using Poisson and inverse gamma distributions, relates to the nonconjugacy in the \"update $z_i$\" section of the appendix. Some reference to equation (3) from the main paper in this section of the appendix would also be helpful for clarification.\n\nA few more sentences providing background on the MFM approach might also help clarify this section. What exactly is the form of H?  \n\nI think giving the specific model here a name would be useful. This would help the authors more easily reference their specific approach in the paper and later allow future researchers to succinctly refer to the work. \n\nThe resolution of the figures is low. Zooming in to see detail is difficult, particularly in 2D. ", "limitations": "One limitation that I believe is more significant than the author's state is the selection of p being constant across subpopulations. The authors briefly mention ways around this issue in the discussion. I think it is important to more directly point this out as a limitation -- groups of neurons across brain regions very clearly do not have the same dimensionality (low level visual areas and higher level decision making areas, e.g.) This model is only able to isolate subpopulations with this specific constrained feature. \n\nNo potential negative social impact of this work.", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 4, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "This paper is technically strong and is an important tool with potential broad use in neuroscience. The addition of a prior over cluster assignment was a non-trivial ML contribution and significant step in understanding increasingly common multi-region datasets with large numbers of neurons. However, there are many pieces to the inference procedure and the appendix has a lot going on. The authors reference many existing approaches and it can be hard to parse the specific approach used in this work.\n\n", "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "bLslKcgePFE", "reviewer": "Reviewer_dDsp", "summary": "This manuscript proposes a new approach called \"Dynamic Poisson Factor Analyzers\" (DPFA) that models neural spike trains.  This methodology is developed to split individual neurals into distinct clusters, based upon the hypothesis that neural spike trains from each neural cluster into distinct populations.  This approach uses fast approximations and modern MCMC approaches to implement an inference scheme, which is applied to synthetic data and high-dimensional real datasets.  Results are compared to essentially a PLDS model (no clustering) and a clustering model based upon anatomical location.", "questions": "What scientific question will this methodology help answer?\n\nIn what real-world datasets could this method significantly outperform the baselines?\n\nAssuming discrete clustering of neurons is a strong assumption.  To my knowledge, most neurons participate in multiple subnetworks of the brain.  Why are the authors considering an explicit clustering of neurons rather than an admixture, which may be more biologically realistic?", "limitations": "No concerns.", "rating": 4, "confidence": 5, "soundness": 4, "presentation": 3, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The biggest strength and contribution of this manuscript is the inference algorithm.  The inference algorithm seems very well-considered and explored, using fast approximations in the form of well-design Laplacian approximations and approximate conjugate distributions.  Combining these tricks with a Mixture of Finite Mixture (MFM) approach, rather than a Dirichlet Process, seems to mix fast and provide robust inference.  This approach seems fairly original and significant, although I would appreciate additional details and discussion on the efficacy and reliability of individual approximations (e.g., how close is the gamma-poisson distribution to the lognormal?).\n\nThe model itself is a reasonable model, but does not seem especially novel or significant given the existence of MixPLDS.  There are changes and advances, but they seem relatively minor.  It would be helpful if the authors could highlight their particular advances more clearly, and explicit discuss how this improves the model.\n\nOne of the areas that would benefit from improvement is clarity on how this relates to a clear scientific challenge.  How the neurons cluster is potentially interesting, but I believe that the manuscript would be greatly improved by focusing on a clear scientific question that this methodology can help answer.\n\nThe results on real data, though, are one of the biggest weaknesses.  In the real data, the proposed method does not fit the data better based on log likelihood than a clustering based on known anatomy, and is barely better than not using a clustering approach at all.  It is difficult to motivate using this approach based on those results for me, especially as other competing approaches such as deep RNNs (e.g., LFADS) do tend to fit the data better.\n", "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "LubRU5yvC5", "reviewer": "Reviewer_WZGP", "summary": "The authors develop a novel model based on a mixture of dynamic Poisson factor analyzers to cluster neurons in large-scale neural recordings. They treat the number of clusters treated as an unknown parameter and propose a novel MCMC algorithm to efficiently sample its posterior distribution. They show in a simulation that their model can accurately recover the true clustering and latent states, and then explore the results of their method on a recording from the Allen Institute Visual Coding Neuropixels dataset.", "questions": "The model is very interesting and well constructed. However, it is a bit hard to grasp, and I feel like the section 2. could be a little more clear. I find that the graph in Fig. 1, C is not super informative, and could include more information, such as the priors and the additional parameters. H is also never introduced before it is mentioned on line 105. I understand that it corresponds to h, but in this case shouldn't there be an additional line in Equation 1 for mu and H? \n\nI would have also liked to gain more intuition on why you chose the different priors. For example, I understand the use of a Geometric distribution as the prior on k, and that the mean of this prior should increase when we want to get more clusters, but an experiment showing the different results and recommended parameter with different datasets, probes, tasks, or number of brain region recording would help a lot. \n\nYou propose to use cross-validation for choosing p. I assume you optimize over the log-likelihood. Won't the log-likelihood necessarily increase as you use higher values for p? What criteria is used? \n\nThe latent states are chosen to evolve linearly over time with Gaussian noise. This makes sense when the animal is not doing any tasks. What happens when the animal change state, receives a stimulus or start doing a task? It would be great to add a visualization of the latent states throughout trials. \n\nIn section 2.3, you mention that \"Though we may evaluate it by a Laplace approximation, but iterating over all potential clusters for each neuron is computationally intensive. To make faster clustering, we approximate the marginal likelihood by utilizing a Poisson-Gamma conjugacy.\" and \"Another possible idea is to approximate the log-likelihood by second-order polynomials, with coefficients determined by Chebyshev polynomial approximation [Keeley et al., 2019]. However, we find that this approximation doesn’t work well in practice when spike counts have a wide range.\". Could you show results and comparisons, in both speed and accuracy for these three approaches? \n\nThe main flaw of the model, to me, is that \"Low firing rates tend to cause confusions in clustering\". Low-firing rate neurons can encode crucial information in the brain, and artificially removing neurons with rates < 1Hz in the experiments does not seem to have any biological justification, but simply a way to get nicer results ass the firing rates will be more uniform. Some neurons can also one strongly associated to a state of the animal and have overall low firing rate as they fire strongly during specific tasks. \nCould this be due to the choice of prior on lambda (and more specifically on delta/mu?)? The prior on delta is not specified (unless I miss it), but this could maybe help and solve this issue? \nMoreover, the results with low-firing rate units are not even included in the paper. It is essential to include them (at least in the appendix) to have a sense of how it contaminates the clustering results. \nWhat proportion of the neurons are discarded when keeping only neurons with FR>1Hz?\n\nThe low firing rate neurons can also be due to artifacts of spike sorting. It is especially true in Neuropixels datasets, where neurons can be split and thus have their firing rate cut in half (or more). There are many quality metrics developed for the output of spike sorting. (Hill et al. (2011), Harris et al. (2001), Schmitzer-Torbert and Redish. J Neurophy (2004), Chung et al. (2017)...). I think using such metrics to discard \"bad\" units would be much more appropriate (as is done in Reproducibility of in-vivo electrophysiological measurements in mice, International Brain Laboratory (2022) for example).\nit would also be interesting to look at how the clusters are impacted by spike sorting artifacts. If you split a neuron in half, are the two half clustered together? How are the clusters impacted by electrode drift? You need to make sure the clustering algorithm is robust to these sources of noise.  \n\nWhat exactly is the \"anatomical cluster model\" mentioned in the result section. Are you running your clustering model on each region independently, using different values for p? (p=[1,11,5,3]). How do you choose the number of clusters per region? Do you use Cross-validation on each region? If you cluster neurons from each region independently, are the neurons clustered together also clustered together when clustering all neurons together? This would be a very interesting experiment to show in the paper.\n\nFinally, MCMC models can be computationally expensive, and your model requires (5-fold) cross-validation. How does it scale with the number of neurons and length of the recording? ", "limitations": "I don't see any potential negative societal impact of this work.", "rating": 6, "confidence": 5, "soundness": 4, "presentation": 3, "contribution": 4, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The strength of the paper lies in the proposed model, that is novel and very well-thought. It is technically sound and well constructed. However, it is a complicated model and the section 2 of the paper could be a little more clear and give more intuition behind the choice of priors and parameters. The results on simulated data are very compelling, but I would have liked to see more experiments on real data (for example, how to choose the parameter of the geometric prior on k depending on the type of probe, or the number of anatomical regions) to gain more intuition. The main weakness of the model lies in the fact that \"Low firing rates tend to cause confusions in clustering\" and are excluded from the data in the analysis. This does not seem very biologically relevant.  ", "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}]

Author rebuttal:

Decision: Accept