Frequency-dependent functional connectivity within resting-state networks: An atlas-based MEG beamformer solution

1The brain consists of billions of interconnected neurons, forming an extremely complex system in which clusters of neurons are organised as functional units with more-or-less specific information processing capabilities (e.g. ). Yet, cognitive functions require the coordinated activity of these spatially separated units, where the oscillatory nature of neuronal activity, and phase relations between units, may provide a possible mechanism. Not only has it been shown that rhythmic activity plays an important role in perception and sensori-motor systems, as well as in higher cognitive functions, but it has also been shown that patterns of resting-state oscillatory activity in patients with neurological disorders differ from those in healthy subjects, and that these differences correlate with cognitive performance.

2Magnetoencephalography (MEG), with its high temporal resolution, can be used to characterise the (resting-state) networks formed by interacting sources of oscillatory activity. However, in many MEG studies the required estimation of functional connectivity is performed at the sensor-level, which impedes comparison with the rapidly growing literature on resting-state functional connectivity using functional Magnetic Resonance Imaging (fMRI; ). Another problem is that multiple recording sites pick up signals from a single source due to both the nature of the induced magnetic flux (see e.g. ) and volume conduction, which can lead to erroneous estimates of functional connectivity when these estimates are based on sensor-level measurements. Moreover, the signals originating from spatially separated brain areas are mixed at the sensor level, which can result in over/underestimation of synchronisation, where the exact effect is dependent on a complex interplay of modulations in source- and noise-power, source interactions, as well as relative position and orientation of the sources (e.g. ).

3These limitations have provoked research into four different directions: i) using estimates of the expected functional connectivity that is due to volume conduction/field spread without true interactions, and subtracting such estimates from the measured functional connectivity, or using such estimates to derive a statistical threshold for the measured interactions (e.g. ); ii) the development of functional connectivity estimators that are insensitive to these confounds of field spread and volume conduction; iii) development of techniques such as DCM (Dynamic Causal Modelling; ) that, given a set of hypotheses, test between different bio-physically motivated source- or network-models based on their model evidence; and iv) investigations into the utility of functional connectivity analysis at the source level. The simplest approach here is to create a source model to project the (unaveraged) sensor data into source space, although the creation of such a source montage requires the availability of averaged evoked data (see though). The construction of a source montage can also be achieved by using a combination of Independent Component Analysis (ICA) and a multivariate autoregressive (MVAR) model, or using Principle Component Analysis (PCA), ICA or MVAR in combination with an inverse estimator. Alternatively, functional networks can be directly estimated at the source level, although this efficient approach is only applicable when the source matrix is obtained though a distributed linear inversion, and when the source coupling can be described as a scalar product between the source signals.

4Various modifications of linear estimators have been used to reconstruct the time-series for a large number of locations, for a set of cortical patches (e.g. ), or for a limited set of a-priori defined Regions-of-Interest (ROIs; e.g. ), where the time-series were subsequently used for the estimation of functional connectivity.

5A drawback of the above approaches is that the spatially smooth estimates of neuronal activity that are obtained contain widespread correlations between reconstructed source elements, so that estimates of functional connectivity between such sources are likely to be, as with sensor-level analysis, erroneous and/or difficult to interpret.

6Here we propose to use a source reconstruction approach that results in sharper 3-dimensional images of neuronal activity, known as beamforming, which has recently been used to map functional connectivity across the entire brain, or to characterise interactions between a few ROIs; see also and for related approaches. We use beamforming to estimate time-series for a set of atlas-based ROIs that cover the brain, where the use of a standard brain-atlas aids the interpretation of our results, gives a robust platform for group-level statistics, and enables a straightforward comparison with results obtained using other modalities. Functional connectivity between these ROIs is then estimated, and we first demonstrate that the effects of volume conduction and biases introduced by the beamformer can be removed by using the Phase Lag Index (PLI) for the estimation of functional connectivity. Applying this approach to resting-state MEG data in healthy controls reveals clearly distinct frequency-band dependent patterns of resting-state functional connectivity. Generally, these patterns corresponded closely to patterns of relative source power, suggesting that the most active brain regions are also the ones that are most-densely connected.

7We used previously analysed MEG data from 13 healthy subjects, where they formed part of studies on Parkinson's disease for which approval was obtained from the medical ethics committee of the VU University Medical Center. In these studies oscillatory power, as well as functional connectivity and network characteristics at the sensor level, were estimated and compared between healthy controls and demented and non-demented patients with Parkinson's disease.

