
==== Front
Hum Brain Mapp
Hum Brain Mapp
10.1002/(ISSN)1097-0193
HBM
Human Brain Mapping
1065-9471
1097-0193
John Wiley & Sons, Inc. Hoboken, USA

10.1002/hbm.70018
HBM70018
Research Article
Research Article
Comparison between EEG and MEG of static and dynamic resting‐state networks
Cho et al.
Cho SungJun https://orcid.org/0000-0002-2787-3739
1 sungjun.cho@ndcn.ox.ac.uk

van Es Mats 1
Woolrich Mark 1
Gohil Chetan 1
1 Oxford Centre for Human Brain Activity, Wellcome Centre for Integrative Neuroimaging, Department of Psychiatry University of Oxford Oxford UK
* Correspondence
SungJun Cho, Oxford Centre for Human Brain Activity, Wellcome Centre for Integrative Neuroimaging, Department of Psychiatry, University of Oxford, Oxford OX3 7JK, UK.
Email: sungjun.cho@ndcn.ox.ac.uk

04 9 2024
9 2024
45 13 10.1002/hbm.v45.13 e7001829 7 2024
08 5 2024
20 8 2024
© 2024 The Author(s). Human Brain Mapping published by Wiley Periodicals LLC.
https://creativecommons.org/licenses/by-nc-nd/4.0/ This is an open access article under the terms of the http://creativecommons.org/licenses/by-nc-nd/4.0/ License, which permits use and distribution in any medium, provided the original work is properly cited, the use is non‐commercial and no modifications or adaptations are made.

Abstract

The characterisation of resting‐state networks (RSNs) using neuroimaging techniques has significantly contributed to our understanding of the organisation of brain activity. Prior work has demonstrated the electrophysiological basis of RSNs and their dynamic nature, revealing transient activations of brain networks with millisecond timescales. While previous research has confirmed the comparability of RSNs identified by electroencephalography (EEG) to those identified by magnetoencephalography (MEG) and functional magnetic resonance imaging (fMRI), most studies have utilised static analysis techniques, ignoring the dynamic nature of brain activity. Often, these studies use high‐density EEG systems, which limit their applicability in clinical settings. Addressing these gaps, our research studies RSNs using medium‐density EEG systems (61 sensors), comparing both static and dynamic brain network features to those obtained from a high‐density MEG system (306 sensors). We assess the qualitative and quantitative comparability of EEG‐derived RSNs to those from MEG, including their ability to capture age‐related effects, and explore the reproducibility of dynamic RSNs within and across the modalities. Our findings suggest that both MEG and EEG offer comparable static and dynamic network descriptions, albeit with MEG offering some increased sensitivity and reproducibility. Such RSNs and their comparability across the two modalities remained consistent qualitatively but not quantitatively when the data were reconstructed without subject‐specific structural MRI images.

Our study shows that medium‐density electroencephalography systems can offer static and dynamic network descriptions of resting‐state networks that are comparable in quality to those identified by high‐density magnetoencephalography systems. This comparability was reproducible within and across two modalities and remained qualitatively consistent when the data were source‐reconstructed without structural magnetic resonance imaging.

electroencephalography
hidden Markov Model
magnetoencephalography
neural dynamics
resting‐state networks
Wellcome Trust 10.13039/100010269 106183/Z/14/Z 203139/Z/16/Z 215573/Z/19/Z NIHR Oxford Biomedical Research Centre 10.13039/501100013373 NIHR203316 New Therapeutics in Alzheimer's Diseases (NTAD), UK Medical Research Council 10.13039/501100000265 Dementias Platform UK 10.13039/501100021686 RG89702 RG94383 source-schema-version-number2.0
cover-dateSeptember 2024
details-of-publishers-convertorConverter:WILEY_ML3GV2_TO_JATSPMC version:6.4.8 mode:remove_FC converted:04.09.2024
Cho, S. , van Es, M. , Woolrich, M. , & Gohil, C. (2024). Comparison between EEG and MEG of static and dynamic resting‐state networks. Human Brain Mapping, 45 (13 ), e70018. 10.1002/hbm.70018
==== Body
pmc1 INTRODUCTION

Resting‐state networks (RSNs) were first observed in functional magnetic resonance imaging (fMRI) and positron emission tomography in the late 20th century and later independently identified in magnetoencephalography (MEG) and electroencephalography (EEG) (Biswal et al., 1995; Shulman et al., 1997). A study by Brookes et al. (2011) provided robust evidence for the electrophysiological basis of RSNs, employing temporal independent component analysis (ICA) on MEG amplitude envelopes to demonstrate that the spatial patterns of brain networks in MEG resemble those of the RSNs previously identified by fMRI (Brookes et al., 2011). While early descriptions of RSNs were essentially static (i.e., time‐averaged), further research has revealed that they can be well described as activating in a transient manner, with different networks activating at distinct time points (Baker et al., 2014; Gohil et al., 2022).

The contribution of electrophysiology to our understanding of RSNs is underscored by their functional significance. MEG and EEG studies have allowed RSNs to be linked to sensorimotor functions at fast subsecond timescales, as well as to a range of cognitive tasks concerning working memory and language comprehension (Barnes et al., 2015; Becker et al., 2020; Quinn et al., 2018). Research on neuropsychiatric disorders has also elucidated the functional relevance of RSNs across various pathological conditions, including Alzheimer's disease, major depressive disorder, and schizophrenia, highlighting their utility in understanding brain dysfunction (Leuchter et al., 2012; Lottman et al., 2019; Sitnikova et al., 2018). From a technical standpoint, electrophysiological analysis of RSNs offers unparalleled insight into the spectral and temporal domains. M/EEG‐derived functional networks allow for the investigation of their spectral contents; empirical evidence suggests that neural activities can be segmented into distinct states of global synchronous fluctuations, each characterised by specific oscillatory frequencies (de Pasquale et al., 2010; Mantini et al., 2007). Furthermore, the superior temporal resolution of electrophysiological approaches facilitates a detailed examination of the transient dynamics of resting‐state and task‐based networks. For example, previous research indicates that M/EEG‐derived RSNs activate on faster timescales (~100–200 ms) than those discerned through fMRI (~10 s) (Baker et al., 2014; Gohil et al., 2022).

To date, studies examining electrophysiological RSNs, particularly those based on the amplitude time‐courses of brain regions, have affirmed a reasonable comparability of EEG RSNs to their MEG or fMRI counterparts (Knyazev et al., 2016; Liu et al., 2017; Siems et al., 2016; Sockeel et al., 2016). However, these studies relied on static analyses (by which we mean that they do not provide a time‐varying description of interactions between brain areas), such as spatial ICA or seed‐based correlation. More recently, Coquelet et al. (2020) evaluated the correspondence of RSNs between 306‐channel MEG and 256‐channel EEG data (Coquelet et al., 2020). This study revealed that while static functional connectivity (FC) patterns are generally consistent across modalities in the resting state, dynamic FC shows less agreement. Their dynamic FC methods, however, used sliding windows and clustering approaches, which are heavily dependent on user‐defined hyperparameters and compromise temporal resolution. Finally, studies to date have typically employed high‐density EEG systems (>100 channels) (Coquelet et al., 2020; Knyazev et al., 2016; Liu et al., 2017), facilitating the computation of RSNs in brain space, where spatial leakage can be more appropriately controlled and the networks are more interpretable (Colclough et al., 2015). Yet, this overlooks the wider availability and clinical utility of lower‐density EEG systems, for which the ability to reliably estimate RSNs in brain space would be a valuable asset.

Acknowledging these limitations, this article aims to establish the efficacy of EEG in providing static and dynamic functional network perspectives akin to those obtained from MEG and fMRI within brain space, with a state‐of‐the‐art dynamic modelling technique—time‐delay embedded hidden Markov model (TDE‐HMM) (Vidaurre, Hunt, et al., 2018) and without the need for high density data. Drawing on existing literature, we hypothesised that both EEG and MEG will provide comparable static and dynamic RSN descriptions, with dynamic features anticipated to exhibit less consistency between modalities. Using MEG RSNs as a benchmark for assessing fast oscillatory networks, we compared static and dynamic network descriptions against those derived from EEG, with the aim of establishing a basis set of EEG RSNs that can be identified with medium‐density EEG data (~61 channels). In addition, we also examined whether M/EEG RSNs, derived from data reconstructed without subject‐specific structural MRI (sMRI) images, can yield network descriptions similar to those obtained using subject sMRI images without loss of quality.

This endeavour seeks to demonstrate that MEG and EEG provide similar and reliable network descriptions, allowing for research that can leverage the accessibility and cost‐effectiveness of EEG in the exploration of brain oscillatory networks in both health and disease. Additionally, this work provides a publicly available set of M/EEG RSNs, complete with scripts for generating precise network inferences, thereby contributing to our comprehensive understanding of the complementary potential of MEG and EEG in brain network research.

2 MATERIALS AND METHODS

2.1 Datasets

In this study, we used the openly available resting‐state EEG and MEG datasets called the Leipzig Study for Mind‐Body‐Emotion Interactions (LEMON) dataset and the Cambridge Centre for Ageing and Neuroscience (Cam‐CAN) dataset, respectively. Both datasets were obtained from healthy adults, and recordings were measured while participants were seated. The two datasets were matched to have an identical age distribution and number of subjects. To specify, each dataset was divided into subsets of a 5‐year age interval, and for imbalances between each subset of two datasets, subjects were randomly subsampled from the larger subset to ensure parity. The age‐matched datasets comprised a total of 96 subjects, consisting of 60 young (aged between 20 and 35 years) and 36 old (aged between 55 and 80 years) participants1 (Figure A1). In the EEG data, the young group comprised 44 females and 16 males, while the old group consisted of 23 females and 13 males. In the MEG data, the young group comprised 25 females and 35 males, while the old group consisted of 19 females and 17 males. The protocols and demographics related to the data collection are outlined in Babayan et al. (2019) for LEMON and Shafto et al. (2014) for Cam‐CAN in detail and are partly repeated here for clarity.

2.1.1 EEG data

The EEG recordings were obtained using a BrainAMP MR plus amplifier with 62‐channel ActiCAP electrodes (Brain Products, Gilching, Germany). These channels consisted of 61 EEG electrodes and 1 vertical electrooculogram channel. A surrogate horizontal electrooculogram channel was additionally added by taking the difference between channels F7 and F8. The channel montage adapted the 10–10 layout, which was referenced and grounded at FCz and the sternum, respectively. The acquisition time of the recordings was 16, with 1 min blocks alternating between eyes‐closed and eyes‐open conditions. Signals were collected with a sampling rate of 2500 Hz and band‐pass filtered between 0.015 Hz and 1 kHz. For coregistration purposes, the T1‐weighted sMRI was collected over 8 min 22 s, using Magnetization‐Prepared 2 RApid Acquisition Gradient Echos (MP2RAGE) sequences with a 3 T MAGNETOM Verio scanner (Siemens Healthcare, Erlangen, Germany) equipped with a 32‐channel head coil.

2.1.2 MEG data

