Quantifying the Test-Retest Reliability of Magnetoencephalography Resting-State Functional Connectivity

Abstract

The coordinated activity of the resting-state brain can be evaluated with magnetoencephalography (MEG) for distinct brain rhythms by performing source reconstruction to estimate the activities of target brain regions and employing one of the many existent functional connectivity (FC) algorithms. Although this procedure has been applied in a great amount of studies both with healthy and pathological populations, the reliability of such FC estimates is unknown, and this impairs the use of resting-state MEG FC at the individual level. In this study, the test-retest reliability of MEG resting FC was evaluated by exploring both within- and between-subject variability in FC in 16 healthy subjects who underwent three resting-state MEG scans. FC was computed after beamforming source reconstruction with four popular FC metrics: phase-locking value (PLV), phase lag index (PLI), direct envelope correlation (d-ecor), and envelope correlation with leakage correction (lc-ecor). Then, test-restest reliability and within- and between-subject agreement were evaluated with the intraclass correlation coefficient (ICC) and Kendall's W, respectively. Reliability was found to depend on the FC metric, the frequency band, and the specific link. As a general trend, greater test-retest reliability was found for PLV in theta to gamma, and for lc-ecor and d-ecor in beta. Further inspection of the ICC distribution revealed that volume conduction effects could be contributing to high ICC in PLV and d-ecor. In addition, stronger links were found to be more reliable. Overall, this encourages the further use of resting-state MEG FC for individual-level studies, especially with PLV or envelope correlation metrics.

Introduction

Since Biswal et al. (1995) observed a highly correlated brain activity between bilateral motor cortices, increasing focus has been shifted to the functioning of the resting-state brain (Deco et al., 2011; van den Heuvel and Hulshoff Pol, 2010). In fact, the resting-state brain activity is highly organized, and separate brain regions have been repeatedly shown to be spontaneously coordinated, accounting for sensory and cognitive functions, such as visual or motor processes, attention, or executive control (Rosazza and Minati, 2011).

This spontaneous brain coordination is usually characterized with functional connectivity (FC), which evaluates statistical dependencies between the activities of distinct brain regions (Pereda et al., 2005). FC can be seen as a mechanism of information transfer between brain areas, and resting FC has been proven to be indicative of an individual's brain functioning in a variety of studies (Baldassarre et al., 2012; Rosazza and Minati, 2011; Yamashita et al., 2015).

Resting-state FC has been predominantly examined with functional magnetic resonance imaging (fMRI) as correlations between blood oxygenation level-dependent (BOLD) signals (van den Heuvel et al., 2009). Although fMRI is of great use in the exploration of human brain function, slowly fluctuating BOLD signals constitute only an indirect measure of brain activity. Therefore, over the past few years, attention has been drawn to the electrophysiological basis of FC. In fact, although magnetoencephalography (MEG) is often overshadowed by the high spatial resolution of fMRI, MEG provides an exceptional insight into the time–frequency dynamics of brain activity. It is indeed a direct measure of neuronal firing and offers an excellent temporal resolution. Resting-state networks (RSNs) have been successfully extracted with MEG over the past years using source-space FC (Brookes et al., 2011a,b; de Pasquale et al., 2010; Hall et al., 2013; Hipp et al., 2012; Luckhoo et al., 2012; Wens et al., 2014). Moreover, resting-state MEG FC has been proven sensible to detect abnormal brain functioning in a variety of diseases, including Alzheimer's disease (Engels et al., 2015; López et al., 2014), multiple sclerosis (Tewarie et al., 2015), and schizophrenia (Bowyer et al., 2015).

Overall, this indicates that FC can be robustly assessed with MEG. However, FC estimates could be affected by measurement noise and the subjects' cognitive or emotional state, and it is unclear to what extent FC estimates are representative of an individual. This prevents the potential use of MEG FC in the characterization of individual subjects or in the proper classification between healthy and pathological populations. Therefore, the reliability of FC estimates needs to be carefully examined. To date, only a handful of studies have tackled test-retest reliability of MEG/electroencephalography (EEG) FC (Deuker et al., 2009; Hardmeier et al., 2014; Jin et al., 2011), and they were restricted to sensor space, therefore failing to provide the spatial resolution necessary for the evaluation of whole-brain FC. To our knowledge, no study has quantified the test-retest reliability of source-space MEG FC estimates.

In this work, we investigated the test-retest reliability of MEG FC. For that, we computed whole-brain resting state FC for 16 subjects, who were scanned thrice with a 1-week test-retest interval. For each subject and session, source reconstruction of the sensor space data was performed with beamforming, and FC was computed between 66 brain regions. FC was evaluated with four phase and envelope synchronization metrics, both with and without correction for volume conduction: envelope correlation (d-ecor), envelope correlation with leakage correction (lc-ecor), phase locking value (PLV), and phase lag index (PLI). The test-retest reliability of the FC estimates was then assessed with the intraclass correlation coefficient (ICC). The effect of the volume conduction on the variability in ICC values across FC metrics and frequency bands was assessed through a linear mixed model (LMM), including the ICC of the beamformer weights correlation. In addition, within- and between-subject agreement of FC maps was computed with Kendall's W.

Materials and Methods

Subjects and MEG acquisition

Sixteen healthy subjects (age 30.4 ± 5.8, 10 women and 6 men) underwent three resting-state MEG sessions (two conditions: eyes open and eyes closed) over 3 consecutive weeks. For each participant, scans were scheduled at the same day of the week and same time of the day (8 subjects in the morning, 8 subjects in the afternoon). The study was approved by the local Ethics Committee. MEG resting-state recordings were acquired at the Center for Biomedical Technology (Madrid, Spain) using an Elekta Vectorview system with 306 sensors (102 magnetometers and 204 planar gradiometers), inside a magnetically shielded room (Vacuumschmelze GmbH). During the measurements, participants sat 4 min with their eyes open staring at a white cross at 2 m distance and 4 min with their eyes closed. They were instructed to remain calm and move as little as possible. A Fastrak Polhemus system digitized each subject's head and four coils were attached to the forehead and mastoids, so that the head position with respect to the MEG helmet was continuously determined. The activity in electrooculogram channels was also recorded to keep track of ocular artifacts. Signals were sampled at 1000 Hz with an online filter of bandwidth 0.1–300 Hz. Maxfilter software (version 2.2., Elekta Neuromag) was used to remove external noise with the temporal extension of the signal space separation (tsss) method with movement compensation (Taulu and Simola, 2006).

Source reconstruction

The continuous resting-state sensor data were scanned for jump, muscle, and ocular artifacts with FieldTrip (Oostenveld et al., 2011) artifact detection procedures, and nonoverlapping artifact-free 6-sec epochs were located in the continuous resting-state recordings. This yielded 24.3 ± 5.7 clean epochs for the eyes open condition and 26.8 ± 4.7 clean epochs for the eyes closed condition. Then, MEG data were filtered into delta (2–4 Hz), theta (4–8 Hz), alpha (8–13 Hz), beta (13–30 Hz), and gamma (30–45 Hz) frequency bands with a finite impulse response filter of order 1000 and downsampled to a 250 Hz sampling rate.