8All subjects gave written informed consent prior to participating. MEG data were acquired in the morning, using a 151-channel whole head MEG system (CTF Systems Inc., Port Coquitlam, Canada), situated in a magnetically shielded room (Vacuum-schmelze GmbH, Hanau, Germany). The data were sampled at 312.5 Hz, with a recording pass-band of 0–125 Hz, and a third-order software gradient was applied. Each session started with an approximately 5 minutes eyes-closed (EC) resting-state recording, followed by an approximately 5 minutes eyes-open (EO) recording. We only analysed the data recorded during the eyes-closed resting-state. Due to technical problems, 1–3 channels were discarded from the analysis (3, 3, and 7 datasets contained 148, 149, and 150 channels, respectively). For the construction of the beamformer weights, the eyes-closed data were band-pass filtered from 0.5 to 48 Hz, and after visual inspection, trials containing artefacts were removed. A time-window of, on average, 264.2 seconds (range: 175–360 s.) was used for the computation of the data covariance matrix. Broadband data were used for the estimation of the beamformer weights as this avoids overestimation of covariance between channels.

9For each subject, an anatomical MRI of the head was obtained at 1 T (Impact, Siemens, Erlangen, Germany), with an in-plane resolution of 1 mm and slice thickness of 1.5 mm. Vitamin E capsules were placed at anatomical landmarks, the pre-auricular points and the nasion, to guide co-registration with the MEG data. In the MEG setting, three head position indicator coils were placed at the same fiducial locations, and these coils were activated at the start of each MEG acquisition. Head position and orientation were computed on the basis of the magnetic fields produced by these coils. Using these two corresponding sets of fiducial markers, the MEG and MRI coordinate systems were matched. The co-registered MRI was subsequently segmented, and the outline of the scalp was used to compute a multi-sphere head model for the calculation of the lead-fields.

10Fig. 1 provides an overview of the analysis framework. The initial step is to project the MEG sensor signals in a meaningful way to a set of time-series of neuronal activation in the brain. To this end, we use a popular beamforming technique, known as Synthetic Aperture Magnetometry (SAM, ). The details of this technique are described elsewhere (, see also for a review), but we will describe its main features below.

11The beamformer output at a target location, for a given orientation of the target source, can be defined as the weighted sum of the output of all (N) signal channels, or mathematically:with V the beamformer output (source strength in nAm), W the 1xN weight vector (units: nAm/T), and B the NxT1 matrix of the magnetic field (in Tesla) at the sensor locations at all (T1) latencies. V is often referred to as a virtual electrode, and has the same temporal resolution as the recorded MEG signals.

12The weights determine the spatial filtering characteristics of the beamformer and are designed to increase the sensitivity to signals from a location of interest whilst reducing the contribution of signals from (noise) sources at different locations. The beamformer weights for a source at a location of interest are completely determined by the data covariance matrix and the forward solution (lead field) for the target source (see ):with Cj the source current covariance matrix, Cb the data covariance matrix, L the lead field1 and T the matrix transpose.

13Differences between various source reconstruction algorithms arise from the different assumptions that are made about the source current covariance matrix (see ). In the case of the beamforming approach it is assumed that all sources are linearly uncorrelated, i.e. Cj is a diagonal matrix, and that each diagonal element in Cj, corresponding to a location θ, can be related to the measured data as follows

14Eq. (3) is the crux of the beamformer algorithm. It is here that the source covariance Cj is estimated based on the data. Combination of the above three equations allows for the computation of the beamformer weights and beamformer output.

15So far we assumed that the source orientation is known. SAM performs a search for the orientation that optimises the normalised beamformer output, defined aswith Σ the sensor noise covariance matrix, Nθ the power of the projected sensor noise, and the pseudo-Z statistic for location and orientation θ.

16Alternatively, an eigendecomposition can be used to determine the optimum source orientation, or one can simply compute the beamformer weights for the two tangential orientation components for each source (or all 3 orthogonal components in the case of electroencephalography (EEG); ).

17The standard approach is to subsequently form a 3-dimensional image of source activity (see for review) that quantifies the (change) in activity for a given time-segment in the MEG recording. Our approach however is to compute beamformer weights for those voxels that contained an anatomical label (see below), and use the complete time-series for all these virtual electrodes for further analysis, notably functional connectivity and network analysis.

