A dynamic generative model can extract interpretable oscillatory components from multichannel neurophysiological recordings

eLife assessment

This method article proposes a valuable oscillation component analysis (OCA) approach, in analogy to independent component analysis (ICA), in which source separation is achieved through biophysically inspired generative modeling of neural oscillations. The empirical evidence justifying the approach's advantage is solid. This work will be of interest to researchers in the fields of cognitive neuroscience, neural oscillation, and MEG/EEG.

https://doi.org/10.7554/eLife.97107.3.sa0

Significance of the findings:

Valuable: Findings that have theoretical or practical implications for a subfield

Landmark
Fundamental
Important
Valuable
Useful

Strength of evidence:

Solid: Methods, data and analyses broadly support the claims with only minor weaknesses

Exceptional
Compelling
Convincing
Solid
Incomplete
Inadequate

During the peer-review process the editor and reviewers write an eLife Assessment that summarises the significance of the findings reported in the article (on a scale ranging from landmark to useful) and the strength of the evidence (on a scale ranging from exceptional to inadequate). Learn more about eLife Assessments

Abstract
Introduction
Theory
Results
Discussion
Materials and methods
Appendix 1
Appendix 2
Appendix 3
Appendix 4
Data availability
References
Article and author information
Metrics

Abstract

Modern neurophysiological recordings are performed using multichannel sensor arrays that are able to record activity in an increasingly high number of channels numbering in the 100s to 1000s. Often, underlying lower-dimensional patterns of activity are responsible for the observed dynamics, but these representations are difficult to reliably identify using existing methods that attempt to summarize multivariate relationships in a post hoc manner from univariate analyses or using current blind source separation methods. While such methods can reveal appealing patterns of activity, determining the number of components to include, assessing their statistical significance, and interpreting them requires extensive manual intervention and subjective judgment in practice. These difficulties with component selection and interpretation occur in large part because these methods lack a generative model for the underlying spatio-temporal dynamics. Here, we describe a novel component analysis method anchored by a generative model where each source is described by a bio-physically inspired state-space representation. The parameters governing this representation readily capture the oscillatory temporal dynamics of the components, so we refer to it as oscillation component analysis. These parameters – the oscillatory properties, the component mixing weights at the sensors, and the number of oscillations – all are inferred in a data-driven fashion within a Bayesian framework employing an instance of the expectation maximization algorithm. We analyze high-dimensional electroencephalography and magnetoencephalography recordings from human studies to illustrate the potential utility of this method for neuroscience data.

Introduction

Human neurophysiological recordings such as scalp electroencephalogram (EEG), magnetoencephalogram (MEG), stereoelectroencephalogram (SEEG), local field potentials, etc., consist of ∼10² of sensors that record mixtures of predominantly cortical network oscillations (Buzsáki, 2006; Buzsáki and Draguhn, 2004; Lopes da Silva, 2013; Wang, 2010; Helfrich and Knight, 2016). The network oscillations have distinct spatio-temporal signatures, based on the functional brain areas involved, their interconnections, and the electromagnetic mapping between the source currents and the sensors. However, the source-to-sensor mixing and the further superposition of measurement noise complicate the interpretation of the sensor-level data and its topography (Schaworonkow and Nikulin, 2022). Given the widespread and growing availability of high-density neural recording technologies (Stevenson and Kording, 2011), there is clearly a pressing need for analysis tools that can recover underlying dynamic components from highly multivariate data.

This problem fits within a larger class of blind source separation (BSS) problems for which there are a plethora of component analysis algorithms that attempt to extract underlying source activity as linear weighted combinations of the sensor activity. The decomposition weights are designed according to some predefined criterion on the extracted time series depending on the application. Independent component analysis (ICA) (Hyvärinen and Oja, 2000; Jung et al., 2000; Cardoso, 1999; Hyvärinen, 2013) has been particularly popular within neuroscience for a number of reasons. It requires no assumption of the original sources except for statistical independence, that is, no or minimal mutual information between the components (Bell and Sejnowski, 1995), and in principle requires little to no user intervention to run. However, a drawback of the method is that it assumes that the samples of each component time trace are independent and identically distributed, which is not generally true in physical or physiological applications. This leads to another major drawback of ICA: it relies on the cumulative histograms of sensor recordings and if these histograms are Gaussian, ICA is in principle unable to separate such sources (Brookes et al., 2011). This includes, for example, physiological signals such as EEG, MEG, sEEG, etc., that are at least approximately Gaussian distributed due to the central limit theorem (Muirhead, 1982) since these signals reflect the linear combination of many sources of activity. Finally, given the lack of assumptions on the structure of the underlying signals, there is no guarantee that ICA will extract relevant oscillatory components or any other highly structured dynamic pattern for that matter.

Given the prominence and ubiquity of oscillations in neurophysiological data, applications of component analysis methods in neuroscience have taken a different approach to emphasize oscillatory dynamics, employing canonical correlation analysis (Robinson et al., 2017), generalized eigenvalue decomposition (Parra and Sajda, 2003), joint decorrelation (de Cheveigné and Parra, 2014), or similar methods to identify a set of spatial weights to maximize the signal-to-noise ratio in the extracted component within a narrow band around a given frequency of interest (Nikulin et al., 2011; de Cheveigné and Arzounian, 2015; Cohen, 2017; Cohen, 2018). The extracted component then inherits the intrinsic oscillatory dynamics around that frequency without the need for a pre-designed narrow-band filter (Yeung et al., 2007). These techniques do acknowledge the inherent temporal dynamics of the oscillatory sources, but do so via nonparametric sample correlation or cross-spectrum matrix estimates that are sensitive to noise or artifacts and that require substantial amounts of data to achieve consistency.

Another problem common to all of these component decomposition methods is that they do not directly estimate the source to sensor mixing matrix. Instead, they estimate spatial filters that extract the independent components, which are not directly interpretable. The source to sensor mixing matrix can only be estimated by solving an inverse problem that requires knowledge of the noise covariance matrix (Haufe et al., 2014), which adds another complication and source of potential error. Finally, the properties of the extracted component time courses are not known a priori and must be assessed after the fact, typically by using nonparametric tests as well as visual inspection.

These difficulties with component selection and interpretation occur in large part because existing BSS methods lack a generative model for the underlying spatio-temporal dynamics. With a probabilistic generative model, it is possible to specify a soft constraint on the dynamic properties of the underlying components, while maintaining other desirable properties such as independence between components. Component selection can be handled automatically within a statistical framework under the model, and interpretation is straightforward in principle if the components can be described by a small number of parameters. Here, we propose a novel component analysis method that uses a clever state-space model (Wiener, 1966; Harvey, 1990; Harvey and Trimbur, 2003) to efficiently represent oscillatory dynamics in terms of latent analytic signals (Bracewell, 2000) consisting of both the real and imaginary components of the oscillation. The observed data are then represented as a superposition of these latent oscillations, each weighted by a multichannel mixing matrix that describes the spatial signature of the oscillation. We estimate the parameters of the model and the mixing matrices using generalized expectation maximization (GEM) and employ empirical Bayes model selection to objectively determine the number of components. In these ways, we address the major shortcomings described above for many component analysis methods. We refer to our novel method as oscillation component analysis (OCA) akin to ICA. In what follows, we describe the model formulation in detail and demonstrate the performance of the method on simulated and experimental data sets, including high-density EEG during propofol anesthesia and sleep, as well as resting-state MEG from the Human Connectome Project.

Theory

State-space oscillator model

Oscillatory time series can be described using the following state-space representation (Wiener, 1966; Harvey, 1990; Harvey and Trimbur, 2003):

\begin{matrix} x_{n} = a R (f) x_{n - 1} + υ_{n}, & υ_{n} \sim N_{2} (0, σ^{2} I); \\ y_{n}^{(l)} = c_{l}^{T} x_{n} + ϵ_{n}, & ϵ_{n} \sim N (0, R), \end{matrix}

where the oscillation state, $x_{n} = [x_{n, 1}, x_{n, 2}]^{⊤}$ , is a two-dimensional state vector. The stochastic difference equation summarizes oscillatory dynamics using random rotation in two-dimensional state space through a deterministic rotation matrix, $R$ , explicitly parameterized by the oscillation frequency, $f$ , the sampling rate, $f_{s}$ , and the damping factor, a, as

R (f) = [\begin{matrix} \cos 2 π f / f_{s} & - \sin 2 π f / f_{s} \\ \sin 2 π f / f_{s} & \cos 2 π f / f_{s} \end{matrix}]

and a stochastic driving noise, assumed to be stationary and Gaussian with variance $σ^{2}$ . These elements of this state vector trace out two time series that maintains an approximate $π / 2$ radian phase difference, and therefore are closely related to the real and imaginary parts of an analytic signal (Bracewell, 2000) in the complex plane (see Appendix 2, section ‘Oscillation states and analytic signals’). The oscillation states are henceforth called analytic signal with minor abuse of notation. An arbitrary fixed projection (i.e., on the real line) of the state vector realizations generates the observed noisy oscillation time series at the sensor. Multiple oscillations can be readily incorporated in this state-space model by simply considering their linear combination. Recently, several investigators (Matsuda and Komaki, 2017a; Beck et al., 2018; Beck et al., 2022) utilized this state-space representation to extract underlying oscillatory time courses from single-channel EEG time traces.

Generalization to multichannel data

In order to represent multichannel neural recordings, we employ this oscillatory state-space representation within a BSS model (Parra, 1998). In the proposed generative model, a sensor array with L sensors records neurophysiological signals produced by superimposition of M distinct oscillations supported by underlying brain networks or circuits, which we will refer to as oscillation sources.

\begin{aligned} (3) & \begin{matrix} [\begin{matrix} x_{t}^{(1)} \\ ⋮ \\ x_{t}^{(M)} \end{matrix}] = [\begin{matrix} a^{(1)} R (f^{(1)}) & \dots & 0 \\ ⋮ & ⋱ & ⋮ \\ 0 & \dots & a^{(M)} R (f^{(M)}) \end{matrix}] [\begin{matrix} x_{t - 1}^{(1)} \\ ⋮ \\ x_{t - 1}^{(M)} \end{matrix}] + [\begin{matrix} υ_{t}^{(1)} \\ ⋮ \\ υ_{t}^{(M)} \end{matrix}] \\ [\begin{matrix} y_{t}^{(1)} \\ ⋮ \\ y_{t}^{(L)} \end{matrix}] = [\begin{matrix} c_{1, 1}^{⊤} & c_{1, 2}^{⊤} & \dots & c_{1, M}^{⊤} \\ c_{2, 1}^{⊤} & c_{2, 2}^{⊤} & \dots & c_{2, M}^{⊤} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ c_{L, 1}^{⊤} & c_{L, 2}^{⊤} & \dots & c_{L, M}^{⊤} \end{matrix}] [\begin{matrix} x_{t}^{(1)} \\ x_{t}^{(2)} \\ ⋮ \\ x_{t}^{(M)} \end{matrix}] + [\begin{matrix} ϵ_{t}^{(1)} \\ ⋮ \\ ϵ_{t}^{(L)} \end{matrix}] \end{matrix} \end{aligned}

Each oscillation source is governed by the above-mentioned state-space representation, resulting in a structured 2M dimensional state space. The structure of the underlying brain networks or circuits governs how each oscillation source is observed at the sensor array, via a $M \times 2$ spatial distribution matrix, $c_{*, m} = [c_{1, m}, c_{2, m}, \dots, c_{L, m}]$ . In other words, the l th electrode in the sensor array observes a specific projection of the underlying mth analytic oscillation given by $c_{l, m} = [c_{l, m, 1}, c_{l, m, 2}]$ , which encodes the amplitude and phase of the mth oscillation time course at that electrode. This probabilistic generative model for multichannel recordings is motivated by a potential biophysical mechanism of electromagnetic traveling wave generation in brain parenchyma (Galinsky and Frank, 2020; see Appendix 1, section ‘Mechanistic origin’).

An illustration

In Figure 1A, middle panel, we depict an oscillation state as an analytic signal, $x_{n}$ (black arrows), in 2D state space that is rotating around the origin according to the given state-space model of left panel: the red dashed line traces the oscillation states over time. The actual amount of rotation between every pair of consecutive timepoints is a random variable centered around $2 π f / f_{s} = 0.22 π$ , with the spread determined by damping parameter, $a = 0.99$ , and process noise covariance. The right panel shows two noisy measurements, reflecting two different projections: one on the real axis (blue traces), another on the line making $π / 4$ radian angle to the real axis (orange traces). Because of this angle between the lines of projection, these two measurements maintain approximately $π / 4$ radian phase difference throughout the time course of the oscillation.

Figure 1

Download asset Open asset

From state-space oscillator model to oscillation component decomposition.

(A) An illustrative example of a multichannel state-space oscillation model: a single oscillation realized as an analytic signal $x_{n}$ is measured as two different projections having radian phase difference. (B) Graphical representation of the probabilistic generative model describing the oscillations as dynamic processes that undergo mixing at the sensors and that are observed with additive Gaussian noise. (C) Graphical representation of the variational Bayes approximation that allows iterative closed-form inference. (D) Oscillation component analysis fitting and reconstruction pipeline for experimentally recorded neurophysiological data. The pipeline exposes a number of methods for ease of analysis, that is, for fitting the oscillation component analysis (OCA) hyperparameters fit () method, which accepts the sensor recordings and an initial sensor noise covariance matrix as input, extracting the oscillation time courses get_sources(), reconstructing a multichannel signal from any arbitrary subset of oscillation components, apply(), getting a final noise covariance estimate from the residuals of OCA fitting, get_noise_covariance() method, etc.

