Hierarchical decoupling of electromagnetic and haemodynamic cortical networks

| Whole-brain neural communication is typically estimated from statistical associations among electromagnetic or haemodynamic time-series. The relationship between functional network architectures recovered from these two types of neural activity remains unknown. Here we map electromagnetic networks (measured using magnetoencephalography; MEG) to haemodynamic networks (measured using functional magnetic resonance imaging; fMRI). We ﬁnd that the relationship between the two modalities is regionally heterogeneous and systematically follows the cortical hierarchy, with close correspondence in unimodal cortex and poor correspondence in transmodal cortex. Comparison with the BigBrain histological atlas reveals that electromagnetic-haemodynamic coupling is driven by laminar differentiation and neuron density, suggesting that the mapping between the two modalities can be explained by cytoarchitectural variation. Importantly, haemodynamic connectivity cannot be explained by electromagnetic activity in a single frequency band, but rather arises from the mixing of multiple neurophysiological rhythms. Correspondence between the two is largely driven by slower rhythms, particularly the beta (15-29 Hz) frequency band. Collectively, these ﬁndings demonstrate highly organized but only partly overlapping patterns of connectivity in MEG and fMRI functional networks, opening fundamentally new avenues for studying the relationship between cortical micro-architecture and multi-modal connectivity patterns.


INTRODUCTION
The structural wiring of the brain imparts a distinct signature on neuronal co-activation patterns.Inter-regional projections promote signaling and synchrony among distant neuronal populations, giving rise to coherent neural dynamics, measured as regional time series of electromagnetic or hemodynamic neural activity [40].Systematic co-activation among pairs of regions can be used to map functional connectivity networks.Over the past decade, these dynamics are increasingly recorded without task instruction or stimulation; the resulting "intrinsic" functional connectivity is thought to reflect spontaneous neural activity.
The macro-scale functional architecture of the brain is commonly inferred from electromagnetic or haemodynamic activity.The former can be measured using electroencephalography (EEG) or magnetoencephalography (MEG), while the latter is measured using functional magnetic resonance imaging (fMRI).Numerous studies -using both MEG and fMRI -have reported evidence of intrinsic functional patterns that are highly organized [5,10,14,17,30,90,119], reproducible [15,26,47,85] and comparable to task-driven co-activation patterns [15,27,101].
How do electromagnetic and haemodynamic networks relate to one another?Although both modalities attempt to capture the same underlying biological process (neural activity), they are sensitive to different physiological mechanisms and ultimately reflect neural activity at fundamentally different time scales [4,52,56,92,93].Emerging theories emphasize a hierarchy of time scales of intrinsic fluctuations across the cortex [43,83,91,98], where unimodal cortex is more sensitive to immediate changes in the sensory environment, while transmodal cortex is more sensitive to prior context [6,23,24,57,64,67].This raises the possibility that the alignment between the relatively slower functional architecture captured by fMRI and faster functional architecture captured by MEG may systematically vary across the cortex.
Previous reports have found some, but not complete, global overlap between the two modalities.Multiple MEG and fMRI independent components -representing spatiotemporal signatures of resting-state intrinsic networks -show similar spatial topography, particularly the visual, somatomotor and default mode components [5,14,17,62].Moreover, fMRI and MEG/EEG yield comparable fingerprinting accuracy, suggesting that they encode common information [29,33,39,95].Finally, global edge-wise comparisons between fMRI networks and electrocorticography (ECoG) [12], EEG [32,117,118] and MEG [44,63,106] also yield moderate correlations.These studies make two tacit assumptions.First, by focusing on global comparisons, they assumed that electromagnetic and haemodynamic connectivity are related in exactly the same way across the brain, precluding the possibility that the coupling between fMRI and MEG connectivity varies from region to region, i.e. is regionally heterogeneous.Second, previous studies typically compared haemodynamic and electromagnetic networks for each band separately, assuming that individual rhythms contribute independently to haemodynamic connectivity.This effectively precludes the possibility that superposition and mixing of elementary electromagnetic rhythms manifests as patterns of haemodynamic connectivity [74,106].Whether and how the two modalities can be mapped to one another at the network level remains unknown.
Here we directly study the link between MEG-and fMRI-based functional networks.We use a linear multifactor model that allows to represent the haemodynamic functional connectivity profile of a given brain region as a linear combination of its electromagnetic functional connectivity in multiple frequency bands.We then explore how the two modalities align across the neocortex.

RESULTS
Data were derived using task-free MEG and fMRI recordings in the same unrelated participants from the Human Connectome Project (HCP [108]; n = 33).We first develop a simple regression-based model to map regional MEG connectivity to regional fMRI connectivity using group-average data.We then investigate how re-gionally heterogeneous the correspondence between the two is, and how different rhythms contribute to this regional heterogeneity.Finally, we conduct extensive sensitivity testing to demonstrate that the results are robust to multiple methodological choices.