18Our aim here is to define a set of regions-of-interest (ROIs) that is consistent across individuals. To this end, we spatially normalised a subject's anatomical MRI, using SPM99, to a template MRI. The nonlinear transformation matrix that is necessary to perform this normalisation was stored.

19The Talairach Daemon Database (TCDB), as built in the MEG/fMRI visualisation package mri3dX (https://cubric.psych.cf.ac.uk/Documentation/mri3dX), was used to label the voxels in the template MRI. We used the labels in the database at the level of description of Brodmann areas, and removed deeper structures, resulting in a set of 68 ROIs in template space (see Appendix A for a list of the labels that were used). The number of voxels per ROI varies, and is dependent on the chosen spatial sampling — voxels with side lengths of 5 mm were used in this study. The inverse of the nonlinear transformation matrix was applied to these ROIs, to create labelled target voxels in a subject's MRI for computation of the beamformer weights.

20Once the beamformer weights were estimated, Eq. (1) was used to reconstruct the time-series for each voxel. These virtual electrode time-series exhibit a non-uniform projection of sensor noise (the weights increase with depth, but the sensor level noise remains constant throughout the volume). In order to compensate for this inherent bias, we therefore normalise each beamformer weight by its vector norm before reconstructing the time-series.

21A ROI in a subject's MRI contains several voxels, all with their own time-series (virtual electrodes). The direction of each estimated virtual electrode is arbitrary (a source pointing inwards with negative amplitude produces the same external magnetic field as a source pointing outwards with positive amplitude), hence the estimated time-series for neighbouring virtual electrodes may have opposite polarities, rendering averaging of time-series across a ROI meaningless. We therefore proceeded to compute the spectrum for each virtual electrode time-series and divided the spectrum into the 5 classical EEG bands (delta (0.5–4 Hz), theta (4–8 Hz), alpha (8–13 Hz), beta (13–30 Hz), and gamma (30–48 Hz)). For each ROI and frequency band separately, we selected the voxel with maximum power in that frequency band, and used the time-series for this voxel for further analysis, resulting in a total of 5 sets of 68 time-series (one for each frequency band). Note that this procedure was carried out for each subject independently, such that the voxels that were selected to represent the ROIs were allowed to vary across subjects, mitigating the effects of co-registration, normalisation, and modelling errors (see e.g. for a similar strategy).

22Functional interactions between sources of oscillatory activity can be captured by quantifying the phase relationship between their time-series (see for a review of coupling measures). Unfortunately, despite the assumptions underlying beamformers, the beamformer reconstructed sources may still show spurious, field spread and volume conduction related, interactions, which manifests itself as locking with zero-phase lags. To show that this is the case, and to demonstrate how this problem can be solved, we use both Phase Coherence (PC) and PLI to estimate functional connectivity between ROIs.

23The Phase Coherence quantifies the phase coupling between two signals as follows:where ΔΦ is the phase difference between the instantaneous phases for the two time-series, defined in the interval [0, 2π], tk are discrete time-steps and S is the number of samples.

24Phase Coherence captures consistent phase differences and is, unlike coherence, not influenced by the amplitude of the signals. Phase Coherence is maximal when the phase difference has a constant value, whatever the value of this phase difference is, and is therefore equally sensitive to both trivial (zero-phase) and true (zero-phase and nonzero-phase) interactions.

25In contrast, the PLI is defined as:where the phase difference is defined in the interval [− π, π] and <> denotes the mean value. The PLI is non-zero when there is an asymmetry in the distribution of the instantaneous phase differences, and therefore only quantifies non-trivial connections, at the expense of potentially discarding true interactions with zero-phase lag.

26For the computation of the functional connectivity, using software developed by one of the authors (CS; Brainwave, version 0.8.92; http://home.kpn.nl/stam7883/brainwave.html), 5 artefact-free data-segments of 4096 samples were selected from the ROI time-series after careful visual inspection.

27For each ROI we computed the mean PLI and Phase Coherence with all other areas. This is also known as the weighted degree or node strength in terms of graph theory, where individual values reflect the importance of nodes in the network, the mean across ROIs indicates the total ‘wiring-cost’, and the distribution of degrees is an important marker of network development and resilience. We then computed the mean of this quantity across trials and subjects to get group mean node strength values per ROI.

28In order to determine the significance of the empirical group mean values, we created 100 sets of null data by phase randomising (whilst maintaining the power spectra) the ROI time-series. For each realisation, this gave rise to 68 new ROI time-series per subject, which were analysed in exactly the same way as the recorded data, giving a mean PLI value per ROI. Taking the maximum mean-PLI value over ROIs on each realisation gave a null distribution of PLI corrected for multiple comparisons across the volume.

29Fig. 2 displays the mean functional connectivity for each ROI with all other ROIs. Note the differences between the maps for PLI and Phase Coherence, where Phase Coherence shows strongest functional connectivity for deeper structures, whereas the strongest PLI is found mainly for superficial areas in the occipital and parietal lobe and the superior temporal gyrus, as well as in the posterior cingulate.

30These differences can be explained by the different sensitivity of these measures to volume conduction/field spread (in the form of correlation between beamformer weights). Fig. 3 shows that there is a close correspondence between the weight correlations (Fig. 3c) and the Phase Coherence (Fig. 3b). This is confirmed in Fig. 3e, where a strong relationship between the weight correlations and Phase Coherence can be seen. This means that it is likely that any observed Phase Coherence could have been caused by the beamformer itself. Such a clear relationship between the mean PLI and the weight correlations is not observable in Fig. 3d.

31Fig. 4 displays the mean PLI for the alpha, beta and gamma bands (the adjacency matrices themselves are given in Supplementary Fig. 1), as well as the mean relative power for these frequency bands. The mean PLI for the delta and theta bands did not reach statistical significance (see Supplementary Fig. 2 for the unthresholded PLI maps and maps of relative power for the delta and theta bands; Supplementary Fig. 3 shows results for the alpha1 and alpha2 bands).

32The most strongly connected, as well as most strongly active, regions in the alpha band were the posterior cingulate, regions in visual cortex and parietal lobe, as well as in the superior and inferior temporal lobe.

33The connectivity and power maps in the beta band showed a relatively restricted pattern of highly connected and strongly activated regions in the sensorimotor cortex, extending into the inferior parietal lobe and the dorsolateral prefrontal cortex (for source power).

34In the gamma band, the most strongly connected regions were in the temporal lobe, sensorimotor cortex and the inferior frontal and parietal lobes. Dominant power was also found in these regions, but covered more of the frontal lobe and also included the visual cortex, with notably less power in the parietal lobe and posterior part of the sensorimotor cortex.

35It is clear from Fig. 4 that the patterns for PLI and power are similar, but that there are also notable exceptions where regions with high power do not have corresponding high PLI values, and vice versa (compare for example the power and PLI values for the occipital pole and parietal lobe in the gamma band). This is further illustrated in Fig. 5, where the relationship between PLI and power for the different frequency bands is shown. For all frequency bands, except the gamma band, there is a significant positive linear relationship between PLI and relative power (for delta, theta, alpha, beta and gamma bands respectively: F(1,66) = 69.98, p < 10− 11; F(1,66) = 8.99, p < 0.01; F(1,66) = 194.31, p < 10− 15; F(1,66) = 167.97, p < 10− 15; F(1,66) = 1.30, p = 0.26). However, the PLI values are not trivially related to the power values, as, for example, the relative power for the beta band and delta bands do not differ significantly (two-sample t(134) = − 1.56, p = 0.12), whereas the mean PLI values are significantly higher in the delta band than in the beta band (two-sample t(134) = 118.66, p < 10− 15). Similarly, the PLI for the theta and alpha bands do not differ significantly (two-sample t(134) = − 1.45, p = 0.15), whereas the mean power values are significantly higher in the alpha band than in the theta band (two-sample t(134) = 10.83, p < 10− 15).

36We have presented a robust method for assessing significant functional connectivity across a group of subjects that is insensitive to volume conduction and statistically well controlled. We first demonstrated that when using Phase Coherence, functional connectivity estimates, even at the source-level, are biased due to the effects of volume conduction and field spread. In contrast, functional connectivity estimates based on PLI are not affected by these biases. Subsequently, we used the method to show significantly higher than chance interactions between regions in three of the five frequency bands studied, revealing distinct frequency-band dependent patterns of functional connectivity across the brain.

37Strong connectivity was observed in the visual cortex and in the parietal and temporal lobes, consistent with the findings by. However, we could not find any direct evidence in the existing literature for the strong resting-state functional connectivity in the posterior cingulate for the alpha band.

38Similarly, the dominant alpha power in the posterior part of the brain is consistent with the established patterns of eyes-closed resting-state alpha activity (e.g. ). Additionally, several studies, using different methodologies and modalities, have reported activations similar to ours for the visual cortex, parietal and temporal lobe as well as the posterior cingulate.

39We found strong resting-state beta-band functional connectivity in the sensorimotor cortex, in agreement with a previous MEG resting-state functional connectivity study. Moreover, the large body of literature on beta-band synchrony in sensorimotor systems involved in co-ordinated movement and posture (e.g. ), suggests that these systems are also connected in this frequency band during rest. Further evidence for this hypothesis comes from single-pulse TMS studies that have shown that TMS synchronises the phase of beta-oscillations, and that stimulating one region of the sensorimotor system can induce activity in another part of this system (e.g. ).

40Our finding of strong resting-state beta power in sensori-motor cortex and inferior parietal lobe is consistent with studies that have revealed event-related beta-power reductions in these regions following somatosensory stimulation or movement (e.g. ), or that showed that the beta-band is the natural frequency of these circuits. Our observed beta band power in the dorsolateral prefrontal cortex can tentatively be related to the reorienting of attention.

41We had naively expected that gamma band PLI (like alpha band PLI) would predominate in the visual cortex. We found little literature to corroborate the resting-state functional connectivity in the temporal lobe, sensorimotor cortex and the inferior frontal and parietal lobe for the gamma band. Recent invasive electrode and non-invasive studies on pitch perception have pointed to a functional role for the gamma band in the region of human auditory cortex. It would make sense that similar functional units are also to some degree engaged during the resting-state.

42In terms of source power however, the observed dominant power in the frontal and temporal lobe is in agreement with previous EEG findings. Our observations are also in agreement with reports of both resting-state measures of gamma power in the primate and event-related changes in gamma power in visual cortex following visual stimulation (e.g. ).

43We found a positive relationship between the patterns of source power and functional connectivity, using a connectivity measure that is not sensitive to power or volume conduction. The fact that we found this relationship for all frequency bands (except the gamma band) suggests that there may be a general mechanism that explains this relationship. The simplest explanation is that a positive relationship between source power and connectivity was introduced by our analysis approach, since the functional connectivity measure that we used relies on accurate estimation of phase differences, and therefore on a high-enough SNR. Consequently, functional connectivity between regions with low source power may be missed or underestimated, introducing a bias towards a positive correlation between power and connectivity. However, we did find regions with strong connections despite low source power (e.g. parietal lobe for the gamma band, Fig. 4), hence it is unlikely that the positive relationship between source power and connectivity can be fully explained by such a methodological bias. One possible physiological mechanism that could provide an alternative explanation has been identified through modelling studies, which have shown that increased mean spiking activity within two connected neuronal populations leads to increased intra- and interregional phase locking, or in other words, to simultaneous increases in power and functional connectivity. This coupling between functional connectivity and mean activity was mainly achieved through a reduction in mean membrane integration times as activity increased, which then introduced a bias towards synchronous firing. Recent modelling work has also shown an amplitude dependency of interregional phase locking, albeit that this dependency was more complex.

44A frequent criticism of sensor-level connectivity/network analysis is that the results are biased by the effects of volume conduction/field spread. Our results show that this problem is not completely solved by going to source-space (Fig. 3), as it reveals itself, in our case, in the form of correlation between beamformer weights. As a solution, we used the PLI, which quantifies functional interactions that are not caused by volume conduction or common sources. The main reason for going to the source-level, in combination with an atlas-based analysis approach, is therefore that it provides a general framework that allows for a direction anatomical interpretation of MEG data, as well as a direct comparison with (functional) connectivity and network studies based on anatomical MRI (Diffusion Tensor Imaging, Voxel-based Morphometry) and functional MRI (e.g. ). We envisage that the understanding of resting-state networks will be much enhanced by such a combination of different modalities. Similarly, our approach enables the integration and direct comparison of data recorded with different MEG systems. Although we have focussed here on the eyes-closed resting-state, our approach can also be applied to compare patterns of oscillatory activity and functional networks for different cognitive states.

45One can consider the choice of using a set of atlas-based ROIs as a compromise between all-to-all connectivity estimates and methods based on a-priori selection of a small number of regions. The use of atlas-based ROIs does certainly compromise the potential spatial resolution, but could be the most efficient level of description given the inter-individual variability. Similarly, the use of predefined anatomical regions does mean that one samples from all sources rather than just those with highest source power (although within each ROI we did select the voxel with highest source power), thereby avoiding the danger that weakly activated, but strongly interacting, sources are missed (compare for example the connectivity and power in the parietal lobe for the gamma band (Fig. 4)).

46We chose to reconstruct the resting-state sources, and their time-series, using a beamforming approach, as this approach does not suffer from the problems associated with linear inverse solutions (widespread correlations between reconstructed source elements and problematic interpretation of reconstructed images of source power).

47A potential limitation of the beamformer-based approach is that the beamformer weights are based on source power and interacting sources with small amplitude could therefore be missed (but see Fig. 4, gamma band power and connectivity). Moreover, the beamformer weights are designed so that the sensitivity to signals from a location of interest is increased, whilst reducing the contribution of signals from (noise) sources at different locations. As a consequence the contribution from sources that are perfectly linearly correlated is cancelled or underestimated (which might in particular be problematic for weak long-distance interactions). At first sight, this property of beamforming seems therefore at odds with our aim to estimate interactions between sources. However, there are several reasons why beamforming can still be applied in studies on functional connectivity: i) it is important to stress that only (zero-lag) linearly correlated sources are problematic for the beamformer. Non-linear measures of functional connectivity should therefore be less affected; ii) even for linearly correlated sources, there is a remarkable tolerance to deviations from the uncorrelatedness-assumption, such that source interdependencies can be accurately reconstructed even when these sources are correlated for 30–40% of the analysis time-window; iii) anatomical and electrophysiological data suggest that the uncorrelatedness-assumption is plausible; and iv) modified beamformer approaches are available to deal with those rare cases where sources with strong linear interactions are present. In addition, the PLI and beamformer blindspots sit comfortably together — the beamformer will potentially mis-localise sources with zero-lag correlation but by using the PLI we ignore these effects.

