Resolving the connectome – Spectrally-specific functional connectivity networks and their distinct contributions to behaviour

In the absence of a task, the brain at rest spontaneously displays activity that reflects features of the underlying neural substrate. Examination of inter-areal coupling of resting state oscillatory activity reveals that the brain’s resting activity is composed of functionally connected networks, which differ depending upon oscillatory frequency, suggesting a role for carrier frequency as a means of creating multiplexed, or functionally segregated, communication channels between brain areas. In a group of 89 participants we examined spectrally-resolved resting-state connectivity patterns derived from MEG recordings to determine the relationship between connectivity intrinsic to different frequency channels and a battery of over a hundred behavioural and demographic indicators, using canonical correlation analysis. We demonstrate that each of the classical frequency bands in the range 1-40Hz (delta, theta, alpha, beta and gamma) delineates a subnetwork that is behaviourally relevant, spatially distinct, and whose expression is predictive of either positive or negative individual traits.


Introduction
Even in the absence of any task, the brain expresses patterns of functional connectivity. Such restingstate activity is thought to reflect intrinsic features of the underlying neural architecture 1 . Thus, spontaneous activity of the brain provides an insight into endogenous features, which are increasingly understood to be substantial contributors to individual differences in behaviour. A number of recent studies have highlighted the extent to which resting state connectivity metrics reveal individual "wiring patterns" of the brain that are significantly predictive of behavioural, and even demographic features. Notably, resting state fMRI (RS-fMRI) connectivity patterns have been shown to be predictive of taskrelated activation patterns 2 and even to constitute uniquely identifying fingerprints that can distinguish between individuals 3 . Such findings suggest that resting-state activity patterns constitute a framework that determine a large degree of individual differences in brain responses to tasks, and of individual differences in behaviour. In a sense, they may be considered the brain's "priors" 4 for determining how it responds to input.
A recent effort to relate whole-brain resting state fMRI connectivity to individual differences in behaviour employed canonical correlation analysis (CCA) to establish how inter-individual differences in connectivity can be related to differences in a broad battery of behavioural and cognitive variables 5 (nonimaging subject measures, SMs), using a sample of 461 participants from the Human Connectome Project 6 . CCA was used to find the maximal correlations between combinations of variables in the connectivity and SMs, revealing a canonical mode of covariation that is highly positively correlated with positive cognitive traits as well as desirable lifestyle indicators. This is a fascinating finding, which highlights the extent to which large-scale neural connectivity of the human brain is functionally relevant. However, this discovery invites further and deeper inquiry: fMRI is rightly prized for its high spatial resolution, but its temporal resolution is inherently limited by the sluggishness of the haemodynamic response that it records. Consequently, fMRI cannot be used to discriminate between neural processes that occur at timescales in the sub-second range. These are, however, of profound relevance to further elucidating the nature of resting state brain activity, and its relationship to behaviour.
A characteristic feature of the human brain's electrical activity at rest is the presence of oscillatory signals. Oscillations are the product of repetitive, or cyclical patterns of brain activity, which are observed to occur at a number of different frequencies. Oscillations exist at multiple spatial scales, from the single neuron, to local field potentials and large-scale magneto-electric brain signals. Different oscillatory frequencies are thought of as representing different channels in which neural activity is communicated and might therefore fulfil different functional roles 7 . Even at a relatively coarse level of spatial resolution, an extensive literature of M/EEG or electrophysiological studies has established how neuronal processes occurring in different frequency bands have different functional roles -from the relatively slow sleep rhythms (e.g. delta waves in NREM sleep 8 ) to memory processes in the hippocampus 9 (mostly in the 4-7 Hz theta range), to processes related to sensory processing (covering the alpha and beta range of 8-30 Hz 10,11 ) or more cognitive functioning (e.g. working memory [12][13][14] and feature binding during visual and other task processing 15,16 , as reflected in higher frequencies, such as the gamma band, i.e. 30Hz and above.
Oscillations can be characterised in terms of frequency, amplitude and phase. Much of the research into the functional relevance of neural oscillations has focused on properties like phase or amplitude, which have been shown to index different states of neuronal excitability [17][18][19][20][21][22][23] and ultimately, to predict behavioural and perceptual performance in a variety of experimental task conditions, domains and sensory systems [24][25][26][27] . While such research has been fruitful, if we are to consider different oscillatory frequencies as carriers of different information, an analysis of the networks in which they occur is necessary.
Synchronised patterns of modulation of activation at separate loci are believed to reflect coupled activation, and consequently to define functionally connected areas, and ultimately, networks. By examining resting state activity in spectrally-distinct frequency ranges, it is thus possible to define largescale brain connectivity networks that are frequency-resolved. Recent efforts have identified large-scale (i.e. whole-brain) electrophysiological networks with topographical and spatial structures comparable to the resting state networks established by fMRI studies. For example, the motor networks or visual components, analogous to those previously described in fMRI 28,29 have been reported in MEG [30][31][32] . Interestingly, these networks, while spatially similar to fMRI-derived RSNs, had distinct spectral properties. That is -they arise when looking at coherent activity in different frequency regimes: the motor network being most pronounced in the beta range (14)(15)(16)(17)(18)(19)(20)(21)(22)(23)(24)(25)(26)(27)(28)(29)(30) whereas networks attributable to visual areas were most visible in the alpha range 8-12 Hz. The finding of different topographies in spectrallyresolved networks hints at the potentially different function these spontaneous networks might serve, but beyond the more intuitively interpretable sensory-and motor-related networks mentioned above, the roles and particular functional relevance of the spatial networks that emerge at different frequencies, remains elusive. At this point, it is worth highlighting that all neural communication is constrained by the underlying anatomical connectivity, and that any network of information transfer that is discovered is necessarily a reflection of intrinsic brain architecture. Thus, any relationship elucidated between functional connectivity patterns and behaviour depend upon neuroanatomy, but they are not solely defined by it -the existence of connectivity in different frequency ranges is highly suggestive of the possibility that the anatomically constrained connections can be exploited to different ends.
The precise role of network communication at different frequencies is the subject of ongoing and increasing interest. One overarching notion is that of "multiplexing", viz. that different oscillatory frequencies comprise different communication channels allowing the same anatomical connections to be used for different functional roles 33 . An influential theory with a very broad scope -'communication through coherence' (CTC) -proposes that the observed frequency-constrained networks, i.e. networks that share a common rhythm, reflect activity of cell assemblies that are in communication with each other -e.g. when binding different sensory modalities together 15 . In particular, coherence in the faster gamma band (40 Hz and above) is well examined and seems to confirm the proposals of CTC theory. CTC also affords the existence of similar mechanisms in other frequency bands, although these are less well explored and less is known about their likely functional role.
Another recent focus of enquiry, underscoring the importance of neural activity in different frequency bands has been how different frequencies reflect the activity of neurons having different roles in the directionality of neuronal communication. For example, it is known that bottom-up processes are preferentially mediated by higher-frequency (e.g. gamma) activity, while top-down processes have principally been associated with rhythms in lower frequencies 34 . In sum, there is compelling evidence to believe that spectrally confined, oscillatory activity may serve different purposes in network communication, be it in directing information flow or in defining the neural assemblies that preferentially operate as coherent ensembles.
Despite these advances, so far, a holistic analysis of the functional role of oscillatory networks, taking into account the full spectral richness of the underlying constituents of human spontaneous brain activity, and its relationship to an extensive range of behavioural factors has been beyond our reach. Here, we seek to expand upon and complement the existing rs-fMRI findings by exploring the relationships between the same set of SMs and spectrally-resolved connectomes derived from rs-MEG, also made available as part of the HCP dataset 35 .
Specifically, we address the following questions. Firstly, can we identify and characterize a global mode linking brain connectivity and behaviour using MEG? Secondly, is this mode comparable to that revealed for rs-fMRI? And thirdly, what additional information can we glean from the spectral richness of MEG afforded by its high temporal resolution? By exploring these issues, we will shed light upon the spatial structure and functional contributions of the spectral components of resting state connectivity and how they relate to a broad spectrum of behavioural indicators.

