Venue: NeurIPS 2021 Conference

Reviews:

final_score: rating: 6: Marginally above the acceptance threshold
review_content: summary: This study proposes a statistical model to jointly estimate the mean firing rate and covariance structure based on the spike count of neural populations, and test the method on mouse head-direction cells. 
main_review: The paper proposes a new descriptive statistical model to fit neural data. By using Gaussian process, the model is able to jointly fit mean firing rate (1st order statistics) and covariance (2nd order statistics) simultaneously, which is  the main novelty of this paper. And the author used this model to fit synthetic data and real neural data, and compared the model performance with other alternative models used in neural data analysis. The study is solid, and figures look clean and great. However, since some key math notations are not defined yet, I am not sure the math used in this study is correct or not although I hope so.

### Writing:
The writing and the definition of math notations can be improved a lot in a revised version. For example, some underlying motivations of proposing a math equation are lacking, and some math notations are not defined at all when using.

- I don’t understand what are observed X and latent Z, and what is the difference between the two variables. From the generative model shown in Fig. 1, the observed X seems to be a latent variable generating recorded spike count Y. Or X is regarded as the presented stimulus to an animal? Please give a concrete example of what X and Z can be referred to in neurophysiology experiments and otherwise it is not easy to understand.
is the motivation to include those two variables.

- Line 100: It is not unclear what $\Pi$ is. It is not defined. The $\Pi$ seems a lumped variable including all observations but I am not sure.
- Eq.8: what is \phi^{-1}? it is not defined.
- Eq.11: it is not clear what F(q) is. What is the difference between $F(q)$ and $q$? Is $F(q)$ the empirical distribution while $q$ the parametric distribution defined in Eq. 2?
- Eq. 12: please add the motivation of including $1/T+ 1/3T^2$ in the equation. 

limitations_and_societal_impact: Yes
review_title: null

final_score: rating: 6: Marginally above the acceptance threshold
review_content: summary: The authors present a very broad framework to model multi-neuronal spike trains (multivariate time series over the nonnegative integers). The model class is sufficiently broad that it can express any arbitrary distribution over N-dimensional vectors containing the integers {0, ..., K}, and thus model the activity of N neurons within a single time bin (assuming an upper bound on the number of spikes, K).
main_review: This paper is well-written, the figures look polished, and I believe the core claims of the paper are correct. Overall, I would like to see this paper published if the authors can address my concerns below, which mostly deal with the interpretation of the model, relation to existing methods, and framing of the paper. In short, I think the model could be useful in cases where one would like to specify large time bin sizes, but the utility of the method is less clear in the common case where one takes the time bins to be very short.

Most existing methods, the authors claim, have a "fundamental underlying assumption of Poisson spiking statistics." The model they propose does not have this assumption if time bins are large and K > 1. However, if we consider the limit of very small time bins (on the order of 1 millisecond) then it is reasonable to set K=1, since the probability of a neuron firing more than one spike is negligible. The multinomial distribution in thier model becomes a binomial distribution. In the limiting case of of time bin duration going to zero, you recover Poisson statistics (see, e.g., sec 5.2 of Kass et al. "Analysis of Neural Data"). In this limit, the model proposed by the authors looks very similar to existing models, like GPFA or GLMs.

If the above analysis is correct, then I think the authors should acknowledge it in the paper and explain that their method is a generalization that retains "universal" modeling capabilities when time bins are not very small and when K > 1. Then, they should explain what the relative advantages and disadvantages are to choosing larger time bins.

There are other parts of this paper that I really like, such as the combination of supervised learning and latent variable modeling to capture representational drift and head direction tuning. However, these modeling innovations don't seem to be sold as the key results of the paper and it doesn't seem like the "universal" spike count model is really necessary to enable this analysis. Again, it seems like one could choose small time bins and use GPFA (w/ poisson noise) + GLMs to arrive at a similar result.

Overall, I thought this paper was very well-written and presented some elegant analysis of real experimental data, so I have scored it as a marginal accept. The modeling innovations that are most heavily emphasized in the manuscript seem to be dependent on having large time bin sizes. If my interpretation here is wrong I am happy to raise my score higher.
limitations_and_societal_impact: The limitations of the model are clear and I do not foresee any negative societal impacts.
review_title: null

final_score: rating: 7: Good paper, accept
review_content: summary: The paper presents a model of spike count distributions in a population of neurons. The model’s spike count distribution depends both on stimulus information and latent variables through a Gaussian process to capture a large range of joint spike count distributions (as the number of factors in the model increases). Moreover, the authors present a set of model fit statistics for evaluating the model fit, and how correlations factor into the model’s predictions.
main_review: The model improves upon previous non-parametric approaches (such as the universal binary count model) to modeling spike count distributions in neural population recordings: importantly, the proposed model includes stimulus information and temporal correlations. The example applications are very thorough in examining the model fit and comparisons to the Poisson and negative binomial models.

The presentation was generally clear. However, I am not convinced I understand the procedure for quantifying correlations given on lines 183-190. I’m certain I’m overthinking things here, but I do think this one portion of the methods could be presented more clearly.

Minor point: one limiting factor not mentioned in the discussion is that the stimulus term in the examples here are low-dimensional compared to stimuli in sensory systems.
limitations_and_societal_impact: Yes
review_title: null

final_score: rating: 6: Marginally above the acceptance threshold
review_content: summary: In this paper authors provide a method for modelling spike counts as a function of input covariates and latent covariates. A Markovin prior is assumed over latent parameters, and the effect of latent parameters and input covariates on spike counts are modelled using a Gaussian process. The authors also suggest a method for calculating Z-scores and tuning index based on their method. The results include analyses of a synthetic dataset and an experimental data, which generally shows the proposed method outperforms considered baselines. 
main_review: The paper is generally clear (but see below). The statistical model is interesting and suitable for modelling the joint effect of latent and observed covariates on the neural activities. (the main text needs to be updated based on the typo mentioned in the supplementary materials). The results are also encouraging and show that the model is improving on baseline methods.

I have the following concerns about the study that require further improvements:

Some of the results are presented in a rather confusing way. For example cvLL results are presented relative to the Poisson model in some panels (e.g,, Fig 2A) but not in other places (e.g., Fig 3 in supple). In Fig 2A for example it's hard to see whether differences with the performance of the Poisson model is significant or not. Also RMSE is reported for synthetic experiments, but not for the experimental data.

The baselines are also different between the experiments; the performance of ANN is not reported for the experimental data. Also, I couldn’t find how the architecture of the neural networks was selected; how many neurons? How were hyper-parameters tuned? This is particularly important since recent approaches report promising results using recurrent neural networks, e.g. below, and the reflective performance of the suggested method w.r.t.to these methods is unclear.

Pandarinath, Chethan, et al. "Inferring single-trial neural population dynamics using sequential auto-encoders." Nature methods 15.10 (2018): 805-815.

Sections 2.4 and 2.5 are rather hard to follow. For example it would be good to expand why the definition of Z in eq 8 is a sensible one and how it is related to the common definition of z score. Below eq 11, for example it's unclear what does data being more/less regular than the data mean.

limitations_and_societal_impact: They are adequately discussed.
review_title: null

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