48Other methods, such as minimum norm based approaches, have no such constraints on source correlation. However, this flexibility comes at the price of poorer noise rejection capability. The work of Ghuman et al. highlights the problems that occur when combining a minimum norm approach with a phase-locking (rather than phase-lagging) approach, in that one has to rely on the subtraction of empty room data in order to try to remove the large number of spurious interactions. A particular worry here is that minimum norm based source reconstruction approaches inherently model all the data and therefore project artefacts (such as heartbeat or from external noise sources) into the source space. In contrast, beamforming approaches only localise those components in the data that match the lead field for a source at a particular location, giving beamformers the ability to reject artefacts (e.g. ). Importantly, these projected artefacts lead to further spurious connectivity estimates when using methods that are sensitive to zero-phase interactions (i.e. methods based on phase-locking, rather than PLI).

49Beamformer implementations exist that use the matrix of connections or higher order statistics, rather than the data covariance matrix, for the computation of the beamformer weights (Eq. (2)). One could also consider the direct replacement of the data covariance matrix by the PLI adjacency matrix when computing the beamformer weights, thereby avoiding problems related to volume conduction and zero-lag phase relations between sources. This may be problematic though, since the PLI is independent of signal amplitude. However, further research could show whether the use of the imaginary coherence, which also minimises the influence of volume conduction, proves to be more fruitful in this type of approach.

