Bayesian Bi-clustering of Neural Spiking Activity with Latent Structures
Authors: Ganchao Wei
Venue: ICLR 2024 poster (2024)
Keywords: neural spikes
spatio-temporal statistics
Bayesian nonparametric
bi-clustering
Primary area: applications to neuroscience & cognitive science

Abstract: Modern neural recording techniques allow neuroscientists to obtain spiking activity of multiple neurons from different brain regions over long time periods, which requires new statistical methods to be developed for understanding structure of the large-scale data. In this paper, we develop a bi-clustering method to cluster the neural spiking activity spatially and temporally, according to their low-dimensional latent structures. The spatial (neuron) clusters are defined by the latent trajectories within each neural population, while the temporal (state) clusters are defined by (populationally) synchronous local linear dynamics shared with different periods. To flexibly extract the bi-clustering structure, we build the model non-parametrically, and develop an efficient Markov chain Monte Carlo (MCMC) algorithm to sample the posterior distributions of model parameters. Validating our proposed MCMC algorithm through simulations, we find the method can recover unknown parameters and true bi-clustering structures successfully. We then apply the proposed bi-clustering method to multi-regional neural recordings under different experiment settings, where we find that simultaneously considering latent trajectories and spatial-temporal clustering structures can provide us with a more accurate and interpretable result. Overall, the proposed method provides scientific insights for large-scale (counting) time series with elongated recording periods, and it can potentially have application beyond neuroscience.