Results
2.1. MEG-derived functional connectivity reveals a global brain-behaviour mode along a positivenegative axis of subject measures.
In order to provide a functional characterization of this mode of maximum covariation between SMs and connectivity, we correlated individual subjects' scores on each of the SMs with the observed canonical variate for the SMs. This procedure reveals the way in which the mode that best aligns the full set of SMs and the connectivity patterns explains the inter-individual variation in each of the SMs. A positive correlation between an SM and the canonical variate suggests that the higher a subject scores on this SM, the larger this subject's score for the canonical variate (and vice versa). Conversely, a negative relationship between an SM and the canonical variate indicates that higher scores on that SM are associated with the subject score of the canonical variate being smaller or more negative.
The behavioural variables most strongly positively and negatively associated with this mode are shown in Figure 1A. The SMs are arranged along a positive-negative axis that is qualitatively the same as that previously described for an fMRI derived CCA mode 5 . It ranges from subject measures that entail tobacco consumption and somatic problems (which may legitimately be interpreted as 'negative') to subject measures indicating high cognitive performance (at the other end of this ranking, interpreted as 'positive'). The nature of the positive relationship between the mode and performance on a working memory task is further illustrated in Figure 1B, where a clear pattern emerges demonstrating that higher individual scores on the SM and connectome canonical variate are predictive of superior list sorting performance. This outcome alone is remarkable -despite using a different imaging modality and having access to a substantially smaller sample (89 vs. 461 data sets) the nature of the mode revealed is highly consistent with that previously demonstrated for the same set of SMs and rs-fMRI connectivity (a descriptive comparison with rs-fMRI, as well as the impact of parcellation scheme and subdivision into the selected frequency bands is presented in Supplementary Figure 3).
Although this MEG-derived canonical mode resembles and partially even reproduces a previously identified fMRI derived brain-behaviour mode, we are at pains to point out that this does not automatically render it the 'universal' brain-behaviour mode. Given a different array of subject metrics, or a qualitatively different array of brain metrics, the principal brain-behaviour mode uncovered by CCA could be different to the one at hand. This possibility is, however, speculative. What we may state with confidence is that, independent of the underlying imaging modality, a canonical brain-behaviour mode, optimally matching spectrally resolved connectomes and subject measure, can be identified that arranges subjects on a positive-negative axis consistent with a study in which connectivity was derived from hemodynamic measures of brain activity.