For each subject, source models were built by spreading dipoles over the cortical surface, with 6 mm spacing. This was achieved by segmenting each subject's T1-weighted MRI with Freesurfer (version 5.1.0; Fischl et al., 2002; Ségonne et al., 2007) and using MNE software (Gramfort et al., 2014) to downsample the corresponding surfaces and realigning to the MEG coordinate system. The forward model was then solved with a three-shell boundary element method, using the inner skull, outer skull, and scalp surfaces extracted from the subject's MRI with NFT software (Acar and Makeig, 2010) and computing leadfield values with MNE.

Source time series were reconstructed with a linearly constrained minimum variance beamformer (Van Veen et al., 1997). For each subject, session, source, and frequency band beamforming filters were built with the magnetometer covariance matrix averaged over trials and a 5% regularization factor. These filters were then applied to the continuous resting-state data, yielding a 4-min time series per source location. To avoid edge artifacts, source time series were only split into clean trials after the FC computation.

FC computation

To compute FC between source time series, four different metrics were used: envelope correlation (d-ecor), envelope correlation with leakage correction (lc-ecor), PLV, and PLI. First, for all source time series \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$x \left( t \right)$$ \end{document} , the phase \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \varphi _x} \left( t \right)$$ \end{document} and the amplitude \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${A_x} \left( t \right)$$ \end{document} were extracted as argument and modulus of the complex analytical signal (Niso et al., 2013; Rosenblum et al., 2001): \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}\zeta \left( t \right) = x \left( t \right) + i{x_H} \left( t \right) = {A_x} \left( t \right) {e^{i{ \varphi _x} \left( t \right) }}. \tag{1}\end{align*} \end{document}

Phase synchronization was computed by evaluating the distribution of phase differences extracted from two source time series \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$x \left( t \right)$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$y \left( t \right)$$ \end{document} , sampled simultaneously at \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \left\{ {{x_k}} \right\} _{k = 1 , 2 \ldots , n}}$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \left\{ {{y_k}} \right\} _{k = 1 , 2 \ldots , n}}$$ \end{document} . First, the PLV measures how much this phase difference differs from a uniform distribution (Lachaux et al., 1999): \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}PLV_{xy} = \frac { 1 } {n} \left| { \sum \limits_{ k = 1 }^n {\rm exp} \left({i\left({{\varphi _{x, k}} - { \varphi_{y, k}}} \right)} \right)} \right|. \tag { 2 } \end{align*} \end{document}

\documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$PL{V_{xy}}$$ \end{document} varies between 1 for a constant phase difference and 0 for a randomly distributed phase difference (within \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$]{-\pi , \pi}]$$ \end{document} ). Second, an alternative phase synchronization measure is the PLI, which has the advantage of escaping volume conduction bias, since it is insensitive to zero and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\pi$$ \end{document} phase differences (Stam et al., 2007): \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}PLI_{xy} = \frac { 1 } { n } \left| {\sum \limits_ { k = 1 } ^n { \rm { sign } } \left[ { { \rm{ sin } } \left( { { \varphi _ { x , k } } - { \varphi _ { y , k }} } \right) } \right] } \right| . \tag { 3 } \end{align*} \end{document}

\documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$PL{I_{xy}}$$ \end{document} also varies between 0 (no synchronization or synchronization centered around 0 or \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\pm \pi$$ \end{document} ) and 1 (maximal synchronization). It quantifies the asymmetry in the phase difference over the upper and lower half-unit circle ( \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$] {0 , \pi } ]$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$] { - \pi , 0} ]$$ \end{document} , respectively). For instance, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$PL{I_{xy}} = 1$$ \end{document} if \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \varphi _x} - { \varphi _y}$$ \end{document} remains in the \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\left[ {0 , \pi } \right]$$ \end{document} interval.

Envelope synchronization can be directly assessed with the correlation coefficient between \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${A_x} \left( t \right)$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${A_y} \left( t \right)$$ \end{document} , which captures simultaneous power increases/decreases between \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$x \left( t \right)$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$y \left( t \right)$$ \end{document} . \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}{\rm d-ecor}_{xy} = \left| {corr \left( {{A_x} \left( t \right) , {A_y} \left( t \right) } \right) } \right| . \tag{4}\end{align*} \end{document}

Envelopes were smoothed with a 0.5 sec moving average filter. Envelope correlations are particularly relevant in MEG RSNs, since they seem to be closely related to fMRI BOLD fluctuations and have successfully reconstructed fMRI RSNs (Brookes et al., 2011a,b; Maldjian et al., 2014). As for the PLV, d-ecor estimates could be affected by volume conduction. The lc-ecor measure accounts for this effect by performing an orthogonalization of the source time series (Brookes et al., 2012) before envelope and correlation computation: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}{y_R} \left( t \right) = y \left( t \right) - \beta x ( t ) , \tag{5}\end{align*} \end{document}

To avoid edge artifacts, envelopes and phases were estimated in the whole 4-min time series, but correlations PLV and PLI were only computed for the previously determined 6-sec clean trials. FC values were averaged over trials to produce a final FC estimate per pair of source locations, condition, frequency band, subject, and session.

Reliability definition and estimation

Reliability is usually considered a relative quantity, which can be defined as follows (Weir, 2005): \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}reliability = { \frac { { \sigma _r } ^2 } { { \sigma _r } ^2 + { \sigma _w } ^2 } } , \tag { 6 } \end{align*} \end{document}

where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \sigma _r}^2$$ \end{document} represents the true variance of the magnitude of interest and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \sigma _w}^2$$ \end{document} the error variance. Reliability ranges therefore from 0 (no reliability) to 1 (perfect reliability). In this study, we focus on test-retest reliability, which is estimated by carrying out various measurement sessions for each subject. Let us call \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${z_{ij}}$$ \end{document} the measured (or observed) score for subject \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$i$$ \end{document} in session \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$j$$ \end{document} . For instance, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${z_{ij}}$$ \end{document} can designate an FC estimate between two sources for subject \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$i$$ \end{document} and session \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$j$$ \end{document} . Then, reliability is defined as follows (Bartlett and Frost, 2008; Scholtes et al., 2011; Weir, 2005): \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}reliability = { \frac { { \sigma _ { z , i } } ^2 } { { \sigma _ { z , i } } ^2 + { \sigma _w } ^2 } } , \tag { 7 } \end{align*} \end{document}

We assume thereby that a true quantity (or, in our case, coupling strength) \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$Z_{true, i}$$ \end{document} exists for each subject \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$i$$ \end{document} , and we call error any deviation between \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${z_{i , j}}$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${z_{true , i}}$$ \end{document} .

The intraclass correlation (or ICC, Shrout and Fleiss, 1979) is the most commonly employed reliability estimate when dealing with interval and ratio values. It evaluates between-subject variance \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \sigma _{z , i}}^2$$ \end{document} and error variance \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \sigma _w}^2$$ \end{document} from a measurement sample, and then uses Equation (7) to quantify reliability. Various ICC types exist, depending on the model used to estimate \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \sigma _{z , i}}^2$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \sigma _w}^2$$ \end{document} . In this text, we focus on the ICC type 1–1 described in McGraw and Wong (1996); Shrout and Fleiss (1979); and Weir (2005), which applies to our study of the test-retest reliability of MEG estimates. This model is based on a one-way random effects model, in which the ordering of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$j$$ \end{document} is assumed to be irrelevant. This is equivalent to assuming that there is no systematic error between observations or subjects do not change consistently across MEG sessions.