50We chose to parcellate the source space on the basis of the Talairach Daemon Database, although alternative atlases are available (e.g. ). Instead, one could also use a parcellation-scheme that is based on the source-sensor geometry in order to obtain a set of maximally independent patches, or perform parcellation in the native MRIs. This latter option would in particular be preferably for patients who's MRI match poorly to a template MRI, for example due to atrophy. In addition, it remains an open question how best to deal with ROIs of unequal size. One could argue that all ROIs should contain an equal number of voxels so that estimates of interdependencies are not affected by differences in ROI size. In contrast, perhaps the size of a ROI should reflect the variations in sensitivity of MEG to neuronal activity in different regions of the brain. Overall, these disadvantages of using a standard brain-atlas are outweighed by the important advantage that its use enables a more direct comparison between data from different modalities.

51Here, we used the beamformer formulation for point sources, i.e. the lead fields for equivalent current dipoles were used for the weight computations. In order to take into account the spatial extent of the ROIs for which representative time-series are estimated, one could use the lead fields for spatially extended sources, use a set of basis functions, or use Singular Value Decomposition (SVD) to define a lower-dimensional representation for each ROI. Similarly, a combination of Signal Space Separation (SSS) and beamforming allows for the estimation of time-series for pre-defined spherical ROIs, although with this approach one would lose the advantages of using a standard brain atlas. Instead, the direct replacement of the SAM beamformer by the recently developed SSS-beamformer would fit more naturally in our proposed analysis framework.

