Venue: NeurIPS 2021 Conference

Reviews:

final_score: rating: 7: Good paper, accept
review_content: summary: This work develops a set of tools for studying higher order interactions in neural activity recordings. Population neural activity is represented as a marked point process, and fit using a time-varying multinomial GLM. Further, the authors present hypothesis testing tools that compare the fit model with a reduced model that models neurons independently. Specifically, the authors present a Deviance measure that can be used in testing whether r-th order interactions are present in data. 
main_review: * This paper gives a comprehensive set of tools to study higher order interactions. The tools would be very valuable for computational neuroscience audience. Most of my comments are about making some practical/application aspects clearer in the text. 

* While there are no apparent inaccuracies in the text, but the main bulk of paper is heavy in mathematical symbols and equations. Improving the writing could greatly improve the readership and appeal.  Defining the notations in a separate subsection in beginning, and stating all the modeling assumptions (time-varying models, sparsity assumptions, smoothening steps) at same place would help. 

* Review of previous works that use marked point process tools for neural activity is missing. 

*  For modeling time history dependance, the authors switch from $n^*_t$ (marked representation) to $n_t$  (original neural activity) for computational tractability. This should be explained more carefully since it might have significant implications on results / interpretation. 

* How sensitive are measures/models to timing jitter? Small timing changes in individual spikes could change the assigned mark. For example, for two neurons and two nearby time bins, an observation of [(0, 0), (1, 1)] could easily become [(1, 0), (0, 1)] by moving the first neuron's spike ahead in time slightly. Since the mapping from actual spike train to marked point process is highly non-linear, this could significantly affect results of the methods.

* While the presented method gives a statistical test to test if there is rth order interaction in a population of C neurons, it does not identify which neurons are interacting. Repeating the test on all subsets of $r$ neurons might not be ideal. Expanding the discussion around this issue would clarify the application for neuroscientific questions. 

* How do the empirical results (especially on human data) change with the smoothening parameter $\beta$ ? Showing dependance on/robustness to $\beta$ would help understand the time-varying model better. 

* In general, since the approach uses a time-varying model, the methods presented in the paper do not identify higher order effects that occur at very slow time scales. For example, for the data in Figure 2, considering very large time scales (large $W$ or $\beta$ ~ 1) would show significant r=2 order interaction for the whole duration of experiment. Choice of these hyper-paramters is important for practical application and should be mentioned in text.

* Large differences between history-dependent and independent models is revealing, suggesting that large order interactions in simple models can be attributed to simple mechanisms in slightly more complex models. This is explained in main text associated with Figure 2. However, making these analyses clearer would help with a broader readership, who want to apply these tools to data. 

limitations_and_societal_impact: no negative societal impact
review_title: null

final_score: rating: 5: Marginally below the acceptance threshold
review_content: summary: The authors describe a time-varying multinomial GLM model which detects higher-order correlations in multi-dimensional spike trains. See main review for further summary.
main_review: It took me a while to understand the contributions made in this paper, but I think they can be summarized as extending the multinomial GLM model of Ba et al. (2014, Frontiers in Comp Neuro) to have time-varying parameters. This extension leans heavily on methods described in Sheikhattar et al. (2015, IEEE Transactions on Signal Processing). The goal of this new model is to capture higher-order spike time correlations between simultaneously recorded neurons. The new method does not require repeated trials.

My biggest concern with this paper is the difficulty I had in understanding the motivations behind various modeling choices and implementation details. I feel like I only began to understand this paper after reading Ba et al. (2014) and Sheikhattar et al. (2015) in careful detail.

For example, a window length W and a forgetting factor \beta control how quickly the model parameters vary over time. If the forgetting factor is close to zero and the window lengths are relatively small, then the parameters can vary rapidly, potentially resulting in overfitting -- how is this avoided? If the forgetting factor is close to one, then the parameters vary only very slowly, potentially resulting in an underfit model. Even if the forgetting factor is well tuned, this model can capture a gradual drift in parameters, but potentially it will have trouble accounting for rare, step-like changes in network state, as happens during propofol induced loss of consciousness (investigated in Fig 3). The paper would be stronger if the dynamics of the time-varying parameters were more clearly discussed and illustrated in practice.