We note that several earlier works Matsuda and Komaki, 2017b; Quinn et al., 2021 have similar multivariate oscillator models, albeit from the perspective of a spatio-spectral eigendecomposition of the companion form of a multivariate autoregressive (MVAR) parameter matrix (Neumaier and Schneider, 2001). As described earlier, we employ this state-space form here in the context of a BSS problem.

Next, we briefly describe how one can infer the hidden oscillator states from observed multivariate time-series dataset given the oscillation state-space model parameters and potentially adjust the oscillation state-space model parameters to the dataset. We defer the mathematical derivations to Appendix 3 to maintain lucidity of the presentation.

Learning algorithm

Priors

We assume simple prior distributions on the measurement noise and sensor-level mixing coefficients to facilitate stable and unique recovery of the oscillation components. To obtain an $M$ -oscillator analysis of $L$ -channel data, we consider a Gaussian prior on the spatial mixing components, $c_{l, m}$ , with precision α, and an inverse-Wishart prior on sensor noise covariance matrix, $R$ , with scale, $Ψ$ , and degrees of freedom, $ν$ :

\begin{array}{ll} p (C ∣ α, M) & = {(\frac{α}{2 π})}^{M L} \exp (- \frac{α}{2} \sum_{l \geq 1, m \geq 1} ‖ c_{l, m} ‖^{2}) \\ p (R ∣ M) & = γ_{ν, L} {| Ψ |}^{ν / 2} {| R |}^{- (ν + L + 1) / 2} \exp - \frac{1}{2} Tr {Ψ R^{- 1}} \end{array}

We treat the oscillation parameters, $(f^{(m)}, a^{(m)}, {(σ^{2})}^{(m)}) := θ^{(m)}$ , and distributional parameters of the assumed priors, $α$ , $Ψ, ν$ , as hyperparameters. Figure 1B shows the probabilistic graphical model portraying oscillations as a dynamical system evolving in time, with the priors as parent nodes.

Variational Bayes inference

Unfortunately, given the priors and the interplay between oscillation source time courses and sensor-level mixing patterns, exact computation of the log-likelihood, and thus the exact posterior, is intractable. We therefore employ variational Bayes (VB) inference (Quinn and Šmídl, 2006), a computationally efficient inference technique to obtain a closed-form approximation to the exact Bayes posterior. Originally introduced by Hinton and van Camp, 1993, VB inference simplifies the inference problem through a restrictive parameterization that reduces the search space for distributions, splitting up the problem into multiple partial optimization steps that are potentially easier to solve (Attias, 1999). The name comes from the negative variational free energy, also known as evidence lower bound (Neal and Hinton, 1998), used to assess the quality of the aforementioned approximation and as a surrogate for the log-likelihood (see Appendix 3, section ‘Negative variational free energy’).

In particular, given the number of oscillatory components $M$ and the corresponding hyperparameters $θ = [θ^{(1)}, θ^{(2)}, \dots, θ^{(M)}]$ and $α, Ψ, ν$ , we use the following VB decoupling (Attias, 1999) of the posterior distribution of mixing matrix $C$ , oscillation states $x_{t}$ , and noise covariance matrix $R$ :

q (C, {x_{t}}, R ∣ α, θ, M) = q (C ∣ M) q ({x_{t}} ∣ M) q (R ∣ M)

This particular choice allows the approximate posteriors of these quantities to be represented in closed form by multivariate Gaussian, multivariate Gaussian, and inverse-Wishart distributions, respectively (see Appendix 3, section ‘Variational Bayes inference’). Figure 1C shows the graphical representation of posterior distribution after the VB decoupling. This essentially allows us to perform an iterative posterior inference procedure, where we cyclically update the posteriors $q (x_{t} ∣ M)$ , $q (C ∣ M)$ , and $q (R ∣ M)$ using the latest sufficient statistics from other two distributions (see Appendix 3, section ‘Variational Bayes inference’ for more details).

Generalized EM to update hyperparameters

Since the parameters of the state-space model, and of the assumed prior distributions, are not known a priori, we then obtain their point estimates using an instance of the GEM algorithm, which utilizes the aforementioned approximate inference in the E-step (Dempster et al., 1977; Neal and Hinton, 1998). We start the learning algorithm by initializing θ, $\bar{Λ}$ , $\bar{C}$ , α to appropriate values (see Appendix 3, section ‘Initialization of parameters and hyperparameters’ for details). We then cyclically update the posteriors $q (x_{t} ∣ M)$ , $q (C ∣ M)$ , and $q (R ∣ M)$ until these iterations stop changing the negative variational free energy. At the end, we update the hyperparameters θ, α from the current inference and proceed to the next inference with the updated hyperparameters. These update rules form an outer update loop that refines the hyperparameters, within which operates an inner update step that iteratively improves the approximate posterior distribution of the model parameters. We defer the update rules to Appendix 3, section ‘Generalized EM algorithm’. The algorithm terminates when the hyperparameter adjustments no longer increase the free energy of the model, that is, we resort to parametric empirical Bayes estimation (see Appendix 3, section ‘Empirical Bayes inference and model selection).

Selecting optimal number of oscillations

The learning algorithm assumes that the number of oscillation sources, $M$ , is known a priori, which is rarely the case for experimentally recorded data. In theory, our probabilistic treatment could assign a log-likelihood for a given dataset to every model structure, that is, the number of oscillation sources, that can be used as a goodness-of-fit score. Unfortunately, the log-likelihood is intractable to evaluate; however, the negative variational free energy can be used as a surrogate.

We formalize this idea by considering that the model structure, M, is drawn from a discrete uniform distribution over a finite contiguous set, $(1, M_{m a x})$ :

\begin{matrix} p (M) = 1 / M_{m a x}, & 1 \leq M \leq M_{m a x}, \end{matrix}

with $M_{m a x}$ being the maximal number of oscillation sources. This allows us to compute the posterior on the model structure, $q (M)$ , within the same variational approximation framework. Using the negative variational free energy expression of the oscillator model, it is easy to show that $q (M)$ is proportional to the exponential of the negative variational free energy of the M-oscillator model (see Appendix 3, section ‘Model structure posterior’). Once the model posteriors are obtained, one can choose to simply select the model with maximum negative variational free energy (empirical Bayes model selection, see Appendix 3, section ‘Empirical Bayes inference and model selection’), or locate the knee or jump point that exhibits a significant local change in $q (M)$ (for such an example, see the DiffBic function; Zhao et al., 2008).

Oscillation component analysis

A preliminary version of this algorithm has been presented previously (Das and Purdon, 2022). In order to analyze experimental data, we devised the a standardized pipeline as demonstrated in Figure 1D closely following the MNE ICA pipeline (Gramfort et al., 2014). We refer to this pipeline as OCA. The pipeline, implemented as Python class OCA, accepts the sensor recordings either as continuous raw data or a set of discrete epochs, each with an initial sensor noise covariance matrix. The hyperparameters for the measurement noise covariance prior are derived from the supplied initial sensor noise covariance matrix and are not updated in this pipeline. The sensor data is standardized (i.e., pre-whitened) against the sensor noise covariance: this step also applies all active signal-space projectors to the data if present such as the average reference (Uusitalo and Ilmoniemi, 1997). The pre-whitened data is then decomposed using PCA. The first n_components are then passed to the parameter learning iterations for hyperparameters θ, α, followed by the joint inference of the mixing matrix, oscillation time courses, and residual noise covariance. The optimal number of oscillations is then selected as described earlier. The final noise covariance is computed from the remaining PCA components that are not supplied to oscillation fitting and the OCA residuals. Once the oscillation hyperparameters and the mixing matrix have been estimated, the oscillation time courses can be estimated from any continuous recordings or discrete epochs. Finally, the pipeline can also reconstruct a multichannel signal from any arbitrary subset of oscillation components. We emphasize here that the main reason for using PCA here is to perform rotation of the data, so that (1) rank-deficient data (i.e., average referenced EEG data) can be handled conveniently by throwing out the zero-variance component. (2) The variance the multichannel data can be uniformly distributed in the remaining PCs. For these reasons, we retain components explaining 99.9% variance of the multichannel recordings, that is, almost all components with nonzero variance. The PCA step essentially performs rotation of the data in the spatial dimension, leaving the temporal relationships untouched.

Results

Finally, we apply OCA on a simulation example and several experimentally recorded datasets. We demonstrate the versatility of OCA analysis using two EEG datasets, one under propofol-induced anesthesia, another during sleep, and one resting-state MEG dataset. The simulation example shows the utility of such time-domain analysis over frequency-domain analysis in the context of multichannel recordings. In the real data applications, we showcase the implications of OCA as a mechanistically grounded tool, suggesting various downstream analysis involving the oscillation time courses and their sensor distributions.

Simulation example

We first apply the proposed oscillation decomposition method on a synthetic dataset to illustrate how our time-domain approach differs from traditional frequency-domain approach. We generated a synthetic EEG dataset using a 64-channel montage, forward model and noise covariance matrix from a sample dataset distributed with the MNE software package (Gramfort et al., 2014). To select the active current sources, four regions with $5 {m m}^{2}$ area were centered on chosen labels in the DKT atlas (Desikan et al., 2006): ‘transversetemporal-lh’, ‘precentral-rh’, ‘inferiorparietal-rh’, and ‘caudalmiddlefrontal-lh’. Three time courses were separately generated at a sampling frequency of 100 Hz according to univariate AR(2) dynamics tuned to generate stable oscillations at 1.6 Hz (slow/delta), 10 Hz (alpha), and 12 Hz (alpha) frequencies. The first two regions were simulated with the same slow/delta time course such that the activity in the ‘precentral-rh’ area (Figure 2A, orange) lags that of ‘transversetemporal-lh’ area (Figure 2A, blue) by 10 ms. The last two regions, ‘inferiorparietal-rh’ and ‘caudalmiddlefrontal-lh’, were simulated with two independent alpha time courses. We projected this source activity to the EEG sensors via a lead-field matrix and add spatially colored noise generated with the covariance structure in Figure 2D to the sensor data. Two 20 s epochs were chosen for OCA and model selection was performed within the range of (2, 3, 4, 5, 6, 8, 10) oscillations. Figure 2B presents the two most common frequency-domain visualizations of multichannel data: the top panel shows the power distribution over the scalp within different canonical EEG frequency bands, while the channel-wise power spectral density (PSD) is shown in the bottom panel. Figure 2C demonstrates how the decomposition obtained by OCA provides an accurate delineation of different oscillation time courses and effective characterization of their spatial distributions in the sensor-level mixing maps. Unlike the frequency-domain techniques, the residual time series that follows OCA extraction provides the temporally unstructured part of the multichannel recording. The estimated covariance matrix from these residuals closely resembles the original sensor noise covariance matrix that was used to corrupt the recordings (see Figure 2E). The covariance matrix estimate is strikingly close to the true covariance matrix (i.e., compare with Figure 2D). OCA provides a measure, $q (M)$ , for model structure selection that objectively identifies the existence of three oscillation time courses with different time dynamics (see Figure 2F).

Figure 2

Download asset Open asset

Simulation study.

(A).Four neural sources carrying three oscillations, where the orange time course is a lagged instance of blue time course, red and green time courses are independent. (B) Power spectral density of the simulated electroencephalogram (EEG) recording and power distribution over the EEG sensors in different frequency bands. (C) Recovered oscillation components and their sensor-level maps. (D) Sensor noise covariance for the simulation. (E) Estimated sensor noise covariance matrix. (F) Model structure selection via model structure posterior $q (m)$ .

How does OCA compare to conventional approaches?

We performed additional analyses to compare OCA with Fourier-based frequency domain method and ICA in a different realization of a similar synthetic EEG dataset (see Figure 3A). The frequency-domain approach applied here is based on the multitaper method: given a bandwidth parameter of 2 Hz, power spectrum density is computed for individual channels. The oscillatory activity can be identified visually as the peaks in the power spectrum plot in Figure 3B, top panel, and its scalp distribution can be obtained by averaging the power within given frequency bands around the peaks in each channel (as shown in Figure 3B, bottom panels). However, these visualizations do very little to identify underlying oscillatory sources and pose additional questions. For example, the left topographical plot hints to two possible sources, but provides no evidence if they are separated spatially, temporally or both. Similarly, the right two topographical plots are almost identical and show no separability between these the oscillatory sources that generate two distinct peaks in spectrum. Regarding ICA, we use the ‘extended infomax’ (Lee et al., 1999) implementation provided by MNE-python 1.2 (Gramfort et al., 2014) to see if ICA can distinguish the three statistically independent sources in these data. The leading four ICA components are shown in Figure 3C. The ICA-identified components mix the underlying independent signals, which is not surprising since all components follow a Gaussian distribution. Meanwhile, OCA is able to recover three distinct oscillatory components consistent with the ground truth. OCA is able to do so because it explicitly models temporal dynamics of the components.

Figure 3

Download asset Open asset

Simulation study (extended).