Published as a conference paper at ICLR 2024
BAYESIAN BI-CLUSTERING OF NEURAL SPIKING AC- TIVITY WITH LATENT STRUCTURES
Ganchao Wei
Department of Statistical Science
Duke University
Durham, NC 27708, USA
ganchao, wei wei@duke, edu
ABSTRACT
Modern neural recording techniques allow neuroscientists to obtain spiking ac- tivity of multiple neurons from different brain regions over long time periods, which requires new statistical methods to be developed for understanding struc- ture the large-scale data. In this paper, we develop a bi-clustering method to cluster the neural spiking activity spatially and temporally, according to their low- dimensional latent structures. The spatial (neuron) clusters are defined by the la- tent trajectories within each neural population, while the temporal (state) clusters are defined by (populationally synchronous local linear dynamics shared with dif- ferent periods. To flexibly extract the bi-clustering structure, build the model non parametrically, and develop an efficient Markov chain Monte Carlo (MCMC) algorithm to sample the posterior distributions of model parameters. Validating our proposed MCMC algorithm through simulations, we find the method can re- cover unknown parameters and true bi-clustering structures successfully. We then apply the proposed bi-clustering method to multi-regional neural recordings un- der different experiment settings, where we find that simultancously considering latent trajectories and spatial-temporal clustering structures can provide with a more accurate and interpretable result. Overall, the proposed method provides scientific insights for large-scale scale- (counting) time series with elongated recording periods, and can potentially have application beyond neuroscience.
INTRODUCTION
In neuroscience, identifying types of neurons a longstanding challenge (Nelson al., 2006: Bota & Swanson 2007 Zeng 2022). Some criteria based on features such as anatomical regions, ge- nomics and synaptic connectivity have been proposed, and there are some Bayesian approaches to integrate these features (Jonas & Kording 2015). On the other hand, the response pattern and inter- actions between neural populations may change over time, especially when the experiment stimuli changes (Pooresmaeili et al., 2014, Oemisch et al. 2015 Ruff & Cohen 2016 Steinmetz et al.) 2019 Cowley et al., (020). However, these complex dynamical observations can often be broken down into simpler units, and it can be appropriate to assume static linear dynamics within chunks of the data. Moreover, it is usually appropriate and helpful to assume similar linear dynamics can be shared by different epochs (Zoltowski et al. 2020 Glaser al. 20206). Here, we consider the problem of how to identify both spatial and temporal clusters of neural spikes.
The modern techniques such as the high-density probes (Jun allaa2211111 2017 Steinmetz et al., 2021| Mar- shall et al. (022) allow us to obtain large-scale multi-clectrode recordings from muluple neurons across different anatomical regions over an elongated session. Several models have been developed to extract shared latent structures from simultaneous neural recordings, assuming that the activity of all recorded neurons can be described through common low-dimensional latent states (Cunningham & Yu 2014 Gao al. 2017). These approaches have proven useful in summarizing and interpret- ing high-dimennsioal population activity. Inferred low diwdmmmnniioal latent states can provide insight into the representation of task variables (Churchland et al.] 2012 Mante et al. [013 Cunningham &
Yu 2014 Saxena & Cunningham (019) and dynamics of the population itself (Vyas et al. 2020).
Many existing approaches are based on linear dynamical system (LDS) model (Macke al.
2011

---

Published as a conference paper at ICLR 2024
which is built on the state-space model and assumes latent factors evolve with static linear dynamies. Although assuming static linear dynamics over time can be valid in some tasks and in small chunks of experiment, the assumption is not generally appropriate. To tackle the nonlinear dynamics, some (SL.DS,Ghahrrmani ee al.] (2008a): Murphy (2012)) and recurrent-SI SLDS (RSLDS, Linderman al. (2017 2019) have been proposed. The parametric Gaussian process factor analysis (GPFA) model (Yu et al. 2009)
variants of the LDS,, such as switching-LDDS (SLDD, Ghahramani & Hinton (2000); Fox (2009); Fox
and its variants provide a more flexible way model non-linaer neural data, although most these methods assume independent GP and doesn't 't allow for interactions between latent factors. Recently, (Cai et al. 2023) proposed the dependent GP method using the kernel convolution framework (KCF, Boyle & Frean (2004); Alvarez & Lawrence 2011); Sofro et al. (0017)) but their method may not scalable for elongated neural recordings. Several methods have been developed and implemented to analyze multiple neural populations and their interactions (Semedo et al. 2019, Glaser et al. 2020b), as the interactions may occur in oow-dimensional subspace (Stavisky et al. 2017] Kaufman et al., 2014 Semedo et al., 2019). But the neural populations are prespecified, and the spatio-temporal clustering structures of neural data haven't I been evaluated systematically in general.
The clustering of neural spikes is an important and long-lasting problem, and people put lot effort into developing methods to uncover patterns from the complex data. The neural spiking activity is essentially a point process, and there are some methods for finding clusters in point process by, such as, Dirichlet mixture of Hawkes process (Xu & Zha 2017). mixture of multi- point process and
(Yin et al., 2021) and group network Hawkes process (Guanhua Fang & Guan 2023). All these methods are general and have wide rage of applications, but the defined "clusters" may not directly related to underlying mechanisms and may lose some scientific insights. To deal with this problem. some methods like miXPLDS (Buesing al., 2014)) and recent mixDPFA method (Wei al. 2022 2023) try to cluster neurons according to their latent structures, by using the mixture of LDS model. These approach provides a more interpretable and accurate way clusters neurons, and may be useful for identifying "functional populations" of neurons. However, these methods assume the static linear dynamies and don't it allow for the interactions between neural populations, which can limit the usage of these methods, and may bias or even fail the detection of neural populations when considering the elongated recordings, especially under different experiment conditions. On the other hand, for the clustering structures in terms of temporal states, most methods are developed based on the SLDS, by modeling the nonlinear dynamics with local linear pattern. Instead of clustering based on linear dynamics. [D' Angelo al.) (2023) recently tried to cluster the experiment periods based on the distributions of spiking amplitude, using a nested formulation of mixture of finite mixtures model (MFMM), i.e., exploiting the generalized MFMM (gFMFMM, Frühhwith-Schhntter et al. 202T prior with common atom model (Denti et al. 2023)
In this research, we develop a bi-clustering method to cluster neural spikes both spatially (to give subject clusters) and temporally (to give state clusters), according to the latent structures of these neurons (TigureTTA). The neural population is defined via private low-dimensional latent trajectories as in miXPLDS (Buesing et al., 2014) or mixDPFA (Wei et al., 2022; 2023). For the state clusters, we assume the linear dynamics can be shared across different chunks and the state clustering struc- tures are defined by local linear manner as in (R)SLDS. Neurons in each population is assume to have private latent trajectories, but all time series are assumed to switch between different states syn- chronously, to use the information from all observations. Besides extending the previous clustering method like mixDPFA to simultancously detect the state cluster, the proposed bi-clustering method also allow for interactions between neural populations and non-stationary dynamics for neural re- sponse, using similar idea from (Glaser et al. 2020a). Simultaneously considering all these effects in the proposed bi-clustering method is necessary, since incorrect population assignments can lead to biased and inconsistent inference on the latent structure (Ventura 2009). On the other hand, these flexibility allows for more accurate estimate of latent trajectories, and hence will lead to a more accurate estimates of the subject clustering structure.
To flexibly infer the bi-clustering structure, we model them non-parametrically to avoid prespecify- ing the number for subject and state clusters. Specifically, the subject clustering structure is modeled by a mixture of finite mixtures model ( (MFMM, Miller Harrison| 2018)) of latent trajectories and the state clustering structure is modeled by a sticky Hierarchical Dirichlet Process Hidden Markov Model (sticky- HDP-HMM, et al. (2008b)) The posteriors of model parameters are sampled us- ing a Markov Chain Monte Carlo (MCMC) algorithm, where the Pólya-Gamma data augmentation technique (Polson et al. 2013) is used to handle the counting observations for neural spikes.
2

---

Published as a conference paper at ICLR 2024
The rest of this paper is structured follows. In Section introduce the bi-clustering method for time series with counting observations, and provide brief explanations of the MCMC algorithm to sample posterior distributions of parameters. After validating the proposed bi-clustering method with a synthetic dataset Section 3 we then apply our method to analyze multi-regional experi- mental recordings from behaving mouse under different experiment settings in Section Finally, in Section 5 we conclude with some final remarks and highlight some potential extensions of our current model for future research.
A
B
clustered
latent
counts
100
200 100
Figure I: Model overview. A. The goal for our proposed model to do clustering both spatially and temporally (i.e. "bi-clustering") for neural spike data (time series data with counting observations), according to their latent structures. The neural spiking counts are determined by low dimensional latent factors, specific to the spatially subject clustering assignment (e.g. green ,blue and red). On the other hand, all neurons are assumed to switch between states synchronously. and are temporally clustered according different states of linear dynamics (e.g. gray and white). B. Graphical model of the proposed bi-clustering model. All prior parameters are summarized as 0, and parameters such d, and e; dropped for simplicity.
BI-CLUSTERING METHOD FOR NEURAL SPIKES
In this section, we introduce a bi-clustering model for neural spiking activity, i.c., the time series data with counting observations. The goal for the proposed model is to cluster neural spikes both spatially (for subject cluster) and temporally (for state cluster), based on the moltiipoppulation and -state latent structures. To flexibly capture the clustering structures, build model non-parametrically, The graphical representation of the model summarized in Figure B. After introducing the model, we briefly describe how a MCMC algorithm is used infer model parameters.
2.1 MULTI-POPULATION AND -STATE LINEAR DYNAMIC MODEL
Assume we can observe spiking activity of N neurons up to recording length T. Denote the number of counts for neuron .,, N) indicator of neuron i as (i.e. the "subject indicator"). Assume the count yı follows a negative- binomial distribution, where the log-mean response modeled by linear combination of subject
at time bin E (1....,77 T}
as Vit E Zzo- and the cluster
(a:) baseline d, population baseline (2:) and p-dimensional latent factor zf(t) E RP (Here we assume all populations have the same latent dimension for convenience). other words, the observation equation is as follows:
Vi NB Mss)
log u4: d; 14 + cxx( (21)
where NB(r, 11) denotes the negative-binomial distribution (NB) with mean and variance μ + 42/r f2/r. and e, ~ (0,, I,). The NB distribution can be replaced by Poisson distribution when it's appropriate assume equi-dispersion, for ease of model inference. If further denote a(()) (1(5) ()) xt fyy and é; = (1, <)y, then log Aix = d, + (x.) specified, refer a((6) "latent factor" for cluster j. which also includes the population baseline.
(&.) To save words and notations, not
3

---

Published as a conference paper at ICLR 2024
Although each neural population is modeled with private latent factors, there usually exist some interactions between clusters (Musall et al. 2019; Stringer et al. 2019), and these interactions can change over time (Ruff& & Cohen|[2111 Steinmetz et al. Sel19Cooleyy al. 2020), especially when the external condition changes. On the other hand, interactions between neural populations and receiving common inputs for all neurons suggest that neurons in different clusters may synchronize the response states over time. Therefore, to allow for the interactions between populations and model the synchronous state switching, we stack the latent factors for all clusters together, and assume all latent factors evolve in conditional linear manner, given the discrete latent states & shared across the cluster, as in Glaser al. (2020b). In other words, the state clustering structure defined by the local linear dynamics of latent factors, by assuming complex dynamies can be decomposed into simple linear unit and the small chunks of the neural response can be sufficiently described by the LDS model. Specifically, assume there are k unique clusters, i.e., (
k. the cluster-stacked latent factors (including population baseline) is denoted as X. = (क) (a))) ................................ a E RR(x++1, To ................................
capture temporal dynamics of the data, we further put AR(1) structure onto the latent factors X,- In other words, given the discrete latent state & as & (i.e. the "state indicator" shared across the subject cluster), X, is assumed evolve linearly with Gaussian noise follows:
X = bb. + Ag.. + cG.
(2)
where Eg. NN Qx.) and X1 NNO,, 1))+))) Here, the dynamic parameters (bg., Ac. Qc,) summarize dynamics (state changes) within and across population. To make the model identifiable, further assume 1) a(5) is zero-centered, i.e. Σr a.(1))
and 2) a (x 111-1 is diagonal, where (03) (z( ..... 2/ R(+11XT Fokoué & Titterington, 2003). The identifiability i:
a/())
issue under pernutttermttaiionnn is handled by alignment to samples early stage of the chain. For more details on model constraints, see Section A.1 In summary, given the neuron belonging to cluster 2; = j. the counting series is generated by negative- binomial linear model M defined in equation as (J.1....,
(()) ~ MMde, C,, x (I:T) ) where the prior for 1:T is denoted as H. The within- and netween-population linear dynamics t-th step are captured by dynamical parameters (b&,, Aç, ). where & is the state indicator at t. To do clustering both spatially (subject cluster)
and temporally (state cluster) in a flexible way, we model each of these two clustering structures non-parametrically as follows.
2.2 SUBJECT CLUSTERING MODEL
Since the number of neural populations is finite but unknown, we put prior on number of subject cluster [ = & as in Wei et al. 2023), which leads to the mixture of the finite mixtures model (MFMM) as follows:
K~ ~ Sk, π " m&) ~ Dirx(r,.....
fk is p.m.f. (1, given K = k, given n
i.d.
Z1.... nN
π
(k) iid.
x
H.
given k
(Jil,...,.r)) .TT) M(di, C, a(x)) given di, Ci x 1 Zi, for I,...,N
where f. denotes the probability mass function. By using the MFMM, we can integrate the field knowledge about the number of clusters into our analysis, by specifying the ft- In the analysis of this paper, we assume & follows geometric distribution, i.e., le Geometric(C) with p.m.f. defied as (k () (1 -)) for k = 1,, 2, .... and
For general use of the proposed method to some problems where the number of subject cluster number can potentially grow to infinity, using the mixture model such the Dirichlet process mixtures model (DPMM) maybe conceptually more appropriate. See Miller Harrison| (2018) for more detailed discussion.
2.3 STATE CLUSTEERIGG MODEL.
For state clustering structure, as the number of states can potentially shoot to infinity, we model the discrete state & by a sticky Hierarchical Dirichlet Process Hidden Markov Model (sticky-HDP-
4

---

Published as a conference paper at ICLR 2024
HMM) proposed by Fox et al. 2008b to encourage the re-occurrence of states as follows:
(T),
vi fid DP(a - m, aB môi
a + m
&1 ~ v2-11
(4)
(b,, A,, Qi) X (be. + As, X. Qei)
isd.
S,
for I = 1,2,-.-7
for 1 = ,
where GEM denotes the stick breaking process (Sethuraman 1994). 6, denotes the indicator func- tion at index i, DP denotes the Dirichlet process and S denotes the normal-inverse-Wisse Wishart prior of (b,, A,, Q1)- - The details of S can be found in the appendix A. when introducing the MCMC algorithm. The sticky HDP-HMM extends the HDP-HMM with "sticky" parameter m 0 to encourage longer state duration, and hence can handle the rapid-switching problem to some degree. Some more careful methods for modeling the state duration and state transition further discussed in the Section 5
2.4 MODEL INFERENCE
We do Bayesian inference on the proposed bi-clustering model by an efficient MCMC algorithm. In each sampling iteration, there are approximately four steps: 1) sample dynamical latent factors (2:) #1:T" 2) sample remaining subject-specific parameters in observation equation equation includ- ing subject baseline d,, factor loading é and dispersion r,. 3) sample the temporal states & and corresponding dynamical parameters (b,, A,, Q:) for each sampled state, and 4) sample the subject cluster indices žj. The details of sampling procedures can be found in the appendix Section|A.. and we briefly introduce the key sampling method for each step here.
In step 1), the full conditional distribution of latent factors =0] is equivalent the posterior distri- binomial dynamic GLM (NB-DGLM), which has no closed form. However, the NB distribution falls within Pólya-Gamma (PG) augmentation scheme (Polson et al. 2013
bution of the negative-binomial
Windle et al.) 2013 Linderman et al. 2016). therefore we can sample them in closed form by in- troducing the PG augmented variables. Conditioning on the auxiliary variables wit the transformed (()
"effective" observations Uie has Gaussian likelihood, and hence we can sample the posterior of a L.T using the forward-fitering--aakwww wew-bbakward--mmlpinn (FFBS. Carter & Kohn (194)); Frrühwirth-Schnatee (994g) algorithm. For Poisson observation model, we can treat the data as coming from the NB distribution, use the samples as proposal and add one more Eatropolis-Hastting (MH) step to accept or reject the proposal. In Poisson case, the dispersion becomes the tuning parameter to achieve desirable acceptance rate (Wei al., 2022 2023).
In step 2), the sampling of d, and é is regular NB regression problem, and we again use the PG data augmentation technique to sample them. The dispersion parameter r, is updated via a Gibbs sampler, using the method described in ZZou et (2012), as the gamma distribution is the conjugate prior to the T3 under the compound Poisson representation.
In step 3), the discrete states & are sampled by weak-limit Gibbs sampler for sticky HDP-HMM as Fox et al. (2008b). The weak-limit sampler constructs finite approximation to the HDP transitions prior with finite Dirichlet distributions, as the infinite limit converges in distribution to true HDP. Given the latent factors a #1.7 and state indicator & we can update dynamical parameters (br, Ai, Qt) for each state separately in closed form.
In step 4), given the parameters in observation equation equation we sample the subject cluster indices using the algorithm for MFMM proposed by Miller & Harrison (2018), which is analogous to the partition-basee algorithm for DPMM (Neal 2000). The label switching issue of z; is handled by the Equivalence Classes Representatives (ECR) algorithm (Papastamoulis & Iliopoulos) 2010), by using the sample from early stage of the chain as pivot allocation. When sampling the subject clustering assignments 2, in such high dimensional time series data with large T, evaluating the full likelihood given samples of c; as in Gaussian MFA (Fokoué & Titterington) 2003) leads to a poor mixing chain. Instead, we evaluate the marginalized likelihood by integrating out the subject-
5

---