The MEG recordings were obtained using a 306‐channel Vector View system (Elekta Neuromag, Helsinki, Finland), consisting of 102 magnetometers and 204 orthogonal planar gradiometers. The acquisition time of the recordings was 8 min 40 s, with the first 20 s discarded. Only eyes‐closed resting‐state data were recorded at a sampling rate of 1 kHz, band‐pass filtered between 0.03 and 330 Hz. The eye‐ and pulse‐related artefacts were monitored using two pairs of biopolar EOG electrodes and one pair of bipolar electrocardiogram (ECG) electrodes, respectively. For coregistration, the T1‐weighted sMRI was recorded for 4 min 32 s, using a Magnetization‐Prepared RApid Gradient Echos sequence with a 3 T TIM Trio scanner (Siemens Healthcare, Erlangen, Germany) equipped with a 32‐channel head coil.

2.2 Data preprocessing and source reconstruction

The preprocessing and source reconstruction processes largely followed the pipeline outlined in Gohil et al. (2023). Detailed steps of the pipeline applied to the M/EEG datasets were slightly different due to the unique data characteristics intrinsic to each dataset. However, they were aligned as closely as possible to mitigate any unknown variance that could arise from discrepancies in these steps. Before applying the pipeline, the MEG data were Maxfiltered using the temporal signal‐space separation method to separate recordings of neural activity within the brain from any external noise sources (Quinn et al., 2018). The entire process of preprocessing and source reconstruction was conducted using the OHBA Software Library (OSL) package (Quinn et al., 2023).

2.2.1 Preprocessing

The LEMON and Cam‐CAN data were first cropped to exclude the first 30 s of the recordings to account for the subject acclimation and stabilisation of measurement devices. These data were band‐pass filtered between 0.5 and 125 Hz using a fifth‐order IIR Butterworth filter and subsequently notch filtered at 50 and 100 Hz to remove mains artefact (and additionally at 88 Hz for Cam‐CAN to remove spike artefacts). The data were then downsampled to 250 Hz. The bad segments and channels with significantly high variance were removed using an automated generalised‐extreme studentised deviate algorithm (Rosner, 1983). For Cam‐CAN, bad segment and channel detection procedures were applied separately to different sensor types (i.e., magnetometers and gradiometers). The data were further denoised by applying an FastICA decomposition (Hyvarinen, 1999) to the M/EEG data, decomposing signals into 64 and 54 components, respectively. Components associated with ocular artefacts (i.e., blinks and saccades) were removed, and ECG‐related artefacts were additionally excluded for Cam‐CAN. Prior to the ICA decomposition, bad segments from EOG recordings were identified to prevent the inadvertent removal of pertinent M/EEG signals due to concurrent noise artefacts in both M/EEG and EOG. Finally, to keep the data dimension consistent across the subjects, any bad channels detected earlier were interpolated from ICA‐cleaned data using spherical spline interpolation (Perrin et al., 1989).

2.2.2 Source reconstruction

The M/EEG data were coregistered to sMRI image and digitised head points using the affine transformation algorithm in OSL. For LEMON, a boundary element model (BEM) with triple layers of scalp, skull, and cortex surfaces was employed as a head model; for Cam‐CAN, a BEM with a single layer of a scalp surface was employed. The surfaces of a scalp, inner skull, and brain were extracted from sMRI data using FSL's BET tool (Jenkinson et al., 2005; Smith, 2002). The nose was not included during the coregistration step, since the sMRI images were defaced.

After coregistration, the preprocessed sensor data were band‐pass filtered between 1 and 45 Hz and reconstructed onto an 8 mm isotropic grid using a linearly constrained minimum variance (LCMV) volumetric beamformer (Veen & Buckley, 1988). The rank of a data covariance matrix used to compute the beamformer weights was set to 50,2 which regularises the covariance estimation and is comfortably below the rank of both the MEG (following Maxfiltering) and EEG data. Voxels were then parcellated into 38 anatomically defined regions. Source reconstructed signals were obtained by applying principal component analysis (PCA) to these voxel‐wise data, wherein a time series of each parcel would be the first principal component (PC) that explains most variance from all voxels within a parcel. The symmetric multivariate leakage reduction technique (Colclough et al., 2015) was applied to minimise spurious correlations between all parcels and reduce source/spatial leakage including so‐called “inherited” or “ghost interactions” (Palva et al., 2018). Finally, the combination of the LCMV volumetric beamformer and the PCA used in the parcel time series calculations means that the sign of the estimated parcel time series varies arbitrarily across subjects. The sign of each channel was therefore matched across subjects using a random flip algorithm described in Vidaurre, Hunt, et al. (2018).

The source reconstruction procedure was repeated using the standard Montreal Neurological Institute (MNI) T1 structural brain image instead of subject‐specific sMRI images. We employed this approach to explore whether similar M/EEG RSN features could be discerned using a standard brain structural image (i.e., without subject sMRI files) and if the comparability between MEG and EEG remains under these conditions. Unless indicated otherwise, all results presented herein are based on the data reconstructed with the subject‐specific sMRI images.

2.2.3 Data organisation

To make the LEMON and Cam‐CAN dataset more comparable, only eyes‐closed segments were extracted from the source reconstructed EEG LEMON data. After extracting the eyes‐closed segments, the mean and standard deviation of the data lengths were 476.4 ± 23.78 s for LEMON and 521.9 ± 59.77 s for Cam‐CAN, with an average data length ratio between the two resulting in 0.9128. Given the data length roughly matched and to retain as much MEG data as possible, we did not further shorten the Cam‐CAN dataset to match LEMON. All data analyses delineated below took place in source space.

2.3 Hidden Markov model

Once the data were source reconstructed, the next step is to identify dynamic RSNs and pinpoint the time periods during which they activate. The TDE‐HMM is a generative model that can partition time series data into a sequence of recurring functional brain networks called states (Vidaurre, Hunt, et al., 2018). Each state represents an RSN characterised by distinct spatiotemporal patterns of spectral content. The basic structure of the TDE‐HMM is depicted in Figure 1.

FIGURE 1 Overview of brain network feature analysis. This figure illustrates the methodology for extracting brain network features from MEG and EEG time series data. The brain network features we used consisted of power spectral density (PSD), spatial power distribution, and functional connectivity maps derived from spectral information within M/EEG recordings (top). The time series are segmented into discrete state activations via a Hidden Markov Model (HMM), resulting in an HMM state time course (middle). This process allows for the differentiation between dynamic network features—obtained from state‐specific descriptions—and static network features, which are computed by time‐averaging across all states. The HMM parameterises functional resting‐state networks using a transition probability matrix and a multivariate Gaussian observation model with distinct means and covariances for each state by minimising variational free energy (bottom). The transition matrix quantifies the likelihood of transition from one state to another over consecutive time points, whereas the observation model elucidates the data generating process based on the latent (i.e., hidden state) representations. When combined with time‐delay embedding, the HMM can identify states with distinct parcel power spectra and coherence networks (not shown).

The HMM typically comprises two primary components: a set of hidden states and an observation model for each state. Our observation model is a multivariate normal (Gaussian) distribution that has a distinct covariance matrix for each state. Crucially, the HMM is fit to TDE data, which augments the parcel time series with extra channels containing time‐lagged versions of the original time series. Using TDE makes the state‐specific covariance matrices sensitive to the power and frequency of oscillations in the original time series and results in states that have distinct power spectra for each parcel and coherence networks (Gohil et al., 2023).

The HMM parameters were inferred using the variational Bayesian (VB) inference (Bishop, 2006; Rezek & Roberts, 2005). This approach recasts the inference process as an optimisation problem, wherein the objective is to minimise the variational free energy. Using the VB method allows the model to learn uncertainty in parameter estimates as probabilistic distributions. To scale the learning process to large datasets, we utilised the stochastic version of the VB method adapted from (Vidaurre, Abeysuriya, et al., 2018). As a result, the HMM outputs a time course of posterior probabilities, which signifies the probability of a state activation at each time point. By extracting the maximum a posteriori probabilities across all states by applying an argmax operation to these probabilities, we can generate a Viterbi path, or a state time course, delineating the most probable, yet mutually exclusive, state at each time point. The training process is described below. For more details on the design and implementation of the TDE‐HMM, readers should refer to Vidaurre, Hunt, et al. (2018) and Gohil et al. (2023).

2.3.1 Data preparation

Before training a TDE‐HMM model, the source reconstructed data were prepared to allow the model to learn states with unique spatio‐spectral patterns. First, the data were TDE with ±7 lags. With a 38‐region parcellation, this embedding resulted in a total of 570 channels. PCA was then applied to reduce the dimensionality down to 80 channels to reduce memory requirements, and the transformed data were subsequently standardised over time. Since we are primarily interested in the transient patterns of spectral events, which are reflected in the state covariance matrices, we fixed the state mean vectors of the observation model to zero during the training process (Gohil et al., 2023). These prepared data served as the input training data, which was then shuffled and batched prior to model training.

2.3.2 Hyperparameters

Training the model involves several key hyperparameters. First, the number of states was set to 6 to facilitate a reasonably low‐dimensional data representation, which tend to be more reproducible and interpretable (see Section 2.3.5 for its validation). Next, we selected 64 as the batch size necessary for the stochastic VB inference. The input sequence length for batching the time series data was 800 samples. This length ensures minimal discontinuities between batched data, leading to less noisy updates of the parameter estimates. Finally, the Adam optimiser (Kingma & Ba, 2014) was used for updating the trainable parameters with a learning rate of 0.01 over a total of 20 training epochs.

2.3.3 Model initialisation

As the optimisation process for the HMM is stochastic, a certain level of variability in the final model parameters, referred to as run‐to‐run variability, is inevitable. To mitigate this, we use a number of random starts, or “initialisations.” This initialisation process was repeated five times, trained for two epochs in each repetition. The initialisation that yielded the lowest free energy was selected for the full model training. This approach has previously been shown to produce reproducible results (Gohil et al., 2023).

2.3.4 Model runs

Despite the initialisation technique above, different model runs may converge to different local optima due to the stochastic nature of the parameter optimisation process. To obtain the best possible parameter estimates, therefore, the HMM model was trained 10 times using the identical training dataset for both LEMON and Cam‐CAN. From these 10 runs, a model with the lowest free energy was selected as the best run. This procedure was again repeated five times,3 yielding a total of 50 runs (and 5 best runs) for each dataset. The best of the five best runs, hereinafter denoted as the optimal run, was chosen to be used for all subsequent analyses. As the order of inferred states is random for each model, the EEG and MEG states from their respective optimal HMM runs were manually aligned by eye4 based on their power and FC spatial patterns.

2.3.5 Model validation

Before the datasets were analysed, we validated that 6 states produce interpretable RSNs without any loss in training quality. For validation purposes, we trained 50 HMM runs with 8 and 10 states, following the method outlined in Section 2.3.4. For each set of dynamic RSNs (for 6, 8, and 10 states), we counted how many times a particular RSN appears across the best runs by qualitatively comparing its PSDs, power maps, and FC networks (see Section 2.4). Results indicated that with six states, the HMM consistently converged to the same six RSNs across the best runs. In contrast, for 8 and 10 states, approximately 6 RSNs demonstrated consistency across the best runs, while the remaining states exhibited greater variability across runs for both EEG and MEG (i.e., the additional states corresponded to less reproducible RSNs). Figure A7 shows 10‐state dynamic M/EEG RSNs, inferred from the optimal model run, as an example. The configuration of six states that guarantees the network features with the lowest model run‐to‐run variability was selected to ensure easier comparison between the modalities.

2.4 Network feature analysis

