Adaptive modeling and inference of higher-order coordination in neuronal assemblies: a dynamic greedy estimation approach

Central in the study of population codes, coordinated ensemble spiking activity is widely observable in neural recordings with hypothesized roles in robust stimulus representation, interareal communication, and learning and memory formation. Model-free measures of synchrony characterize coherent pairwise activity but not higher-order interactions, a limitation 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 repeated trials 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 statistical inference framework to identify significant higher-order coordination. In the course of developing a precise statistical test, we 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 from human and rat cortical assemblies, our proposed methods provide new insights into the dynamics underlying localized population activity during transitions between brain states.


Introduction
Synchronous neuronal ensemble activity is central in the study of neural population codes.Coordinated ensemble spiking has been observed in a variety of brain areas, prompting a variety of hypotheses about its role in cognitive function.For instance, studies have documented synchronous spiking at all levels of the mammalian visual pathway [1][2][3].Synchronized thalamic population activity has also been widely observed, a phenomenon to which visual cortical neurons have been found sensitive, suggesting the importance of synchronized neuronal activity in thalamocortical communication [4,5].
Synchronized spiking has, more broadly, been hypothesized to influence inter-areal communication and the flow of neural information [6][7][8][9][10].The study of coordinated neural activity is also closely tied to oscillatory activity and memory.Synchronized hippocampal and hippocampal-cortical activity are thought to have significant roles in memory formation, working memory tasks, and encoding information for spatial navigation [11][12][13].Coordinated ensemble spiking has additionally been postulated to be mediated by oscillations in local field potentials [14][15][16].
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 quantifying spiking synchrony.Perhaps the most intuitive model-free metric is the pairwise correlations of spike trains smoothed by a Gaussian (or exponential) kernel [17,18].Other model-free measures include a range of spike train distance metrics that also perform pairwise comparisons [19].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.
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) [20,21].
Maximum entropy models describe the state of the neural population only in terms of its instantaneous correlational structure [22,23].Models are estimated to match observed firing rates and all pairwise (and potentially higher order) correlations simultaneously.The suitability of the maximum entropy model formulation for analyzing coordinated spiking has motivated several extensions.For instance, Bayesian state-space filtering algorithms have been developed to capture dynamics in the strength of higher-order spiking interactions [24,25].A stimulus-dependent maximum entropy model has also been proposed to address potential synchrony-modulating factors [26].
Extensions also include efforts to address the computational complexity associated with analyzing higher-order interactions amongst large neuronal assemblies [27,28].
Point process GLMs are a common alternative to maximum entropy models for ensemble spiking [29,30] that can characterize the influence of past population activity and other relevant covariates.Though useful in estimating functional connectivity [31,32], each neuron must be assumed conditional independent due to regularity conditions that prohibit simultaneous spiking events [33][34][35].This can be circumvented by using an equivalent marked point processes (MkPP) representation that explicitly models each disjoint simultaneous spiking event [34].MkPP representations of ensemble activity have also been utilized to analyze neuronal population coding in unsorted spiking data [36,37].A related approach models disjoint simultaneous spiking events as log-linear combinations of point process models that permits an intuitive representation of excess or suppressed synchrony [15,35].
The aforementioned statistical models enable the analysis of higher-order coordination in ensemble spiking, though each with their respective limitations.
Dynamics in the correlational structure of maximum entropy models may be tracked with state-space filtering algorithms and credible intervals can be constructed to assess October 16, 2023 2/27 the statistical significance of correlations; however, the influence of past population activity on the ensemble state is neglected, and assumptions on the relevance of higher-order interactions are typically imposed for tractability of model estimation.
Log-linear point process models can track dynamics in coherent spiking while incorporating the effects of temporal dynamics of population activity, and confidence intervals may be approximated for statistical inference; the necessity of a priori assumptions on the relevance of higher-order interactions still, however, remains.
Additionally, model-fitting in both approaches requires multiple repeated trials to capture dynamics in correlational structure, thus limiting their applicability to continuously recording spiking data.A discretized MkPP model is capable of capturing greater detail in the effects of past population activity on coordinated spiking, though these effects are assumed to be static.To our knowledge, the tractability of the MkPP model for ensemble spiking has not been addressed, and a corresponding statistical inference framework is lacking.
We address these gaps by proposing an adaptive greedy filtering algorithm based on the discretized MkPP formulation in [34] to model dynamics in higher-order spiking coordination in single-trial recordings while capturing the influence of past ensemble activity.Incorporating similar data-driven restrictions on modeled interactions as in [27], we also address the question of tractability of the discretized MkPP formulation.Furthermore, we build on recent theoretical results related to Adaptive Granger Causality (AGC) analysis [31] to provide a precise statistical framework to detect significant coordinated spiking activity of arbitrary order.We demonstrate our proposed method's utility in tracking dynamics in synchronous activity with statistical confidence on simulated ensemble spiking.Applying our method to continuous multi-electrode recordings of human cortical assemblies during anesthesia and to rat cortical assemblies during sleep provides novel insights into coordinated spiking dynamics that underlie transitions between brain states.