52Recently, have proposed a nulling-beamformer, which removes potential cross-talk between ROIs by incorporating additional (nulling-) constraints for these ROIs into the beamformer design. It is not clear though how the reduction in the degrees of freedom due to the use of many ROIs in a standard atlas (i.e. many nulling-constraints would be needed) would affect both the ability to reject noise and the accuracy of the reconstructed time-series.

53We used a multi-sphere head model for the computation of the lead fields, which provides an accurate approximation of the volume conductor for MEG, although more complex numerical models may provide increased accuracy in certain situations. Inaccuracies in the volume conductor model may lead to underestimation of source power. However, given that PLI is independent of amplitude, it is unlikely that connectivity biases were introduced, even if there were regions for which source power was underestimated.

54The length of the data covariance window that was used for the weight computations was determined by the amount of data that was available from each recording session, the duration of which was set at what was common practice for a resting-state MEG session at the time of recording (2003/2004). Simulation studies have shown that the required number of samples for beamforming depends, among other factors, on the frequency-band of interest, sampling rate and source power. The results by Brookes et al. suggest that even for our worst case scenario (small band width of 3.5 Hz for the delta band and smallest co-variance window of 175 s) the errors in the estimation of the data covariance, and therefore in source power estimates, were minimal (less than 10% underestimation in source power).