The CCA Mode is Composed of Spatio-Spectrally Segregated Subcomponents
In a second step, we characterize the connectivity components of the identified brain-behaviour mode, for each of the five frequency bands and their connectomes. In order to do this, the connectivity patterns across subjects in each band were correlated with the individual connectome weights of the significant CCA mode. This reveals which connections (hereafter referred to as edges) of the connectomes are more, or less, related to the mode, and therefore more or less explanatory of the individual differences over SMs. The outcome of this analysis, thresholded to show only the top 0.5% most positively and the 0.5% most negatively related edges is shown in Figure 2. The correlation coefficients resulting from this analysis will be further referred to as edge modulations (EMs). The patterns of within-band edge modulations (EMs) significantly related to the mode are clearly differently distributed over the brain, and have distinctly negative or positive relationships with the mode in each of the analysed frequency bands. The observed patterns, i.e. the EMs in each band, do not merely mirror the group-average functional connectivity (see Supplementary Figure 4), but are different with respect to the spatial characteristics (i.e. when comparing peak connectivity vs. peak edge modulations, quantified in Supplementary Figure  5).

Figure 2.
Edge modulations (i.e. edges whose connectivity on a subject-by-subject level covaries with behaviour) for each frequency band, from delta, theta, alpha and beta to gamma band. Edge modulations of frequency bands are colour coded (2 colours per band, positive and negative, encoding the correlation coefficients, in a range from -.5 to .5). Edges rendered on the brain templates are thresholded at the 0.5 th bottom and 99.5 th top percentile for visualization. Each of the bands shows a strong preference for either positive or negative relationship to the mode, but not a mixture of both. For a more complete picture, the histograms in the bottom row depict the distributions of the (unthresholded) correlation coefficients, with comparison to null distributions generated by a permutation test (n=10000, in grey). Figure 3A shows the top 50 edge modulations across all five frequency bands, providing an overview of the diversity of connectivity patterns in each of the frequency ranges analysed (we will further abbreviate edge modulations by 'EMs', or 'EM maps' when referring to whole-brain maps of EMs). To reveal the relative functional importance of each node in terms of its relationship to the mode, averaging was performed over all the EMs in which the respective node is involved. We thus derive maps of accumulated EMs, which are shown in Figure 3B. As described above, the distributions of EMs across bands is primarily unipolar, each of the five frequency bands preferentially showing either a positive or a negative relationship with the mode, but not both ( Figure 3B). In order to further characterize the five band-specific networks of EMs that we revealed (shown in Figure  2), we used the NeuroSynth database to establish the meta-analytic correlation between each of the band-specific sets of EMs and a set of 50 topics derived from the database. We carried this out for maps of the implicated nodes ( Figure 4) and the accumulated EMs (Supplementary Figure 6). This analysis enabled us to obtain a broad quantitative functional characterization, based upon spatial distribution, of what subset of cognitive functions each of the networks most closely reflect.

Connectivity of spectrally-confined brain activity is preferentially 'good' or 'bad', but not mixed
The discovery of frequency-specific connectivity components that correlate either negatively or positively with the CCA-derived mode invites the conclusion that some connectivity patterns are "good", while others are "bad". This interpretation is somewhat consistent with previous reports from rs-fMRI data 5 , which also showed that certain subcomponents of the whole-brain connectome are positively (especially areas pertaining to the default mode network), and others negatively, associated with the positive-negative axis exposed by CCA (mostly visual areas).