To compare the RSNs of the two modalities, we had to begin by extracting their static and dynamic network features. These features included PSDs, power maps, and FC networks, which were obtained by either averaging over the entire recording duration or computing state‐specific descriptions based on the state time courses (Figure 1). These features were qualitatively compared across modalities and evaluated for their efficacy in capturing age‐related effects at the group level. Additionally, for the age classification task (see Section 2.7), we specifically computed TDE‐PCA covariance matrices and summary statistics as input network features. The network feature computations outlined in this subsection were largely adopted from Gohil et al. (2022) and conducted using the OSL Dynamics toolbox (Gohil et al., 2023).

When interpreting these network features, it is important to note that the approach we used for source reconstruction (i.e., unit noise gain invariance LCMV beamformer) implicates the loss of absolute amplitude information from the sensor recordings. This step was necessary to address the depth bias in source reconstruction, but it precludes us from making any direct comparison of absolute network features across modalities. Furthermore, we standardised (i.e., z‐transformed) the data across subjects before computing network features (see below), which eliminates amplitude information across individuals. This procedure allows us to compare the relative network features in each frequency band across modalities and subjects. Therefore, all network features employed in this article are measures relative to the time‐average for each sensor and subject, not absolute measures.

2.4.1 General linear model

We employed a group‐level general linear model (GLM) to study age‐related effects and measure the differences in various network features between the young and old groups. The value of a network feature (e.g., power for a parcel, or FC for a parcel‐pair, at a particular frequency) as it varied of subjects was fed in as the dependent variable to the group GLM, and the GLM was fitted independently to each network feature. The group‐level GLM design matrix (Figure A2) contained two regressors that modelled the group means of the young and old cohorts separately. Sex and head size were included as covariates to remove any confounding effects introduced by these variables. In particular, the GLM contained two sex regressors, one for male and the other for female, and one head size regressor, all of which were demeaned prior to model fitting. From the fitted GLMs we computed the mean group differences as a contrast of parameter estimates (COPEs). Group‐averaged values of each network feature were simply taken from the regression estimates of the fitted GLMs. For a comprehensive overview of the GLM and its application to neuroimaging data, readers should consult Quinn et al. (2024).

2.4.2 Power spectral densities

First, we computed subject‐level PSDs to examine the distribution of power across distinct frequency ranges. For the static analysis, we employed Welch's method (Welch, 1967) on the original parcel time course to calculate static PSDs of size N×M×P, where N is the number of subjects, M is the number of parcels, and P is the number of frequency bins. Each PSD was calculated using a 2‐s Hann window with 50% overlap across a 1.5–45 Hz frequency range, and each parcel time course was standardised (to have zero mean and unit variance) at the individual subject level before calculation. To compare the static PSDs between the young and old cohorts, PSDs were fitted to a GLM, producing the group‐averaged PSDs of size M×P.

For the dynamic analysis, the power spectrum for each subject, region, and state was separately calculated using the multi‐taper method, which is the typical approach used with TDE‐HMMs (Vidaurre et al., 2016). As in the static analysis, the data were standardised at the individual subject level. PSDs for each state were estimated using the source data and inferred state time course, which helped identifying time points corresponding to the active states. For each subject and state, a multi‐taper spectrogram was generated over the entire time points during which a single state was activated, using seven tapers with a time‐half bandwidth of 4. By averaging a spectrogram over time, state‐specific multi‐taper spectra (i.e., dynamic PSDs) of size N×M×P were obtained. In addition, state‐specific coherence spectra of size N×M×M×P were calculated using the multi‐taper cross‐spectral densities. For visualisation purposes, the mean PSD across all states, weighted by the duration of state activations, was subtracted from the state‐specific multi‐taper spectra at the subject level, and the group‐averaged dynamic PSDs were computed by averaging over the subjects in each group.

2.4.3 Power maps

Power maps were generated by averaging PSDs across specified frequency ranges of 1.5–20 Hz (wide‐band), 1.5–4 Hz (delta), 4–8 Hz (theta), 8–13 Hz (alpha), and 13–20 Hz (beta). These frequency bands were manually defined, rather than by using the non‐negative matrix factorisation method employed in Gohil et al. (2023), due to the failure of EEG data to converge during factorisation. For static analyses, power maps were derived from static PSDs; for dynamic analyses, we used dynamic, state‐specific multi‐taper spectra. Mean group differences in static power maps were quantified as COPEs by fitting the power maps to GLMs. For visualisation purposes, power values across brain regions were voxel‐weighted and interpolated. No threshold was applied to the plots. In dynamic power maps, the mean power map across all states, weighted by the duration of state activations, was subtracted from the state‐specific power maps at the subject level.

2.4.4 FC networks

FC was measured as the amplitude envelope correlation (AEC) for static analyses and coherence for dynamic analyses. The preference for AEC in static analyses was due to the method's known consistency in capturing stationary connectivity measures with relatively high test–retest reliability (Colclough et al., 2016). Amplitude envelopes were computed by applying the Hilbert transform to the source data, which were band‐pass filtered between 1.5 and 20 Hz and subsequently standardised. The static AEC networks were estimated by taking pairwise correlations of these amplitude envelopes across brain regions for each subject.

Conversely, the dynamic coherence networks were derived by averaging a coherence matrix for each state across 1.5–20 Hz. Note that it is conventional to use AEC for the static analysis and coherence for the dynamic analysis. Since the HMM infers coherence networks given the TDE data, using coherence is better suited to the dynamic case (Gohil et al., 2023). Moreover, applying AEC to dynamic FC usually involves segmenting signals according to state time courses and concatenating these segments. This process, however, can introduce edge effects and phase discontinuities, which may be detrimental to the Hilbert transform and its amplitude envelope. On the other hand, we opted not to use coherence for static analyses because it is reported to perform poorly in static connectivity estimation, showing low test–retest reliability (Colclough et al., 2016).

A threshold corresponding to the 95th percentile was applied to the FC networks for visualisation. In dynamic analyses, FC networks averaged across all states, weighted by the duration of state activations, were subtracted from the state‐specific FC maps.

2.4.5 Summary statistics of network dynamics

In the age group classification task (outlined in Section 2.7), we incorporated summary statistics detailing the dynamics of RSNs as part of the input features for our classifier. Four summary metrics were computed to quantify the temporal characteristics of the HMM states for each subject (Baker et al., 2014; Gohil et al., 2023). These metrics were calculated using the state time courses and include: Fractional occupancy: the proportion of the total time spent in each state.

Mean lifetime: the average amount of time spent in each state before transitioning into another state.

Mean interval: the average amount of time elapsed between visits to the same state.

Switching rate: the number of state visits per second.

2.4.6 TDE‐PCA covariance matrix

Another input feature we used for the classification task was static and dynamic TDE‐PCA covariance matrices. These features summarise the spatiotemporal properties of the source data. Subject‐specific static TDE‐PCA covariance matrices were generated by calculating the covariance of the TDE parcel data for each subject, performing a PCA to reduce the data to a rank of 80 and then standardising the resulting PCs. Note that this corresponds precisely to the steps taken in the data preparation to compute what we refer to as the standardised PCA‐reduced TDE data (cf., Section 2.3.1).

The HMM explicitly infers dynamic TDE‐PCA covariance matrices. However, these values are group‐level estimates. To calculate subject‐specific estimates of the dynamic TDE‐PCA covariance, we used the dual‐estimation method (Huang et al., 2024; Vidaurre et al., 2021), where we fixed the state probabilities in the group‐level HMM and re‐estimated the state TDE‐PCA covariance values using the subject‐wise prepared source data.

2.5 Statistical analysis

Having extracted the static and dynamic RSN features, we next examined the statistical differences between the two age groups to probe how age‐related effects are represented by each modality. We employed non‐parametric permutation tests for all pertinent brain network features, using t‐statistics as the permuted statistic (Holmes et al., 1996; Nichols & Holmes, 2002). For PSDs, we determined group‐level significance using a two‐tailed cluster permutation test. Clusters were formed over frequencies, with a cluster‐forming threshold set to 3. For power maps, we utilised a two‐tailed max‐t permutation test to evaluate group‐level significance.

For both types of permutation tests, a metric of interest was first fitted to a group‐level GLM as previously described. From the fitted GLMs, we computed COPEs (i.e., mean group differences) and their variances, from which t‐statistics could be calculated (Quinn et al., 2024). The group regressors were permuted 5000 and 10,000 times to construct a null distribution for the cluster and max‐t permutation tests, respectively. We then selected a significance threshold from this null distribution for statistical testing.

Due to the multiplicity of different non‐parametric permutation tests, we adjusted for multiple comparisons across frequency bands and states using Bonferroni correction. For the number of bands and states we accounted for, please refer to the figure legends for further details. This correction was applied within the MEG or EEG dataset; tests repeating across the modalities were not treated as part of the same statistical family.

2.6 Reproducibility of static and dynamic RSNs

To gauge whether the static RSNs derived from the EEG and MEG data are reproducible, we split each dataset into two halves and inferred networks on each half separately (Gohil et al., 2022). For each half, we compared static PSDs, power maps, and FC networks both within and across modalities. PSDs were averaged across subjects and parcels, while power maps and FC networks were averaged over a frequency range of 1.5–20 Hz and across subjects.

In Section 2.3.4, it was noted that the HMM is prone to converging to local optima, resulting in the model potentially arriving at similar variational free energies with different parameter estimates. Therefore, in addition to the initialisation method and the selection of the optimal runs, we implemented an additional validation step that assessed reproducibility across split‐half datasets. Specifically, the LEMON and Cam‐CAN datasets were each divided into two halves, with separate HMMs trained on each subset. For each half, a model was trained 10 times.

We assessed the reproducibility of the network dynamics through three measures. First, we examined the spatial distribution of power within and across the modalities by comparing the state‐specific power maps derived from each of the split‐half data. Next, we quantified the similarity of group‐level state covariance matrices both within and across modalities by computing RV coefficients for pairs of matrices (Robert & Escoufier, 1976). Finally, we evaluated the consistency of inferred transition probability matrices within and across modalities using the Jensen–Shannon (JS) distance (Lin, 1991). This metric was chosen for its symmetric nature and finite values, making it suitable for comparing probabilistic distributions. The JS distance was calculated for each row of transition probability matrices, and these values were averaged over all rows. This approach treats each row as a distinct probability distribution, viewing the transitions from one state to every other as separate entities and considering the overall similarity across these distributions.

It should be noted that power maps and RV coefficients were calculated using the optimal HMM run with the lowest variational free energy. The JS distances, on the other hand, were computed using all 10 runs. For each comparison within and across modalities, a distribution of JS distances was obtained from permuted pairs of 10 HMM runs across two split‐half datasets. We additionally tested whether the JS distances within EEG and MEG are drawn from the same distribution using the Wilcoxon rank‐sum test, after checking the assumptions of normality and homogeneity of variances.

2.7 Age group classifications

Finally, to evaluate the degree to which the network descriptions derived from MEG and EEG are meaningful, we assessed their utility in predicting age groups by employing a logistic regression classifier with L2 regularisation (Hoerl & Kennard, 1970; Ng, 2004). The task incorporated three distinct categories of input features: static, dynamic, and a combination of both. Static features comprised static TDE‐PCA covariances, while dynamic features included dynamic TDE‐PCA covariances alongside summary statistics representing network dynamics. These input features were concatenated before training a classifier.