Published as a conference paper at ICLR 2024
specific loading Ct, similar to Wei al. (023). The marginalized likelihood is evaluated by Laplace approximation.
The Python implementation of the NB and Poisson bi-clustering model is available in https: 7/github. com/ weigcdsb/bi httpaaabbbiiuuusteeeeeeeeeeeeeee and supplementary material. The additional details for MCMC sampling can be found in appendix Section A.1
3 SIMULATION
To validate and illustrate the proposed bi-clustering method, we simulate neural spikes from the NB bi-clustering generative model defined in observation equation and system equation 2 In this simulation, we generate clusters with 10 neurons in each cluster (NV = 30 total). The recording (5) length is T 500 and the dimension for 21-r are all p p= 2. For each neuron, the individual baseline is generated by d; ~ (0,, 0.5²), the factor loading is generated by G₂ ~ NN,, I2) and dispersion are all r; 10. For latent factors of these three clusters { {= they are generated from two discrete states, and the state indicator & 1, 2 is generated from a semi- Markov chain (Sansom & Thomson 200T] Yu 2010), to encourage longer state duration. These states correspond two sets of linear dynamics: I) independent state, where A E R is diagonal and 2) interactive state, where A a random rotation of an orthogonal matrix, and hence there are interactions between clusters. The bias term is b = 0 and noise covariance is Q Ig 10-2 for both states.
We then apply the proposed bi-clustering methods to the simulated data, by setting the maximum number of states be 10 for the weak-limit sampler. Here, we run MCMC chain for 10,000 itera- tions, starting from subject cluster and 10 uniformly distributed temporal states. The results shown here summarize the posterior samples from iteration 2500 to 10,000. The trace plots for several pa- rameters are shown in appendix (Figure 4. and they don't suggest significant convergence issues. First, to evaluate the inferred clustering structure, (Figure A and subject cluster (Figure 2 The entry (4,j) )) for a similarity matrix is the posterior probability that data points i and j belong to the same cluster. Both subject and state similarity matrices are sorted according to true clustering assignments, and hence if the algorithm can recover the simulated cluster structures, the diagonal blocks will have high posterior probability, which is the pattern shown in Figure Z and 2 The histograms of posterior samples (Figure [C))how that our method can successfully recover the number of subject and state cluster. To represent the clus- tering structure in temporal state more intuitively, we also provide the single point estimates of & (Figure 2) by maximizing the posterior expected adjusted Rand index (maxEPAR, Fritsch & Ick- stadi (2009), which performs better than other points estimates such as the MAP estimates. All these results show that we can successfully recover the clustering structures on spatial and temporal dimension, including the number of clusters for each.
check the similarity matrices for both state
One the other hand, the advantage for proposed bi-clustering method is that it can simultaneously provide unbiased estimates of the latent trajectories for each cluster, which can be very helpful for scientific interpretation and insights. Here, we show the latent trajectories for the cluster that most neuron 1-10 belong to in Figure 2 (subject 1-10 also forms the maxPEAR subject cluster). The samples for a! (excluding n&)) in Figure 2 are rotated by single matrix, to match the posterior mean to ground truth (the rotation matrix is found least square between posterior mean and ground truth). The trace plots for L2/ Frobenius norms of latent trajectories (show both raw and rotated traces for x!) for each detected cluster is shown in Figure 4C-D. The fitting results show that si- multancously considering subject and state clustering structure is necessary for estimation of latent structures. We also show the performance of Poisson bi-clustering model (i.e. replace the NB distri- bution in observation equation by Poisson distribution). The clustering results are summarized by similarity matrices (Figure 2 and maxPEAR estimate of state (the third bar in Figure 2D). Since in this simulation example, the over-dispersion is not severe (r; = 10), the Poisson version can also recover the true clustering structure, but with some bias in the estimation of latent trajectories. How- ever, for data with large over-dispersion (which is common for real data), the wrong assumption on equi-dispersiie will hugely influence the clustering structures, and it would be necessary to use the more flexible NB bi-clustering model.

---

Published as a conference paper at ICLR 2024
A
B
state cluster atateaarrdd
subject cluster
maxPEAR states
Poisson model
cluster
123
NB
Prissom
390
361
latent trajectories of cluster nt 11-110
postertor
subject cluster
Xm
2000
stalt clusler
e sutjjct chasler
Figure 2: Simulations. Here, we show the results for posterior samples from iteration 2500 to 10,000 for each MCMC chain on the simulated dataset. These results are from NB bi-clustering model not specified. A. The posterior similarity matrix for temporal states are ordered according to ground true states, representing the inferred clustering structures relative to ground truth. B. Spa- tially, the similarity matrix for subject are ordered according true subject clusters. C. The histograms of posterior samples number of state cluster (true 2) and subject cluster (true 3). D. The max- PEAR estimates (point estimates) the discrete states for NB and Poisson bi-clustering model, comparing to the true temporal states. E. The inferred latent trajectories for the detected cluster that most subject 1-10 belong to. Here μr denotes /1::: X denotes i-th row of ø1.7, and samples of (a) transformed by single matrix, to match posterior mean and ground truth (by least square). The black lines are truths, blue lines are posterior means and shaded light blue regions are 95% highest posterior density (HPD) regions. F. The similarity matrices of state and subject cluster for Poisson bi-clustering model, sorted using the same order as in panel A and respectively.
APPLICATION
We then apply our bi-clustering method to Allen Institute Visual Coding Neuropixels dataset. The dataset contains neural spiking activity from multiple brain regions of an awake mouse, under dif- ferent visual stimuli. See SSiegle et al. (202T) for more detailed data description. Here, we use electrophysiology session 719161330 to investigate the bi-clustering structures of neurons from three anatomical sites, under three consccutive experimental epochs. After excluding neurons having less than IHz response rate, 78 neurons are contained in the following analysis. Among these neu- rons, 37 originate from the hippocampal CA1, 20 from the lateral posterior nucleus of the thalamus (L.),, and 21 from the primary visual cortex (VISp). The neural spikes are recorded when the mouse exposed three consecutive visual stimuli spontaneous (S, durates 30.025s), natural movie (N, durates 300.251s) and again spontancous (S, durates 30.025s). These three epochs are chosen, according to the hypothesis that spontaneous connectivity patterns in visual cortex are shaped by ha- bitual coactivations produced by natural visual stimulation (Fiser et al., 2010; Harmelech & Malach, 2013 Sadaghiani & Kleinschmidd, 2013), and the similarity in connectivity patterns between spon- taneous and natural movie are found in human visual cortex (Wilf et al.) 2017). Here, we rebin the data with 500ms, and hence T 720. For formal application, we may need smaller bin size for the higher resolution. The binned spiking counts for these 78 neurons are shown in Figure AA.
Then, we the data with both NB bi-clustering and Poisson bi-clustering model, and run two independent chains for each. The results from all four chains can be found in the Section A.3 Figure 6 Although the formal analysis requires us to tune some parameters such as latent dimension p and the sticky parameter in equation 4 here
10,000 iterations using p 2 and simply to illustrate the usage of proposed method on real data. Since these neurons come from three brain regions, set the prior for number of subject cluster Geometric(O.415) such that
10,

---

