Dynamic Analysis of Higher-Order Coordination in Neuronal Assemblies via De-Sparsified Orthogonal Matching Pursuit
Year: 2021

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 pointprocess 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.

heading: 1 Introduction
text: Coordinated ensemble spiking has been observed in a variety of brain areas, prompting a range of hypotheses about its role in cognitive function. Studies have documented synchronous spiking at all levels of the mammalian visual pathway [1, 2, 3]. Coordinated neural activity has additionally been hypothesized to influence interareal communication and the flow of neural information [4, 5, 6, 7, 8], and been postulated to be mediated by oscillations in local field potentials [9, 10, 11]. The study of synchrony is also closely tied to memory [12, 13, 14].
The prevalence of coordinated spiking and its functional implications for a range of neural processes have motivated both model-free and model-based approaches to its characterization. An intuitive model-free metric is the pairwise correlations of spike trains smoothed by a 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 in or determine the significance of pairwise coherence without repeated trials.
35th Conference on Neural Information Processing Systems (NeurIPS 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) [18, 19]. Maximum entropy models describe the state of the neural population only in terms of its 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 to 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, as derived in [29] and expounded upon in [27]. MkPP representations of ensemble activity have also been utilized to analyze neuronal population coding in unsorted spiking data [30, 31]. A related approach models disjoint simultaneous spiking events as log-linear combinations of point process models, hence permitting an intuitive representation of excess or suppressed synchrony [10, 28].
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. Log-linear point process models address this shortcoming, 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 are required to capture dynamics in correlational structure and to evaluate the statistical significance of coordinated spiking.
We address these limitations by proposing an adaptive greedy filtering algorithm based on the discretized MkPP formulation in [27] to model dynamics in 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 utility in tracking dynamics in coordinated spiking with statistical confidence on simulated ensemble spiking. Applying our 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.

heading: 2 Preliminaries
text: In this section, we review the discrete-time marked process representation and two corresponding likelihood models. First, we present a summary of key notation used throughout the subsequent sections in Table 1.

heading: 2.1 Marked Process Representation of Ensemble Spiking
text: To characterize coordinated spiking, we utilize the discrete-time marked process (MkPP) representation of ensemble neuronal activity [27, 28]. For an ensemble of C neurons, the C-variate spiking process, binned with small bin size ∆, at time bin index t is denoted by nt := [n (1) t , n (2) t , . . . , n (C) t ] ′, 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 nt as multivariate Bernoulli observations. The spiking process nt is mapped to a C∗-variate process n∗t := [n ∗ t (1), n∗t (2), . . . , n∗t
(C∗)]′, which are the binned observations of a marked point process whose marks count the number of exactly one ofC∗ := 2C−1 disjoint non-zero spiking events; we refer to n∗t as the marked Bernoulli process, distinguishing it from the multivariate Bernoulli process nt. We define the mark space K := {1, . . . , C∗} [26]. Fig. 1 shows an example of mapping the activity of C = 3 neurons to C∗ = 7 marked processes. At each time tj such that ntj 6= 0, the sole non-zero element of n∗tj indicates the mark. We also define the binned ground process n(g)t that takes value 1 at each such tj and is zero otherwise [26]; the ground process indicates the occurrence of any spiking event and is represented by n(g)t := ∑C∗ m=1 n ∗ t (m).
The marked process representation is not unique, but can be defined in a convenient fashion: treating the components ofnt as the bits of aC-bit binary number, the mark indexed by the decimal equivalent of a particular realization of nt will corresponded to that realization. By the disjointness of the marked representation, the spiking process of the cth neuron can be recovered as the sum of all marked process whose index, in binary, takes value 1 at the cth bit. For instance, in Fig. 1, the spiking data of neuron 3 (in blue) is the sum of simultaneous 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. 1, bottom panel). To this end, we next describe two useful likelihood models for simultaneous spiking.