(A) Four neural sources carrying three oscillations, similar to Figure 2. (B) Power spectral density of the simulated electroencephalogram (EEG) recording and power distribution over the EEG sensors in the frequency bands (red overlay) around visually identifiable peaks. (C) Recovered independent component analysis (ICA) components (left, middle, and right columns show topographic maps, power spectrum density, and time courses, respectively). (D) Recovered oscillation component analysis (OCA) components (the topographic maps show the magnitude [left] and phase [right], while line plots show power spectrum density [left] and time courses [right], respectively).

We repeated similar analyses using real human EEG recording and found a similar result where ICA produced components that mixed slow and alpha band signals, whereas OCA identified distinct oscillatory components (see Appendix 4, section ‘Comparison of OCA to traditional approaches in experimental EEG data’).

EEG recording during propofol-induced unconsciousness

Next, we demonstrate utility of OCA on experimental data using EEG recordings from a healthy volunteer undergoing propofol-induced unconsciousness, previously described in Purdon et al., 2013. For induction of unconsciousness, the volunteer underwent a computer-controlled infusion of propofol to achieve monotonically increasing levels of effect-site concentration in steps of $1 µ g {m L}^{- 1}$ . Each target effect-site concentration level was maintained for 14 min. The EEG was recorded using a 64-channel BrainVision MRI Plus system (Brain Products) with a sampling rate of 5000Hz, bandwidth 0.016–1000 Hz. The volunteers were instructed to close their eyes throughout the study. Here, we investigated EEG epochs during maintenance of effect-site concentrations of $0 µg {m L}^{- 1}$ (i.e., baseline, eyes closed), $2 µg {m L}^{- 1}$ , and $4 µg {m L}^{- 1}$ for OCA. We selected 10 clean 3.5 s epochs corresponding to those target effect-site concentration (see ‘Materials and methods’).

We fitted OCA models with 20, 25, 30, 35, 45, 50, 55, 60 oscillation components and selected the model with highest negative variational free energy as the best OCA model. Following these empirical Bayes criteria, we selected OCA models with 30, 30, 50 components for these three conditions, respectively. When the center frequencies were grouped within the canonical frequency bands (Buzsáki, 2006), we found 17, 17, 29 slow/delta (0.1–4 Hz) oscillation components and 17, 9, 15 alpha (8–13 Hz) oscillation components for the three conditions, respectively. Figure 4A–C shows PSDs of reconstructed EEG activity within each band, defined as the cumulative projection of these grouped oscillation components to the sensors. Since each of these oscillations have different sensor distributions, the increasing number of oscillation components can be associated with fragmentation of neuronal networks that support these oscillations.

Figure 4

Download asset Open asset

Oscillation component analysis (OCA) of the electroencephalogram (EEG) from a healthy volunteer undergoing propofol-induced unconsciousness.

Conditions of target effect-site concentration of (A, D) 0 (i.e., baseline), (B, E) $2 µg {m L}^{- 1}$ , and (C, F) $4 µg {m L}^{- 1}$ are analyzed. Panels (**A–C**) show the power spectral densities (PSDs) of reconstructed EEG activity within each canonical band. Panels (**D–F**) show the three dominant (in terms of sensor wide power) alpha component: the topographic maps show the magnitude (left) and phase (right) distribution of sensor-level mixing, the time courses are 1 s representative example of the extracted oscillations from the selected epochs. The black bars on the right display the coherency measure within alpha band. OCA correctly identifies that the spatial mixing sensor maps of the alpha waves (8–12 Hz) are oriented posteriorly at baseline, but gradually become frontally dominant under propofol. The sensor weights are scaled to have maximum value 1. So, the units of time series traces can be considered to be in µV.

We also quantified the coherency of the reconstructed EEG activity within these commonly used frequency bands as the ratio of the power of the strongest component to the total power of the oscillation components within the bands. For the data analyzed here, slow + delta band coherency decreases gradually while the alpha band coherency increases with the increasing propofol effect-site concentration. To investigate this observation further, we visualized three dominant alpha components from each conditions in order of their strengths in Figure 4D–F. In each column, the left sub-panels consist of two topographic plots: the left one showing the distributions of the strength (scaled between 0 to 1) and the right one showing the relative phase of the estimated analytic oscillation mixing, $c_{i, j}$ . The right sub-panels display examples of the extracted analytic oscillations, that is, estimated real and imaginary components of 1 s long oscillation segments. The rightmost black bars show the coherency measure within alpha band. Clearly, during the $4 µg {m L}^{- 1}$ effect-site concentration, the amplitude of the leading alpha oscillation component is much bigger than the others, which aligns with our observation of the coherency measure. This finding previously reported analyses showing a rise of slow + delta oscillation power, rapid fragmentation of a slow + delta oscillation network, and emergence of global alpha network during propofol-induced anesthesia (Lewis et al., 2012; Murphy et al., 2011). Further, as the propofol effect-site concentration increases, the temporal dynamics increasingly resemble a pure sinusoidal signal, as evident from the PSD of the alpha band reconstruction that becomes more concentrated around the center frequency. Lastly, the topographic plots in Figure 4D–F illustrate how alpha activity shifts from posterior electrodes to frontal electrode with increasing propofol concentration (Cimenser et al., 2011; Purdon et al., 2013).

EEG recordings under different sleep stages

The HD-EEG data shown in Figure 5 is from a young adult subject (age 28, female) undergoing a sleep study, recorded using the ANT Neuro eego mylab system. After preprocessing (see ‘Materials and methods’), we visually selected 10 clean 5 s segments of data that exhibited strong alpha oscillations, each from relaxed wakeful (eyes closed) and rapid eye movement (REM) sleep to apply OCA.

Figure 5

Download asset Open asset

Oscillation component analysis (OCA) of the electroencephalogram (EEG) from a healthy young volunteer to compare between (A, C) wakeful resting state and (B, D) rapid eye movement (REM) sleep.

Panels (**A, B**) show the power spectral densities (PSDs) of reconstructed EEG activity within each canonical band. Panels (**C, D**) show the three dominant (in terms of sensor-wide power) alpha component: the topographic maps show the magnitude (left) and phase (right) distribution of sensor-level mixing, the time courses are 1 s representative example of the extracted oscillations from the selected epochs. The rightmost black bars display the coherency measure within alpha band. The contrasting topographic distribution of the alpha components (8–12 Hz), the shape of the oscillation power spectrum and alpha coherence hints at a distinct generating mechanism for alpha waves during stage 2 REM sleep compared to awake eyes closed alpha wave. The sensor weights are scaled to have a maximum value 1 so that the units of time series traces are in µV.

We fitted OCA models with 20, 25, 30, 35, 45, 50, 55, 60 oscillation components and selected the model structure with the highest negative variational free energy as the optimum OCA model. The empirical Bayes criterion selected 55 and 50 OCA components for relaxed wakeful (eyes closed) and REM sleep, respectively. Figure 5 follows a format similar to Figure 4: panels A nd B show the power spectrum density of the recorded EEG and the reconstructed sensor data from the OCA components grouped in slow, theta and alpha bands; while panels C and D visualize three of the extracted oscillation components within the alpha band, ordered according to their strengths. The topographical maps show the scalp distribution of these components, while the blue and red traces show 1-s-long oscillation time-series pair corresponding to each component.

Activity within the alpha band (8–13 Hz) could be summarized by a few oscillation components: for example, only seven and three OCA components had their center frequency within alpha band during relaxed wakeful (eyes closed) condition and REM sleep, respectively. These alpha oscillation components during relaxed wakeful (eyes closed) condition and REM sleep appeared to be distinct in terms of their temporal dynamics: relaxed wakeful alpha components are more regular (i.e., closer to sinusoidal activity) than REM alpha. But REM alpha coherency (0.4307) is moderately higher than that of relaxed wakeful alpha (0.3195). The spatial distribution of the REM alpha component appeared to be confined primarily within the posterior channels while the awake alpha is more distributed. To quantify the degree of similarity or difference between the spatial distributions, we calculated the principal angle between the spatial mixing maps of the awake eyes closed and REM alpha components. The principal angle measures the similarity between two subspaces, where an angle of 0° indicates that one subspace is a subset of the other and an angle of 90° indicates that at least one vector in a subspace is orthogonal to the other (Bjorck and Golub, 1973). The spatial mixing maps for awake eyes closed alpha and REM sleep had a principal angle of 44.31°, suggesting that the subspaces, and thus the underlying cortical generators, were substantially different.

Resting-state MEG recording

Finally, we demonstrate the utility of OCA on a MEG dataset from the Human Connectome Project.The MEG data used here is the resting-state recording from subject 104 012_MEG, session 3.

We fit OCA on 34 clean 5s epochs, and from among the OCA with 20, 25, 30, 35, 40, 45, 50, 55 components, empirical Bayes criterion selected the OCA with 30 components. Out of these, 13 were within the slow + delta band and 10 were alpha oscillations. For this example, we sought to highlight a different downstream analysis of the extracted oscillation components. We chose the three leading slow oscillation components and the three leading alpha oscillation components from the set of identified oscillation components.

For each oscillation, OCA extracts a pair of time traces, that is, one real time trace and one imaginary time trace (see Appendix 2, section ‘Oscillation states and analytic signals’). These ‘real’ and ‘imaginary’ indices can be utilized to compute the instantaneous phase and instantaneous amplitude of individual oscillation components at every timepoint, without having to rely on Hilbert transform (Wodeyar et al., 2021).

A_{t}^{(*)} = \sqrt{{x_{t, 1}^{(*)}}^{2} + {x_{t, 2}^{(*)}}^{2}}, P_{t}^{(*)} = \arctan (\frac{x_{2, t}^{(*)}}{x_{1, t}^{(*)}})

Further, we pick any one of the slow oscillation components and one of the alpha oscillation components and consider the ordered pair of slow component phase and alpha component amplitude at a single timepoint as a sample drawn from the joint distribution of the ‘phase–amplitude’ of these two components. These samples then can be used to quantify cross frequency phase–amplitude coupling between these components via either nonparametric (Cohen, 2008; Cheng et al., 2018; Martínez-Cancino et al., 2019) or parametric modeling (Soulat et al., 2022; Perley and Coleman, 2022).

Figure 6 demonstrates a possible way of investigating the coupling between the slow oscillation component phase and alpha oscillation component amplitude: the conditional mean of the alpha amplitude given slow/delta phase is higher at some specific phases between all three slow waves and OSC012, but overall flat for other alpha oscillation components. This example illustrates how the OCA algorithm can be used to identify specific oscillatory components that show phase–amplitude coupling.

Figure 6

Download asset Open asset

Cross–frequency phase–amplitude coupling in oscillation component analysis (OCA) components extracted from resting-state magnetoencephalogram (MEG) recording.

The black traces show the conditional mean of a selected alpha component (8–12 Hz) amplitude given another selected slow/delta component (0–4 Hz) phase. The three slow oscillations and three alpha oscillations that explained the highest variance were selected for demonstration purposes. The topographic maps show the magnitude (left) and phase (right) distribution of sensor-level mixing of the selected components.

Discussion

OCA is a novel approach to the multichannel component decomposition problem that can identify a number of independent spatio-temporal components, but does so in the context of underlying temporal dynamics specified by a state-space model. In the state-space representation used here, the dynamics of an elemental oscillation are parameterized by a central frequency, f, a damping parameter, a, and the second-order statistics of stochastic driving noise, $σ^{2}$ . Meanwhile, an associated spatial mixing pattern at the sensor level quantifies the contribution of the elemental oscillation to each observed channel.

The OCA learning algorithm uses an instance of the GEM algorithm to iteratively match the parameters of the oscillation state-space parameters to the second-order statistics of the M/EEG data. The goodness-of-fit for the data-driven matching procedure is defined within a Bayesian framework, as is the inference for the oscillation time courses, their sensor-level mixing patterns, and the measurement noise covariance matrix. This same goodness-of-fit metric is further used to determine the number of oscillations present in the multichannel data via empirical Bayes model selection. Once the number of oscillations, oscillation state-space parameters, their sensor-level mixing patterns, and the measurement noise covariance matrix are estimated, the elemental oscillatory activities can be extracted as a pair of time courses, that is, the ‘real’ and ‘imaginary’ parts of the analytical signal, from any given observation.

The oscillator state-space parameters discovered in OCA are akin to the ad hoc parameters, for example, peak frequency, full-width-half-maxima bandwidth, and peak oscillation power, etc., encountered in Fourier/wavelet-based frequency-domain nonparametric time-series analysis (Donoghue et al., 2020). However, OCA circumvents a major drawback of frequency-domain methods when searching for neural oscillations from multichannel recordings. Generally, the source-to-sensor mixing complicates the interpretation of the topographies generated by signal processing tools that treat each channel individually (Schaworonkow and Nikulin, 2022). In a nonparametric frequency-domain approach, the only way to discover this spatially correlated structure across sensors is to perform eigenvalue analysis on the cross-spectral density on a frequency-by-frequency basis, also known as global coherence analysis (Cimenser et al., 2011; Weiner et al., 2023; Mitra and Bokil, 2007). If the eigenvalue distribution exhibits sufficient skewness, the associated frequency is deemed coherent and the leading eigenvectors are identified as the ‘principal’ sensor networks at that frequency; in other words, a proxy for a strong oscillatory component. As the number of sensors increases, this analysis becomes increasingly intractable: in a sensor array with ∼10² sensors, the cross-spectral density matrix has ∼10⁴ entries. Alternatively, the state-space representation used here in OCA provides a convenient and compact way to perform multichannel time-domain modeling and analysis (Nise, 2011). The OCA modeling approach decouples the oscillatory dynamics and the spatial mixing pattern effectively by estimating one set of parameters for each discovered oscillation and the associated spatial mixing matrix, thereby avoiding the need for cross-spectral density matrices and frequency parameters for individual channels altogether (Gunasekaran et al., 2023). Another major advantage of OCA over global coherence analysis is that OCA automatically identifies the dominant coherent oscillations across the channels based on its probabilistic generative model and Bayesian learning approach.