Published as a conference paper at ICLR 2024
P(k- 3) 0.8. The trace plots for several parameters provided in appendix (Figure 5 which don't 't suggest significant convergence issues.
For subject clustering structure, the NB bi-clustering model detects around 10 clusters (histogram and trace in Figure(5A-B), and the posterior similarity matrix sorted by maxPEAR estimate is shown in Figure CC- Generalluy, the method detects large neural population with high "confidence", with several weak clusters. We further sort the similarity matrix according to the anatomical labels, to examine the relationship between subject clustering results and anatomy (Figure Bii). The re- sorted result show that most neurons of the detected largest cluster come from CAT, while some neurons in LP and VISp are also included. Moreover, although most identified subject clusters are neurons from the same anatomical area, there are some mismatches between these two criteria. Es- pecially, some neurons in CA1 I are grouped into the same cluster with neurons VISp, and also be- tween LP and VISp. This may imply that there are some "functional interactions" between CA1 1 and VISp, and LP and VISp. We also compare the subject clustering results from Poisson bi-clustering model and mixDPFA, and sort the similarity matrices using the same order as in Figure C-ii When assuming the equi-dispersion and fit the Poisson bi-clustering model, there are more mismatches between the detected clusters and anatomy (Figure BC-iii,, which may suggest some spurious inter- actions are detected when ignoring the dver-dispersion. The mixDPFA assumes Poisson distributed spikes and further ignores the potential state changes along the time. The detected clustering struc- tures are more noisy, and especially can hardly distinguish between VISp and IP (Figure (iv- The results are consistent with previous finding of the mixDPFA which shows that the neural pop- ulation may change under different experimental settings, if the static dynamics is assumed (Wei al., 2023). Overall, these results suggest that it is necessary to consider the over-dispersion and time-varying nonlinear dynamics, to obtain unbiased estimate of clustering structures.
For the state clustering structures (histogram and trace in Figure 5 the algorithm detects around 13 clusters, and we show the similarity matrix (Figure 33)) and the maxPEAR estimate (Figure Bp). The stage changes faster for spontaneous (S) than natural movie (N), but the temporal states don't show a clear pattern as in subject cluster in NB- IB-bbclusterrig model. This may suggest similarity between spontancous and natural visual stimulation for temporal response, in addition to spatial connectivity pattern (Wilf et al.) 2017). On the other hand, the Poisson bi-clustering model suggests the earlier and later neurons have different stages. But since the bin size is relatively large in current implementation, the over-dispersion may be significant and the results here may be spurious by ignoring the over-dispersion. Further application to data with "ground-truth states, such as sleep/ wake, may help us study effects of over-dispersion in different models.
Finally, we also show the details of the largest maxPEAR subject cluster (traces of 1.2 norm for latent trajectoreis in Figure 5C). The largest maxPEAR cluster has 21 neurons, which contains 9 neurons from CAL, 7 neurons from LP and 5 neurons from VISp. The spiking counts of these neurons (Figure 3 may suggest periodic pattern, i.e., alternating strong and weak response, in the middle portion of the natural movie epoch. The observed pattern is captured by the latent trajectories that most of these neurons belong to (Figure B)).
5 DISCUSSION
In this paper, we introduce a Bayesian nonparametric method to cluster the neural spiking activity spatially and temporally, according to the latent structures (trajectories) for each neural population. Compared to other clustering method for time series (e.g. distance-based methods), the clustering structures defined by latent trajectories can provide scientific insights for the large-scad scale complicated time series data. Moreover, simultaneously consider the subject and state clustering structures can provide us unbiased and consistent estimates of latent structures (e.g. trajectories) and vice versa.
Although the proposed method can successfully bi-cluster the neural spikes, there are some potential improvements. First, the subject clustering structures are modeled by MFMM, which consider the nature for number of neural populations. However, the uncertainty of clustering results can be large in some cases, and hence it may be better to consider the generalized MFMM (gMFMM), which can provide greater efficiency in the cluster estimation (Frühwirth -ruhhwith--cchntter et al.
202T). Moreover, the common atom specification of gMFMM (Denti al. [023] D Angelo et al. 2023) provide
flexiblity in partitions estimations, resolve the degeneracy isuues and more importantly can allow us borrow information from different neural populations. Second, we currently pre-specify and
8

---

Published as a conference paper at ICLR 2024
A
neural spikes
neuron
D
maxPEAR states
PMN CAR, Nao
CAAT aerori 10
Experiment
180
290
330
304
901
042
CAL
largest maxPEAR cluster
322
402
nNFal sgiene
B
state clustor- NB
Now = 19
MITEAA undaary 14
Intert Injectories
013
64
-2
233
423
802
125
Figure 3: Application in multi-regional neural data under different experiment epochs. A. Here, we apply our method to multi-regional Neuropixels data, which contains neural spikes from 3 regions (CAI, LP and VISp) across periods with different visual stimuli: spontancous (S), natural movie (N) and spontaneous (S). The results from iteration 2500 to 10,000 for each chain are shown here. B. The similarity matrix of state cluster for NB bi-clustering model. C. The similarity matrices of neuron cluster sorted by maxPEAR estimate for NB bi-clustering model (NB- --BB-maaEAAR upper-left). The clustering results sorted by both MB-maxPEAR and anatomical sites for three different clustering models (NB bi-clustering, Poisson bi-clustering and mixDPFA) are shown here for comparison. D. The maxPEAR estimates of the discrete states for NB and Poisson bi-clustering model. E. The largest maxPEAR cluster contains 9 neurons from CAL. 7 neurons from LP and neurons from VISp. The upper panel shows the observed neural spikes. The lower panel shows the latent trajectories that most these neurons belong to, where the blue lines are posterior means and shaded light blue regions are 95% HPD regions.
assume all clusters share the same dimension of latent factors p for convenience. However, this assumption may be inappropriate for real data application, and the method can be more flexible to infer p at the same time. Previously, Wei et al. (2023) sample the latent dimension by a birth-and- death MCMC (BDMCMC) (Fokoué & Titterington 2003, [Stephens] 2000) with the marginalized likelihood, which requires very little mathmatical sophistication and is easy for interpretation. Some other methods, such as putting multiplicative Gamma process prior (Bhattacharya & Dunson 2011
multiplicative exponential process prior (Wang et al. 2016) and Beta process priopr (Paisley &
Carin 2009 Chen et al. 2010) and Indian Buffet process prior (Knowles & Ghahramani] [200]]
2011 Ročková & George 2016) may also be useful. Third, when clustering the temporal state, we
tried to the Hiicky-HDP-HMM to handle the rapid-swit issue. However, the method restrict to geometric state duration and doesn't it allow for learning state-specific duration information. When applying to the Neuropixels data, the state looks change fast. This may suggest that we need to model the state duration more carefully, for example, by HDP-HSMM (Johnson & Willsky| 2013)) Moreover, neither HiCKy-HDP-HMM nor HDP-HSMM allow the transition of discrete latent state & to depend on latent trajectories H. Therefore, it may be possible to combine idea of recurrent HMM (Linderman et al. 2017) with HDM-HSMM, which may lead to some method like HDP-recurrent- HSMM, for instance. Finally, although the MCMC algorithm developed here is quite efficient, deterministic approximation of MCMC, such as variational inference may be more computationally efficient and can be more attractive for scientific application.
To sum up, as the scale of neural spiking data becoming large both spatially and temporally, un- derstanding the latent structures of multiple populations under different conditions can be a major statistical challenge. Here, we provide a way to extract opatio-temporal clustering structure, accord- ing to their lwwddimensional latent trajectories. Although the proposed bi-clustering method is to resolve problems in neuroscience, this method can be potentially useful to extract insightful latent structures (bi-clustering and trajectories) from general large-scale (counting) time series.

---

Published as a conference paper at ICLR 2024
REFERENCES
Mauricio A. Álvarez and Neil D. Lawrence. Computationally efficient convolved multiple output
gaussian processes. Journal of Machine Learning Research, 12(41):1459-1500 2011. URL
httt:://mm...oogg/apprrss / /v12/
httpeeesssy////aavelllleae wwvarezlla.....oo htmm
A. Bhattacharya and D. B. Dunson. Sparse Bayesian infinite factor models. Biometrika, 98(2):291,
2011. ISSN 00063444, doi: 10. 10.1099//10MET//RRRR
Mihail Bota and Larry W Swanson. The neuron classification problem. Brain Research
Reviews, 56(1):79-88) 2007. ISSN 0165-0173. doi: http:///bb.orr//00011..j..
2007.05.005. URL https:/ /tttps:/www.... wtw.:s/ceeeeeeeeee.cccccoceeeees
5016501 7307000663
Phillip Boyle and Marcus Frean. Dependent Gaussian Processes. In L Saul, Y Weiss, and L Bottou (eds.), Advances in Neural Information Processing Systems, volume 17. MIT Press, 2004. URL fttpps://ppoceeednggs
cwtp://ng..nniiii......co.ooocce http:ppppppeeillee/////////// fil File/50 File/59eb5dd3e http:://cccceae...........ccescp
Lars Buesing, Timothy A Machado, John P Cunningham, and Liam Paninski. Clustered factor analysis of multineuronal spike data. In Advances in Neural Information Processing Systems, volume 27, 2014.
Jiachen Cai, Robert J. B. Goudie, Colin Starr, and Brian D. M. Tom. Dynamic factor analysis with dependent gaussian processes for high-dimensional gene expression trajectories, 2023.
C. K. Carter and R. Kohn, On Gibbs sampling for state space models. Biometrika, 81(3):441-553, 1994. ISSN 0006-3444, doi: 10. 1w...///iiom///1.......
Minhua Chen, Jorge Silva, John Paisley, Chunping Wang, David Dunson, and Lawrence Carin. Compressive sensing on manifolds using nonparametric mixture of factor analyzers: Algorithm and performance bounds. IEEE Transactions on Signal Processing, 58(12)61140-6155, 2010.
doi: 10. 10.1009999222000220777966
Mark M Churchland, John P Cunningham, Matthew T Kaufman, Justin D Foster, Paul Nuyujukian, Stephen Ryu, and Krishna V Shenoy. Neural population dynamics during reaching. Nature, 487 (7405):51-55) -6, 2012. ISSN 1476-4687.
doi: 10. 10.1t08//aturel1 1229. https:::: :ttetss://dee w/tp://o...ot. org/
10.1033/naturel112
Benjamin Cowley, Adam C Snyder, Katerina Acar, Ryan Williamson, Byron M Yu, and Matthew A Smith. Slow Drift of Neural Activity Signature of Impulsivity in Macaque Visual and Prefrontal Cortex. Neuron, 108(3):551-567..8 2020. ISSN 0896-6273. doi:
021. sst.aaeee/zrrccae///////86222223
John Cunningham and Byron M Yu. Dimensionality reduction for large-scale neural recordings.
Nature Neuroscience, 17(11):10001100, 2014. ISSN 1546-1726. doi: 10.1038/nn.3776 URL
https:://di..rr///1
Laura D' Angelo. Antonio Canale, Zhaoxia Yu, and Michele Guindani. Bayesian nonparametric analysis for the detection of spikes in noisy caleium imaging data. Biometries, 79(2):1370-1335 jun 2023. ISSN 0006-341X. doi: https://doi.org org/10..00/1000 1/biomm. 13626. URL hhttps://doo
org/ /0/ 1111/biom. 13626
Francesco Denti, Federico Camerlenghi, Michele Guindani, and Antonietta Mira. A Common Atoms Model for the Bayesian Nonparametric Analysis of Nested Data. Journal of the American Statistical Association. (ttttaaioonn((54))))) 118(541):405-116 jan 2023. ISSN 0162- 1459. doi: 10. 1080/01621459 0800
2021. 1933499. URL URL.ht::p::
2021 1933499
József Fiser, Pietro Berkes, Gergö Orbán, and Máté Lengyel. Statistically optimal perception and learning: from behavior to neural representations. Trends in Cognitive Sciences. 14(3):119-130, 2010. ISSN 1364-6613, doi: http:///d...gg///....0.......s.
URL https: //ww.
sciencedi rect. com/ sci http:///cw.ctt.....o....aaaaaaaa
10

---

