Spatial (mis)match between EEG and fMRI signal patterns revealed by spatio-spectral source-space EEG decomposition

This study aimed to directly compare electroencephalography (EEG) whole-brain patterns of neural dynamics with concurrently measured fMRI BOLD data. To achieve this, we aim to derive EEG patterns based on a spatio-spectral decomposition of band-limited EEG power in the source-reconstructed space. In a large dataset of 72 subjects undergoing resting-state hdEEG-fMRI, we demonstrated that the proposed approach is reliable in terms of both the extracted patterns as well as their spatial BOLD signatures. The five most robust EEG spatio-spectral patterns not only include the well-known occipital alpha power dynamics, ensuring consistency with established findings, but also reveal additional patterns, uncovering new insights into brain activity. We report and interpret the most reproducible source-space EEG-fMRI patterns, along with the corresponding EEG electrode-space patterns, which are better known from the literature. The EEG spatio-spectral patterns show weak, yet statistically significant spatial similarity to their functional magnetic resonance imaging (fMRI) blood oxygenation level-dependent (BOLD) signatures, particularly in the patterns that exhibit stronger temporal synchronization with BOLD. However, we did not observe a statistically significant relationship between the EEG spatio-spectral patterns and the classical fMRI BOLD resting-state networks (as identified through independent component analysis), tested as the similarity between their temporal synchronization and spatial overlap. This provides evidence that both EEG (frequency-specific) power and the BOLD signal capture reproducible spatio-temporal patterns of neural dynamics. Instead of being mutually redundant, these only partially overlap, providing largely complementary information regarding the underlying low-frequency dynamics.

1 Introduction

Finding a reliable link between electroencephalography (EEG) and functional magnetic resonance imaging (fMRI) data measured during resting brain activity seems to be just as difficult as fulfilling the “to rest” instruction for the volunteering subject. Yet, there is a variety of strong reasons for attempting to fuse the two modalities. First, finding a relation may help in understanding the origin and character of each type of signal. Second, by using simultaneous EEG-fMRI measurements, we can mitigate the limitations of each modality when applied separately, allowing for a more accurate study of both healthy brains and brain disorders. The third motivation is to search for aspects of the EEG signal that could reproduce the well-documented features of the BOLD signal - the resting-state network (RSN) dynamics. Discovering such relationships could potentially enable the replacement of costly fMRI measurements with more affordable EEG experiments in a greater number of scenarios in the future. This could ultimately provide more affordable diagnostics and treating methods (Vega et al., 2022), increase the ecological validity of experiments as well as enable a wider range of experiment designs.

Since Logothetis et al. (2001) experimentally showed a delayed hemodynamic response of the Blood-oxygen-level-dependent (BOLD) fMRI signal to the measured local field potential (LFP) signal envelope, questions concerning the details of the relation of electrophysiological signal and fMRI have been around and are still far from satisfactorily answered, particularly when one moves outside of the well-controlled situation of a controlled task experiment and localized intracortical EEG recording. The standard non-invasive EEG signal can be naturally described in terms of three domains: temporal, spatial, and frequency domains. Initial resting-state studies aimed to find the BOLD correlate of the most prominent resting-state EEG feature, that is, the occipital alpha activity (Berger, 1929). In those studies, the spatial domain mainly was reduced by selecting the occipital subset of electrodes. Similarly, the frequency domain was limited to the alpha frequency band by utilizing a band-limited power (BLP) regressor for BOLD signal prediction (Goldman et al., 2002; Moosmann et al., 2003; Feige et al., 2005). Other studies varied in terms of the used electrode subsets (Laufs et al., 2003a) or frequency bands (Laufs et al., 2003b), and also the potential use of multiple frequency specific regressors in the same statistical general linear model (GLM) (Tyvaert et al., 2008; de Munck et al., 2009). Despite a very similar experimental design, the initial studies generally identified two distinct spatial BOLD fMRI correlation patterns associated with the alpha regressor. Namely, a widespread bilateral frontoparietal correlation pattern was reported by Laufs et al. (2003a), while (Goldman et al., 2002; Moosmann et al., 2003) observed an occipitoparietal pattern. Thus, while both patterns showed a negative relationship between the EEG alpha BLP regressor and the BOLD signal (which was in line with the intuitive interpretation of alpha as the idling rhythm corresponding to low metabolism requirements), the localization of these anticorrelations differed substantially between the studies, creating confusion concerning which brain areas were actually most involved in this idling. Later, by reanalyzing their data, Laufs et al. (2006) showed that both these alpha correlates can be observed within the same dataset when considering different subsets of subjects. Similarly, Gonçalves et al. (2006) showed that the alpha correlate shows a high level of inter- as well as intrasubject variability.

Other studies moved the focus of the field to other frequency bands than alpha: Scheeringa et al. (2012) showed that frontal theta rhythm power fluctuation positively correlates with hubs of the default mode network (DMN). Further, Mantini et al. (2007) found that the BLP time series within standard EEG frequency bands, calculated as an average signal across all electrodes, correlates with the time series of multiple BOLD independent components (ICs), specifically resting-state networks (RSNs). This, along with previous findings, suggests that the dynamics of RSNs may be reflected in EEG as a combination of multiple frequency and spatial patterns. The first EEG-fMRI data-driven method that did not impose any a priori assumptions on the spatial, temporal, or frequency domains was introduced in a study by Miwakeichi et al. (2004) and later in a simultaneous EEG-fMRI study by Martınez-Montes et al. (2004), where the three-dimensional (3D) spatio-temporal-frequency array was factorized using parallel factor analysis (PARAFAC) (Harshman and Lundy, 1994) to obtain spatio-temporal-frequency components. In line with the study by Goldman et al. (2002), the authors discovered a single alpha band occipital component that exhibited a statistically significant temporal correspondence with the BOLD signal in the occipital lobe. The authors also reported a frontal theta component in accordance with the study by Scheeringa et al. (2012), which did not show a statistically significant correlation with the BOLD signal. A similar analysis was conducted by Mareček et al. (2017); Marecek et al. (2016), where several components representing brain rhythms and artifacts were identified. Although these trilinear decomposition methods represent a truly data-driven approach for identifying EEG spatio-temporal-frequency patterns, they may struggle to identify spatially and frequency-wise distinct BLP fluctuations. This is because each component in the model has only one spatial signature across all frequency bands and vice versa.

Therefore, it appears more convenient to use a model that permits different frequency bands (or their combinations) to exhibit distinct spatial signatures within a single component. One such method, called spatio-spectral decomposition (Bridwell et al., 2013), is based on a concatenation of spatial and frequency domains. A matrix decomposition, such as independent component analysis (ICA) (Hyvärinen and Oja, 2000), is then applied to this joint spatio-spectral dimension. Bridwell et al. (2013) introduced this idea for simultaneous EEG-fMRI data and associated the resulting components with the BOLD RSNs. Bridwell et al. (2018) further validated several blind source separation algorithms (BSS) using both realistic and simulated data. Labounek et al. (2018) reported robust spatio-spectral patterns (SSPs) across three datasets from different paradigms, and later (Labounek et al., 2019) established a temporal relationship between SSPs and the BOLD signal in various hubs of BOLD resting state networks (RSNs). Additionally, a statistically significant relationship between several SSPs and the stimulus types across different paradigms was also found by Labounek et al. (2021). A potential drawback of this method is that the spatial domain is defined by the electrode-space, making the spatial signatures difficult to interpret or directly compare with fMRI BOLD maps. Readers will notice that we review only studies within the so-called EEG-informed fMRI integration approach (Abreu et al., 2018) as well as studies that primarily derive BLP features from EEG data to predict the BOLD signal. For a broader overview of integration approaches and EEG feature extraction methods, readers are encouraged to refer to the comprehensive review articles by Jorge et al. (2014); Abreu et al. (2018); Murta et al. (2015).

Recently, there has been an effort to make high-density EEG (hdEEG) a true neuroimaging tool, especially when compared to standard low-density EEG, where the few recording electrodes hardly provide an image of brain activity in the computer science sense; the methods applied correspond much more to its time series character. In a recent resting-state study, Liu et al. (2017) obtained a comprehensive set of RSNs using precisely source-localized BLP ICA decomposition, similar to what is typically performed on the BOLD signal. In a study by Liu et al. (2018), the authors highlighted the significance of accurately conducted source localization steps, such as precise electrode localization (Marino et al., 2016), a multilayer individual head model to establish a dependable forward model (Taberna et al., 2021), and the density of electrode coverage. In a study by Marino et al. (2019), the authors recently obtained EEG source-reconstructed DMN that was spatially and temporally related to the BOLD-derived DMN from simultaneous resting-state EEG-fMRI recordings. In contrast, a different analysis aimed at obtaining DMN characteristics from EEG was presented in a study by Prestel et al. (2018), where the authors highlight that the component with a high temporal correspondence to the DMN was associated with eye movement artifacts.

Besides the above-mentioned studies, there are other recent studies investigating the link between resting-state EEG and fMRI in the EEG source-reconstructed space. In a study by Sockeel et al. (2016), the authors concatenated temporal and frequency domains to perform a matrix decomposition to obtain separate temporal dynamics for each of the five frequency bands. They reported not only a correspondence between EEG and fMRI RSNs but also substantial mismatches. A similar analysis was also performed by Li et al. (2018). In a study by Abreu et al. (2020), despite having a dataset with only 10 subjects, the authors report a high spatial correspondence between the EEG and fMRI RSNs and also a statistically significant match between the EEG and fMRI dynamical functional connectivity (dFC) states. Another study by Yuan et al. (2016) aimed to obtain the full set of EEG RSNs and compare them spatially and temporally with BOLD RSNs. They found a spatial correspondence but, at the same time, very low temporal correspondence between EEG and fMRI RSNs. In a study by Meyer et al. (2013), the authors claim that the correlation between EEG BLPs and BOLD RSNs is unstable over time.

The goal of this study is to integrate the literature on spatio-spectral decomposition with reliable source localization techniques and to introduce a spatio-spectral decomposition in the source-reconstructed space that leverages the data-driven nature of spatio-spectral decomposition and cortical sources in the spatial domain. This approach allows us to directly examine the relationship between these modalities in a way that was previously impossible. In a large simultaneous EEG-fMRI dataset, we explore the stability of decomposition in both EEG source- and electrode-spaces, as well as how the EEG patterns correspond to BOLD activation patterns and BOLD RSNs. These are all questions currently at the center of hot debate in the literature.

Throughout this manuscript, we emphasize using robust evaluation methods that offer an unbiased perspective on intermodal relationships. Specifically, we start by examining the stability of EEG decomposition itself in Sections 2.6, 3.1.1 to ensure the algorithmic stability of the EEG decomposition. Next, we provide a biased perspective on EEG feature reproducibility testing in the Section 2.9.1 and advocate for more advanced evaluation methods that take into account not only EEG patterns but also their relationship with fMRI BOLD data. To address this, we introduce EEG-fMRI integration reproducibility analysis (Sections 2.9.2, 3.1.2), which assesses the reproducibility of derived EEG patterns and their association with fMRI data. Next, we identify the most reproducible EEG-fMRI patterns (Sections 2.10, 3.2). Finally, employing the source localization methods, we directly test three hypotheses regarding spatio-temporal relationships between EEG and BOLD data: (1) Are EEG patterns co-localized with their corresponding fMRI BOLD activation maps? (Methods Section ??); (2) Do EEG patterns that explain a greater portion of BOLD data variability exhibit higher colocalization? (Methods Section ??); and (3) Are EEG patterns spatio-temporally associated with BOLD resting-state networks? (Methods Section 2.11.3). The results of these hypothesis tests are presented in Section 3.3. We believe that addressing these questions is essential for developing a more reliable understanding of the EEG/fMRI relationship.

2 Methods

2.1 Participants and experimental design

For this study, we used two datasets. The primary dataset includes simultaneous EEG-fMRI recordings along with various types of structural images essential for creating individual models of the brain. The second dataset comprises out-of-scanner EEG recordings without functional MRI data, and it illustrates its robustness concerning changes in the dataset, preprocessing, and the absence of MRI artifacts. Both are described in the following subsections.

2.1.1 EEG-fMRI dataset

We used a dataset from a single study obtained over a 4-year period at the Central European Institute of Technology (CEITEC) in Brno, Czech Republic. We analyzed 72 healthy participants (mean age: 31.4, range: 18.3–50.6; 36 males and 36 females) with consistent data acquisition parameters. All participants underwent a 20-min eyes-closed resting-state EEG-fMRI recording session with instructions to lie still, avoid falling asleep, and not think about anything in particular. The study was approved by the ethics committee of Masaryk University and was conducted in accordance with the Declaration of Helsinki. All participants provided written informed consent to participate in the study.

2.1.2 Out-of-scanner EEG dataset

To relate the EEG patterns extracted from the combined EEG-fMRI dataset to the patterns obtainable when applying common processing steps to out-of-scanner EEG, we utilized the following out-of-scanner EEG dataset. The EEG data from healthy subjects were acquired at the National Institute of Mental Health (NIMH, Klecany, Czech Republic) as part of the multimodal prospective database of the first episodes of psychotic illness project. The dataset consists of 50 healthy participants (mean age: 29.9, range: 19.1–43.0; 21 males and 29 females). Each participant underwent a 15-min resting-state paradigm (first 5 min with eyes open and the last 10 min with eyes closed). The 10 min of resting-state recording with eyes closed were preprocessed and analyzed. The study received approval from the local ethics committee of the NIMH and was carried out in accordance with the Declaration of Helsinki. All participants provided written informed consent to participate in the study. The differences in acquisition parameters and processing steps between this out-of-scanner dataset and the previously mentioned EEG-fMRI dataset are briefly discussed at the end of each of the following subsections. In both research and clinical practice, the same acquisition and processing steps typically cannot be guaranteed. Therefore, some differences should be expected in most subsequent studies that aim to apply the proposed methodology for EEG SSP extraction to datasets from other EEG experiments. Thus, this out-of-scanner dataset is not meant for testing replicability under optimal conditions; instead, it serves as an illustrative example that clarifies the degree to which the results observed in the EEG-fMRI dataset are robust when applied to another dataset with established standards, albeit with slightly different processing.

2.2 Data acquisition

