Mapping electromagnetic networks to haemodynamic networks in the human brain

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 find 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, potentially reflecting patterns of laminar differentiation, recurrent subcortical input and neurovascular coupling. Correspondence between the two is largely driven by slower rhythms, particularly the delta (2-4 Hz) and beta (15-29 Hz) frequency band. Moreover, haemodynamic connectivity cannot be explained by electromagnetic activity in a single frequency band, but rather arises from the mixing of multiple neurophysiological rhythms. Collectively, these findings 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 [32]. 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.
How do electromagnetic and haemodynamic networks relate to one another? Although both modalities attempt to capture the underlying biological process (neural activity), they are sensitive to different physiological mechanisms and ultimately reflect neural activity at fundamentally different time scales [3,42,66]. Emerging theories emphasize a hierarchy of time scales of intrinsic fluctuations across the cortex [35,61,65,70], where unimodal cortex is more sensitive to immediate changes * bratislav.misic@mcgill.ca in the sensory environment, while transmodal cortex is more sensitive to prior context [5,17,18,45,47,50]. 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. First 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 [4,9,12]. Moreover, fMRI and MEG/EEG yield comparable fingerprinting accuracy, suggesting that they encode common information [22,26,31,68]. Finally, global edge-wise comparisons between fMRI networks and electrocorticography (ECoG) [8], EEG [25] and MEG [36] 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 [54]. 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 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 [13]) 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 (pair-wise 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). 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 [77]; 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 regionally 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 ( Fig. 1 [79]). 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 func-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 [55]. 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 [85]. 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. Intrinsic networks: vis = visual; sm = somatomotor; da = dorsal attention; va = ventral attention; lim = limbic; fp = frontoparietal; dmn = default mode. tional 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 [13]; 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 multiplexed to give rise to regionally heterogeneous haemodynamic connectivity. Note that the exponential inter-regional Euclidean distance trend is removed from all connectivity matrices before fitting any model, given the distance-dependence of both fMRI and MEG functional connectivity data (Fig. S1). All findings are reported for distance-corrected data.
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.1 to r = 0.3 ( Fig. 1b; r denotes Pearson correlation coefficient). The regional multi-linear model fits range from adjusted R 2 = −0.01 to R 2 = 0.57 (R 2 denotes coefficient of determination; hereafter we refer to adjusted R 2 as R 2 ), suggesting a close correspon-dence in some regions and poor correspondence in others (Fig. 1c,e). 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.1 (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 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 functional architecture of the cerebral cortex. To assess this hypothesis, we used the NeuroMaps toolbox [57] to compare the cross-modal R 2 map with the unimodal-transmodal principal functional gradient, estimated using diffusion 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 [2,14]. (a) The overall contribution of each frequency band is depicted for the regional model (box plots). Slower bands (delta and beta) 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. (c) The most dominant predictor (frequency band) is depicted for each brain region, confirming overall higher contributions from slower rhythms across the cortex. Faster rhythm appear to contribute to model fit in medial prefrontal regions and primary visual cortex. (d) Frequency band contribution to the regional cross-modal correspondence is shown separately for different rhythms across the cortex (95% intervals). map embedding on the group-average fMRI connectivity [53,55]. The two are significantly anti-correlated (Spearman rank correlation coefficient r s = −0.50, p spin = 0.0001), suggesting strong cross-modal correspondence in unimodal sensory cortex and poor correspondence in transmodal cortex. Finally, we stratify regions by their affiliation with macro-scale intrinsic networks and computed the mean R 2 in each network [85]. 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.

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 [2,14]. 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 and 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 regionally 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 delta and beta bands, contribute the most to model prediction across the cortex, while fast rhythm connectivity patterns appear to be mainly contributing to the model fit in medial prefrontal regions and primary visual 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 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) [44]. 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). (c) MEG connectivity is computed using an alternative, phase-based connectivity measure (phase locking value; PLV [52,60]), rather than the amplitude envelope correlations (AEC) (Fig. S2). The two FC measures yield similar cross-modal correspondence maps (Fig. 2a). (d) We repeated the analysis by fitting the multi-linear model on MEG and fMRI connectivity data without removing the exponential inter-regional distance effect. The results are largely consistent with and without distance correction (Fig. 2a). All brain maps are shown at 95% intervals. 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) [44]. 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.76, p spin = 0.0002). 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.84, p spin = 0.0002). 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.
Second, we compute MEG connectivity using an alternative, phase-based connectivity measure (phase locking value; PLV [52,60]), rather than the amplitude envelope correlations (AEC). The two FC measures yield similar cross-modal correspondence maps ( Fig. 4c; r = 0.77, p spin = 0.0001). Finally, we also consider how the results would look if we did not remove the inter-regional distance effect. Although we consider the results with the distance effect removed to be a more direct test of our hypothesis, we nevertheless include these results for completeness. As before, the results are largely similar with and without distance correction ( Fig. 4d; r = 0.69, p spin = 0.0001).

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 unimodal-transmodal cortical hierarchy. 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 [50,55,58]. Indeed, similar spatial variations are observed for multiple micro-architectural features, such as gene expression [15,34,44], T1w/T2w ratio [49], laminar differentiation [82] and neurotransmitter receptor profiles [33,40]. 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 [35,61,70,80,84]. In particular, multiple studies have reported faster timescales of intrinsic neural activity in unimodal sensory cortex, and slower timescales in transmodal cortex [35,61,65]. 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 [5,17,18,45,47,50]. 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 [27].
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, 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 [6,38,75]. Patterns of laminar differentiation -which define the unimodal-transmodal hierarchy [41,46,82] -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 [28,29,43,48,59,76]. As a result, patterns of haemodynamic connectivity may systematically decouple from the underlying electromagnetic connectivity patterns. Another related possibility is that unimodal cortex receives greater recurrent subcortical input [21,83]. This may result in more sustained and stable oscillations, and a greater correspondence between electromagnetic and haemodynamic connectivity. In contrast, transmodal regions have weaker recurrent connectivity with subcortex, yielding more flexible and less sustained electromagnetic dynamics, resulting in decoupling of haemodynamic and electromagnetic connectivity.
Throughout the present report, we find that fMRI networks are best explained as arising from the multiplexing 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 multiplexed model that takes into account all bands simultaneously. Previous work on cross-frequency interactions [30] and on multi-layer MEG network organization [11] 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 multiplexed 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 delta and beta. Multiple authors have reported that -since it captures slow haemodynamic coactivation -fMRI network connectivity would be mainly driven by slower rhythms [12,25,30,54,66]. 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 [16,73]? 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 [79]. Second, the present results show how the two modalities are related in a task-free resting state, but what is the relationship between fMRI and MEG connectivity during cognitive tasks [51]? 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 [62,81]. Second, to ensure that the analyses were performed in the same participants 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 multiplexing 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 [77]). The data includes resting state scans of about 6 minutes long (sampling rate = 2034.5 Hz; 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 ( [74]; http://neuroimage.usc.edu/brainstorm).
The MEG recordings were registered to the structural MRI scan of each individual using the anatomical transforma-tion matrix provided by HCP for co-registration, following the procedure described in Brainstorm's online tutorials for the HCP dataset (https://neuroimage.usc.edu/ brainstorm/Tutorials/HCP-MEG). The pre-processing was 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. 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 time-series were then used to estimate functional connectivity with an amplitude-based connectivity measure (amplitude envelope correlation; AEC [13]) and a phase-based connectivity measure (phase locking value; PLV [52,60]) from Brainstorm. AEC and PLV 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) frequency bands). Finally, the vertex-level connectivity data were parcellated into 400 cortical regions using the Schaefer-400 atlas [69]. Group-average MEG functional connectivity matrices were constructed as the mean functional connectivity across all individuals for each frequency band.

Resting state functional MRI
The functional MRI data were pre-processed using HCP minimal pre-processing pipelines [37,77]. Detailed information regarding data acquisition and preprocessing is available elsewhere [37,77]. 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) [24]. 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 [24,67]. The pre-processed timeseries were then parcellated into 400 cortical areas using Schaefer-400 parcellation [69]. 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.

Multi-linear model
A multiple linear regression model was used to assess regional associations between haemodynamic (fMRI) and electromagnetic (MEG) functional connectivity ( Fig. 1 [79]). 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 (F C(δ) i , 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). Also note that all multi-linear models are fit on distancecorrected data throughout the manuscript.

Region-level cross-validation
Region-level cross-validation was performed to assess out-of-sample model performance [44]. Briefly, we pseudorandomly split the connectivity profile of a given region (e.g. node i) into train and test sets based on spatial separation. 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 [2,14] (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,56]. We used the Schaefer-400 atlas in the HCP's fsLR32k grayordinate space [69,77]. 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 [71,78,79]. 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.

Code and data availability
Code used to conduct the reported analyses is available on GitHub (https://github.com/netneurolab/ shafiei_megfmrimapping). Data used in this study were obtained from the Human Connectome Project (HCP) database (available at https://db.humanconnectome. org/). Figure S1. Exponential relationship between functional connectivity and Euclidean distance | MEG functional connectivity (AEC) and fMRI functional connectivity are plotted against inter-regional Euclidean distance. Functional connectivity decreases exponentially with increased Euclidean distance. An exponential model is fitted to each connectivity matrix separately. The model parameters are then used to remove the distance effect from the data. Figure S2. Phase-based MEG functional connectivity | MEG connectivity is estimated using an alternative, phase-based connectivity measure (phase locking value; PLV [52,60]), rather than the amplitude envelope correlations (AEC; Fig. 1) for different frequency bands.