Published as a conference paper at ICLR 2024
Ernest Fokoué and D. M. Titterington. Mixtures of factor analysers. Bayesian estimation and infer- ence by stochastic simulation. Machine Learning, 50(1-2):73 944, 2003. ISSN 08856125. doi:
10.1023/A.10000998288222.
Emily Fox, Erik Sudderth, Michael Jordan, and Alan Willsky. Nonparametric Bayesian Learning of Switching Linear Dynamical Systems. In D Koller, D Schuurmans, Y Bengio, and L Bottou (eds.), Advances in Neural Information Processing Systems, volume 21. Curran Associates, Inc., 2008a. URL hattts::/proceeeeeee
www.:/t..wiii.ccoo..o.....o cc/paper fi Files/paper/2008/ paf
Emily B Fox, Erik B Sudderth, Michael I Jordan, and Alan S Willsky. An HDP-HMM for Systems with State Persistence. In Proceedings of the 25th International Conference Ma- chine Learning, ICML "08, pp. 312-319, New York, NY, USA, 2008b. Association for Com-
puting Machinery. ISBN 9781605582054, doi: (1114//1390115.139019999 URL https:
Emily Beth Fox. Bayesian nonparametric learning of complex dynamical phenomena. PhD thesis,
Massachusetts Institute of Technology, 2009.
Arno Fritsch and Katja Ickstadt. Improved criteria for clustering based on the posterior similarity matrix. Bayesian Analysis, 4(2):367-391, 2009. doi: 10.1214/09-BA414, URL https://
doi ooi. //1 BA414
S Frühhwitt--Shhnatter,, G Malsiner-Wallh and Grün. Generalized Mixtures of Finite Mixtures and Telescoping Sampling. Bayesian Analysis, 16(4):1279-1348)) 1307, 2021. doi: 10.1214/21-ba 1944.
URL ottp://////o.oo/1............
Sylvia Seathwrtth-chhattee Data augmentation and dynamic linear models. Journal of Time Series Analysis, 15:183-202, 1994.
Peiran Gao, Eric Trautmann, Byron Yu, Gopal Santhanam, Stephen Ryu, Krishna Shenoy, and Surya Ganguli. A theory of multineuronal dimensionality, dynamics and measurement. bioRxiv, pp. 214262, jan 2017. doi: 10.1 101/2142..... URL htt httt::::/bbooxww.ooo///ooooo 2017 2017/11/12 2421222222222... 22t.22112222....aacacco.a
Z Ghahramani and G E Hinton. Variational Learning for Switching State-Space Models. Neural Computation, 12(4):831-64, 2000. ISSN 0899-7667 VO 12. doi: 10.1162/
0899766003000156199
Joshua Glaser, Matthew Whiteway, John P Cunningham, Liam Paninski, and Scott Linder- man. Recurrent switching dynamical systems models for multiple interacting neural pop- ulations. In Advances in Neural Information Processing Systems, volume 33, pp. 14867- 14878, 2020a. URL httt:p:://proceediing wwurin......iss.s.eeee.........s httpssccppppee//2////// aalf5f7 3322 02t.:322aaa0deeeeeeeeeeeeee 785aaf Bapeeaaaaafaaaaaaaaaaaa pdf
Joshua Glaser, Matthew Whiteway, John P Cunningham, Liam Paninski, and Scott Linder- man. Recurrent Switching Dynamical Systems Models for Multiple Interacting Neural Pop- ulations. In Advances in Neural Information Processing Systems, volume 33, pp. 14867- 14878, 2020b. URL Ba1f5f7332 7ba 40 40d44404771e5555555 85aaf Paper. pdf
Haochen Xu Xuening Zhu Guanhua Fang, Ganggang Xu and Yongtao Guan. Group network hawkes process. Journal of the American Statistical Association, 0(a):1-78, 2023. doi: 10.
URL
1000110..100111111 01 20223
2023.
2257889
Tal Harmelech and Rafael Malach. Neurocognitive biases and the patterns of spontaneous cor- relations in the human cortex. Trends in Cognitive Sciences, 17(12))606-115 dec 2013. ISSN 1364-6613, doi: 11001//. 2013 09.014
Matthew J Johnson and Alan Willsky. Bayesian Nonparametric Hidden Semi-Markov Models. J. Mach. Learn. Res., 14(1)6673-701, feb 2013. ISSN 1532-4435.
11

---

Published as a conference paper at ICLR 2024
Eric Jonas and Konrad Kording. Automatic discovery of cell types and microcireuitry from ral comnectomics. eLife, 4zc04250. 2015. ISSN 2050-084X. doi: 10.7554/eLife.042550 URL
https://doi org/ /00 7554 04250
James Jun, Nicholas A. Steinmetz, Joshua H. Siegle, et al. Fully integrated silicon probes for high-
density recording of neural activity. Nature 2017 531:7679, 551(7679):232-233, 2017. ISSN
1476-4687. doi: 10.1058/mature24663.
URL https: www nature.tue........ http://wwwwwtatee..cmaaaaaattttt wwt.///aaaateeeeeeeeeee
nature24636
Matthew Kaufman, Mark M Churchland, Stephen I Ryu, and Krishna Shenoy. Cortical activity in the null space: permitting preparation without movement. Nature Neuroscience, 17(3):440- 448, 2014. ISSN 1546-1726. doi: 10.1038/mm.3643 URL htttps://cci...... http://///o..r.......o 3643
David Knowles and Zoubin Ghahramani. Infinite sparse factor analysis and infinite independent components analysis. In Independent Component Analysis and Signal Separation, pp. 381-388, Berlin, Heidelberg, 2007. ISBN 978-3-540-77444--8.
David Knowles and Zoubin Ghahramani. Nonparametric Bayesian sparse factor models with appli- cation to gene expression modeling. The Annals of Applied Statistics, 5(2B):1534 1552, 201 doi: 10.1214/1--AOASAA
URL /ttpttttttppo...ooeettottttttttt oroii.org// 12141111111110000 A0AS435
Scott Linderman, Ryan P Adams, and Jonathan W Pillow. Bayesian latent structure discovery from multi-neuron recordings. Advances in Neural Information Processing Systems, volume 29, 2016.
Scott Linderman, Matthew Johnson, Andrew Miller, Ryan Adams, David Blei, and Liam Paninski. Bayesian Learning and Inference in Recurrent Switching Linear Dynamical Systems. In Pro- ceedings of the 20th International Conference on Artificial Intelligence and Statisties, volume 54, pp. 914-922, 2017. URL
mlr www.pp..s//a.........
html
Scott Linderman, Annika Nichols, David Blei, Manuel Zimmer, and Liam Paninski. Hierarchical recurrent state space models reveal discrete and continuous dynamics of neural activity in jemaC. elegans;/em,, bioRxiv, pp. 621540, jan 2019, doi: 10. 10.1101/621540 http:::://biirxx botr://ooteee//eerrr2//////////2
abstract]
Jakob H. Macke, Lars Buesing, John Cunningham, Byron M. Yu, Krishna Shenoy, and Maneesh Sahani. Empirical models of spiking in neural populations. Advances in Neural Information Processing Systems, 24, 2011.
Valerio Mante, David Sussillo, Krishna Shenoy, and William T Newsome. -ontexttdependent computation by recurrent dynamics in prefrontal cortex. Nature, 503(7474):78 84, 2013. doi: 10.1038/nature12744 URL https:// http:////o..rr//111111111
ISSN 1476-4687.
na nature12742]
Najja Marshall, Joshua I Glaser, Eric M Trautmann, Elom A Amematsro, Sean M Perkins, Michael N Shadlen, L Abbott, John Cunningham, and Mark M Churchland. Flexible neural control of motor units. Nature Neuroscience, 2022. ISSN 1546-1726. doi: 10.1038/ bintttss://di..rgggg/ http://...101///// 1031/ 022-5599--22--1--101
941593-022--1166 165-8.
Jeffrey W. Miller and Scott L. Carter. Inference in generalized bilinear models, 2020. URL https:
7/ //arxiv. ort/ /2010.200000 04895
Jeffrey W. Miller and Matthew T. Harrison. Mixture models with a prior on the number of
components. Journal of the American Statistical Association, 113(521):340, 2018. doi: 10.10800006224599206665555
Kevin Patrick Murphy. Machine learning: a probabilistic perspective. MIT press, 2012.
Simon Musall, Matthew Kaufman, Ashley L Juavinett, Steven Gluf, and Anne K Churchland. Single-tring trial neural dynamics dominated by richly varied movements. Nature Neuroscience, 17700 2019. ISSN 1546-1726. doi: 10.10338/41599--119-000222 URL https: 7/doi. hawploorg/100.0 -0-
12

---