Separate classifiers were trained for each category of input features with a nested cross‐validation approach, which involved two levels of cross‐validation. In this framework, the inner cross‐validation focused on hyperparameter optimisation, while the outer cross‐validation estimated the generalisation error across multiple test sets. In the outer loop, input features were randomly shuffled and subjected to a fivefold cross‐validation. The training set partitioned from this outer loop underwent another fivefold cross‐validation within the inner loop. For the inner cross‐validation, we standardised the data and applied dimensionality reduction (i.e., PCA with whitening) before fitting any classifiers to enhance predictive accuracy. The optimal number of PCs was determined via a grid search on the inner training folds over 10, 20, and 30 PCs.

The best‐performing classifier from the inner loop was then selected and evaluated on the test set in the outer loop to compute its overall predictive accuracy. Here, the test set was processed using the same standardisation and PCA parameters fitted on the training set to ensure consistency across data splits and maintain fair model evaluation. This entire nested cross‐validation procedure was repeated 10 times. We calculated and reported the mean and standard deviation of the accumulated test accuracies to provide a comprehensive overview of the classifier's performance on the M/EEG network features.

To determine if the test accuracies exceeded the level anticipated by random chance alone (i.e., 50%), a one‐tailed, one‐sample t test was used. We employed a two‐tailed, two‐samples t test to assess whether there are significant differences in the prediction scores between MEG and EEG, as well as to examine if the prediction scores varied depending on how the data were reconstructed (i.e., the use of subject sMRI data in source reconstruction). The significance of all tests was evaluated against a Bonferroni‐corrected threshold, adjusted for three categories of input features. The assumptions of normality and homogeneity of variances were verified prior to conducting these tests. When the assumption of equal variances was violated, we employed Welch's t test as an alternative to the standard two‐sample t test.

2.8 Use of artificial intelligence assistance

While preparing this article, the authors used ChatGPT to refine the styles and grammar of the parts of the written text. This tool was never used to create and write new content. The authors reviewed and edited any changes made and take full responsibility for the content of the publication.

2.9 Data and code accessibility

For the original LEMON and Cam‐CAN datasets, the readers should refer to Babayan et al. (2019) and Shafto et al. (2014), respectively. All codes for analyses and data visualisations, including scripts for preprocessing, source reconstruction, and model training, were written in Python and shared at https://github.com/OHBA-analysis/Cho2024_MEEG_RSN. Some of the scripts were built upon the OSL Dynamics software (Gohil et al., 2023). The inferred HMM parameters and computed RSN features are shared in the same repository. Further requests for prepared data or concerns regarding the codes should be directed to and will be fulfilled by the authors.

3 RESULTS

3.1 Static analysis

3.1.1 Static RSN features are qualitatively similar across the modalities

As discussed in the Introduction, the comparability between static functional brain network features of MEG and EEG has been established by previous literature. To ascertain this similarity with a medium‐density EEG (61 channels), we examined how static RSN features qualitatively compare between these modalities (Figure 2). The subject‐averaged, wide‐band (1.5–20 Hz) power maps for EEG and MEG both revealed pronounced activation in the frontal, parieto‐occipital, and sensorimotor cortices. MEG exhibited a higher signal‐to‐noise ratio, with its power values higher than EEG. The subject‐averaged, wide‐band FC maps also indicated similarities between the modalities, in which the connections were strongest at the parieto‐occipital region. Here, EEG showed stronger connection strengths, whereas MEG FC was more concentrated on the parietal cortex. Likewise, M/EEG PSDs averaged across subjects and parcels displayed comparable spectral distributions with a prominent unimodal peak in the alpha frequency band. However, their aperiodic activities diverged, as EEG showed a steeper decay following a 1/f pattern. All these trends were also evident in the M/EEG data that had been reconstructed using a standard brain structural image instead of subject‐specific sMRI images (Figure A3).

FIGURE 2 Qualitative comparison of static resting‐state networks between EEG and MEG. (a) Brain surface map of static power for each brain parcel, calculated by first averaging PSDs across the frequency range of 1.5–20 Hz at the individual level and then across all subjects (left). On the top right, an FC network of static AEC averaged across subjects, with only the top 5% of correlation strengths shown for better visualisation. AEC values were calculated on envelopes obtained from subject‐wise parcel time courses, band‐pass filtered from 1.5 to 20 Hz. On the bottom right, the static PSD averaged across both brain parcels and subjects, with standard errors across parcels indicated by grey shading. (b) The subject‐averaged static power map, functional connectivity network, and PSD generated using the MEG data following the same analyses as in (a).

To ascertain whether the static RSNs derived from M/EEG and their overall comparability are reproducible, we extracted the same static PSDs, power maps, and FC networks from the EEG and MEG datasets divided into halves. We then assessed their reproducibility within and across the modalities (Figure 3). A qualitative comparison illustrated that static spectral, power, and connectivity distributions remained consistent within each modality. In both EEG and MEG, each half highlighted identical brain regions with high power and connectivity values, while also demonstrating similar PSD shapes (Figure 3a,b). Across the modalities, any comparisons of two halves mirrored the observations from the full dataset. That is, MEG showed higher power than EEG, whereas EEG exhibited stronger FC than MEG, with MEG accentuating the parietal cortex relatively more. To this extent, the comparability of the static RSN features was higher within each modality than between the two modalities.

FIGURE 3 Resting‐state network reproducibility of static network features within and across modalities. (a) The subject‐averaged, wide‐band (1.5–20 Hz) power maps and FC networks, along with the PSDs averaged across both brain parcels and subjects, computed from the first half (top panel) and the second half of the resting‐state EEG data (bottom panel). Each half comprises data from 48 subjects. The power maps depict lateral surfaces at the top and medial surfaces at the bottom, and the FC networks show only the top 5% of correlation strengths. For the PSDs, standard errors across parcels are indicated by grey shading. (b) The same analysis as in (a) is conducted using the MEG data.

In summary, EEG and MEG demonstrated qualitative correspondence across all static network features, despite some discrepancies. This comparability was noted in the data that were reconstructed using a standard brain structural image as well, with which we could observe the same discrepancies between the two modalities. These network features also displayed a good level of reproducibility, both within and across the modalities, when qualitatively contrasted.

3.1.2 MEG and EEG reveal comparable age effects in static network features

To investigate whether MEG and EEG provide equally meaningful static network descriptions, we examined age‐related effects captured by their static RSN features. Both modalities demonstrated comparable age‐related effects in PSDs, identifying statistically significant differences between the age groups across the delta, theta, and beta bands (Figure 4a,f). However, slight discrepancies were noted at higher frequencies (30–45 Hz). The topographic maps highlight the alpha peak PSD differences, the power differences in lower frequencies, and the power differences in the beta band; these topographical patterns additionally showed qualitative similarities between the two modalities.

FIGURE 4 EEG and MEG reveal comparable age‐related effects in static network features. (a) The between‐group PSD difference (red dotted line), averaged over parcels and subjects, is depicted alongside green lines denoting the subject‐averaged PSD differences for each parcel (the inset in the top‐right corner indicates the shade of green used for each brain parcel) of the EEG data. Frequency ranges with significant PSD differences (p < .05; horizontal red bars) are identified via a cluster permutation test on parcel‐averaged PSDs of the two age groups. Additionally, the top topographical maps show the peak/trough group‐level PSD differences found in the alpha frequency band. The bottom topographies show group‐level PSD differences integrated separately across the lower (1.5–8 Hz) and beta (13–20 Hz) frequency bands, represented by the grey shading. (b) On the left, a brain surface map illustrating delta band power differences between the old and young groups is plotted for the EEG data, with lateral surfaces at the top and medial surfaces at the bottom. On the right, the max‐t permutation test (p < .05; Bonferroni‐corrected, n = 4 frequency bands) was applied to the power map to highlight the age‐related effects, with significant parcels accentuated. Red and blue colours indicate higher power in old and young participants, respectively, and non‐significant findings are labelled as n.s. The same analysis as in (b) is shown for the (c) theta, (d) alpha, and (e) beta band EEG power maps. (f) The same analysis as described in (a) is shown for the MEG data. (g–j) The same analysis as described in (b) is shown for the delta, theta, alpha, and beta band power differences, respectively, using the MEG data.

Further analysis of age‐related effects in narrow‐band (delta, theta, alpha, beta) power spatial maps indicated that the general trend of power differences between the age groups (whether young participants exhibited higher power than old participants, or vice versa) was statistically significant across two modalities (Figure 4b–e,g–j; right). Nonetheless, their specific spatial distributions still demonstrated slight discrepancies. In the delta power map, the young cohort exhibited higher power in the frontal and temporal cortices, and specifically in the sensorimotor region in MEG (Figure 4b,g). In the theta power map, the young group showed increased power in the frontal and sensorimotor regions, with only EEG displaying higher power in the parietal region (Figure 4c,h). Significant age‐related alpha effects were exclusively observed in MEG within the occipital cortex, whereas EEG did not identify any regions with notable age effects (Figure 4d,i). In the beta power map, age‐related effects were more prominent in MEG, which identified effects in the parietal, frontal, and sensorimotor cortices, albeit both modalities recognised effects in the parietal region (Figure 4e,j). The general trend of between‐group static power differences was more evident in the qualitative representation of age effects (Figure 4b–e,g–j; left).

The observed age‐related effects remained consistent in the M/EEG data that were reconstructed using a standard brain structural image (Figure A4). However, notable differences were identified, including minor deviations in spatial patterns within the parietal and temporal cortical regions in the EEG beta power map (Figure A4e), along with the absence of age‐related effects in the MEG alpha power map (Figure A4i).

In summary, although some of the spatial details of the age‐related effects varied, MEG and EEG often captured group‐level differences of a similar nature. Their comparability persisted even when the data were reconstructed with a standard brain structural image. Nonetheless, there was a tendency for MEG to be slightly more sensitive to age‐related changes.

3.2 Dynamic analysis

3.2.1 Dynamic RSN features are qualitatively similar across MEG and EEG

Following the static analysis, we examined how dynamic RSN features compared between the two modalities. Figure 5 presents the subject‐averaged PSDs and wide‐band (1.5–20 Hz) power map and FC network for each state, demonstrating their alignment between EEG and MEG. Notably, the spatial patterns of power distributions showed a marked resemblance between the modalities. Nonetheless, distinctions arose particularly in states 2 (the visual network) and 6 (the sensorimotor network), where there was relatively less correspondence. These states exhibited divergent areas of stronger activation or deactivation between the modalities; for instance, state 2 displayed elevated occipital activity in EEG, while more pronounced deactivation in non‐occipital areas was observed with MEG.

FIGURE 5 Qualitative comparison of dynamic resting‐state networks between EEG and MEG. (a) Each box shows the subject‐averaged power map (left), FC network (middle), and parcel‐averaged PSD (right) of each HMM state for 96 EEG subjects. Each state is considered as a distinct RSN. The power maps show lateral surfaces at the top and medial surfaces at the bottom. The FC networks illustrate connections with the top 5% coherence values (irrespective of sign). The shaded areas of the PSDs represent the standard error of the mean across the parcels. The power maps, FC networks, and PSDs are visualised relative to their average across all states. (b) The plots follow the same format as (a), showing the RSNs computed from 96 MEG subjects.

Similarly, the state‐specific FC networks and PSDs largely maintained a high degree of comparability between EEG and MEG, alongside more detailed discrepancies. Such discrepancies were most pronounced in states 1 (a putative anterior default mode network) and 2 (the visual network), in which the top 5% coherence values and PSDs highlighted different brain regions and spectral patterns, respectively, when visualised relative to the average values across all states.