The magnetic resonance imaging was performed using a 3T Siemens Prisma magnetic resonance scanner equipped with a 64-channel radiofrequency (RF) receiving head coil (Siemens Healthineers, Erlangen, Germany). The functional magnetic resonance imaging was performed with a multiband multiecho two-dimensional (2D) echo-planar imaging (EPI) sequence with the following parameter settings: Multiband factor: 6; number of echos: 3; PAT factor: 2 (PAT) factor: 2; repetition time (TR) = 650 ms; echo time (TE) = 14.60/33.56/52.52 ms; 48 axial slices with 3 mm slice thickness; slice: 64 × 64 matrix, 194 × 194 mm; number of volumes: 1,840; and flip angle (FA) = 30°. Two structural MRI images were acquired. The first one was obtained without an electrode net with the T1 Magnetization Prepared – RApid Gradient Echo (MPRAGE) sequence and the following parameters setting: TR = 2,300 ms; TE = 2.34 ms; inversion time (TI) = 900 ms; 240 sagittal slices with 1 mm slice thickness; slice: 260 × 256 matrix, 260 × 256 mm; FA = 8°; PAT factor: 2. The second structural MRI image with T1 MPRAGE sequence was acquired with the electrode net put on for precise electrode localization and the parameter setting: TR = 2,300 ms; TE = 2.33 ms; TI = 900 ms; 240 sagittal slices with 1 mm slice thickness; slice: 224 × 224 matrix, 224 × 224 mm; FA = 8°; and PAT factor: 7. The EEG data with a sampling frequency of 1,000 Hz were recorded with an EGI Hydrocell MR-compatible 256-channel high-density electrode net plugged into the EGI GES 400 signal amplifier (Electrical Geodesics, Inc., Eugene, Oregon, USA). To obtain the ECG signal, one additional channel was recorded. A breathing belt was also attached to the participant's chest to record breathing cycles.

The out-of-scanner EEG data acquisition was performed with the very same recording set up with the only difference of having the non-magnetic resonance (MR)-compatible version of the 256-channel high-density electrode net. The individual structural T1 MRI MPRAGE image was also provided with the following acquisition parameters: Voxel size of 1 × 1 × 1 mm, 224 sagittal slices, TE = 4.63 ms, TR = 2,300 ms, TI = 900 ms, FA = 10°, and time of acquisition (TA) = 5:30, and field of view (FOV) = 256 mm.

2.3 EEG data preprocessing

The raw EEG data were preprocessed by a fully automated pipeline introduced in Marino et al. (2019); Liu et al. (2017, 2018). The preprocessing pipeline utilizes built-in and in-house MATLAB (MathWorks, Natick, MA, USA) functions as well as SPM (Penny et al., 2011), Fieldtrip (Oostenveld et al., 2011), and EEGLAB (Delorme and Makeig, 2004) toolboxes. The preprocessing steps are summarised in the following paragraph.

The first step involved gradient artifact removal using the FMRI Artifact Slice Template Removal (FASTR) method (Niazy et al., 2005) in EEGLAB, followed by the removal of ballistocardiogram artifacts through the adaptive optimal basis set method introduced in Marino et al. (2018). The channels with poor signal quality were identified based on low correlation with all other channels in across a frequency band (1—80 Hz) and the variance in EEG non-physiological frequency band 200—250 Hz. The latter metric serves as a noise variance indicator. In case when at least one criterion marked an outlier in distribution across channels, the channel time course was subsequently interpolated by the time courses of neighboring channels based on the weighted average (electrode distances) scheme implemented in Fieldtrip (Oostenveld et al., 2011). Subsequently, EEG data were filtered in the 1—80 Hz frequency band. ICA was applied to remove movement and other biological artifacts including electrooculographic (EOG) and electromyographic (EMG) artifacts from the EEG recordings. For that purpose, FastICA (Hyvärinen and Oja, 2000) algorithm based on a deflation approach and hyperbolic tangent as the contrast function was applied and artifactual components were detected based on three parameters, namely the correlation values between ICs time course and reference EOG and EMG signals, the similarity of ICs power spectrum with a 1/f function, and kurtosis of ICs timecourse (Mantini et al., 2008). The artifact-suppressed EEG data were subsequently filtered into several frequency bands: delta (δ, 1—4 Hz), theta (θ, 4—8 Hz), alpha (α, 8—12 Hz), low beta (β₁, 12—15 Hz), middle beta (β₂, 15—18 Hz), high beta (β₃, 18—30 Hz), and gamma (γ, 30—44 Hz). Although the EEG data were comprehensively preprocessed, we decided to exclude all cheek electrodes and the two lowest layers of the neck electrodes from further analysis since there is a concern in current literature that sensors placed at those electrode sites contain disproportionally more artifacts (Vorderwülbecke et al., 2020). Therefore, for all further analyses we utilized 195 out of 257 electrodes. As the last step of the preprocessing, we rereferenced the EEG data to average reference.

The out-of-scanner EEG data were preprocessed manually by the expert with the following steps: EEG data were preprocessed in BESA version 7.0 software (MEGIS, Munich, Germany) by removing noisy epochs followed by a semi-automated ECG and eye movement-related artifact correction by the signal-space projection method (Uusitalo and Ilmoniemi, 1997).

2.4 MRI data preprocessing

The functional MRI data were preprocessed in MATLAB utilizing a combination of SPM12 (Penny et al., 2011) pipelines and in-house scripts. At first, the fMRI volumes were spatially realigned, followed by a fusion of three echos by a weighted averaging based on a temporal signal-noise-ratio (tSNR). To regress out cardiac and breathing artifactual signals from the BOLD signal, the RETROICOR method (Glover et al., 2000) informed by the ECG and breathing signals was employed. Then, a coregistration of the structural MRI was followed to the average functional image, and a structural MRI was normalized to MNI space based on an image segmentation of gray and white matter. Functional data were subsequently transformed to MNI space using the combined transformation matrices from the previous step. All volumes were also resampled into 3 × 3 × 3 mm isotropic voxels. Individual fMRI volumes were then spatially smoothed by a Gaussian filter (full width at half maximum [FWHM] 5 mm). To further mitigate the risk of spurious correlations that can cause common artifacts appearing over the entire volume and having a non-physiological nature, the signal from the gray matter was orthogonalized to the following proxies of contributing artifactual signals: a bank of sinusoidal signals with low frequencies (periods slower than 128 s), the first PCA component from the voxel signals of white matter mask, the first PCA component from the voxel signals of CSF mask, set of 24 rotation and translation parameters of estimated head motion, namely the estimated motion parameters themselves, their first differences, squares, and squared first differences. Finally, we applied a low-pass filter to the BOLD data with a cut-off frequency 0.09 Hz, which is a typical filtering step for resting-state connectivity analyses. The resulting preprocessed data were utilized for subsequent analyses.

2.5 EEG source localization

To estimate the sources of brain activity as reliably as possible, we implemented a source localization pipeline using individual-level data. Each step is briefly described in the following paragraphs.

Head tissue segmentation of structural T1 images without the electrode net was performed by the automated 12-compartment segmentation tool (Taberna et al., 2021), which consists of image preprocessing, tissue probability mapping, and tissue segmentation steps. Subsequently, we followed a source localization pipeline based on the Fieldtrip toolbox (Oostenveld et al., 2011) functions. Based on all 12 compartments of the segmented T1, the hexahedral mesh was generated, and compartment conductivities were assigned based on Liu et al. (2018) to define a head model. A precise electrode localization for the head model, provided by the Multimodal and Functional Imaging Laboratory of the Central European Institute of Technology (Brno, Czech Republic), was obtained semi-automatically from the acquired T1 image with the electrode net on, see Section 2.2. Electrode locations were manually marked on rendered head surface where the electrode artifact was prominent. After the electrode location definition, structural T1 image with electrode net was coregistered to the second structural T1 without electrode net and resulting transformation matrix was also applied to the electrode positions. Electrodes were projected to the closest surface point of the head model. The source model was generated in the brain gray matter and cerebellar gray matter compartments with a grid size of 6 mm. The finite element method (FEM) SimBio (Vorwerk et al., 2018) solver was used to compute the leadfield matrix. The sources were estimated by the eLORETA inverse algorithm (Pascual-Marqui et al., 2011). Furthermore, we established a template source model utilizing the MNI-template anatomy from SPM (Penny et al., 2011) and the methodology stated above. This template source model was subsequently designated as the reference for all subsequent analyses. Results obtained from the source space of the EEG-fMRI dataset, out-of-scanner EEG dataset, and fMRI were interpolated to this reference model using the nearest neighbor method. The source localization pipeline for out-of-scanner EEG data generally followed the same structure with the following differences: The source model resolution was 10 mm instead of 6 mm, the inverse warped source model grid from the MNI to the individual space was used instead of a regular grid directly generated in the individual space, the standard 5-layer head model in the Fieldtrip toolbox (Oostenveld et al., 2011) was used instead of the advanced 12-layer model in the MRTIM toolbox (Taberna et al., 2021), an electrode template was co-registered with an individual head structure and projected onto a head surface instead of individual electrode positioning, as discussed in Section 2.1.2. Here, we followed a well-established Fieldtrip toolbox (Oostenveld et al., 2011) pipeline, mimicking a typical standard processing setup.

2.6 Electrode-space spatio-spectral decomposition

In the present study, we assume that the electrical activity expressed as the EEG signal envelope measured at the electrode locations on the scalp is a linear mixture of source signals representing brain activity. We further assume that the brain electrical activity power may have a different spatial profile across different EEG frequency bands, which is a substantial difference between the spatio-spectral models (allowing thus estimation of a much richer structure) and the trilinear models (Marecek et al., 2016; Mareček et al., 2017) where the spatial mode is assumed to be the same across the whole frequency range, only with a different strength. The schematic flowchart of the spatio-spectral decomposition methodology is depicted in Figure 1. At first, we computed a signal envelope for each electrode and frequency band EEG time series by applying the Hilbert transform. Then we downsampled the data to the fMRI BOLD TR (650 ms) (Brookes et al., 2011), obtaining thus a time series of band-limited power (BLP). Subsequently, we implemented an outlier correction procedure using Tukey's fences (k = 3) on each time series. Outliers were identified and replaced with the mean value. Since the data were acquired in resting-state condition, we further band-pass filtered all BLP time series in a standard way to 0.008–0.09 Hz as for the resting-state fMRI BOLD data. Next, we concatenated the (space × time) BLP matrices of all frequency bands along the spatial domain, forming thus a joint spatio-spectral domain on the single-subject level. Each single-subject BLP time series was then z-score normalized to handle inter-individual variability. All single-subject normalized BLP matrices were subsequently concatenated along the temporal domain, forming a group BLP matrix X. Assuming that X is a linear mixture (by a mixing matrix A) of temporally independent BLP sources S, we can express the linear data generation process in a matrix form as

$\begin{array}{ll}X=AS,& (1)\end{array}$

where the matrix X has size (B ·N_e) x (P ·T). Here, B is the number of frequency bands, that is, 7, N_e is the number of electrodes, that is, 195, P stands for the number of subjects, that is, 72, and finally T is the number of single-subject BLP time points downsampled to TR, that is, 1,800. We applied a temporal (group-level) ICA by RUNICA algorithm (Makeig et al., 1995) with a PCA dimensionality reduction set to C = 30 components. Several information criteria, including both the Bayesian and the Akaike information criteria, were evaluated to ascertain the optimal number of components. Due to significant discrepancies between the criteria, the final number was chosen to align with the typical number of components in BOLD resting-state ICA analyses, see Section 4.1 for further discussion. Thus, we obtain a matrix S with dimensions C x (P ·T). Each row represents a single component time series considered a temporal signature of the given component. Matrix A, called mixing matrix with dimensions (B ·N_e) x C represents SSPs of the components (one in each column). Since ICA decomposition is based on iterative algorithms with random initialization, we utilized the ICASSO tool (Himberg et al., 2004) to investigate the algorithmic stability of the ICA decomposition itself by performing the ICA decomposition 20 times with different initial conditions. The cluster centroid time series were then considered the most reliable estimate of components and, Therefore, utilized for the subsequent analyses. Furthermore, SSPs were not obtained directly from the matrix A but via correlating each component time series (ICASSO cluster centroids) with each time series in BLP matrix X. The correlation coefficients were subsequently transformed by a Fisher Z-transformation. Obviously, each of the C = 30 SSPs consists of B spatial patterns (one for each frequency band), and each of C = 30 temporal signatures in S consists of the concatenation of P individual time series (one for each subject). We also performed the spatio-spectral decomposition on the first and second half of the dataset for reproducibility analyses and statistical evaluation. To handle temporal discontinuities in the case of the out-of-scanner dataset, we treated each data epoch separately up to a point before concatenation in the time domain. After that, the epochs within the subject were concatenated and subsequently across subjects. Finally, the ICA decomposition was performed in a similar way as for the previous dataset.

Figure 1

Figure 1. A schematic flowchart of the spatio-spectral decomposition in both electrode-space and source space with highlighted main steps and examples of SSPs. Please note that electrode- and source-space data are subjected to the ICA decomposition separately. The label electrodes/dipoles is meant to denote in what dimension both approaches differ. Also note that the actual size of matrices in Dipoles dimension (N_d, source-space version) is much higher compared to Electrodes dimension (N_e, electrode-space version).

2.7 Source–space spatio-spectral decomposition

The source–space spatio-spectral decomposition generally follows the workflow visualized in Figure 1 and described in Section 2.6 with the following differences specific to the source-projected EEG data. At each source model position, the dipole moment is expressed as three time series, one for each direction of the Cartesian coordinate system. Under the assumption that the direction of the dipole moment is not fixed but may rotate freely in time, we estimate the (scalar) signal power (amplitude) time series as in Liu et al. (2017)

$\begin{array}{ll}s(t)=\sqrt{j_x^2(t)+j_y^2(t)+j_z^2(t)},& (2)\end{array}$

where s(t) is a BLP time series and j_x, j_y, j_z are dipole moments in each of x, y, and z axes of the Cartesian coordinate system. Instead of N_e spatial points or electrodes, in the source–space case we obtain N_d spatial points. These points correspond to the individual source model positions. As mentioned in Section 2.5, individual BLP time series were transformed to the template source model positions to allow concatenation across subjects. Decomposing the source BLP matrix enables us to compare the SSPs spatially with the statistical maps of the explained BOLD signal or with the BOLD-derived RSN spatial signatures. This was carried out by also interpolating the BOLD voxel positions to the template source model positions (rather than interpolating source model positions to the positions of the BOLD voxels), which was advantageous in terms of the computational complexity of the subsequent analyses. The out-of-scanner SSPs were also transformed into the template source model positions.

2.8 BOLD signatures of EEG spatio-spectral patterns