Published as a conference paper at ICLR 2024
Radford M Neal. Markov Chain Sampling Methods for Dirichlet Process Mixture Models. Journal of Computational and Graphical Statistics, 9()):249-265, 2000.
Sacha B Nelson, Ken Sugino, and Chris Hempel. The problem of neuronal cell types: a physio-
logical genomics approach. Trends in Neurosciences, 29(6):339 345, 2006. ISSN 0166-2236. doi:
200.0.0......22000000000000..... URL https: wwww http://wwww.ciieccrreeee....
com/ http://cceearttteee///////622222
Mariann Oemisch, Stephanie Westendorff, Stefan Everling, and Thilo Womelsdorf. Interareal Spike-Train Correlations of Anterior Cingulate and Dorsal Prefrontal Cortex during Atten- tion Shifts. The Journal of Neuroscience, 35(38):13076 LP 13089, sep 2015. doi: 10. 1523/INEUROSCI.12222 38/13076.abstr..
URL http:/ http:///wwwwwwwo........c
John Paisley and I Lawrence Carin. Nonparametric factor analysis with beta process priors. In Proceedings of the 26th Annual International Conference on Machine Learning, ICML '09, pp. 777-784, New York, NY. USA, 2009. ISBN 9781605585161. doi: 10.1145/153374415554444
URL htttp::::///o..coo org 10. 11..1...1113377771 1113447454
Panagiotis Papastamoulis and George Iliopoulos. An Artificial Allocations Based Solution to the Label Switching Problem in Bayesian Analysis of Mixtures of Distributions. Journal of Computational and Graphical Statistics, 19(2):313-331, jan 2010. ISSN 1061-8600. doi: /tt.g//1..111 11t.::..11.......0/.0 09008
1..11998jcss22010000000..
Panagiotis Papastamoulis and loannis Nızoufras. On the identifiability Bayesian factor ana- lytic models. Statistics and Computing, 32(2):23, 2022. ISSN 1573-1375 doi: 10.1007/
111222-022-1 022-002222000--4.. URL
111-11222222222000-----
Nicholas G. Polson, James G. Scott, and Jesse Windle. Bayesian inference for logistic models using pólya-gamma gamma latent variables. Journal of the American Statistical Association, 108(504):3339 1349, 2013. doi: 10. 1021100001621459.2133.322222
URL httpps://cci.orr////
01 01621459.20m 20144999211... 8290011
Arezoo Pooresmaeili. Jasper Poort, and Pieter R Roelfsema. Simultaneous selection by object-based attention in visual and frontal cortex. Proceedings of the National Academy of Sciences, 111(17): 6467 6477, apr 2014. doi: 10. 2..1073/paa..3111111 pnas. 1316181111
UtLhttps::/ /ellattpps://c..orr//// 10777/111111
Raquel Prado and Mike West. Time Series: Modeling, Computation, and Inference. Chapman & Hall/CRC, 1st edition, 2010. ISBN 1420093363.
Carlos E Rodríguez and Stephen G Walker. terministic Relabeling Strategies. Journal of Computational and Graphical Statistics, 23(1):
Label Switching in Bayesian Mixture Models: De-
45. jan 2014. ISSN 1061-8600. doi: 10. 00.10000106116002222.36666/ //doi.org/ 10
URL https:
2012. 72112.73562
Veronika Ročková and Edward I. George. Fast bayesian factor analysis via automatic rota- tions to sparsity. Journal of the American Statistical Association, 111(516): 1622, 2016. doi: 10. 0100000162145.2015.5
URL http::://a.i..cgg//0000000//////
2015. 1100620
Douglas A Ruff and Marlene R Cohen. Attention Increases Spike Count Correlations between Vi- sual Cortical Areas. The Journal of Neuroscience, 36(28):7523 LP - 7534, jul 2016. doi: 10.
1523/INEUROSCI.6110 16.2016. URL http::// http:::/wwwwwweeureeecoooo jweurosci................ wttp://w.tg...oooeeeeoooomooooee
28/7523 28/7523.abstractt ract
Sepideh Sadaghiani and Andreas Kleinschmidt. Functional interactions between intrinsic brain ac- tivity and behavior. Neurolmage, 80:379-386 2013. ISSN 1053-8119. doi: https://doi.rr//10
N016/i.euroomagee2013.000 100. URL peatpps:/wweee www.eiiiiineieieettt............ wttp://ett..oocccccccccccc/ aiticle/pii/s1053811
John Sansom and Peter Thomson. Fitting Hidden Semi-Markov Models to Breakpoint Rainfall
Data. Journal of Applied Probabiliry, 38:142-157, sep 2001. ISSN 19029002. URL http:
7/ jat.r..rrrro... htwp.//stor.rr//ttbe//////////
13

---