The comparability of dynamic network features across the two modalities was reflected in the MEG and EEG RSNs derived from the data reconstructed using a standard brain structural image (Figure A5). Importantly, the discrepancies previously observed persisted, with states 2 and 6 displaying reduced correspondence in power maps and states 1 and 2 showing lesser alignment in FC networks and PSDs.

A possible source of the between‐modality discrepancies might stem from variability among subjects and between different model runs. To ensure that the overall consistency of the dynamic RSNs derived from EEG and MEG are replicable despite these minor differences, we evaluated the reproducibility of HMM states within and across the modalities (Figure 6). We started by comparing the state‐specific power maps generated from the EEG and MEG datasets split into two halves. Our qualitative comparison revealed that spatial power distributions remained consistent within each modality (Figure 6a,b). There was also a notable similarity in power maps across the modalities, with state 6 being the sole exception to this trend.

FIGURE 6 Resting‐state network reproducibility of dynamic network features within and across modalities. (a) The wide‐band (1.5–20 Hz) power maps for resting‐state EEG data, averaged across subjects, are shown for the first half of the dataset (top) and the second half of the dataset (bottom), with each half comprising data from 48 subjects. These maps show lateral surfaces at the top and medial surfaces at the bottom, visualised relative to their average across all states. (b) The same analysis as in (a) is conducted using the MEG data. (c) For each state, the similarities between inferred covariance matrices within and across modalities are measured as RV coefficients (higher value indicates more reproducible). (d) Similarities between inferred transition probability matrices within and across modalities are measured as the Jensen–Shannon (JS) distances and depicted in bar graphs (lower value indicates more reproducible). The statistical difference between JS distances within EEG and MEG was tested using the Wilcoxon rank‐sum test (U = 9.053, p = 1.394e‐19). Significance levels are indicated by ***: p < .05, **: p < .01, *: p < .001.

The reproducibility of HMM states was also measured with respect to the inferred HMM parameters. In particular, we computed the RV coefficients, a generalisation of the squared Pearson correlation coefficient, between the covariance matrices of two split‐half M/EEG datasets (Figure 6c). As anticipated, RV coefficients for comparisons within the same modality were typically higher (more reproducible) than those for comparisons across the modalities. For all states except state 6, RV coefficients exceeded 0.65 for any given comparison, although the variability of this measure over different types of comparisons was relatively higher for states 1 (a putative anterior default mode network) and 6 (the sensorimotor network). Notably, RVs within MEG were generally higher than those within EEG.

Subsequently, we quantified the difference between transition probability matrices of two split‐half datasets as the JS distance (Figure 6d). This distance was lower (more reproducible) for the comparison within MEG than that within EEG or across the modalities. In particular, the JS distances within MEG were significantly smaller than those within EEG. The comparisons across the modalities had median JS distances similar to the median value within EEG, although the latter exhibited a larger variance, including relatively lower distances. Such similarity partly serves as additional evidence that EEG and MEG RSNs are comparable as they pertain to inferred state transition probabilities.

In summary, EEG and MEG were broadly similar across the dynamic network features, while also showing some discrepancies. Their comparability, as well as discrepancies, persisted when the data were reconstructed using a standard brain structural image. Furthermore, we measured the reproducibility of these network features, as well as of the HMM parameters, both within and across the modalities. This showed that state 6 (the sensorimotor network) was the least reproducible and that EEG was less reproducible than MEG.

3.2.2 MEG and EEG reveal comparable age effects in dynamic network features

The age‐related effects in dynamic PSDs and power maps derived from MEG and EEG are illustrated in Figure 7. Note that for each dynamic measure we compute the “dynamic contribution to the state‐specific” measure, resulting in any dynamic effects above the static effects. For example, the power in a state in a particular frequency band is the power in that state in addition to the static power in that frequency band.

FIGURE 7 EEG and MEG reveal comparable age‐related effects in dynamic network features. (a) The state‐specific, dynamic PSDs of the EEG data averaged across grouped subjects and parcels, depicted for old (dotted) and young (solid) participants. For each HMM state, a cluster permutation test is conducted on parcel‐averaged PSDs to detect age‐related dynamic effects, with frequencies exhibiting statistical significance (p < .05; Bonferroni‐corrected, n = 6 states) highlighted in red. An inset at the top‐right corner displays the t‐statistics, signifying the magnitude of age‐related effect sizes between the groups. () The state‐specific power maps with significant between‐group differences for the EEG data. These differences are validated through max‐t permutation tests (p < .05; Bonferroni‐corrected, n = 6 states), with significant parcels accentuated. Red and blue colours indicate higher power in old and young participants, respectively, and non‐significant findings are labelled as n.s. (c) The same analysis as in (a) performed on the MEG data, marking significant age‐related dynamic effects in state‐specific PSDs. (d) The same analysis as in (b) performed on the MEG data, illustrating significant age‐related dynamic effects in state‐specific wide‐band power maps. Note that in all cases we plot the “dynamic contribution to the state‐specific” feature, resulting in dynamic effects that are above the static effects.

For both modalities, the age‐related effects for the dynamic PSDs (Figure 7a,c) broadly appeared in the same frequency bands as for the static PSDs. Specifically, (Figure 7a,c) shows significant spectral differences between the age groups in the delta and theta bands across all states for both MEG and EEG. However, consistency in these differences was less evident in the beta band and higher frequencies. While states 4, 5, and 6 exhibited overlapping effects in these frequencies, states 1, 2, and 3 did not. In state‐specific, wide‐band (1.5–20 Hz) power maps, all states except state 1 displayed roughly equivalent age‐related effects, with the young cohort showing higher power compared to the old group (Figure 7b,d). However, states 2, 4, and 6 presented overlapping patterns across brain regions between the two modalities, whereas states 3 and 5 revealed significant power differences in regions that varied across the modalities.

The observed age‐related effects remained consistent in the M/EEG data that were reconstructed using a standard brain structural image without subject sMRI images, although to a lesser degree compared to the static analysis (Figure A6). In particular, the age effects in higher frequencies altered in the PSDs of EEG states 2 and 4 (Figure A6a). The differences were more pronounced on state‐specific power maps. In general, significant power spatial distributions exhibited different patterns when a standard structural image was used. While EEG and MEG state 1 highlighted significant age‐related effects in the occipital and parietal cortices, respectively, significance in MEG state 6 diminished under this condition (Figure A6b,d).

In summary, MEG and EEG demonstrated broadly comparable age‐related effects in dynamic network features (specifically, the dynamic effects above the static effects), albeit to a lesser extent compared to static features. This comparability persisted when the data were reconstructed using a standard brain structural image, although it was not maintained at a level similar to the static analysis.

3.3 Age group classification task

Our comparative analysis of how MEG and EEG capture age‐related effects in static and dynamic RSN features showed that both modalities offer broadly comparable network descriptions. To extend this analysis further in a quantitative manner, we aimed to estimate the predictive power of these network features in accurately classifying subjects into distinct age groups. Utilising three distinct sets of input features—static, dynamic, and a combination of both—we employed a logistic regression classifier with an L2 regularisation penalty (see Section 2.7).

We evaluated prediction scores against test datasets across 10 repeated classification tasks, with Figure 8 illustrating the resulting accuracies in age prediction. Both modalities showed predictions that exceeded what would be expected by random chance alone (i.e., 50%); the mean prediction accuracies, averaged across 10 repetitions, were significantly higher for MEG than EEG when employing static and dynamic input features (Figure 8a). Although this statistical difference was not observed with combined features, there was still a clear trend towards better prediction for MEG compared to EEG.

FIGURE 8 Predictive accuracy in age group classification using EEG and MEG resting‐state network features. (a) Mean prediction scores for three different RSN feature types derived from the EEG (purple) and MEG (pink) data across 10 repeated runs of classification tasks. The bars denote mean scores, with error bars reflecting standard deviations. A dotted line represents the score expected by random chance. Features were obtained from data reconstructed using individual sMRI images. Prediction scores were evaluated against random chance with a one‐sample independent t test (for EEG, static: T = 19.01, p = 7.085e‐9; dynamic: T = 24.81, p = 6.736e‐10; combined: T = 30.11, p = 1.203e‐10; for MEG, static: T = 44.05, p = 4.003e‐12; dynamic: T = 26.29, p = 4.027e‐10; combined: T = 28.71, p = 1.840e‐10) and compared between EEG and MEG using a two‐sample independent t test (static: T = −18.47, p = 3.801e‐13; dynamic: T = −10.53, p = 4.016e‐9) and a Welch's t test (combined: T = −12.99, p = 8.974e‐9). (b) Mean prediction scores for three distinct RSN feature types derived from the EEG data that were source reconstructed with (grey) and without (orange) individual sMRI images. Bar and error bar representations match (a). Differences due to sMRI image use were assessed with a two‐sample independent t test (dynamic: T = 3.156, p = 5.471e‐3; combined: T = 4.076, p = 7.088e‐4). (c) The same analysis as in (b) using the MEG data (static: T = 3.320, p = 3.807e‐3). Significance levels are indicated by */♦: <0.05, **/♦♦: <0.01, ***/♦♦♦: <0.001 (before Bonferroni‐corrected, n = 3 features), with non‐significant findings unmarked.

Next, we conducted the same analysis using the data that were reconstructed without subject‐specific sMRI images. For EEG, the prediction scores significantly decreased when using dynamic and combined features, while a significant decrease was observed with static features for MEG (Figure 8b,c). Irrespective of statistical significance, however, all three types of input features exhibited lower prediction scores for both EEG and MEG when only the standard MNI template was provided for reconstruction.

4 DISCUSSION

Using the MEG RSNs from the Cam‐CAN dataset as a benchmark for evaluating brain networks of fast oscillatory activity, we investigated EEG‐derived RSNs to ascertain if both modalities can yield meaningful network descriptions and establish how comparable they are to one another. The results indicate qualitative and quantitative comparability between MEG and EEG in representing functional brain networks and suggest some advantages of HMM over traditional network inference methods for EEG (e.g., microstates (Michel & Koenig, 2018), sliding window method (de Pasquale et al., 2010)). Our findings highlighted several key points.

First, consistent with previous research (Coquelet et al., 2020; Siems et al., 2016), we found that MEG and EEG exhibit analogous static RSN features (Figure 2), which were reproducible within and across the MEG and EEG datasets (Figure 3). Age‐related effects in static PSDs and power maps also proved comparable across the modalities. Nonetheless, there was a tendency for MEG to be more sensitive to age‐related changes than EEG (Figure 4). These network descriptions and age effects persisted even when the data were reconstructed using a standard brain structural image (Figures A3 and A4).

Second, MEG and EEG demonstrated broadly comparable dynamic RSN features (Figure 5), with age‐related effects in state‐specific PSDs and power maps displaying marked similarity, albeit to a lesser degree compared to the static RSNs (Figure 7). These network features and age‐related effects remained largely consistent when a standard brain structural image was used to reconstruct the data (i.e., when subject‐specific sMRI images were not available) (Figures A5 and A6). We also used a split‐half approach to quantitatively measure the extent to which the dynamic network features were reproducible within and across the MEG and EEG datasets (Figure 6). Overall, this approach revealed a good level of reproducibility, with state 6 (the sensorimotor network) found to be the least reproducible. We deemed a network description reproducible if the network features qualitatively resembled one another across all states, or at least the vast majority of states. We found that EEG was generally less reproducible than MEG. In sum, the similarity of dynamic RSNs across the modalities was substantial when contrasted with a previous sliding‐window‐based analysis (Coquelet et al., 2020) that reported discordant dynamic RSN FC between MEG and EEG.