Relating haemodynamic and electromagnetic connectivity
To relate fMRI and MEG functional connectivity patterns, we apply a multi-linear regression model [110] (Fig. 1).The model is specified for each brain region separately, attempting to predict a region's haemodynamic connectivity profile from its electromagnetic connectivity Figure 2. Regional model fit | (a) Spatial organization of fMRI-MEG correspondence is depicted for the regional model fit (95% interval).The cross-modal correspondence of connectivity profiles of brain regions is distributed heterogeneously across the cortex, representing regions with low or high correspondence.Strong cross-modal correspondence is observed in sensory areas whereas poor correspondence is observed for higher order regions.(b) Spatial organization of the cross-modal correspondence is compared with the functional hierarchical organization of cerebral cortex [75].The two are significantly anti-correlated, confirming poor fMRI-MEG correspondence in connectivity profile of higher-order, transmodal areas compared to strong correspondence for sensory, unimodal regions.(c) Regions are stratified by their affiliation with macro-scale intrinsic networks [119].The distribution of R 2 is depicted for each network, displaying a systematic gradient of cross-modal correspondence with the highest correspondence in the visual network and lowest correspondence in the default mode network.(d) The model fit is related to the cytoarchitectural variation of the cortex, estimated from the cell staining intensity profiles at various cortical depths obtained from the BigBrain histological atlas [2,88].Bigger circles denote statistically significant associations after correction for multiple comparisons by controlling the false discovery rate (FDR) at 5% alpha [11].The peak association between cross-modal correspondence and cytoarchitecture is observed approximately at cortical layer IV that has high density of granule cells.Staining intensity profiles are depicted across the cortex for the most pial, the middle and the white matter surfaces.rs denotes Spearman rank correlation.Intrinsic networks: vis = visual; sm = somatomotor; da = dorsal attention; va = ventral attention; lim = limbic; fp = frontoparietal; dmn = default mode.
profile.The dependent variable is a row of the fMRI functional connectivity (FC) matrix and the independent variables are the corresponding rows of MEG FC matrices for six canonical electrophysiological bands, estimated using amplitude envelope correlation (AEC [18]) with spatial leakage correction (see Methods for more details).For a model fitted for a given node i, the observations in the model are the connections of node i to the other j = i regions (Fig. 1a).The model predicts the fMRI FC profile of node i (i.e.i-th row) from a linear combination of MEG FC profiles of node i in the six frequency bands (i.e.i-th rows of MEG FC matrices).Collectively, the model embodies the idea that multiple rhythms could be superimposed to give rise to regionally heterogeneous haemodynamic connectivity.Indeed, we find that the relationship between haemodynamic and electromagnetic connectivity is highly heterogeneous.Band-limited MEG connectivity matrices are moderately correlated with fMRI connectivity, ranging from r = −0.06 to r = 0.36 (Fig. 1b; r denotes Pear-son correlation coefficient).The regional multi-linear model fits range from adjusted-R 2 = −0.002 to adjusted-R 2 = 0.72 (R 2 denotes coefficient of determination; hereafter we refer to adjusted-R 2 as R 2 ), suggesting a close correspondence in some regions and poor correspondence in others (Fig. 1c,e; see Fig. S1 for bandspecific regional model fits).For comparison, a single global model is fitted to the data, predicting the entire upper triangle of the fMRI FC matrix from a linear combination of the upper triangles of six MEG FC matrices.The global model, which simultaneously relates whole-brain fMRI FC to the whole-brain MEG FC, yields an R 2 = 0.15 (Fig. 1d,e).Importantly, the global model clearly obscures the wide range of correspondences, which can be considerably greater or smaller for individual regions.