heading: 2.2 Two Likelihood Models of Simultaneous Spiking
text: In the discrete formulation, the conditional intensity functions (CIFs) of nt and n∗t are approximated by the probabilities of observing an event at time bin t given the ensemble’s spiking history. That is,
λ (c) t ∆ = P[n (c) t = 1|Ht], λ∗t (m)∆ = P[n∗t (m) = 1|Ht], (1)
for c = 1, . . . , C and m = 1, . . . , C∗. We can relate λ(c)t ∆ to λ ∗ t (m)∆ in the same manner as n(c)t to n∗t (m), and obtain the CIF of the ground process λ∗t (g)∆ = ∑C∗ m=1 λ ∗ t (m)∆.
The marked process permits a generative description of simultaneous spiking events: ensemble spiking events are characterized by the ground process, occurring with probability λ∗t
(g)∆; the event is then assigned to the mth mark (i.e. the mth simultaneous spiking outcome) with conditional probability λ ∗ t (m)∆
λ∗t (g)∆
. Thus, at time t the likelihood of ensemble event n∗t is given by:
p(n∗t ) = C∗∏ m=1
( λ∗t (m)∆
λ∗t (g)∆
)n∗t (m) ( λ∗t (g)∆ )n(g)t ( 1− λ∗t (g)∆ )1−n(g)t . (2)
The likelihood in (2) is used to form a multinomial generalized linear model (mGLM) with multinomial 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 contemporaneous spiking, permitting compact parameterization by µt = [µ (1) t , µ (2) t , . . . , µ (C∗) t ]
′. Defining the baseline firing parameter for the mth mark to be
µ (m) t := log
( λ∗t (m)∆
1− λ∗t (g)∆
) , m = 1, 2, · · · , C∗, (3)
or equivalently λ∗t (m)∆ = e
µ (m) t
1+ ∑C∗ j=1 e µ (j) t , the log-likelihood can be rewritten as a linear function of
n∗t , resembling the maximum entropy model [20, 21]:
log p(n∗t ) = µ ′ tn ∗ t − ψ(µt), where ψ(µt) := log
( 1 +
C∗∑ m=1 eµ (m) t
) . (4)
The second, more general version utilizes the ensemble history as covariates in the mGLM. Letting the covariate vector xt be the ensemble history up to some fixed lag at time t (augmented by a constant element of 1), the model is parameterized by ωt := [ ω (1)′ t ,ω (2)′ t , . . . ,ω (C∗)′ t
]′, where the parameters for themth mark ω(m)t := [ µ (m) t ,θ (m)′ t
]′ consists of an ensemble history-modulation vector θ
(m) t in addition to the baseline firing parameter. Thus, the log-likelihood in this case admits a similar
form to (4), by simply replacing µ(m)t with x′tω (m) t .