Delta-band and alpha-band EMs show negative relation to mode
Alpha and Delta-band connectivity were almost exclusively negatively associated with the CCA mode. The parcels with the highest incidence of significant connectivity-behaviour relationships in alpha are principally concentrated in visual-posterior areas, consistent with the usual distribution of alpha power at rest. Here, after having regressed out spectral power variation across subjects, we may conclude that higher connectivity in the alpha band, absent the biasing effects of the typically higher power in the alpha band compared to the rest of the spectrum, is predictive of lower cognitive performance and higher incidence of negatively associated health and behaviour indicators. Similarly, delta-band connections correlate negatively with behaviour, albeit with the majority of implicated nodes spatially very differently distributed compared to alpha-band EMs, incorporating bilateral inferior posterior and anterior temporal lobes and inferior and orbitofrontal areas.
The negative relationship of the alpha-band posterior EMs to the canonical mode is comparable to the accumulated EMs known from fMRI derived CCA analysis 5 . Here, by means of spectrally resolved CCA analysis, we can now assign these negative-posterior EMs a neuro-spectral profile, identifying alpha oscillations as being the most relevant across the whole-spectrum connectome. The components residing in the delta-band are less easy to match to previous results. In general, these findings are suggestive of increased low-frequency neural coupling. Usually, pronounced delta oscillations are found during sleep, in the awake human brain it is visible mostly in pathological situations such as in the presence of lesions or tumours. Excessive delta coupling may also be a sign of slowed theta oscillations, a phenomenon often observed in pathological cognitive decline, such as Alzheimer's disease, caused by cholinergic loss 37 . While we do not propose that the young, healthy cohort examined here are suffering from any neurological pathologies, it is intriguing to note the fact that slow oscillatory activation is not typically associated with optimal brain function during wakefulness.

Theta-, and beta-band EMs show positive relationship to mode
Beta, Theta and Gamma-band EMs exhibit a substantially different relationship to the canonical mode, compared to those of Alpha and Delta. This is the case both in terms of the polarity of the relationship, and the spatial distribution of the relevant nodes. For the observed canonical mode, subjects scoring high in frontal and prefrontal beta-band connectivity also score high in the identified canonical mode, indicating high cognitive performance. The patterns that arise from beta-band accumulated EM maps, especially in medial areas, show good correspondence to previous results for fMRI connectivity 5 . For theta EMs and accumulated maps, spatially we see a more pre-central pattern with a functionally similar interpretation -the more pronounced connectivity is in these edges, the higher the subject scores in the canonical mode. The properties of gamma-band EMs are particularly noteworthy. While the overall gamma EM network is somewhat patchy, featuring what appear to be a greater number of short-range connections, its connectivity is still relevant for the global mode of behaviour. It also seems to have fewer sinks or hubs where connections emerge or terminate, resulting in the absence of significant accumulated EMs, in contrast to the more clustered network patterns in frequency bands at lower frequencies. Considering that gamma in general is known to be an index of local activation, offering often relatively fine spatial resolution (e.g. down to mapping the representation of individual digits 38 ), this is not too surprising, given the more global character of the identified brain-behaviour mode.
This striking dichotomy of either good or bad connectivity -when looking at spectrally resolved connectomes -is both novel and unexpected. Previous results -using fMRI connectivity agnostic of the underlying frequencies of neuronal communication -have so far only shown patterns of mixed positive or negative relationship to the observed brain-behaviour mode. We tested whether this segregation does arise from the CCA trivially by creating EMs from surrogate data, by permuting the subject-labels and recomputing the EMs. The resulting null-distributions of correlation coefficients are shown in Figure 2 (bottom row); these do not exhibit the polarity preference and skewedness of the actual EMs.

Nodes that are relevant to the mode differ from other nodes in their envelope dynamics
The functional connectivity metrics employed here do not provide insights into the temporal dynamics of the signal. These can be accessed in a variety of ways, and we chose to examine the modulation spectra of the envelope time-series that constitute the functional connectivity examined. This approach captures local fluctuations in amplitude of a given frequency band in a manner that is consistent with notions of neural information. There is increasing evidence that variability in neural signals may be a determiner of cognitive performance 39 , and we explored the difference between the relevant and irrelevant nodes of the above-described EM networks in terms of these spectral dynamics.
In theta, beta and gamma bands, nodes with strong, positive links show more modulation power (Supplementary Figure 7), i.e. more variability in their envelope dynamics than with less positive links (i.e. the nodes only showing weak links to the mode). Conversely, in the alpha band, nodes with the strongest negative links to the mode show more variability than nodes with a less negative link to the mode (which in case of alpha means a weaker link). These results show that far from strength of connectivity alone being the determiner for behavioural relevance of an edge, the local resting state spectral dynamics of the nodes also play a role.

