Concurrent EEG- and fMRI-derived functional connectomes exhibit linked dynamics

Large-scale functional connectivity of the human brain, commonly observed using functional Magnetic Resonance Imaging (fMRI), exhibits a whole-brain spatial organization termed the functional connectome. The fMRI-derived connectome shows dynamic reconfigurations that are behaviorally relevant. Due to the indirect nature of fMRI, it is unclear whether such topographic changes reliably reflect modulation in neuronal connectivity patterns. Here, we directly compared concurrent fMRI-derived and electrophysiological connectivity dynamics on a connection-wise basis across the whole connectome. Dynamic whole-brain functional connectivity (dFC) was assessed during resting-state in two independent concurrent fMRI-electroencephalography (EEG) datasets (42 subjects total) using a sliding window approach. FMRI- and EEG-derived dFC shared significant mutual information in all canonical EEG frequency bands. Notably, this was true for virtually all connections. Across all EEG frequency bands, connections with the strongest link between EEG and fMRI dynamics tied the default mode network (DMN) to the rest of the brain. Beyond this frequency-independent multimodal dFC, fMRI connectivity covaried with EEG connectivity in a frequency-specific manner in two distributed sets of connections for delta and gamma bands, respectively. These results generalized across the two datasets. Our findings promote the DMN as a universal hub of dynamics across frequencies, but also show that spatial distribution of fMRI and EEG dFC differ across the canonical EEG-frequency bands. This study reveals a close relationship between time-varying changes in whole-brain connectivity patterns of electrophysiological and hemodynamic signals. The results support the value of EEG for studying the whole-brain connectome and provide evidence for a neuronal basis of fMRI-derived dFC.


Introduction
To date, our knowledge about the topography of neural communication in the human brain is largely derived from fMRI, an indirect measure of neural activity. Since the discovery that brain activity is correlated between distant brain regions in resting-state functional magnetic resonance imaging (rs-fMRI) (1) extensive efforts have been undertaken towards establishing the neural origin of this functional connectivity (FC) (2)(3)(4)(5)(6)(7)(8). Using fMRI, the brain has been characterized in terms of different intrinsic connectivity networks (ICNs) (9)(10)(11). Collectively, the connectivity within and between ICNs can be represented as a whole-brain connectivity graph, or connectome (12). This line of research has traditionally focused on static properties of FC by averaging brain activity across the entire recording period (13). Though dominated by fMRI-based research, static whole-brain FC has been shown to have a spatial organization comparable across EEG and fMRI (2,5) and MEG and fMRI (4,6) and intracranial EEG.

Static connectivity relationship across EEG and fMRI
We used resting-state concurrent EEG-fMRI data from two independent datasets (primary n=26, generalization n=16). Preprocessed fMRI signal was averaged to the 68 cortical regions of the Desikan atlas (39,40). Preprocessed EEG signals band-passed to each canonical frequency band (δ,θ,α,β,γ) were source-reconstructed (41)(42)(43) to the same atlas (Fig. 1a). First, we sought to confirm that the previously reported static relation between fMRI-derived and EEG-derived connection-wise connectivity strength (3,8) holds true for the two datasets of this study. To this end, we assessed connectivity averaged across the total duration of EEG and fMRI data. Connectivity in source-projected EEG was quantified as band-limited phase synchrony using imaginary coherence, and connectivity in fMRI was quantified as Pearson's r for the same whole-brain parcellation atlas. In line with our prior work (8), static FC averaged of all EEG bands were correlated to FC of fMRI (main dataset fMRI vs. δ/θ/α/β/γ: r=0.34/0.34/0.33/0.36/0.29; generalization dataset fMRI vs. δ/θ/α/β/γ: r=0.34/0.33/0.36/0.41/0.39, taking the average FC across all subjects). This observation reaffirms the link between the spatial organization of connectivity across modalities and confirms sufficient quality of EEG source localization to the whole-brain parcellation. Fig. 1: a) Construction of EEG and fMRI connectomes. EEG was source reconstructed and fMRI signal was averaged for the 68 cortical regions of the Desikan atlas (39). Pearson's correlation of fMRI timecourses and imaginary coherence of bandlimited EEG source signal were used to build connectomes. b) Dynamic FC was derived from a 1min window sliding at 2s (= repetition time of fMRI). Mutual information between fMRI-dFC and EEG-dFC was statistically compared to a null model consisting of the mutual information between fMRI-dFC and phase scrambled EEG-dFC (8,44).