heading: 3 Adaptive Estimation of the History-Dependent mGLM
text: 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-independent model of [32] and extend it to the history-dependent mGLM.
We assume conditional independence across time bins and that the parameters ωt admit piece-wise constant dynamics and are constant over consecutive windows of length W . The ensemble history up to lag p defines the covariates as xt := [1, n (1) t−1, . . . , n (1) t−p, . . . , n (C) t−1, . . . , n (C) t−p]
′. The set of history covariate vectors at the ith window are denoted byXi = [x1+i(W−1), . . . ,xiW ]′. Since the mapping from n∗t to nt is injective, the influence of past spiking activity can be equivalently captured by defining history covariates in terms of either; however, using nt reduces the dimensionality of ωt and quantifies the influence of past spiking activity directly rather than through categorical variables. Let n ∗(m) i = [n ∗(m) 1+W (i−1), . . . , n ∗(m) iW ]
′ denote the sequence of outcomes of the mth mark in the ith window. The log-likelihood of the ith window is thus given by
`i(ωi) := C∗∑ m=1 n∗i (m)′Xiω (m) i − iW∑ j=1+(i−1)W log
( 1 +
C∗∑ m=1 exj ′ω (m) i
) . (5)
Motivated by the RLS objective function [33], a forgetting factor mechanism is utilized to combine the log-likelihoods up to the kth window, capturing the dynamics in each mark’s rates. For a forgetting factor 0 ≤ β < 1, the adaptively-weighted log-likelihood at window k is thus defined as:
`βk(ωk) := (1− β) k∑ i=1 βk−i`i(ωk). (6)
Parameter estimation can be performed by solving a sequence of maximum likelihood problems:
ω̂k := arg max ωk
`βk(ωk), k = 1, 2, · · · ,K. (7)
Two issues arise when considering large ensembles. First, the dimensionality of µk grows exponentially 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 K̄ = {m ∈ K : ∑ t n ∗ t
(m) > Nthr} for some pre-defined constant Nthr > 0, and treating the rates of the remaining marked processes as 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 a generalized Orthogonal Matching Pursuit (OMP) [35] [36]. The OMP iteratively identifies the parameter support set (the non-zero components) of fixed size (a hyperparameter for which we cross-validate) to encourage sparsity, thus capturing the inherent sparsity of network interactions based on past ensemble activity [37, 38, 39]. Moreover, the greedy estimation over a 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.1. The key element of AdOMP is efficient evaluation of the gradient ∇ω`βk(ωk) at the lth iterate ω̂(l),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 is crucial for the algorithm to operate in an online fashion. To this end, we utilize a recursive update rule to compute the gradient at the kth window, generalizing the adaptive filtering techniques employed in [40] for Bernoulli observations to a multivariate setting.

heading: 4 Statistical Inference of Higher-Order Coordination
text: 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 a 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 a nested hypothesis test; second, we generalize the de-sparsifying procedure for `1-regularized maximum-likelihood estimators established by [41] to the AdOMP; and third, we use the latter to establish a 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 rth-Order Coordinated Spiking
We characterize the significance of all rth-order simultaneous spiking events for the history-dependent mGLM. The corollary for the history-independent model is addressed in Appendix C.4. The significance of r-wise simultaneous spiking for r ≥ 2 is tested by considering the two alternatives:
H0 : r th-order simultaneous spikes occur as frequently as they would between independent units, given ensemble spiking history H1 : r
th-order simultaneous spikes occur at a significantly different rate than they would between independent units, given ensemble spiking history
(8)
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 as 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 rth-order interactions are chance occurrences by constraining the base rate parameters for each rth-order mark. For the mth mark, let u(m)t := xt′ω (m) k = µ (m) k + x̄ ′ tθ (m) k . We decompose the base rate parameter as µ (m) k = µ (m) 0,k + γ (m) k , where µ(m)0,k is rate under the null hypothesis and γ (m) k is analogous to the additional multiplicative factor in [28] that captures potential exogenous effects after conditioning on ensemble spiking history.
We thus estimate the reduced model ω̂(R)k := arg maxω(R)k ` β k(ω (R) k ), where the base rate parameters of rth-order events are constrained to those under the null hypothesis. That is, for each m ∈ Kr := {m ∈ K : ∑C c=1mc = r}, where mc is the c th least significant bit of m in binary, we fix µ(m)k to µ (m) 0,k and optimize the remaining parameters. To explicitly obtain the constraints, first recall that xt′ω(m)k is the log-odds of n∗(m)t = 1 versus n (g) t = 0 given ensemble spiking history. Under the assumption that the neurons are independent, the probabilities of each event is given, respectively, by
P[n∗t (m) = 1|Ht] = ∏ ca:mca=1 ( λ (ca) t ∆ ) ∏ cb:mcb=0 ( 1− λ(cb)t ∆ ) , and P[n(g)t = 0|Ht] = C∏ c=1 ( 1− λ(c)t ∆ ) . (9)
Taking the ratio evaluated at the full model estimate ω̂k, we obtain u (m) 0,t := ∑ c:mc=1 log ( λ̂ (c) t ∆
1−λ̂(c)t ∆
) .
Assuming the difference between u(m)t and u (m) 0,t is due only to exogenous factors, we estimate the corresponding term at the kth window as γ̂(m)k = 1 W ∑kW t=(k−1)W+1 ( u (m) t − u (m) 0,t ) . Thus, for the reduced model, we fix µ(m)k at µ̂ (m) k − γ̂ (m) k for m ∈ Kr. The hypotheses are then quantitatively stated as:
H0 : ωk = ω̂ (R) k , H1 : ωk 6= ω̂ (R) k . (10)
To control the possible abrupt variations of γ̂(m)k across windows, we apply Kalman forward/backward smoothing to the exogenous factor and use the smoothed values, γ̃(m)k , in lieu of γ̂ (m) k . The procedure is summarized in Algorithm 2 in Appendix A.3.

heading: 4.2 De-Sparsifying AdOMP Estimates
text: To test the hypotheses defined above, it is necessary to 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 a paucity of work to this end for greedily-estimated high-dimensional sparse models. In the closely-related problem of `1-regularized maximum-likelihood estimation, a set of elegant results [41, 42, 43] have established techniques to de-sparsify parameter estimates and construct confidence intervals. In particular, we extend the de-sparsification technique of [41], based on close inspection the Karush-Kuhn-Tucker conditions, in the greedy high-dimensional setting. In Appendix A.2, we derive the de-sparsified AdOMP parameters following s∗ iterations as
ŵk := ω̂(s∗),k − ( ∇2`βk(ω̂(s∗),k) )−1 ( ∇`βk(ω̂(s∗),k) ) . (11)
Next, the asymptotic normality of the de-sparsified AdOMP estimates is established. While a 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 `1-regularized maximum likelihood estimation. Hence, we adapt the treatment in [41] for AdOMP to establish the following result:
Theorem 1. Consider the maximization of the total data log-likelihood `βk(ωk) at the k th window, where the true parameter ωk ∈ Rd is (s, ξ)-compressible with ξ < 12 . Let ω 0 k be the maximum likelihood estimate and ω̂k be the AdOMP estimate after O(slog(s)) iterations. If conditions (C1)– (C6) are met, the de-sparsified AdOMP estimate ŵk satisfies√
1 + β 1− β ( ŵk − ω0k ) = Vk + oP(1) · 1,
where, as β → 1, Vk d−→ N (0, I−1k ) with I −1 k = −Σ −1 k 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 a detailed proof. Based on Theorem 1, 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 rth-Order Coordinated Spiking
Classical results on likelihood ratio tests between two nested hypotheses [44, 45] have established the use of the deviance difference as a common procedure. However, they are ill-suited in our setting due to the highly-dependent covariates and forgetting-factor mechanism in the data log-likelihood. These issues are addressed in a related context [24] for the inference of Granger-causal links by defining the adaptive de-biased deviance difference and characterizing its limiting distribution under presence and absence of Granger-causal links. We similarly utilize the adaptive de-biased deviance difference
D (r) k,β
( ω̂
(F ) k , ω̂ (R) k
) := ( 1 + β
1− β
)[ 2 ( `βk(ω̂ (F ) k )− ` β k(ω̂ (R) k ) ) − ( B (F ) k −B (R) k )] (12)
as the test statistic, where B(F )k and B (R) k are the respective biases of the full and reduced models. As we show in Appendix A.1, the full and reduced log-likelihoods can also be computed in an online fashion, in a similar manner as the gradients.
The limiting distributions of the adaptive de-biased deviance difference for the greedily-estimated 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-sparsified AdOMP estimates from Theorem 1:
Theorem 2. Let ω̂(F )k and ω̂ (R) k respectively be the full and reduced greedily-estimated mGLM parameters at window k, where ω̂(R)k assumes conditionally independent r th-order simultaneous spiking. Then, as β → 1,
i) if rth-order coordination matches independent rth-order interactions given ensemble spiking history, then D(r)k,β ( ω̂ (F ) k , ω̂ (R) k ) d−→ χ2(M (r)), i.e. chi-square, and ii) if rth-order coordination diverges from independent rth-order interactions given ensemble
spiking history, and assuming the base rate parameters of rth-order interactions scale at least as O (√
1−β 1+β ) , then D(r)k,β ( ω̂ (F ) k , ω̂ (R) k ) d−→ χ2(M (r), ν(r)k ), i.e. non-central chi-square, where ν(r)k is the non-centrality parameter at window k and the degrees of freedom M
(r) := |Kr| is the difference in the cardinalities of the full and reduced support sets.
A detailed proof is provided in Appendix C.3. In order to fully characterize the limiting distribution of D(r)k,β under H1, we must estimate the non-centrality parameter for each window. Assuming the parameter evolves smoothly in time, we use a state-space smoothing algorithm [24] to estimate it from the observed D(r)k,β values. This not only allows us to identify significant coordination, but to also quantify the degree of significance using Youden’s J-statistic
J (r) k := 1− α− Fχ2(M(d),ν̂(r)k ) ( F−1 χ2(M(d)) (1− α) )
(13)
for significance level α, where F (·) denotes the CDF. Values of Jk close to 1 imply that the rejection of the null is a stronger indication of coordination than for smaller values of Jk. Thus, the J-statistic characterizes the test in terms of both type I and type II errors. By convention, we take Jk = 0 when H0 is not rejected at the kth 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 rth-order coordination and using a signed J-statistic J
(r) k ·sgn (∑ m∈Kr γ̂ (m) k ) . The full procedure for identifying significant rth-order coordinated spiking
is summarized by Algorithm 3 in Appendix C.