Finally, we confirmed that both MEG and EEG can provide meaningful static and dynamic network descriptions by testing the efficacy of M/EEG network features in predicting age groups. Although EEG reported lower test accuracy than MEG, the performance scores of both modalities were above random chance expectations (Figure 8a). However, one caveat is that, while MEG outperformed EEG in the age class prediction task, differences between the two datasets may not solely reflect differences between EEG and MEG. Other factors, such as age range and participant heterogeneity, may have also influenced the results.

Interestingly, the predictive power was similar for the static and dynamic RSN features and did not significantly improve when they were combined. This lack of improvement suggests that, for the simple task of classifying binary age groups with a logistic classifier, using either static or dynamic network descriptions suffices. The complex analysis combining static and dynamic features does not appear to be beneficial for this specific task, possibly due to the overfitting of the data. In addition, these test accuracies varied across different brain spaces. For EEG, performance metrics for dynamic and combined input features showed significant variation depending on whether the data were reconstructed with subject‐specific sMRI images or the standard MNI template (Figure 8b). For MEG, significant variation was observed in performance for static input features (Figure 8c). Although not always statistically significant, performance metrics for all three input feature types decreased for both EEG and MEG when we used the standard MNI template.

4.1 Basis set of EEG RSNs

As noted at the outset of this article, the current study provides a basis set of RSNs derived from the EEG LEMON dataset, alongside their counterparts from the MEG Cam‐CAN data. This basis set was constructed using medium‐density EEG, thereby rendering it as data applicable for subsequent studies employing similar EEG systems. In particular, we offer scripts designed for the preprocessing, source reconstruction, and analysis of M/EEG data, as well as the inferred model parameters of the optimal HMM runs (e.g., state covariances, transition probabilities, time course of state probabilities). These data and resources are made publicly accessible, with detailed guidelines for software installation and data downloading available in the GitHub repository (see Section 2.9).

The basis sets of static and dynamic EEG RSNs are illustrated in Figures 2 and 5, respectively. The designation of these EEG networks as a basis set is justified for several reasons. Notably, in both static and dynamic analyses, EEG RSNs not only resembled those derived from the MEG data but also demonstrated age‐related effects comparable across the two modalities. In addition, their network features were reproducible within and across the modalities, although with the exception of state 6 in the dynamic RSNs. Finally, akin to MEG, network descriptions in the EEG basis set demonstrated robust performance in age prediction, significantly surpassing that of random chance.

Despite the comparability between EEG and MEG, it is important to acknowledge that MEG demonstrated better reproducibility, enhanced predictive accuracy, and greater sensitivity to age‐related effects compared to EEG. This discrepancy is arguably anticipated, given the higher sensor count of MEG and consequently its superior spatial resolution. Previous studies (Coquelet et al., 2020; Knyazev et al., 2016; Liu et al., 2017) suggest that EEG datasets with a greater number of sensors can also achieve a level of comparability with MEG data. Nonetheless, it remains to be elucidated how EEG data, varying in sensor count relative to the EEG LEMON dataset, compares with MEG, especially in terms of dynamic RSNs as derived from state‐of‐the‐art methodologies. Among several potential applications, this basis set could be used for investigating RSNs in high‐ and low‐density EEG systems, developing methods for effectively combining M/EEG with fMRI, and comparing basis sets of RSNs inferred from various modalities such as optically pumped magnetometers (OPM).

4.2 Limitations

One of the main limitations in our study was the relatively small sample size (n = 96 subjects). This limitation meant that while subjects across the two datasets could be age‐matched, further matching based on sex or Mini‐Mental State Examination scores was not feasible. Likewise, demographic attributes, recording sites, and other cognitive or physiological states that have not matched may introduced unaccounted variance into our effects of interest. In the context of our classification task, the small sample size restricted us to a binary categorisation of age groups (young vs. old), as it was difficult to subdivide our sample suitably under the k‐fold cross‐validation framework. Because ages were categorised in intervals within the LEMON dataset, predicting numerical, continuous age values had to be precluded from the outset. Other limitations included the use of random subsampling during the age‐matching process, which may have brought additional subject variability to our comparison between LEMON and Cam‐CAN.

It is also worth noting that the different data acquisition methods for EEG and MEG may have contributed to certain disparities in their RSN descriptions and associated age effects. The EEG data were recorded in 1‐min blocks alternating between eyes‐closed and eyes‐open conditions, while the MEG data were recorded with participants keeping their eyes closed throughout the session. Although the EEG analysis was limited to the eyes‐closed periods, the intermittent eye openings might have altered the brain state during EEG acquisition (e.g., by forcing participants to be more alert on average). On the contrary, the brain state in the second half of the MEG session may differ from the first half due to potential drowsiness. These differences in data acquisition likely impacted the underlying brain dynamics in EEG and MEG, leading to some disparities in the effects observed in the two modalities.

Another notable aspect of our preprocessing pipeline is that we used a parcellation atlas with 38 anatomical regions for reconstructing both EEG and MEG, which has previously been shown to work effectively on MEG data (Colclough et al., 2015). However, it is plausible that using fewer or coarser parcels could enhance the comparability between EEG and MEG. Even in this case, some differences could be preserved, for example, if localised age effects align with the adopted parcellation. Similarly, altering parcel shapes could also possibly increase comparability between the modalities if the parcel borders are defined in such a way that differences between EEG and MEG are reduced. Exploring the effects of different parcellation atlases on the EEG and MEG RSNs would be an intriguing direction for future research.

4.3 Outlooks and future works

Previous research has indicated marginal comparability between MEG and EEG in capturing the temporal dynamics of brain networks (Coquelet et al., 2020). They proposed that the reasons for this poor correspondence are primarily twofold. From a physiological standpoint, it was surmised that MEG and EEG detect distinct components of transient neural activities, particularly in fast oscillatory fluctuations above 1 Hz, originating from different neural populations responsible for generating magnetic and electric fields. From a methodological perspective, it was acknowledged that the sliding‐window approach, compounded by the lack of spatial leakage correction, could result in large statistical errors (Hutchison et al., 2013). Hence, the estimation of dynamic RSN FC at shorter timescales may not benefit from smoothing effects typically associated with static analysis techniques, such as low‐pass filtering or amplitude envelope analysis.

Although we found the methodological limitations plausible, our findings did not support the physiological explanation. When we applied the TDE‐HMM (Vidaurre, Hunt, et al., 2018), a more advanced dynamic modelling technique validated across various MEG datasets, we discovered a previously unobserved similarity between the modalities in both static and dynamic RSNs, even with a lower‐dimensional 61‐channel EEG system and a 38 parcellated brain regions. Consequently, whereas it was believed that MEG and EEG produce complementary views on the dynamic functional network organisations of the human brain, our results indicate that both modalities can identify the same network profiles, at least in the resting state.

Therefore, instead of assuming they measure distinct aspects of brain functionality (Cichy & Pantazis, 2017; Mahmutoglu et al., 2022), our data suggest that the same functional networks underpin brain activity, irrespective of the modalities. This indicates that the RSNs may be less a characteristic of the specific measurement modality and more a neural substrate correlated with cognition and behaviours that is common across neuroimaging modalities. It is conceivable that different modalities may highlight specific aspects of RSNs, as evidenced by statistical testing on age effects; however, the overarching spatial patterns remained consistent. With the application of robust statistical methods capable of enhancing spatiotemporal precision (McFadyen & Dolan, 2023), such as dynamic network modes (Gohil et al., 2022), we may be able to see more clearly the same RSNs in both MEG and EEG.

Our results thus support the notion that previously known properties of MEG and fMRI RSNs may translate into EEG. Furthermore, while both MEG and EEG provide unprecedented temporal resolution over fMRI, their comparability allow us to consider prioritising EEG in future functional brain network studies, given its cheap cost and broader accessibility. If high‐density EEG systems are leveraged in conjunction with advanced statistical modelling methodologies, we anticipate an even greater alignment in the dynamic network features identifiable across these modalities.

To this extent, future works should aim to verify whether cognitive, behavioural, and clinical phenotypes associated with MEG and fMRI RSNs can also be identified in EEG. Given that dynamic RSN phenotypes of neuropsychiatric disorders have predominantly been reported in fMRI and MEG studies, it would be crucial to explore dynamic EEG RSNs in this regard. The established correspondence between MEG and fMRI RSNs (Brookes et al., 2011) implies that biomarkers identified through these modalities could potentially be detected within EEG data. An additional avenue of research could focus on the use of simultaneous M/EEG recordings to more accurately assess how dynamic RSN features and the age‐related effects in them vary between modalities. This approach would offer insights into the efficacy and necessity of multimodal electrophysiological techniques. Furthermore, examining whether the comparability of RSNs across MEG, EEG, and fMRI extends to other measurement devices, such as wearable functional near‐infrared spectroscopy (Ferrari & Quaresima, 2012) or OPM‐based MEG systems (Brookes et al., 2022), would substantiate the validity of our findings.

5 CONCLUSION

In this article, we demonstrated that MEG and EEG produce comparable static and dynamic network descriptions, which were shown to be meaningful in highlighting age‐related effects and predicting age groups. MEG, however, was more reproducible and sensitive in detecting age‐related effects than EEG. This comparability and concomitant discrepancies between the modalities were also observed when the datasets were reconstructed using a standard brain structural image and without requiring subject sMRI images, albeit to a lesser degree. Our results suggest that MEG and EEG share the same functional brain networks, likely paralleling those detected by other modalities like fMRI.

AUTHOR CONTRIBUTIONS

SungJun Cho: Conceptualisation, data curation, software, methodology, formal analysis, writing—original draft, writing—review and editing. Mats van Es: Conceptualisation, software, writing—review and editing, supervision. Mark Woolrich: Conceptualisation, methodology, writing—review and editing, supervision. Chetan Gohil: Conceptualisation, data curation, software, methodology, writing—review and editing, supervision.

CONFLICT OF INTEREST STATEMENT

The authors declare no competing financial interests or personal relationships that could be perceived as influencing the findings reported in this article.

Supporting information

FIGURE A1. Age distributions of participants in the EEG LEMON and MEG Cam‐CAN datasets. Histograms depict the distributions of ages for the young (blue) and old (red) group. Note that the EEG and MEG data have identical distributional shapes, as the subjects in two datasets are age‐matched. To ensure comparable sample sizes between LEMON and Cam‐CAN, a subset of participants was sampled from each age group for each dataset. Both datasets have a sample size of 60 young and 36 old participants.

FIGURE A2. Group‐level design matrices and contrasts for the EEG and MEG datasets. (A) The design matrix for the EEG data with regressors in columns and individual subjects in rows. The matrix consists of the condition (Group) and covariate (Sex, Headsize) regressors, with covariate values demeaned. A colour bar indicates the values for each regressor. The table below outlines the contrast (i.e., mean group difference), denoted as C, including the weightings associated with each regressor. (B) The design matrix for the MEG data visualised in the same structure as in (A).