Hierarchical organization of cross-modal connectivity
We next consider the spatial organization of fMRI-MEG correspondence.Fig. 2a shows the spatial distribution of Figure 3. Dominance analysis | Dominance analysis is performed for each regional multi-linear model to quantify how MEG connectivity at different rhythms contribute to regional patterns of cross-modal correspondence [3,19].(a) The overall contribution of each frequency band is depicted for the regional model (box plots).Slower bands contribute the most to model fit whereas faster bands (low and high gamma) contribute the least.(b) The mean contribution of different rhythms is estimated for the intrinsic networks.Consistent with the overall contributions depicted in panel (a), the greatest contribution is associated with slower frequency bands, specifically the beta band.(c) The most dominant predictor (frequency band) is depicted for each brain region, confirming overall higher contributions from slower rhythms across the cortex.(d) Frequency band contribution to the regional cross-modal correspondence is shown separately for different rhythms across the cortex (95% intervals).
regional R 2 values, representing regions with low or high correspondence.Regions with strong cross-modal correspondence include the visual, somato-motor and auditory cortex.Regions with low cross-modal correspondence include the posterior cingulate, lateral temporal and medial prefrontal cortex.
Collectively, the spatial layout of cross-modal correspondence bears a resemblance to the unimodaltransmodal cortical hierarchy observed in large-scale functional and microstructural organization of the cerebral cortex [67].To assess this hypothesis, we first used the neuromaps toolbox [76] to compare the crossmodal R 2 map with the principal functional gradient [70,75] (Fig. 2b).The two are significantly anticorrelated (Spearman rank correlation coefficient r s = −0.69,p spin = 0.0001), suggesting strong cross-modal correspondence in unimodal sensory cortex and poor correspondence in transmodal cortex.We then stratify regions by their affiliation with macro-scale intrinsic networks and computed the mean R 2 in each network [119] (Fig. 2c).Here we also observe a systematic gradient of cross-modal correspondence, with the strongest correspondence in the visual network and poorest correspondence in the default mode network.
Finally, we relate the cross-modal R 2 map to the cytoarchitectural variation of the cortex (Fig. 2d).We use the BigBrain histological atlas to estimate granular cell density at multiple cortical depths [2,88].Cell-staining intensity profiles were sampled across 50 equivolumetric surfaces from the pial surface to the white matter surface to estimate laminar variation in neuronal density and soma size.Fig. 2d shows the correlation between MEG-fMRI correspondence and cell density (y-axis) at different cortical depths (x-axis).Interestingly, the model fit is associated with cytoarchitectural variation of the cortex, with the peak association observed approximately at cortical layer IV that has high density of granular cells and separates supra-and infra-granular layers [86,87,113].

Heterogeneous contributions of multiple rhythms
How do different rhythms contribute to regional patterns of cross-modal correspondence?To address this question, we perform a dominance analysis for every regional multi-linear model [3,19].The technique estimates the relative importance of predictors by constructing all possible combinations of predictors and quantifying the relative contribution of each predictor as additional variance explained (i.e. gain in adjusted-R 2 ) by adding that predictor to the models.Fig. 3a shows the global dominance of each frequency band.Overall, we observe the greatest contributions from slower bands (delta to beta) and smallest contributions from faster bands (low and high gamma).
Zooming in on individual regions and intrinsic networks, we find that the dominance pattern is also re-Figure 4. Sensitivity analysis | (a) A regional cross-validation was performed by pseudorandomly splitting the connectivity profile of a given region into train and test sets based on spatial separation, such that 75% of the closest regions to a random region are selected as the train set and the remaining 25% of the regions are selected as test set (399 repetitions) [54].The multi-linear model is then fitted on the train set and is used to predict the connection strength of the test set for each region and each split.The mean regional model performance across splits is depicted for train and test sets, displaying consistent results between the two (scatter plot).The out-of-sample model performance is stronger in the sensory, unimodal areas compared to transmodal areas, consistent with original findings (Fig. 2).(b) A subject-level cross-validation was performed using a leave-one-out approach.The regional multi-linear model is trained using data from n − 1 subjects and is tested on the held-out subject for each region separately.The mean regional model performance is shown for train and test sets, displaying consistent results between the two (scatter plot).The out-of-sample model performance is stronger in the sensory, unimodal areas compared to transmodal areas, consistent with original findings (Fig. 2).The analysis is also repeated for various processing choices: (c) after regressing out interregional Euclidean distance from connectivity matrices, (d) using MEG connectivity data without spatial leakage correction, (e) using another MEG source reconstruction method (standardized low resolution brain electromagnetic tomography; sLoreta [89]), (f) using a phase-based MEG connectivity measure (phase-locking value; PLV [69,82]), and (g) at a low resolution parcellation (Schaefer-200 atlas [96]).The results are consistent across all control analyses, identifying similar cross-modal correspondence maps as the original analysis (Fig. 2a).All brain maps are shown at 95% intervals.rs denotes Spearman rank correlation.gionally heterogeneous.Namely, the make-up and contribution of specific MEG frequencies to a region's fMRI connectivity profile varies from region to region.Fig. 3b shows the dominance of specific rhythms in each intrinsic network.Fig. 3c shows the most dominant predictor for every brain region.Fig. 3d shows the dominance of specific rhythms separately for each region across cortex.Overall, we observe that slow rhythm connectivity patterns, specifically from beta band, contribute the most to model prediction across the cortex.