heading: 5 Applications

heading: 5.1 Simulated Ensemble Spiking Data
text: We validate our proposed methods in a simulated example, performing complementary historyindependent and history-dependent analyses. Let the base rate parameter and exogenous effect for the history-independent model be denoted by µk and γk; and the same for the history-dependent model by µk,H and γk,H, with history-modulation parameter θk. Then, the reduced model constraints imply γk = γk,H+ x̄ ′ tθk. If the observed rate of higher-order events is equal to that of independent neurons, γk = 0; however, higher-order interactions may still be coordinated, i.e. γk,H = −x̄′tθk 6= 0. Conversely, the observed rate of higher-order events may differ from that of independent neurons, i.e., γk 6= 0. If γk,H = 0, 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 in Eq. (2) with an average rate of ∼ 0.1 spikes per bin. In the first and third simulated epochs, 4th-order spiking events were excited by amplifying the default history-modulation parameters. In the second epoch, the base rate parameter was increased to induce 3rd-order spiking events. These adjustments respectively reflected simultaneous spiking induced by ensemble history and by an unobserved process. Figure 2–A shows the simulated spiking activity, from which no obvious coordination is observable. The aggregate rth-order marks are visualized in Fig. 2–B, with apparent increased rates of 3rd- and 4th-order spiking events.
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
variability due to self-exciting dynamics while smaller ratios indicate regularity in spiking (i.e. globally coordinated spiking). Both measures are computed over non-overlapping windows of 250 samples to track dynamics. The third measure is the average difference between rth-order mark CIFs and probabilities of rth-order 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 rth-order coordination for r = 2, . . . , 5 using the history-independent model (W = 10; β = 0.975) reveals facilitated 3rd-order coordination during the second epoch, indicated by large positive values of the J-statistics (Fig. 2–C). Facilitated 4th-order coordination is detected during the first and third epochs. Ensemble spiking was also analyzed using the history-dependent model (W = 10; β = 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. 2–F) indicates Poisson spiking statistics rather than coordinated activity. The 3rd- and 4th-order mark CIF differences (Fig. 2–
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–G, insets). We further examine the effect of hyperparameters W and β in Appendix B.1, inspect the model goodness-of-fit in Appendix B.2, and provide a complementary simulation study with more complex dynamics in Appendix B.3.