55It is feasible to apply noise-regularisation during the computation of the beamformer weights (not used in this study), which would lead to increased signal-to-noise ratio (SNR) for the estimated time-series, and would also mitigate against the effects of using a limited set of (atlas-based) voxels in template space (activation could be missed due to co-registration and normalisation errors — see e.g. for a similar strategy). However, regularisation comes at the expense of a decreased ability to distinguish spatially separate sources. The optimal trade-off between this temporal and spatial accuracy for studies that aim to compute the topography and topology of functional networks, has yet to be determined for empirical data.

56To deal with the issue of arbitrary sign for the orientation of the source at each voxel, which makes straightforward averaging of source waveforms across a ROI impossible, we selected the voxel with maximum source power within a ROI. It has been shown previously that the time-series estimated at local maxima best describe the underlying source activity. For datasets with large artefacts, these artefacts could leak into the reconstructed time-series and potentially bias the selection of the voxel with maximum power. Alternative approaches include: i) performing a check on (and adjustment of) the polarity of the time-series of neighbouring voxels before averaging time-series across a ROI. This assumes that the source orientation varies smoothly when moving through the source space; ii) using SVD to find the eigenvectors (time-series) that best represent the time-series of the ROI (see e.g. ); and iii) using the time-series that most strongly correlates with the time-series for the other voxels in the ROI.

57A consequence of selecting the voxels with maximum power to represent the ROIs is that, particularly when source activation spreads over multiple ROIs, the peak voxels for neighbouring ROIs can be close together (Supplementary Fig. 5), i.e. such voxels share (almost) the same signal. For such cases, the source reconstruction approach could not unambiguously determine whether the activity is coming from one or the other ROI (or both). Importantly, PLI is insensitive to the spurious zero-lag interactions that could exist between voxels that are close together, hence the voxel selection will not lead to overestimates of (local) connectivity.

58Alternatively, one could, for each voxel, extract the power modulations of the time-series using the Hilbert transform and average these envelopes across a ROI. This approach has already demonstrated interesting relationships between functional networks constructed on the basis of fluctuations in MEG band-limited power and those based on low-frequency modulations in BOLD fMRI time-series. A limitation of this approach is that differences in ROI size might result in biases due to differences in SNR.