FIGURE A3. Static resting‐state networks inferred from datasets source reconstructed using a standard brain structural image. (A) A subject‐averaged, wide‐band (1.5–20 Hz) power map displayed for the EEG LEMON dataset (left). On the top right, a subject‐averaged, wide‐band (1.5–20 Hz) FC network, with only the top 5% of correlation strengths shown for visualisation purposes. On the bottom right, a parcel‐ and subject‐averaged PSD, with standard errors across parcels indicated by grey shading. (B) The subject‐averaged power map, FC network, and PSD are generated using the MEG Cam‐CAN dataset following the same analyses described in (A). For both EEG and MEG, the data were source reconstructed using a standard MNI brain, not individual sMRI files.

FIGURE A4. Age‐related effects in static network features extracted using a standard brain structural image. (A) For the EEG data, the between‐group PSD difference (red dotted line), averaged over parcels and subjects, and its standard error (pink shade) are depicted alongside green lines denoting the parcel‐wise, subject‐averaged PSD differences. An inset in the top‐right corner marks the brain parcel locations associated with each PSD difference line. Frequency ranges with significant (p < 0.05) PSD differences between the two age groups are indicated by red bars. Additionally, the group‐level peak differences in the alpha (8–13 Hz) band and group‐level PSD differences integrated across the lower (1.5–8 Hz) and beta (13–20 Hz) bands are depicted in topographical maps at the top and bottom, respectively. (B) On the left, a brain surface map illustrating delta band power differences between the old and young groups is plotted for the EEG data. On the right, the brain regions with significant (p < 0.05; Bonferroni‐corrected, n = 4 frequency bands) between‐group differences are highlighted. Red and blue colours indicate higher power in old and young participants, respectively, and non‐significant findings are labelled as n.s. The same analysis as in (B) was replicated for the (C) theta, (D) alpha, and (E) beta band EEG power maps. (F) The same analysis as described in (A) is conducted on the MEG data. (G)‐(J) The same analysis as described in (B) is applied to delta, theta, alpha, and beta band power differences using the MEG data, respectively. For both EEG and MEG, the data were source reconstructed using a standard MNI structural brain image.

FIGURE A5. Dynamic resting‐state networks inferred from datasets source reconstructed using a standard brain structural image. (A) Each box shows the subject‐averaged power map (left), FC network (middle), and parcel‐averaged PSD (right) of each HMM state for 96 EEG subjects. The power maps show lateral surfaces at the top and medial surfaces at the bottom. The FC networks illustrate connections with the top 5% coherence values. The shaded areas of the PSDs represent the standard error of the mean. The power maps, FC networks, and PSDs are visualised relative to their average across all states. (B) The plots follow the same format as (A), showing the RSNs computed from 96 MEG subjects. For both EEG and MEG, the data were source reconstructed using a standard MNI brain and not individual sMRI files.

FIGURE A6. Age‐related effects in dynamic network features extracted using a standard brain structural image. (A) The state‐specific PSDs of the EEG data are averaged across grouped subjects and parcels and depicted for old (dotted) and young (solid) participants. For each HMM state, frequencies exhibiting statistically significant (p < 0.05; Bonferroni‐corrected, n = 6 states) age‐related dynamic effects are highlighted in red. An inset at the top‐right corner displays the t‐statistics, signifying the magnitude of age‐related effect sizes between the groups. (B) The state‐specific power maps with significant (p < 0.05; Bonferroni‐corrected, n = 6 states) between‐group differences are shown for the EEG data. Red and blue colours indicate higher power in old and young participants, respectively, and non‐significant findings are labelled as n.s. (C) The same analysis as in (A) is performed on the MEG data, marking significant age‐related dynamic effects in state‐specific PSDs. (D) The same analysis as in (B) is performed on the MEG data, illustrating significant age‐related dynamic effects in state‐specific power maps. For both EEG and MEG, the data were source reconstructed using a standard MNI structural brain image.

FIGURE A7. 10‐state dynamic resting‐state networks of EEG and MEG. Dynamic RSNs of 96 subjects are shown for the (A) EEG and (B) MEG datasets, identical to Figure 5, but with 10 HMM states inferred from the data. Note that source reconstruction was performed using individual sMRI images for both EEG and MEG.

ACKNOWLEDGEMENTS

This research was supported by the National Institute for Health Research (NIHR) Oxford Health Biomedical Research Centre. The Wellcome Centre for Integrative Neuroimaging is supported by core funding from the Wellcome Trust (203139/Z/16/Z). M.E. and M.W. are supported by the Wellcome Trust (106183/Z/14/Z and 215573/Z/19/Z), the New Therapeutics in Alzheimer's Diseases (NTAD) study supported by UK MRC, and the Dementia Platform UK (RG94383/RG89702). M.W. is further supported by the NIHR Oxford Health Biomedical Research Centre (NIHR203316). C.G. is supported by the Wellcome Trust (215573/Z/19/Z). The views expressed are those of the author(s) and not necessarily those of the NIHR or the Department of Health and Social Care.

DATA AVAILABILITY STATEMENT

In this study, we used the openly available resting‐state EEG and MEG datasets called the Leipzig Study for Mind‐Body‐Emotion Interactions (LEMON) dataset (Babayan et al., 2019) and the Cambridge Centre for Ageing and Neuroscience (Cam‐CAN) dataset (Shafto et al., 2014), respectively.

1 These age ranges were selected for the young and old cohorts based on the age distribution of the LEMON dataset (i.e., how the two groups were divided in the original dataset) (Babayan et al., 2019).

2 It is worth considering that, in matching the rank of both the MEG and EEG data to 50, we may have lost some of the additional spatial resolution available in MEG, as the MEG data has a higher rank than EEG data at the sensor level, even after maxfiltering. Therefore, this matched rank could have artificially increased the similarity between the EEG and MEG descriptions beyond what would naturally occur. In this article, however, we decided that the most appropriate strategy is to align the preprocessing steps as closely as possible across both modalities, to ensure that any observed differences between the two are not due to trivial variations arising from preprocessing and source reconstruction steps.

3 This procedure applied exclusively to analyses involving data reconstructed with subject sMRI images. For data reconstructed with the standard MNI template, the procedure was not reiterated. Instead, the optimal run was chosen from a total of 10 runs.

4 Typically, state alignments in multiple HMM runs are conducted automatically through the linear sum assignment algorithm (Crouse, 2016). This approach tends to yield satisfactory results only when HMMs are trained on datasets of identical modalities. Aligning states across HMM runs trained on datasets of different modalities thus necessitated manual intervention due to the algorithm's decreased efficacy under such conditions. However, it should be noted that while manual alignment by eye was reliable for matching six states in our case, this process may still have introduced some human‐induced bias.
==== Refs
REFERENCES