heading: 5.2 Real Data Example: Anesthesia Data
text: We next present our analysis of human cortical neuronal assemblies during the transition into propofolinduced general anesthesia. The data were retrieved in a 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. 3–A, their ensemble spiking activity is shown, aligned to the loss of consciousness (LOC) at 0 s. Ensemble activity recovered after ∼250 s, when propofol was re-administered. The decomposition of ensemble spiking into rth-order events (Fig. 3–B) highlights lower rates of higher-order spiking events.
History-independent analysis (W = 10; β = 0.99) of rth-order spiking revealed high rates of 3rd-, 4th-, and 5th-order events during consciousness (Fig. 3–C). Additionally, 2nd-order coordination was suppressed during consciousness. However, history-dependent analysis (W = 10; β = 0.995) did not identify conditional coordination during consciousness (Fig. 3–D). 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 3rd-order 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 at ∼250 s. 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 are poorly reflected by the Pearson correlation and spiking regularity, both computed over windows of 250 samples (resp., Fig. 3–E and Fig. 3–F). Prior to LOC, pairwise correlations may have been high but spiking was not globally synchronized. Hence, spiking statistics seem 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.

heading: 6 Concluding Remarks
text: The proposed modeling and statistical inference algorithms constitute a novel approach to studying coordinated neuronal spiking. In contrast to previous model-based approaches, the proposed method is 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 to sparsity-regularized 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 or 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. Additionally, 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 to address current gaps in understanding the local mechanisms underlying the emergence of different brain states.

heading: Acknowledgements
text: This work is supported in 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 [47]) recorded from human cortical assemblies during anesthesia.

Review summaries:
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.