Materials and methods
We first summarize essential components of the proposed methods for inferring latent higher-order spiking coordination in order to contextualize our results.The primary contribution of this work is a framework for the dynamic and statistically precise inference of latent coordinated spiking in neuronal assemblies using their simultaneous spiking representation (Fig.

Discretized marked point process likelihood model
To characterize coordinated spiking, it is necessary to use an appropriate representation of neuronal ensemble spiking.Because multivariate point processes as defined in literature [33] do not permit simultaneous events at arbitrarily small time scales, point process models of ensemble spiking treat neurons as conditionally independent elements of a multivariate process.To avoid this assumption, Solo introduced a theoretical model for simultaneous spiking events based on marked point processes (MkPP) that explicitly model each possible ensemble spiking event [38].Due to the convenience of its disjoint October 16, 2023 3/27 CIF of the ground process, t , . . ., µ Base rate parameters of mark events , . . ., ω Model parameters of history-dependent model Log-odds of m th mark event vs. no spiking event Log-odds of m th mark event vs. no spiking event (reduced model) Exogenous factor for m th mark β Forgetting factor, 0 < β < 1 W Window length decomposition, MkPP representations have been used by Kass et al. in [35] and Ba et al. in [34] to model simultaneous spiking in neuronal assemblies.
Here, we utilize a discrete-time MkPP representation of ensemble neuronal activity.For an assembly of C neurons, the C-variate spiking process, binned with small bin size ∆, at time bin index t is denoted by n , 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 n t as multivariate Bernoulli observations.The spiking process n t is mapped to a C * -variate process n * which are the binned observations of a marked point process whose marks count the number of exactly one of C * := 2 C − 1 disjoint non-zero spiking events; we refer to n * t as the marked Bernoulli process, distinguishing it from the multivariate Bernoulli October 16, 2023 4/27 process n t .We define the mark space K := {1, . . ., C * } [33].The example in Fig. 1 shows the activity of C = 3 neurons mapped to a marked process with C * = 7 marks.
At each time t j such that n tj = 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 t j and is zero otherwise [33]; the ground process indicates the occurrence of any spiking event and is represented by n The marked process representation is not uniquely defined, but can be done so in a convenient fashion: treating the components of n t as the bits of a C-bit binary number, the mark indexed by the decimal equivalent of a particular realization of n t will corresponded to that realization.By the disjointness of the marked representation, the spiking process of the c th neuron can be recovered as the sum of all marked process whose index, in binary, takes value 1 at the c th bit.For instance, in Fig. 1, the spiking activity of neuron 3 (in blue) is the sum of simultaneous spiking event processes 4-7.
The conditional intensity functions (CIFs) of n t and n * t are approximated by the probabilities of observing an event at time bin t given ensemble spiking history.That is, for c = 1, . . ., C and m = 1, . . ., C * .We can relate λ 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 m th mark (i.e. the m th simultaneous spiking outcome) with conditional probability Thus, at time t the likelihood of ensemble event n * t is given by: The likelihood in Eq. ( 2) is used to form a multinomial generalized linear model (mGLM) with multinomial logistic link function of which we consider two versions.
The first, more general version utilizes the ensemble history as covariates in the mGLM.Letting the covariate vector x t be the ensemble history up to some fixed lag at time t (augmented by a constant element of 1), the model is parameterized by , . . ., ω , where the parameters for the m th mark .The CIF of the m th mark is hence defined as: equivalently, the log-odds of the m th mark occurring versus no spiking event occurring The second version of the mGLM makes the simplifying assumption that there is no history dependence.The resulting model depends only on contemporaneous spiking, permitting compact parameterization by the baseline firing parameters t , . . ., µ ] .In the history-independent model, the log-odds of the m th mark occurring versus no spiking event occurring are: October 16, 2023 5/27 The log-likelihood can be expressed in a fashion resembling the maximum entropy model [22,23].For the history-independent case, the log-likelihood of n * t is log p(n * t ) = µ t n * t − ψ(µ t ), where ψ(µ t ) := log 1 + while the history-dependent formulation admits a similar expression by simply replacing Adaptive estimation of marked point process models In contrast with conventional mGLM models, here the parameters are allowed to change in time.Hence, we use adaptive algorithms to capture the dynamics of the history-dependent and history-independent models.However, analyzing large neuronal assemblies raises the issue of tractability since the number of parameters to be estimated scales exponentially with C. Since it is likely that some marks will not contain any events, we employ a thresholding rule similar to [27], considering only "reliable interactions", i.e. the subset of the mark space K = {m ∈ K : t n * t (m) > N thr } for some pre-defined constant N thr > 0, and treating the rates of the remaining marked processes as negligible due to their infrequency.For generality and clarity in notation, subsequent discussions are in terms of the full mark space K.
The parameters of the history-dependent discretized MkPP model are obtained by solving a sequence of maximum likelihood problems, formulated as follows.We assume that the parameters ω t admit piece-wise constant dynamics and are constant over consecutive windows of length W ; additionally, observations are assumed conditionally independent across time bins.The ensemble history up to lag p defines the covariates as t−1 , . . ., n denote the sequence of outcomes of the m th mark in the i th window.The log-likelihood of the i th window is thus given by Motivated by the RLS objective function [39], a forgetting factor mechanism is utilized to combine the log-likelihoods up to the k th 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: Parameter estimation is hence performed by solving the sequence of maximum likelihood problems: To efficiently solve the sequence of problems in Eq. ( 8) in an online fashion, we use the Adaptive Orthogonal Matching Pursuit (AdOMP) [40], an adaptive version of the October 16, 2023 6/27 Orthogonal Matching Pursuit (OMP) [41] [42].The AdOMP and OMP both iteratively identify 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 [43][44][45].
However, AdOMP permits the support set to change between windows.Moreover, greedy estimation over a sparse subset of parameters mitigates the intractability of the estimation problem for large neuronal assemblies where regularization-based constraints still require optimization over all parameters.The AdOMP relies on efficient evaluation of the gradient ∇ ω β k (ω k ) at the l th 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 k th window, generalizing the adaptive filtering techniques employed in [46] for Bernoulli observations to a multivariate setting.A complete derivation of the recursive update rule and the utilization of the gradient in AdOMP is provided in supporting information (S1 Appendix).
The sequence of maximum likelihood problems that must be solved to obtain the history-independent model takes a similar form as in Eqs. ( 6) -( 8) under the same assumption of piece-wise constant dynamics of µ t .Defining n * i := 1 W iW j=(i−1)W +1 n * j , the equivalent of the i th window log-likelihood in Eq. ( 6) is given by . Hence, the adaptively-weighted log-likelihood at the k th window is: The sequence of maximum likelihood estimates μk = arg max are obtained by gradient descent.The gradient of the history-independent log-likelihood in Eq. ( 9), while involving a weighted summation of the terms n * i for which a simple recursion is derived, can be computed directly.The procedure for computing the maximum-likelihood estimates of the history-independent model is detailed in supporting information (S1 Appendix).