Babayan, A. , Erbey, M. , Kumral, D. , Reinelt, J. D. , Reiter, A. M. F. , Röbbig, J. , Schaare, H. L. , Uhlig, M. , Anwander, A. , Bazin, P. L. , Horstmann, A. , Lampe, L. , Nikulin, V. V. , Okon‐Singer, H. , Preusser, S. , Pampel, A. , Rohr, C. S. , Sacher, J. , Thöne‐Otto, A. , … Villringer, A. (2019). A mind‐brain‐body dataset of MRI, EEG, cognition, emotion, and peripheral physiology in young and old adults. Scientific Data, 6 (1 ), 180308. 10.1038/sdata.2018.308 30747911
Baker, A. P. , Brookes, M. J. , Rezek, I. A. , Smith, S. M. , Behrens, T. , Probert Smith, P. J. , & Woolrich, M. (2014). Fast transient networks in spontaneous human brain activity. eLife, 3 , e01867. 10.7554/elife.01867 24668169
Barnes, J. J. , Woolrich, M. W. , Baker, K. , Colclough, G. L. , & Astle, D. E. (2015). Electrophysiological measures of resting state functional connectivity and their relationship with working memory capacity in childhood. Developmental Science, 19 (1 ), 19–31. 10.1111/desc.12297 25782537
Becker, R. , Vidaurre, D. , Quinn, A. J. , Abeysuriya, R. G. , Parker Jones, O. , Jbabdi, S. , & Woolrich, M. W. (2020). Transient spectral events in resting state MEG predict individual task responses. NeuroImage, 215 , 116818. 10.1016/j.neuroimage.2020.116818 32276062
Bishop, C. M. (2006). Pattern recognition and machine learning. Springer.
Biswal, B. , Zerrin Yetkin, F. , Haughton, V. M. , & Hyde, J. S. (1995). Functional connectivity in the motor cortex of resting human brain using echo‐planar MRI. Magnetic Resonance in Medicine, 34 (4 ), 537–541. 10.1002/mrm.1910340409 8524021
Brookes, M. J. , Leggett, J. , Rea, M. , Hill, R. M. , Holmes, N. , Boto, E. , & Bowtell, R. (2022). Magnetoencephalography with optically pumped magnetometers (OPM‐MEG): The next generation of functional neuroimaging. Trends in Neuroscience, 45 (8 ), 621–634. 10.1016/j.tins.2022.05.008
Brookes, M. J. , Woolrich, M. , Luckhoo, H. , Price, D. , Hale, J. R. , Stephenson, M. C., Barnes, G. R., Smith, S. M., & Morris, P. G. (2011). Investigating the electrophysiological basis of resting state networks using magnetoencephalography. National Academy of Sciences of the United States of America, 108 (40 ), 16783–16788. 10.1073/pnas.1112685108
Cichy, R. M. , & Pantazis, D. (2017). Multivariate pattern analysis of MEG and EEG: A comparison of representational structure in time and space. NeuroImage, 158 , 441–454. 10.1016/j.neuroimage.2017.07.023 28716718
Colclough, G. , Brookes, M. , Smith, S. , & Woolrich, M. (2015). A symmetric multivariate leakage correction for MEG connectomes. NeuroImage, 117 , 439–448. 10.1016/j.neuroimage.2015.03.071 25862259
Colclough, G. L. , Woolrich, M. W. , Tewarie, P. K. , Brookes, M. J. , Quinn, A. J. , & Smith, S. M. (2016). How reliable are MEG resting‐state connectivity metrics? NeuroImage, 138 , 284–293. 10.1016/j.neuroimage.2016.05.070 27262239
Coquelet, N. , De Tiège, X. , Destoky, F. , Roshchupkina, L. , Bourguignon, M. , Goldman, S. , Peigneux, P., & Wens, V. (2020). Comparing MEG and high‐density EEG for intrinsic functional connectivity mapping. NeuroImage, 210 , 116556. 10.1016/j.neuroimage.2020.116556 31972279
Crouse, D. F. (2016). On implementing 2D rectangular assignment algorithms. IEEE Transactions on Aerospace and Electronic Systems, 52 (4 ), 1679–1696. 10.1109/TAES.2016.140952
de Pasquale, F. , Della Penna, S. , Snyder, A. Z. , Lewis, C. , Mantini, D. , Marzetti, L. , Belardinelli, P. , Ciancetta, L. , Pizzella, V. , Romani, G. L. , & Corbetta, M. (2010). Temporal dynamics of spontaneous MEG activity in brain networks. National Academy of Sciences of the United States of America, 107 (13 ), 6040–6045. 10.1073/pnas.0913863107
Ferrari, M. , & Quaresima, V. (2012). A brief review on the history of human functional near‐infrared spectroscopy (fNIRS) development and fields of application. NeuroImage, 63 (2 ), 921–935. 10.1016/j.neuroimage.2012.03.049 22510258
Gohil, C. , Huang, R. , Roberts, E. , van Es, M. W. J. , Quinn, A. J. , Vidaurre, D. , & Woolrich, M. W. (2023). osl‐dynamics: A toolbox for modelling fast dynamic brain activity. eLife, 12 , RP91949. 10.7554/eLife.91949.2
Gohil, C. , Roberts, E. , Timms, R. , Skates, A. , Higgins, C. , Quinn, A. , Pervaiz, U. , van Amersfoort, J. , Notin, P. , Gal, Y. , Adaszewski, S. , & Woolrich, M. (2022). Mixtures of large‐scale dynamic functional brain network modes. NeuroImage, 263 , 119595. 10.1016/j.neuroimage.2022.119595 36041643
Hoerl, A. E. , & Kennard, R. W. (1970). Ridge regression: Biased estimation for nonorthogonal problems. Technometrics, 12 (1 ), 55. 10.2307/1267351
Holmes, A. P. , Blair, R. C. , Watson, J. D. G. , & Ford, I. (1996). Nonparametric analysis of statistic images from functional mapping experiments. Journal of Cerebral Blood Flow & Metabolism, 16 (1 ), 7–22. 10.1097/00004647-199601000-00002 8530558
Huang, R. , Gohil, C. , & Woolrich, M. (2024). Modelling variability in dynamic functional brain networks using embeddings. bioRxiv. 10.1101/2024.01.29.577718
Hutchison, R. M. , Womelsdorf, T. , Allen, E. A. , Bandettini, P. A. , Calhoun, V. D. , Corbetta, M. , Della Penna, S. , Duyn, J. H. , Glover, G. H. , Gonzalez‐Castillo, J. , Handwerker, D. A. , Keilholz, S. , Kiviniemi, V. , Leopold, D. A. , de Pasquale, F. , Sporns, O. , Walter, M. , & Chang, C. (2013). Dynamic functional connectivity: Promise, issues, and interpretations. NeuroImage, 80 , 360–378. 10.1016/j.neuroimage.2013.05.079 23707587
Hyvarinen, A. (1999). Fast and robust fixed‐point algorithms for independent component analysis. IEEE Transactions on Neural Networks, 10 , 626–634. 10.1109/72.761722 18252563
Jenkinson, M. , Pechaud, M. , & Smith, S. (2005). BET2: MR‐based estimation of brain, skull and scalp surfaces. In Eleventh annual meeting of the organization for human brain mapping, 17:167. Toronto, Canada.
Kingma, D. P. , & Ba, J. (2014). Adam: A method for stochastic optimization. arXiv. p. 1412.6980. 10.48550/arXiv.1412.6980
Knyazev, G. G. , Savostyanov, A. N. , Bocharov, A. V. , Tamozhnikov, S. S. , & Saprigyn, A. E. (2016). Task‐positive and task‐negative networks and their relation to depression: EEG beamformer analysis. Behavioural Brain Research, 306 , 160–169. 10.1016/j.bbr.2016.03.033 27001453
Leuchter, A. F. , Cook, I. A. , Hunter, A. M. , Cai, C. , & Horvath, S. (2012). Resting‐state quantitative electroencephalography reveals increased neurophysiologic connectivity in depression. PLoS One, 7 (2 ), e32508. 10.1371/journal.pone.0032508 22384265
Lin, J. (1991). Divergence measures based on the Shannon entropy. IEEE Transactions on Information Theory, 37 (1 ), 145–151. 10.1109/18.61115
Liu, Q. , Farahibozorg, S. , Porcaro, C. , Wenderoth, N. , & Mantini, D. (2017). Detecting large‐scale networks in the human brain using high‐density electroencephalography: Imaging brain networks with high density EEG. Human Brain Mapping, 38 (9 ), 4631–4643. 10.1002/hbm.23688 28631281
Lottman, K. K. , Gawne, T. J. , Kraguljac, N. V. , Killen, J. F. , Reid, M. A. , & Lahti, A. C. (2019). Examining resting‐state functional connectivity in first‐episode schizophrenia with 7T fMRI and MEG. NeuroImage: Clinical, 24 , 101959. 10.1016/j.nicl.2019.101959 31377556
Mahmutoglu, M. A. , Rupp, A. , & Baumgärtner, U. (2022). Simultaneous EEG/MEG yields complementary information of nociceptive evoked responses. Clinical Neurophysiology, 143 , 21–35. 10.1016/j.clinph.2022.08.005 36087398
Mantini, D. , Perrucci, M. G. , Del Gratta, C. , Romani, G. L. , & Corbetta, M. (2007). Electrophysiological signatures of resting state networks in the human brain. National Academy of Sciences of the United States of America, 104 (32 ), 13170–13175. 10.1073/pnas.0700668104
McFadyen, J. , & Dolan, R. J. (2023). Spatiotemporal precision of neuroimaging in psychiatry. Biological Psychiatry, 93 , 671–680. 10.1016/j.biopsych.2022.08.016 36376110
Michel, C. M. , & Koenig, T. (2018). EEG microstates as a tool for studying the temporal dynamics of whole‐brain neuronal networks: A review. NeuroImage, 180 , 577–593. 10.1016/j.neuroimage.2017.11.062 29196270
Ng, A. Y. (2004). Feature selection, L1 vs. L2 regularization, and rotational invariance. In Proceedings of the 21st International Conference on Machine Learning. Banff, Canada.
Nichols, T. E. , & Holmes, A. P. (2002). Nonparametric permutation tests for functional neuroimaging: A primer with examples. Human Brain Mapping, 15 , 1–25. 10.1002/hbm.1058 11747097
Palva, J. M. , Wang, S. H. , Palva, S. , Zhigalov, A. , Monto, S. , Brookes, M. J. , Schoffelen, J. M. , & Jerbi, K. (2018). Ghost interactions in MEG/EEG source space: A note of caution on inter‐areal coupling measures. NeuroImage, 173 , 632–643. 10.1016/j.neuroimage.2018.02.032 29477441
Perrin, F. , Pernier, J. , Bertrand, O. , & Echallier, J. (1989). Spherical splines for scalp potential and current density mapping. Electroencephalography and Clinical Neurophysiology, 72 , 184–187. 10.1016/0013-4694(89)90180-6 2464490
Quinn, A. J. , Atkinson, L. Z. , Gohil, C. , Kohl, O. , Pitt, J. , Zich, C. , Nobre, A. C. , & Woolrich, M. W. (2024). The GLM‐spectrum: A multilevel framework for spectrum analysis with covariate and confound modelling. Imaging Neuroscience, 2 , 1–26. 10.1162/imag_a_00082
Quinn, A. J. , van Es, M. W. , Gohil, C. , & Woolrich, M. W. (2023). OHBA software library in python (OSL) (0.1.1). Zenodo. 10.5281/zenodo.6875059
Quinn, A. J. , Vidaurre, D. , Abeysuriya, R. , Becker, R. , Nobre, A. C. , & Woolrich, M. W. (2018). Task‐evoked dynamic network analysis through hidden Markov modeling. Frontiers in Neuroscience, 12 , 603. 10.3389/fnins.2018.00603 30210284
Rezek, I. , & Roberts, S. (2005). Ensemble hidden Markov models with extended observation densities for biosignal analysis. In D. Husmeier , R. Dybowski , & S. Roberts (Eds.), Probabilistic modeling in bioinformatics and medical informatics (pp. 419–450). Springer.
Robert, P. , & Escoufier, Y. (1976). A unifying tool for linear multivariate statistical methods: The RV‐coefficient. Applied Statistics, 25 (3 ), 257. 10.2307/2347233
Rosner, B. (1983). Percentage points for a generalized ESD many‐outlier procedure. Technometrics, 25 , 165–172. 10.2307/1268549
Shafto, M. A. , Tyler, L. K. , Dixon, M. , Taylor, J. R. , Rowe, J. B. , Cusack, R. , Calder, A. J., Marslen‐Wilson, W. D., Duncan, J., Dalgleish, T., Henson, R. N., Brayne, C., Matthews, F. E., & Cam‐CAN. (2014). The Cambridge Centre for Ageing and Neuroscience (Cam‐CAN) study protocol: A cross‐ sectional, lifespan, multidisciplinary examination of healthy cognitive ageing. BMC Neurology, 14 , 204. 10.1186/s12883-014-0204-1 25412575
Shulman, G. L. , Fiez, J. A. , Corbetta, M. , Buckner, R. L. , Miezin, F. M. , Raichle, M. E. , & Petersen, S. E. (1997). Common blood flow changes across visual tasks: II. Decreases in cerebral cortex. Journal of Cognitive Neuroscience, 9 (5 ), 648–663. 10.1162/jocn.1997.9.5.648 23965122
Siems, M. , Pape, A. , Hipp, J. F. , & Siegel, M. (2016). Measuring the cortical correlation structure of spontaneous oscillatory activity with EEG and MEG. NeuroImage, 129 , 345–355. 10.1016/j.neuroimage.2016.01.055 26827813
Sitnikova, T. A. , Hughes, J. W. , Ahlfors, S. P. , Woolrich, M. W. , & Salat, D. H. (2018). Short timescale abnormalities in the states of spontaneous synchrony in the functional neural networks in Alzheimer's disease. NeuroImage: Clinical, 20 , 128–152. 10.1016/j.nicl.2018.05.028 30094163
Smith, S. M. (2002). Fast robust automated brain extraction. Human Brain Mapping, 17 (3 ), 143–155. 10.1002/hbm.10062 12391568
Sockeel, S. , Schwartz, D. , Pélégrini‐Issac, M. , & Benali, H. (2016). Large‐scale functional networks identified from resting‐state EEG using spatial ICA. PLoS One, 11 (1 ), e0146845. 10.1371/journal.pone.0146845 26785116
Veen, B. V. , & Buckley, K. (1988). Beamforming: A versatile approach to spatial filtering. IEEE ASSP Magazine, 5 (2 ), 4–24. 10.1109/53.665
Vidaurre, D. , Abeysuriya, R. , Becker, R. , Quinn, A. J. , Alfaro‐Almagro, F. , Smith, S. M. , & Woolrich, M. W. (2018). Discovering dynamic brain networks from big data in rest and task. NeuroImage, 180 , 646–656. 10.1016/j.neuroimage.2017.06.077 28669905
Vidaurre, D. , Hunt, L. T. , Quinn, A. J. , Hunt, B. A. E. , Brookes, M. J. , Nobre, A. C. , & Woolrich, M. W. (2018). Spontaneous cortical activity transiently organises into frequency specific phase‐ coupling networks. Nature Communications, 9 (1 ), 2987. 10.1038/s41467-018-05316-z
Vidaurre, D. , Llera, A. , Smith, S. M. , & Woolrich, M. W. (2021). Behavioural relevance of spontaneous, transient brain network interactions in fMRI. NeuroImage, 229 , 117713. 10.1016/j.neuroimage.2020.117713 33421594
Vidaurre, D. , Quinn, A. J. , Baker, A. P. , Dupret, D. , Tejero‐Cantero, A. , & Woolrich, M. W. (2016). Spectrally resolved fast transient brain states in electrophysiological data. NeuroImage, 126 , 81–95. 10.1016/j.neuroimage.2015.11.047 26631815
Welch, P. (1967). The use of fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms. IEEE Transactions on Audio and Electroacoustics, 15 (2 ), 70–73. 10.1109/TAU.1967.1161901