Sensitivity analysis
Finally, we note that the present report goes through several decision points that have equally-justified alternatives.Here we explore the other possible choices.First, rather than framing the report from an explanatory perspective (focusing on model fit), we instead derive an equivalent set of results using a predictive perspective (focusing on out-of-sample prediction).We perform cross-validation at both the region-and subjectlevel (Fig. 4a,b).For region-level cross-validation, we pseudorandomly split the connectivity profile of a given region into train and test sets based on spatial separation (inter-regional Euclidean distance), such that 75% of the closest regions to a random region are selected as the train set and the remaining 25% of the regions are selected as test set (399 repetitions; see Methods for more details) [54].We then train the multi-linear model using the train set and predict the connection strength of the test set for each region and each split.The mean regional model performance across splits is consistent for train and test sets (Fig. 4a; r = 0.78, p spin = 0.0001).For subject-level cross-validation, we use leave-one-out-cross validation, wherein we train the regional multi-linear models using data from n − 1 subjects and test each one on the held-out subject.The mean regional model performance is consistent for train and test sets (Fig. 4b; r = 0.90, p spin = 0.0001).Altogether, both analyses give similar, highly concordant results with the simpler model fit-based analysis, identifying strong cross-modal correspondence in unimodal sensory regions and poor correspondence in transmodal areas.
To consider the effect of spatial proximity on the findings, we remove the exponential inter-regional Euclidean distance trend from all connectivity matrices before fitting any model.The results are consistent with and without distance correction (Fig. 4c; correlation with functional hierarchy: r s = −0.53,p spin = 0.0001; correlation with original R 2 : r s = 0.67, p spin = 0.0001).We also obtain consistent findings when we repeat the analysis without accounting for spatial leakage effect in estimating MEG connectivity with AEC (Fig. 4d; correlation with functional hierarchy: r s = −0.60,p spin = 0.0001; correlation with original R 2 : r s = 0.84, p spin = 0.0001).Next, we use another source reconstruction method (standardized low resolution brain electromagnetic tomography; sLoreta [89]) instead of LCMV beamformers, as previous reports suggest that sLoreta improves source localization accuracy [59,60].We then estimate MEG connectivity with AEC and repeat the multi-linear model analysis, identifying similar results as before (Fig. 4e; correlation with functional hierarchy: r s = −0.80,p spin = 0.0001; correlation with original R 2 : r s = 0.85, p spin = 0.0002).Next, we compute MEG connectivity using an alternative, phase-based connectivity measure (phase locking value; PLV [69,82]), rather than the AEC.The two FC measures yield similar cross-modal correspondence maps (Fig. 4f; correlation with functional hierarchy: r s = −0.53,p spin = 0.0022; correlation with original R 2 : r s = 0.66, p spin = 0.0001).We also repeat the analysis using a low resolution parcellation (Schaefer-200 atlas [96]) to ensure that the findings are independent from the choice of parcellation.As before, the results demonstrate similar cross-modal correspondence map (Fig. 4g; correlation with functional hierarchy: r s = −0.70,p spin = 0.0001).Finally, to assess the extent to which the results are influenced by MEG source localization error, we compare the cross-modal correspondence pattern to peak localization error estimated using cross-talk function (CTF) [58-60, 71, 81].No significant association is observed between R 2 pattern and CTF for LCMV (Fig. S2a; r s = −0.14, p spin = 0.6) and sLoreta (Fig. S2b; r s = −0.04,p spin = 0.9) source reconstruction solutions.