Discussion
The concept of spectrally resolved connectivity, which can also be considered neuronal communication at multiple time-scales or channels, is not new. Some models of how neuronal assemblies communicate have focused on activity confined to certain frequency bands 40 , while more recent perspectives also incorporate multi-channel or multi-frequency scenarios 15 . Furthermore, empirical studies, mostly MEG, have also increasingly shown connectivity and networks at multiple frequencies or showing features or states with specific spectral features or profiles 41,42,43,44 ). What was missing hitherto was the functional interpretation of these spectrally resolved features or parameters, due to the absence of an attested link to behaviour. Here, to the best of our knowledge, for the first time, we have demonstrated the nature of the relationship between a spectrally resolved connectivity estimate to not just one, but an entire battery of behavioural and other subject variables.
Although the classification of brain activity into predefined frequency bands might be an oversimplification of the richness of neuronal dynamics in the brain, it can be a useful one. The fact that using these predefined, canonical, spectral ranges and the corresponding connectomes reveals a global brain-behaviour mode with little overlap between bands lends support to the idea of spectrally-confined functionally meaningful networks. Further studies looking into cross-frequency coupling or more complex types of neuronal communication might be needed to fully uncover the role of interactions in neural transient dynamics.
Further, the strong preference for either positive or negative relationships with the mode seems to be supportive of the existence of different types or 'channels' of neuronal communication, characteristic of and confined to these frequency bands. Furthermore, we do not find that a single frequency band completely dominates the CCA derived brain-behaviour mode. Rather, analyses of both the functionally relevant patterns of connectivity -the EMs -as well as the accumulated EM maps, point to a relatively distributed and balanced contribution of network components indicative of behaviour.
Another point of interest is the fundamental and qualitative difference of the observed EMs in alpha and beta-band. While both rhythms, at rest, express functional connectivity network patterns that are relatively amenable to interpretation -they resemble visual or motor related networks -the resulting EMs, for the observed mode at least for beta, are much more frontal or prefrontal localized and can be less easily assigned directly to sensorimotor functions. In addition, higher beta-band connectivity is beneficial, while higher alpha connectivity (in the predominantly visual nodes implicated in the CCA mode) is apparently deleterious for cognitive performance linked to the observed mode. This implies that more connectivity is not always better, but its impact on the mode (and thus behaviour) appears to be a function of frequency.
This fundamental issue of 'good' and 'bad' brain patterns is also related to another line of research, which has investigated how variability predicts healthy ageing and cognitive performance 39,45 . While the approach employed here emphasised connectivity rather than variability per se, the two concepts are not totally unrelated. For example, neuronal variability (or activity in the widest sense) is a necessary pre-condition for meaningful connectivity. What directly links our results to the concept of variability, is that the nodes implied in our functionally relevant networks (as seen in the accumulated EMs, for example) differ from others in their degree of envelope modulation at rest, which may be a sign of increased information transmission or communication. What is particularly noteworthy is that this level of modulation is associated with this node being important or not; it does not predict whether this variability is good or bad -this is rather dependent on the frequency band the network is operating at.
The difference in the modulation of envelope dynamics and the absence of clear spectrally confined effects also touches upon another recent discussion -the debate of dominant neuronal time scales in the human brain. One of the studies supporting such a view suggests that there is evidence for the existence of neuronal processes with life-times on the order of 200 milliseconds 46 , using computational modelling of the neuronal responses underlying the relatively sluggish fMRI time-series. Other studies, using temporally highly resolved MEG data used statistical models such as Hidden-Markov-Modeling to infer life-times of neuronal states hidden in resting state time-series 47,48 . Both sources, while indirectly, seem to suggest the idea of dominant time-scales clearly below a second. In our case, when examining the more directly accessible envelope dynamics of the different band-specific EM networks, as identified by CCA, independent of whether they are beneficial or deleterious, we do not see evidence for dominant time scales. Rather, we observe very broad changes in the modulation spectra not bound to specific frequencies. In general, it is interesting to note that there is no marked peak in the relationship between behaviour and the cerebral signal for any of the carrier frequencies (i.e. the five bands employed here) or modulation frequencies. Of course, these findings are not directly comparable to and thus cannot exclude or contradict the presence of dominant time-scales as identified by other approaches. Yet the demonstration of the spectrally broad distribution of meaningful connectivity may be insightful for our understanding of (large-scale) neuronal communication.
Assuming the existence of multiple, functionally-relevant frequencies allows for the possibility of multiplexing, which is an elegant concept of how communication and information transfer may be realized in the human brain 33,49 . While we do find patterns of edges implicated in the canonical brainbehaviour mode, these are for the most part spatially distinct. This spatial separation of the EM networks at each oscillatory frequency band implies that with respect to the overall spectrally-delineated connectivity patterns, the behaviourally relevant components may not necessarily require multiplexed communication. This, of course, does not exclude that multiplexing occurs, it rather indicates that the spatial segregation of functional networks is relevant.
In conclusion, we have provided evidence that electrophysiologically derived and spectrally resolved connectivity present in MEG resting state data can be used to index the ranking of subjects across a large range of subject measures. It exhibits highly structured patterns that are functionally relevant and provides sensitivity comparable to rs-fMRI derived networks in explaining variability of subject measures across a wide range of different features and domains. Furthermore, the additional dimension afforded by spectrally resolving connectivity measures opens up new avenues into a better and more holistic understanding of the roles and contributions of brain rhythms at rest that ultimately will help cast light on the intrinsic features of the brain that determine positive and negative cognitive and behavioural traits. Given recent breakthroughs in using frequency-tuned transcranial electrical stimulation techniques to enhance cognition 50 , outlining the roles of frequency-specific connectivity networks is all the more timely, and provides a basis for considering potential network-level targets for novel neuromodulatory interventions.