We reconstructed the EEG spatio-spectral component time series S_p of individual subjects by splitting S to P segments of length T. To assess the single-subject relationship between EEG spatio-spectral component time courses and the BOLD signal in each voxel, the general linear model (GLM) was utilized. Note that several resting-state studies showed a substantial HRF variability of band-limited power regressor both when derived from a subset of electrodes (de Munck et al., 2007), using a bilinear (Labounek et al., 2019) or a trilinear (Marecek et al., 2016) band-limited power decomposition approach. To take into consideration such variability of the hemodynamic response function (HRF) across subjects, brain areas, and components, we included three regressors into the general linear model for each component separately: a component time series convolved with the canonical hemodynamic response function as the first regressor β₁, and a component time series convolved with the first (temporal) and the second (dispersion) derivative of the canonical hemodynamic response function, as the regressors β₂ and β₃—a procedure proposed and examined in previous event-related studies (Lindquist et al., 2009; Friston et al., 1998). In contrast to Labounek et al. (2019); Marecek et al. (2016), where authors employed F-statistics test inference for the GLM, our aim was to utilize such inference method specifically to discern the direction (positive/negative correlation) of the relationship between the component BLP time series and the BOLD signal. Therefore, we implemented a formula proposed in Calhoun et al. (2004):

$\begin{array}{ll}H=sign({\widehat{\beta }}_1)\sqrt{{\widehat{\beta }}_1^2+{\widehat{\beta }}_2^2+{\widehat{\beta }}_3^2},& (3)\end{array}$

where H is an amplitude combining ${\widehat{\beta }}_1$, ${\widehat{\beta }}_2$, ${\widehat{\beta }}_3$ absolute values, and the directionality is determined by the sign of ${\widehat{\beta }}_1$, that is, the regression coefficient of component time series convolved with the canonical HRF. This combined voxel-wise beta coefficient H was subsequently used in a group-level analysis. Group-level statistics was performed by a one-sample t-test at each voxel, and a cluster-based permutation statistics (Maris and Oostenveld, 2007) was applied to correct for multiple comparisons (α = 0.05, α_cluster = 0.05). Those fMRI BOLD activation maps (later referred to as BOLD signatures), as well as source–space signatures, were visualized with the BrainNet Viewer tool (Xia et al., 2013).

2.9 Spatio-spectral EEG-fMRI integration evaluation

One of the challenges in resting-state EEG-fMRI research is evaluating the reliability of the correspondence between both modalities. On the one hand, there are methodological studies analyzing relatively small EEG-fMRI datasets in terms of a number of subjects (Abreu et al., 2020; Sockeel et al., 2016; Yuan et al., 2016); thus, statistical power is low. On the other hand, there are studies examining reliability and reproducibility with a relatively larger amount of subjects (Bridwell et al., 2018) or across more datasets, and also across different paradigms (Labounek et al., 2019). Since we do not possess two equivalent EEG-fMRI datasets, we opted for checking the reproducibility of observed patterns and relationships by splitting the sample - in our case, we split our dataset into two equal subsets, that is, 36 subjects each. The choice to divide the data into two equal subsets, was guided by practical considerations aimed at balancing statistical power with the need for subgroup analysis. This particular split allowed us to maintain a sufficiently large sample size in each subset, which is essential for observing typically weak levels of EEG-BOLD correlations while facilitating comparisons between the two groups. We want to highlight that the methods described in the following parts can be broadly applied to other EEG-fMRI datasets acquired under different acquisition parameters. However, generalizing the findings and results presented in this study should be approached with caution since they originate from one yet comprehensive EEG-fMRI dataset.

2.9.1 Between-group reproducibility of spatio-spectral EEG patterns

For an initial check of the robustness of the spatio-spectral patterns, we ran the SSP identification algorithm on each subset separately. We then computed the Spearman correlation between each SSP from the first subset and each SSP obtained from the second subset, forming a square matrix with dimensions C × C.

A natural challenge in evaluating the reproducibility of EEG or fMRI components is that the components provided have the in-principle arbitrary order and are not thus directly comparable. This can be solved by sorting them by selecting the best match to each from a predefined list of templates (template matching), or by simultaneously optimizing pairwise matches between all components of both subsets. We chose this latter approach because it does not require selecting a template. In particular, we first transformed this matrix by calculating absolute values and then multiplying by –1 to create a cost function for the Munkres algorithm (Bourgeois and Lassalle, 1971), and then applied a MATLAB implementation of Munkres algorithm to find C unique pairs of SSPs between subsets. The algorithm satisfies a condition of maximum similarity and ensures that each component appears exactly once.

Note that establishing statistical significance for these maximum values is highly non-trivial, due to both the intricate sample dependencies, multiple testing, and ultimately the maximization procedure involved in selecting the best matches. Hence, we present the median and interquartile range of similarity in the results section. The results of the between-group reproducibility of SSPs evaluation will be presented in Section 3.1, along with the reproducibility analysis of spatio-spectral EEG-fMRI integration. Indeed, while both the template matching and the global optimization procedures have previously been commonly applied to provide unique matching, the obtained best matches naturally suffer from upward bias related to the problem of overfitting by optimizing the matches, implicitly assuming a one-to-one mapping, that might be in many scenarios unrealistic, or carry biases in the selection of the templates (if used). Thus, wherever practical, we developed alternative reliability testing procedures that avoid explicit template matching. In the following sections, we describe all statistical procedures used in this study for evaluation of the proposed integration approach.

2.9.2 Reproducibility of spatio-spectral EEG-fMRI integration

While the reproducibility of the Spatio-spectral EEG decomposition per se is not straightforward to establish (see previous section), our main interest here lies in the EEG decomposition for fusion with fMRI. Thus, we introduce a novel approach, where we focus on the reproducibility of the link between an EEG component and its BOLD correlates, that is, of the EEG-fMRI patterns. By EEG-fMRI pattern, we mean the spatio-spectral and temporal signatures of a given EEG BLP IC together with the BOLD signature of this EEG BLP IC, that is, the statistical GLM map obtained when using the EEG BLP IC time course as a regressor for concurrently measured voxel-wise BOLD data. We work with the same subsets as in the previous case. For both subsets, a separate source- as well as electrode-space spatio-spectral decomposition was performed in the same manner as described in Sections 2.6, 2.7, respectively. We test the hypothesis that components with similar spatio-spectral signatures (between the test and retest decomposition) also have a similar pattern in explained BOLD signal, that is, BOLD signatures. To this end, a correlation matrix was computed between the SSPs of the first and the second subsets. Besides the correlation matrix of the whole length of SSPs, separate band-wise correlation matrices were also computed to test band-specific reproducibility. After that, a group-level statistical GLM matrix with beta coefficients as columns for each component was created for both subsets, and a correlation matrix expressing similarity between components was computed. If the EEG components and their correspondence with the BOLD signal in both subsets are similar, then those two correlation matrices should be more similar than by chance, that is, the strength of the match between the EEG SSPs should predict the strength of the match between their BOLD signatures. To test this hypothesis, we implemented permutation statistical testing. During each of 1,000 iterations, columns and rows of spatio-spectral correlation matrix were randomly permuted, and vectorized forms of the permuted SPP correlation matrix and the GLM map correlation matrix were correlated to generate a permutation null distribution. Statistical significance was determined based on a percentile of original not permuted similarity between spatio-spectral and statistical GLM correlation matrices (α = 0.05, right-sided test). A schematic flowchart of the proposed testing of EEG-fMRI pattern reproducibility is shown in Figure 2. The results of this reproducibility test are presented in Section 3.1.

Figure 2

Figure 2. A schematic flowchart of a permutation statistical testing of EEG-fMRI patterns reproducibility between subsets of the dataset.

2.10 Identification of reproducible EEG-fMRI patterns

Following the assessment of overall reproducibility in EEG-fMRI integration, our goal is to identify the most reproducible EEG-fMRI patterns across subsets that would call for more detailed neuroscientific interpretation. To achieve this, we introduce a heuristic metric to identify EEG-fMRI patterns that significantly contribute to the observed reproducibility. This metric considers both the similarity of EEG BLP spatio-spectral maps and the similarity of the corresponding BOLD signatures. For each pair of spatio-spectral components within and across both subsets, we computed the following multiplicative criterion to quantify their similarity:

$\begin{array}{ll}{R}_{M_{x,y}}=\{\begin{array}{ll}\sqrt{R_{I{C}_{x,y}}{R}_{GL{M}_{x,y}}},\hfill & if{R}_{I{C}_{x,y}}{R}_{GL{M}_{x,y}}>=0.\hfill \\ 0,\hfill & otherwise.\hfill \end{array}& (4)\end{array}$

where _RMx,y is the similarity (by multiplicative criterion) between components x and y, _RICx,y is the Spearman correlation coefficient between spatio-spectral signatures of components x and y, and R_GLMx,y is the Spearman correlation coefficient between BOLD signatures of components x and y. This similarity metric was further transformed to a distance metric D = 1 − R_Mx,y for the application of a hierarchical clustering algorithm to form an agglomerative hierarchical cluster tree with the shortest distance method to compute a distance between clusters. We aimed to report only the first five most reproducible EEG-fMRI patterns (which corresponded to the threshold for defining clusters D_thresh = 0.4 for electrode and source spaces). To associate a single whole dataset EEG-fMRI pattern to a pair of EEG-fMRI patterns from subsets 1 and 2, we applied the same criterion as in Equation 4 to both whole dataset EEG-fMRI patterns with defined EEG-fMRI pattern from subsets 1 and 2. The best-matching whole dataset pattern was determined to have the highest product value for these two criteria. The resulting EEG-fMRI patterns are presented in Section 3.2.

Up to this point, reproducibility testing and identifying the most reproducible EEG-fMRI patterns were performed solely on the simultaneous EEG-fMRI data. To investigate whether the observed spatio-spectral components can also be found in the EEG data outside the MRI environment, we performed the following analysis: For each spatio-spectral pattern of the five most reproducible EEG-fMRI patterns, the most similar SSP in the out-of-scanner dataset–based on the highest correlation of the spatio-spectral signatures-was determined and reported. It is important to stress that reporting the most similar patterns from the out-of-scanner dataset does not provide a strict validation or generalization of the results to other datasets due to the challenging problem of overfitting through the matching procedure. However, it demonstrates the level of variability to expect in resting-state EEG data.

2.11 Source–space EEG patterns spatio-temporal relation to BOLD signatures and BOLD RSNs

So far, no one has directly spatio-temporally compared BLP-based features (in our case, SSPs) obtained from EEG with their activation maps in BOLD (BOLD signatures). Even though it seems (at first sight) natural that temporally related EEG and BOLD features should also be similar spatially, there is no strong evidence in the literature that this holds. Furthermore, whether BLP-based features exhibit spatio-temporal correspondence with BOLD RSNs remains unclear. Therefore, we introduce three novel statistical testing procedures, directly corresponding to the three hypotheses regarding spatio-temporal relationships between EEG and BOLD data: (1) whether EEG SSPs are spatially co-localized with their corresponding BOLD activation maps, (2) whether EEG SSPs explaining more BOLD variability exhibit stronger spatial alignment, and (3) whether EEG SSPs are spatio-temporally associated to BOLD RSNs.

2.11.1 Spatial similarity of EEG spatio-spectral and BOLD signatures

Here, for the five most reproducible patterns, we tested whether the similarity of their spatio-spectral signatures to their BOLD signatures is higher than to BOLD signatures of other EEG-fMRI patterns, indicating a potential spatial correspondence between spatio-spectral signatures and temporally coactivated areas observed in the BOLD data. To be able to compare the EEG spatio-spectral signatures and, in this case, purely spatial BOLD signatures, we separately correlated each band-wise part of SSP with a given BOLD signature. Then we computed a weighted average (weighted by the contribution of each band to the overall SSP weights) of the absolute values of Spearman correlations. In that manner, we computed a complete (weighted average) correlation matrix of all SSPs and BOLD signature combinations. We compared the mean value on the diagonal (i.e., the correct correspondence between SSPs maps and BOLD signatures) to the null distribution obtained by computing the same statistic for each of 1,000 random permutations of columns of the matrix, representing random assignment of different BOLD signatures to SSPs. The significance level α was set to 0.05 (right-sided test). We further evaluated each component separately and how similar it is to its corresponding GLM map compared to other BOLD signatures (each row of the correlation matrix). On the significance level α = 0.05 (right-sided test, uncorrected for a number of components), we report the most similar components. A schematic flowchart of the statistical testing is included in the Supplementary Figure S1. The results are presented in Section 3.3.

2.11.2 spatio-temporal similarity of EEG spatio-spectral and BOLD signatures

Furthermore, we tested whether the EEG components that are more temporally synchronized with the BOLD signal also tend to be more spatially similar to their BOLD signatures. It is reasonable to assume that multiple of the 30 EEG components may not correspond to strong stable neuronal dynamics but rather to some transients or residual artifacts. To the extent that both BOLD and EEG reflect predominantly local neuronal activity (an assumption commonly used in brain activity modeling, based on experimental observations such as those presented by Logothetis et al. (2001) and others), one would expect that the reproducible EEG patterns would have BOLD signatures co-localized to the EEG spatial patterns, while the weaker, less robust, or artifactual components not temporally related to BOLD signal, would also have BOLD signature unrelated to the spatial EEG pattern.

To test this hypothesis, we again use a permutation-based test. Similarly to the previous subsection, we compute the Spearman correlation between each of the 30 EEG spatio-spectral components and its BOLD signature, resulting in a vector. Then, for each EEG SSP, we correlate its time series with the BOLD signature time series. Since the BOLD signature is a statistical map representing spatial activation, it does not inherently possess a time series. Therefore, to establish a comparable time series for each BOLD signature, we computed the average BOLD time series. This average was computed using a weighted approach, where the weights corresponded to the BOLD signature's spatial pattern. This means we prioritized voxels that exhibited a stronger correlation with the EEG SSP time series when calculating this weighted average. This process also resulted in a vector of temporal similarities, with a length corresponding to the number of components, specifically 30. If the described spatio-temporal similarity is valid, these two vectors should exhibit greater similarity than what would occur by chance. The null distribution was obtained by randomly permuting elements of both vectors (1,000 iterations) and the original not permuted absolute value of Spearman correlation was tested against the permutation distribution at the significance level α = 0.05 (right-sided test). A schematic flowchart of the statistical testing is included in the Supplementary Figure S2. The results are presented along with the previous test results in Section 3.3.

2.11.3 Source–space EEG patterns spatio-temporal relation to BOLD RSNs

In previous EEG-only studies (Liu et al., 2017, 2018; Sockeel et al., 2016) and EEG-fMRI studies (Yuan et al., 2016; Abreu et al., 2020; Marino et al., 2019) the authors aimed to link EEG-derived RSNs and BOLD RSNs both spatially and temporally. We have decided to investigate the spatial and temporal correspondence between EEG SSPs and BOLD RSNs involving a broad range of RSNs in the same unbiased spirit as for the previous tests. At first, we performed a spatial group ICA decomposition of the BOLD data. To this end, we combined all of the single-subject z-score normalized datasets across the temporal domain. The dimensionality was then reduced to 30 PCA components before running spatial ICA decomposition by the FastICA algorithm (Hyvärinen and Oja, 2000). Subsequently, ICs representing RSNs were identified based on a functional network atlas obtained from (Shirer et al., 2012). The criterion was based on the Spearman correlation between the spatial distribution of the IC weights and atlas masks. Subsequently, statistical testing was carried out in an analogous fashion for the spatio-temporal similarity of EEG SSPs and BOLD signatures. In particular, the (weighted across-band-averaged) correlation matrix of spatio-spectral signatures in the source-reconstructed space on one side with the spatial signatures of BOLD ICs spatial maps on the other side was computed. After that, temporal signatures of SSPs convolved with the canonical HRF were correlated with the BOLD ICs temporal signatures. If there was spatial and temporal correspondence between EEG-derived SSPs and BOLD IC components representing RSNs, then those two correlation matrices should be more similar than by chance. Similar to the approach in Section 2.9.2, we implemented permutation statistical testing. In each of the 1,000 iterations, columns and rows of spatial patterns correlation matrix were randomly permuted (by two different permutations), and vectorized forms of the permuted spatial patterns correlation matrix and temporal patterns correlation matrix (absolute values) were correlated to generate a permutation null distribution. The statistical significance was determined based on the percentile of the permutation distribution occupied by the original data similarity (Spearman correlation) between spatial and temporal patterns correlation matrices (α = 0.05, right-sided test). A schematic flowchart of statistical testing of EEG and BOLD ICs similarity testing is in Figure 3. The results are presented along with the two previous test results in Section 3.3.

Figure 3

Figure 3. A schematic flowchart of permutation statistical testing of the correspondence of spatial and temporal similarities between the EEG source–space spatio-spectral components and the BOLD ICs.

3 Results

The results section is divided into three subsections. In Section 3.1 we report EEG-fMRI patterns reproducibility for source- and electrode-space approaches, in Section 3.2, we show the reproducible EEG-fMRI patterns in both spaces, in Section 3.3, a summary from spatial and temporal correspondence testing of EEG SSPs, their BOLD signatures, and BOLD ICs is reported.

3.1 Stability of SSPs and reproducibility of EEG-fMRI integration

3.1.1 Numerical SSP decomposition stability

As described in Section 2.6, first, the spatio-spectral components stability testing in terms of ICA decomposition was performed utilizing the ICASSO tool. Based on a similarity index and visual assessment from the 2D canonical correlation analysis (CCA) projection provided by the ICASSO tool, we conclude that the ICA decomposition is algorithmically stable for both source- and electrode-space approaches as well for both subsets and the whole dataset. The ICASSO outputs are included in the Supplementary Figure S3.

3.1.2 EEG-fMRI pattern reproducibility

The EEG-fMRI patterns proved generally reproducible in the test–retest setting, with the more matching EEG SSPs typically corresponding to those with a more similar BOLD signature. This was true for both EEG spaces, with marginally better reproducibility for the source–space approach, in particular only two statistically non-significant results were observed for the electrode-space band-wise patterns in δ (r_s = 0.03) and γ (r_s = 0.04) bands. The highest reproducibility was observed for the multiband (whole) SSPs—(r_s = 0.29 for the source–space and r_s = 0.25 for the electrode-space). The results are summarized in Table 1.

Table 1

Table 1. Reproducibility of EEG-fMRI patterns: match of EEG SSPs corresponds to match of their BOLD GLM signatures.

Of note, the SSPs reproducibility test described in Section 2.9.1 and based on components matching by Munkres algorithm clearly showed a biased result with median Spearman correlation 0.46 (IQR= 0.42, 0.58). This approach has two shortcomings: (1) It is difficult to test the level of similarity between matched SSPs statistically and (2) It does not test BOLD signature reproducibility. These results are further elaborated in the Discussion Section 4.1.

3.2 Reproducible EEG-fMRI patterns

We identified the five most reproducible EEG-fMRI patterns based on the similarity measure and clustering described in Section 2.9.2. The hierarchical tree is illustrated in Supplementary Figure S4. The reproducibility is meant in the sense of similarity of SSPs and also of BOLD signatures. The five most reproducible source- as well as electrode-space EEG-fMRI patterns are depicted in Figures 4, 5, respectively (ordered based on descending similarity defined in Section 2.10).

Figure 4

Figure 4. Visualization of the most reproducible source–space EEG-fMRI patterns (A–E) in terms of the band-wise Euclidean norm of an SSP, SSPs, and the BOLD signatures (from left to right). For each EEG-fMRI pattern, only selected band-wise SSPs are visualized based on the weight distribution and their mutual similarity. Pattern (A) is similar between all frequency bands, pattern (B) is the most prominent in the α band, pattern (C) is similar in all β subbands and has opposite weights in γ band, pattern (D) in changes gradually from δ to γ) band. Finally, the pattern (E) changes its spatial distribution from δ to γ band.