The iterative parameter estimation procedure of OCA is clearly more computationally burdensome compared to conventional frequency-domain nonparametric methods. However, once the OCA parameters are estimated, OCA can decompose any given multichannel recording segments into the oscillation time-series pairs. In that sense, OCA parameter learning can be viewed as iterative estimation of a pair of optimal spatio–temporal filters (Anderson and Moore, 2005) for each oscillation, parameterized by the estimated oscillation state-space parameters and spatial maps. These optimal spatio-temporal filters are then applied to the multichannel data to extract the oscillation time courses. As a result, the extracted narrowband activity is endogenous to the data, rather than being imposed by an arbitrary narrowband filter (Yeung et al., 2007).

This behavior of OCA is similar to oscillatory component extraction methods based on data-driven spatial filtering, where the extracted component inherits the intrinsic oscillatory dynamics around the specified frequency (Nikulin et al., 2011; de Cheveigné and Arzounian, 2015; Cohen, 2017; de Cheveigné and Arzounian, 2015). In fact, OCA employs the same philosophy as spatio-temporal source separation (Cohen, 2018), where the extracted narrow-band signal and its time-shifted versions are again projected to a temporal coordinate space to further enhance the power within the given frequency band. However, spatio-temporal source separation establishes the temporal constraint via nonparametric sample correlation estimates that are sensitive to noise or artifacts and that require substantial amounts of data for estimation. Another important distinction is that spatio-temporal source separation determines the sensor weights first and then obtains the temporal filter kernel from the projection of the multichannel data. In contrast, OCA updates the sensor weights and the parameters for the temporal filter iteratively within the Bayesian formulation of the state-space framework, thus jointly estimating the spatial and temporal components simultaneously from the data.

The class of ICA-based methods, on the other hand, assumes that the component time courses share no mutual information, that is, are statistically independent. Without any explicit assumption on the temporal dynamics of the generative process, the properties of the identified components may be difficult to interpret and may also be unreliable. In fact, the identified ICA components typically require subjective visual inspection by experts for their possible interpretation. A number of investigators have recognized these limitations with ICA, in particular the assumption of temporal independence, and have proposed generalizations of the ICA framework to incorporate auto-regressive modeling of source dynamics (Pearlmutter and Parra, 1996; Parra, 1998). Despite demonstrating improved source separation performance in naturalistic signals like music (Pearlmutter and Parra, 1996), adoption of these methods for the analysis of neurophysiological data has been slow. This may be due in part to the lack of interpretability of the higher order auto-regressive model structures. Alternatively, Brookes et al., 2011 use a combination of Hilbert envelope computation within predefined frequency bands and temporal downsampling to identify meaningful temporally independent time signals, but stop short of trying to disambiguate different oscillatory sources. Thus, they only assess the connectivity pattern within predefined bands, that is, how different areas of the brain are harmonized through modulation of the oscillations or vice versa inside those predefined bands. The spatial maps recovered from anatomically projected resting-state MEG data by this method resemble spatial patterns of fMRI resting-state networks. OCA is close to the first variant of ICA, but describes the identified component sources using a 2D state-space vector that efficiently represents oscillatory dynamics (Matsuda and Komaki, 2017a) in a manner that is easy to interpret. Matsuda and Komaki, 2017b described a similar state-space model for multivariate time series, but did not apply the model to neurophysiological data, perhaps due to its high dimensionality. The OCA algorithm differs by way of its learning and model selection algorithms, which allow it to select the number of components to extract in a statistically principled and data-driven manner. Overall, OCA delivers the best of the previously mentioned eigen-analysis-based source separation methods and ICAs: a decomposition method that can identify independent components, but where each component represents an underlying oscillatory dynamical system. The generative link between each component and the underlying dynamical system makes interpretation of the components straightforward. Here, we focus specifically on analyzing neurophysiological data that are exemplified by highly structured oscillatory dynamics. But the methods we describe apply to more general, arbitrary dynamics that can be approximated by linear state-space models, and can be equally useful so long as the state-space model is interpretable.

There is another important distinction between OCA and the aforementioned BSS techniques. These methods directly estimate the sensor weights for a weighted linear combination of sensor recordings that produce the individual temporal components (Haufe et al., 2014). In other words, these BSS methods provide a backward decoding algorithm. The linear weights are not directly interpretable since they do not specify how the individual temporal components map to the observed data. That uniquely interpretable information is provided by a forward mapping which must be estimated in a separate step. When the number of linear components matches with the number of sensors, this forward mapping can be calculated via simple matrix inversion. However, when there is a reduced set of components the forward mapping can only be obtained by solving an inverse problem that requires knowledge of the source and sensor sample covariance matrix. In contrast, the OCA model involves only the (forward) spatial distribution matrix and our algorithm directly estimates it. In fact, the backward extraction model of OCA involves nonzero weights for time-shifted signals and is never explicitly computed. In that sense, OCA avoids an extraneous transformation step that could inadvertently introduce errors to the spatial mixing pattern estimates.

Another popular time-domain approach for oscillatory signal discovery is the multichannel extension of empirical mode decomposition. Empirical mode decomposition models the time-series data as linear combination of intrinsic oscillations called intrinsic mode functions (IMFs). Identification of IMFs depends critically on finding the local mean of the local extrema, which is not well-defined in the context of multichannel recordings. Rehman and Mandic, 2010 chose to take real-valued projections of the $n$ -channel data, with the projections taken along direction vectors uniformly sampled on a unit spherical surface in an $n$ -dimensional coordinate space. The local extrema of each of the projected signals are extrapolated to obtain a set of multivariate signal envelopes. Clearly, this random sampling in a high-dimensional space is computationally demanding, and the fully nonparametric formulation requires a substantial amount of data, making the procedure sensitive to noise. On the contrary, OCA employs a parametric model that represents oscillatory dynamics in a manner that is statistically efficient and resilient to noise.

Conclusion

In summary, starting from a simple probabilistic generative model of neural oscillations, OCA provides a novel data-driven approach for analyzing multichannel synchronization within underlying oscillatory modes that are easier to interpret than conventional frequency-wise, cross-channel coherence. The overall approach adds significantly to existing methodologies for spatio-temporal decomposition by adding a formal representation of dynamics to the underlying generative model. The application of OCA on simulated and real M/EEG data demonstrates its capabilities as a principled dimensionality reduction tool that simultaneously provides a parametric description of the underlying oscillations and their activation pattern over the sensor array.

Materials and methods

Key resources table

Reagent type (species) or resource	Designation	Source or reference	Additional information
Software, algorithm	MNE-python 1.2	https://mne.tools/stable/index.html; Gramfort et al., 2014
Software, algorithm	Eelbrain 0.37	https://eelbrain.readthedocs.io/en/stable/; Brodbeck et al., 2023
Software, algorithm	purdonlabmeeg	This paper; Das, 2024	Available at https://github.com/proloyd/purdonlabmeeg

M/EEG preprocessing

All prepossessing were performed using MNE-python 1.2 (Gramfort et al., 2014) and Eelbrain 0.37 (Brodbeck et al., 2023), with default setting of the respective functions.

EEG recording during propofol-induced unconsciousness

Request a detailed protocol

The EEG was bandpass filtered between 0.1 Hz to 40Hz, followed by downsampling to 100 Hz and average referencing prior to performing OCA. The 0.1 Hz highpass filtering was done to filter out slow drifts in the EEG recordings, which can adversely affect the OCA fitting procedure. All selected epochs had peak-to-peak signal amplitude less than 1000 mV. The prior on the noise-covariance was estimated from the same EEG recording, after high-pass filtering above 30 Hz.

Sleep EEG recording

Request a detailed protocol

The EEG data was first downsampled to 100 Hz after bandpass filtering within 1–40 Hz. The flat and noisy channels are first identified upon visual inspection. We computed neighborhood correlation for the rest of the channels and marked channels with median correlation, $ρ < 0.4$ as bad channels. These channels are dropped for the subsequent analysis.

MEG recording from the Human Connectome Project

Request a detailed protocol

We divided the entire resting-state MEG recording from subject 104 012_MEG, session 3 into 5 s epochs, and selected 34 epochs with peak-to-peak amplitude less than 4000 fT for OCA. Prior to fitting OCA, we used signal-space projection (Uusitalo and Ilmoniemi, 1997) to remove heartbeat and eye movement artifacts from MEG signals. The repaired MEG recordings are then downsampled to 100 Hz after applying appropriate anti-aliasing filter.

Code availability

Request a detailed protocol

An implementation of the state inference and parameter learning algorithms described in this article is available as part of purdonlabmeeg Python library on GitHub at https://github.com/proloyd/purdonlabmeeg, copy archived at Das, 2024 under MIT license. The simulation script to generate Figure 2, included in the repository, demonstrates how to run OCA on preprocessed EEG/MEEG data and visualize the results. Codes for analyzing experimental EEG and MEG recordings exactly follow the simulation script, i.e., is essentially a single function call with the experimental data (preprocessed according to previous section) as an argument and the results are to also be visualized as done in the simulation script. Thus they were not included as a separate repository, but can be made available to individuals upon request.

Appendix 1

Mechanistic origin

Here, we detail how the multichannel multi-oscillator state-space model can be understood in terms of a biophysical model of brain wave generation.

Charge-continuity and electromagnetic brain waves

Maxwell’s equation provides the following charge-continuity equation that governs the spontaneous spatial-temporal evolution of the potential, $Φ$ , in brain tissues:

\frac{\partial}{\partial t} \nabla^{2} Φ = - \nabla \cdot Σ \cdot Φ .

This differential equation is satisfied by traveling waves of the form, $Φ \sim \exp [- j (k \cdot r + Ω t)]$ , where $r$ and $t$ denote the spatial coordinate and time, respectively, with the wave number, $k$ , and the complex frequency, $Ω$ , satisfying following dispersion relation: $D (Ω, k) = j Ω | k |^{2} - Σ_{i j} k_{i} k_{j} - j \partial_{i} Σ_{i j} k_{j}$ , in tensor notation. Separating the real and imaginary parts of the complex frequency, $Ω = ω + j γ$ , we arrive at the following expressions for decaying and oscillatory frequencies:

γ = - \frac{Σ_{i j} k_{i} k_{j}}{| k |^{2}}, ω = \frac{\partial_{i} Σ_{i j} k_{j}}{| k |^{2}} .

Clearly, wave solutions $Φ$ satisfying $| γ / ω | ≪ 1$ , that is, the oscillatory frequency is much larger than the decaying frequency, can give rise to electromagnetic oscillations in brain parenchyma, governed by the brain composition and architecture. Detailed treatment of such a model, considering the anisotropy and inhomogeneity of the brain tissue, can be found in Galinsky and Frank, 2020.

Traveling wave solution to the state-space model

Here, we consider one such oscillatory electric potential being recorded at $L$ sensors placed at coordinates $r_{1}, r_{2}, \dots, r_{L}$ (w.r.t. a reference electrode placed at infinity). The potentials sampled (sampling frequency, $f_{s} = 1 / Δ t$ ) at those electrodes at time index $n$ are given as

Φ (l, n) = ℜ {A_{l} \exp [- j (k \cdot r_{l} + Ω n Δ t)]}

with $A_{l}$ representing the amplitude observed at lth electrode. The separability of the spatial part $ϕ^{s p a t i a l} (l) = A_{l} \exp [- j (k \cdot r_{l})]$ and the temporal part $Φ^{t e m p o r a l} (n) = \exp [- j (Ω n Δ t)]$ of the general solution is crucial for the following development since only the temporal part, $Φ^{t e m p o r a l} (n + 1)$ , evolves in time,

Φ^{t e m p o r a l} (n + 1) = \exp [j Ω Δ t] Φ^{t e m p o r a l} (n),

while the spatial part $ϕ^{s p a t i a l} (l)$ remains constant. Adapting the following real-valued matrix-vector notation that consists of the real and imaginary parts of the temporal and spatial components, $Φ^{t e m p o r a l} (n) = x_{n, 1} + j x_{n, 2}$ and $Φ^{s p a t i a l} (l) = m_{l, 1} + j m_{l, 2}$ , we arrive at:

[\begin{matrix} x_{n + 1, 1} \\ x_{n + 1, 2} \end{matrix}] = \exp [- γ Δ t] [\begin{matrix} \cos (ω Δ t) & - \sin (ω Δ t) \\ \sin (ω Δ t) & \cos (ω Δ t) \end{matrix}] [\begin{matrix} x_{n, 1} \\ x_{n, 2} \end{matrix}],