DISCUSSION
In the present report we map electromagnetic functional networks to haemodynamic functional networks in the human brain.We find two principal results.First, the relationship between the two modalities is regionally heterogeneous but systematic, reflecting the unimodaltransmodal cortical hierarchy and cytoarchitectural variation.Second, haemodynamic connectivity cannot be explained by electromagnetic connectivity in a single band, but rather reflects mixing and superposition of multiple rhythms.
The fact that the association between the two modalities follows a gradient from unimodal to transmodal cortex resonates with emerging work on cortical hierarchies [67,75,79].Indeed, similar spatial variations are observed for multiple micro-architectural features, such as gene expression [21,42,54], T1w/T2w ratio [66], laminar differentiation [113] and neurotransmitter receptor profiles [41,48,55].Collectively, these studies point to a natural axis of cortical organization that encompasses variations in both structure and function across micro-, meso-and macro-scopic spatial scales.
This spatial variation in structure is thought to manifest as spatial variation in intrinsic dynamics [43,83,98,111,116].In particular, multiple studies have reported faster timescales of intrinsic neural activity in unimodal sensory cortex, and slower timescales in transmodal cortex [43,83,91].This hierarchy of timescales is thought to support a hierarchy of temporal receptive windows: time windows in which a newly arriving stimulus will modify processing of previously presented, contextual information [6,23,24,57,64,67].As a result, unimodal cortex needs to adapt to rapid, uncertain changes in sensory input, while transmodal cortex registers slower prediction error signals, resulting in greater sensitivity to contextual information [34].
Interestingly, we find the closest correspondence between fMRI and MEG functional connectivity in unimodal cortex (including the visual and somatomotor networks) and the poorest correspondence in transmodal cortex (default mode, limbic, fronto-parietal and ventral attention networks).In other words, the functional architectures of the two modalities are consistent early in the cortical hierarchy, presumably reflecting activity related to instantaneous changes in the external environment.Conversely, as we move up the hierarchy, there is a gradual separation between the two architectures, suggesting that they are differently modulated by endogenous inputs and contextual information.How the two types of functional connectivity are related to ongoing task demand is an exciting question for future research.
Why is there systematic divergence between the two modalities?One possibility is that neurovascular coupling is heterogeneous across the cortex [7,13,46,104].Patterns of laminar differentiation -which define the unimodal-transmodal hierarchy [49,61,113] -result in regionally heterogeneous haemodynamic response functions, such that the same stimulus may elicit bloodoxygen-level-dependent (BOLD) responses with different amplitudes and latencies across the cortex [36,37,53,65,80,105].As a result, patterns of haemodynamic connectivity may systematically decouple from the underlying electromagnetic connectivity patterns.
A related possibility is that greater coupling in unimodal cortex is driven by the underlying cytoarchitecture, reflecting greater neuron density, stronger recurrent subcortical inputs and specialized local processing in sensory regions compared to transmodal cortex [28,115].Indeed, we find that the correspondence between haemodynamic and electromagnetic connectivity is associated with cytoarchitectural variation across the cortex, such that regions with higher density of granular cells have higher cross-modal correspondence, and vice versa.This is consistent with the notion that neural oscillations and the BOLD response mirror cytoarchitectural organization, as a result of distinct feedforward and feedback projections between different cortical layers and subcortical regions [20,35,51,97].For example, previous studies of animal electrophysiological recordings demonstrated that visual and frontal cortex gamma activity can be localized to superficial cortical layers (supragranular layers I-III and granular layer IV), whereas alpha and beta activity are localized to deep, infragranular layers (layers V-IV) [8,9,20,72,73,100].Similar findings have been reported in humans using EEG and laminar-resolved BOLD recordings, demonstrating that gamma and beta band EEG power are associated with superficial and deep layer BOLD response, respectively, whereas alpha band EEG power is associated with BOLD response in both superficial and deep layers [97].Another study using human MEG recordings found that source activity differs when deep and superficial surface models are used for source reconstruction analysis [107].
Throughout the present report, we find that fMRI networks are best explained as arising from the superposition of multiple band-limited MEG networks.Although previous work has focused on directly correlating fMRI with MEG/EEG networks in specific bands, we show that synchronized oscillations in multiple bands could potentially combine to give rise to the well studied fMRI functional networks.Indeed, and as expected, the correlation between any individual band-specific MEG network and fMRI is substantially smaller than the multi-linear model that takes into account all bands simultaneously.Previous work on cross-frequency interactions [38] and on multi-layer MEG network organization [16] has sought to characterize the participation of individual brain regions within and between multiple frequency networks.Our findings build on this literature, showing that the superimposed representation may additionally help to unlock the link between MEG and fMRI networks.
It is noteworthy that the greatest contributions to the link between the two modalities came from slower oscillations, such as beta band.Multiple authors have reported that -since it captures slow haemodynamic coactivation -fMRI network connectivity would be mainly driven by slower rhythms [17,32,38,74,92].Here we empirically confirm that although all frequency bands contribute to the emergence of fMRI networks, the greatest contributions come from relatively slower frequencies and the lowest contributions come from relatively faster frequencies.
The present results raise two important questions for future work.First, how does structural connectivity shape fMRI and MEG functional networks [22,102,118]?Given that they capture distinct neurophysiological mechanisms, it is possible that haemodynamic and electromagnetic architectures have a different relationship with structural connectivity and this could potentially explain why they systematically diverge through the cortical hierarchy [110].Second, the present results show how the two modalities are related in a taskfree resting state, but what is the relationship between fMRI and MEG connectivity during cognitive tasks [68]?Given that the two modalities become less correlated in transmodal cortex in the resting state, the relationship between them during task may depend on demand and cognitive functions required to complete the task.
Finally, the present results should be interpreted in light of several methodological considerations.First, although we conduct extensive sensitivity testing, including multiple ways of defining functional connectivity, there exist many more ways in the literature to estimate both fMRI and MEG connectivity [84,112].Second, to ensure that the analyses were performed in the same par-ticipants using both resting state fMRI and MEG, and that the participants have no familial relationships, we utilized a reduced version of the HCP sample.Third, in order to directly compare the contributions of multiple frequency bands, all were entered into the same model.As a result however, the observations in the linear models (network edges) are not independent, violating a basic assumption of these statistical models.For this reason, we only use model fits and dominance values to compare the correspondence of fMRI and MEG across a set of nodes, each of which is estimated under the same conditions.
Despite complementary strengths to image spatiotemporal brain dynamics, the links between MEG and fMRI are not fully understood and the two fields have diverged.The present report bridges the two disciplines by comprehensively mapping haemodynamic and electromagnetic network architectures.By considering the contributions of the canonical frequency bands simultaneously, we show that the superposition and mixing of MEG neurophysiological rhythms manifests as highly structured patterns of fMRI functional connectivity.Systematic convergence and divergence among the two modalities in different brain regions opens fundamentally new questions about the relationship between cortical hierarchies and multi-modal functional networks.