The type 1–1 ICC is: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}ICC = { \frac { MSR - MSW } { MSR + 2MSW } } , \tag { 9 } \end{align*} \end{document}

where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$n = 16$$ \end{document} is the number of subjects, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$k = 3$$ \end{document} is the number of MEG sessions per subject, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \overline x}_{i.}$$ \end{document} is the session average for subject \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$i$$ \end{document} , and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \overline x}{..}$$ \end{document} is the subject and session average. We note that, although reliability is a nonnegative magnitude, in an experiment with no true test-retest reliability, ICC estimates fluctuate around zero. It is then possible to obtain \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$MSR < MSW$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$ICC < 0$$ \end{document} .

Reliability of whole-brain FC

The test-retest FC reliability was first assessed for all links connecting 66 cortical regions of interest (ROIs). These ROIs were defined in the individual's T1 volume, using Freesurfer software (version 5.1.0) and all 33 cortical regions per hemisphere of the Desikan–Killiany atlas (Desikan et al., 2006). For each ROI, a representative time series was chosen: it corresponded to the source time series that presented the highest average correlation with the other sources within the ROI. FC between ROIs was then computed using all four FC metrics described in the FC Computation section, yielding a 66 × 66 connectivity matrix per FC metric and frequency band. ICC was then computed for each link, FC metric, and frequency band following Equation (9).

The ICC quantifies the test-retest reliability of our FC estimates, but cannot separate which fraction of the reliability is due to the underlying true FC strength and which fraction is due to artifacts in the FC estimate, such as volume conduction. In fact, we can expect that volume conduction presents lower within- than between-subject variability, and therefore may contribute to higher test-retest reliability in FC metrics such as d-ecor or PLV, which are sensible to volume conduction. As a further step, we attempted therefore to elucidate the effect of volume conduction in the ICC. We note that this is a tricky task, since there are no direct measures of volume conduction for a given MEG recording. The correlation between beamformer weights is, however, indicative of source leakage (Hillebrand et al., 2012). In fact, if source leakage exists between two sources, their beamfomer weights are similar (or highly correlated) and their reconstructed time series are highly coupled. In the opposite case, if the beamfomer weight correlation between two sources is low, it is unlikely that high FC between these two sources result from volume conduction. We accordingly computed 66 × 66 beamfomer weight correlation matrices for each frequency band, subject, and session as follows: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}BW{C_{pq}} = corr \left( {{w_p} , {w_q}} \right) , \tag{12}\end{align*} \end{document}

where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$corr \left( \ldots \right)$$ \end{document} denotes Pearson correlation and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${w_p}$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${w_q}$$ \end{document} are \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$1 \times Nchannels$$ \end{document} beamformer weights for the representative sources in ROIs \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$p$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$q$$ \end{document} , respectively. Next, the test-retest reliability of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$BWC$$ \end{document} , which will be referred to as \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$ICCBWC$$ \end{document} , was computed for each link and frequency band with the ICC.

To assess the effect of beamformer weight correlation (BWC) in the ICC of FC values, the following LMM was employed: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}ICCFC = X \beta + {z_0} + { \varepsilon _0} , \tag{13}\end{align*} \end{document}

where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$X = [iscorrd, \ iscorrleak, \ isPLV, \ isPLI , \ iscorrd \cdot ICCBWC , \ iscorrleak \cdot ICCBWC , \ isPLV \cdot ICCBWC , \ isPLI \cdot \ ICCBWC ]$$ \end{document} is the fixed effect design matrix, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\beta = { \left[ {{a_{d - ecor}}{a_{lc - ecor}}{a_{PLV}}{a_{PLI}}{b_{d - ecor}}{b_{lc - ecor}}{b_{PLV}}{b_{PLI}}} \right] ^T}$$ \end{document} is the vector of fixed effect coefficients, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${z_0}$$ \end{document} is the random effect coefficient, and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \varepsilon _0}$$ \end{document} is the error term. The coefficients \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$b$$ \end{document} represent the strength of the effect of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$ICCBWC$$ \end{document} in \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$ICCFC$$ \end{document} for the corresponding FC metric, and the coefficients \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$a$$ \end{document} represent the strength of all other fixed effects in ICCFC different than ICCBWC, for the corresponding FC metric. The model was fit (for each frequency band separately) using Matlab's fitlmematrix function and a restricted maximum likelihood estimation method.

Finally, to explore a possible relation between ICCFC and FC strength, Spearman correlations were computed across links between both variables, for each FC metric and frequency band separately. Furthermore, partial correlation between ICCFC and FC, controlling for the effect of ICCBWC, was computed to account for volume conduction effects. Bonferroni correction was applied to control for the multiple comparisons problem.

Within- and between-subject agreement of whole-brain FC maps

Reliability, as computed above, depends both on within- and between-subject variability. To inspect both effects separately, we evaluated the within- and between-subject agreement of FC maps. For a given condition, FC metric, and frequency band, FC maps were sorted across links (ROI pairs) \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$l$$ \end{document} and transformed into ranks \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${r_{ijl}}$$ \end{document} ( \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$1 \le {r_{ijl}} \le {N_{links}} , Nlinks = 66*65 / 2 = 2145$$ \end{document} ) for every subject \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$i$$ \end{document} and session \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$j$$ \end{document} separately. Then, the agreement between two sessions \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${j_1}$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${j_2}$$ \end{document} (or within-subject agreement) was assessed with Kendall's W (Legendre, 2005): \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} { W_w } = \frac { 12 \sum \nolimits_ { l = 1 } ^ { { N_ { links } } } \left( { A_ { il } } - { { { \overline A } _i } } \right) ^2 } { { 2^3 } \left( { N_ { links } } ^ { 3 } - { N_ { links } } \right) , } \tag { 14 } \end{align*} \end{document}

Reliability of RSNs

To further investigate the spatial characteristics of the reliability of resting-state FC, we focused on seven well-known RSNs: visual, sensorimotor, auditory, DMN (default mode network), left FP (frontoparietal), right FP, and frontoinsula networks. Each RSN was defined as a set of 2 or 3 nodes, following their spatial extent described previously in the literature (Laird et al., 2011; Smith et al., 2009). The nodes were defined as a 1 cm-radius sphere around a central coordinate defined in MNI space (Table 1), which was then transformed into subject space with a homogeneous transformation extracted from Freesurfer's cortical segmentation process.

Table 1.

MNI Coordinates of the Nodes Forming the Resting-State Networks

		MNI coordinate
Network	Node	x	y	z	BA
Visual	Left visual cortex	−41	−77	3	19
	Right visual cortex	41	−72	1	19
Sensorimotor	Left primary somatosensory cortex	−38	−27	52	1
	Right primary somatosensory cortex	40	−25	50	1
Auditory	Left primary auditory cortex	−55	−21	7	41
	Right primary auditory cortex	57	−20	7	41
DMN	Precuneus	1	−57	28	31
	Left inferior parietal	−45	−65	30	39
	Right inferior parietal	52	−60	26	39
	Anterior cingulate	2	42	7	32
Left FP	Left angular gyrus	−40	−58	56	39
	Left middle frontal gyrus	−42	28	23	9