59In this work we contrasted PLI with Phase Coherence as measures of functional connectivity between ROIs. As expected, the Phase Coherence suffers from spurious correlations between ROIs due to correlations between their beamformer weights (Fig. 3). PLI on the other hand is a conservative measure that is insensitive to such spurious interactions, albeit at the expense that true zero-lag correlations are also missed. Zero-lag correlations are most-likely to be short-range (e.g. ), but see, and functional networks constructed on the basis of PLI could therefore have a topology for which the clustering and/or modularity are underestimated. Note that this is a general issue for all MEG/EEG connectivity studies, and not a problem that is specific for the proposed analysis framework. In fact, when improved connectivity estimators become available, then they can easily be incorporated in our analysis framework. have recently described a promising modification of PLI with reduced underestimation of connectivity between sources with small-lag interactions, as well as a reduced estimator bias, although this comes at the expense of introducing an arbitrary bias favouring large phase differences and mixing of the estimation of consistency of phase differences with the estimation of the magnitude of the phase difference. An alternative approach that would avoid these issues, and at the same time reduce the effects of noise on PLI-based connectivity estimates (particularly when interactions occur with almost zero-phase lag) would be to ignore phase differences within a small window around zero and around ± π. However, although the use of such an offset could potentially be useful when considering small numbers of trials, noise-related counts of positive and negative phase lags (regardless of magnitude) should not bias the statistics.

60We should note that the time window for the connectivity analysis was fixed here at 13.1 s, and that we analysed 5 of such artefact-free data-segments. Clearly, the choice of time window determines not only the expected period of stationarity of the interactions, but also the bandwidth over which interactions can occur. Further investigations may reveal whether different forms of interactions are highlighted for different choices of time-frequency parameters. Indeed, empirical observations have revealed that relatively short epochs (~ 10 s) are quasi-stationary, and that such short epochs can already be representative of a subject's cognitive fingerprint (see ). Moreover, the work by has shown that functional networks at these time-scales are stable yet dynamic (on longer time-scales the functional networks resemble the underlying (static) anatomical network, whereas on shorter timescales the functional networks are highly variable). Additionally, it has been shown that using 5 epochs of 10 s resting-state data results in fairly good to good levels of test–retest reliabilities, depending on the frequency band and metrics analysed. It should be noted that Lin et al. analysed magnetometer data at the sensor level and computed complex network metrics, whereas in our study third-order gradiometer data was analysed using a beamformer approach, both resulting in improved SNR, and that we used simple metrics to characterise the functional networks; these factors all contribute to an increased stability of the estimated resting-state network properties. Finally, using the same number of epochs as we routinely use in our clinical studies, which has proven to give stable estimates of resting-state activity/network parameters (e.g. ), renders our developed methodology and the results from the current study directly relevant to our clinical work.

61We looked at the average overall connectivity for a ROI, not, as is typical in fMRI, specifically at connections between certain ROIs. For example, the average connectivity for a ROI could be relatively low, but the ROI could still have strong connections with (only) some ROIs in, let's say, the default mode network. A more in depth comparison with the fMRI literature would be an interesting topic for future research. Indeed, recent work has shown the correspondence between beamformer estimated power envelope correlations and fMRI defined resting-state networks. It will be interesting to see the make-up of the fMRI defined networks in terms of the relative contributions of electrical interactions over dimensions of power, phase and frequency.

62For the creation of the phase-randomised surrogates we did not take into account any jumps at the boundaries, which could have introduced high-frequency artefacts in our surrogate data. However, any biases introduced in our statistics will have been minimal, as we subsequently filtered the surrogate data to relatively low frequencies (maximal 48 Hz for the gamma band).

63We have introduced a general MEG analysis framework for the reconstruction of frequency-dependent profiles of source power and functional connectivity, which is robust to artefactual connectivity estimates caused by volume conduction, due to the use of the PLI; robust to artefactual connectivity estimates caused by physiological (e.g. heartbeat) and environmental (e.g. power-line) noise, due to the use of the beamformer approach; robust to co-registration errors, due to the use of unconstrained source orientation at each voxel; and robust to modelling errors (e.g. introduced by co-registration errors or by ignoring source extent), due to the use of coarse spatial sampling (defined by the ROIs). Finally, the results are based on non-parametric statistics with few underlying assumptions.

64The analysis framework contains two important elements: i) activity is reconstructed for an atlas-based set of ROIs in order to facilitate interpretation and comparison with results obtained with other modalities; and ii) effects of volume conduction/field spread on estimated interactions between ROIs are removed using PLI, a measure that is insensitive to these effects. Using this framework we have revealed distinct frequency-dependent patterns of source power and source interactions. We envisage that this approach will be used to further elucidate the patterns of resting-state activity in health and disease.

65The following are the supplementary materials related to this article.