Dataset: Human Connectome Project (HCP)
Resting state magnetoencephalography (MEG) data of a set of healthy young adults (n = 33; age range 22-35 years) with no familial relationships were obtained from Human Connectome Project (HCP; S900 release [108]).The data includes resting state scans of about 6 minutes long (sampling rate = 2034.5Hz; anti-aliasing filter lowpass filter at 400 Hz) and noise recordings for all participants.MEG anatomical data and 3T structural magnetic resonance imaging (MRI) data of all participants were also obtained for MEG pre-processing.Finally, we obtained functional MRI data of the same n = 33 individuals from HCP dataset.All four resting state fMRI scans (two scans with R/L and L/R phase encoding directions on day 1 and day 2, each about 15 minutes long; T R = 720 ms) were available for all participants.

Resting state magnetoencephalography (MEG)
Resting state MEG data was analyzed using Brainstorm software, which is documented and freely available for download online under the GNU general public license ( [103]; http://neuroimage.usc.edu/brainstorm).The MEG recordings were registered to the structural MRI scan of each individual using the anatomical transformation matrix provided by HCP for co-registration, follow-ing the procedure described in Brainstorm's online tutorials for the HCP dataset (https://neuroimage.usc.edu/brainstorm/Tutorials/HCP-MEG).The pre-processing performed by applying notch filters at 60, 120, 180, 240, and 300 Hz, and was followed by a highpass filter at 0.3 Hz to remove slow-wave and DCoffset artifacts.Bad channels were marked based on the information obtained through the data management platform of HCP for MEG data (ConnectomeDB; https: //db.humanconnectome.org/).The artifacts (including heartbeats, eye blinks, saccades, muscle movements, and noisy segments) were then removed from the recordings using automatic procedures as proposed by Brainstorm.More specifically, electrocardiogram (ECG) and electrooculogram (EOG) recordings were used to detect heartbeats and blinks, respectively.We then used Signal-Space Projections (SSP) to automatically remove the detected artifacts.We also used SSP to remove saccades and muscle activity as low-frequency (1-7 Hz) and highfrequency (40-240 Hz) components, respectively.The pre-processed sensor-level data was then used to obtain a source estimation on HCP's fsLR4k cortex surface for each participant.Head models were computed using overlapping spheres and the data and noise covariance matrices were estimated from the resting state MEG and noise recordings.Linearly constrained minimum variance (LCMV) beamformers method from Brainstorm was then used to obtain the source activity for each participant.We performed data covariance regularization to reduce the effect of variable source depth, where the estimated source variance was normalized by the data covariance matrix.Source orientations were constrained to be normal to the cortical surface at each of the 8,000 vertex locations on the fsLR4k surface.Source-level timeseries were then parcellated into 400 regions using the Schaefer-400 atlas [96], such that a given parcel's time series was estimated as the first principal component of its constituting sources' time series.Parcellated timeseries were then used to estimate functional connectivity with an amplitude-based connectivity measure from Brainstorm (amplitude envelope correlation; AEC [18]).An orthogonalization process was applied to correct for the spatial leakage effect by removing all shared zero-lag signals [25].AEC functional connectivity were derived for each participant at six canonical electrophysiological bands (i.e., delta (δ: 2-4 Hz), theta (θ: 5-7 Hz), alpha (α: 8-12 Hz), beta (β: 15-29 Hz), low gamma (lo-γ: 30-59 Hz), and high gamma (hi-γ: 60-90Hz)).Group-average MEG functional connectivity matrices were constructed as the mean functional connectivity across all individuals for each frequency band.For comparison, band-limited group-average AEC matrices were also estimated without correcting for spatial leakage effect.
We also processed the MEG data using additional methodological choices.First, the LCMV source reconstructed and parcellated time-series were used to estimate functional connectivity with an alternative, phasebased connectivity measure (phase locking value; PLV [69,82]) for each frequency band.Second, another source reconstruction method (standardized low resolution brain electromagnetic tomography; sLoreta [89]) was used instead of LCMV beamformers to obtain sourcelevel time-series, given that previous reports suggest that sLoreta improves source localization accuracy [59,60].Source-level time-series, obtained by sLoreta, were then parcellated into 400 regions and were used to estimate AEC matrices with spatial leakage correction for the six frequency bands.Third, to ensure that the findings are independent from choice of parcellation, a low resolution atlas ) was used to parcellate the original LCMV time-series to 200 cortical regions and obtain spatial leakage corrected AEC connectivity matrices.Finally, we estimated MEG source localization errors for LCMV and sLoreta source reconstruction solutions using cross-talk functions (CTF) [58-60, 71, 81].CTF of a given source i is a measure of how activity from all other sources contributes to the activity estimated for the i-th source.Following guidelines from Brainstrom [103] and MNE-Python software packages [50], we used CTF to calculate peak localization error of a given source i as the Euclidean distance between the peak location estimated for source i and the true source location i on the surface model [59,81].