Right FP	Right angular gyrus	53	−50	43	39
	Right middle frontal gyrus	45	28	26	9
Frontoinsular	Left insula	−38	16	2	13
	Right insula	38	16	2	13
	Median cingulate	−2	12	40	32

MNI coordinates are obtained from the functional magnetic resonance imaging RSNs released by Laird et al. (2011) and Smith et al. (2009). Labels are assigned with FSL atlas tools (Desikan et al., 2006; Eickhoff et al., 2005; Jenkinson et al., 2012).

BA, Broadmann area; DMN, default mode network; FP, frontoparietal.

First, within-network reliability was quantified as the ICC of FC between all source pairs connecting given RSNs nodes, for each FC metric, RSN, and frequency band separately. Second, the reliability of seed-based FC maps was also quantified, using all nodes in Table 1 as seeds and computing the ICC of the FC between the seed and the remaining source locations. To enable within-subject comparisons, individual source meshes were transformed into a template mesh in MNI space with \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${N_{sources}} = 4739$$ \end{document} sources generated with Freesurfer from the subjects' cortical segmentations. FC values were then spatially smoothed with a 15 mm moving average filter.

Results

Reliability of whole-brain FC maps

ICCFC distribution

The distribution of ICCFC values, which quantify the test-retest reliability of eyes closed resting-state FC, is represented in Figure 1A. ICCFC was greatly dependent both on the FC metric and the frequency band. Highest values were reached for the PLV in alpha, beta, and gamma bands (averages 0.74–0.82). For d-ecor and lc-ecor, ICC values were generally smaller, and peaked in the beta band (average 0.23–0.37). For PLI, average ICC was lower (<0.1) for all frequency bands. The same effects across frequency bands and FC metric were obtained for the eyes open condition: results are displayed in Supplementary Figure S1 (Supplementary Data are available online at www.liebertpub.com/brain). Given that eyes open and closed conditions were also found to behave similarly in subsequent analyses, the following sections will for simplicity only describe the eyes closed results.

FIG. 1.

Distribution of ICC values for whole-brain FC during resting-state eyes closed. (A) Distribution of ICCFC values across links, for each FC metric and frequency band separately. The colored boxes indicate the position of the first, second, and third quartile, and the central lines span the range from the 5th to the 95th percentile. (B) Estimated a coefficients in the LMM described in the Reliability of whole-brain FC section. (C) Estimated b coefficients in the LMM. For (B) and (C), the colored boxes span the coefficients' 99% confidence interval. Stars indicate that the corresponding coefficient is significantly larger than zero (p < 0.01, Bonferroni corrected). FC, functional connectivity; ICC, intraclass correlation coefficient; LMM, linear mixed model. Color images available online at www.liebertpub.com/brain

BWC effect

The contribution of BWC to the previously described ICCFC was explored by introducing both ICCFC and ICCBWC in an LMM, as detailed in the Reliability of whole-brain FC section. The 99% confidence intervals of the LMM coefficients are displayed in Figure 1B and C. On the one hand, the a coefficients, which weigh ICCBWC-unrelated fixed effects on ICCFC, were significantly greater than zero (p < 0.01, Bonferroni corrected) for PLV for all frequency bands, lc-ecor in alpha and beta, d-ecor in beta, and PLI in alpha. On the other hand, the b coefficients, which weigh the contribution of ICCBWC to ICCFC, were significantly greater than zero (p < 0.01, Bonferroni corrected) for d-ecor and PLV for all frequency bands, while there were close to zero for leakage-corrected FC metrics (PLI and d-ecor).

Reliability was, however, unevenly distributed across links. Figure 2 represents the 5% of links with highest \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$ICCFC - b \cdot ICCBWC$$ \end{document} (i.e., the links with highest ICCFC after subtracting the contribution of BWC). The topography of these links differed across FC metrics and frequency bands. As a general trend, the following areas were found to be of especial relevance: anterior cingulate and temporal ROIs in delta, temporal, occipital, and posterior ROIs in theta, temporal, occipital, and anterior ROIs in alpha, central and temporal ROIs in beta, and parietal and frontal ROIs in gamma.

FIG. 2.

Spatial distribution of links with highest (ICCFC − b·ICCBWC). For each FC metric and frequency bands, the 5% links with highest (ICCFC − b·ICCBWC) are represented on a template brain, viewed from the top (frontal regions up, posterior regions down). For better visualization, the template brain surface was inflated with Caret software (Van Essen et al., 2001). BWC, beamformer weight correlation. Color images available online at www.liebertpub.com/brain

Dependence between ICC and FC strength

The dependence between ICC and FC was investigated with Spearman correlations. The results are listed in Table 2. ICC and FC were positively correlated, especially for d-ecor and PLV (r = 0.50–0.78). Correlations for lc-ecor and PLI were significant in alpha and beta bands (r = 0.14–0.46), and smaller (ABS|r| < 0.08) in delta, theta and gamma. When computing partial correlation between ICC and FC, while controlling for ICCBWC, correlation coefficients for d-ecor and PLV were much smaller, although remained significant for all frequency bands (r = 0.29–0.57). On the other hand, for lc-ecor and PLI, as expected, correlation coefficients between ICCFC and FC remained unchanged after controlling for ICCBWC.

Table 2.

Correlation Between ICC and FC

	Correlation ICCFC-FC				Partial correlation ICCFC-FC (ICCBWC)
	d-ecor	lc-ecor	PLV	PLI	d-ecor	lc-ecor	PLV	PLI
Delta	0.54^*	−0.04	0.78^*	−0.01	0.30^*	−0.04	0.57^*	−0.02
Theta	0.63^*	0.00	0.71^*	0.00	0.36^*	0.00	0.48^*	0.00
Alpha	0.70^*	0.16^*	0.62^*	0.27^*	0.40^*	0.16^*	0.43^*	0.27^*
Beta	0.70^*	0.46^*	0.50^*	0.14^*	0.56^*	0.46^*	0.30^*	0.13^*
Gamma	0.69^*	0.08^*	0.51^*	−0.02	0.43^*	0.08^*	0.29^*	−0.02

In the left side, the Spearman correlation coefficients between ICC and FC across links are listed, for all FC metrics and frequency bands separately, for the resting-state eyes closed condition. In the right side, Spearman partial correlation coefficients between ICCFC and FC, controlling for BWC ICC are listed. Asterisks denote significant correlations (Bonferoni corrected, p < 0.01).

BWC, beamformer weight correlation; FC, functional connectivity; ICC, intraclass correlation coefficient; PLI, phase lag index; PLV, phase locking value.

Within- and between-subject agreement in whole-brain FC maps

The within- and between-subject agreement of FC maps was computed with Kendall's W from pairwise maps of a single and distinct subjects. The distribution of W values for every frequency band and FC metric is displayed in Figure 3. First, highest within- and between-subject average W was found for PLV (0.59–0.98). Moreover, within-subject W was significantly higher than between-subject W for PLV (20% increment on average), and this difference was higher for high-frequency bands (22% increment in beta and gamma vs. 16% in delta). Second, high average W values (0.56–0.82) were also found for d-ecor, and within-subject W was higher than between-subject W for all frequency bands (8.4% on average). Third, for lc-ecor, within- and between-subject W were generally smaller (0.42–0.63) and similar, although for the beta band, a higher W (0.42–0.79) was obtained, and within-subject W significantly exceeded between-subject W. Fourth, for PLI both between- and within-subject W values were smaller (0.48–0.60), and within-subject W was significantly greater than between-subject W for alpha and beta bands.