Figure 5

Figure 5. Visualization of the most reproducible electrode-space EEG-fMRI patterns (A–E) in terms of the band-wise Euclidean norm of an SSP, SSPs, and the BOLD signatures (from left to right). For each EEG-fMRI pattern, only selected band-wise SSPs are visualized based on the weight distribution and their mutual similarity. Pattern (A) is the most prominent in the α band, pattern (B) is similar between all frequency bands, pattern (C) is similar for θ band and neighboring frequencies and differs in the γ band, pattern (D) is less prominent in δ and θ bands and more in high (β₂ and β₃) frequency bands. Finally, the pattern (E) changes its spatial distribution from δ to γ bands.

All these reproducible source- as well as electrode-space EEG-fMRI patterns have statistically significant correlation (cluster-based permutation test, α = 0.05, α_cluster = 0.05) with the BOLD signal (see BOLD signatures in Figures 4, 5). Moreover, they differ substantially in frequency as well as band-wise spatial distributions. Furthermore, some of the source–space patterns are noticeably similar to electrode–space patterns in terms of their BOLD signatures as well as ICs weights distribution.

3.2.1 EEG-fMRI patterns common to source and electrode-space

Specifically, the first EEG-fMRI pattern (Figure 4A) is neither frequency- nor spatially specific and is represented by positive weights across all frequency bands and source positions, that is, widespread positive weights from δ to γ bands focused more frontally and medially and less occipitally and parietally. This global pattern does not correlate negatively with the BOLD signal but it shows a positive correlation pattern with the BOLD signal in temporal (Heschl's, superior temporal gyri), frontal (subcentral, frontal middle, frontal superior medial gyri, and supplementary motor area), parietal (supramarginal gyrus and precuneus), and occipital (lingual, occipital inferior gyri) lobes. Positive clusters also occupy the insular and cingulate cortex and some subcortical brain areas such as the amygdala, caudate nucleus, thalamus, and putamen. In several parts of the cerebellum (not visualized), positive clusters can also be observed. Noticeably, the second electrode-space EEG-fMRI pattern (Figure 5B) can also be called global since its weight distribution is broad and very similar across all frequency bands. This pattern correlates only positively with the BOLD signal in very similar brain areas as the mentioned pattern in Figure 4A, although in contrast to the source–space global pattern, the electrode-space positive clusters were not observed in middle frontal gyrus and were in fusiform gyrus.

The following pattern (Figure 4B) has weights distributed mostly around the alpha band and spatially located positively across occipital, parietal, and partly temporal lobes. This pattern has widespread negative correlation with the BOLD signal across the whole cortex and frontally focused positive activations. Negative clusters can be found almost over the whole occipital lobe (calcarine, lingual, inferior, middle, superior gyri, and cuneus), parts of the parietal (postcentral, supramarginal gyri, parietal superior cortex, and paracentral lobe), and temporal lobes (middle, superior temporal, Heschl's, and fusiform gyri). In frontal lobe negative clusters are located in subcentral, precentral, inferior frontal gyri, supplementary motor area, rolandic and inferior frontal operculum. Positive clusters can be observed in anterior cingulate as well as parts of orbitofrontal cortices. This source-space pattern is highly similar to the electrode-space EEG-fMRI pattern shown in Figure 5A. This pattern is spatially distributed over the whole scalp with the maximum in the parieto-occipital sites and also partly frontally. The statistical GLM map is again very similar to its source–space counterpart; only the electrode-space component does not have positive clusters in the frontal lobe, while in some parts of the cerebellum, positive clusters were observed.

3.2.2 Source–space specific EEG-fMRI patterns

The third EEG-fMRI pattern (Figure 4C) seems to be source–space specific and is characterized by mostly frontoparietal and medial positive weights across all three β subbands and negative bilateral (mostly temporal) γ band negative weights. This pattern correlates positively as well as negatively with the BOLD signal. The negative cluster can be observed in the parietal (postcentral, supramarginal gyri, inferior cortex, and paracentral lobe) lobe, temporal superior gyrus, and frontal (precentral, subcentral gyri, and inferior operculum) lobe. The positive cluster is located in the precuneus as well as again in several parts of the cerebellum.

The next EEG-fMRI pattern (Figure 4D) is characterized by negative parietal and occipital weights for all frequency bands and positive bilateral temporal and frontal weights from β₁ to γ bands. This pattern shows only negative BOLD correlate clusters in the parietal (postcentral, supramarginal gyri, inferior cortex, and paracentral lobe) as well as in frontal (precentral, middle, superior, subcentral gyri, and inferior operculum) lobes.

Finally, the last EEG-fMRI pattern in the source–space (Figure 4E) is characterized by negative weights in low frequencies in bilateral temporal and occipital sites and, at the same time, positive γ weights again on bilateral temporal and occipital cortex sites. This pattern correlates exclusively negatively in parietal (postcentral, supramarginal gyri, inferior cortex, paracentral lobe, and precuneus) and frontal (precentral, superior, middle and inferior gyri, supplementary motor area, and frontal inferior operculum) lobes as well as in middle part of the cingulate cortex.

3.2.3 Electrode-space specific EEG-fMRI patterns

The third most reproducible electrode-space EEG-fMRI pattern (Figure 4C) is characterized by whole scalp negative weights in low frequencies around θ band in contrast to positive weights in γ band spread mostly on parieto-occipito-temporal electrode sites. This SSP correlates positively with the BOLD signal in temporal (Heschl's, superior, and middle temporal gyri), frontal (supplementary motor area, rolandic, frontal inferior operculum, frontal inferior triangularis, frontal superior, frontal superior medial, and precentral gyri), parietal (precuneus, postcentral gyrus), and lingual gyrus in occipital lobe. Anterior and middle cingulate, as well as insular cortices and parts of the cerebellum, also show a negative correlation with the BOLD signal.

The fourth scalp SSP (Figure 5D) is again not spatially specific and demonstrates changes in the band-limited power between low and high frequencies. Therefore, we observe gradually changing polarity of weights from the lowest frequencies δ to β₃ widely across the whole scalp. For this pattern, specifically, α and γ frequency bands are the least important. The correlation pattern of the BOLD signal is exclusively negative and very similar to the statistical GLM map of α band pattern (Figure 5A). Apart from α band pattern, negative clusters can also be observed in parietal inferior, insular, middle cingulate cortices, and frontal superior gyrus. A positive cluster is only in some parts of the cerebellum.

The last electrode-space EEG-fMRI pattern (Figure 5E) is relatively complex. From δ to θ band, the spatial pattern is mostly positive in frontal, central, and temporal electrode sites. From α to β₃ bands also, negative parietal and occipital weights are present, and the γ band is represented by a reversed pattern than for all other frequencies. This complex EEG-fMRI pattern has an exclusively positive BOLD signature, and it is very similar to previous negative BOLD signatures in Figures 5A, D, respectively. In contrast to those, clusters in the parietal and frontal lobes have relatively higher T values compared to the occipital ones.

3.2.4 Variability of the five most reproducible EEG-fMRI patterns

Last but not least, to capture the variability of components between subsets of the same dataset as well as between different datasets recorded in a different environment (out-of-scanner dataset), in Figures 6, 7 we show the SSPs from both subsets and the closest SSP in out-of-scanner dataset for each of five the most reproducible EEG-fMRI patterns for source- as well as electrode-space approaches.

Figure 6

Figure 6. Visualization of corresponding subset 1 and 2 SSPs together with out-of-scanner SSP for each of the five most reproducible source–space EEG-fMRI patterns (A–E).

Figure 7

Figure 7. Visualization of corresponding subset 1 and 2 SSPs together with out-of-scanner SSP for each of the five most reproducible electrode-space EEG-fMRI patterns (A–E).

Please note that in this case statistical testing is not provided. Therefore, both figures serve only as an example of SSP's variability that can be expected when applying the methodology to a different dataset. Also, note that the most similar component of a global SSP from the out-of-scanner dataset was not determined based on the highest correlation with a global SSP from the whole dataset but rather by selecting a component having all weights with the same sign. This feature cannot be measured with a similarity measure based on Spearman correlation. Also, the visualized out-of-scanner SSP in Figure 6C is the second most similar component to the whole dataset component since the most similar component was the same as for Figure 6E.

3.3 Source–space EEG patterns spatio-temporal relation to BOLD signatures and BOLD RSNs

We report the findings of the three hypotheses testing regarding spatio-temporal relationships between EEG and BOLD data: (1) whether EEG SSPs are spatially co-localized with their corresponding BOLD activation maps (Methods Section ??), (2) whether EEG SSPs that explain more BOLD variability exhibit stronger spatial alignment (Methods Section ??), and (3) whether EEG SSPs are spatio-temporally associated to BOLD RSNs (Methods Section 2.11.3). The results are summarized in Figure 6.

We found for the first 5 most reproducible EEG-fMRI patterns weak but statistically significant (R = 0.14, p = 0.031) spatial relation between their spatio-spectral signatures and BOLD signatures, suggesting spatial colocalization between them, see Figure 8A. In post-hoc testing, we observed that the similarity was mostly driven by spatio-spectral patterns B and D, see Figure 8B.

Figure 8

Figure 8. Spatial (spatio-temporal) relation of EEG spatio-spectral components and BOLD. (A) Spatial similarity of the EEG spatio-spectral signatures and their corresponding BOLD signatures: mean value across the five most robust components (orange bar), the permutation test statistic (blue histogram). (B) Spatio-temporal similarity of the EEG spatio-spectral signatures and their corresponding BOLD signatures: mean values across all components (orange bar), the permutation test statistic (blue histogram). (C) Spatio-temporal similarity of the EEG spatio-spectral signatures and BOLD RSNs: mean values across all components (orange bar), permutation test statistic (blue histogram). (D) Spatial similarity of all robust source–space spatio-spectral patterns (A-E) all BOLD signatures (violin plots, dots) with denoted actual BOLD GLM map (square point). (E) Pair-wise similarity (Spearman correlation) of EEG spatio-spectral and BOLD signatures spatial (scatter plot x-axis) and temporal (scatter plot y-axis) signatures. Colored points correspond to EEG-fMRI patterns from (D). (F) Pair-wise similarity (Spearman correlation) of EEG spatio-spectral and BOLD RSNs spatial (scatter plot x-axis) and temporal (scatter plot y-axis) signatures. Colored points correspond to EEG-fMRI patterns from (D). Three points for alpha EEG-fMRI pattern are labeled as follows: (1) most spatially similar to Precuneal network (PN), (2) most temporally similar to Visuospatial network (VSN), and (3) most spatio-temporally similar to Primary visual network (PVN).

The suggested spatio-temporal test of EEG spatio-spectral components and BOLD signatures shows statistically significant relation between the levels of their temporal and spatial correlation (R = 0.42, p = 0.009), suggesting that EEG spatio-spectral components having a higher level of correlation in the temporal domain with its statistical GLM map pair are also more similar spatially, see Figures 8B, E.

Statistical testing of the correspondence between EEG-derived spatio-spectral components and BOLD-derived ICs in terms of spatial and temporal similarities (R = 0.03, p = 0.255) did not reveal any similarity, see Figure 8C, F. Specifically to demonstrate this on the most reliable alpha pattern, we picked 3 pairs of alpha-RSNs as: (1) The most spatially similar (PN), (2) The most temporally similar (VSN), and (3) The most spatio-temporally similar (PVN) (up-right in scatter plot) pair alpha pattern-RSNs, see Figure 8F. In Figure 9, the alpha EEG-fMRI pattern together with those 3 RSNs is visualized. These results are commented and put in the context in the Discussion Section 4.3.

Figure 9

Figure 9. Example of spatial mismatch between alpha EEG-fMRI pattern and BOLD RSNs.

4 Discussion

In the present study, we, for the first time, introduce a fully data-driven method of spatio-spectral decomposition in the source-reconstructed space applied to hdEEG data recorded simultaneously with the fMRI. We directly compare it to the electrode-space spatio-spectral decomposition on the same dataset. Rather than solely testing against other traditional methods, our study is focused on addressing intriguing questions regarding the EEG-fMRI relationship. Therefore, we apply a novel, unbiased approach to evaluate the reproducibility of the observed EEG-fMRI patterns and, maybe even more importantly, of their spatio-temporal relation to their BOLD signatures and BOLD RSNs. A comparison with the SSPs obtained from the out-of-scanner data and with the results of electrode-space decomposition is also provided.

4.1 Stability of SSPs and their relation to fMRI

The unbiased reproducibility testing based on the split-half analysis of the large EEG-fMRI dataset clearly showed that the SSPs, and at the same time, their correspondence to the BOLD signal, are reproducible between subsets for both approaches (source- and electrode-space). All frequency-specific (single-band) EEG spatial signatures are reproducible across subsets (besides electrode-space δ and γ bands, see Table 1). The whole SSPs robustness suggests that the spatio-spectral components are formed across multiple frequency bands, as in the case of most electrode-space (Figures 5B–E) as well as source space (Figures 4A, C–E) SSPs. Additionally, (single-band) SSP reproducibility does not exclude the possibility that some SSPs might be narrow (defined by a single frequency band) and that most frequency bands have at least a few unique components.

It is also interesting to note that the highest reproducibility is achieved by the full (multiband) spatio-spectral signatures compared to band-wise signatures for both spaces; see Table 1. This suggests that the pattern across joint spatio-spectral space is important, and performing only single-band decomposition might miss some spectrally wide patterns.

In this study, the concurrent measurement with BOLD provided that we could simultaneously test SSPs and their BOLD signatures. Previously, when assessing the stability/reproducibility of EEG patterns between datasets or subjects, many studies utilized template matching procedures to pair the most similar components to BOLD-derived RSNs templates (Liu et al., 2017, 2018; Marino et al., 2019) or component clustering algorithms applied only to EEG signatures (Labounek et al., 2018, 2019). However, such matching approaches suffer from artificially inflated match quality (as the optimal match is sought for across many components, even a random set of candidate maps, if sufficiently large, can provide decent “optimal” matches). For illustrative purposes and to maintain some comparability with results of previous studies (despite many other methodological differences), we also applied such straightforward matching procedure in the evaluation of SSPs reproducibility (the best matching SSPs between two subsets), where we reached a relatively higher mean level of similarity (R = 0.46).

A potential drawback of our reproducibility analysis is the number of PCA components retained before the ICA decomposition, that is, model selection. It is because overestimating components may result in a component splitting (Abou-Elseoud et al., 2010). This could decrease the level of reproducibility. It is difficult to adequately estimate the number of PCA components, and various model selection criteria often provide different estimates. Some authors estimate the number of PCA components on BLP matrices based on information theory criteria such as minimum description length (MDL) (Liu et al., 2017, 2018; Marino et al., 2019). Others Labounek et al. (2018) arbitrarily chose a number of ICs and tested the stability of recovered ICs by metrics obtained from toolboxes such as ICASSO software (Himberg et al., 2004) or setting the number of ICs to be around the typical number from fMRI studies (Li et al., 2018). A combination of the last two mentioned model selection strategies was utilized in the present study. A different approach might be to estimate the number of non-random components via a data-modeling and testing against surrogates (Vejmelka et al., 2015). Generally, the number of ICs throughout the literature fluctuates from 20 to 50; therefore, we chose the number of components in the current study to be 30. Another improvement might be a generalization of the reproducibility evaluation between different datasets or even different paradigms, not only resting-state (Labounek et al., 2018, 2019).

4.2 Reproducible EEG-fMRI patterns

While we chose, in line with the literature, a generous number of 30 ICA components to be extracted, unlike several previous studies we present only a relatively modest number of 5 of them. Two factors may play a role here. First, we pose an additional requirement on the presented components, that not only they have a counterpart in the decomposition of an independent EEG dataset, that is, the EEG SSP patterns is reproducible, but moreover, we require that its BOLD signature is also conserved between the test and retest dataset. This focuses our analysis only on those EEG patterns that do have some robust BOLD signature, but perhaps even more importantly provides a control for false matches between random EEG patterns, which are likely to occur due to the sheer number of available components for the matching. Secondly, the fact that we search for multiband spectral patterns that do not require the same spatial signature for each frequency means that the approach presented here can encompass richer patterns within a single component, rather than dividing them into individual frequencies.

Thus, the resulting number of components described in detail in both source- and electrode-space are relatively smaller than those reported by other authors. Previous studies usually concentrated on the single-band source–space BLP decomposition and directly focused on spatial and temporal correspondence. They typically report almost a complete set of recovered EEG-derived RSN patterns. In EEG-only studies (Liu et al., 2017, 2018), the authors reported finding a complete set of 14 EEG-BLP patterns spatially similar to BOLD-derived RSNs. Conversely, Marino et al. (2019) focused on extracting a single DMN pattern from the source–space EEG BLP matrix.

In Bridwell et al. (2013), the authors reported 10 robust EEG SSPs performing electrode-space spatio-spectral decomposition. Five of those components were specifically spectrally focused around the α band. Unlike our approach, the authors had a finer spectral resolution of 0.5 Hz, which, in combination with ICA decomposition along the spatio-spectral domain, may result in the decomposition of one alpha rhythm into several components. All five alpha components exhibited remarkably similar occipital SSPs in their study, and all demonstrated a negative correlation with the BOLD signal in the occipital, parietal, and certain regions of the frontal lobe. This suggests that they may relate to the over-splitting of a single, relatively spatio-temporally consistent alpha component. The authors also report a widespread positive correlation across the cortex for the rest of the components, where three were in δ–θ and two in β₂–γ frequency bands.

In Labounek et al. (2018) and Labounek et al. (2019), the authors found 14 SSPs (across three different datasets) having substantial similarity with components from Bridwell et al. (2013). We present the five most reliable EEG-fMRI patterns that are reproducible in terms of both SSPs and their BOLD signatures. Furthermore, it is noteworthy that some EEG-fMRI patterns between spaces (electrode and source level) intuitively correspond to each other. In the following paragraphs, we will discuss all reported patterns regarding interpretation and connect them to existing studies.

The most reproducible electrode-space EEG-fMRI pattern (Figure 5A) resembles the alpha components from a couple of previous studies (Bridwell et al., 2013; Labounek et al., 2018). It is a data-driven derived occipital alpha correlate known from many previous studies (Goldman et al., 2002; Moosmann et al., 2003; Feige et al., 2005; Laufs et al., 2003a,b; Tyvaert et al., 2008; Gonçalves et al., 2006). In agreement with Bridwell et al. (2013), we found widespread negative correlation within a middle frontal, superior temporal, inferior occipital, lingual gyri, and cuneus; see statistical GLM map in Figure 5A. Note that we also observed positive clusters in the anterior cingulate and parts of the orbitofrontal cortices. Positive correlations are not widely reported across studies except (Laufs et al., 2003b), where the authors reported a positive correlation with the BOLD signal in the anterior cingulate cortex in a fixed effects group analysis. Intuitively, this electrode-space component has its source–space counterpart as EEG-fMRI pattern in Figure 4B. This source–space pattern has very similar statistical GLM map patterns corresponding to sensory areas and weights distribution, mostly in α band parieto-occipital cortical sites. Since the source–space SSPs were not studied before, linking them with the current literature is difficult. Nevertheless, compared to the TFICA approach on EEG-only dataset in Li et al. (2018), this SSP may be related to the group of components called Default by the authors. It is also interesting that the BOLD signature of those components resembles a more typical EEG-fMRI occipital α band correlate rather than the bilateral frontoparietal pattern reported and discussed in the studies by Laufs et al. (2003b) and Laufs et al. (2006).

The second pattern in the electrode-space shown in Figure 5B has (as in the previous case) its intuitive counterpart in the source space; see Figure 4B. Those global patterns in both spaces also share very similar BOLD signatures when having exclusively positive clusters in all lobes and several subcortical structures such as the thalamus, amygdala, or putamen. To further elucidate the existence of this pattern, it probably represents a global BLP mean time series since they appear to be correlated in space and across all frequency bands. This pattern may be notable in studies where authors derive BLP regressors as average BLP time series across all electrodes (Mantini et al., 2007). Our results show that the average BLP time series is one of the most stable (assessed by ICASSO tool) BLP-independent components. Furthermore, we demonstrated that the average component belongs to a set of robust EEG components, significantly contributing to EEG-fMRI integration. Even though it is difficult to interpret such an unspecific EEG-fMRI pattern, this pattern might represent both neurophysiological activity and/or globally expressing artifactual patterns. The former (neurophysiological activity) explanation can be supported by the study Magri et al. (2012). Despite substantial methodological differences (resting-state simultaneous LFP and BOLD data in monkey visual cortex), the authors made notable conclusions and remarks linking to this global as well as the previous alpha pattern. First, they found strong positive correlations between LFP at frequencies above 20 Hz and much weaker negative correlations between low-frequency LFP (<20 Hz) and the BOLD signal. Second, they report positive covariation of α LFP power and BOLD with total LFP power. Furthermore, they hypothesize that when significant fluctuations in total LFP power take place, the correlation between alpha power and the BOLD signal will seem positive (as reported by some authors), reflecting their shared variability with total LFP power. It seems that utilizing our method, the alpha pattern (anticorrelating with the BOLD signal) can be clearly distinguished from total power fluctuation (global pattern) allowing us to study those phenomena separately. Our method is, therefore, superior to ones considering only temporal averages of BLP across preselected electrodes/ROIs.

The next source–space EEG-fMRI pattern in Figure 5C represents an interplay between β and γ frequency bands. The SSP for all β subbands can be described as primarily frontoparietal medial positive weights compared to negative bilateral temporal and parietal patterns of the γ band. Notably, the BOLD signature of this pattern spatially resembles its SSP in terms of the positive cluster in the precuneus. Furthermore, parietal and temporal negative activations can be spatially linked to the γ band as well. As the positive cluster in the precuneus (the central spatial component of the DMN) together with negative clusters (mainly SMN and AN functional networks), this SSP might represent temporal anticorrelation structure between DMN and so-called task-positive networks as suggested in Fox et al. (2005) and later discussed in many other studies.

Another source–space EEG-fMRI pattern in Figure 4D reflects a broadband pattern from δ to γ band where negative weights can be observed through all bands in parietal and occipital lobes, whereas bilateral positive weights in temporal and frontal lobes are observed in higher β₁ to γ frequency bands. Based on its BOLD signature, this pattern might capture active somatosensory brain areas processing information; thus, positive correspondence between the BOLD signal and β₁—γ weights. At the same time, a decrease of BLP in parieto-occipital cortex sites might be associated with typical α band suppression and probably in combination with inactivity of visual processing since there is no positive cluster in occipital areas. Therefore, this pattern could be considered somatosensory specific and could distinguish between visual and other sensory processing.

The last not mentioned source–space EEG-fMRI pattern relating low frequencies around θ and γ band is shown in Figure 4E. Negative, mostly occipital and midline low-frequency band weights are included together with γ band bilateral frontal, temporal, and parietal positive weights. This pattern exhibits correlation with the BOLD signal in sensory areas such as postcentral and supramarginal gyri, parietal inferior cortex, and paracentral lobe. This pattern might, along with other EEG-fMRI patterns, reflect somatosensory processing.

The electrode–space pattern in Figure 5C also represents mutually anticorrelated BLP patterns between negative weights in low frequencies with a maximum around θ band and positive γ band with a maximum in right parietal as well as left temporal electrode sites. Again, this pattern is not reported in previous studies. It is interesting to note that its BOLD signature differs from other EEG-fMRI patterns by having negative clusters in several parts of the frontal lobe, mainly in the superior frontal and medial superior frontal gyri, which might have a connection to typical sources of the midline theta rhythm (Asada et al., 1999).

The EEG-fMRI pattern in Figure 5D is again a frequency-wide but spatially not specific pattern showing a global anticorrelation pattern of BLP fluctuations between high and low frequencies. Interestingly, both α and γ bands are the least important frequency bands for this pattern. This pattern is not observable in spatio-spectral decomposition studies (Bridwell et al., 2013; Labounek et al., 2018, 2019). The rationale might be that it is unlikely that ICA decomposition performed along a spatio-spectral domain will retrieve such type of frequency-wide pattern. In contrast, performing ICA decomposition along a temporal domain might reveal such types of spatio-spectral signature patterns. A more related frequency-wide component may be found in the trilinear PARAFAC decomposition, such as in Marecek et al. (2016) and Mareček et al. (2017), even though it is a different type of decomposition. Despite this component differing in spatio-spectral signature from SSP in Figure 5A, the fMRI correlation pattern exhibits a statistically significant negative relation when changing from high to low frequencies. This SSP also may be interpreted reversely such that large weights in lower frequencies δ and θ bands positively correlate with the BOLD signal across sensory brain areas in occipital and parietal lobes, which is in correspondence with low-frequency components in Bridwell et al. (2013); Labounek et al. (2018).

Finally, the last not discussed electrode-space EEG-fMRI pattern is a broadband and spatially specific pattern in Figure 5E. Its SSP is represented by (1) frontal, central, and temporal positive weights for low frequencies; (2) α band negative parieto-occipital weights, and (3) reversed pattern in γ band. The statistical GLM map of this pattern is very similar (sign reversed) to one in Figure 5A, except that clusters in parietal somatosensory areas have relatively higher T-values compared to occipital ones.

Examples of SSP's variability in Figures 6, 7 show that to some extent similar components can be obtained not only from different subsets of the same dataset but even from a different dataset recorded out of the MRI scanner environment. Even though not statistically tested, this observation suggests that a similar methodology may also be applied to out-of-scanner datasets, which potentially increases the number of experiment types that the proposed methodology can investigate. Subjectively, the global and the alpha SSPs are the most robust for both source- and electrode-spaces. As highlighted in several sections of our manuscript, the out-of-scanner dataset differed from the simultaneous one in terms of how we acquired and analyzed the data. We acknowledge that these differences could have a meaningful impact on the obtained SSPs and, consequently, on the similarity levels of matched SSPs between the datasets. However, delving into a detailed exploration of the influence of these analysis differences is beyond the scope of our current study. This aspect might become a subject of investigation in future research.

In this work, we focus on the five most robust EEG-fMRI patterns. Indeed, further, particularly less prominent patterns could be extracted by alternative clustering setups (e.g., less strict thresholding) or a combination of dimensionality reduction and clustering algorithms such as applied in Piorecky et al. (2019). A particularly interesting direction would be to consider a possible hierarchical splitting of SSPs into more components and other aspects, such as linking between electrode-space and source–space components. Such an approach might also reveal differences in observed components between data spaces (electrode vs. source) and that way help better understand a source–space version of the spatio-spectral decomposition approach.

The use of source localization before the EEG power decomposition provides a direct anatomical link for spatio-spectral components, which we consider a major benefit of this method. In this study, a source–space spatio-spectral decomposition approach did reveal source–space unique patterns compared to the electrode-space. We observed spatial correspondence between parts of SSP and its BOLD signature for some source–space patterns. In contrast, we should be aware of the ill-posed nature of the source localization and the resulting limitations of this technique (Hallez et al., 2007).

Last but not least, we want to make a comment on inter-individual variability. In this work, we applied a relatively conservative approach to the z-score normalization of each BLP time series. Nevertheless, there may be meaningful variability across subjects, as well as across bands/electrodes, that we here suppressed. This may be an objective of a future study.

4.3 EEG spatio-spectral components relation to BOLD signatures and BOLD RSNs

Based on the results, we observed a far-from-perfect but statistically significant spatial relation between EEG spatio-spectral components and BOLD signatures for the five most reproducible EEG-fMRI patterns. This finding is in line with the common assumption that the BLP fluctuations have spatially (almost) the same source as the BOLD signal. In contrast, by looking at the most similar source–space spatio-spectral signature and its BOLD signature, which is the alpha pattern (Figure 4B), we can still notice that the spatial match is far from perfect. The similarity is mostly driven by both maps' tendency to be more parieto-occipital rather than frontal. Besides all the potential sources of mismatch mentioned in the following paragraphs, the low levels of similarity can be caused by the ill-posed nature of the inverse problem. Indeed, one should consider that even in the case the fMRI and EEG power were perfectly co-localized, the EEG spatial distribution in the brain volume reconstructed from the electrode signals may achieve far-from-perfect recovery of this EEG-BOLD spatial match. The level of this reconstruction-based bias could be at least partially elucidated by simulation studies with realistic signals and forward and backward projection.

Furthermore, a statistically significant linear relationship between temporal and spatial correspondence between EEG spatio-spectral patterns and their BOLD signatures further suggests that there exist EEG components more informative for the BOLD signal, which is (at the same time) more spatially similar to their BOLD signatures. Suppose we assume that this correspondence arises from physiological brain activity. In that case, it means that the more of the same brain activity visible in both EEG and BOLD, the more spatially similar they are and vice versa. On the other side, our statistical evaluation could not confirm the spatial and temporal correspondence between BLP-derived components and the BOLD ICs. This negative finding might contrast with studies (Liu et al., 2017, 2018; Marino et al., 2019). Via a group ICA approach, obtained components represent the structure of the whole dataset. Of note, our approach does not require any pattern matching of components between subjects, and avoids thus potential overfitting.

Our approach is, of course, vulnerable to the ubiquitous problem of selecting the number of ICs recovered and the previously mentioned problem of component splitting, which may cause a problem for the selected similarity statistical method. Nevertheless, our study suggests that while one may find for known BOLD ICA components some EEG SSP patterns that would be spatially similar (as shown previously in the literature), there is not a clear spatio-temporal agreement between main components of EEG BLP and the BOLD ICA components. More research is thus needed to clarify the spatial (spatio-spectral) and temporal correspondence between EEG BLP patterns and BOLD ICs.

For an illustrative example, as demonstrated in Figure 9, we see that even for the most reproducible EEG-fMRI alpha pattern (based on several previously described criteria), there is no clear RSN counterpart that would be the most similar spatially and at the same time temporally. Furthermore, those related networks (precuneal network, visuospatial network, and primary visual network) are associated with a different type of brain processing. Therefore, when performing an EEG-only study, we should be careful when interpreting the resting-state results based on only spatial profiles of BLP features in the sense of BOLD functional networks. In any case, the direct spatio-temporal relation of EEG BLP and RSNs does not seem trivial, and one should keep in mind that we cannot use knowledge from one modality to interpret results from other modality solely on spatial colocalization assumption and without further discussion.

The statistical methods used to test the reproducibility of EEG-fMRI patterns, as well as the relationship between EEG SSPs and BOLD signatures/RSNs, involved multiple levels of construction and processing of descriptive statistics that were only at the end subjected to statistical hypothesis testing. This resulted in an increased abstraction level of the final results. For instance, when testing the spatial colocalization hypothesis in Section 2.11.1, we began with derived EEG spatio-spectral patterns, which are ICs representing band-limited power fluctuations in EEG. We believed that the spatio-spectral signature of a given component represents a spatial pattern (one for each frequency band), and its temporal signature depicts the evolution of this spatial pattern over time. To obtain a BOLD-related map, it is straightforward to use a temporal signature of a given EEG component as an explanatory variable in a typical GLM analysis as is Section 2.8 and obtain BOLD signatures (a spatial map of beta coefficients). In the ideal case of spatial colocalization, the BOLD signature as a spatial map would perfectly spatially correspond to the spatio-spectral signature of a given EEG component. We can compute the level of spatial similarity of those two maps by (in our case) Spearman correlation. We applied a weighted average to address the issue of separate spatial maps for each EEG frequency band. Additionally, given that both types of spatial maps (EEG SSPs and BOLD signatures) exhibit spatial smoothness, assessing the level of statistical significance was challenging. Therefore, we employed permutation statistics to generate a meaningful null distribution as a final step. In summary, by the spatial colocalization, we tested whether the spatial similarity between a given EEG SSP and its BOLD signature (assessed by the weighted average of band-wise spatial map similarity with BOLD signature) was, on average, higher than between the same EEG SSP and BOLD signatures of other components. In other cases, as in Sections 2.9.2, ??, 2.11.3, we utilized a very similar methodology, customizing this multistep approach to suit each specific hypothesis. We aimed to consistently compare comparable derived measures (spatial against spatial and temporal against temporal) and construct a meaningful null distribution for statistical analysis. Of note, the thresholding of the various spatial maps prior to statistical analyses was omitted because those analyses focused on comparing overall pattern similarity between the maps rather than interpreting individual map values. This omission did not introduce bias into the analysis but might introduce some noise levels.

A possible explanation for not observing spatial/spatio-spectral and temporal correspondence between EEG BLP patterns and BOLD ICs might lie in brain activity phenomena that are assumed to give rise to these two signals. While EEG signal is considered to rise from temporally and spatially synchronous neuron activity, specifically excitatory postsynaptic potentials/currents, the BOLD signal generation is dependent on the rather complex vasodilation effect, which is controlled by multiple different pathways as summarized in Drew (2019). It is also a matter of ongoing research which (and how much) cell types in the central nervous system (CNS) contribute to the overall captured BOLD signal. In an optogenetic mouse model stimulation study Vazquez et al. (2018), the authors found that stimulating interneurons leads to a larger increase in blood flow than stimulating pyramidal neurons; thus, they hypothesize that interneurons might be the primary driver of the BOLD response. Other mouse model studies suggest that also astrocyte responses might contribute to the BOLD signal (Tran et al., 2018) and that astrocyte activity might not modulate neuronal activity (Takata et al., 2018). In contrast to human resting-state EEG-fMRI studies, in many animal studies, it was observed that typically gamma band-limited power (Logothetis et al., 2001; Mateo et al., 2017; Winder et al., 2017; Schölvinck et al., 2010) and multiunit spiking activity (Ma et al., 2016; Mateo et al., 2017; Winder et al., 2017) correlates with the BOLD signal. It is also discussed in (Drew, 2019), that observed low correspondence between EEG and fMRI resting-state signals (correlation coefficient typically lower than 0.3) might not be caused by methodological issues but rather by ongoing vasculature dynamics that might be independent of neural activity. This is argued by evoked activity studies (Logothetis et al., 2001; Winder et al., 2017; Lima et al., 2014), which report much higher levels of correlation (even higher than 0.9), thus a high level of explained variance. Findings in those animal studies, as well as other non-neural signals that are superposed on true neural activity and hemodynamic responses, might all contribute to not directly observing spatio-temporal relation between the proposed EEG- and BOLD-derived patterns. These findings and theories follow mostly the assumption that EEG- and BOLD-derived fluctuations are spatially co-located. Nevertheless, since generators of EEG are generally considered EPSCs, and we assume long-range communication between brain areas, the active brain area (observed in BOLD signal) causes observable EEG activity at a different part of the brain. Therefore, perfect spatial colocalization of EEG and fMRI might not likely emerge. In that case, studying causal relationships might further contribute to a better understanding those modalities. It is also important to note that there are other EEG features distinct from the amplitude of band-limited EEG oscillations relating EEG and BOLD signals, such as reported BOLD correlates of infraslow EEG fluctuations (ISF) Hiltunen et al. (2014); a signal that is known to be correlated to amplitude of band-limited EEG oscillations Monto et al. (2008), although the relation of the BOLD correlates of ISF EEG and BLP power oscillations has not yet been explicitly studied. Finally, we must also acknowledge the impact of EEG data quality acquired within the MRI environment. Despite employing state-of-the-art methods for EEG preprocessing, effectively suppressing the two most prominent artifacts (pulse and gradient) remains challenging, partly due to their mutual interaction (Steyrl and Müller-Putz, 2019). One such improvement of the proposed methodology would be to evaluate the quality of EEG recorded during fMRI with the in-scanner EEG recorded prior to fMRI acquisition. Further development of the preprocessing methods is necessary primarily because the mentioned artifacts partially overlap with the EEG frequency bands of interest.

5 Conclusion

We comprehensively compared and integrated EEG whole-brain patterns of neural dynamics with concurrently measured fMRI BOLD data. We derived EEG patterns based on a spatio-spectral decomposition of band-limited EEG power in the source-reconstructed space for that purpose. The decomposition has proven reproducible regarding the similarity of the EEG spatio-spectral signatures with the correlation patterns of the BOLD signal, and we further illustrated that the SSPs obtained from the source-reconstructed space are anatomically interpretable. Furthermore, we observed statistically significant but weak spatial correspondence between EEG spatial profiles and their BOLD signatures, supporting both plausibility in mutual validation, while challenging the view of perfect spatial co-localization of EEG power and BOLD fluctuations. Moreover, we did not observe substantial spatio-temporal correspondence between the EEG and BOLD components. Overall, the introduced data-driven method with enhanced spatial resolution may be helpful for a deeper understanding of the link between EEG and fMRI, namely, human brain networks. Specifically, this approach might identify complex spatio-spectral dynamics of interest not only in rest but also in task studies of the dynamics of brain activity-both in health and in disease states.

From the data analysis perspective, we showed that EEG spatio-spectral decomposition in the source-reconstructed space is a viable dimensionality reduction technique that provides direct information about an anatomical region related to a specific spatio-spectral pattern and, therefore, provides EEG components that are directly spatially comparable to BOLD spatio-temporal patterns. The presented results shed more light on resting-state EEG-fMRI literature when understanding the link between those two modalities. From an application perspective, the extracted EEG spatio-spectral patterns may serve as a spatio-spectral filter for various applied studies (i.e., provide well-informed dimension reduction to time series of 5 interpretable SSP patterns) or as a basis for other methodological EEG-fMRI work studying the link between those modalities.

To further explore the spatio-spectral space by ICA or similar methods, advanced clustering analysis may contribute to potentially detecting subsequent-weaker or more detailed-patterns. Different types of more subject-individual decomposition methods might be beneficial in understanding more specific patterns. A good trade-off between single-subject ICA decomposition and the group ICA method presented in this direction might be an independent vector analysis (IVA) (Lee et al., 2008). There also exist other methods from the trilinear decomposition family, such as Tucker trilinear decomposition methods (Tucker, 1966) or block-term decomposition (Rontogiannis et al., 2021). Those are far less constrained approaches (compared to PARAFAC) which might provide valuable insights into interactions within EEG data. The proposed methodology may also serve for computational modeling purposes as an EEG feature extraction method for whole-brain network computational models (Sanz Leon et al., 2013; Cakan et al., 2021; Schirner et al., 2018).

Data availability statement

The data analyzed in this study is subject to the following licenses/restrictions: The data that support the findings of this study are available from the corresponding author upon reasonable request. Requests to access these datasets should be directed to Jaroslav Hlinka, aGxpbmthQGNzLmNhcy5jeg==.

Ethics statement

The studies involving humans were approved by the Ethical Committee of Masaryk University (Brno, Czech Republic) for the EEG-fMRI dataset and by the Ethical Committee of NIMH (Klecany, Czech Republic) for the out-of-scanner dataset. These studies were conducted in accordance with local legislation and institutional requirements. The participants provided their written informed consent to participate in these studies.

Author contributions

SJ: Data curation, Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing – original draft, Writing – review & editing. VK: Conceptualization, Formal analysis, Investigation, Methodology, Resources, Supervision, Validation, Writing – original draft, Writing – review & editing. DM: Software, Writing – original draft, Writing – review & editing, Methodology. RM: Data curation, Methodology, Writing – original draft, Writing – review & editing, Software. JH: Conceptualization, Funding acquisition, Investigation, Methodology, Project administration, Resources, Supervision, Validation, Writing – original draft, Writing – review & editing.

Funding

The author(s) declare that financial support was received for the research and/or publication of this article. The publication was supported by the ERDF-Project Brain dynamics, No. CZ.02.01.01/00/22_008/0004643, the Czech Technical University Internal Grant Agency - grant number SGS23/120/OHK3/2T/13, the Ministry of Health Czech Republic - DRO 2021 (National Institute of Mental Health - NIMH, IN: 00023752), the Czech Science Foundation project No. 21-32608S, and the Czech Science Foundation project No. 23-07578K.

Acknowledgments

We acknowledge the core facility MAFIL supported by the Czech-BioImaging large RI project (LM2018129 funded by MEYS CR) for their support with obtaining scientific data presented in this study. We would like to express our gratitude to Marco Marino, Gaia Taberna, and Jessica Samogin for their contributions to EEG preprocessing and forward modeling software development.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The author(s) declared that they were an editorial board member of Frontiers, at the time of submission. This had no impact on the peer review process and the final decision.

Generative AI statement

The author(s) declare that no Gen AI was used in the creation of this manuscript.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnins.2025.1549172/full#supplementary-material

References

Abou-Elseoud, A., Starck, T., Remes, J., Nikkinen, J., Tervonen, O., and Kiviniemi, V. (2010). The effect of model order selection in group pica. Hum. Brain Mapp. 31, 1207–1216. doi: 10.1002/hbm.20929

PubMed Abstract | Crossref Full Text | Google Scholar

Abreu, R., Leal, A., and Figueiredo, P. (2018). EEG-informed fMRI: a review of data analysis methods. Front. Hum. Neurosci. 12:29. doi: 10.3389/fnhum.2018.00029

PubMed Abstract | Crossref Full Text | Google Scholar

Abreu, R., Simões, M., and Castelo-Branco, M. (2020). Pushing the limits of EEG: estimation of large-scale functional brain networks and their dynamics validated by simultaneous fMRI. Front. Neurosci. 14:323. doi: 10.3389/fnins.2020.00323

PubMed Abstract | Crossref Full Text | Google Scholar

Asada, H., Fukuda, Y., Tsunoda, S., Yamaguchi, M., and Tonoike, M. (1999). Frontal midline theta rhythms reflect alternative activation of prefrontal cortex and anterior cingulate cortex in humans. Neurosci. Lett. 274, 29–32. doi: 10.1016/S0304-3940(99)00679-5

PubMed Abstract | Crossref Full Text | Google Scholar

Berger, H. (1929). Über das elektroenkephalogramm des menschen. Archiv für psychiatrie und nervenkrankheiten 87, 527–570. doi: 10.1007/BF01797193

Crossref Full Text | Google Scholar

Bourgeois, F., and Lassalle, J.-C. (1971). An extension of the munkres algorithm for the assignment problem to rectangular matrices. Commun. ACM 14, 802–804. doi: 10.1145/362919.362945

Crossref Full Text | Google Scholar

Bridwell, D. A., Rachakonda, S., Silva, R. F., Pearlson, G. D., and Calhoun, V. D. (2018). Spatiospectral decomposition of multi-subject EEG: evaluating blind source separation algorithms on real and realistic simulated data. Brain Topogr. 31, 47–61. doi: 10.1007/s10548-016-0479-1

PubMed Abstract | Crossref Full Text | Google Scholar

Bridwell, D. A., Wu, L., Eichele, T., and Calhoun, V. D. (2013). The spatiospectral characterization of brain networks: fusing concurrent EEG spectra and fMRI maps. Neuroimage 69, 101–111. doi: 10.1016/j.neuroimage.2012.12.024

PubMed Abstract | Crossref Full Text | Google Scholar

Brookes, M. J., Woolrich, M., Luckhoo, H., Price, D., Hale, J. R., Stephenson, M. C., et al. (2011). Investigating the electrophysiological basis of resting state networks using magnetoencephalography. Proc. Nat. Acad. Sci. 108, 16783–16788. doi: 10.1073/pnas.1112685108

PubMed Abstract | Crossref Full Text | Google Scholar

Cakan, C., Jajcay, N., and Obermayer, K. (2021). neurolib: a simulation framework for whole-brain neural mass modeling. Cogn. Comput. 15, 1132–1152. doi: 10.1007/s12559-021-09931-9

Crossref Full Text | Google Scholar

Calhoun, V. D., Stevens, M. C., Pearlson, G. D., and Kiehl, K. A. (2004). fMRI analysis with the general linear model: removal of latency-induced amplitude bias by incorporation of hemodynamic derivative terms. Neuroimage 22, 252–257. doi: 10.1016/j.neuroimage.2003.12.029

PubMed Abstract | Crossref Full Text | Google Scholar

de Munck, J. C., Gonçalves, S. I., Huijboom, L., Kuijer, J. P., Pouwels, P. J., Heethaar, R. M., et al. (2007). The hemodynamic response of the alpha rhythm: an EEG/fMRI study. Neuroimage 35, 1142–1151. doi: 10.1016/j.neuroimage.2007.01.022

PubMed Abstract | Crossref Full Text | Google Scholar

de Munck, J. C., Gonçalves, S. I., Mammoliti, R., Heethaar, R. M., and Da Silva, F. L. (2009). Interactions between different EEG frequency bands and their effect on alpha-fMRI correlations. Neuroimage 47, 69–76. doi: 10.1016/j.neuroimage.2009.04.029

PubMed Abstract | Crossref Full Text | Google Scholar

Delorme, A., and Makeig, S. (2004). EEGlab: an open source toolbox for analysis of single-trial EEG dynamics including independent component analysis. J. Neurosci. Methods 134, 9–21. doi: 10.1016/j.jneumeth.2003.10.009

PubMed Abstract | Crossref Full Text | Google Scholar

Drew, P. J. (2019). Vascular and neural basis of the bold signal. Curr. Opin. Neurobiol. 58, 61–69. doi: 10.1016/j.conb.2019.06.004

PubMed Abstract | Crossref Full Text | Google Scholar

Feige, B., Scheffler, K., Esposito, F., Di Salle, F., Hennig, J., and Seifritz, E. (2005). Cortical and subcortical correlates of electroencephalographic alpha rhythm modulation. J. Neurophysiol. 93, 2864–2872. doi: 10.1152/jn.00721.2004

PubMed Abstract | Crossref Full Text | Google Scholar

Fox, M. D., Snyder, A. Z., Vincent, J. L., Corbetta, M., Van Essen, D. C., and Raichle, M. E. (2005). The human brain is intrinsically organized into dynamic, anticorrelated functional networks. Proc. Nat. Acad. Sci. 102, 9673–9678. doi: 10.1073/pnas.0504136102

PubMed Abstract | Crossref Full Text | Google Scholar

Friston, K. J., Fletcher, P., Josephs, O., Holmes, A., Rugg, M., and Turner, R. (1998). Event-related fMRI: characterizing differential responses. Neuroimage 7, 30–40. doi: 10.1006/nimg.1997.0306

PubMed Abstract | Crossref Full Text | Google Scholar

Glover, G. H., Li, T.-Q., and Ress, D. (2000). Image-based method for retrospective correction of physiological motion effects in fMRI: retroicor. Magn. Reson. Med. 44, 162–167. doi: 10.1002/1522-2594(200007)44:1<162::AID-MRM23>3.0.CO;2-E

PubMed Abstract | Crossref Full Text | Google Scholar

Goldman, R. I., Stern, J. M., Engel, J. Jr, and Cohen, M. S. (2002). Simultaneous EEG and fMRI of the alpha rhythm. Neuroreport 13:2487. doi: 10.1097/00001756-200212200-00022

Crossref Full Text | Google Scholar

Gonçalves, S. I., De Munck, J. C., Pouwels, P. J., Schoonhoven, R., Kuijer, J. P., Maurits, N. M., et al. (2006). Correlating the alpha rhythm to bold using simultaneous EEG/fMRI: inter-subject variability. Neuroimage 30, 203–213. doi: 10.1016/j.neuroimage.2005.09.062

PubMed Abstract | Crossref Full Text | Google Scholar

Hallez, H., Vanrumste, B., Grech, R., Muscat, J., De Clercq, W., Vergult, A., et al. (2007). Review on solving the forward problem in EEG source analysis. J. Neuroeng. Rehabil. 4, 1–29. doi: 10.1186/1743-0003-4-46

PubMed Abstract | Crossref Full Text | Google Scholar

Harshman, R. A., and Lundy, M. E. (1994). Parafac: parallel factor analysis. Comput. Stat. Data Anal. 18, 39–72. doi: 10.1016/0167-9473(94)90132-5

Crossref Full Text | Google Scholar

Hiltunen, T., Kantola, J., Elseoud, A. A., Lepola, P., Suominen, K., Starck, T., et al. (2014). Infra-slow EEG fluctuations are correlated with resting-state network dynamics in fMRI. J. Neurosci. 34, 356–362. doi: 10.1523/JNEUROSCI.0276-13.2014

PubMed Abstract | Crossref Full Text | Google Scholar

Himberg, J., Hyvärinen, A., and Esposito, F. (2004). Validating the independent components of neuroimaging time series via clustering and visualization. Neuroimage 22, 1214–1222. doi: 10.1016/j.neuroimage.2004.03.027

PubMed Abstract | Crossref Full Text | Google Scholar

Hyvärinen, A., and Oja, E. (2000). Independent component analysis: algorithms and applications. Neural Netw. 13, 411–430. doi: 10.1016/S0893-6080(00)00026-5

PubMed Abstract | Crossref Full Text | Google Scholar

Jorge, J., Van der Zwaag, W., and Figueiredo, P. (2014). EEG-fMRI integration for the study of human brain function. Neuroimage 102, 24–34. doi: 10.1016/j.neuroimage.2013.05.114

PubMed Abstract | Crossref Full Text | Google Scholar

Labounek, R., Bridwell, D. A., Mareček, R., Lamoš, M., Mikl, M., Bednařík, P., et al. (2019). EEG spatiospectral patterns and their link to fMRI bold signal via variable hemodynamic response functions. J. Neurosci. Methods 318, 34–46. doi: 10.1016/j.jneumeth.2019.02.012

PubMed Abstract | Crossref Full Text | Google Scholar

Labounek, R., Bridwell, D. A., Mareček, R., Lamoš, M., Mikl, M., Slavíček, T., et al. (2018). Stable scalp EEG spatiospectral patterns across paradigms estimated by group ica. Brain Topogr. 31, 76–89. doi: 10.1007/s10548-017-0585-8

PubMed Abstract | Crossref Full Text | Google Scholar

Labounek, R., Wu, Z., Bridwell, D. A., Brázdil, M., Jan, J., and Nestrašil, I. (2021). Blind visualization of task-related networks from visual oddball simultaneous EEG-fMRI data: spectral or spatiospectral model? Front. Neurol. 12:644874. doi: 10.3389/fneur.2021.644874

PubMed Abstract | Crossref Full Text | Google Scholar

Laufs, H., Holt, J. L., Elfont, R., Krams, M., Paul, J. S., Krakow, K., et al. (2006). Where the bold signal goes when alpha EEG leaves. Neuroimage 31, 1408–1418. doi: 10.1016/j.neuroimage.2006.02.002

PubMed Abstract | Crossref Full Text | Google Scholar

Laufs, H., Kleinschmidt, A., Beyerle, A., Eger, E., Salek-Haddadi, A., Preibisch, C., et al. (2003a). EEG-correlated fMRI of human alpha activity. Neuroimage 19, 1463–1476. doi: 10.1016/S1053-8119(03)00286-6

PubMed Abstract | Crossref Full Text | Google Scholar

Laufs, H., Krakow, K., Sterzer, P., Eger, E., Beyerle, A., Salek-Haddadi, A., et al. (2003b). Electroencephalographic signatures of attentional and cognitive default modes in spontaneous brain activity fluctuations at rest. Proc. Nat. Acad. Sci. 100, 11053–11058. doi: 10.1073/pnas.1831638100

PubMed Abstract | Crossref Full Text | Google Scholar

Lee, J.-H., Lee, T.-W., Jolesz, F. A., and Yoo, S.-S. (2008). Independent vector analysis (iva): multivariate approach for fMRI group study. Neuroimage 40, 86–109. doi: 10.1016/j.neuroimage.2007.11.019

PubMed Abstract | Crossref Full Text | Google Scholar

Li, C., Yuan, H., Shou, G., Cha, Y.-H., Sunderam, S., Besio, W., et al. (2018). Cortical statistical correlation tomography of EEG resting state networks. Front. Neurosci. 12:365. doi: 10.3389/fnins.2018.00365

PubMed Abstract | Crossref Full Text | Google Scholar

Lima, B., Cardoso, M. M., Sirotin, Y. B., and Das, A. (2014). Stimulus-related neuroimaging in task-engaged subjects is best predicted by concurrent spiking. J. Neurosci. 34, 13878–13891. doi: 10.1523/JNEUROSCI.1595-14.2014

PubMed Abstract | Crossref Full Text | Google Scholar

Lindquist, M. A., Loh, J. M., Atlas, L. Y., and Wager, T. D. (2009). Modeling the hemodynamic response function in fMRI: efficiency, bias and mis-modeling. Neuroimage 45, S187–S198. doi: 10.1016/j.neuroimage.2008.10.065

PubMed Abstract | Crossref Full Text | Google Scholar

Liu, Q., Farahibozorg, S., Porcaro, C., Wenderoth, N., and Mantini, D. (2017). Detecting large-scale networks in the human brain using high-density electroencephalography. Hum. Brain Mapp. 38, 4631–4643. doi: 10.1002/hbm.23688

PubMed Abstract | Crossref Full Text | Google Scholar

Liu, Q., Ganzetti, M., Wenderoth, N., and Mantini, D. (2018). Detecting large-scale brain networks using EEG: impact of electrode density, head modeling and source localization. Front. Neuroinform. 12:4. doi: 10.3389/fninf.2018.00004

PubMed Abstract | Crossref Full Text | Google Scholar

Logothetis, N. K., Pauls, J., Augath, M., Trinath, T., and Oeltermann, A. (2001). Neurophysiological investigation of the basis of the fMRI signal. Nature 412, 150–157. doi: 10.1038/35084005

PubMed Abstract | Crossref Full Text | Google Scholar

Ma, Y., Shaik, M. A., Kozberg, M. G., Kim, S. H., Portes, J. P., Timerman, D., et al. (2016). Resting-state hemodynamics are spatiotemporally coupled to synchronized and symmetric neural activity in excitatory neurons. Proc. Nat. Acad. Sci. 113, E8463–E8471. doi: 10.1073/pnas.1525369113

PubMed Abstract | Crossref Full Text | Google Scholar

Magri, C., Schridde, U., Murayama, Y., Panzeri, S., and Logothetis, N. K. (2012). The amplitude and timing of the bold signal reflects the relationship between local field potential power at different frequencies. J. Neurosci. 32, 1395–1407. doi: 10.1523/JNEUROSCI.3985-11.2012

PubMed Abstract | Crossref Full Text | Google Scholar

Makeig, S., Bell, A., Jung, T.-P., and Sejnowski, T. J. (1995). “Independent component analysis of electroencephalographic data,” in Advances in Neural Information Processing Systems, 8.

Google Scholar

Mantini, D., Franciotti, R., Romani, G. L., and Pizzella, V. (2008). Improving meg source localizations: an automated method for complete artifact removal based on independent component analysis. Neuroimage 40, 160–173. doi: 10.1016/j.neuroimage.2007.11.022

PubMed Abstract | Crossref Full Text | Google Scholar

Mantini, D., Perrucci, M. G., Del Gratta, C., Romani, G. L., and Corbetta, M. (2007). Electrophysiological signatures of resting state networks in the human brain. Proc. Nat. Acad. Sci. 104, 13170–13175. doi: 10.1073/pnas.0700668104

PubMed Abstract | Crossref Full Text | Google Scholar

Mareček, R., Lamoš, M., Labounek, R., Bartoň, M., Slavíček, T., Mikl, M., et al. (2017). Multiway array decomposition of EEG spectrum: implications of its stability for the exploration of large-scale brain networks. Neural Comput. 29, 968–989. doi: 10.1162/NECO_a_00933

PubMed Abstract | Crossref Full Text | Google Scholar

Marecek, R., Lamos, M., Mikl, M., Barton, M., Fajkus, J., Brazdil, M., et al. (2016). What can be found in scalp EEG spectrum beyond common frequency bands. EEG-fMRI study. J. Neural Eng. 13:046026. doi: 10.1088/1741-2560/13/4/046026

PubMed Abstract | Crossref Full Text | Google Scholar

Marino, M., Arcara, G., Porcaro, C., and Mantini, D. (2019). Hemodynamic correlates of electrophysiological activity in the default mode network. Front. Neurosci. 13:1060. doi: 10.3389/fnins.2019.01060

PubMed Abstract | Crossref Full Text | Google Scholar

Marino, M., Liu, Q., Brem, S., Wenderoth, N., and Mantini, D. (2016). Automated detection and labeling of high-density EEG electrodes from structural mr images. J. Neural Eng. 13:056003. doi: 10.1088/1741-2560/13/5/056003

PubMed Abstract | Crossref Full Text | Google Scholar

Marino, M., Liu, Q., Koudelka, V., Porcaro, C., Hlinka, J., Wenderoth, N., et al. (2018). Adaptive optimal basis set for BCG artifact removal in simultaneous EEG-fMRI. Sci. Rep. 8, 1–11. doi: 10.1038/s41598-018-27187-6

PubMed Abstract | Crossref Full Text | Google Scholar

Maris, E., and Oostenveld, R. (2007). Nonparametric statistical testing of EEG-and MEG-data. J. Neurosci. Methods 164, 177–190. doi: 10.1016/j.jneumeth.2007.03.024

PubMed Abstract | Crossref Full Text | Google Scholar

Martınez-Montes, E., Valdés-Sosa, P. A., Miwakeichi, F., Goldman, R. I., and Cohen, M. S. (2004). Concurrent EEG/fMRI analysis by multiway partial least squares. Neuroimage 22, 1023–1034. doi: 10.1016/j.neuroimage.2004.03.038

PubMed Abstract | Crossref Full Text | Google Scholar

Mateo, C., Knutsen, P. M., Tsai, P. S., Shih, A. Y., and Kleinfeld, D. (2017). Entrainment of arteriole vasomotor fluctuations by neural activity is a basis of blood-oxygenation-level-dependent “resting-state” connectivity. Neuron 96, 936–948. doi: 10.1016/j.neuron.2017.10.012

PubMed Abstract | Crossref Full Text | Google Scholar

Meyer, M. C., Janssen, R. J., Van Oort, E. S. B., Beckmann, C. F., and Barth, M. (2013). The quest for EEG power band correlation with ica derived fMRI resting state networks. Front. Hum. Neurosci. 7:315. doi: 10.3389/fnhum.2013.00315

PubMed Abstract | Crossref Full Text | Google Scholar

Miwakeichi, F., Martınez-Montes, E., Valdés-Sosa, P. A., Nishiyama, N., Mizuhara, H., and Yamaguchi, Y. (2004). Decomposing EEG data into space-time-frequency components using parallel factor analysis. Neuroimage 22, 1035–1045. doi: 10.1016/j.neuroimage.2004.03.039

PubMed Abstract | Crossref Full Text | Google Scholar

Monto, S., Palva, S., Voipio, J., and Palva, J. M. (2008). Very slow EEG fluctuations predict the dynamics of stimulus detection and oscillation amplitudes in humans. J. Neurosci. 28, 8268–8272. doi: 10.1523/JNEUROSCI.1910-08.2008

PubMed Abstract | Crossref Full Text | Google Scholar

Moosmann, M., Ritter, P., Krastel, I., Brink, A., Thees, S., Blankenburg, F., et al. (2003). Correlates of alpha rhythm in functional magnetic resonance imaging and near infrared spectroscopy. Neuroimage 20, 145–158. doi: 10.1016/S1053-8119(03)00344-6

PubMed Abstract | Crossref Full Text | Google Scholar

Murta, T., Leite, M., Carmichael, D. W., Figueiredo, P., and Lemieux, L. (2015). Electrophysiological correlates of the bold signal for EEG-informed fMRI. Hum. Brain Mapp. 36, 391–414. doi: 10.1002/hbm.22623

PubMed Abstract | Crossref Full Text | Google Scholar

Niazy, R. K., Beckmann, C. F., Iannetti, G. D., Brady, J. M., and Smith, S. M. (2005). Removal of fMRI environment artifacts from EEG data using optimal basis sets. Neuroimage 28, 720–737. doi: 10.1016/j.neuroimage.2005.06.067

PubMed Abstract | Crossref Full Text | Google Scholar

Oostenveld, R., Fries, P., Maris, E., and Schoffelen, J.-M. (2011). Fieldtrip: open source software for advanced analysis of meg, EEG, and invasive electrophysiological data. Comput. Intell. Neurosci. 2011:156869. doi: 10.1155/2011/156869

PubMed Abstract | Crossref Full Text | Google Scholar

Pascual-Marqui, R. D., Lehmann, D., Koukkou, M., Kochi, K., Anderer, P., Saletu, B., et al. (2011). Assessing interactions in the brain with exact low-resolution electromagnetic tomography. Philos. Trans. R. Soc. A 369, 3768–3784. doi: 10.1098/rsta.2011.0081

PubMed Abstract | Crossref Full Text | Google Scholar

Penny, W. D., Friston, K. J., Ashburner, J. T., Kiebel, S. J., and Nichols, T. E. (2011). Statistical Parametric Mapping: The Analysis of Functional Brain Images. New York: Elsevier.

Google Scholar

Piorecky, M., Koudelka, V., Strobl, J., Brunovsky, M., and Krajca, V. (2019). Artifacts in simultaneous hdEEG/fMRI imaging: a nonlinear dimensionality reduction approach. Sensors 19:4454. doi: 10.3390/s19204454

PubMed Abstract | Crossref Full Text | Google Scholar

Prestel, M., Steinfath, T. P., Tremmel, M., Stark, R., and Ott, U. (2018). fMRI bold correlates of EEG independent components: spatial correspondence with the default mode network. Front. Hum. Neurosci. 12:478. doi: 10.3389/fnhum.2018.00478

PubMed Abstract | Crossref Full Text | Google Scholar

Rontogiannis, A. A., Kofidis, E., and Giampouras, P. V. (2021). Block-term tensor decomposition: model selection and computation. IEEE J. Sel. Top. Signal Process. 15, 464–475. doi: 10.1109/JSTSP.2021.3051488

Crossref Full Text | Google Scholar

Sanz Leon, P., Knock, S. A., Woodman, M. M., Domide, L., Mersmann, J., McIntosh, A. R., et al. (2013). The virtual brain: a simulator of primate brain network dynamics. Front. Neuroinform. 7:10. doi: 10.3389/fninf.2013.00010

PubMed Abstract | Crossref Full Text | Google Scholar

Scheeringa, R., Petersson, K. M., Kleinschmidt, A., Jensen, O., and Bastiaansen, M. C. (2012). EEG alpha power modulation of fMRI resting-state connectivity. Brain Connect. 2, 254–264. doi: 10.1089/brain.2012.0088

PubMed Abstract | Crossref Full Text | Google Scholar

Schirner, M., McIntosh, A. R., Jirsa, V., Deco, G., and Ritter, P. (2018). Inferring multi-scale neural mechanisms with brain network modelling. Elife 7:e28927. doi: 10.7554/eLife.28927

PubMed Abstract | Crossref Full Text | Google Scholar

Schölvinck, M. L., Maier, A., Ye, F. Q., Duyn, J. H., and Leopold, D. A. (2010). Neural basis of global resting-state fMRI activity. Proc. Nat. Acad. Sci. 107, 10238–10243. doi: 10.1073/pnas.0913110107

PubMed Abstract | Crossref Full Text | Google Scholar

Shirer, W. R., Ryali, S., Rykhlevskaia, E., Menon, V., and Greicius, M. D. (2012). Decoding subject-driven cognitive states with whole-brain connectivity patterns. Cerebral cortex 22, 158–165. doi: 10.1093/cercor/bhr099

PubMed Abstract | Crossref Full Text | Google Scholar

Sockeel, S., Schwartz, D., Pélégrini-Issac, M., and Benali, H. (2016). Large-scale functional networks identified from resting-state EEG using spatial ica. PLoS ONE 11:e0146845. doi: 10.1371/journal.pone.0146845

PubMed Abstract | Crossref Full Text | Google Scholar

Steyrl, D., and Müller-Putz, G. R. (2019). Artifacts in EEG of simultaneous EEG-fMRI: pulse artifact remainders in the gradient artifact template are a source of artifact residuals after average artifact subtraction. J. Neural Eng. 16:016011. doi: 10.1088/1741-2552/aaec42

PubMed Abstract | Crossref Full Text | Google Scholar

Taberna, G. A., Samogin, J., and Mantini, D. (2021). Automated head tissue modelling based on structural magnetic resonance images for electroencephalographic source reconstruction. Neuroinformatics 19, 585–596. doi: 10.1007/s12021-020-09504-5

PubMed Abstract | Crossref Full Text | Google Scholar

Takata, N., Sugiura, Y., Yoshida, K., Koizumi, M., Hiroshi, N., Honda, K., et al. (2018). Optogenetic astrocyte activation evokes bold fMRI response with oxygen consumption without neuronal activity modulation. Glia 66, 2013–2023. doi: 10.1002/glia.23454

PubMed Abstract | Crossref Full Text | Google Scholar

Tran, C. H. T., Peringod, G., and Gordon, G. R. (2018). Astrocytes integrate behavioral state and vascular signals during functional hyperemia. Neuron 100, 1133–1148. doi: 10.1016/j.neuron.2018.09.045

PubMed Abstract | Crossref Full Text | Google Scholar

Tucker, L. R. (1966). Some mathematical notes on three-mode factor analysis. Psychometrika 31, 279–311. doi: 10.1007/BF02289464

PubMed Abstract | Crossref Full Text | Google Scholar

Tyvaert, L., LeVan, P., Grova, C., Dubeau, F., and Gotman, J. (2008). Effects of fluctuating physiological rhythms during prolonged EEG-fMRI studies. Clin. Neurophysiol. 119, 2762–2774. doi: 10.1016/j.clinph.2008.07.284

PubMed Abstract | Crossref Full Text | Google Scholar

Uusitalo, M. A., and Ilmoniemi, R. J. (1997). Signal-space projection method for separating MEG or EEG into components. Med. Biol. Eng. Comput. 35, 135–140. doi: 10.1007/BF02534144

PubMed Abstract | Crossref Full Text | Google Scholar

Vazquez, A. L., Fukuda, M., and Kim, S.-G. (2018). Inhibitory neuron activity contributions to hemodynamic responses and metabolic load examined using an inhibitory optogenetic mouse model. Cerebral Cortex 28, 4105–4119. doi: 10.1093/cercor/bhy225

PubMed Abstract | Crossref Full Text | Google Scholar

Vega, M. L. B., Michel, C. M., Saxena, S., White, T., and Valdes-Sosa, P. A. (2022). Neuroimaging and global health. Neuroimage 260:119458. doi: 10.1016/j.neuroimage.2022.119458

PubMed Abstract | Crossref Full Text | Google Scholar

Vejmelka, M., Pokorná, L., Hlinka, J., Hartman, D., Jajcay, N., and Paluš, M. (2015). Non-random correlation structures and dimensionality reduction in multivariate climate data. Clim. Dyn. 44, 2663–2682. doi: 10.1007/s00382-014-2244-z

Crossref Full Text | Google Scholar

Vorderwülbecke, B. J., Carboni, M., Tourbier, S., Brunet, D., Seeber, M., Spinelli, L., et al. (2020). High-density electric source imaging of interictal epileptic discharges: how many electrodes and which time point? Clin. Neurophysiol. 131, 2795–2803. doi: 10.1016/j.clinph.2020.09.018

PubMed Abstract | Crossref Full Text | Google Scholar

Vorwerk, J., Oostenveld, R., Piastra, M. C., Magyari, L., and Wolters, C. H. (2018). The fieldtrip-simbio pipeline for EEG forward solutions. Biomed. Eng. Online 17, 1–17. doi: 10.1186/s12938-018-0463-y

PubMed Abstract | Crossref Full Text | Google Scholar

Winder, A. T., Echagarruga, C., Zhang, Q., and Drew, P. J. (2017). Weak correlations between hemodynamic signals and ongoing neural activity during the resting state. Nat. Neurosci. 20, 1761–1769. doi: 10.1038/s41593-017-0007-y

PubMed Abstract | Crossref Full Text | Google Scholar

Xia, M., Wang, J., and He, Y. (2013). Brainnet viewer: a network visualization tool for human brain connectomics. PLoS ONE 8:e68910. doi: 10.1371/journal.pone.0068910

PubMed Abstract | Crossref Full Text | Google Scholar

Yuan, H., Ding, L., Zhu, M., Zotev, V., Phillips, R., and Bodurka, J. (2016). Reconstructing large-scale brain resting-state networks from high-resolution EEG: spatial and temporal comparisons with fMRI. Brain Connect. 6, 122–135. doi: 10.1089/brain.2014.0336

PubMed Abstract | Crossref Full Text | Google Scholar