Dynamic connectivity relationship across EEG and fMRI
Next, we used a sliding window approach to assess whether time-varying changes in the whole-brain connectivity pattern derived from fMRI are linked to dynamics in band-limited phase synchrony across modalities (Fig. 1b). We used the information-theoretic measure of mutual information, which assesses the relationship between the modalities without assuming linearity. We tested these codynamics against a null model that spatially randomizes connections of the EEG connectome while preserving comparable modular organization (8,44). Dynamic FC in δ, θ, α, β and γEEG showed high mutual information with fMRI connectome dynamics significantly outperforming the null model in virtually all region-pairs: in 100% of connections for fMRI vs. δEEG-fMRI, θEEG-fMRI, αEEG-fRMI, and γEEG-fMRI, and in 99.91% of connections for βEEG-fMRI (Bonferroni corrected for 2248 connections, p<2.19*10 -5 ). This strong cross-modal relationship was confirmed in the generalization dataset, although it was significant in a smaller number of connections especially for βEEG-fMRI and γEEG-fMRI: fMRI vs. δEEG-fMRI=98.68%, θEEG-fMRI=91.83% αEEG-fRMI=87.97%, βEEG-fMRI=54.48%, γEEG-fMRI=31.12% (Bonferroni corrected for 2248 connections, p<2.19*10 -5 ). The relative reduction in effect size in the replication dataset is in line with the smaller sample size and shorter recording duration. Fig. 2 show the distribution of mutual information of one randomly selected subject per dataset.  Additionally, we tested for generalizability at a connection-wise level. The connection-wise strength of mutual information averaged across all subjects was strongly correlated across primary and replication datasets for fMRI compared to δ/θ/αEEG (r=0.48/0.46/0.34) although no correlation was observed for fMRI vs. β/γEEG (r=0.04/-0.02, Fig. 3). A split-half approach indicated that the lack of connection-wise generalization for β and γ bands is due to the lower signal-to-noise ratio compared with the slower frequencies (SI Results and Table S1). Fig. 3: Connection-wise comparison of the EEG-fMRI relationship across the two datasets. Each data point reflects the mutual information between EEG-dFC and fMRI-dFC for a given connection averaged across all subjects of the respective datasets. Scatter plots are provided for the relation of fMRI-dFC to a) δEEG, b) θEEG, c) αEEG, d) βEEG, and e) γEEG-dFC. The whole-brain distribution of mutual information correlated across the two datasets for δ, θ and αEEG.
In summary, fMRI signal correlations were linked to slow modulations of oscillatory phase synchrony across the vast majority of the connectome's region-pairs. This was true for all canonical EEG frequency bands across both datasets. The strength of the EEG-fMRI relationship was correlated on a connectionwise basis between primary and generalization datasets in the δ, θ, and α bands.