FIG. 3.

Distribution of within- and between-subject agreement of FC maps. Each pair of same-color boxplots represents the within- (left) and between-subject (right) agreement for given a frequency band and FC metric, for the resting-state eyes closed condition. Stars indicate significantly higher within- than between-subject agreement (Wilcoxon rank sum test, p < 0.01, Bonferroni corrected). Color images available online at www.liebertpub.com/brain

Reliability of RSNs

Within-network reliability

By exploring the test-retest reliability of whole-brain FC in the Reliability of Whole-Brain FC Maps section, it was found that ICC values have a great spatial variability and some of the links with higher test-retest reliability bear some similarity with the RSNs described in the literature. The test-retest reliability of seven RSNs was therefore explored. First, the reliability of the FC of links connecting nodes of the same RSN was assessed with the ICC. Results are listed in Table 3. ICC values were dependent both on the FC metric and frequency band. For all RSNs, highest ICC values were found for PLV, usually ranging from 0.7 to 0.9 in alpha, beta, and gamma bands and smaller values from 0.4 to 0.7 in delta and theta bands. High ICC values (>0.7) were, however, also found in delta and theta bands for the sensorimotor and visual networks. d-ecor and lc-ecor presented greater variability across frequency bands than PLV. The ICC of both FC metrics topped in the beta band (ICC = 0.46–0.78 for d-ecor and ICC = 0.43–0.66 for lc-ecor), while low to medium ICC (<0.5) was found for delta, theta, and gamma bands. Although ICC was rather low in the alpha band, high values (>0.6) were found with d-ecor in the visual and sensorimotor networks. With PLV and envelope correlation metrics, higher ICC values were found for the sensory networks (especially sensorimotor and visual), than for DMN, FP, and frontoinsular networks. Finally, ICC values were medium to low for the PLI in all frequency bands and networks and did not exceed 0.44.

Table 3.

ICC of the Within-Network FC for the Resting-State Eyes Closed Condition

Network		Visual	Sensorimotor	Auditory	DMN	Left FP	Right FP	Frontoinsular
d-ecor	Delta	0.53	0.29	0.23	0.18	0.05	0.20	0.05
	Theta	0.69	0.49	−0.03	0.22	0.13	0.32	0.24
	Alpha	0.61	0.73	0.64	0.32	0.38	0.66	0.29
	Beta	0.69	0.76	0.58	0.42	0.46	0.78	0.59
	Gamma	0.34	0.46	0.43	0.32	0.20	0.10	0.15
lc-eco r	Delta	0.20	−0.03	0.10	0.10	−0.08	−0.22	−0.07
	Theta	0.07	−0.15	0.02	0.01	0.19	−0.06	−0.02
	Alpha	0.51	0.25	0.17	0.32	0.02	0.25	0.11
	Beta	0.66	0.54	0.52	0.45	0.43	0.55	0.47
	Gamma	−0.02	0.15	0.20	0.00	0.00	−0.32	−0.04
PLV	Delta	0.76	0.82	0.72	0.41	0.56	0.53	0.49
	Theta	0.85	0.91	0.67	0.59	0.67	0.71	0.72
	Alpha	0.86	0.95	0.71	0.75	0.87	0.82	0.67
	Beta	0.89	0.91	0.73	0.84	0.88	0.86	0.78
	Gamma	0.90	0.87	0.68	0.88	0.78	0.92	0.84
PLI	Delta	0.07	−0.09	−0.03	−0.06	0.19	0.17	−0.08
	Theta	−0.10	0.05	−0.19	0.06	−0.14	−0.07	−0.05
	Alpha	0.40	−0.10	0.03	0.08	0.11	0.12	0.07
	Beta	−0.14	−0.10	0.22	0.07	−0.02	−0.10	0.02
	Gamma	−0.05	−0.15	−0.01	−0.10	−0.21	0.11	−0.07

ICC values were computed for every frequency band, FC metric, and RSN separately. For RSNs with more than two nodes (DMN and frontoinsular), the average ICC across links within the network is listed.

RSNs, resting-state networks.

Reliability of seed-based FC maps

FC was computed from each seed to the remaining source location for each condition, frequency band, and FC metric separately. Then, the ICC of the FC values was computed for each link. The distributions of ICC values over all seed locations are displayed in Supplementary Figure S2. These ICC values were lower than the corresponding within-network ones for most seeds for PLV and d-ecor across frequency bands, lc-ecor in alpha, and beta and PLI in alpha indicating greater within- than between-network reliability. Moreover, the variability over frequency bands and FC metrics was similar than in the within-network analysis. First, highest values (ICC >0.6) were found for PLV, especially in alpha, beta, and gamma bands for all seed locations, although slightly higher ICC was found for sensory seeds and precuneus. Second, ICC values for d-ecor and lc-ecor were small for delta, theta, and gamma bands (third quartile <0.3), but higher values were obtained for alpha and especially beta bands (ICC 0.3–0.6). ICC values also varied somewhat across seeds, and were higher in the beta band in sensorimotor and left middle frontal seeds with lc-ecor and in the alpha band in the sensory and precuneus seeds for both lc-ecor and d-ecor. Third, ICC was poor for PLI (third quartile <0.2).

ICC values were unevenly distributed across brain regions. As an example, beta band ICC and average FC maps for three selected seeds are displayed for lc-ecor in Figure 5 and PLV in Figure 4 (corresponding figures for d-ecor and PLI can be found in the Supplementary Figures S3 and S4). Average FC was generally higher between regions belonging to the same network. For instance, seeds in sensory networks (visual, somatosensory, and auditory cortices) were strongly connected to their homologous contralateral areas. This was less evident for DMN, FP, and frontoinsular networks. ICC maps resembled average FC maps, showing high ICC for links between homologous sensory areas. DMN and FP networks presented lower ICC than sensory networks for all FC metrics. Moreover, FC increased along with ICC values for d-ecor, lc-ecor, and PLV, as revealed by the joint distribution of both magnitudes. This relation was, however, not linear, especially for FC metrics without zero-lag correction (PLV and d-ecor), for which the ICC increased first rapidly along with FC values and then saturated for medium to high FC (Fig. 5 and Supplementary Fig. S3).

FIG. 4.

FC and ICC for beta band resting-state eyes closed and lc-ecor. For three selected seeds (left primary somatosensory, left primary auditory, and precuneus), average FC maps (first row), ICC (second row), and the normalized joint histogram of ICC and FC values (third row) are displayed. FC maps are averaged over all sessions and subjects. The seed location is indicated with a white circle. Color images available online at www.liebertpub.com/brain

FIG. 5.

FC and ICC for beta band resting-state eyes closed and PLV. For three selected seeds (left primary somatosensory, left primary auditory, and precuneus), average FC maps (first row), ICC (second row), and the normalized joint histogram of ICC and FC values (third row) are displayed. FC maps are averaged over all sessions and subjects. The seed location is indicated with a white circle. PLV, phase locking value. Color images available online at www.liebertpub.com/brain