Resting state functional MRI
The functional MRI data were pre-processed using HCP minimal pre-processing pipelines [45,108].Detailed information regarding data acquisition and preprocessing is available elsewhere [45,108].Briefly, all 3T functional MRI time-series were corrected for gradient nonlinearity, head motion using a rigid body transformation, and geometric distortions using scan pairs with opposite phase encoding directions (R/L, L/R) [31].Further pre-processing steps include co-registration of the corrected images to the T1w structural MR images, brain extraction, normalization of whole brain intensity, high-pass filtering (> 2000s FWHM; to correct for scanner drifts), and removing additional noise using the ICA-FIX process [31,94].The pre-processed timeseries were then parcellated into 400 cortical areas using Schaefer-400 parcellation [96].The parcellated timeseries were used to construct functional connectivity matrices as Pearson correlation coefficients between pairs of regional time-series for each of the four scans and each participant.A group-average functional connectivity matrix was constructed as the mean functional connectivity across all individuals and scans.

BigBrain histological data
To characterize the cytoarchitectural variation across the cortex, cell-staining intensity profile data were obtained from the BigBrain atlas [2,88].The BigBrain is a high-resolution (20 µm) histological atlas of a post mortem human brain and includes cell-staining intensities that are sampled at each vertex across 50 equivolumetric surfaces from the pial to the white matter surface using the Merker staining technique [2,78].The staining intensity profile data represent neuronal density and soma size at varying cortical depths, capturing the regional differentiation of cytoarchitecture [2,87,88,113,114].Intensity profiles at various cortical depths can be used to approximately identify boundaries of cortical layers that separate supragranular (cortical layers I-III), granular (cortical layer IV), and infragranular (cortical layers V-VI) layers [88,113,114].The data were obtained on fsaverage surface (164k vertices) from the BigBrainWarp toolbox [88] and were parcellated into 400 cortical regions using the Schaefer-400 atlas [96].

Multi-linear model
A multiple linear regression model was used to assess regional associations between haemodynamic (fMRI) and electromagnetic (MEG) functional connectivity (Fig. 1 [110]).A separate multi-linear model is applied for each brain region from the parcellated data, predicting the region's fMRI functional connectivity profile from its band-limited MEG functional connectivity.The dependent variable is a row of the fMRI connectivity matrix and the independent variables (predictors) are the corresponding rows of MEG connectivity for the six canonical electrophysiological bands.The linear regression model for each brain region i is constructed as following: where the dependant variable F C i is the set of fMRI connections of node i to the other j = i regions and the predictors are sets of MEG connections of node i to the other j = i regions for the six frequency bands The regression coefficients b 1 , ..., b 6 and the intercept b 0 are then optimized to yield maximum correlation between empirical and predicted fMRI connectivity for each brain region.Goodness of fit for each regional model is quantified using adjusted-R 2 (coefficient of determination).

Region-level cross-validation
Region-level cross-validation was performed to assess out-of-sample model performance.Given the spatial autocorrelation inherent to the data, random splits of brain regions into train and test sets may result in out-ofsample correlations that are inflated due to spatial proximity [77].To take this into account, we used a distancedependant cross-validation approach where we pseudorandomly split the connectivity profile of a given region (e.g.node i) into train and test sets based on spatial separation [54].We used inter-regional Euclidean distance to select 75% of the closest regions to a randomly selected source region as the train set and the remaining 25% of the regions as test set.The random source region can be any of the 399 regions connected to node i; hence, the connectivity profile of node i is split into 399 unique train and test sets.We then train the multi-linear model using the train set and predict functional connectivity of the test set for each region and each split.Finally, the model performance is quantified using Pearson correlation coefficient between empirical and predicted values.The cross-validated regional model performance is then estimated as the mean correlation coefficient between empirical and predicted values across splits for each brain region.

Subject-level cross-validation
Leave-one-out cross-validation was performed to assess model performance on held-out subjects.Briefly, the regional multi-linear model is trained using the groupaverage data from n − 1 subjects.The trained model is then used to predict fMRI connectivity profile of each region on the held-out subject (test set).The model performance is quantified as the Pearson correlation coefficient between empirical and predicted connectivity of each region.The analysis is repeated for all subjects and the regional model performance is averaged across individuals.