Spatial topography of the dynamic relationship
Next, we sought to characterize the spatial topography of co-dynamics beyond the above-described all-encompassing relation between EEG and fMRI dFC. Specifically, we assessed how the connections for which fMRI dFC was most strongly linked to EEG dFC were distributed over canonical neurocognitive ICNs (11) and across EEG frequency bands. To this end, the contrast of EEG-fMRI mutual information in real against null data was subjected to network based statistics (NBS) (45) separately for each EEG frequency band. The initial element-wise t-threshold was chosen so that NBS would identify the connected set of 200 connections with strongest mutual information between EEG and fMRI dynamics. The connection-wise threshold and the ensuing NBS cluster-corrected significance corresponding to the top-200 cutoff are listed in Table 1. For both datasets the top 200 connections are consistently drawn from connections that have high mutual information (see SI Results).  We observed no significant difference between datasets in this permutation analysis (n=100,000 iterations, p>0.0018 for all EEG frequency bands, Bonferroni corrected for multiple comparisons at 28 network pairs). Thus, while we observed connection-wise generalization for δ/θ/α but not β/γ bands (Fig. 3), the results in Fig. 4 support generalization for all bands at the coarser ICN-wise resolution.  Interestingly, for all EEG frequencies, this distributed network of strongest cross-modal dynamics aligned with the Default Mode Network (DMN). We established that this distribution of the top 200 connections and the ensuing DMN-dominance were not driven by the number of ICN nodes or other potential biases. For each network pair (e.g. DMN-VIS) we tested whether the number of connections was significantly higher than chance by randomly selecting (n=100,000) 200 connections from the main dataset (Tables S2-S6 In summary, while our first analysis showed a widespread significant relationship between fMRI and EEG for virtually all connections and irrespective of oscillation frequency, the connections of the strongest (top 200) cross-modal relationship showed a dominant role of the DMN.

Frequency-specificity of the dynamic relationship
Finally, we sought to directly and statistically corroborate the frequency specificity of the cross-modal relationship on a connection-wise basis. To this end, we combined the EEG-fMRI mutual information matrices for all EEG bands into an ANOVA (5 levels for 5 frequency bands). Statistical testing indicated a main effect of EEG band (F-test, p=0.041, NBS corrected). Exploratory post-hoc t-tests revealed that both delta-and gamma-band dFC are organized in a frequency specific network of increased mutual information to fMRI dynamics relative to all other frequencies (p<0.05, NBS-corrected, Table 2). Fig. 5a visualizes this frequency-specific network for delta and gamma bands in physical brain space, and Fig.  5b maps the connections according to canonical ICNs (we report the connected set of top 100 connections). Confirming the observation described above (Fig. 4b), both networks of frequency specific cross-modal dynamics consisted predominantly of DMN connections to the rest of the brain; the delta-specific set showed a strong preference for Limbic-DMN connections whereas the gamma specific set shared a preference between DMN-Ventral-attention and DMN-Limbic connections (Fig.  5c).  In the generalization dataset the F-test and NBS analysis likewise showed a significant main effect of frequency (F<0.0002). Exploratory post-hoc analysis replicated a set of connections with significantly higher mutual information for δEEG-fMRI and γEEG-fMRI, and additionally for βEEG-fMRI (Table 2). To conclude, when contrasted directly across frequency bands, frequency-specificity of the EEG-fMRI relationship was confirmed with a DMN dominance.

Discussion
Dynamics of intrinsic brain activity are necessary for human behavior, which is inherently flexible. While these dynamics are most commonly studied using fMRI, the relationship between BOLD dynamics and their underlying neurophysiological basis is still unclear. In this study, we demonstrated that dFC derived from fMRI shares mutual information with EEG-derived dFC in all canonical frequency bands. Importantly, this applies to the whole connectome. Further, the relative contribution of each EEG frequency band to BOLD connectivity dynamics varies across space. These results shed new light on the relationship between electrophysiological and BOLD connectivity dynamics and in particular provide strong evidence that BOLD-derived dFC is directly linked to neurophysiological dynamics with corresponding spatial organization.

Static connectivity
In line with prior work, the spatial pattern of static EEG and fMRI connectomes were significantly correlated. Only two prior studies have investigated whole-brain connectomes in concurrent EEG and fMRI (3,8), and the current study confirms this important observation. However, the two prior studies were restricted to static relationships, while the true advantage of concurrent as opposed to separately recorded EEG and fMRI connectivity lies in its potential to reveal dynamics co-fluctuating across modalities.

Dynamic connectivity
Dynamic FC has been shown to exist both on slow hemodynamic (14) and fast electrophysiological timescales (31,32). Regarding the relationship between EEG and fMRI, Chang et al. (33) and Allen et al. (35) showed that BOLD FC dynamics are linked to band-limited EEG global field power. The current study advances to spatially resolved EEG connectomes, demonstrating the link between slow dynamics of fast electrophysiological connectivity and the slow BOLD connectivity dynamics.
This finding has implications for the interpretation of fMRI-derived dFC. Since physiological and nonphysiological noise heavily contribute to fMRI-derived dFC, the degree to which such dynamics reflect changes in neural communication is difficult to assess; choosing an appropriate null model can be challenging (15,22). Instead of comparing fMRI dFC to a null model, we chose to compare fMRI dynamics to direct measures of neural dynamics such as MEG or EEG. Our approach of concurrently assessing dFC in EEG provides evidence that veridic dynamics can indeed be derived from fMRI, as they are directly related to slow changes in the underlying electrophysiological connectivity. Another important advantage of our approach is the independence of outcomes from both EEG-and fMRIrelated artefacts (e.g. eye movements in EEG or magnetic field inhomogeneities in fMRI), since spurious dynamics due to random noise in each modality will cancel each other out in the joint model we utilized. Comparing a global EEG parameter to local fMRI connections is an approach prone to spurious cross-modal links stimming from a connection-unspecific global shift (e.g. breathing (46)). In contrast, our approach comparing both modalities on a connection-wise level is unlikely to be impacted by averaged global patterns.
Our results also have important implications for neurophysiological connectivity dynamics with respect to methodological reliability and observed timescales. Due to concerns regarding the ill-posed nature of EEG source localization, connectivity approaches in whole-brain parcellation space are underused (30). The close connection-wise relationship to fMRI-derived dFC provides strong support for the relevance of source-localized EEG to the study of the whole-brain functional connectome. With respect to timescales, prior MEG-based whole-brain investigations have established fast dynamics in interregional connectivity at ~50-100ms (31,32). Connectivity dynamics at these fast timescales in EEG have been shown to correlate with the slow changes observed in fMRI, albeit in EEG sensor space rather that reconstructed brain parcellations (47,48). Concurrent EEG-fMRI studies investigating neurophysiological dynamics at infraslow speeds (typically defined as < 0.1Hz (49,50)) have been either limited to a coarse global field average of connectivity across EEG electrode pairs (27,28), or have focused solely on signal fluctuations (Hiltunen et al. (49)). Extending beyond these important studies, we show that connectivity derived from fast oscillation phase synchrony in EEG exhibits meaningful fluctuations in the infraslow range.

Spatial distribution of the cross-modal dynamics
We found the strongest relationship between EEG and fMRI dynamics in cortico-DMN connections. The dominance of DMN interactions with other ICNs has been previously demonstrated for dynamic connectivity derived from MEG (51). Similarly, Vidaurre et al. (32) observed increased MEG withinnetwork coherence in the DMN as opposed to visual and sensorimotor areas. More generally, the DMN has been proposed to form a central hub that integrates multisensory input from different brain regions (52). As such, it is no surprise that we observed connections between DMN and the visual and motor systems to show some of the strongest relations between EEG and fMRI.
Interestingly, connections with the strongest interrelation between EEG and fMRI connectivity predominantly spanned the canonical intrinsic networks. While we found that the vast majority of connections were dynamically linked between EEG and fMRI, the cross-modal relation was weaker within as compared to across networks with the exception of DMN-DMN connectivity. This finding may seem counterintuitive since the investigation of amplitude fluctuations (as opposed to connectivity) shows a strong link between EEG and fMRI within ICNs (53,54), suggesting a strong static withinnetwork connectivity. Similarly, de Pasquale et al. (51) observed strong within-network connectivity dynamics for αand βbands alongside between-network connectivity dynamics. One possible explanation for our results is that within-network connections are less dynamic than connections between different intrinsic networks. This is in line with the observation that DMN regions are among the most dynamic, as measured by fMRI-derived dFC (15).

Frequency specificity of the cross-modal dynamics
We observed both a delta dominance in DMN-limbic connections and γ-dominance in VA-DMN connections (Fig. 4). The relationship between static MEG and fMRI connectivity has been shown to vary across connections depending on neurophysiological frequency (5). Likewise, spatially distinct patterns of frequency-specific connectivity dynamics have been observed in MEG (32). As such, Vidaurre et al. (32) have suggested a division of the DMN into a low-frequency and a high-frequency subnetwork. Although the network subdivisions used by Vidaurre et al. (32) do not directly map to the spatial distribution of the canonical ICN networks that we chose, our results for delta vs. γ-are in line with a low and high frequency-specific contribution of connections encompassing the DMN. This suggests a differential organization of frequencies in the DMN, which are synced between EEG and fMRI dynamics.

Methodological considerations and limitations
The relatively low spatial resolution (number of regions) of the selected Desikan atlas was imposed by the number of 64 EEG electrodes. While it has been shown that a parcellation adapted to the actual EEG montage improves the quality of the source reconstruction (55), it is unclear if the fMRI signal would suffer from a parcellation scheme imposed by the EEG montage. A future approach could be extending the MR-compatible EEG setup to 128 or 256 electrodes (56) to gain data quality comparable to MEG recordings (57). Regarding temporal resolution, the sliding window approach has been previously criticized for assuming slowly changing dynamics as opposed to fast instantaneous switches (see (21) for review). While the timescale of observable dynamics depends on the choice of window length, a comparative analysis of the optimal windowing parameters to derive dFC was beyond the scope of this study. Importantly, we show that results from the chosen parameter set generalize across two independent datasets.
Our focus on common patterns generalizing across datasets resulted in a stringent statistical threshold (cf. Table 1) that comes at the cost of being insensitive to small effects. This was primarily a result of the limited power of the generalization dataset (smaller sample size and shorter recordings). In particular, β-and γ-connectivity seem to be influenced by the low signal-to-noise ratio of the generalization dataset (with 45.52% and 68.88% of βEEG-and γEEG-derived connections sharing no significant dynamics with fMRI, as opposed to virtually no insignificant connections in the primary dataset). Beyond the signal-to-noise limitations in human concurrent EEG-fMRI, there is evidence that γ-and β-band connectivity may contain information complementary to fMRI-derived connectivity. A weaker relation to fMRI connectivity has been reported for β and low-γ compared to other bands in intracranial electrophysiological recordings in humans (38,58) and in animals (59,60). The weaker relationship for β-and low γ-bands likely reflects a general property of the electrophysiology-fMRI relationship. This interpretation is also in line with our previous finding that γ-connectivity shows a unique relationship to structural connectivity not shared by fMRI-derived connectivity (8). Importantly, at the coarser ICN-wise (Fig. 4) as compared to connection-wise resolution (Fig. 3), the results generalized across all bands including β and γ. The generalization of effects in this study is especially supportive of the robustness of the EEG-fMRI dFC relationship in light of substantial differences across the two datasets (3T vs. 1.5T MRI field strength, eyes-closed vs. eyes-open resting-state and differences in fMRI sequences and subject demographics).

Conclusion
We observed a link between electrophysiological and fMRI-derived dynamic functional connectivity which demonstrates, on a connection-wise level across the whole brain, that fMRI-derived connectivity entails slow dynamics of fast electrophysiological connectivity. While this link exists across all canonical electrophysiological frequency bands, the strength of the cross-modal relationship varies over connections in a frequency-specific manner, especially for δ-and γ-bands. In conclusion, this study provides strong multimodal evidence of slow time varying connectivity dynamics of intrinsic brain activity.

Methods
We analyzed a primary dataset (n=26) and tested for generalizability across a different sample by using an independent dataset (n=16) from a different site. The generalization dataset is openly available at https://osf.io/94c5t/ and described in detail in Deligianni et al. (3,61).

Primary Dataset Subjects
We recruited 26 healthy subjects (8 females, mean age 24.39, age range 18-31) with no history of neurological or psychiatric illness. Ethical approval has been obtained from the local Research Ethics Committee (CPP Ile de France III) and informed consent has been obtained from all subjects.

Data Acquisition
We acquired three runs of 10 minutes eyes-closed resting-state in one concurrent EEG-fMRI session (Tim-Trio 3T, Siemens The acquisition was part of a study with two additional naturalistic film stimulus of 10 minutes not analyzed in the current study, and acquired after runs 1 and 2 of the resting state as described in Morillon et al. (63). The three runs resulted in a total length of 30 minutes of resting-state fMRI per subject. Subjects wore earplugs to attenuate scanner noise and were asked to stay awake, avoid movement and close their eyes during resting-state recordings. In three subjects, one of three rest sessions each was excluded due to insufficient EEG quality.

EEG
EEG was corrected for the gradient artefact induced by the scanner using the template subtraction and adaptive noise cancelation followed by lowpass filtering at 75Hz, downsampling to 250Hz (65) and cardiobalistic artefact template subtraction (66) using EEGlab v.7 (http://sccn.ucsd.edu/eeglab) and the FMRIB plug-in (https://fsl.fmrib.ox.ac.uk/eeglab/fmribplugin/). Data then was analyzed with Brainstorm software (43), which is documented and freely available under the GNU general public license (http://neuroimage.usc.edu/brainstorm, version 10 th August 2017). Bandpass-filtering was carried out at 0.3-70 Hz. Data was segmented according to TR of the fMRI acquisition (2s epochs). Epochs containing head motion artifacts in EEG were visually identified after semi-automatically preselecting epochs where signal in any channel exceeded the mean channel timecourse by 4 std. These segments were excluded from the analysis. Electrode positions and T1 were coregistered by manually moving the electrode positions onto the electrode artifacts visible in the T1 image. Using the OpenMEEG BEM model, a forward model of the skull was calculated based on the individual T1 image of each subject (42,67).

Joint motion scrubbing
For all analyses, both fMRI volumes and EEG epochs were excluded for time periods where motion was identified in either modality. Time periods with motion were defined as volumes exceeding the framewise displacement threshold FD=0.5 in fMRI (69), and by visual inspection in EEG as described above. Additionally, for sliding window connectivity (see section Sliding window connectivity below), windows with more than 10% of their datapoints (>3 fMRI volumes or >3 EEG epochs) removed by this motion scrubbing procedure were excluded from dynamic connectivity analysis. The joint motion scrubbing approach resulted in a mean of 544 out of 870 sliding windows (range 262-813) for the main dataset and 216 out of 272 sliding windows (range 112-259) for the generalization dataset.

Static connectivity
Static connectivity was estimated for fMRI data by calculating Pearson's correlation of the BOLD timecourse between each region pair over the duration of each run and averaged across the 3 runs. For EEG, the connectivity (imaginary coherence) calculated for each 2s epoch was averaged across all runs (Fig. 1a).

Sliding window connectivity
dFC matrices were calculated using a rectangular sliding window of 1 min (using imaginary coherence for EEG and Pearson's correlation for fMRI, this resulted in 30 datapoints per window). The window length was chosen as a tradeoff between maximizing the number of datapoints without discarding relevant dynamic BOLD frequencies (19,64,70) while also taking into account the theoretical limitations of shorter window lengths to reliably detect dFC (24,71). Most importantly we show that the chosen parameter set reliably replicated across two independent datasets.  (72,73) was calculated for the resulting EEG and fMRI dFC matrices, building a new connectivity matrix of joint EEG-fMRI, based on mutual information strength between the modalities (for each EEG frequency band). In contrast to linear measures such as correlation, mutual information is an information theoretic measure which is able to also capture cross-modal relationships in connectivity dynamics without assuming linearity or Gaussian constraints (74). Mutual information has previously been shown to be helpful when combining EEG and fMRI (75).

Null model
For each connection, mutual information was then compared to a null model of EEG-fMRI mutual information using spatially phase-randomized EEG matrices (8,44). In brief, we applied the approach of phase-randomization proposed by Prichard and Theiler (76) to whole-brain connectomes by extracting the Eigenvectors of the EEG connectome at each sliding window. Subsequently each Eigenvector was Fourier transformed and the phases of this transformation were then randomly shifted. The result of the phase-shift then was back-transformed using the inverse Fourier transform, and the phase shifted Eigenvectors were used to reverse the Eigen-decomposition. This approach has been shown to generate pseudo-connectomes uncorrelated to the original matrix while keeping a similar spatial structure (such as interhemispheric connections and ordering in ICN-like networks) and constant global mean comparable to the original connectome (44) (Fig. 1b).
An explicit test for dynamics is not needed as spurious dynamics are implicitly excluded by combining the two modalities: in the scenario that the previously established static link across the modalities (8) is driving a window-wise constant connectivity in both modalities the entropy of the cross-modal relation -dominated by random noise -would be low, resulting in spurious mutual information close to zero. Spurious dynamics in one modality resulting e.g. from physiological artifacts or sampling error are unlikely to co-occur in the other modality. Exceptions, such as head motion that may affect both modalities at the same time, likely have a spatially distributed impact on many connections of the connectome, and are controlled for by our null model that would preserve such "mass" connectivity changes. The randomization process was carried out 50 times for the EEG connectome, and mutual information to the unaltered fMRI connectome was calculated for each iteration. For final statistical comparison of this null model to the original EEG-fMRI mutual information, we calculated an average mutual information matrix from the 50 iterations for each subject.
To assess the connected set with the strongest EEG-fMRI relationship in each EEG frequency band, we additionally subjected the connection-wise tests against the null model to Network Based Statistics (NBS https://sites.google.com/site/bctnet/comparison/nbs, Version 1.2, correcting for multiple comparisons; (45)). NBS controls the family-wise error rate of the mass-univariate testing at every connection. This method is a non-parametric cluster-based approach to finding connected sets of nodes that significantly differ across thresholded connectivity matrices.

Frequency specific analysis
To test for frequency specificity, we included frequency-specific mutual information matrices of all bands (δEEG-fMRI, θEEG-fMRI, etc.) and for all subjects in an ANOVA (frequency band as 1 factor with 5 levels), while discarding the pseudo-matrices generated for the first statistical analysis. Using NBS an F-Test was carried out to determine if the EEG-bands contributed differentially to the mutual information with fMRI-derived dFC. Posthoc t-tests between one band vs. all the other bands were carried out to explore if any EEG band expressed a network of stronger mutual information than observed in the other bands.

Intrinsic network analysis
To further interpret outcomes in the context of neurocognitive networks, we mapped the extracted 200 connections to 7 canonical intrinsic networks (Visual, Somato-Motor, Default Mode, Fronto-Parietal, Ventral Attention (largely corresponding to Cingulo-opercular (Dosenbach et al. 2006)), Dorsal Attention and Limbic) as described in Yeo et al. (11). The number of connections falling into each network pair were counted (e.g. DMN to Visual). To assure that the observed connectivity pattern did not arise from random sampling into the different networks, we also created 100,000 random networks of 200 connections to derive the probability that a connection randomly falls into one of the network pairs.

Generalization dataset Subjects
This dataset comprises 17 healthy adults. Ethical approval has been obtained from the UCL Research Ethics Committee (project ID:4290/001) and informed consent has been obtained from all subjects. One subject was excluded as T1 data quality was not sufficient to run the Freesurfer recon-all command, resulting in a final group of 16 subjects (6 females, mean age: 32.41, range 22-53).

Data Acquisition
We

Data processing
The fMRI data was processed as described for the primary dataset with the exception that no slicetime correction was carried out (in accordance with the original processing in Deligianni (3)). EEG was corrected for the gradient artefact using the template subtraction and adaptive noise cancelation followed by a downsampling to 250Hz and cardiobalistic artefact template subtraction using the Brain Vision Analyzer 2 software (Brain Products, Gilching, Germany). Due to apparent low frequency drift artefacts in several subjects, EEG data was high pass filtered at 0.05Hz instead of the 0.03Hz used in the primary dataset. Because of the differing TR the sliding window for 1 minute was now consisting of 28 volumes (28*2.16s = 60.48s). EEG data processing was equivalent to the primary dataset, with the epochs being 2.16s instead of 2s to match the fMRI TR.
All following analysis steps were identical to the primary dataset.

Mutual information strength
The top-200 connections (Fig. 4)  We further confirmed that the MI magnitude difference between the top 200 connections and the rest of the connections is consistent across the two datasets. We performed an analogous test between the two datasets by masking the generalization data with the top-200 connections of the main dataset. MI was generally increased as compared to the rest of the connections for MI between fMRI and δ/θ/α but not increased for MI between fMRI and β/γ (one-sided ttest, fMRI vs. δ/θ/α/β/γ: p=5.00*10 -7 /3.71*10 -8 /0.012/0.053/0.428, p<0.05).
When assessing generalization at a connection-wise resolution, we observed correlation of connection-wise EEG-fMRI mutual information for δ, θ, and α bands but not for β and γEEG (Fig. 3). To assess whether this discrepancy was due to lower signal to noise ratio in the higher frequencies, we investigated whether replication within each dataset would show a similar pattern. We divided the primary dataset (respectively the generalization dataset) into two groups of 8 subjects each (taking only the first 16 subjects of the primary dataset). Indeed, we observed a within dataset correlation of fMRI vs. δ/θ/α connection-wise mutual information, whereas no correlation between fMRI vs. β/γ was observed (Table S1).

Mutual information strength of frequency specific networks
Within the subnetworks selected by NBS (Fig. 5), we further tested for connection-wise correlation of mutual information strength between the two datasets. When masked by the top-100 connections of the 'δ > other bands' contrast of the main dataset (Table 2), the correlation between mutual information (averaged across all subjects) across the two datasets was: r=0.56 (p=2.0*10 -9 ), a value qualitatively higher than for the whole brain (r=0.48). Contrarily, masking the top-100 connections of the 'γ > other bands' contrast r=-0.0785 was not significant (p=0.44, in line with the missing relationship between connection-wise MI correlation in gamma; Fig. 3).

Generalization of MI distribution across canonical intrinsic networks
We established that this distribution of the top 200 connections and the ensuing DMN-dominance were not driven by the number of ICN nodes or other potential biases. For each network pair (e.g. DMN-VIS) we tested whether the number of connections was significantly higher than chance by randomly selecting (n=100,000 iterations) 200 connections from the main dataset (Table S2-S6).