Discussion

In this work, we assessed the test-retest reliability of resting-state FC by exploring within- and between-subject variability of FC values in a cohort of 16 subjects. Following previous resting-state MEG literature (Brookes et al., 2011b; Hillebrand et al., 2005; Schoffelen and Gross, 2009), source time series were estimated with beamforming. Then, source space FC was computed between source locations using four distinct FC metrics: d-ecor, lc-ecor, PLV, and PLI. Within- and between-subject variability was employed to assess both reliability (with the ICC) and agreement (with Kendall's W) of the FC maps. First, overall trends over FC metrics and frequency bands were inspected for whole-brain connections, exploring all links connecting 66 cortical ROIs. Average ICCFC was highest in the PLV in alpha, beta, and gamma bands (0.74–0.82), medium for d-ecor and lc-ecor in the beta band (average 0.23–0.37), and considerably low (<0.1) for PLI across frequency bands. An LMM, including ICCBWC revealed that volume conduction, indirectly assessed through BWC, could be contributing to high test-retest reliability in PLV and d-ecor. When accounting for this effect, the overall test-retest reliability of whole-brain FC was found to be significant for PLV for all frequency bands, lc-ecor in alpha and beta, d-ecor in beta, and PLI in alpha. When evaluating within- and between-subject agreement separately with Kendall's coefficient of concordance W, medium to high values were obtained (0.5–0.97), depending on the frequency band and the FC metric. W was highest with PLV and lowest with lc-ecor. In addition, within-subject W was significantly higher than between-subject FC for d-ecor and PLV throughout frequency bands, PLI in alpha and beta, and lc-ecor in beta.

In addition, FC metric and frequency band were not the only important factors influencing test-retest reliability. Instead, the FC test-retest reliability was highly dependent on the specific link considered. First, ICCFC correlated positively with the FC strength, especially in alpha and beta bands. Second, when focusing on well-known RSNs, within-network connections were particularly reliable. Within each network, similar FC metrics and frequency band trends were found than for whole-brain FC: highest values found for PLV in alpha, beta, and gamma bands and for d-ecor and lc-ecor in beta band.

To our knowledge, this constitutes the first test-retest reliability study of MEG source space FC. Test-retest reliability of FC has, however, already been explored with both sensor space MEG/EEG and fMRI. Deuker et al. (2009) and Jin et al. (2011) employed the ICC to estimate the reliability of resting-state graph metrics derived from MEG sensor space mutual information. In the first study, Deuker et al. (2009) obtained high ICC (0.5–0.8) for strength, clustering, minimum path length, synchronizability, global efficiency, cost efficiency, and small world index in alpha band, and low to medium ICC (0–0.5) in other frequency bands. Interestingly, the ICC of these graph metrics during a working memory task was found to be higher than during resting state. In the latter study, Jin et al. (2011) found medium ICC (0.4–0.65) for nodal efficiency, nodal degree, and betweeness centrality in alpha and beta bands and lower ICC for theta and gamma bands. Using EEG and sensor space PLI, Hardmeier et al. (2014) obtained good reliability (ICC range: 0.5–0.8) for clustering coefficient, average path length, and small-world index in theta, alpha, and beta bands. Test-retest studies of sensor space graph metrics have therefore obtained divergent results, although they cannot be easily compared because of differences in the processing methodology, the choice of FC metric, and the test-retest interval. Using resting-state fMRI, Shehzad et al. (2009) could directly inspect the test-retest reliability of RSNs by computing FC between ROIs selected from the fMRI literature. They overall obtained low reliability (average ICC = 0.10–0.40), although higher ICC was found when considering significant correlations exclusively. This is in agreement with our results and suggests that stronger connections are more reliable. Although not focusing on reliability, Wens et al. (2014) studied short-term within-subject variability in MEG FC with four subjects who underwent 20 scan sessions in a single day. After selecting three seeds in visual, auditory, and sensorimotor cortices, they evaluated the spatial similarity between FC maps. They finally recommended between 2 and 11 sessions to produce a robust subject average, depending on the network. Moreover, in line with our work, they obtained that the sensorimotor network presented the highest within-subject similarity and hence required the least sessions to produce good subject-representative average FC maps.

An interesting result in our work is the frequency specificity of the reliability of RSNs. With PLV, RSNs reached medium to good reliability for all frequency bands, but highest ICC values were consistently found in alpha, beta, and gamma bands. In the previously commented graph metric studies (Deuker et al., 2009; Hardmeier et al., 2014; Jin et al., 2011), reliability was generally highest for alpha, reasonably good for theta and beta, and quite low for gamma. Although we cannot establish any direct analogy between our approach and sensor space graph metrics, given that high frequency oscillations are supposed to be more local than low-frequency oscillations (Buzsáki and Draguhn, 2004), we could speculate that we found higher gamma band reliability because gamma network properties are not well captured with sensor space FC. When using d-ecor and lc-ecor, reliability varied greatly among frequency bands and was higher in beta, followed by alpha. This finding is not particularly surprising, since MEG RSNs have been predominantly studied with envelope correlation in the alpha and beta bands. In fact, although the seminal work of (Brookes et al., 2011a,b; de Pasquale et al., 2010; Hipp et al., 2012) investigated RSNs with envelope correlation at several frequency bands, they concluded that most networks were best extracted in the beta band and alpha band in de Pasquale et al. (2010) and in Brookes et al. (2011b) for the DMN. In this study, we found that the test-retest reliability was quite high across networks for alpha and beta band with envelope correlation measures, although higher values were generally found in sensory networks and DMN.

In addition, we found that the choice of FC metric greatly influenced reliability. ICC was highest for PLV, medium for d-ecor and lc-ecor, and lowest for PLI, indicating that PLV estimates are the most reliable ones. The most striking difference was found between PLV and PLI: PLV presented consistently higher test-retest reliability than PLI. This could be due to the fact that both metrics capture different aspects of phase synchronization: while PLV evaluates how much the distribution of phase differences between brain oscillations deviates from a uniform distribution, PLI quantifies asymmetries in the phase differences over \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$] {0 , \pi } ] $$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$] { - \pi , 0} ]$$ \end{document} . The differences in test-retest reliability could be therefore caused by very different effects. On the one hand, it is possible that zero-lag synchronization represents meaningful and reliable FC, since various animal studies with intracranial EEG have proven the existence of true zero-phase lag synchronization between distinct areas (Fell and Axmacher, 2011; Frien et al., 1994; Roelfsema et al., 1997). On the other hand, it is also possible that volume conduction, which presents higher within- than between-subject variability, spuriously increased the magnitude of the test-retest reliability of PLV (and d-ecor) estimates. Recent literature of MEG RSNs has, however, predominantly employed d-ecor and lc-ecor. Envelope correlation has in fact produced FC maps, which resemble the ones obtained with fMRI (Brookes et al., 2011b; de Pasquale et al., 2010, 2012). Furthermore, lc-ecor refines d-ecor by eliminating zero-lag effects, using either time domain (Brookes et al., 2012) or frequency domain (Hipp et al., 2012; Wens et al., 2014) orthogonalization. The average FC maps obtained in this study were similar to the ones obtained in the previously commented studies: with lc-ecor, homologous areas in sensory networks were strongly connected and, contrary to PLV or d-ecor, the neighborhood of the seed did not present artifactual high values. Our results suggest nonetheless, that, although lc-ecor produces meaningful and unbiased FC maps, PLV estimates present higher test-retest reliability. Furthermore, phase and amplitude synchronization evaluate different mechanisms of information transfer, so that the choice of FC metric should depend on the underlying FC model. On the one hand, the phase of the oscillations is supposed to represent the degree of excitability of neuronal populations, and phase synchronization has been thought to contribute both to precise spike timing, especially in the gamma band, and to neural communication and plasticity (Fell and Axmacher, 2011). On the other hand, the amplitude of the oscillations represents the intensity of coherent neural activity, and amplitude correlation reflects the simultaneous activation of distinct brain areas.