Dominance analysis
Dominance Analysis was used to quantify the distinct contributions of resting state MEG connectivity at different frequency bands to the prediction of resting state fMRI connectivity in the multi-linear model [3,19] (https://github.com/dominance-analysis/dominance-analysis).Dominance analysis estimates the relative importance of predictors by constructing all possible combinations of predictors and re-fitting the multilinear model for each combination (a model with p predictors will have 2 p − 1 models for all possible combinations of predictors).The relative contribution of each predictor is then quantified as increase in variance explained by adding that predictor to the models (i.e. gain in adjusted-R 2 ).Here we first constructed a multiple linear regression model for each region with MEG connectivity profile of that region at six frequency bands as independent variables (predictors) and fMRI connectivity of the region as the dependent variable to quantify the distinct contribution of each factor using dominance analysis.The relative importance of each factor is estimated as "percent relative importance", which is a summary measure that quantifies the percent value of the additional contribution of that predictor to all subset models.

Null model
To make inferences about the topographic correlations between any two brain maps, we implement a null model that systematically disrupts the relationship between two topographic maps but preserves their spatial autocorrelation [1,77].We used the Schaefer-400 atlas in the HCP's fsLR32k grayordinate space [96,108].The spherical projection of the fsLR32k surface was used to define spatial coordinates for each parcel by selecting the vertex closest to the center-of-mass of each parcel [99,109,110].The resulting spatial coordinates were used to generate null models by applying randomly-sampled rotations and reassigning node values based on the closest resulting parcel (10,000 repetitions).The rotation was applied to one hemisphere and then mirrored to the other hemisphere.
ology, 106(3):1125.Figure S1.Band-specific regional model fit | Separate regional regression models were applied to map MEG functional connectivity (AEC) to fMRI functional connectivity at each frequency band.Distributions of adjusted-R 2 are depicted for band-specific regional model fits and for the multiband model fit obtained by the original analysis.The multi-linear regional that combines MEG connectivity at multiple rhythms to predict regional fMRI connectivity profiles performs better than the band-specific models.[58-60, 71, 81].CTF is used to calculate peak localization error of a given source i as the Euclidean distance between the peak location estimated for source i and the true source location i on the surface model [59,81].No significant association is observed between the cross-modal correspondence R 2 map and peak localization error for LCMV and sLoreta.

Figure 1 .
Figure 1.Relating haemodynamic and electromagnetic connectivity | (a) A multi-linear regression model was applied to predict resting state fMRI connectivity patterns from band-limited MEG functional connectivity.The model is specified for each brain region separately, attempting to predict a region's haemodynamic connectivity profile from its electromagnetic connectivity profile.The dependent variable is a row of the fMRI FC matrix and the independent variables are the corresponding rows of MEG amplitude envelope correlation (AEC [18]) matrices for six canonical electrophysiological bands.For a model fitted for a given node i, the observations in the model are the connections to the other j = i regions.The regression coefficients b are optimized to yield maximum correlation between empirical and predicted fMRI connectivity profile of each brain region.(b) The overall relationship between fMRI and MEG functional connectivity is estimated by correlating the upper triangle of fMRI FC with the upper triangles of band-limited MEG FC, suggesting moderate relationship between the two across frequency bands.(c) Regional multi-linear model shown in panel (a) is used to predict fMRI FC from band-limited MEG FC for each brain region (i.e.row) separately.The empirical and predicted fMRI FC are depicted side-by-side for the regional model.The whole-brain edge-wise relationship between the empirical and predicted values is shown in the scatter plot.Each grey dot represents an edge (pairwise functional connection) from the upper triangles of empirical and predicted fMRI FC matrices.(d) A global multi-linear model is used to predict the entire upper triangle of fMRI FC from the upper triangles of the MEG FC matrices.The empirical and predicted fMRI FC are depicted side-by-side for the global model.The whole-brain edge-wise relationship between the empirical and predicted values is shown in the scatter plot.Each grey dot represents en edge from the upper triangles of empirical and predicted fMRI FC matrices.(e) The distribution of regional model fit quantified by R 2 is shown for regional model (grey histogram plot).The global model fit is also depicted for comparison (pink line).

Figure S2 .
Figure S2.Source localization error | MEG source localization error is estimated for (a) LCMV and (b) sLoreta source reconstruction solutions using cross-talk functions (CTF)[58-60, 71, 81].CTF is used to calculate peak localization error of a given source i as the Euclidean distance between the peak location estimated for source i and the true source location i on the surface model[59,81].No significant association is observed between the cross-modal correspondence R 2 map and peak localization error for LCMV and sLoreta.