Subjects and data
Resting state MEG data from the Human Connectome Project were used for this study. 89 subjects had complete resting state recordings (i.e. 3 recording sessions of about 6 mins). Mean age: 29 +/-4 years, Sex distribution: 41 female / 48 male. Subjects were recorded in supine position and instructed to remain relaxed, with eyes open.

Subject and behavioural measures.
Following along the lines of a previous study showing CCA derived brain-behaviour mode with fMRI 5 , we extracted a large number of subject measures from the HCP data set. We excluded gender related subject measures and structural measures, but kept all other measures (cognitive, emotional, sensory performance tests, psychiatric and personality tests, family history of mental or other disorders, consumption of alcohol, tobacco or other drugs, in-scanner (fMRI) task performance). From those, we excluded variables that did not fulfil the following criteria: These exclusion criteria resulted in 182 behavioural variables that were used for further analysis.
We also deconfounded, i.e. regressed out, effects of age, height, weight, blood pressure (systolic / diastolic) and brain volume for all subjects. All variables were further subject to rank-based inverse normal transformation and demeaned.
MEG pre-processing and parcellation, source leakage correction.
MEG acquisition. Magnetoencephalography data was recorded by a whole-head MAGNES 3600 (4D NeuroImaging, San Diego, CA). Data were acquired in three runs of resting state recordings, with 6 minutes each. The MEG recorded from 248 magnetometer sensors, with 23 reference channels. Sampling rate was 2035 Hz and data were down-sampled for further processing to 200 Hz. The standard pre-processing pipeline as offered by the HCP was used which applied independent component analysis to remove potential artefacts from ocular, muscular or cardiac sources. For source reconstruction, single shell volume models were used, based on the individual anatomical MRI T1 images. Linearly constrained minimum variance (LCMV) beamforming was employed, resulting in regular 3D grids in normalized MNI source space with a resolution of 8mm 3 . After source reconstruction, which resulted in a grid of 5798 source voxels, we parcellated the data into 100 source parcels using a semi-data driven parcellation approach. Parcels were derived by starting out from an initial 246-regions anatomical brain atlas 51 covering the two hemispheres (including subcortical cerebral areas) and reducing them to 100 parcels. Reduction of numbers of parcels was achieved by a k-means clustering of group-(and session) average correlation matrices of the original resulting 3 x 246 time-series. Parcels that clustered (because of similar absolute correlation coefficients) were merged into common parcels. Parcellation at higher dimensionality, e.g. with N=200, lead to failure of the multivariate source leakage correction needing full rank covariance matrices. A visualization is provided in Supplementary Figure 1, for a list of areas see Supplementary Table 1).
We chose to use this data-driven customized parcellation and not the fMRI based group-ICA as for the fMRI data, with the rationale that the fMRI-derived parcellation is a spatial ICA capitalizing on the large amount of spatial information available in fMRI data, resulting in sometimes fine-grained components, e.g. in visual areas. However, since this type of fine-grained spatial resolution is not always given in MEG data, similarly to the N=200 BrainNetome 51 derived parcellation from above, rank deficiency of the parcel-time series can render the source leakage correction infeasible.
To obtain parcel-time series, the first principal component over all voxels belonging to any given parcel was extracted. After extraction, resulting parcel time-series underwent reduction of source leakage by using a multivariate leakage reduction approach 52 , effectively orthogonalising all parcel time-series and removing zero-lag correlations.

Analysis of functional connectivity
After source-leakage correction, we estimated functional connectivity (FC), separately for each individual and session, and on a group level. To do so, we first computed envelope dynamics in 20 contiguous frequency bins, with a frequency range of 0.5 to 40 Hz, and a resolution of 2 Hz, using the resulting envelopes in each bin after Hilbert transform of the narrow-band filtered time-series. Within each bin, full linear Pearson correlation was computed between all possible pairing of parcel time-series, resulting in a 100 x 100 connectivity matrix for each resting state run, frequency bin and subject. Connectivity matrices were based on the full correlation matrices computed in each frequency bin. For estimation of functional connectivity, we used the OSL toolbox (Oxford Software Library), freely available at https://github.com/OHBA-analysis/osl-core. After FC estimation, all three runs were then averaged to increase stability of FC estimations.
We observed that functional connectivity in spectrally resolved MEG data is not always independent of spectral power or frequency, but correlated with group-average correlations, being higher in lower frequency bands. To ensure that further findings are not simply confounded by these effects, we removed this systematic variation by pair-wise regressing out of spectral power from connectivity (correlation coefficients) from all possible pairings (i.e. edges). This resulted in subject-specific connectivity estimates without low-frequency bias. Subsequently, the resulting connectivity matrices were normalized to zero mean unit variance.
Analysis of brain-behaviour mode by means of canonical correlation analysis (CCA).
In contrast to the fMRI data set -which is missing the dimension of frequency -we here have 5 frequency bands that are taken into account. First, PCA is performed on the concatenated data sets (i.e. on the 24750 x 89 concatenated connectivity matrix, see Supplementary Figure 2 for visualised approach). By retaining 22 principal components we obtain a ratio of features to subjects (or samples) of 1 to 4, avoiding overfitting when using CCA. This results in a 22 x 89 connectivity matrix. Subsequently, PCA is performed on the behavioural data, reducing the 182 x 89 matrix of subject and behavioural measures into a 22 x 89 matrix as well. CCA is then performed on these two matrices, finding the best linear weighting of connectivity and behavioural features to maximize brain and behavioural covariation across subjects. Here, CCA results in a maximum of 22 canonical modes, of which statistical significance is established by permutation testing (n = 10000). For permutation testing, subject labels were swapped while considering family structure in the data. For running CCA and permutation testing we used an adapted version of the scripts made available in a previous CCA study 5 .