Finally, it is important to note that test-retest reliability does not evaluate the goodness of the measurement (in this case, FC estimation). The test-retest reliability is instead a relative quantity, which depends on within- and between-subject variability. Any within-subject variability produces a reliability decrease, independent of the origin: instrumentation (e.g., MEG sensors deviating from their ideal working point), measurement situation (e.g., varying magnetic fields interfering in the MEG measurement, acoustic noise that could induce interfering brain activity in the subject) or the subject itself (true change in the coupling strength between two regions within a single individual). We can therefore not conclude that an FC estimate with low test-retest reliability (small ICC) contains a great amount of error or is inaccurate: it could as well be possible that it captures a true within-subject change. In this line, higher within-subject agreement in FC maps does not necessarily imply that FC is more stable within individuals: it could as well be due to volume conduction bias in the FC metric. Reliability relates instead on how well subjects can be distinguished in a given group. It is therefore indicative of whether FC estimates can be employed to characterize individual subjects within the group sample. We should also indicate that reliability is a population-specific quantity: for a fixed within-subject variability, the more heterogeneous the subject sample, the higher the reliability.

Conclusions

In conclusion, we provided the first source space test-retest reliability report of MEG resting-state FC. We computed both whole-brain FC between 66 ROIs and seed-based FC (focusing on RSNs nodes) using four commonly employed FC metrics: d-ecor, lc-ecor, PLV, and PLI. We then evaluated the test-retest reliability of the FC estimates with the ICC, and the within- and between-subject agreement of the FC maps with Kendall's W. ICC values were greatly dependent on the FC metric, the frequency band, and the source locations. Overall, highest ICC values were obtained for PLV, even reaching 0.8 or higher for broad regions, especially in alpha, beta, and gamma bands. d-ecor and lc-ecor presented similar reliability, and generally yielded low to medium ICC (<0.5), except for alpha and beta bands, where ICC values of 0.6 were reached. PLI estimates yielded the lowest ICC throughout frequency bands, and ICC values rarely exceeded 0.4. Joint inspection of the ICCFC (test-retest reliability of FC) and the ICCBWC (test-retest reliability of beamformer weights correlation) revealed that volume conduction (assessed through the BWC) could be partly responsible for the high ICCFC obtained in PLV and d-ecor estimates. Across FC metrics and frequency bands (but especially in alpha and beta bands), ICCFC values correlated positively with FC strength, and within-network connections presented higher ICCFC than between-network connections. Both findings point out a link between stronger connections and higher test-retest reliability. Reliability was defined in this study as a relative quantity (Weir, 2005), which represents the fraction of the measurement variability that is accounted for by between-subject variability. Our results are therefore encouraging for future studies, which would aim at using FC strength in RSNs to characterize individual subjects, especially when using PLV in theta, alpha, beta, and gamma bands, or envelope correlation measures in the alpha and beta bands.

Footnotes

Acknowledgments

The authors are grateful to the research participants for their participation in this study. Research by P.G. has been supported by a PICATA predoctoral contract of the Moncloa Campus of International Excellence (UCM–UPM).

Author Disclosure Statement

No competing financial interests exist.

References

Acar

, Makeig

. 2010. Neuroelectromagnetic forward head modeling toolbox. J Neurosci Methods, 190:258–270.

Baldassarre

, Lewis

, Committeri

, Snyder

, Romani

, Corbetta

. 2012. Individual variability in functional connectivity predicts performance of a perceptual task. Proc Natl Acad Sci USA, 109:3516–3521.

Bartlett

, Frost

. 2008. Reliability, repeatability and reproducibility: analysis of measurement errors in continuous variables. Ultrasound Obstet Gynecol, 31:466–475.

Biswal

, Zerrin Yetkin

, Haughton

, Hyde

. 1995. Functional connectivity in the motor cortex of resting human brain using echo-planar MRI. Magn Reson Med, 34:537–541.

Bowyer

, Gjini

, Zhu

, Kim

, Moran

, Rizvi

, et al. 2015. Potential biomarkers of schizophrenia from MEG resting-state functional connectivity networks: Preliminary data. J Behav Brain Sci, 05:1–11.

Brookes

, Hale

, Zumer

, Stevenson

, Francis

, Barnes

, et al. 2011a. Measuring functional connectivity using MEG: Methodology and comparison with fcMRI. Neuroimage, 56:1082–1104.

Brookes

, Woolrich

, Luckhoo

, Price

, Hale

, Stephenson

, et al. 2011b. Investigating the electrophysiological basis of resting state networks using magnetoencephalography. Proc Natl Acad Sci USA, 108:16783–16788.

Brookes

, Woolrich

, Barnes

. 2012. Measuring functional connectivity in MEG: A multivariate approach insensitive to linear source leakage. Neuroimage, 63:910–920.

Buzsáki

, Draguhn

. 2004. Neuronal oscillations in cortical networks. Science, 304:1926–1929.

10.

Deco

, Jirsa

, McIntosh

. 2011. Emerging concepts for the dynamical organization of resting-state activity in the brain. Nat Rev Neurosci, 12:43–56.

11.

de Pasquale

, Della Penna

, Snyder

, Lewis

, Mantini

, Marzetti

, et al. 2010. Temporal dynamics of spontaneous MEG activity in brain networks. Proc Natl Acad Sci USA, 107:6040–6045.

12.

de Pasquale

, Della Penna

, Snyder

, Marzetti

, Pizzella

, Romani

, Corbetta

. 2012. A cortical core for dynamic integration of functional networks in the resting human brain. Neuron, 74:753–764.

13.

Desikan

, Ségonne

, Fischl

, Quinn

, Dickerson

, Blacker

, et al. 2006. An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest. Neuroimage, 31:968–980.

14.

Deuker

, Bullmore

, Smith

, Christensen

, Nathan

, Rockstroh

, Bassett

. 2009. Reproducibility of graph metrics of human brain functional networks. Neuroimage, 47:1460–1468.

15.

Engels

, Stam

, van der Flier

, Scheltens

, de Waal

, van Straaten

. 2015. Declining functional connectivity and changing hub locations in Alzheimer's disease: An EEG study. BMC Neurol, 15:145.

16.

Fell

, Axmacher

. 2011. The role of phase synchronization in memory processes. Nat Rev Neurosci, 12:105–118.

17.

Fischl

, Salat

, Busa

, Albert

, Dieterich

, Haselgrove

, et al. 2002. Whole brain segmentation. Neuron, 33:341–355.

18.

Frien

, Eckhorn

, Bauer

, Woelbern

, Kehr