The assumption of sparsity on the GLM parameters is not clearly motivated or discussed. Why not use l2 regularization instead of greedy sparse estimates? The Orthogonal Matching Pursuit algorithm is not reviewed, and terms like "support set" (refering to the subset of nonzero parameters) are used without a clear definition. On lines 171-172 the authors say they apply a Kalman forward filter and backward smoother control "abrupt variations" but the details of this smoothing are not present in the main or supplemental text.

The application of the model to real data is discussed only in a short section (5.2). How were the \beta parameters selected in this analysis? How much sparsity was imposed on the parameter estimates, and how was this tuned?

Despite these misgivings, I think the model addresses an interesting topic and I am open to raising my score if the other reviewers were better able to understand the details of the model.


Other Comments
--------------

- On line 318 there is a typo, "estimate" should be "estimated"

- In figures 2 + 3, please label the blue-to-red colorbar

- On lines 167-169, are the "u" variables meant to be "\mu"?
limitations_and_societal_impact: The limitations of the work and motivation behind modeling choices could be more clearly and concisely discussed. I foresee no potential negative societal impacts.
review_title: null

final_score: rating: 7: Good paper, accept
review_content: summary: The current manuscript presented an algorithm based on marked point processes to distill higher-order interactions in neuron ensembles. The authors provided thorough theory for the algorithm and demonstrated the efficacy of the algorithm in detecting both time-varying spiking rates and coordinated spiking components in simulated and real neuron ensembles. 
main_review: In general, the authors presented a novel approach for analyzing higher-order interactions in a neuron ensemble. The simulation results suggested that the algorithm successfully identify the endogenous components of higher-order coordinated spikes. In application of real neuron ensemble, the algorithm identified brain state modulated higher-order interactions. The paper was well written, and the experiments were carefully designed to demonstrated the power of the current algorithm. However, I would concern about the scope of application of this model. First, the anesthesia data has relatively simple structure, that single factor (anesthesia) generated large effects on neuron spiking rates. And the authors only tested the model on eight neurons with the highest spiking rate from the ensemble. It is unclear whether the model would converge for a data set with more complication structure and lower spiking rate. Second, the identified higher-order interactions and the history-dependence did not readily interpretable, as the modeling components that generated this knowledge did not have a biology meaning. In many real-world applications, users would not be able to tell whether identified higher-order interactions and history dependence were real or meaningful. 
limitations_and_societal_impact: The authors did not discuss the limitation of current work. I think one of the major limitations of the current model is lack of functional foundation or biology meaning of model components.
review_title: null

Meta-review: This paper presents a model for multi-neuronal spike trains called the mGLM. The observation in time bin $t$ is a binary vector of length $N$, where $N$ is the number of neurons, and it is modeled as a categorical random variable that can take on one of $2^N$ values. The categorical probabilities are modeled via a GLM, with past spiking and external inputs as covariates. The authors then allow for time-varying parameters and propose an OMP-based fitting procedure. 

The reviewers were generally favorable of this paper. However, as Reviewer KgmN pointed out, the core idea of a multinomial GLM to capture dependencies in instantaneous counts is one that has been explored previously (Ba et al, 2014), so the main technical advance here is in the time-varying parameters and the OMP fitting procedure. I'd also note that many have recognized the limitations of conditional independence assumptions in standard GLMs for neural data; indeed this is a main motivation for latent variable models like the Poisson LDS (see Macke et al, NeurIPS 2011, e.g.) and the recurrent linear model (Pachitariu, NeurIPS 2011). I'm also skeptical of the exponential complexity in model parameters in the mGLM, and the time-varying aspect seems to exacerbate the problem. 

Despite my reservations, I will go along with the reviewers and recommend acceptance.
Decision: Accept (Poster)