Visualisation of post-CCA results
For characterization of the observed CCA mode(s) with respect to the underlying behavioural variables and connectivity, the resulting behavioural CCA weights (i.e. one per subject, n=89 in total) were correlated with the full, deconfounded behavioural variables (i.e. 182 x 89), and the mode-derived connectivity weights (n=89), in turn, were correlated with the full deconfounded individual FC networks (i.e. the 5 x 4950 x 89 connectivity matrix). For visualization of connectivity and the resulting edge modulations, we used the freely available Circos software suite 53 (http://circos.ca/software/download/circos/). We also derived accumulated EMs from the data where we collapsed one dimension in the EM matrices by averaging. For visualizing the results, we used the parcellation functionality of OSL (Oxford Software Library), available at https://github.com/OHBAanalysis/osl-core .

Meta-analytic Characterisation of EM network parcels using NeuroSynth
In order to describe the functional systems in which the five behaviourally-relevant connectivity parcellations -i.e. the EMs -are embedded, we carried out a meta-analysis using the NeuroSynth database 30 , similar to recent efforts to characterize functional gradients in brain connectivity 54 . We carried out analyses using the "decode" tool provided by NeuroSynth, to evaluate the relationship between a list of topic terms (derived from the database) and the maps revealed by correlating the subject weights of the behavioural mode with the connectivity values in each of the five frequency bins. We conducted two analyses in each of the five frequency bands. In the first, we collated all the nodes that were significantly (using a threshold including the top and bottom 0.5 th percentile of all correlation values) involved in any number of connections, shown in Figure 2. The second analysis examined the aggregated EMs, thresholded at the 10 th and 90 th percentile, as displayed in Figure 3B. For both of the two representations, parcels containing significant nodes were recorded and an aggregate binary mask was composed reflecting all the parcels shown to be implicated, for each frequency band. The binary maps were then evaluated for their relationship with a list of 50 topic terms derived from the NeuroSynth database 55 , 56 . Of the 50 terms, a number were excluded from the analysis a priori based on a lack of interpretability. The full term list is provided in Supplementary table 2.

Characterisation of EM network parcels by analysis of modulation spectra
To understand what renders the edges relevant in the context of the observed CCA mode different to other edges we performed by analysis of the modulation spectra of the dynamics of the band-specific envelopes underlying the function connectivity used for CCA and following analyses. Using a cutoff at the 10 th and 90 th percentile, we average modulation spectra of the 10 most important edges, i.e. with the largest accumulated EMs, compared against the averaged modulation spectra of the 10 edges with the smallest influence (i.e. lowest EMs).

Acknowledgments
Funded by the Swiss National Science Foundation (grant: PP_163726 awarded to A.H-A). Data were provided by the Human Connectome Project, WU-Minn Consortium (Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657) funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at Washington University. We are grateful to Steve Smith and colleagues for sharing the code of their HCP-fMRI CCA study 5 online. The authors are also grateful to Daniel Margulies for assistance in preparing the NeuroSynth meta-analysis and for sharing crucial elements of the NeuroSynth database.

Supplementary Methods
Comparison of MEG CCA mode to other parametrisations, parcellations and comparison to fMRIderived CCA modes We analysed different types of MEG and fMRI derived connectivity to obtain an overview about the importance of the chosen approach and how they compare to the one approach chosen for this study. Estimation of functional connectivity was varied in the following ways: Approach 1: This is the approach presented and used throughout the study, as described in Methods.
Approach 2: Originally, we used a higher frequency resolution and binning, i.e. functional connectivity estimation was based on 20 frequency bins (with 2Hz steps each). Remaining steps were identical to approach 1.
Approach 3: Same as approach 2, but without any source leakage correction performed. Usually it is recommended to perform source leakage correction of some sort to avoid spurious correlations.
Approach 4: Instead of 5 bands or 20 bins, here we only used broadband-defined connectivity, that is the linear correlation of Hilbert transformed envelopes in the frequency band up to 40 Hz.
Approach 5: No actual connectivity was used here, instead we used the average resting state spectral power in each of the 5 frequency bands as feature to feed into the CCA.
Approach 6: This approach used functional connectivity estimation as in 2, with a different parcellation.
Here we used an fMRI derived 50 parcel parcellations (called 'PTN parcellation', as offered by the HCP platform as part of the HCP S1200 release on https://db.humanconnectome.org/). Approach 7: This is using fMRI connectivity instead of MEG (as used elsewhere in this study), as extracted from the HCP S1200 release, for the 89 subjects used in the present study, using the 100 parcel fMRI derived HCP parcellation.
For these different approaches of estimating functional connectivity, we only varied the part of FC estimation, the following steps as explained in the Methods section above were the same -using the same subjects, using the same number of principal components retained, using CCA in the same manner (testing against the same behavioural variables) and using the same permutation tests for establishing significance. Results are given in Supplementary Figure 3, preserving the order of the analyses described above.
We wanted to understand and characterize better the spatial network structures emerging as being relevant to the CCA mode, and how these structures relate to the underlying functional connectivity. For this we performed two analyses: First, we linearly correlated both group-average FCs and (group-based) EMs separately, across frequency bands to estimate the spatial autocorrelation across frequency boundaries. Then, we correlated FCs and EMs with each other, both for corresponding frequency bands and also across frequency bands. Results are shown in Supplementary Figure 5.

Supplementary Figures
Supplementary Figure 1. Visualisation of the used parcellation (n = 100). Parcel identities can be found in Supplementary Table 1, with anatomical labels.
Supplementary Figure 2. Illustration of analysis approach. Functional connectivity as defined by envelope correlations within five conventional frequency bands are extracted from n=89 subjects. For each subject, these connectivity features are concatenated and subjected to group PCA, retaining the first 22 principal components per subject. The same approach is used for reducing dimensionality of the subject measures of all subjects, retaining the first 22 principal components. These two matrices are then subject to canonical correlation analysis (CCA), which identifies within each matrix the optimal linear weighting to maximize correlation of features between the two sets of data (brain vs behavior). Rows represent modes where the linear combination of connectivity and behavior correlated most strongly (first mode here is indicated by thin black line).
Supplementary Figure 3. Comparison of CCA results for different types of functional connectivity in MEG and fMRI for the 89 subjects in this study, derived from different types of MEG connectivity (1-6) or from fMRI (7), for details see Supplementary Methods. For each approach, correlation coefficients and resulting significance levels are indicated. Abbreviations & explanation of analyses: 1: Analysis with MEG connectivity as used throughout the study. 2: Same data, but instead of using 5 bands we used the complete frequency resolution available (20 2Hz-wide bins from 0.5 -40Hz), 3: Same analysis, but without source leakage correction (SLC), 4: broadband conn = connectivity derived from broadband activity, 5: only using power values from the resting state spectra in the corresponding frequency bands, 6: As 2, but with parcellations from fMRI derived PTN parcellation ('Parcels, time series and netmats' HCP release), 7: power = CCA applied to spectral power only (no connectivity), 8: fMRI connectivity extracted from the HCP S1200 release. All analyses have been performed with the same 89 subjects. Levels of significance as tested by permutation testing: ** p<=0.01, *** p<=0.001, n.s. = not significant.
Supplementary Figure 4. Group-average functional connectivity of all subjects, in order of frequency bands -delta, theta, alpha, beta and gamma-band. Colour codes for connectivity, i.e. the linear correlation coefficient of envelope time-series in each frequency band (see also Methods). Please note, that while we use the same colour scheme throughout the paper, here, all bands, including delta and alpha-band FCs, have positive correlations, and thus only the colour schemes encoding the positive correlations in each band are used.
Supplementary Figure 5. The spatial structure of average functional connectivity is more correlated across frequency bands (left panel) than the spatial structure of the observed resulting edge modulations (middle panel). Furthermore, also the correlation -within each frequency band -between average functional connectivity and group edge modulations (rightmost panel) is showing relatively little overlap. The strongest relationship between average functional connectivity and edge modulations is in the alpha band, however the relationship is negative (diagonal in right-most panel).  Figure 3B) and 50 topics from the NeuroSynth database. Results are comparable to the non-accumulated EMs (as depicted in Figure2), with gamma-band accumulated EMs not reaching threshold and thus being left out.
Supplementary Figure 7. Modulation spectra of the 10 edges with strongest accumulated edge modulations vs the 10 edges with lowest or most negative accumulated edge modulations (averaged over subjects and edges, error bar depicts standard error of the mean over subjects). Apart from the delta-band edges, we see more modulation for edges that are more relevant which is in theta, beta, gamma edges with positive relationship to the mode, while in alpha the relevant edges exert a negative influence on the mode. While there are clear differences between these modulation spectra, they tend to show broad and almost scale-free effects rather than spectrally confined differences.