. 1994. Stimulus-specific fast oscillations at zero phase between visual areas V1 and V2 of awake monkey. Neuroreport, 5:2273–2277.

19.

Gramfort

, Luessi

, Larson

, Engemann

, Strohmeier

, Brodbeck

, et al. 2014. MNE software for processing MEG and EEG data. Neuroimage, 86:446–460.

20.

Hall

, Woolrich

, Thomaz

, Morris

, Brookes

. 2013. Using variance information in magnetoencephalography measures of functional connectivity. Neuroimage, 67:203–212.

21.

Hardmeier

, Hatz

, Bousleiman

, Schindler

, Stam

, Fuhr

. 2014. Reproducibility of functional connectivity and graph measures based on the phase lag index (PLI) and weighted phase lag index (wPLI) derived from high resolution EEG. PLoS One, 9:e108648.

22.

Hillebrand

, Barnes

, Bosboom

, Berendse

, Stam

. 2012. Frequency-dependent functional connectivity within resting-state networks: An atlas-based MEG beamformer solution. Neuroimage, 59:3909–3921.

23.

Hillebrand

, Singh

, Holliday

, Furlong

, Barnes

. 2005. A new approach to neuroimaging with magnetoencephalography. Hum Brain Mapp, 25:199–211.

24.

Hipp

, Hawellek

, Corbetta

, Siegel

, Engel

. 2012. Large-scale cortical correlation structure of spontaneous oscillatory activity. Nat Neurosci, 15:884–890.

25.

Jin

S-H

, Seol

, Kim

, Chung

. 2011. How reliable are the functional connectivity networks of MEG in resting states?. J Neurophysiol, 106:2888–2895.

26.

Lachaux

, Rodriguez

, Martinerie

, Varela

. 1999. Measuring phase synchrony in brain signals. Hum Brain Mapp, 8:194–208.

27.

Laird

, Fox

, Eickhoff

, Turner

, Ray

, McKay

, et al. 2011. Behavioral interpretations of intrinsic connectivity networks. J Cogn Neurosci, 23:4022–4037.

28.

Legendre

. 2005. Species associations: The Kendall coefficient of concordance revisited. J Agric Biol Environ Stat, 10:226–245.

29.

López

, Bruña

, Aurtenetxe

, Pineda-Pardo

JÁ

, Marcos

, Arrazola

, et al. 2014. Alpha-band hypersynchronization in progressive mild cognitive impairment: A magnetoencephalography study. J Neurosci, 34:14551–14559.

30.

Luckhoo

, Hale

, Stokes

, Nobre

, Morris

, Brookes

, Woolrich

. 2012. Inferring task-related networks using independent component analysis in magnetoencephalography. Neuroimage, 62:530–541.

31.

Maldjian

, Davenport

, Whitlow

. 2014. Graph theoretical analysis of resting-state MEG data: Identifying interhemispheric connectivity and the default mode. Neuroimage, 96:88–94.

32.

McGraw

, Wong

. 1996. Forming inferences about some intraclass correlation coefficients. Psychol Methods.

33.

Niso

, Bruña

, Pereda

, Gutiérrez

, Bajo

, Maestú

, del-Pozo

. 2013. HERMES: Towards an integrated toolbox to characterize functional and effective brain connectivity. Neuroinformatics, 11:405–434.

34.

Oostenveld

, Fries

, Maris

, Schoffelen

J-M

. 2011. FieldTrip: Open source software for advanced analysis of MEG, EEG, and invasive electrophysiological data. Comput Intell Neurosci, 2011:156869.

35.

Pereda

, Quiroga

, Bhattacharya

. 2005. Nonlinear multivariate analysis of neurophysiological signals. Prog Neurobiol, 77:1–37.

36.

Roelfsema

, Engel

, König

, Singer

. 1997. Visuomotor integration is associated with zero time-lag synchronization among cortical areas. Nature, 385:157–161.

37.

Rosazza

, Minati

. 2011. Resting-state brain networks: Literature review and clinical applications. Neurol Sci, 32:773–785.

38.

Rosenblum

, Pikovsky

, Kurths

, Schäfer

, Tass

. 2001. Neuro-Informatics and Neural Modelling, Handbook of Biological Physics. New York, NY: Elsevier.

39.

Schoffelen

, Gross

. 2009. Source connectivity analysis with MEG and EEG. Hum Brain Mapp, 30:1857–1865.

40.

Scholtes

, Terwee

, Poolman

. 2011. What makes a measurement instrument valid and reliable?. Injury, 42:236–240.

41.

Ségonne

, Pacheco

, Fischl

. 2007. Geometrically accurate topology-correction of cortical surfaces using nonseparating loops. IEEE Trans Med Imaging, 26:518–529.

42.

Shehzad

, Kelly

AMC

, Reiss

, Gee

, Gotimer

, Uddin

, et al. 2009. The resting brain: Unconstrained yet reliable. Cereb Cortex, 19:2209–2229.

43.

Shrout

, Fleiss

. 1979. Intraclass correlations: Uses in assessing rater reliability. Psychol Bull, 86:420–428.

44.

Smith

, Fox

, Miller

, Glahn

, Fox

, Mackay

, et al. 2009. Correspondence of the brain's functional architecture during activation and rest. Proc Natl Acad Sci USA, 106:13040–13045.

45.

Stam

, Nolte

, Daffertshofer

. 2007. Phase lag index: Assessment of functional connectivity from multi channel EEG and MEG with diminished bias from common sources. Hum Brain Mapp, 28:1178–1193.

46.

Taulu

, Simola

. 2006. Spatiotemporal signal space separation method for rejecting nearby interference in MEG measurements. Phys Med Biol, 51:1759–1768.

47.

Tewarie

, Schoonheim

, Schouten

, Polman

, Balk

, Uitdehaag

BMJ

, et al. 2015. Functional brain networks: Linking thalamic atrophy to clinical disability in multiple sclerosis, a multimodal fMRI and MEG Study. Hum Brain Mapp, 36:603–618.

48.

van den Heuvel

, Hulshoff Pol

. 2010. Exploring the brain network: A review on resting-state fMRI functional connectivity. Eur Neuropsychopharmacol, 20:519–534.

49.

van den Heuvel

, Mandl

RCW

, Kahn

, Hulshoff Pol

. 2009. Functionally linked resting-state networks reflect the underlying structural connectivity architecture of the human brain. Hum Brain Mapp, 30:3127–3141.

50.

Van Veen

, van Drongelen

, Yuchtman

, Suzuki

. 1997. Localization of brain electrical activity via linearly constrained minimum variance spatial filtering. IEEE Trans Biomed Eng, 44:867–880.

51.

Weir

. 2005. Quantifying test-retest reliability using the intraclass correlation coefficient and the SEM. J Strength Cond Res, 19:231–240.

52.

Wens

, Bourguignon

, Goldman

, Marty

, Op de Beeck

, Clumeck

, et al. 2014. Inter- and intra-subject variability of neuromagnetic resting state networks. Brain Topogr, 27:620–634.

53.

Yamashita

, Kawato

, Imamizu

. 2015. Predicting learning plateau of working memory from whole-brain intrinsic network connectivity patterns. Sci Rep, 5:7622.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.42 MB

0.53 MB

0.38 MB

0.39 MB