Dynamic Analysis of Higher-Order Coordination in Neuronal Assemblies via De-Sparsified Orthogonal Matching Pursuit
Authors: Shoutik Mukherjee
Behtash Babadi
Venue: NeurIPS 2021 Poster (2021)
Keywords: neural coding
generalized linear models
synchrony
orthogonal matching pursuit
high-dimensional inference

Abstract: Coordinated ensemble spiking activity is widely observable in neural recordings and central in the study of population codes, with hypothesized roles including robust stimulus representation, interareal communication of neural information, and learning and memory formation. Model-free measures of synchrony characterize the coherence of pairwise activity, but not higher-order interactions; this limitation is transcended by statistical models of ensemble spiking activity. However, existing model-based analyses often impose assumptions about the relevance of higher-order interactions and require multiple repeated trials in order to characterize dynamics in the correlational structure of ensemble activity. To address these shortcomings, we propose an adaptive greedy filtering algorithm based on a discretized mark point-process model of ensemble spiking and a corresponding precise statistical inference framework to identify significant coordinated higher-order spiking activity. In the course of developing the statistical inference procedures, we also show that confidence intervals can be constructed for greedily estimated parameters. We demonstrate the utility of our proposed methods on simulated neuronal assemblies. Applied to multi-electrode recordings of human cortical ensembles, our proposed methods provide new insights into the dynamics underlying localized population activity during transitions between brain states.