y_{n + 1}^{(l)} = [\begin{matrix} A_{l} m_{l, 1} & A_{l} m_{l, 2} \end{matrix}] [\begin{matrix} x_{n + 1, 1} \\ x_{n + 1, 2} \end{matrix}],

where $y_{n + 1}^{(l)}$ is the recorded potential at electrode $l$ at time index $(n + 1)$ . Essentially, the complex-valued oscillatory function, $Φ^{t e m p o r a l} (n)$ , is an analytic signal whose real and imaginary parts constitute the ordered pair, $x_{n} = [x_{n, 1}, x_{n, 2}]$ , which we call the oscillation state. Equation A1.5 describes a circular motion of the oscillation state $x_{n}$ in a 2D state space with frequency $ω$ at every time step through the damped rotation matrix $\exp [- γ Δ t] R (ω Δ t)$ , while a projection, dictated by the electrode location (Equation A1.6), is recorded at the EEG electrodes as an oscillation (see Figure 1a). This deterministic model forms the basis for the probabilistic linear Gaussian state-space model in Equation 1 which introduces process and observation noise, $υ_{n}$ and $ϵ_{n}$ , respectively, with redefined spatial components as $c_{l, *} = A_{l} m_{l, *}$ , and the decay term $\exp [- γ Δ t]$ as a damping parameter, $a$ . We also represented the frequency of oscillation, f, as f = ω/2π by replacing $ω Δ t = 2 π f / f_{s}$ .

Appendix 2

Oscillation states and analytic signals

Here, we consider real-valued observations, $y_{n} \in R$ , admitting a state-space representation in complex plane, that is, generated from complex-valued states, $z_{n} \in C$ :

\begin{matrix} z_{n} = ρ x_{n - 1} + q_{n} \\ y_{n} = \frac{1}{2} (z_{n} + {\bar{z}}_{n}) + r_{n} \end{matrix}

where ${\bar{z}}_{n}$ denotes complex conjugate of ${\bar{z}}_{n}$ , q_n and r_n are complex-valued and real-valued random variables. In steady state, that is, $t \to \infty$ , the states, z_n, admit Fourier representation given by

Z (e^{- j ω}) = \frac{Q (e^{- j ω})}{1 - ρ e^{- j ω}}

It is clear that if q_n is an analytic signal, that is, $Q (e^{- j ω})$ has no negative frequency component, the same will be true for $Z (e^{- j ω})$ .

Further, since q_n is an analytic signal, there exists a real-valued signal s_n, such that $q_{n} = s_{n} + j (h * s)_{n}$ , where discrete time Fourier transform of h_n is given by