Statistical inference of higher-order coordination
Coordinated spiking can indicate relationships between components of a neuronal assembly and, potentially, effects of unobserved processes.However, independent neurons can also spike concurrently by chance, necessitating statistical inference to distinguish between excessive (or suppressed) and chance simultaneous spiking.To this end, we quantify the two alternatives as nested hypotheses and prove that an adaptive de-biased deviance test used for identifying significant Granger-causal influences [31] is applicable to our setting, thus establishing a precise statistical inference framework.
Here, we focused on characterizing the significance of r th -order simultaneous spiking and have formulated the hypothesis test accordingly; however, similarly constructed null hypotheses can be used to test the significance of any set of simultaneous spiking events using the same inference procedure (see supporting information S1 Appendix).For cogency, we focus on the statistical inference procedure for history-dependent MkPP models; differences for the history-independent model are addressed in supporting October 16, 2023 7/27 information (S1 Appendix).However, the complementary nature of the two models are summarized here.Theoretical results pertaining to the precise inference framework are summarized here, and comprehensively described in supporting information (S1 Appendix).
Formulating nested hypotheses to test for r th -order coordination The significance of r-wise simultaneous spiking for r ≥ 2 is tested by considering the two alternatives: H 0 : r th -order simultaneous spikes occur as frequently as they would between independent units, given ensemble spiking history H 1 : r th -order simultaneous spikes occur at a significantly different rate than they would between independent units, given ensemble spiking history (11) A similar formulation was used in [35] to determine whether one mark occurs at a significantly different rate than expected.The likelihood of the mark was modeled as the product of marginal likelihoods and an additional multiplicative factor.Noting that the additional factor takes value 1 if the neurons are truly independent, the null hypothesis was quantified accordingly.Instead, we estimate a reduced model that assumes r th -order interactions are chance occurrences by constraining the base rate parameters for each r th -order mark.For the m th mark, let We decompose the base rate parameter as 0,k is rate under the null hypothesis and γ (m) k is analogous to the additional multiplicative factor in [35] that captures potential exogenous effects after conditioning on ensemble spiking history.
The reduced model is estimated by solving the sequence of maximum likelihood problems ω(R) k ) with the base rate parameters of r th -order events constrained to their values under the null hypothesis.Let where m c is the c th least significant bit of m in binary.Then, to obtain the reduced model for any m ∈ K r , we fix µ (m) k to µ (m) 0,k and optimize the remaining parameters.To explicitly obtain the constraints, first recall that versus n (g) t = 0 given ensemble spiking history.Under the assumption that the neurons are independent, the probabilities of each event are given, respectively, by Evaluating the log-ratio of these two probabilities at the full model estimate ωk , we obtain that October 16, 2023 8/27 Attributing the difference between u (m) t and u (m) 0,t only to exogenous factors, we estimate the exogenous factors at the k th window as γ(m) Thus, for the reduced model, the value of µ (m) The hypotheses are then quantitatively stated as: 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 given in more details in supporting information (S1 Appendix).