Dynamic Analysis of Higher-Order Coordination in Neuronal Assemblies via De-Sparsified Orthogonal Matching Pursuit
Shoutik Mukherjee and Behtash Babadi Department of Electrical and Computer Engineering
Institute for Systems Research University of Maryland, College Park
College Park, MD 20742 (smukher2, behtash) Centtashh@umd... edu
Abstract
Coordinated ensemble spiking activity is widely observable in neural recordings and central in the study of population codes, with hypothesized roles including robust stimulus representation, interareal communication of neural information, and learning and memory formation. Model-free measures of synchrony characterize the coherence of pairwise activity, but not higher-order interactions; this limitation is transcended by statistical models of ensemble spiking activity. However, existing model- base analyses often impose assumptions about the relevance of higher-order interactions and require multiple repeated trials in order characterize dynamics the correlational structure of ensemble activity. To address these shortcomings, we propose an adaptive greedy filtering algorithm based on discretized mark point- process model of ensemble spiking and corresponding precise statistical inference framework to identify significant coordinated higher-order spiking activity. the course of developing the statistical inference procedures, we also show that confidence intervals can be constructed for greedily estimated parameters. We demonstrate the utility of our proposed methods on simulated neuronal assemblies. Applied to multi-electrode recordings of human cortical ensembles, our proposed methods provide new insights into the dynamics underlying localized population activity during transitions between brain states.
Introduction
Coordinated ensemble spiking has been observed in a variety of brain areas, prompting a range of hypotheses about its role cognitive function. Studies have documented synchronous spiking at all
levels of the mammalian visual pathway [ 2 3 Coordinated neural activity has additionally been
hypothesized to influence interareal communication and the flow of neural information [45 6) 7 g.
and been postulated to be mediated by oscillations in local field potentials [9] TO 11
The study of
synchrony also closely tied to memory [12 14.
The prevalence of coordinated spiking and its functional implications for range of neural processes have motivated both model-free and model-based approaches its characterization. An intuitive model-free metric is the pairwise correlations of spike trains smoothed by Gaussian (or exponential) kernel [15 16]; several pairwise distance metrics have been proposed [17] as alternatives. Though the coherence of pairwise activity can be described, such measures do not capture higher-order coordination, and are limited in the ability to model dynamics or determine the significance of pairwise coherence without repeated trials.
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

Statistical models of neuronal ensemble activity transcend the limitation of model-free metrics to pairwise comparisons. Two widely-used approaches are the maximum entropy models and point process generalized linear models (GLM) [8 19 Maximum entropy models describe the state of the neural population only in terms of instantaneous correlational structure [20] 21]. Models are estimated to match observed firing rates and all pairwise (and potentially higher-order) correlations simultaneously. Alternatively, point process GLMs for ensemble spiking [22] 23
characterize the influence of past population activity, or other relevant covariates. Though useful in estimating functional connectivity [24 25], each neuron must be assumed conditionally independent due regularity conditions that prohibit simultaneous spiking events [26] 27 28) This can be circumvented by using an equivalent marked point processes (MkPP) representation that explicitly models each disjoint simultaneous spiking event, derived in 22 and expounded upon [27]. MkPP representations coding in unsorted spiking data 30 31] related approach models disjoint simultaneous spiking events as log-linear combinations of point process models, hence permitting an intuitive representation of excess suppressed synchrony [10
ensemble activity have also been utilized to analyze neuronal population
Though statistical models can capture higher-order neural coordination, existing approaches face key limitations. Maximum entropy models can track dynamics in coordination using state-space filtering algorithms, but neglect the influence of past population activity on the ensemble state. point process models address this shorteoming, but still share two shortcomings with maximum entropy models. First, assumptions on the relevance of higher-order interactions are typically imposed for tractable model estimation. Second, multiple repeated trials required to capture dynamics in correlational structure and to evaluate the statistical significance of coordinated spiking.
We address these limitations by proposing adaptive greedy filtering algorithm based on the discretized MkPP formulation (27) to model dynamics coordinated spiking within continuous recordings while capturing the influence of past ensemble activity. Furthermore, we build on recent theoretical results related to Adaptive Granger Causality analysis [24] to provide a precise statistical framework to detect significantly coordinated activity of arbitrary order. We demonstrate our proposed method's 's utility tracking dynamics in coordinated spiking with statistical confidence on simulated ensemble spiking, Applying method to continuous multi-electrode recordings of human cortical assemblies during anesthesia provides novel insights into coordinated spiking dynamics that underlie transitions between brain states.
2 Preliminaries
In this section, we review the discrete-time marked process representation and two corresponding likelihood models. First, we present summary of key notation used throughout the subsequent sections in Table
Notation
Definition
Ensemble spiking observation at time bin t of C neurons
A
Conditional Intensity Function (CIF) of th neuron
Marked observations at time bin of C = 2 - marks
A (m)
CIF of mh mark
(9)
Ground process, D T0
X A0(0)
CIF of the ground process, I X=(m)
Base rate parameters of mark events
W.
(ay
Model parameters of history-dependeer model
(m)
Log-odds of mh mark event vs. spiking event Log-odds of mn mark event vs. no spiking event (restricted model) Exogenous factor for mhth mark Forgetting factor, Window length
14
(m)
10.0
(m))
[m]
up.
19
<
W
Table 1: Summary of key notation

---

2.1 Marked Process Representation of Ensemble Spiking
To characterize coordinated spiking, we utilize the discrete-time marked process (MkPP) representa- tion of ensemble neuronal activity [27] 28). For an ensemble of C' neurons, the C-variate spiking process, binned with small bin size A, at time bin index denoted by n. == [n, ne ..... n. 1 where each component is the spiking process of one neuron. Conventional discrete point process models treat the components as conditionally independent Bernoulli observations. Given our interest in simultaneous spikes, we instead treat nu as multivariate Berouulii observations. The spiking pro- (C"))) 'I which are the binned ..... "4
(2)
cess is mapped to -"-variate process n. == [n; (1), (2))) observations a marked point process whose marks count the number of exactly one of C := 2 disjoint non-zero spiking events; we refer to nł as the marked Bernoulli process, distinguishing from the multivariate Bernoulli process n2e We define the mark space := {l,.... C*) 26] Fig. shows example of mapping the activity of C = 3 neurons to C* 7 marked processes. At each time ty such that nu, + O, the sole non-zero element of mi, indicates the mark. We also define the
ni
binned ground process nt (x) that takes value 1 at each such 1 and is zero otherwise 226 ; the ground process indicates the occurrence of any spiking event and is represented by n
Em=1 ni (m)
Spiking Data
Simu Inneous Spiking Events
Spiking Coordination
Latent Highec-Order
inile netron
7
II
F1
I
IB
o p-lrreer
Neuson
Slcueon
Dynamic Higher-Orde Coordination Analysis
Netson
Coordination Strength: J-Statistic
Dynamic Inference of --order Coordination
Adaptive Discretized
MkPP Model
Figure 1: Ensemble spiking is mapped to disjoint representation of simultaneous spiking events. The proposed
method is used infer the strength of higher-order coordination amongst neurons in dynamic fashion. The marked process representation is not unique, but can be defined in convenient fashion: treating the components n, as the bits of C'-bit binary number, the mark indexed by the decimal equivalent of a particular realization of will corresponded to that realization. By the disjointness of the marked representation, the spiking process of the th neuron can recovered as the sum of all marked process whose index, in binary, takes value at the ct bit. For instance, in Fig. the spiking data of neuron (in blue) is the sum of simultancous spiking event processes 4-7.
Our main contribution in this work is the dynamic and statistically precise inference of latent coordinated spiking of C neurons using their simultaneous spiking representation (Fig. bottom panel). To this end, we next describe two useful likelihood models for simultaneous spiking.
2.2 Two Likelihood Models of Simultancous Spiking
In the discrete formulation, the conditional intensity functions (CIFs) of n, and nž approximated by the probabilities of observing an event at time bin i given the ensemble's spiking history. That is,
A (c) A 1P(r) (c) 1/Hr),
X (m) A
PPnn((m)
1|Ht],
(c) for c = I,.... and m = 1, 1,.... C*. We can relate X A to Xp (m)) A in the same manner as n. to nf (m) and obtain the CIF of the ground process X;
A > (m)
A.
The marked process permits a generative description of simultancous spiking events: ensemble spiking events are characterized by the ground process, occurring with probability >; 2² A; the event then assigned to the th mark (i.e. the mhth simultaneous spiking outcome) with conditional
probability (m)m) A Thus, at time / the likelihood of ensemble event nł is given by:
(m))
p(n)) = X (m) A
())()) m(9) (1-17 1- 1-n(0)
II
X [
m=1
3

---

The likelihood in (2) is used to form multinomial generalized linear model (mGLM) with multi- nomial logistic link function of which we consider two versions. The first makes the simplifying assumption that there is no history dependence; the resulting model depends only on contempora-
neous spiking, permitting compact parameterization by μr 144 : .....P .r '1. Defining the baseline firing parameter for the mh
(2) (C*)
mark to be
X (m)
H4 (m)
log
1,2,
(3)
oPt (m)
or equivalently (m A
the log-likelihood be rewritten linear function of
1+ e"!
n., resembling the maximum entropy model 20
log p(n.) = p!n! - (pp)), where ((4)) == log 1
Σ ea!
(4)
m=11
The second, more general version utilizes the ensemble history covariates in the mGLM. Letting the covariate vector F; be the ensemble history up to some fixed lag at time (augmented by a constant element of 1), the model is parameterized by w, := w !2 ".....,("")
where the parameters for the m²th mark w²!") := [A4Tm] "("""' consists of an ensemble history modulation vector 0!m) in addition the baseline firing parameter. Thus, the log-likelihood in this case admits similar
(C*)) I1,
form to (4), by simply replacing pi (m)) with z/w(""
Adaptive Estimation of the History- Deppedeet mGLM
Unlike conventional mGLM models, here the parameters are allowed to change in time. To capture their dynamics, we take a similar approach to the dynamic history-independenr model of 32 and extend it to the history-dependent mGLM.
We assume conditional independence across time bins and that the parameters w, admit piece- wise constant dynamics and are constant over consecutive windows of length W. The ensemble history up history covariate vectors at the ;th window are denoted by X; = (1) Øaw]'. Since the mapping from nł to n, is injective, the influence of past spiking activity be equivalently captured by defining history covariates in terms of either, however, using n, reduces the dimensionality quantifies the influence of past spiking activity directly rather than through categorical variables. Let (m)) ( ( n,w "(m)) denote the sequence of outcomes of the m mark in the 7th window. n;
lag p defines the covariates as Φ; := [1, n - -4-1: T2 (C) The
(1)
(C)
w; and
The log- likelihodo of the 5th window is thus given by
WW log
C
(((w;) = Σ (m) n,
X,w
(m)
Σ
+
Σ m=1
(m)
(5)
m==
(= 1JW
Motivated by the RLS objective function 33 forgetting factor mechanism is utilized to combine the log-likelle liiikkikkhhods up to the htth window, capturing the dynamics each mark's rates. For forgetting factor 0 < B < 1, the adaatively-weiihhee log- ligeeikli at window k is thus defined as:
(ww) = (1 B) Σ g (((w))
(6)
Parameter estimation can be performed by solving a sequence of maximum likelihood problems:
W == arg max (wxx). k 1,2, K.
we
Two issues arise when considering large ensembles. First, the dimensionality of μu grows exponen- tially with C: second, it is likely that some marks will not contain any events. To address this, we employ a thresholding rule similar to [34,, considering only "reliable interactions". i.e. the subset of the mark space & {m E K. E, n (m) No.r } for some pre-defined constant Ner > 0, and treating the rates of the remaining marked processes negligible due to their infrequency.

---

To efficiently solve the sequence of problems in 7 in an online fashion, we develop an adaptive greedy
approach based on generalized Orthogonal Matching Pursuit (OMP) 33 36.. The OMP iteratively identifies the parameter support set (the non-zero components) of fixed size (a hyperparameter for which we (ross-validatt) to encourage sparsity, thus capturing the inherent sparsity of network interactions based on past ensemble activity [37 38 39] Moreover, the greedy estimation over sparse subset of parameters mitigates the intractability of the estimation problem for large ensembles where regularization-based constraints still require optimization over all parameters. The proposed adaptive OMP (AdOMP) algorithm, so named because the support set is permitted to change between windows, is detailed in Algorithm 1 in Appendix A.I. The key element of AdOMP is efficient evaluation of the gradient V (wxx) at the 7th iterate (())) k- to determine the next addition to the parameter support set and to solve the new maximization problem via gradient descent. Hence, its recursive computation crucial for the algorithm to operate an online fashion. To this end, we
utilize a recursive update rule to compute the gradient at the htth window, generalizing the adaptive filtering techniques employed in 40 for Bernoulli observations multivariate setting.
4 Statistical Inference of Higher-Order Coordination
Coordinated spiking can indicate relationships between components of a neuronal ensemble and, potentially, effects of unobserved processes. However, simultaneous spiking events can still occur by chance amongst independent neurons, necessitating test of significance to distinguish between excessive (or suppressed) and chance simultaneous events. In this section, we detail such a framework: first, we quantify the two alternatives by constructing nested hypothesis test; second, we generalize the de-sparsifying procedure for f--regularized maximum- likeeloooo estimators established by [41] to the AdOMP; and third, we use the latter to establish precise statistical inference framework by proving the applicability of an adaptive de-biased deviance test, used for identifying significant Granger-causal influences (24), to our setting.
4.1 Hypothesis Test Formulation for -th-Order Coordinated Spiking
We characterize the significance of all -onder sinultaneous spiking events for the history-dependent mGLM. The corollary for the history-independee model is addressed in Appendix C.4. The sig- nificance of r-wise simultaneous spiking for r 2 2 is tested by considering the two alternatives:
Ho -order simultaneous spikes occur as frequently as they would between independent units, given ensemble spiking history
(8)
H1
th -order simultaneous spikes occur at a significantly different rate than they would between independent units, given ensemble spiking history
A similar formulation is used in (28) to determine whether one mark occurs at a significantly different rate than expected. The likelihood of the mark is modeled the product of marginal likelihoods times an additional multiplicative factor. Noting that the additional factor takes value 1 if the neurons are truly independent, the null hypothesis is quantified accordingly. To account for all marks of order r, we instead estimate a reduced model that assumes yhhh ooder interactions are chance occurrences by constraining the base rate parameters for each yTth -oder mark. For the mhth mark, (m))
(mn) +76 let uf ((m) := xz. wo (m) = Us = + _ (M) We decompose the base rate parameter as pa (m) = Po, (m)
where Po is rate under the null hypothesis and 7/ is analogous the additional multiplicative factor in (28 that captures potential exogenous effects after conditioning on ensemble spiking history.
We thus estimate the reduced model a/ == arg max, E (w] (H) where the base rate parameters of th events are constrained to those under the null hypothesis. That is, for each m E K,, == {m E K: E.= m. = r). where m. is the th least significant bit of m in binary, we fix u to Po and optimize the remaining parameters. To explicitly obtain the constraints, first recall that x (u("") is the log-odds of n (m)) versus n[(9) given ensemble spiking history. Under the assumption that the neurons are independent, the probabilities of each event is given, respectively, by
P[n (mm)
1|24d]
II
(
II
(1
((a) A). and P(n 0/71.]
II c=1
(1 (1) A). (9)
GA =1
=0
5

---

Taking the ratio evaluated at the full model estimate Wy. obtain up,t Eem.-1 log (1~A"A
(m))
Assuming the difference between u. (mn) and (m) is due only to exogenous factors, we estimate the corresponding term at the beth window as %/(m) = w Ei 21=(t-1)w+1 (u(")) ). Thus, for the reduced
model, we fix (m) at Mk
(m) (m) for m E Ky. The hypotheses are then quantitatively stated as:
Ho Ho: =
H1 Wk + 26 (R)
(10)
To control the possible abrupt variations of m/(m) across windows, apply Kalman forward/backward smoothing to the exogenous factor and use the smoothed values, ((()) in lieu of </m") The procedure is summarized in Algorithm 2 in Appendix A.3.
4.2 De-Sparsitying AdOMP Estimates
To test the hypotheses defined above, is necessary be able to construct confidence intervals for the parameter estimates. The procedure is well-established for unrestricted or unregularized linear regression models, but there is paucity of work this end for greedily-estimatee high-dimensional sparse models. In the closcly-related problem of ( f1-regularized maximum maxiimun-lkkcliood estimation, set of elegant results 41
have established techniques de-sparsify parameter estimates and construct confidence intervals. In particular, we extend the de se-ssssfiicatiin technique of [41], based on close inspection the Karush- Kuhn-Tuche conditions, the greedy high-dimensional setting. In Appendix A.2, we derive the de- sparsiisif AdOMP parameters following s" iterations as
Wh == (1*)) - (F(
(Ve) ((0))))))))
(1I)
Next, the asymptotic normality of the de-sparsified AdOMP estimates is established. While related result is established in [41], the independence of each realization of the covariates and observations is assumed; additionally, several conditions involved are tailored for é é-regularized maximum likelihood estimation. Hence, we adapt the treatment in [11 for AdOMP to establish the following result:
Theorem 1. Consider the maximization of the total data log- aigliikllhood (we) at the h window, where the true parameter We E IRd is (s. comppesssibee with & < 1. Let w be the maximum likelihood estimate and Ws be the AdOMP estimate after Osslog (1)) iterations. 1f conditions (CI)- (C6) are met, the de-sparsified AdOMP estimate dd satisfies
1 1+ 8
1- (wr, wh) = Vv op() 1,
V
where, as 3 1, V AN 2T with IF² - the inverse of the Fisher information matrix.
For brevity, the technical conditions (C1)-(C6) are omitted here, but are presented in Appendix C.2 along with detailed proof. Based on Theorem confidence intervals for the AdOMP estimates can be constructed by adapting the recursive procedure of [40] to our setting.
4.3 Deviance Difference Test for tth-Order Coordinated Spiking
Classical results on likelihood ratio tests between two nested hypotheses 444 45 have established the use of the deviance difference as a common procedure. However, they are ill-suited in our setting due covariates and forgetting-factte mechanism in the data lig-liiklihood. These
to the highly- issues are addressed in related context (24) for the inference Granger-causal links by defining the adaptive de-biased deviance difference and characterizing its limiting distribution under presence and absence of Granger-cause links. We similarly utilize the adaptive de- deviance difference
Di (
1++3 a0 = 1 1-
(6 {((
(2((((
(12)
(R)) are the respective biases of the full and reduced models. As show in Appendix A.I, the full and reduced ligeliill can also be computed in an online
as the test statistic, where ((F)) and (
fashion, in similar manner as the gradients.
The limiting distributions of the adaptive de-biased deviance difference for the joint model under both the null and alternative hypotheses are characterized in a similar fashion to 24 by utilizing the asymptotic normality of the de- searssiiiie AdOMP estimates from Theorem 1:

---

Theorem 2. Let ( and (i respectively be the full and reduced greediy-estimated mGLM parameters at window k. where (R) assumes conditionally independent 1thh -order simultaneous spiking. Then, as B + 1,
i) if thth -order coordination matches independent -order interactions given ensemble spiking history, then D (00 (F) ) 4 x2 i.e. chi-square, and
ū) if yth-order coordination diverges from independent thh - interactions given ensemble spiking history, and assuming the base rale parameters thh-dder interactions scale at least as Ol then D k.B (ook) (,(()) xx((1)) 1Vx {²)), ), i.e. non-central chi-square,
(r) is the non-centrality parameter at window and the degrees of freedom M == [/G]] where Uk the difference in the cardinalities of the full and reduced support sets.
| is
detailed proof is provided in Appendix C.3. In order to fully characterize the limiting distribution of D( under H1, we must estimate the non-centrality parameter for each window. Assuming the
parameter evolves smoothly in time, we use a state-space space smoothing algorithm (24) to estimate it
from the observed D/ values. This not only allows us to identify significant coordination, but to
also quantify the degree of significance using Youden's 's J-statistic
(r) = 1 - a FFFFFFF2aaaa XM () (F x3( (1-ai) a))
(13)
for significance level where F(-) denotes the CDF. Values of Jk close to 1 imply that the rejection of the null is stronger indication of coordination than for smaller values of Jg. Thus, the J-statistic characterizes the test in terms of both type I and type II errors. By convention, we take Jk = when Hg not rejected at the 1.th window. Under the alternative, it is possible to observe either significant excess or suppressed coordination; this can be reflected in the J-statistic by incorporating the net exogenous effect on thth -order coordination and using a signed J-statistic vk( -snn( Emek, 1mm)). The full procedure for identifying significant -oddr coordinated spiking summarized by Algorithm 3 in Appendix C.
Applications
5.1 Simulated Ensemble Spiking Data
We validate our proposed methods in simulated example, performing complementary history- independent and history-dependen analyses. Let the base rate parameter and exogenous effect for the history-independent model be denoted by l4 and yx: and the same for the
model by MR,H and Yk. H. with moduulaion parameter Og. Then, the reduced model constraints imply T = 7 + +/00p. If the observed rate of higher-order events is equal to that of independent neurons, Thk 0; however, higher-order interactions may still be coordinated, i.c. Yk.H -²,0, o O. Conversely, the observed rate of higher-order events may differ from that of independent neurons, i.e., 7 + 0. If 7 O. observed coordination can be attributed to the effects of ensemble history; otherwise, observed coordination was driven by an unobserved process. A MATLAB implementation of both algorithms is provided in supplementary material.
Ensemble spiking of five neurons was generated by a marked Bernoulli process as described Eq. (2 with an average rate of ~ 0.1 spikes per bin. In the first and third simulated epochs, -rder spiking events were excited by amplifying the default epoch, the base rate parameter was increased to induce -order spiking events. These adjustments respectively reflected simultancous spiking induced by ensemble history and by an unobserved process. Figure 2 shows the simulated spiking activity, from which no obvious coordination is observable. The aggregate th order marks are visualized in Fig. B, with apparent increased rates of 3rd_ and th order spiking events.
parameters. In the second
For comparison, we also used three single-trial measures of coordinated spiking. The first is the average Pearson correlation between smoothed spiking responses. The second, is the spiking regularity, quantified by average coefficient of variation (ratio of the standard deviation to the mean inter-spike interval) [46]. A ratio close to 1 indicates Poisson statistics; larger ratios indicate greater
7

---

variability due to self-exciting dynamics while smaller ratios indicate regularity in spiking (i.e. globally coordinated spiking), Both measures are computed over non-ovellapping windows of 250 samples to track dynamics. The third measure is the average difference between rthh -order mark CIFs and probabilities of yTh -odder independent interactions, generalizing the measure employed in [27] to higher-order simultaneous spiking. Other model-based analyses require multiple trial repetitions and are thus unsuited to our single-trial simulation setting.
Statistical analyses of -ohder coordination for ......55 using the history-independent model (W = 10: / 0.975) reveals facilitated -oddrr coordination during the second epoch, indicated by large positive values of the J J-statistics (Fig. 2). Facilitated 4th -oderr coordination is detected during the first and third epochs. Ensemble spiking was also analyzed using the distere-dependent model (W = 10; 3 = 0.99) (Fig. 2-D). Conditional facilitation of 3rd -order coordination was correctly detected during the second epoch and 4th-order coordination was correctly conditioned out. The history-dependent analysis also detected conditional suppression of 2nd -order coordination.
In contrast, the three control measures are unable to capture the underlying dynamics. Significant pairwise correlations (Fig. 2-E) are stably indicated throughout the simulation, insensitive to changes in coordinated spiking across epochs. Similarly, the spiking regularity (Fig. (2FF) indicates Poisson spiking statistics rather than coordinated activity. The 3rd_ and 4th -order mark CIF differences (Fig.
A
Simulated Spike Raster
nomnm indes
B
Morder Event Raster
00 1 IUIH
111 111
111
C
_ Hllstory-Independen
D
History-Dependeen Dependeen
E
Avg. Correlation
F
Avg. Spiking Regularity
G
Avg. Mark CIF Difference
0.15
3rd and order
0.15
4000
12000
Time Index
Figure 2: Analysis of ensemble spiking with non-overlapping epochs of 3rd_ and th torerr coordination. A. Simulated ensemble spiking of five neurons. B. Sum of the Th ""dooddrr sinultaneous spiking events for r= 2,3, 4, Spiking coordination varies across 3 epochs, demarcated by vertical dashed lines. C. Significant "hh-order coordination neglecting ensemble history. Significant th -order coordination based on history- dependent ensemble spiking model. Statistical testing in C-D performed at level a 0.001. E. Average Pearson correlation with 95% confidence interval. F, Average spiking regularity: coefficient of variation 12 SEM. G. Average mark CIF differences of 3d (green) and 4 orrder (tcal) spiking interactions +2 SEM.

---

A
Spike Raster
Timee
500
Deuron index
B
-order Event Raster
100
C
c History-Inddepedeent
D
x() History- Dependent
E
Avg. Correlation
F
Avg. Spiking Regularity
Figure 3: Coordinated spiking in human cortical neurons is exogenously induced during unconsciousness. A.
Raster of 8 cortical neurons aligned to loss of consciousness (LOC) O The anesthetic administered twice: 0 and at 250 B. Sum of the -ororr simultaneous spiking events for ........ C. Significant th -order coordination neglecting ensemble history. Significant ythh -order coordination based on history- dependent analysis. Statistical testing in C-D performed at level a = 0.001. E. Average Pearson correlation with 95% confidence interval. F. Average spiking regularity: coefficient of variation #2 SEM.
G) weakly reflect the underlying dynamics, but closer inspection reveals the oscillatory nature of this sample-to-sample measure that diminishes its reliability (Fig. 2 insets). We further examine the effect of hyperparameters W and in Appendix B.I, inspect the model goodness-of-ftt in Appendix B.2, and provide a complementary simulation study with more complex dynamics in Appendix B.3.
5.2 Real Data Example: Anesthesia Data
We next present our analysis of human cortical neuronal assemblies during the transition into propofol- induced general anesthesia. The data were retrieved in fully anonymized format with permission from the authors in [47] who obtained written consent from the participants in compliance with the institutional review board (please refer to 47 for details). We employed the proposed algorithms to analyze higher-order coordination and compared them against the average Pearson correlation and spiking regularity. The CIF-based measure is omitted given its highly oscillatory nature, rendering its interpretation uncertain. We analyzed spiking data from one subject, selecting the 8 neurons with the highest average firing rate. In contrast to the simulation in Fig. 2 the average firing rate was 0.05 spikes per bin. In Fig. 3AA, their ensemble spiking activity is shown, aligned to the loss of consciousness (LOC) at 0 Ensemble activity recovered after ~50 when propofol was re- re-administered. The decomposition of ensemble spiking into -tt-order events (Fig. 3 3BB) highlights lower rates of higher-order spiking events.
History-independen analysis (W 10; 0.99) of otddrr spiking revealed high rates of 3rd.. 4th and 5th -odder events during consciousness (Fig. 3 C). Additionally. 2nd-order coordination was suppressed during conseiousness. However, vistory-dependees analysis (W 10; 3 = 0.995) did not identify conditional coordination during consciousness (Fig. 3DD Together, these suggest that rates of higher-order simultaneous spiking events diverged from the rate of such interactions amongst independent neurons, but coordination during consciousness is attributable to ensemble history. The

---

high J-statistic values at the start of the recording were transients related to the initial convergence of the adaptive filters.
After LOC, the rates of simultaneous events matched those amongst independent neurons (Fig. 3-C). However, significant 2nd and -odder conditional coordination were detected during anesthesia. Third-order coordination was persistently and exogenously facilitated; 2nd-order coordination changed from a state of initial suppression to facilitation ~50 Our results show local neuronal networks during both conscious and unconscious states exhibited coordinated spiking, but the underlying mechanisms differed between states, and the differences manifested rapidly in concurrence with LOC.
These dynamics poorly reflected by the Pearson correlation and spiking regularity, both computed over windows of 250 samples (resp., Fig- 3-R and Fig. 3-F). Prior to L.OC, pairwise correlations may have been high but spiking was globally synchronized. Hence, spiking statistics to match Poisson spiking and neurons seem uncorrelated because of irregular low-order coordinated spiking. Slight increases in correlation during anesthesia weakly indicate coordinated spiking. The fluctuating spiking regularity leaves the nature of higher-order coordination indeterminate.
Concluding Remarks
The proposed modeling and statistical inference algorithms constitute a novel approach to studying coordinated neuronal spiking. In contrast to previous model- -ased approaches, the proposed method tailored for the analysis of continuous recordings of neuronal data. We demonstrated that the framework can capture both time-varying spiking rates and the influence of spiking history, and thus can detect endogenously or exogenously induced coordinated spiking.
In developing this framework, we showed that confidence intervals can be constructed around greedily estimate parameters in similar fashion sparsity-regularizes parameter estimates. We found this to be a noteworthy gap in existing literature, as theoretical analyses of greedy algorithms focused instead on guarantees of model recovery. This result enabled us to develop a precise statistical inference framework in which the statistical strength of discovered synchronous spiking can be quantified.
Simulation studies demonstrated the efficacy of our framework in detecting suppressed facilitated coordinated spiking activity. Moreover, in application to spontaneous ensemble spiking during the transition into propofol-induced anesthesia, our proposed method provided greater detail about the correlation structure of local neuronal networks in both the conscious and unconscious states. Addi- tionally, our results reflected the abruptness of the transition between network states by characterizing dynamics in coordinated spiking. The ability to track transitions in higher-order network interactions through adaptive filtering techniques can be used address current gaps in understanding the local mechanisms underlying the emergence of different brain states.
10

---

Acknowledgements
This work is supported part by the National Science Foundation awards no. 1552946 and 2032649 and the National Institutes of Health award no. 1U19NS107464.
The authors declare no competing financial interests.
We thank Emery N. Brown and Patrick L. Purdon for providing neural data (published in recorded from human cortical assemblies during anesthesia.
References
[1] M. Meister, R. O. Wong, D. A. Baylor, and C. J. Shatz, "Synchronous bursts of action potentials in ganglion cells of the developing mammalian retina," Science, vol. 252, no. 5008, pp. 939-943, 1991.
[2] M. J. Schnitzer and M. Meister, "Multineuronal firing patterns in the signal from eye brain," Neuron, vol. 37, no. 3, pp. 199-511, 2003.
[3] W. M. Usrey and R. C. Reid, "Synchronous activity in the visual system," Annual Review of Physiology, vol. 61, no. 1, pp. 435-456, 1999.
[4] E. Salinas and T. J. Sejnowski, "Correlated neuronal activity and the flow of neural information," Nature Reviews Neuroscience, vol. 2, 2001.
[5] M. Diesmann, M.-O. Gewaltig, and A. Aertsen, "Stable propagations of synchronous spiking in cortical neural networks." " Nature, vol. 402, 1999.
[6] C. Rossant, Leijon, A. K. Magnusson, and R. Brette, "Sensitivity of noisy neurons to coincident inputs," Journal of Neuroscience, vol. 31, no. 47, pp. 17193-17206, 2011.
[7] S. Moldakarimov, M. Bazhenov, and T. J. Sejnowski, "Feedback stabilizes propagation of synchronous spiking in cortical neural networks," Proceedings of the National Academy of Sciences, vol. 112, no. 8, pp.
2015.
[8] M. M. Tran, L. Y. Prince, D. Gray, L. Saad, H. Chasiotis, J. Kwag, M. M. Kohl, and B. A. Richards, "Neocortical inhibitory interneuron subtypes display distinct responses to synchrony and rate of inputs," bioRxiv, 2019.
[9] M. J. Jutras and E. A. Buffalo, "Synchronous neural activity and memory formation," Current Opinion in Neurobiology, vol. 20, no. 2, pp. 150-155, 2010.
[10] P. Zhou, S. D. Burton, A. C. Snyder, M. A. Smith, N. N. Urban, and R. E. Kass, "Establishing a statistical link between network oscillations and neural synchrony." PLOS Comput Biol, vol. 11, c1004549, 10 2015.
[11] M. Denker, S. Roux, H. Lindén, M. Diesmann, A. Riehle, and S. Grün "The local field potential reflects surplus spike synchrony," Cerebral Cortex, vol. 21, pp. 2681-2695, 04 2011.
[12] M. C. Zielinski, J. D. Shin, and S. P. Jadhav, "Coherent coding of spatial position mediated by theta oscillations in the hippocampus and prefrontal cortex," Journal of Neuroscience, vol. 39,
no. 23, pp. 4550-4565, 2019.
[13] E. Boran, T. Fedele, Klaver, P. Hilfiker, Stieglizz, T. Grunwald, and J. Sarnthein, "Persistent hippocampal neural firing and hippocampal-cortical coupling predict verbal working memory load," Science Advances, vol. 5, no. 3, 2019.
[14] L. Meshulam, J. L. Gauthier, C. Brody, D. W. Tank, and W. Bialek, "Collective behavior of place and non-place neurons in the hippocampal network," Neuron, vol. 96, no. 5, pp. 1178 1191.e4. 2017.
[15] S. Schreiber, J. M. Fellous, D. Whitmer, Tiesinga, and T. J. Sejnowski, "A new correlation- based measure of spike timing reliability," Neurocomputing, vol. 52-54, pp. 925-931 2003.
[16] S. Haas and J. A. White, "Frequency selectivity of layer ii stellate cells in the medial entorhinal cortex," Journal of Neurophysiology, vol. 88, no. 5, pp. 2422-2429, 2002.
[17] Kreuz, S. Haas, A. Morelli, H. Abarbanel, and Politi, "Measuring spike train synchrony," Journal of Neuroscience Methods, vol. 165, no. pp. 151 161, 2007.

---

[18] W. Truccolo, L. R. Hochberg, and J. P. Donoghue. "Collective dynamics in human and monkey sensorimotor cortex: predicting single neuron spikes," Nature Neuroscience, vol. 13. pp. 105- 111, 2010.
[19] Y. Roudi, B. Dunn, and J. Hertz, neuronal activity and functional connectivity in cell assemblies," Current Opinion in Neurobiology, vol. 32. pp- 38 44, 2015.
[20] S. ichi Amari, "Information geometry on hierarchy of probability distributions," IEEE Transac- tions on Information Theory, vol. 47, pp. 1701 1711, July 2001.
[21] E. Schneidman, M. J. Berry, R. Segev, and W. Bialek, "Weak pairwise correlations imply strongly correlated network states in a neural population," " Nature, vol. 440, pp. 1007-1012, 2006.
[22] J. W. Pillow, J. Schlens, L. Paninski, A. Sher, A. M. Litke, E. J. Chichilnisky, and E. P. Simoncelli, "Spatio-temporal correlations and visual signalling in a complete neuronal population." Nature, vol. 454, pp. 995-999, 2008.
[23] W. Truccolo, "From point process observations collective neural dynamies: Nonlinear hawkes process glms, low-dimensional dynamics and coarse graining," Journal of Physiology-Pure vol. 110, no. 4, pp. 336-347, 2016.
[24] A. Sheikhattar, S. Miran, J. Liu, J. B. Fritz, S. A. Shamma, P. O. Kanold, and B. Babadi, "Extracting neuronal functional network dynamics via adaptive Granger causality analysis," Proceedings of the National Academy of Sciences, vol. 1 115, no. 17, pp. E3869 E3878, 2018.
[25] S. Kim, D. Putrino, S. Ghosh, and E. N. Brown, "A granger causality measure for point process models of ensemble neural spiking activity," PLOS Computational Biology, vol. 7, P. e1001110. 03 2011.
[26] D. J. Daley and D. Vere-Jones, An Introduction to the Theory of Point Processes, vol. 1, New York, NY: Springer, 2nd ed., 2003.
[27] D. Ba, S. Temereanca, and E. N. Brown, "Algorithms for the analysis of ensemble neural spiking activity using simultaneous-event multivariate point-process models," Frontiers in Computational Neuroscience, vol. 8, no. 6, 2014.
[28] R. E. Kass, R. C. Kelly, and W.-L. Loh, "Assessment of synchrony in multiple neural spike trains using loglinear point process models," Ann. Appl. Stat., vol. 5, pp. 1262-1292, 06 2011.
[29] V. Solo, "Likelihood functions for multivariate point processes with coincidences," in Proceed- ings of IEEE Conference on Decision and Control, vol. 3, pp. 4245 4250, 2007.
[30] F. Kloosterman, S. P. Layton, Z. Chen, and M. A. Wilson, "Bayesian decoding using unsorted spikes in the rat hippocampus." " Journal of Neurophysiology, vol. 111, no. pp. 217-227, 2014. [31] X. Deng, D. F. Liu, K. Kay, L. M. Frank, and U. T. Eden, "Clusterless decoding of position from multiunit activity using a marked point process filter," " Neural Computation, vol. 27,
pp. 1438-1460, 07 2015.
[32] S. Mukherjee and B. Babadi, "A statistical approach to dynamic synchrony analysis of neuronal ensemble spiking," " in Proceedings of the 2019 Asilomar Conference on Signals, Systems, and Computers, Nov. 3-6, Pacific Grove, CA, 2019.
[33] S. S. Haykin, Adaptive Filter Theory. Upper Saddle River, NJ: Prentice Hall, 1996.
[34] E. Ganmor, R. Segev, and F. Schneidman, "Sparse low-order interaction network underlies highly correlated and learnable neural population code," Proceedings of the National Academy of Sciences, vol. 108, pp. 9679-96988 9684, June 2011.
[35] T. Zhang, "Sparse recovery with orthogonal matching pursuit under RIP." IEEE Transactions on Information Theory, vol. 57, pp.
September 2011.
[36] A. Kazemipour, M. Wu, and B. Babadi, "Robust estimation of self-exciting generalized linear models with application to neuronal modeling," IEEE Transactions on Signal Processing. vol. 65, pp. 1733-3748, July 2017.
[37] A. Litwin-Kumal K. D. Harris, R. Axel, H. Sompolinsky, and I.. Abbott, "Optimal degrees of synaptic connectivity." Neuron, vol. 93, no. 5, pp. 1153-1164.c7, 2017.
[38] D. Pfau, F. A. Pnevmatikakis, and I.. Paninski, "Robust learning of low-dimensional dynamics from large neural ensembles." in Advances Neural Information Processing Systems, vol. 26, Curran Associates, Inc., 2013.
12

---

[39] A. M. Thomson and A. P. Bannister, "Interlaminar connections in the neocortex," Cerebral Cortex, vol. 13. pp. 5-14, 01 2003.
[40] A. Sheikhattar, J. B. Fritz, S. Shamma, and B. Babadi, "Recursive sparse point process regres- sion with application to spectrotemporal receptive field plasticity analysis," IEEE Transactions on Signal Processing, vol. 64, pp. 2026-2039, April 2016.
[41] S. van de Geer, P. Buhlmann, Y. Ritov, and R. Dezeure, "On asymptotically optimal confidence regions and tests for nigh-dimensional models," Ann. Starist., vol. 42, no. 3, pp. 1166-1202, 2014.
[42] A. Javanmard and A. Montanari, "Confidence intervals and hypothesis testing for high- dimensional regression," Journal of Machine Learning Research, vol. 15, no. 82, pp. 2869-2909, 2014.
[43] C.-H. Zhang and S. S. Zhang, "Confidence intervals for low dimensional parameters in high dimensional linear models," Journal of the Royal Statistical Society. Series B (Statistical Methodology), vol. 76, no. 1, pp. 217-242, 2014.
[44] S. S. Wilks, "The large-sample distribution of the likelihood ratio for testing composite hypothe- ses," Ann. Math. Statist., vol. 9, no. 1, pp. 60-62, 1938.
[45] A. Wald, "Tests of statistical hypotheses concerning several parameters when the number of observations is large," Transactions of the American Mathematical Society, vol. 54, no. 3, pp. 426-482, 1943.
[46] P. Martinez-Cañada T. V. Ness, G. T. Einevoll, T. Fellin, and S. Panzeri, "Computation of the electroencephalogram (EEG) from network models of point neurons," PLOS Computational Biology vol. 17, no. 4, pp. 1 1-41, 2021.
[47] L. D. Lewis, V. S. Weiner, E. A. Mukamel, J. A. Donoghue, E. N. Eskandar, J. R. Madsen, W. Anderson, L. R. Hochberg, S. S. Cash, E. N. Brown, and P. Purdon, "Rapid fragmentation of neuronal networks at the onset of propofol-induced unconsciousness." Proceedings of the National Academy of Sciences, vol. 109, no. 49, pp. E3377-E3386, 2012.
13

---

Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's
contributions and scope? [Yes]
(b) Did you describe the limitations of your work? [Yes] We discussed the challenges of
tractability in modeling synchronous neuronal activity.
(c) Did you discuss any potential negative societal impacts of your work? [N/A] Our results provide a computational method for analyzing neuronal data. As such, no negative societal impacts are noted.
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes]
2. If you are including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [Yes] See Appendix C for the details of assumptions of our theoretical results.
(b) Did you include complete proofs all theoretical results? [Yes] See Appendix C for complete proofs of our theoretical results.
3. If you ran experiments..
(a) Did you include the code, data, and instructions needed to reproduce the main exper- imental results (either in the supplemental material or URL)? [Yes] Code and
simulated data are included as supplementary material; the experimental data acquired from human subjects
available upon request from the authors in [47]
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)? [Yes] Hyperparameters used to generate the results are all specified.
(c) Did you report errоr bars (c.g., with respect to the random seed after running exper- iments multiple times)? [Yes] Confidence intervals and type and II error rates are quantified and reported wherever applicable.
(d) Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, or cloud provider)? [N/A] The proposed algorithms are not computationally intensive. As such, no notable computing resources were used.
4. If you using existing assets (e.g., code, data, models) curating/releusing
assets... (a) If your work uses existing assets, did you cite the creators? [Yes] The experimental data was obtained with permission, and the original study cited.
(b) Did you mention the license of the assets? [N/A]
(c) Did you include any new assets either in the supplemental material or as URL? [Yes] Custom MATLAB code and simulated data are included as supplemental material.
(d) Did you discuss whether and how consent was obtained from people whose data you're using/curating? [Yes] See Section 5.2: human cortical spiking data were acquired with permission from the authors [47].
(e) Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content? [Yes] See Section 5.2: human cortical spiking data were acquired in fully anonymized format.
5. If you used crowdsourcing or conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, applicable? [N/A]
(b) Did you deseribe any potential participant risks, with links to Institutional Review Board (IRB) approvals, applicable? [N/A]
(c) Did you include the estimated hourly wage paid to participants and the total amount spent on participant compensation? [N/A]
14


Reviews:
[{"review_id": "Yyux9ooMg-b", "reviewer": "Reviewer_T3r8", "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. ", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "* 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. \n\n* 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. \n\n* Review of previous works that use marked point process tools for neural activity is missing. \n\n*  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. \n\n* 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.\n\n* 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. \n\n* 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. \n\n* 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.\n\n* 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. \n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "no negative societal impact", "review_text": "", "extra_scores": {"ethical_concerns": "no ethical concerns", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "6 hours"}, "extra_text": {}}, {"review_id": "TL8-ziGYCLA", "reviewer": "Reviewer_nqVp", "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. ", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "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. ", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "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_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5"}, "extra_text": {}}, {"review_id": "Kl0vp6u181C", "reviewer": "Reviewer_KgmN", "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.", "questions": "", "limitations": "", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "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.\n\nMy 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.\n\nFor 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.\n\nThe 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.\n\nThe 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?\n\nDespite 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.\n\n\nOther Comments\n--------------\n\n- On line 318 there is a typo, \"estimate\" should be \"estimated\"\n\n- In figures 2 + 3, please label the blue-to-red colorbar\n\n- On lines 167-169, are the \"u\" variables meant to be \"\\mu\"?", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "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_text": "", "extra_scores": {"ethical_concerns": "n/a", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}]

Author rebuttal:

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)