Published as a conference paper at ICLR 2024
Shreya Saxena and John P Cunningham. Towards the neural population doctrine. Current Opin-
ion in Neurobiology, 55:103-111, 2019. ISSN 0959-4388. doi:
2019.02.002. URL https:/ wwwp//wwwww..o// www.sccennediiect.roo com http://wt....c/seeeeeeeerrtll/aa
509594388 8300990
João D. Semedo, Amin Zandvakili, Christian K. Machens, Byron M. Yu, and Adam Kohn. Cor- tical areas interact through a communication subspace. Neuron, 102(():249-255 2019. ISSN 0896-6273. doi: https://doi.orr//00. w0w......oo.......ooo...m.
.026. URL httpp: //ww:/
SCi scciencedinn rect www.co///cc.oooo htt:/scienccce httm:/cienneeearriiccio / 9300534
Jayaram Sethuraman. A CONSTRUCTIVE DEFINITION OF DIRICHLET PRIORS. Statistica
Sinica, 4(2):639-650, nov 1994. ISSN 10170405, 19968507. URL http://////// http://wwwr jstor.
arg/stable/2430555
Joshua H. Siegle, Xiaoxuan Jia, Séverine Durand, et al. Survey of spiking in the mouse visual system reveals funetional hierarchy. Nature 202/ 592:7852, 592(7852):86-92, 2021. ISSN 1476- 4687. doi:
nat wwt. wtt..///t.a.ttee/c// URL https: http://www.aaatee...mrrr.....ae.
151586-020-03177
A' A'yunin Sofro, Jian Qing Shi, and Chunzheng Cao. Regression analysis for multivariate dependent count data using convolved gaussian processes, 2017.
Sergey D Stavisky, Jonathan Kao, Stephen I Ryu, and Krishna Shenoy. Motor Cortical
Visuomotor Feedback Activity Is Initially Isolated from Downstream Targets in Output-Null Neural State Space Dimensions. Neuron, 95(1):195-2000 2017. ISSN 0896-6273. doi:
StRLhhhts:://// www. www.ncedirennee..cc.c....... com/
8ci Ssteenee//rrttcee///////66666712
Nicholas A Steinmetz, Peter Zatka-Haan Matteo Carandini, and Kenneth D Harris. Distributed coding of choice, action and engagement across the mouse brain. Nature, 57667786):26--73, 2019. ISSN 1476-4687. doi:
http:t/:::/ /tttps::/doi... /ttt:::://ooi..//11... /10
Nicholas A. Steinmetz, Cagatay Aydin, Anna Lebedeva, Michael Okun, Marius Pachitariu, Marius Bauza, Maxime Beau, Jai Bhagat, Claudia Böhm, Martijn Broux, Susu Chen, Jennifer Colonell, Richard J. Gardner, Bill Karsh, Fabian Kloosterman, Dimitar Kostadinov, Carolina Mora-Lopez, John O'Callaghan, Junchol Park, Jan Putzeys, Britton Sauerbrei, Rik J. J. van Daal, Abra- ham Z. Vollan, Shiwei Wang, Marleen Welkenhuysen, Zhiwen Ye, Joshua T. Dudman, Barundeb Dutta, Adam W. Hantman, Kenneth D. Harris, Albert K. Lee, Edvard Moser, John O'Keefe, Alfonso Renart, Karel Svoboda, Michael Häusser, Sebastian Haesler, Matteo Carandini, and Timothy D. Harris. Neuropixels 2.0: A miniaturized high-density probe for stable, long-term brain recordings. Science, 37266539), 2021. doi: 10. 00.1126sscienceaab445588 URL https: http://wwwwwww.cee//tttt.
www.nce.eccc........o.... wtt. wwc.iinne.org/bb//as////o. 1ew.ll1266cc6 sci
Matthew Stephens. Bayesian analysis of mixture models with an unknown number of compo- nents-an an alternative reversible jump methods. The Annaly of Statistics, 28(1))40 74, 2000. doi: (0.1214/a0s/001612000
Carsen Stringer, Marius Pachitariu, Nicholas Steinmetz, Charu Bai Reddy, Matteo Carandini, and Kenneth Harris. Spontaneous behaviors drive multidimenssinal, brainwide activity. Science, 364464437)eaav7893,, apr 2019, doi: 10. httl///ieneee.aay//B 78893.
1126/sciencn 7a93
Valérie Ventura. Traditional waveform based spike sorting yields biased rate code estimates. Pro- ceedings of the National Academy of Sciences of the United States of America, 106(1 (06(17)66221- 6926, 2009. ISSN 00278424. doi: 10. (0070//NNA..001111111
PDF. URL wwws. org. http//gi///o// wtwggg/////////eeooooooooeooeeee http://///////ooottee/t..... full/
Saurabh Vyas, Matthew D Golub, David Sussillo, and Krishna V Shenoy. Computation Through Neural Population Dynamics. Annual Review of Neuroscience, 43(1):249-211 275,, 2020. ISSN 0147- 006X. doi: hwwp:mmuwuuuuuuuuuuum09999991999
URL https: //cci http//ii..ooogo///
Enmurev-reuron (0te::o00991119---
14

---

Published as a conference paper at ICLR 2024
Ye Wang, Antonio Canale, and David Dunson. Scalable geometric density estimation. In Pro- ceedings of the 19th International Conference on Artificial Intelligence and Statistics, vol- ume 51 of Proceedings of Machine Learning Research, pp. 857-865, Cadiz, Spain, 2016. URL pittps:/prrceest:
ml m press v51 /aagg be html
Ganchao Wei, Ian H Stevenson, and Xiaojing Wang. Bayesian clustering of neural spik- ing activity using a mixture of dynamic poisson factor analyzers. In S. Koyejo, S. Mo- hamed, A. Agarwal, D. Belgrave. Cho, and Oh (eds.), Advances in Neural In- formation Processing Systems, volume 35, pp. 19499-19510. Curran Associates, Inc., 2022. URL. htttps://proceedd P File/ http://9994511222233eecccc333333 P
neuri ct/ wwtpscpppppppceeeeeeceeeeeeeeeee http://perrrfllessppppr/////////
Ganchao Wei, lan H. Stevenson, and Xiaojing Wang. A flexible bayesian clustering of dynamic subpopulations in neural spiking activity, 2023.
Meytal Wilf, Francesca Strappini, Tal Golan, Avital Hahamy, Michal Harel, and Rafael Malach. Spontaneously Emerging Patterns in Human Visual Cortex Reflect Responses Naturalis- tic Sensory Stimuli. Cerebral Cortex, 27(1):750-763, jan 2017. ISSN 1047-3211. doi: 10. 10.1099/cercor/bb227. URL httph:tts//eoooo org.
cercor/bbo225
Jesse Windle, Carlos M. Carvalho, James G. Scott, and Liang Sun. Efficient Data Augmentation in Dynamic Models for Binary and Count Data. 2013. URL /KK/http://arrziv.rr///bbb
0774
Hongteng Xu and Hongyuan Zha. dirichlet mixture model hawkes processes for event sequence clustering. In I. Guyon, U. Von Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett (eds.), Advances in Neural Information Processing Systems, volume 30. Cur- ran Associates, Inc., 2017. URL wtt.:///......teee/////peeee.ee.
wwtpsper.pppppeeeeeeeeeeeeeeeeee files/ 25207700//1111dee777aaaaaaaee000 pap
Lihao Yin, Ganggang Xu, Huiyan Sang, and Yongtao Guan. Row-clustering of a point process- valued matrix. In M. Ranzato, A. Beygelzimer, Y. Dauphin, P.S. Liang, and J. Wortman Vaughan
(eds.), Advances in Neural Information Processing Systems, volume 34, pp. 20028-20039. Cur- ran Associates, Inc., 2021. URL /ttp/t/pss//prreeeeen...eeeeeee
cc/paper files/
Byron M. Yu, John P. Cunningham, Gopal Santhanam, Stephen Ryu, Krishna V. Shenoy, and Ma-
neesh Sahani. paussisn-process factor analysis for low-dimensional single-trial analysis of neural
population activity. Journal of Newrophysiology, 102(1)):14-635, 2009. ISSN 15221598. doi:
URL https:
Shun-Zheng Yu. Hidden semi-Markov models. Artificial Intelligence, 174(2):215 2010.
ISSN 0004-3702. doi: http:///o.aan////..000o/mmmmmm.0 URL https://www.
Hongkui Zeng. What is a cell type and how to define it? Cell, 185((15)223999775553 2022. ISSN
0092-8674. doi: 10.1016fj.cll12202200... 2022.06 031
Mingyuan Zhou, Lingbo Li, David Dunson, and Lawrence Carin. Lognormal and Gamma Mixed Negative Binomial Regression. In Proceedings of the 29th International Coference on Interna- tional Conference on Machine Learning, ICML' 12. pp. 859-866, 2012. ISBN 978 1450312851.
David M Zoltowski, Jonathan W Pillow, and Scott W Linderman. A General Recurrent State Space Framework for Modeling Neural Dynamics during Decision-Making Proceedings of the 37th International Conference Machine Learning, ICML'20. JMLR.org. 2020.
15

---

Published as a conference paper at ICLR 2024
A APPENDIX
A.1 MCMC UPDATES
The posteriors of the model parameters are sampled by MCMC algorithm. To illustrate the sam- pling steps for each iteration, we expand them into six steps and provide the details as follows.
A.1.1 UPDATE LATENT FACTORS
Denote the p.m.f. for (B, as
fNBly [1,))= T(r v) (r)
9
r+
p
then the full conditional distribution of X1.r is as follows:
P P(XIT R.TT (yr,.r, de (br, Ar, (1):
N T
T
(N de. dis
- 11
|
II
1 exp (-26% Qe.1 Qc,11-s)))
exp
2
1=1
where Pat - exp (d: + c. and 8: X, As, - The full conditional distribution has no closed form, and we sample it via PG augmentation technique, i... by introducing the auxil- iary PG variable, the transformed "effective" observation dar performs like Gaussian, and hence the (FFBS) can be implemented to sample XI.r- Specifically, we can sample XI.r from full conditionals two steps:
regular sampling algorithm like
Step I: sample the auxiliary PG parameter wa and calculate the transformed response ga- For i = 1,.... N and L = 1.... T, sample the wa from Pólya-Gamma distribution as Wit ~ PG re + Yie. d; + E, - log rit Then we can calculate the transformed response as jar Wit 2 Ku. where K2t (ya, r2)/2 + wallogr; - dd,)
Step 2: sample the X1.T by FFBS. Since the transformed response Da N(d,ā,, w ) we can use the regular FFBS algorithm for Gaussian state-space model, the detailed algo- rithm can be found Chapter 4 in Prado & West (2010).
For data doesn't tt have severe over-dispersion, it can be useful to assume Poisson distributed obser- vations (i.e. assume Vir ~ Poisson(μ4z)), for ease of the model inference. The Poisson distribution doesn't 't fall within the PG augmentation scheme, and we cannot use the method described here. However, motivated by the fact that NB distribution approximates to Poisson distribution as goes to infinity, we can use the sample from NB model as the proposal, and add one more Metropolis- Hasting step to accept or reject it as in (Wei et al. 2023).
To ensure the model identifiability, we project the posterior samples to the constraint space, such (j)
that 1) a(() is zero-centered, i.c. ΣT a( = 0 and 2) x #1:7"217 L:T is diagonal, where x 11.7 =
=0) (())
(a)() ab² E (+++))) (Fokoué& Titterington, 2003). However, the model still not identi- For deterministic algorithm,
fiable in terms of the signed-permutatiin on singular value decomposition (SVD) as in Miller & Carter (2020). but this is not appropriate for the sampling algorithm. To resolve this, simply search for the higned-permuteen the closest Euclidean distance to reference trajectories (in this implementation, we use sample in the 100-th iteration. Before that, align current step to previous). Instead of fixing the reference trajectories manually, we may find the optimized value recursively using the Varimax Rotation- Sign- Pegmmeermttation Warimax-RSPI algorithm developed by (Papastamoulis & Ntzoufras 2022).
can put constraints based
that has
A.1.2 UPDATE SUBJECT BASELINE AND FACTOR LOADING
The sampling of subject baseline d, and factor loading c; from full conditional distribution is regular NB regression problem for each neuron, which is again updated by PG augmentation technique. The idea is the same in sampling latent factors, i.e., transform the spikes yir to be Gaussian like
16

---

Published as a conference paper at ICLR 2024
by introducing the augmented parameters wie ~ PG (rs + yien Hi (())) + (1, (1 (d,, ")) - -logr:) Therefore, Di Wit Ka 2f (d, ci), wi) where Kig (ya r:)/2 wa w1oglog7
H4 (a:)).
A.1.3 UPDATE DISPERSION
The dispersion for each neuron r, is updated via a Gibbs sampler, since the gamma distribution is the conjugate prior to it, under the compound Poisson representation. Specifically, let pat Pa (M + ri), then the conditional posterior of is Gamma ao + Σi= T,, 1/(1 Σ log(1 pst)) where L ~ Poisson --r; log(1 - -pa)). Refer to Zhou et al. 2012) for more technical details such as how sample Lt.
A.1.4 UPDATE DISCRETE STATES
To update the discrete states Çt, we use the weak-limit Gibbs sampler for sticky HDP-HMM in (Fox al. 2008a), by constructing finite approximation to the HDP transitions with finite Dirichlet distribution. This is motivated by the fact the infinite limit of hierarchical mixture model converges in distribution to true HDP as M + oc, such that
B ~ Dir(n/M,..... n/M),
1/1 Dir(a aBt).
By using this weak limit approximation, we can update the states by an efficient, blocked sampling algorithms. Refer to Fox et al. (2008b) for more details.
A.1.5 UPDATE LINEAR DYNAMICS
Because the Gaussian assumption in the model, updating the linear dynamics for each state is Bayesian multivariate linear regression problem. Specifically, let Xi = (1r, X1:r), then for state I assume the conjugate priors for (b), A,, Qi) as
Q: W (Vo 70), vec((b,, AD) NVvee((B) Qi 1
where W-1 denotes the inverse-Wishart distribution. The priors are set as Vn 0.01 1 7o (p+ 1) + and B, (0, rY. If there are Q-chunks latent factors having & I, and the denote the time steps for q-th chunk as Tg kg Tq. Then, the full conditional distribution are:
Q: XIT W- 7n), vec((bi, Ai)
X1.7 N(B,, Qi Q 1
where
V, Vo
S SyS (B, B...o(B, Bo),
Sq X, r
ho +1:,
T, -kyi Bn
7n = 70 +
kig
1=1
B, r-1
X
- + Bo)
9=1
Q
r, Σ
x -
T'
To
9=1
17

---

Published as a conference paper at ICLR 2024
A.1.6 UPDATE SUBJECT CLUSTER ASSIGNMENTS
To update the subject cluster labels 2,S, use partition based algorithm, similarly Miller & Harrison (0118). Let C denote a partition of neurons, and CV denote the partition obtained by removing neuron from
1. Initialize C and (c)
C) (e.g., one cluster for all neurons in our simulation).
2. In each iteration, for = ,,...,.. N: remove neuron from C and place it:
(a) in c E Cii with probability x (lc) + y Mc(.....T)) where y is the hyperparameter of the Dirichlet distribution in equation 3 and )MM, denotes the marginalized likelihood of neuron i in cluster c, when integrating the loading &n out.
(b) in a new cluster e* with probability x
Va((a) VV,s)
Va(++11 Mainr) where s is the number
(w])ent (!), with x(m) a<x + 1) ...................... (r + - 1), =(m) = x(x - 1) (x: m + 1), x0()) = and 2(0) 1.
of partitions by removing the neuron and Vn(ss) Di
The update is an adaptation of partition-baae based algorithm for DPM (Neal 2000), but with two sub- stitutions: 1) replace [cg] by |c. + 7 and 2) replace a by (V I)/V, (t). See more details and
discussions in (Miller & Harrison 2018).
Instead of evaluating the full likelihood, we integrate the subject-specific factor loading é; ob- tain the marginalized likelihood (M(.i.1.)) to achieve better mixing for high dimensional situation. The marginalized likelihood is evaluated by the Laplace approximation as follows:
T
Mo(MIITT)
X
INNNI
ri, ri, (T
(c) Pit
d, + Pit E,x
where 2 means "approximately proportional to",, é and E., are MLE estimates and corresponding variance estimates (inverse of the negative Hessian at 6;) from NB regression on C₁-
Since we need to align the latent trajectories to reference value for each cluster, it's necessary to handle the label switching problem for cluster assignment. Here, we implement the Equivalence Casses Representative (ECR) alforithm (Papastamoulis & Iliopoulos] 2010)), by setting the cluster assignment at the 100-th iteration (before 100-th step, align current allocation to previous step) as the pivot allocation. To relax the dependency setting pivot manually, can further consider iterative versions by /Rodriguez & (21k))
A.2 SUPPLEMENTARY RESULTS FOR SIMULATIONS
Here, we show trace plots (Figure 4 of several parameters for NB bi-clustering model in simulation (section B and Figure 2) Visual inspcctoo and Geweke diagnostics don't 't show significant issues with convergence.
A.3 SUPPLEMENTARY RESULTS FOR NEURAL PIXELS
In this section, we first provide histograms and trace plots for parameters in the MCMC chain for NB bi-clustering model (chain of NB bi-clustering model) in shown in Figure Again, visual inspection and Geweke diagnostics don't show significant issues with convergence.
We then provide results from two independent chains for both NB and Poisson bi-clustering model (i.e., four chains in total). Specifically, show the similarity matrices of subject clusters sorted by maxPEAR estimates (Figure (A) and by maxPEAR with anatomical sites (Figure 63). These orders obtained from first chain of NB model, which the same order as used in Figure 3C. We further show the similarity matrices for state cluster for these four chains. Overall, these results show that NB model provides different clustering results compared the Poisson one.
18

---

Published as a conference paper at ICLR 2024
# of subject cluster
B
# of state cluster
10.0
7.5
5.0
2.5
2000
4000
6000
BO00
10000
2000
4000
6000
8000
10000
iteration
iteration
C
trace of ||
D
trace of ||X]|r
now rolated
20
10
-20
50
10
0
n21 n30
10
2000
4000
6000
8000
10000
2000
4000
6000
8000
10000
iteration
iteration
Figure 4: Trace plots for simulations. Here we show trace plots of several parameters for NB bi-clustering model in simulation (section and Figure including number of subject cluster (A)
and number of state cluster (B). We further show the traces of || llz (C) and || (D) for the 3 detected clusters (equivalent to true subject cluster). The way we define detected cluster is (5)
illustrated in section and Figure 2E. For ||=|// [|x, the traces for both the raw (blue) and rotated
(orange) samples are shown here
of neuron cluster
of state cluster
largest maxPEAR cluster
trace of ||
1500
2000
7.5
1500
1000
1000
2.5
500
trace of [|xaella
10
15
10
15
15
15
trace of [Xxl
10
2000
6000
10000
2000
6000
10000
iteration
iteration
iteration
Figure 5: Histograms and trace plots for chain of NB-biclustering. A. Histograms for number of neuron/ state cluster, using the samples from iteration 2500 to 10,000. The corresponding trace plots are shown in panel We then show the L2 norm for each latent trajectory, i.e. μ,, X(1) and ((2), for the largest maxPEAR cluster, corresponding to latent trajectories in Figure 3E. Here μ denotes PI. (o) denotes i-th row of 1
19

---

Published as a conference paper at ICLR 2024
A
cluster maxPEAR
B
cluster
state cluster
maxPEAR anatomy
Possen-tckole
Noarn
Polsor clsty chiin
Paisson Pelionbeluioo chain
Figure 6: Supplementary results for the application to Neuropixels data. A. similarity matrices for neuron cluster, for independent MCMC chains each on both NE and Poisson bi-clustering models. These matrices sorted by maxPEAR estimates of chain for NB model (also used in main text, Figure i). B. We further sorted the matrices by anatomical sites, which is the same order used in the main text, Figure BC ii-iv. C. The similarity matrices for state cluster, for independent chains on each model.
20


Reviews:
[{"review_id": "nu9zXTbB8o", "reviewer": "Reviewer_96hg", "summary": "This paper proposes a nonparametric Bayesian biclustering algorithm to simultaneously form spatial (neuron) clusters and temporal clusters for multiple (count) time series. The algorithm is carried out through efficient MCMC sampling and is shown to be able to recover all model parameters in simulation studies. The proposed algorithm is then applied to a real data set to show its effectiveness.", "questions": "1. The authors only reviewed some work in time series clustering. However, the neural spiking activity data is originally a point process and one has to aggregate them using small bins into count time series. There is some existing work on finding clusters in point process literature as well, which in my opinion should be reviewed for relevance. For example, [1], [2], [3], and references therein. \n\n2. The simulation study is a little bit too simple, with only one parameter setting, and there are no comparisons to the state-of-the-art methods.   I would suggest investigating three things: (1) the sensitivity of the bin sizes; (2) how the estimation performances change when the number of nodes and the Time length increases; (3) compared to some existing methods in terms of the predictive performances (one may consider using cross-validations.)\n\n3. In the real data analysis, the authors argue that it is necessary to use the biclustering algorithm. However, no comparison to any existing method is given.  Similar to the simulation studies, I suggest comparing the predictive performance of the proposed algorithm to some existing methods.\n\n4. In the simulation study, there is no histogram in Figure 3(c). Could you please clarify?\n\n[1] Xu, H., & Zha, H. (2017). A dirichlet mixture model of hawkes processes for event sequence clustering. Advances in neural information processing systems, 30.\n\n[2] Yin, L., Xu, G., Sang, H., & Guan, Y. (2021). Row-clustering of a Point Process-valued Matrix. Advances in Neural Information Processing Systems, 34, 20028-20039.\n\n[3] Fang, G., Xu, G., Xu, H., Zhu, X., & Guan, Y. (2023). Group network Hawkes process. Journal of the American Statistical Association, (just-accepted), 1-78.", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The algorithm is clearly described, reasonable assumptions are imposed, and it has a potentially large range of applications.\nThe overall writing is good and the presentation is clear.", "weaknesses": "The numerical results, including simulation studies and real data application, are not enough convincing.", "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": "", "extra_scores": {}}, {"review_id": "o7w1IISG30", "reviewer": "Reviewer_xUve", "summary": "Motivated by modeling of neural spiking activities, the paper proposes a bi-clustering approach at both spatial and temporal levels. Individual count data are modeled through negative binomial regression. The spatial (subject) clusters are modeled through mixture of finite mixtures. The temporal (state) clusters are modeled through hierarchical Dirichlet process. The dependency of negative binomial regression on the temporal dimension is given by a hidden Markov model. A MCMC algorithm is designed for sampling from the posterior distribution. Both synthetic numerical experiments and real applications are provided.", "questions": "The temporal clusters are captured in the temporal dynamic in the AR(1) model. What is the motivation of such modeling from an application perspective, compared to simpler modeling (e.g. directly model temporal clusters of \\tilde{X}_t through change points, etc)?", "limitations": "", "rating": 8, "confidence": 5, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The problem is well-motivated with a meaningful and important application.\n2. The paper is mostly well-written and clearly presented.\n3. Sufficient background and preliminaries are provided.\n4. Details on the derivation of the MCMC algorithms are provided.\n5. The challenging goal of conducting full Bayesian inference for a complex clustering task is of itself great importance.", "weaknesses": "1. While some constraints required for identifiability are provided at the end of section 2.1, I am not convinced that these are sufficient conditions. A theoretical proof of the model identifiability along with all necessary conditions seems important here, given the vast number of parameters.\n2. The MCMC algorithm has incorporated all the most modern efficient MCMC techniques, including the Polya-Gamma augmentation, Miller&Harrison sampler for mixture of finite mixtures, FFBS for state space models, etc. The effort here is worth being recognized. However, for both the simulation and application, the MCMC sampler is run for only 1000 iterations. I am concerned of whether the MCMC sampler has really converged given its complexity and the vast amount of model parameters. If it has indeed converged, it would be nice to present trace plots of parameters and also conduct MCMC convergence diagnosis (e.g. using Gelman-Rubin or Geweke statistics).\n3. There are some typos in the paper and some terms undefined, e.g. \n    - the distribution \\mathcal{S} in equation (4) is only given in the appendix\n    - in the first paragraph of Sec 2.1, where \\log\\mu_{i, t} = d_i + \\tilde{c}' \\tilde{x}_t^{(z_i)} has missing subscript i in \\tilde{c}.\n\nThe general idea of this paper and the designed MCMC sampler are both quite nice. I would be willing to raise the rating by one or two levels if the points 1 and/or 2 can be fully addressed.", "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": "", "extra_scores": {}}, {"review_id": "XTk0Ep4BDR", "reviewer": "Reviewer_DWGp", "summary": "This paper proposes a spatio-temporal clustering method to analyze multiple neural populations. This method could find the clusters in neurons (spatial) and different states over time (temporal). To flexibly extract the bi-clustering structure, the authors model the spike data in a non-parametrical way, where the subject clustering structure is modeled by a mixture of finite mixtures model and the state clustering structure is modeled by a sticky Hierarchical Dirichlet Process Hidden Markov Model. The inference is performed by MCMC with the Polya-Gamma technique. In the experiments, the authors evaluate their method on both simulated data and neural recordings.", "questions": "*  What's the time complexity (concerning the number of iterations, trials, neurons, and time points) of the proposed model with the efficient MCMC algorithm? \n\n* What kind of desirable neural data could be better analyzed by such spatiotemporal clustering structures rather than LDS-based models (e.g., SLDS, r-SLDS)? \n\n* Could you please discuss the relationship of this method to the PP-Seq model [1]? PP-Seq is an unsupervised non-parametrical model to detect neural sequences in high dimensional neural recordings, which could be considered as a clustering approach to finding spatiotemporal neural patterns.\n\n[1] Williams, Alex, et al. \"Point process models for sequence detection in high-dimensional neural spike trains.\" Advances in neural information processing systems 33 (2020): 14350-14361.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "* Different from previous LDS methods, this paper explores multi-region neural data from a new perspective: the authors try to understand multi-neural populations by spatiotemporal clustering structures.\n\n* Non-parametrically model the neural data so that there is no need to prespecify the number for subject and state clusters.", "weaknesses": "* No analysis of the scalability of this method. For both syntactic data and real neural data, the number of neurons is small (e.g., 30, 60). Could this method generalize to a large neural recording? Such as a larger number of neurons and a longer time stamps.\n\n* No comparison of the proposed model with other latent variable models like SLDS and rSLDS.\n\n* Some typos, e.g., the \"cite\" doesn't refer to a paper in section 2.3, the figure index should be 2 rather than 3 in section 3, and the figure index should be 3 rather than 4 in section 4.", "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": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (poster)