Precise statistical inference using deviance differences
The use of the deviance difference test statistic has been established in classical statistical methodology [47,48] as a common procedure for likelihood ratio tests between two nested hypotheses.However, such a test is ill-suited to our setting due to the highly-dependent covariates and forgetting-factor mechanism in the data log-likelihood.These issues were addressed in a related context [31] 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, as the test statistic, where are the respective biases of the full and reduced models.The full and reduced model log-likelihoods can be efficiently computed online in a similar manner as the gradients (see supporting information S1 Appendix).We precisely characterize the limiting behavior of the deviance difference in Eq. ( 19) under both the null and alternative hypotheses, showing that as β → 1: i) if r th -order coordination matches independent r th -order interactions given ensemble spiking history, then ), i.e. chi-square, and ii) if r th -order coordination diverges from independent r th -order interactions given ensemble spiking history, and assuming the base rate parameters of r th -order interactions scale at least as In order to fully characterize the limiting distribution of D k,β values.Hence, in addition to identifying significant coordination, we also quantify the degree of significance using Youden's J-statistic for significance level α, where F (•) denotes the CDF.Values of J k .Thus, the J-statistic characterizes the test in terms of both type I and type II errors.By convention, we take J (r) k = 0 when H 0 is not rejected at the k 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 r th -order coordination and using a signed J-statistic The full procedure for identifying significant r th -order coordinated spiking is summarized by Algorithm 1.
Algorithm 1 Dynamic inference of r th -order spiking coordination for m ∈ K r do 6: Evaluate {u end for

9:
Estimate ω(R) k , and end if 14: end for 15: Estimate {ν In establishing the limiting behavior of D k,β , we also proved a result of independent interest; namely, we generalized asymptotic properties of de-sparsified 1 -regularized estimates [49] to de-sparsified greedy estimates by showing that the de-sparsified AdOMP estimate behaves asymptotically like the maximum likelihood estimate.
Crucially, this allows for the construction of confidence intervals around greedily estimated parameters thus enabling precise statistical inference.

Complementary characterizations of higher-order coordination by history-independent and history-dependent models
The history-independent model is a special case of the history-dependent model; hence, the statistical inference procedure summarized by Algorithm 1 can be appropriately October 16, 2023 10/27 modified.This specialization is expounded in supporting information (S1 Appendix).
However, the interpretation of statistically significant results obtained using history-independent analysis is distinct from but complementary to the interpretation of statistically significant results obtained using history-dependent model.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 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 = 0. Conversely, the observed rate of higher-order events may differ from that of independent neurons, i.e., γ k = 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.Thus, history-independent analysis reveals if the observed rate of simultaneous spiking events deviates from the expected rate in a group in independent neurons, while history-dependent analysis reveals if the observed rate of simultaneous events cannot be attributed to endogenous network effects.

Results
The proposed methods for analyzing higher-order spiking coordination were validated through simulations.First, we empirically verified our theoretical results characterizing the limiting behavior of the adaptive debiased deviance difference on simulated spiking.Next, we demonstrated the utility of the analyses on simulated ensemble spiking data, showing that the ground-truth latent dynamics in higher-order coordinated spiking can be recovered with statistical confidence and with greater accuracy than existing single-trial metrics.Then, the proposed methods were applied to continuous multi-electrode recordings of human cortical assemblies during anesthesia and to rat cortical assemblies during sleep in order to infer latent dynamics in coordinated spiking during transitions between brain states.

Empirical Validation of the Limiting Behavior of Deviance Differences
We first validated the proposed statistical test for r th -order coordinated spiking by empirically verifying the limiting distributions of the adaptive debiased deviance difference derived under the null and alternative hypotheses.We simulated and analyzed 50 realizations of a 5-neuron ensemble spiking process; the spiking activity was generated by a marked Bernoulli process, as described by Eq. ( 2).Each realization was 4000 samples long with dynamics in 3 rd -order spiking coordination.Namely, a step function was used to exogenously facilitate 3 rd -order spiking during the second half of each realization.
The simulated spiking data were analyzed by applying Algorithm 1 to each realization; restricted models were only computed for 3 rd -order spiking coordination.A threshold of N thr = 1 was used to pruned marked events that occurred no more than once per realization; this excluded 5 th order events from the estimated model.The window size over which parameters were assumed constant was set to W = 10 in order to enable stable estimation at each window while still allowing for fast changes.The forgetting factor was set to β = 0.99.For the purpose of validating the limiting distributions, a forgetting factor close to 1 was desirable; however, β = 0.99 is also a practical choice for the forgetting factor that serves to illustrate the utility of the limiting distributions when analyzing physiological data.Indeed, the hyperparameter October 16, 2023 11/27 choices (W, β) = (10, 0.99) were used in the applications to physiological data presented later.Hyperparameter choices for analysis are discussed in the context of two additional simulations.
The limiting distribution under the null hypothesis was validated by compiling adaptive de-biased deviance differences that were computed during the first half of each realization and corresponded to small estimated non-centrality parameter values.Their distribution, depicted by the blue histogram in Fig. 2, closely matched the theoretical distribution of deviance differences under the null hypothesis, i.e. a χ 2 distribution with M = 10 degrees of freedom.

Simulated Ensemble Spiking: Example 1
We next demonstrated the utility of the proposed methods in application to two sets of simulated ensemble spiking data with dynamic latent higher-order coordination.In the first example, spiking activity of 5 neurons was simulated for 16000 samples.Spiking activity (Fig. 3A) included 3 rd -order events exogenously facilitated and suppressed in alternation by a square wave and 4 th -order events induced through endogenous effects of ensemble spiking history throughout the simulated duration.The latent spiking coordination was evident when visualizing the sums of all r th -order events in Fig. 3B.
Four epochs of the simulation were defined by the periods of 3 rd -order facilitation and suppression, shown in Fig. 3C, and are indicated by vertical dashed lines common across all panels.
History-independent and history-dependent analyses of higher-order coordination were applied to the simulated spiking data using the same hyperparameters.A conservative threshold of N thr = 1 was chosen to prune marked events that occurred unreliably over the simulated duration.The window size over which parameters were assumed constant was set to W = 10 in order to enable stable estimation at each window while still allowing for fast changes.Fixing W , several candidate values for the forgetting factor were considered to obtain the most appropriate effective integration window, N eff = O( W 1−β ) [46], and was set to β = 0.975.The effect of varying the forgetting factor are presented in supporting information (S1 Appendix).In summary, for simulations the best choice of β corresponded to N eff ≈ τ 10 , where τ denotes the duration of the shortest latent state.For the present simulation, τ = 4000, the half-period of the square wave that defined alternating states of 3 rd -order facilitation and suppression; hence, N eff = 400.Statistical tests were performed at level α = 0.01.
October 16, 2023 13/27 History-independent inference of higher-order synchrony (Fig. 3D) accurately characterized the periods of facilitated 3 rd -order coordinated spiking and correctly assessed the rates of 4 th -order spiking events to be significantly higher than expected amongst independent neurons, indicated by J-statistic values close to +1.However, 3 rd -order suppression was not detected as a statistically significant, likely reflecting that the expected rate of 3 rd -order was low to begin with and therefore difficult to distinguish.In complement, history-dependent inference of higher-order coordination (Fig. 3E) correctly attributed 4 th -order spiking, which occurred at a statistically significantly high rate, to endogenous network effects captured by ensemble spiking history regressors while detecting that 3 rd -order spiking was exogenously facilitated.
For comparison, three single-trial measures of coordinated spiking were utilized.The first is the average Pearson correlation between smoothed spiking responses.The second is the spiking regularity, quantified by the average coefficient of variation (ratio of the standard deviation to the mean inter-spike interval) [50].A ratio close to 1 indicates Poisson spiking 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 were computed over non-overlapping windows of 200 samples to track dynamics, which while not identical are of a similar order of magnitude as N eff .The third measure is the average difference between r th -order mark CIFs and probabilities of r th -order independent interactions, generalizing the measure employed in [34] to higher-order simultaneous spiking.Other existing model-based analyses require multiple trial repetitions and were thus unsuited to the single-trial simulation setting.
In application to simulated ensemble spiking data, the three control measures were unable to capture the latent dynamics in spiking coordination.Significant pairwise correlations (Fig. 3F) were detected throughout the simulated duration, indicating only that several pairs of neurons were spiking concurrently, but were insensitive to the changes between facilitative and suppressive states of the exogenous process.Similarly, the spiking variability measure (Fig. 3G) indicated Poisson-like spiking statistics throughout the simulation without reflecting any latent dynamics.The average mark CIF differences of 3 rd -order events (Fig. 3H) weakly reflected the dynamics of the exogenous process, but closer inspection (Fig. 3H inset) reveals the oscillatory nature and wide confidence intervals of this sample-by-sample measure which pose challenges in interpreting the analysis.

Simulated Ensemble Spiking: Example 2
The second simulated example utilized an autoregressive process instead of a square wave to induce exogenous 3 rd -order coordinated spiking in a 5-neuron assembly.
Ensemble spiking was simulated for 12000 samples (Fig. 4A) with 3 rd -order events exogenously induced by one realization of an autoregressive process.Additionally, 4 th -order events were induced through endogenous effects for the first and last 4000-samples periods of the simulated duration, but occurred with chance-level probability otherwise.The sums of all r th -order events (Fig. 4B) reflected the latent spiking coordination.Coordinated 3 rd -order spiking was most evidently facilitated during an interval when the exogenous variable had value greater than 2 (Fig. 4C); the interval is indicated by vertical dashed lines common across all panels.
Both history-independent and history-dependent analyses were applied to the second simulated spiking data set using the same hyperparameters.The mark space was again pruned to include only events that occurred more than N thr = 1 times; parameters were assumed constant over windows of W = 10 samples; and statistical tests were performed at level α = 0.01.The forgetting factor was varied as previously described (supporting information S1 Appendix) and set to β = 0.95, which corresponded to an effective  integration window N eff = 200.We noted that the exogenous autoregressive process 451 October 16, 2023 15/27 most persistently facilitated 3 rd -order coordinated spiking for a duration of ∼ 4000 samples; and within that duration, two subintervals of ∼ 2000 samples separated at time index ∼ 8000 can be discerned upon visual inspection (Fig. 4C).Hence, taking τ = 2000, the choice of β = 0.95 satisfies the criterion that N eff ≈ τ 10 .
The history-independent analysis of higher-order coordination correctly detected statistically significant 3 rd -order coordination during the interval in which the exogenous variable was greater than 2 and 4 th -order coordination when they were induced by ensemble spiking history (Fig. 4D).The history-dependent analysis also correctly identified the exogenous facilitation of 3 rd -order events while attributing 4 th -order coordination to endogenous effects (Fig. 4E).For comparison, the average Pearson correlation, average spiking regularity, and average mark CIF differences were computed in identical fashion as for the first simulation.The average Pearson correlation exhibited indicated significant pairwise correlations concurrently with both exogenously induced 3 rd -order events and endogenously induced 4 th -order events (Fig. 4F); however, these two facets of latent higher-order coordination could not be disambiguated.In contrast, the average spiking regularity did not exhibit any dynamics; Poisson-like spiking statistics were indicated throughout the simulated duration (Fig. 4G).The average mark CIF differences for 3 rd -and 4 th -order marks both weakly indicated the latent higher-order coordination (Fig. 4H).In addition to previously issues concerning large confidence intervals and oscillatory nature, deviant average mark CIF differences for 3 rdand 4 th -order events appear identical despite being induced in different manners.This illustrates that the average mark CIF differences only indicate when rates of r th -order events deviate from the expected rate and cannot further address latent structure.

Ensemble Spiking in Anesthetized Humans
In the first application to recorded spiking data, we analyzed microelectrode recordings of human cortical neurons during the transition into propofol-induced general anesthesia.Commonly used in surgical procedures, general anesthesia is a drug-induced neurophysiological state of sedation and unconsciousness.In a study of the transition into unconsciousness, simultaneous recordings of single-neurons, LFP, and electrocorticograms were acquired to analyze changes to neural activity and functional connectivity over multiple spatial scales (full details of the experimental procedure are described in [51]).To complement previous analysis of pairwise spiking correlations, we employed the proposed methods for characterizing higher-order coordinated spiking.
Spiking data from the microelectrode recordings of one subject were analyzed, focusing specifically on the 8 neurons with the highest average firing rate over the 1000 second recording.Multi-unit spike recordings were originally oversampled at 1kHz, but downsampled by a factor of 50 to reduce computational complexity.Hence, the definition of simultaneous spiking in this analysis was taken to be the occurrence of spiking events across multiple neurons within at least 50ms of each other.Ensemble spiking activity is shown in Fig. 5A, aligned to the loss of consciousness (LOC) at 0s when propofol was first administered; the effect was evident from the rapid decrease in spiking.Spiking activity recovered and after 250s propofol was administered again.In order to analyze higher-order coordination with the proposed methods, the mark space of C * = 2 8 − 1 possible simultaneous spiking events, K, was pruned to the set of reliable interactions K that occurred more than N thr = 15 times; that is, simultaneous spiking events with average rates less than 0.015Hz were treated as negligible.The cardinality of the set of reliable interactions, defined with a conservative threshold, was The sums of all r th -order events (Fig. 5B) show that up to 4 th -order coordinated spiking occurred reliably, though less frequently after LOC.
History-independent and history-dependent analyses were performed using the same hyperparameters.The window over which parameters were assumed constant was set to October 16, 2023 16/27 W = 10.The forgetting factor, β = 0.99, was selected so that N eff ≈ τ 5 ; here, we used 503 τ = 5000, the approximate number of samples between the two administrations of 504 propofol.For applications to recorded data, choosing β such that N eff ≈ τ 10 (as was 505 done for simulated data) yielded inferred higher-order coordination that was statistically 506 weak (as quantified by J-statistics) and transient, resembling simulated examples where 507 N eff was mismatched to the duration of latent states (supporting information S1 508 Appendix).We speculate that a shorter effective integration was appropriate in 509 simulations because the assemblies were comprised of 5 neurons with similar firing rates, 510 which facilitated tracking latent dynamics.This contrasts with the variability in firing 511 rates that can be observed in Fig. 5A.Finally, statistical inference was performed at 512 level α = 0.01.

513
Applying history-independent higher-order coordination analysis revealed sustained 514 significantly high rates of 2 nd -, 3 rd -, and 4 th -order events prior to LOC (Fig. 5C).

515
Moreover, conditioning on ensemble spiking history indicated that 3 rd -and 4 th -order 516 events were exogenously facilitated, while 2 nd -order events were exogenously suppressed 517 (Fig. 5D).This latent structure was disrupted immediately following LOC; as spiking 518 activity diminished, no higher-order coordination was detected.However, as spiking 519 activity recovered, 3 rd -and 4 th -order events (but not 2 nd -order events) occurred at 520 significantly high rates.As the second administration of propofol again diminished 521 spiking activity, the rate of 4 th -order events became insignificantly different from the 522 October 16, 2023 17/27 expect rate amongst independent neurons and did not recover.However, transient 3 rd -order spiking after the second administration continued that the history-independent analysis detected as statistically significant.Third-order spiking was sustained at a significantly high rate once ensemble spiking activity recovered.Notably, none of the higher-order coordinated spiking after LOC was exogenously induced.
The dynamics in higher-order spiking coordination described by the proposed methods were poorly reflected by the average Pearson correlation and average spiking regularity.Both measures were computed over windows of 200 samples in order to track changes during the transition into anesthesia.Average correlations seemed to be significantly greater than zero for longer intervals after LOC than during consciousness, but trends in the average correlation were difficult to distinguish (Fig. 5E).The average spiking regularity measure indicated Poisson-like spiking statistics throughout, contrasting the dynamics of higher-order coordination described by the proposed analyses.Average spiking regularity was ill-suited to analyzing dynamics after LOC due to the reduced spiking activity; this was reflected by abrupt changes and wide confidence intervals (Fig. 5E).
In summary, history-independent and history-dependent analyses of ensemble spiking during the transition into anesthesia revealed the rapid onset of differences in latent higher-order coordination that distinguished consciousness from anesthesia.
Specifically, comparisons between the history-independent and history-dependent results suggest that exogenous influences on the higher-order interactions of small neuronal assemblies during consciousness are disrupted during anesthesia.These results are corroborated by previous analyses of these data [51] that indicated a rapid state change in which local network interactions were preserve but spatially distant network interactions were disrupted during anesthesia.Previous studies have shown that propofol acts by enhancing GABAergic circuits whose recurrent dynamics contribute to inducing synchronized slow-wave oscillatory activity [51][52][53][54][55]. Ensemble spiking history regressors likely accounted for these recurrent dynamics in the history-dependent model so that no exogenous effects were detected.

Ensemble Spiking in Sleeping Rats
We additionally analyzed ensemble spiking data recorded from rat cortical neurons during sleep.Sleep consists of cyclical transitions between brain states that maintain homeostatic neural activity distinct from waking states; however, both the purpose and mechanisms of these transitions remain unclear.We analyzed large-scale spike recordings from frontal and motor cortices during sleep obtained to study the effects of different sleep stages on the firing rate dynamics of putatively excitatory (pE) pyramidal neurons and putatively inhibitory (pI) interneurons [56,57].By examining neuronal activity recorded during several instances of rapid eye movement (REM), non-REM (nREM), and microarousal states over multiple sleep cycles, the study sought to address homeostatic effects of sleep.Instead, we sought to use the proposed analyses of higher-order spiking coordination to study the dynamics during transitions into sleep and between REM and nREM states in one sleep cycle.
We analyzed spiking data during one 182s long sleep cycle from one animal in which at least 10 pE and pI neurons were identified, selecting the 10 neurons of each class with the highest average firing rate.Recordings were originally oversampled at 20kHz, but downsampled to 200Hz to reduce computational complexity.Simultaneous spiking in this analysis hence equated to the occurrence of spiking events across multiple neurons within at least 5ms of each other.Ensemble spiking activity of pE and pI neurons were analyzed separately; the activity of each population during the sleep cycle is shown in For tractable analysis, the mark spaces of C * = 2 10 − 1 possible simultaneous spiking events, K, of both populations were pruned to the set of reliable interactions K that occurred more than N thr = 10 times; that is, simultaneous spiking events with average rates less than 0.055Hz were treated as negligible.The cardinality of the set of reliable interactions amongst pE neurons was | KpE | ≈ 0.017 • C * and amongst pI neurons was | KpI | ≈ 0.058 • C * .The sums of all r th -order events (Fig. 6B,H) show that up to 3 rd -order coordinated spiking occurred reliably amongst pI neurons while only up to 2 nd -order interactions occurred reliably amongst pE neurons.The same effective integration windows were used for history-independent and history-dependent analyses of both neuronal populations; with W = 10, the forgetting factor β = 0.99 so that N eff = τ 5 , where τ ≈ 5000 was the duration of the second nREM interval.Statistical inference was performed at level α = 0.01.
Applying the history-independent and history-dependent analyses of higher-order coordination to the ensemble spiking of pE neurons in concert identified intervals of significantly higher rates of 2 nd -order events that could be attributed to effects of ensemble spiking history (Fig. 6C-D).Most of the detected intervals were not sustained during either REM or nREM sleep; rather, they were aligned to the transitions between states.However, pI neurons exhibited more structured higher-order coordination.History-independent analysis of pI neurons revealed that 2 nd -order events had significantly higher rates during two intervals; the first was during the wake-period, and the second started at the end of the first wake-period and ending at the transition from REM to nREM sleep (Fig. 6I).While during the first of these intervals the October 16, 2023 19/27 facilitation of 2 nd -order events could largely be attributed to ensemble history effects, there was a shift in the exogenous effects on 2 nd -order during the second interval (Fig. 6J).That is, after the transition from the wake state to REM sleep, the exogenous suppression of 2 nd -order events gradually shifted to exogenous facilitation by the middle of the REM period that persisted into the nREM period.Exogenous 2 nd -order coordination was no longer detected strongly after the first half of the nREM period, but exogenous suppression emerged again in the second wake-period.
In addition to dynamics in 2 nd -order coordination, pI neurons also exhibited significant 3 rd -order coordination.The rate of 3 rd -order events was significantly high during the first wake period and REM sleep; though significantly higher at the start of nREM sleep, only statistically weak and transiently high rates were detected during the middle and end of nREM sleep.However, in the second wake period, the rate of 3 rd -order coordinated events again became significantly high (Fig. 6I).Notably, the high rate of 3 rd -order events during REM was distinctive because it was exogenously facilitated, whereas 3 rd -order events during other periods occurred at significantly high rates because of endogenous effects (Fig. 6J).
In contrast to the proposed analyses, neither the average Pearson correlation nor average spiking regularity, computed over windows of 200 samples, reflected similar latent dynamics of higher-order coordination amongst pE or pI neurons.For pE neurons, pairwise correlations were close to 0 for much of the sleep cycle with the exception of a few windows (Fig. 6E).However, the spiking regularity was significantly less than 1 for much of the sleep cycle (Fig. 6F); the implication of globally coordinated ensemble spiking is at odds with the absence of reliably occurring higher-order spiking events amongst pE neurons.For pI neurons, the average correlation was significantly higher than 0 during the first and second wake periods, mirroring the significantly high rates of higher-order events during these intervals; however, excepting a few windows, the average correlation did not significantly differ from 0 during REM sleep (Fig. 6K), presenting an inconsistency with the rates of higher-order events.Meanwhile, the average spiking regularity did not differ significantly from 1 for most of the sleep cycle, indicating Poisson-like spiking activity (Fig. 6L); this contrasts sharply from the reliable occurrence of higher-order events.
Summarily, applying the history-independent and history-dependent analyses of higher-order spiking coordination revealed distinctive latent dynamics amongst pE and pI neurons during the same sleep cycle.Intervals of significant 2 nd -order spiking coordination amongst pE neurons were attributable to the effects of ensemble spiking history and occurred around the transitions between arousal states rather than being sustained during the arousal states, possibly relating to a hypothesis that transition periods are themselves distinct states [58].In contrast, 2 nd -and 3 rd -order spiking events amongst pI neurons were detected to be exogenously coordinated, especially during sleep states.The observed changes in higher-order coordination of pI neurons during REM sleep are consistent with previous results that have shown excitation of pI neuronal activity and coordination during REM sleep [59,60].Additionally, the detected exogenous influences on pI neurons may be explained by studies that have indicated signatures of REM sleep can be found in hippocampal neurons prior to cortical neurons [61,62].

Discussion and concluding remarks
Relations to other models of coordinated spiking activity The proposed algorithms integrate some notable functionalities of existing maximum entropy model variations with the GLM framework, and are tailored for the analysis of October 16, 2023 20/27 continuously acquired neuronal data.As Truccolo's comparisons in [20] suggest, GLMs account for temporal dynamics explicitly in modeling ensemble spiking, and thus are arguably more predictive than maximum entropy models.Within the context of the MkPP mGLM we utilized, temporal dynamics of neuronal spiking were modeled as relevant covariates in the estimation of ensemble spiking events.Such a model can be simplified to exclude spiking history, as demonstrated by the history-independent model; and can be expanded to model the influence of stimuli, as previously addressed for maximum entropy models [26].
Due to the large number of possible interactions, challenges in the tractability of synchrony analyses are inherent, particularly when modeling the effects of relevant covariates.Incorporating the emphasis on reliable interactions, as proposed in [27], model complexity may be managed in a data-driven fashion.The proposed adaptive greedy filtering algorithm for sparse model estimation ensures only the salient effects of covariates are captured.The adaptive filtering algorithm also characterizes dynamics in network correlational structure, analogous to Bayesian state space filtering algorithms [24,25], and is thus applicable in the analysis of non-stationary neuronal processes.In lieu of constructing credible intervals around the aforementioned Bayesian estimates, we utilize a statistical test for which the test statistic's limiting distribution is precisely characterized.Unlike existing analyses, the proposed statistical tests do not require repeated trials of data to detect coordinated spiking activity, and are thus suitable for the analysis of continuous recordings of ensemble neuronal spiking.
Extending previous results in high-dimensional statistics, we have shown in Theorem 1 that the elegant procedure of [49] for LASSO estimation may be adapted to de-sparsify OMP estimates, and that de-sparsified estimates are asymptotically normal.
In reviewing the existing literature, we noted a paucity in work on variable selection algorithms concerning the construction of confidence intervals.The OMP has been shown to have similar consistency properties as LASSO regression under appropriate conditions [41,42]; however, in settings with large quantities of data, the latter becomes intractable.The result established by Theorem 1 enables the construction of confidence intervals around OMP-estimated parameters in order to provide analogous methods of statistical inference as LASSO for an algorithm suitable in settings with large data sets, addressing this gap in the high-dimensional statistics literature.

Novel Insights into Coordinated Network Activity
The proposed modeling and statistical inference framework constitute a novel approach to studying coordinated neuronal spiking by enabling the adaptive analysis of continuously acquired or single-trial data.The applications presented in this work demonstrate that tracking dynamics in higher-order coordinated spiking activity provides new perspectives on the latent structure of neuronal assemblies during fast state transitions.
Analyzing spontaneous ensemble spiking recorded during the transition into propofol-induced anesthesia provided greater detail about the correlation structure of local neuronal networks during both the conscious and unconscious states, revealing that simultaneous spiking activity is induced differently in each state.Similarly, analyzing two populations of neurons in rat motor cortex during one sleep cycle revealed distinctions in latent higher-order coordinated structure between them.While periods of endogenously facilitated simultaneous spiking in excitatory neurons aligned mainly with state transitions, higher-order interactions between inhibitory neurons were most prominently exogenously coordinated during sleep.The ability to track dynamics and detect exogenous influences on ensemble spiking with statistical confidence provides a new approach to probing the neural mechanisms underlying transitions between and characteristics of arousal states.October 16, 2023 21/27

Extensions
The proposed statistical inference framework was developed to test for significant coordination of r th -order spiking events, and the presented results demonstrated its efficacy.Specifically, Theorem 2 characterized the limiting distributions for the adaptive de-biased deviance difference test statistic under both outcomes of a nested hypothesis test in which the null hypothesis restricted parameters to impose conditionally independent r th -order spiking.However, a nested null hypothesis can, in principle, be constructed to impose different assumptions.An immediate extension of the proposed analysis could include spatial information, for example, so that a null hypothesis assumes r th -order spiking amongst a spatially localized subset of a recorded neuronal assembly is conditionally independent.The proposed inference framework was hence established to readily extend to any nested hypothesis test in Corollary 2.2.
An important consequence of this corollary result is that it provides a theoretical foundation for adaptive Granger causality using greedy algorithms.Since the proposed methods utilize a multinomial extension of generalized linear models, Corollary 2.2 establishes the asymptotic result in [31] for greedy parameter estimates in the limiting case of a single-neuron model.Notably though, Corollary 2.2 also implies that a nested hypothesis test can be formulated to determine if exogenous signals, such as sensory stimuli or concurrent activity in other brain regions, have Granger-causal effects on a neuronal network or its subsets.Thus, the methods proposed in the present study can be extended to investigate the local network effects of global neural dynamics.

H 0 :Fig 1 .
Fig 1.Ensemble spiking is mapped to a disjoint representation of simultaneous spiking events.The proposed methods (grey box) are used to infer the strength of higher-order coordination amongst C neurons in a dynamic fashion.
t consists of an ensemble history-modulation vector θ (m) t and the baseline firing parameter, µ (m) t k ), i.e. non-central chi-square, where ν (r) k is the non-centrality parameter at time k that depends only on the true parameters, and M (r) = |K r | is the difference in the cardinalities of the full and reduced support sets.Our theoretical results are comprehensively discussed in supporting information (S1 Appendix).
under H 1 , we must estimate the non-centrality parameter for each window.Assuming it evolves smoothly in time, we use a state-space smoothing algorithm[31] to estimate ν (r) that the rejection of the null is a stronger indication of coordination than for smaller values of J (r)

Fig 2 .Fig. 2 ,
Fig 2.Theoretical versus empirical distributions of the adaptive de-biased deviance difference under the null (blue) and alternative (red) hypotheses.Empirical distributions of adaptive de-biased deviance differences were compiled from history-dependent analyses of 50 realizations of a discretized marked point process simulating the ensemble spiking of 5 neurons with exogenously induced 3 rd -order spiking dynamics.The distributions were compared to the probability distribution functions (PDF) of a χ 2 distribution with M = 10 degrees of freedom (left), and a non-central χ 2 distribution with M = 10 degrees of freedom and non-centrality parameter ν = 327 (right), corresponding respectively to the limiting distribution under the null and alternative hypotheses.

Fig 3 .
Fig 3. Analysis of ensemble spiking with 3 rd -order coordination induced exogenously by a square wave. A. Simulated ensemble spiking of five neurons.B. Sum of the r th -order simultaneous spiking events for r = 2, 3, 4, 5. C. Spiking coordination varies across 4 epochs defined by the exogenous process that alternatingly facilitated and suppressed 3 rd -order events.These are demarcated by vertical dashed lines.D. Significant r th -order coordination neglecting ensemble history.E. Significant r th -order coordination based on history-dependent ensemble spiking model.Statistical testing in D-E performed at level α = 0.01.F. Average Pearson correlation with 95% confidence interval.G. Average spiking regularity: coefficient of variation ±2 SEM.H. Average mark CIF differences of 3 rd -order (green) spiking interactions ±2 SEM.

Fig 4 .
Fig 4.Analysis of ensemble spiking with 3 rd -order coordination induced exogenously by an autoregressive process.A. Simulated ensemble spiking of five neurons.B. Sum of the r th -order simultaneous spiking events for r = 2, 3, 4, 5. C.An autoregressive process was used to exogenously induce 3 rd -order spiking coordination.The effect was most evident when the exogenous variable had value larger than 2 over an interval is demarcated by vertical dashed lines.D. Significant r th -order coordination neglecting ensemble history.E. Significant r th -order coordination based on history-dependent ensemble spiking model.Statistical testing in D-E performed at level α = 0.01.F. Average Pearson correlation with 95% confidence interval.G. Average spiking regularity: coefficient of variation ±2 SEM.H. Average mark CIF differences of 3 rd -(green) and 4 th -order (teal) spiking interactions ±2 SEM.

Fig 5 .
Fig 5. Higher-order spiking coordination analysis of human cortical neurons during anesthesia.A. Ensemble spiking of 8 neurons aligned to loss of consciousness (LOC) at 0s induced by administering propofol.A second administration of propofol occurred at ∼ 250s.B. Sum of the r th -order simultaneous spiking events for r = 2, . . ., 8. C. Significant r th -order coordination neglecting ensemble history.D. Significant r th -order coordination based on history-dependent ensemble spiking model.Statistical testing in C-D performed at level α = 0.01.E. Average Pearson correlation with 95% confidence interval.F. Average spiking regularity: coefficient of variation ±2 SEM.

Table 1 .
1, bottom panel); the following subsections describe each analysis stage.Key notation used throughout the remainder of the paper is summarized in Algorithm development and theoretical results are comprehensively

Table 1 .
Summary of key notation t , . . ., n (C) t Ensemble spiking observation at time bin t of C neurons λ (c) t ∆ Conditional Intensity Function (CIF) of c th neuron