H (e^{- j ω}) = {\begin{cases} - j & for 0 < ω \leq π \\ 0 & for ω = 0 \\ j & for - π < ω < 0 \end{cases} .

Considering s_n to be a white noise with variance $σ^{2}$ , it is easy to verify that $h_{n} * s_{n}$ is also white noise with variance $σ^{2}$ and uncorrelated to s_n. This makes q_n to be a ‘analytic’ white noise, that is power spectrum density of q_n is flat for $0 > ω > π$ , and zero in $- π > ω > 0$ . We can then derive the following state-space representation, with real valued but 2D states as in Equation 1, but with additional constraint that $υ_{n, 1}$ and $υ_{n, 2}$ are related by $υ_{n, 2} = (h * υ_{1})_{n}$ . But since such constraint is not easy to enforce during the inference or model learning from experimental data, so we relaxed the constraint on the state-noise covariance to be the one stated in Equation 1.

Extraction of instantaneous amplitude and phase

However, we can still compute the instantaneous amplitude and phase of the signal from this state-space representation from the following observation. The oscillation states assume their successive values in a way that traces limit cycle-like maps around origin (see Figure 1a). This allows us to define the instantaneous phase as the angle the oscillation state makes from a fixed reference direction at any given timepoint (Rosenblum et al., 1997, section A.1 and 2). For simplicity, we consider the positive direction on the real line as the reference. The amplitude is simply given by the distance from the origin.

A_{t}^{(*)} = \sqrt{{x_{t, 1}^{(*)}}^{2} + {x_{t, 2}^{(*)}}^{2}} P_{t}^{(*)} = \arctan (\frac{x_{2, t}^{(*)}}{x_{1, t}^{(*)}})

In this sense, the oscillator state-space representation provides a generalization of phasor concept, similar to analytic signals (i.e., by allowing for time-varying amplitude, phase, and frequency, in contrast to invariant amplitude, phase, and frequency of phasor) but with relatively relaxed constraint on the relation between real and imaginary parts.

Appendix 3

Model parameter estimation and model selection

Negative variational free energy

For an L-channel M/EEG recording $y_{t}, t = 1, 2, \dots, T$ , $M$ -oscillator probabilistic state-space oscillator model admits to the following distribution:

\begin{array}{ll} p ({x_{t}, y_{t}} | C, R, F, Q, M) = \prod_{t = 1}^{T} \frac{1}{\sqrt{| (2 π) R |}} \exp - \frac{1}{2} {‖ y_{t} - C x_{t} ‖}_{R^{- 1}}^{2} \\ \frac{1}{\sqrt{| (2 π) Q |}} \exp - \frac{1}{2} {‖ x_{t} - F x_{t - 1} ‖}_{Q^{- 1}}^{2} \end{array}

Since $F$ and $Q$ are parameterized using hyperparameter $θ$ , we will re-parameterize left-hand side of Equation A3.1 as $p ({x_{t}, y_{t}} ∣ C, R, θ, M)$ . To simplify our exposition, we define $X = [x_{1}, x_{2}, \dots, x_{N}]$ and use $X$ and ${x_{t}}$ interchangeably (similarly for $Y$ ). We also extensively use $vec (\circ)$ notation and associated terminology introduced in Muirhead, 1982.

We note that computation of ensemble likelihood of the presented model requires marginalizing over $X, C, R$ , and $M$ given the priors in Equations 4 and 6:

\begin{array}{ll} p ({y_{t}}) & = \sum_{M = 1}^{M_{m a x}} \int p (X, Y, C, R | α, θ, M) p (M) d {X, C, R} \\ = \sum_{M = 1}^{M_{m a x}} \int p (X, Y | C, R, α, θ, M) p (C ∣ α, M) p (R) p (M) d {X, C, R} \end{array}

However, this involves an intractable integration that cannot be performed analytically. We avoid this intractability by using a lower bound to the ensemble likelihood as a surrogate for the same. Specifically, we invoke Neal–Hinton representation theorem (Neal and Hinton, 1998) to obtain a lower bound on the ensemble log-likelihood as

\begin{array}{ll} L & = \log p ({y_{t}}) \\ = \log \sum_{M = 1}^{M_{m a x}} \int p (X, Y, C, R ∣ α, θ, M) p (M) d {X, C, R} \\ = \log \sum_{M = 1}^{M_{m a x}} \int \frac{p (X, Y, C, R ∣ α, θ, M) p (M)}{q (X, C, R ∣ M) q (M)} q (X, C, R ∣ M) q (M) d {X, C, R} \\ \geq \sum_{M = 1}^{M_{m a x}} q (M) \int \log \frac{p (X, Y, C, R ∣ α, θ, M) p (M)}{q (X, C, R, ∣ M) q (M)} q (X, C, R ∣ M) d {X, C, R} \\ = \sum_{M = 1}^{M_{m a x}} q (M) [\log \frac{p (M)}{q (M)} + \int \log \frac{p (X, Y, C, R, α, θ ∣ M)}{q (X, C, R ∣ M)} q (X, C, R ∣ M) d {X, C, R}] \\ = {⟨ \log \frac{p (M)}{q (M)} + {⟨ \log \frac{p (X, Y, C, R, α, θ ∣ M)}{q (X, C, R ∣ M)} ⟩}_{q (X, C, R ∣ M)} ⟩}_{q (M)} := F, \end{array}

which holds for any arbitrary conditional distribution $q$ and can be computed analytically if $q$ is chosen carefully. This lower bound is known as negative variational free energy (Quinn and Šmídl, 2006) and attains the actual ensemble log-likelihood only when $q$ is the exact Bayes posterior distribution (Attias, 1999). In the last line, we use ${⟨ \circ ⟩}_{q (X, C, R ∣ M)}$ to denote the average of the expression within $⟨ ⟩$ with respect to the model posterior $q (X, C, R ∣ M)$ .

Since negative variational free energy as a lower bound on the ensemble log-likelihood, maximization of negative variational free energy will tend to maximize of the ensemble log-likelihood. Even though global maximization of the ensemble log-likelihood is not guaranteed, the ensemble log-likelihood at the point of maximum negative variational free energy is guaranteed to at least greater than the negative variational free energy. On the other hand, the negative variational free energy can be shown to be proportional to the negative Kullback–Leibler divergence between the exact posterior and approximate posterior (Quinn and Šmídl, 2006). So, the maximization shall result in closer approximation exact posterior. In the next section, we work with the intuition that maximization of negative variational free energy leads to better approximation to the exact posterior and increment of ensemble log-likelihood.

Variational Bayes inference

However, the fact that the ensemble log-likelihood cannot be found analytically makes computation of the exact Bayes posterior distribution intractable. We thus employ an efficient VB inference procedure that approximates the exact Bayes posterior with a distribution of form:

q ({x_{t}}, C, R ∣ M) = q ({x_{t}} ∣ M) q (C ∣ M) q (R ∣ M),

that is, that factorizes over ${x_{t}}, C, R$ , conditional on the model structure, $M$ . This particular choice of functional constraint leads to tractable conditional distributions when the trailing term of Equation A3.2 is subjected to maximization w.r.t. functional form:

q (X ∣ M) = \frac{1}{\sqrt{| (2 π) Σ_{x} |}} \exp - \frac{1}{2} {‖ v e c (X) - μ_{x} ‖}_{Σ_{x}^{- 1}}^{2}

q (C ∣ M) = \frac{1}{\sqrt{| (2 π) Σ_{c} |}} \exp - \frac{1}{2} {‖ v e c (C) - μ_{c} ‖}_{Σ_{c}^{- 1}}^{2}

q (R ∣ M) = γ_{ρ, L} {| Ω |}^{ρ / 2} {| R |}^{- (ρ + L + 1) / 2} \exp - \frac{1}{2} T r {Ω R^{- 1}}

Here, we omit the trivial derivation of the functional form for brevity of the presentation; curious readers are encouraged to expand the trailing term of Equation A3.2 and follow Attias, 1999, section 2.4 as an exercise to obtain the optimal functional forms.

Once the abovementioned optimal functional forms of the approximate posterior distribution have been identified, the next task is to find the optimal quantities that parameterize those functional forms. We show that the parameters of any of these distribution can be expressed in terms of the M/EEG data and moments of other two distributions. Generally speaking, the expressions can be found by collecting relevant terms in the expansion of trailing term of Equation A3.2, followed by simple algebraic manipulation. For example, the parameters of $q (R ∣ M)$ are given by

Ω = Ψ + \sum_{t = 1}^{T} {⟨ y_{t} - C x_{t} ⟩}_{q (X, C ∣ M)} a n d ρ = ν + L .

Expression for the parameters of the Gaussian distribution is a bit involved if one attempts collecting terms and try completing squares. Instead, we use the following facts about Gaussian log-posteriors: (1) since the mode and mean of a multivariate Gaussian distribution coincide, the mean of Gaussian distributions can be found by maximizing the approximate log-posterior, and (2) the negative Hessian of the log-posterior is the inverse of covariance matrix (Fahrmeir, 1992). It is worth mentioning here that the log-posterior maximization can be carried out in closed form by solving a high-dimensional system of linear equations (first-order optimality criterion, Boyd and Vandenberghe, 2004), which involves Hessian matrix inversion.

In case of the latent oscillator states, q( $X ∣ M)$ , the expressions of inverse covariance matrix and the mean vector is as follows:

\begin{matrix} Σ_{x}^{- 1} = [\begin{matrix} D & S & 0 & 0 & \dots & 0 & 0 \\ S^{⊤} & D & S & 0 & \dots & 0 & 0 \\ 0 & S^{⊤} & D & S & \dots & 0 & 0 \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ 0 & 0 & 0 & 0 & \dots & S^{⊤} & D \end{matrix}] & \begin{matrix} where \\ D = {⟨ C^{⊤} R^{- 1} C ⟩}_{q (C, R ∣ M)} \\ + {\hat{F}}^{⊤} {\hat{Q}}^{- 1} \hat{F} + {\hat{Q}}^{- 1} \\ S = - {\hat{F}}^{⊤} {\hat{Q}}^{- 1} \end{matrix} \\ μ_{x} = Σ_{x} \times v e c (B y_{t}) & B = {⟨ C^{⊤} R^{- 1} ⟩}_{q (C, R ∣ M)} \end{matrix}

We exploit the block tridiagonal structure of the hessian matrix to realize a fast and stable inverse (Asif and Moura, 2000; Jain et al., 2007). This facilitates fast computation of the mean and covariance matrix of the approximate log-posterior.

Before tackling the case of the mixing matrix, $C$ , we note that each latent oscillator state is a pair of two coordinates and these two oscillator coordinates correspond to two consecutive columns of the mixing matrix, that is, contribution of any given latent oscillator state on the observation can be viewed as the inner product, that is, $y_{t}^{(l)} = c_{l, m} x_{t}^{(m)}$ , so that the column pairs $(C_{:, 2 l - 1 : 2 l})$ and latent oscillator states $(X_{2 l - 1 : 2 l, :})$ are only unique up to a scaling and a rotation. We resolve this ambiguity by simply fixing the first row of the mixing matrix column pair to $[\begin{matrix} 1 & 0 \end{matrix}]$ for all the oscillators. With this modification, the parameters of $q (C ∣ M)$ are given as

\begin{array}{ll} Σ_{c}^{- 1} & = Λ_{2 : L, 2 : L} \otimes T P_{x x} + α I_{2 (L - 1) M \times 2 (L - 1) M} \\ μ_{c} & = Σ_{c} (T vec ({(Λ_{2 : N} P_{y x})}^{⊤}) - λ_{1} \otimes (P_{x x} m_{1}^{⊤})) \end{array}

where

Λ = {⟨ R^{- 1} ⟩}_{q (R ∣ M)}, P_{x x} = \frac{1}{T} \sum_{t = 1}^{T} {⟨ x_{t} x_{t}^{⊤} ⟩}_{q (X ∣ M)}, P_{y x} = \frac{1}{T} \sum_{t = 1}^{T} {⟨ y_{t} x_{t}^{⊤} ⟩}_{q (X ∣ M)} .

Using the properties of Kronecker products (Muirhead, 1982), the inversion of the negative Hessian matrix and solution of the system of the linear equations can be carried out in a fast and efficient manner.

Evidently, this interconnectedness of the distributional parameters of the marginals of ${x_{t}}, C, R$ , given model structure $M$ , is a testament to the inter-dependence of these three variables. Since joint posterior distribution of is assumed to factorize over the marginals, that is, ${x_{t}}, C, R$ , given model structure $M$ and observations $Y$ are assumed to be independent, inter-dependence has been captured by the interrelated distributional parameters. In fact, the interrelated parameters pass information among them using the sufficient statistic (SS) of these distributions that are required to compute the averages in the update equations. The VB approach for posterior inference therefore leads to an iterative scheme for each model structure, $M$ , when the hyperparameters $(θ, α)$ are specified: we cyclically update the marginal posterior parameters of ${x_{t}}, C, R$ according to Equations A3.8, A3.9, and A3.7, , respectively. Each time a distribution (i.e., its parameters) is updated, the SS associated with the distribution is also recomputed to pass to the next update. The iterations are repeated until a stopping criteria is met: stopping criteria could be a fixed number of iterations or when the objective, that is, variational free energy stabilizes. The marginal posterior parameters obtained at the last iteration are the optimal distributional parameters providing the VB inference.

Generalized EM algorithm

The VB inference of form presented in Equation A3.3 still requires the hyperparameters $θ$ and $α$ to be specified. The hyperparameters that maximize the negative variational free energy are the second best choice for the hyperparameters, right after the ones that maximize the ensemble likelihood. Since the direct maximization w.r.t. the hyperparameters is impractical, we use an generalized version of expectation maximization (EM) algorithm that uses the variational posterior inference instead of the exact posterior to compute the expectation of the complete log-likelihood expression (Dempster et al., 1977; Fahrmeir, 1992; Shumway and Stoffer, 1982) to learn these hyperparameters from the M/EEG recording. One requires the following averages computed w.r.t. the approximate posterior $q (X ∣ M)$ ,

\begin{matrix} U^{(m)} = \sum_{t = 1}^{T} {⟨ x_{t - 1}^{(m)} {x_{t - 1}^{(m)}}^{⊤} ⟩}_{q (X ∣ M)}, \\ V^{(m)} = \sum_{t = 1}^{T} {⟨ x_{t}^{(m)} {x_{t - 1}^{(m)}}^{⊤} ⟩}_{q (X ∣ M)}, \\ W^{(m)} = \sum_{t = 1}^{T} {⟨ x_{t}^{(m)} {x_{t}^{(m)}}^{⊤} ⟩}_{q (X ∣ M)}, \end{matrix}

to update the hyperparameters $θ^{(m)} = (f^{(m)}, a^{(m)}, {(σ^{2})}^{(m)})$ as follows:

f^{(m)} = \frac{1}{2 π} \arctan {\frac{r t {V^{(m)}}}{t r {V^{(m)}}}},

a^{(m)} = \frac{\sqrt{tr {V^{(m)}}^{2} + rt {V^{(m)}}^{2}}}{t r {U^{(m)}}}

{σ^{2}}^{(m)} = \frac{1}{2 π} (tr {W^{(m)}} - {a^{(m)}}^{2} tr {U^{(m)}})

On the other hand, α update takes a relatively simple form:

α = \frac{2 (L - 1) M}{{⟨ {‖ C_{2 : L, :} ‖}^{2} ⟩}_{q (C ∣ M)}} .

In brief, we start with an initial set of the hyperparameters θ and α, and keep alternating between variational posterior inference (E-step) and hyperparameter update (M-step) until the hyperparameters or the variational free energy stabilize. Once the suitable hyperparameters are obtained, the variational inference is run for the last time, providing us with the final oscillation components.

Evaluation of variational free energy

Considering update rules of Equations A3.7, A3.13, , and Equation A3.10 we can simplify the second term in Equation A3.2 as

\begin{aligned} {⟨ \log \frac{p (X, Y, C, R, α, θ ∣ M)}{q (X, C, R ∣ M)} ⟩}_{q (X, C, R ∣ M)} = - \frac{(L + 2 M) T}{2} \log 2 π \\ - \frac{T}{2} \log | Q | - \frac{2 M (L - 1)}{2} - ⟨ \log q {(X ∣ M)}_{q (X ∣ M)} ⟩ \\ + \frac{2 M (L - 1)}{2} \log \frac{α}{2 π} - \frac{2 M T}{2} - ⟨ \log q (M ∣ m) ⟩ \\ (2) & + \frac{ν}{2} \log | Ψ | + γ_{ν, n} - (\frac{ν + T}{2} \log | Ω | + γ_{ν + T, n}) . \end{aligned}

After noting the following property of Gaussian random variable, $Z \sim N_{p} (μ_{z}, Σ_{z})$ :

{⟨ \log p (Z) ⟩}_{p (Z)} = - \frac{1}{2} \log | (2 π) Σ_{z} | - \frac{p}{2}

Equation A3.14 can be further simplified as:

\begin{aligned} {⟨ \log \frac{p (X, Y, C, R, α, θ ∣ M)}{q (X, C, R ∣ M)} ⟩}_{q (X, C, R ∣ M)} & = - \frac{L T}{2} \log 2 π - \frac{T}{2} \log | Q | + \frac{1}{2} \log | Σ_{x} | \\ (3) & + \frac{1}{2} \log | α Σ_{c} | + \frac{ν}{2} \log | Ψ | & + γ_{ν, n} - (\frac{ν + T}{2} \log | Ω | + γ_{ν + T, n}) \end{aligned}

We note here that the log-determinant of $Σ_{x}$ could be computed while solving Equation A3.8 exploiting the tridiagonal structure of ${Σ_{x}}^{- 1}$ . Similarly, the log-determinant of $Σ_{c}$ could be computed while solving Equation A3.9 using the properties of Kronecker products. Matrices $Φ$ and $Ω$ are in the sensor space dimension, thus log-determinant computation is relatively cheaper. The remaining matrix $Q$ is a diagonal matrix, thus its log-determinant can be computed easily. As a whole, the variational free energy for model structure $M$ can be computed as a by-product of the iterative updates.

Initialization of parameters and hyperparameters

Because the ensemble log-likelihood (hence the negative variational free energy) is not concave in general, initial values for the parameters and hyperparameters must be carefully chosen. We use the interpretation of the MVAR processes as a superimposition of multiple oscillatory components (Matsuda and Komaki, 2017b; Quinn et al., 2021; Neumaier and Schneider, 2001).

First of all, we estimate $Ψ$ using the fact that for single-channel noisy recording, $y^{(l)}$ where $S_{y^{(l)} y^{(l)}} (f)$ denotes the PSD of the signal:

10 \log_{10} \frac{R}{f_{s}} = lim_{f \to \infty} S_{y^{(l)} y^{(l)}} (f) .

We generalize Equation A3.16 to multichannel recording by high-pass filtering the M/EEG recording at $f_{s} - Δ f$ , ( $Δ f > 0$ ), and estimating inverse-Wishart parameters from the sample covariance of the high-pass filtered data, ${\tilde{y}}_{t}$ . We set $ν = T$ , and $Ψ = \sum_{t = 1}^{T} x_{t} x_{t}^{⊤}$ . We also initialize the noise covariance matrix as $R = \sum_{t = 1}^{T} x_{t} x_{t}^{⊤} / T$ .

We then fit an MVAR model of order $p$ on the multichannel data. Order $p$ is chosen according to Akaike information criteria (Cavanaugh and Neath, 2019). We perform eigendecomposition of the companion form of the MVAR parameter matrix (Neumaier and Schneider, 2001) to yield $λ$ , $V$ , and $W$ as eigenvalues, right and left eigenvectors. We choose the eigenvalues with $ℑ {λ^{(m)}} > 0$ and corresponding right eigenvectors $w^{(m)}$ and collect them as $\hat{λ}$ and $\hat{W}$ (assume there are $M$ such eigenvalues). The frequency and damping parameters are initialized as

f^{(m)} = arg (λ^{(m)}) \frac{f_{s}}{2 π} a^{(m)} = | λ^{(m)} | .

For ${σ^{(1)}}^{2}, {σ^{(m)}}^{2}, \dots, {σ^{(M)}}^{2}$ , we use the linear combination of the theoretical PSD of the oscillatory components that approximate the multitaper PSD (Babadi and Brown, 2014) of the leading time series the best. In short, we compute the multitaper PSD of the leading time series and the theoretical PSDs of the oscillatory components (Soulat et al., 2022, Supplementary Information) on the same frequency grid of spacing $f_{s} / (2 K)$ . The multitaper PSD, $ρ$ , is a $1 \times K$ vector, whereas the theoretical PSDs, $ϕ^{(m)}$ , form a $M \times K$ matrix $Φ$ . We initialize the ${σ^{(1)}}^{2}, {σ^{(m)}}^{2}, \dots, {σ^{(M)}}^{2}$ as the solution of the following linear equation:

[\begin{matrix} {σ^{(1)}}^{2} & {σ^{(m)}}^{2} & \dots & {σ^{(M)}}^{2} \end{matrix}] Φ = ρ .

In order to select a subset of the discovered M oscillatory components, we sort the oscillatory components (and columns of the $\hat{W}$ ) according to their theoretical contribution (most to least):

{σ^{2}}^{(m)} ‖ {\hat{w}}^{(m)} / {\hat{w}}_{1}^{(m)} ‖_{2} \sum_{f} ϕ^{(m)} (f),

and select the leading oscillatory parameters (and columns of the $\hat{W})$ .

With Equation A3.17, the spatial mixing components for the leading time series become $c_{l, m} = [1 0]$ . The rest of the elements of spatial mixing matrix, $C$ , are initialized from $\hat{W}$ as

c_{l, m} = [\begin{matrix} ℜ (\frac{{\hat{w}}^{(m)}_{l}}{{\hat{w}}^{(m)}_{1}}) & - ℑ (\frac{{\hat{w}}^{(m)}_{l}}{{\hat{w}}^{(m)}_{1}}) \end{matrix}]

Lastly, α is initialized by $2 (L - 1) M / {‖ C_{2 : L} ‖}_{F r o}$ with the abovementioned $C$ .

Model structure posterior

Finally, it is easy to verify that the following discrete distribution maximizes Equation A3.2:

q (M) = \frac{p (M) \exp - {⟨ \log \frac{p (X, Y, C, R, α, θ ∣ M)}{q (X, C, R ∣ M)} ⟩}_{q (X, C, R ∣ M)}}{\sum_{M^{'} = 1}^{M_{m a x}} p (M^{'}) \exp - {⟨ \log \frac{p (X, Y, C, R, α, θ ∣ M^{'})}{q (X, C, R ∣ M^{'})} ⟩}_{q (X, C, R ∣ M^{'})}},

which can be computed readily once the final variational inferences are made for $M = 1, \dots, M_{max}$ . The negative variational free energy of the modeling framework is given by

F = \sum_{M^{'} = 1}^{M_{m a x}} p (M^{'}) \exp - {⟨ \log \frac{p (X, Y, C, R, α, θ ∣ M^{'})}{q (X, C, R ∣ M^{'})} ⟩}_{q (X, C, R ∣ M^{'})} .

Empirical Bayes inference and model selection

Since the model order is not known, we estimate the model order, $M$ , by maximizing the approximate posterior distribution of model order, $q (M)$ , and treat $M$ as if it were known to be equal to this estimate. The same procedure is followed for the hyperparameters, for which we set them to the values that maximize the negative variational free energy. This way of determining model order or estimating unknown parameters comes under the umbrella of empirical Bayes framework (Robbins, 1964; Efron and Morris, 1973; Efron and Morris, 1975; Morris, 1983). Briefly, it provides a first-order approximation to the posterior mean of these quantites and neglects the uncertainty (i.e., second-order statistics) of the estimated quantities.

Appendix 4

Comparison of OCA to traditional approaches in experimental EEG data

We compare OCA with an oscillation finding approach based on channel-wise PSD visualization and ICA in a subset of real human EEG recording, taken from Figure 4 to provide empirical justification for OCA. For Appendix 4—figure 1, we used a 3.5-s-long EEG segment during the maintenance of 2 µg effect-site concentration of propofol. We expect strong frontal alpha activity and overall increase in slow activity as per existing literature (Purdon et al., 2013; Cimenser et al., 2011).

The frequency-domain approach examined here is based on the popular multitaper method (de Cheveigné and Parra, 2014): given a bandwidth parameter, power spectrum density is computed for individual channels and a butterfly plot is made (see Appendix 4—figure 1A, left panels). The oscillatory alpha activities are then identified visually (red overlay) as the peaks in these power spectrum plots, and their scalp distribution are obtained by averaging channel-wise power within given frequency bands around the peaks as shown in Appendix 4—figure 1A, right panels. We chose three different bandwidth, 1 Hz, 2 Hz, and 5 Hz, to demonstrate the inconsistency of frequency domain methods in ‘producing’ peaks in the line plots of the spectrum. Spectrum with 1 Hz bandwidth (top row) exhibits a lot of spurious peaks, whereas spectrum with 5 Hz bandwidth (bottom row) exhibits wide peaks, possibly encompassing multiple peaks. This inconsistency of ad hoc visual identification of peaks makes it subjective. Similarly, the topographic plots in Appendix 4—figure 1A demonstrates similar inconsistency: as estimation bandwidth increases (i.e., spectral resolution decreases), the topographic maps becomes less susceptible to spurious power leakage until the estimation bandwidth matches the bandwidth of underlying processes (compare middle row to top row). As the estimation bandwidth further increases, a severe spectral leakage renders the topographic maps largely uninformative of underlying separable oscillatory sources (bottom row).

ICA decomposition on this multichannel recordings identifies components that are neither slow oscillation nor alpha oscillation but combination of multiple oscillations as one would expect based on the observation that cumulative histograms of EEG recordings is approximately Gaussian (Brookes et al., 2011). We use ‘extended infomax’ (Lee et al., 1999) implementation provided by MNE-python 1.2 (Gramfort et al., 2014), the leading four ICA components are shown in Appendix 4—figure 1B.

Lastly, OCA is able to identify several oscillatory components (in both slow and alpha band) from this multichannel data, and extract their individual time series and topographic distribution. OCA is able to do so due to its explicit modeling of temporal dynamics in the form of the structured state-space representation. Appendix 4—figure 1B shows leading four components whose estimated center frequencies are within alpha band (9–13 Hz). This also demonstrates that OCA is bandwidth-free, that is, it adjusts estimation bandwidth to match the bandwidths of the underlying oscillations in a data-driven way as evident from the spectrum plots of the estimated oscillation time- courses.

Appendix 4—figure 1

Download asset Open asset

Empirical justification for oscillation component analysis (OCA) in analyzing real data.

A 3.5-s-long electroencephalogram (EEG) recording during propofol-induced anesthesia (effect-site concentration of 2 μg) is considered for this demonstration. (A) Power spectral density using multitaper method (for varying time-bandwidth product) of the EEG recording and power distribution over the EEG sensors in the marked frequency bands (red overlay) around visually identifiable peaks. (B) Four leading independent component analysis (ICA) components (left, middle, and right columns show topographic maps, power spectrum density, and time courses, respectively). (C) Four leading OCA components within alpha band (the topographic maps show the magnitude [left] and phase [right], while line plots show power spectrum density [left] and time-courses [right], respectively).

Data availability

The current manuscript is a computational study, so no data have been generated for this manuscript. An implementation of the state inference and parameter learning algorithms described in this article is available as part of purdonlabmeeg Python library on GitHub at https://github.com/proloyd/purdonlabmeeg (copy archived at Das, 2024) under MIT license. The simulation script to generate Figure 2, included in the repository, demonstrates how to run OCA on preprocessed EEG/MEEG data and visualize the results. Codes for analyzing experimental EEG and MEG recordings exactly follow the simulation script, i.e., is essentially a single function call with the experimental data (preprocessed according to previous section) as an argument and the results are to also be visualized as done in the simulation script. Thus they were not included as a separate repository, but can be made available to individuals upon request.

The following previously published data sets were used

1. Van Essen DC
2. Smith SM
3. Barch DM
4. Behrens TEJ
5. Yacoub E
6. Ugurbil K
(2013) ConnectomeDB
ID 104012MEG. MEG data in Section Resting state MEG recording.

https://db.humanconnectome.org

References

Book
1. Anderson BDO
2. Moore JB
(2005)
Optimal Filtering: Dover Books on Engineering

Dover Publications.
- Google Scholar
Conference
1. Asif A
2. Moura JMF
(2000) Inversion of block matrices with block banded inverses: application to Kalman-Bucy filtering
2000 International Conference on Acoustics, Speech and Signal Processing. pp. 608–611.

https://doi.org/10.1109/ICASSP.2000.862055
- Google Scholar
Conference
1. Attias H
(1999)
Inferring parameters and structure of latent variable models by variational bayes

Proceedings of the Fifteenth Conference on Uncertainity in Artificial Intelligence. pp. 21–30.
- Google Scholar
1. Babadi B
2. Brown EN
(2014) A review of multitaper spectral analysis
IEEE Transactions on Bio-Medical Engineering 61:1555–1564.

https://doi.org/10.1109/TBME.2014.2311996
- PubMed
- Google Scholar
Conference
(2018) State space oscillator models for neural data analysis
Annual International Conference of the IEEE Engineering in Medicine and Biology Society. pp. 4740–4743.

https://doi.org/10.1109/EMBC.2018.8513215
- Google Scholar
Preprint
1. Beck AM
2. He M
3. Gutierrez RG
4. Hotan GC
5. Purdon PL
(2022) An iterative search algorithm to identify oscillatory dynamics in neurophysiological time series
bioRxiv.

https://doi.org/10.1101/2022.10.30.514422
- Google Scholar
1. Bell AJ
2. Sejnowski TJ
(1995) An information-maximization approach to blind separation and blind deconvolution
Neural Computation 7:1129–1159.

https://doi.org/10.1162/neco.1995.7.6.1129
- PubMed
- Google Scholar
1. Bjorck A
2. Golub GH
(1973) Numerical methods for computing angles between linear subspaces
Mathematics of computation 27:579.

https://doi.org/10.2307/2005662
- Google Scholar
Book
1. Boyd SP
2. Vandenberghe L
(2004) Convex optimization
Cambridge University Press.

https://doi.org/10.1017/CBO9780511804441
- Google Scholar
Book
1. Bracewell RN
(2000)
McGraw-hill series in electrical and computer engineering

In: Bracewell RN, editors. The Fourier Transform and Its Applications. McGraw Hill. pp. 361–363.
- Google Scholar
1. Brodbeck C
2. Das P
3. Gillis M
4. Kulasingham JP
5. Bhattasali S
6. Gaston P
7. Resnik P
8. Simon JZ
(2023) Eelbrain, a Python toolkit for time-continuous analysis with temporal response functions
eLife 12:e85012.

https://doi.org/10.7554/eLife.85012
- PubMed
- Google Scholar
1. Brookes MJ
2. Woolrich M
3. Luckhoo H
4. Price D
5. Hale JR
6. Stephenson MC
7. Barnes GR
8. Smith SM
9. Morris PG
(2011) Investigating the electrophysiological basis of resting state networks using magnetoencephalography
PNAS 108:16783–16788.

https://doi.org/10.1073/pnas.1112685108
- PubMed
- Google Scholar
1. Buzsáki G
2. Draguhn A
(2004) Neuronal oscillations in cortical networks
Science 304:1926–1929.

https://doi.org/10.1126/science.1099745
- PubMed
- Google Scholar
Book
1. Buzsáki G
(2006) Rhythms of the Brain
Oxford University Press.

https://doi.org/10.1093/acprof:oso/9780195301069.001.0001
- Google Scholar
1. Cardoso JF
(1999) High-order contrasts for independent component analysis
Neural Computation 11:157–192.

https://doi.org/10.1162/089976699300016863
- PubMed
- Google Scholar
1. Cavanaugh JE
2. Neath AA
(2019) The Akaike information criterion: Background, derivation, properties, application, interpretation, and refinements
WIREs Computational Statistics 11:1460.

https://doi.org/10.1002/wics.1460
- Google Scholar
1. Cheng N
2. Li Q
3. Wang S
4. Wang R
5. Zhang T
(2018) Permutation mutual information: a novel approach for measuring neuronal phase-amplitude coupling
Brain Topography 31:186–201.

https://doi.org/10.1007/s10548-017-0599-2
- PubMed
- Google Scholar
(2011) Tracking brain states under general anesthesia by using global coherence analysis
PNAS 108:8832–8837.

https://doi.org/10.1073/pnas.1017041108
- PubMed
- Google Scholar
1. Cohen MX
(2008) Assessing transient cross-frequency coupling in EEG data
Journal of Neuroscience Methods 168:494–499.

https://doi.org/10.1016/j.jneumeth.2007.10.012
- PubMed
- Google Scholar
1. Cohen MX
(2017) Multivariate cross-frequency coupling via generalized eigendecomposition
eLife 6:e21792.

https://doi.org/10.7554/eLife.21792
- PubMed
- Google Scholar
1. Cohen MX
(2018) Using spatiotemporal source separation to identify prominent features in multichannel data without sinusoidal filters
The European Journal of Neuroscience 48:2454–2465.

https://doi.org/10.1111/ejn.13727
- PubMed
- Google Scholar
Conference
1. Das P
2. Purdon PL
(2022) Extracting common oscillatory time-courses from multichannel recordings: oscillation component analysis
56th Asilomar Conference on Signals, Systems, and Computers. pp. 602–606.

https://doi.org/10.1109/IEEECONF56349.2022.10052084
- Google Scholar
Software
1. Das P
(2024) Purdonlabmeeg, version swh:1:rev:c7e69979d99174348bac7cdd42815ff5815bee52
Software Heritage.

https://archive.softwareheritage.org/swh:1:dir:308d24c1c964363d635456bf74cb6828df80f896;origin=https://github.com/proloyd/purdonlabmeeg;visit=swh:1:snp:2f7581043381ea7d2527b54f4935dcd6509d633e;anchor=swh:1:rev:c7e69979d99174348bac7cdd42815ff5815bee52
1. de Cheveigné A
2. Parra LC
(2014) Joint decorrelation, a versatile tool for multichannel data analysis
NeuroImage 98:487–505.

https://doi.org/10.1016/j.neuroimage.2014.05.068
- PubMed
- Google Scholar
1. de Cheveigné A
2. Arzounian D
(2015) Scanning for oscillations
Journal of Neural Engineering 12:066020.

https://doi.org/10.1088/1741-2560/12/6/066020
- PubMed
- Google Scholar
(1977) Maximum likelihood from incomplete data via the EM Algorithm
Journal of the Royal Statistical Society Series B 39:1–22.

https://doi.org/10.1111/j.2517-6161.1977.tb01600.x
- Google Scholar
1. Desikan RS
2. Ségonne F
3. Fischl B
4. Quinn BT
5. Dickerson BC
6. Blacker D
7. Buckner RL
8. Dale AM
9. Maguire RP
10. Hyman BT
11. Albert MS
12. Killiany RJ
(2006) An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest
NeuroImage 31:968–980.

https://doi.org/10.1016/j.neuroimage.2006.01.021
- PubMed
- Google Scholar
1. Donoghue T
2. Haller M
3. Peterson EJ
4. Varma P
5. Sebastian P
6. Gao R
7. Noto T
8. Lara AH
9. Wallis JD
10. Knight RT
11. Shestyuk A
12. Voytek B
(2020) Parameterizing neural power spectra into periodic and aperiodic components
Nature Neuroscience 23:1655–1665.

https://doi.org/10.1038/s41593-020-00744-x
- PubMed
- Google Scholar
1. Efron B
2. Morris C
(1973) Stein’s estimation rule and its competitors—an empirical bayes approach
Journal of the American Statistical Association 68:117–130.

https://doi.org/10.1080/01621459.1973.10481350
- Google Scholar
1. Efron B
2. Morris C
(1975) Data analysis using stein’s estimator and its generalizations
Journal of the American Statistical Association 70:311–319.

https://doi.org/10.1080/01621459.1975.10479864
- Google Scholar
1. Fahrmeir L
(1992) Posterior mode estimation by extended kalman filtering for multivariate dynamic generalized linear models
Journal of the American Statistical Association 87:501–509.

https://doi.org/10.1080/01621459.1992.10475232
- Google Scholar
1. Galinsky VL
2. Frank LR
(2020) Universal theory of brain waves: from linear loops to nonlinear synchronized spiking and collective brain rhythms
Physical Review Research 2:023061.

https://doi.org/10.1103/PhysRevResearch.2.023061
- PubMed
- Google Scholar
(2014) MNE software for processing MEG and EEG data
NeuroImage 86:446–460.

https://doi.org/10.1016/j.neuroimage.2013.10.027
- PubMed
- Google Scholar
(2023) Characterizing endogenous delta oscillations in human MEG
Scientific Reports 13:11031.

https://doi.org/10.1038/s41598-023-37514-1
- PubMed
- Google Scholar
Book
1. Harvey AC
(1990) Forecasting, structural time series models and the kalman filter
Cambridge University Press.

https://doi.org/10.1017/CBO9781107049994
- Google Scholar
1. Harvey AC
2. Trimbur TM
(2003) General model-based filters for extracting cycles and trends in economic time series
Review of Economics and Statistics 85:244–255.

https://doi.org/10.1162/003465303765299774
- Google Scholar
1. Haufe S
2. Meinecke F
3. Görgen K
4. Dähne S
5. Haynes J-D
6. Blankertz B
7. Bießmann F
(2014) On the interpretation of weight vectors of linear models in multivariate neuroimaging
NeuroImage 87:96–110.

https://doi.org/10.1016/j.neuroimage.2013.10.067
- PubMed
- Google Scholar
1. Helfrich RF
2. Knight RT
(2016) Oscillatory dynamics of prefrontal cognitive control
Trends in Cognitive Sciences 20:916–930.

https://doi.org/10.1016/j.tics.2016.09.007
- Google Scholar
Conference
1. Hinton GE
2. van Camp D
(1993) Keeping the Neural Networks Simple by Minimizing the Description Length of the Weights
Proceedings of the Sixth Annual Conference on Computational Learning Theory - COLT ’93. pp. 5–13.

https://doi.org/10.1145/168304.168306
- Google Scholar
1. Hyvärinen A
2. Oja E
(2000) Independent component analysis: algorithms and applications
Neural Networks 13:411–430.

https://doi.org/10.1016/s0893-6080(00)00026-5
- PubMed
- Google Scholar
1. Hyvärinen A
(2013) Independent component analysis: recent advances
Philosophical Transactions of the Royal Society A 371:20110534.

https://doi.org/10.1098/rsta.2011.0534
- Google Scholar
Book
1. Jain J
2. Li H
3. Cauley S
4. Koh CK
5. Balakrishnan V
(2007)
Numerically stable algorithms for inversion of block tridiagonal and banded matrices

Department of Electrical and Computer Engineering Technical Reports.
- Google Scholar
1. Jung TP
2. Makeig S
3. Humphries C
4. Lee TW
5. McKeown MJ
6. Iragui V
7. Sejnowski TJ
(2000) Removing electroencephalographic artifacts by blind source separation
Psychophysiology 37:163–178.

https://doi.org/10.1111/1469-8986.3720163
- PubMed
- Google Scholar
(1999) Independent component analysis using an extended infomax algorithm for mixed subgaussian and supergaussian sources
Neural Computation 11:417–441.

https://doi.org/10.1162/089976699300016719
- PubMed
- Google Scholar
1. Lewis LD
2. Weiner VS
3. Mukamel EA
4. Donoghue JA
5. Eskandar EN
6. Madsen JR
7. Anderson WS
8. Hochberg LR
9. Cash SS
10. Brown EN
11. Purdon PL
(2012) Rapid fragmentation of neuronal networks at the onset of propofol-induced unconsciousness
PNAS 109:E3377–E3386.

https://doi.org/10.1073/pnas.1210907109
- PubMed
- Google Scholar
1. Lopes da Silva F
(2013) EEG and MEG: relevance to neuroscience
Neuron 80:1112–1128.

https://doi.org/10.1016/j.neuron.2013.10.017
- PubMed
- Google Scholar
(2019) Measuring transient phase-amplitude coupling using local mutual information
NeuroImage 185:361–378.

https://doi.org/10.1016/j.neuroimage.2018.10.034
- PubMed
- Google Scholar
1. Matsuda T
2. Komaki F
(2017a) Time series decomposition into oscillation components and phase estimation
Neural Computation 29:332–367.

https://doi.org/10.1162/NECO_a_00916
- PubMed
- Google Scholar
1. Matsuda T
2. Komaki F
(2017b) Multivariate time series decomposition into oscillation components
Neural Computation 29:2055–2075.

https://doi.org/10.1162/NECO_a_00981
- PubMed
- Google Scholar
Book
1. Mitra P
2. Bokil H
(2007) Observed brain dynamics
Oxford University Press.

https://doi.org/10.1093/acprof:oso/9780195178081.001.0001
- Google Scholar
1. Morris CN
(1983) Parametric empirical bayes inference: theory and applications
Journal of the American Statistical Association 78:47–55.

https://doi.org/10.1080/01621459.1983.10477920
- Google Scholar
Book
1. Muirhead RJ
(1982) Aspects of multivariate statistical theory
Wiley.

https://doi.org/10.1002/9780470316559
- Google Scholar
1. Murphy M
2. Bruno M-A
3. Riedner BA
4. Boveroux P
5. Noirhomme Q
6. Landsness EC
7. Brichant J-F
8. Phillips C
9. Massimini M
10. Laureys S
11. Tononi G
12. Boly M
(2011) Propofol anesthesia and sleep: a high-density EEG study
Sleep 34:283–91A.

https://doi.org/10.1093/sleep/34.3.283
- PubMed
- Google Scholar
Book
1. Neal RM
2. Hinton GE
(1998) Learning in graphical models
In: Jorden MI, editors. A view of the em algorithm that justifies incremental, sparse, and other variants. Springer Netherlands. pp. 355–368.

https://doi.org/10.1007/978-94-011-5014-9_12
- Google Scholar
1. Neumaier A
2. Schneider T
(2001) Estimation of parameters and eigenmodes of multivariate autoregressive models
ACM Transactions on Mathematical Software 27:27–57.

https://doi.org/10.1145/382043.382304
- Google Scholar
(2011) A novel method for reliable and fast extraction of neuronal EEG/MEG oscillations on the basis of spatio-spectral decomposition
NeuroImage 55:1528–1535.

https://doi.org/10.1016/j.neuroimage.2011.01.057
- PubMed
- Google Scholar
Book
1. Nise NS
(2011)
Control systems engineering

Wiley.
- Google Scholar
Book
1. Parra LC
(1998) Temporal Models in Blind Source Separation
Springer.

https://doi.org/10.1007/BFb0054001
- Google Scholar
1. Parra L
2. Sajda P
(2003) Blind source separation via generalized eigenvalue decomposition
Journal of Machine Learning Research 4:1261–1269.

https://doi.org/10.5555/945365.964305
- Google Scholar
Conference
1. Pearlmutter B
2. Parra L
(1996)
Maximum likelihood blind source separation: a context-sensitive generalization of ICA

Advances in neural information processing systems. pp. 613–619.
- Google Scholar
Conference
1. Perley A
2. Coleman TP
(2022) A mutual information measure of phase-amplitude coupling using high dimensional sparse models
44th Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC). pp. 21–24.

https://doi.org/10.1109/EMBC48229.2022.9871816
- Google Scholar
1. Purdon PL
2. Pierce ET
3. Mukamel EA
4. Prerau MJ
5. Walsh JL
6. Wong KFK
7. Salazar-Gomez AF
8. Harrell PG
9. Sampson AL
10. Cimenser A
11. Ching S
12. Kopell NJ
13. Tavares-Stoeckel C
14. Habeeb K
15. Merhar R
16. Brown EN
(2013) Electroencephalogram signatures of loss and recovery of consciousness from propofol
PNAS 110:E1142–E1151.

https://doi.org/10.1073/pnas.1221180110
- PubMed
- Google Scholar
Book
1. Quinn A
2. Šmídl V
(2006)
The Variational Bayes Method in Signal Processing

Springer NetLibrary, Inc.
- Google Scholar
(2021) Delineating between-subject heterogeneity in alpha networks with Spatio-Spectral Eigenmodes
NeuroImage 240:118330.

https://doi.org/10.1016/j.neuroimage.2021.118330
- PubMed
- Google Scholar
1. Rehman N
2. Mandic DP
(2010) Multivariate empirical mode decomposition
Proceedings of the Royal Society A 466:1291–1302.

https://doi.org/10.1098/rspa.2009.0502
- Google Scholar
1. Robbins H
(1964) The empirical bayes approach to statistical decision problems
The Annals of Mathematical Statistics 35:1–20.

https://doi.org/10.1214/aoms/1177703729
- Google Scholar
Conference
(2017) Canonical correlation analysis of EEG for classification of motor imagery
2017 IEEE International Conference on Systems, Man and Cybernetics (SMC). pp. 2317–2321.

https://doi.org/10.1109/SMC.2017.8122967
- Google Scholar
(1997) Phase synchronization in driven and coupled chaotic oscillators
IEEE Transactions on Circuits and Systems I 44:874–881.

https://doi.org/10.1109/81.633876
- Google Scholar
1. Schaworonkow N
2. Nikulin VV
(2022) Is sensor space analysis good enough? Spatial patterns as a tool for assessing spatial mixing of EEG/MEG rhythms
NeuroImage 253:119093.

https://doi.org/10.1016/j.neuroimage.2022.119093
- PubMed
- Google Scholar
1. Shumway RH
2. Stoffer DS
(1982) An approach to time series smoothing and forecasting using the EM algorithm
Journal of Time Series Analysis 3:253–264.

https://doi.org/10.1111/j.1467-9892.1982.tb00349.x
- Google Scholar
1. Soulat H
2. Stephen EP
3. Beck AM
4. Purdon PL
(2022) State space methods for phase amplitude coupling analysis
Scientific Reports 12:15940.

https://doi.org/10.1038/s41598-022-18475-3
- PubMed
- Google Scholar
1. Stevenson IH
2. Kording KP
(2011) How advances in neural recording affect data analysis
Nature Neuroscience 14:139–142.

https://doi.org/10.1038/nn.2731
- PubMed
- Google Scholar
1. Uusitalo MA
2. Ilmoniemi RJ
(1997) Signal-space projection method for separating MEG or EEG into components
Medical & Biological Engineering & Computing 35:135–140.

https://doi.org/10.1007/BF02534144
- Google Scholar
1. Wang XJ
(2010) Neurophysiological and computational principles of cortical rhythms in cognition
Physiological Reviews 90:1195–1268.

https://doi.org/10.1152/physrev.00035.2008
- PubMed
- Google Scholar
1. Weiner VS
2. Zhou DW
3. Kahali P
4. Stephen EP
5. Peterfreund RA
6. Aglio LS
7. Szabo MD
8. Eskandar EN
9. Salazar-Gomez AF
10. Sampson AL
11. Cash SS
12. Brown EN
13. Purdon PL
(2023) Propofol disrupts alpha dynamics in functionally distinct thalamocortical networks during loss of consciousness
PNAS 120:e2207831120.

https://doi.org/10.1073/pnas.2207831120
- PubMed
- Google Scholar
Book
1. Wiener N
(1966)
Nonlinear problems in random theory

MIT press.
- Google Scholar
1. Wodeyar A
2. Schatza M
3. Widge AS
4. Eden UT
5. Kramer MA
(2021) A state space modeling approach to real-time phase estimation
eLife 10:e68803.

https://doi.org/10.7554/eLife.68803
- PubMed
- Google Scholar
(2007) Theta phase resetting and the error-related negativity
Psychophysiology 44:39–49.

https://doi.org/10.1111/j.1469-8986.2006.00482.x
- PubMed
- Google Scholar
Conference
1. Zhao Q
2. Xu M
3. Fränti P
(2008) Knee point detection on bayesian information criterion
20th IEEE International Conference on Tools with Artificial Intelligence (ICTAI). pp. 431–438.

https://doi.org/10.1109/ICTAI.2008.154
- Google Scholar

Article and author information

Author details

Proloy Das

Department of Anesthesiology, Perioperative and Pain Medicine, Stanford University, Stanford, United States

Contribution
Conceptualization, Resources, Data curation, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing - original draft, Writing – review and editing

For correspondence
proloy@stanford.edu

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-8807-042X
Mingjian He
1. Department of Anesthesiology, Perioperative and Pain Medicine, Stanford University, Stanford, United States
2. Department of Psychology, Stanford University, Stanford, United States
Contribution
Resources, Data curation, Visualization

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-6688-8693
Patrick L Purdon
1. Department of Anesthesiology, Perioperative and Pain Medicine, Stanford University, Stanford, United States
2. Department of Bioengineering, Stanford University, Stanford, United States
Contribution
Conceptualization, Resources, Supervision, Funding acquisition, Investigation, Methodology, Project administration, Writing – review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0003-0080-3340

Funding

National Institutes of Health (R01AG054081-01A1)

Patrick L Purdon

Tiny Blue Dot Foundation

Patrick L Purdon

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

This work was supported by the National Institutes of Health (grant no. R01AG054081-01A1) and Tiny Blue Dot Foundation. As requested by the Human Connectome Project, the following text is copied verbatim: MEG data in sSection ‘Resting -state MEG recording’ ‘were provided by the Human Connectome Project, WU-Minn Consortium (Pprincipal Iinvestigators: David Van Essen and Kamil Ugurbil; 1U54MH091657) funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at Washington University.

Version history

Preprint posted: February 8, 2024
Sent for peer review: February 29, 2024
Reviewed Preprint version 1: April 18, 2024
Reviewed Preprint version 2: July 12, 2024
Version of Record published: August 15, 2024

Cite all versions

You can cite all versions using the DOI https://doi.org/10.7554/eLife.97107. This DOI represents all versions, and will always resolve to the latest one.

Copyright

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.