Resting state MEG oscillations show long-range temporal correlations of phase synchrony that break down during finger movement

The capacity of the human brain to interpret and respond to multiple temporal scales in its surroundings suggests that its internal interactions must also be able to operate over a broad temporal range. In this paper, we utilize a recently introduced method for characterizing the rate of change of the phase difference between MEG signals and use it to study the temporal structure of the phase interactions between MEG recordings from the left and right motor cortices during rest and during a finger-tapping task. We use the Hilbert transform to estimate moment-to-moment fluctuations of the phase difference between signals. After confirming the presence of scale-invariance we estimate the Hurst exponent using detrended fluctuation analysis (DFA). An exponent of >0.5 is indicative of long-range temporal correlations (LRTCs) in the signal. We find that LRTCs are present in the α/μ and β frequency bands of resting state MEG data. We demonstrate that finger movement disrupts LRTCs correlations, producing a phase relationship with a structure similar to that of Gaussian white noise. The results are validated by applying the same analysis to data with Gaussian white noise phase difference, recordings from an empty scanner and phase-shuffled time series. We interpret the findings through comparison of the results with those we obtained from an earlier study during which we adopted this method to characterize phase relationships within a Kuramoto model of oscillators in its sub-critical, critical, and super-critical synchronization states. We find that the resting state MEG from left and right motor cortices shows moment-to-moment fluctuations of phase difference with a similar temporal structure to that of a system of Kuramoto oscillators just prior to its critical level of coupling, and that finger tapping moves the system away from this pre-critical state toward a more random state.


Introduction
The human MEG is a rich time-varying signal that provides an important window into ongoing brain processing (Cohen, 1968;Riehle and Vaadia, 2004;Buzsaki, 2006;Hansen et al., 2010). In the frequency domain the MEG is characterized by spectral peaks overlaying 1 f noise, such that its power at each frequency diminishes in an inverse relationship (Buzsaki, 2006). The power of the spectral peaks varies with sensor location and task performed, but in general terms there are peaks of spectral power in the θ , α/µ, β, and γ bands (Berger, 1929). Underlying these peaks is oscillatory neuronal activity due to network synchronization (Buzsaki, 2006). In the motor system there is particular interest in the behavior of spectral peaks in the α, µ, and β bands whose oscillatory activity dominates in resting state MEG recorded over the sensory-motor cortex and which is diminished by transient movement of the contralateral muscle groups (Pfurtscheller, 1977;Pfurtscheller and da Silva, 1999). During bimanual movement and motor learning there are changes in EEG synchronization between different motor areas including left and right motor cortices (Andres et al., 1999). During sustained muscle contraction there is βband synchronization between motor cortex EEG/MEG and the contralateral EMG (Farmer, 1998;Halliday et al., 1998;James et al., 2008). Importantly neural synchronization is weak and varies with time even during a sustained task. These fluctuations of synchrony can be detected using a variety of analytical techniques: time-varying coherence, wavelet coherence, and optimal spectral tracking (Brittain et al., 2009;Fries, 2009). The observation that oscillatory brain signals may come in and out of synchronization with each other spontaneously and with task, and idea that synchronization lends salience to a signal, has led to the hypothesis of "communication through coherence" (CTC). Implicit is the idea that the variability of the effective coupling is in itself of significance (Fries, 2009).
Because neural signals exhibit moment-to-moment fluctuations in their amplitude and phase it is possible to characterize the temporal correlations within such time series. The order within amplitude fluctuations of bandpass-filtered MEG and EEG time series has been previously quantified through estimation of the Hurst exponent (Linkenkaer- Hansen et al., 2001Hansen et al., , 2004. The exponent of a power law relationship between scale and fluctuation magnitude within that scale, characterizes the presence of long-range temporal correlations in the time series. An exponent of 0.5 indicates white noise with no temporal correlation, whereas exponents of >0.5 or <0.5 indicate either long-range or an anti-correlated time series, respectively. These data have been interpreted within the framework of the "critical brain hypothesis" in which the presence of spatio-temporal correlations in brain signals at single unit level, at LFP level and at surface signal level (such as MEG and EEG oscillations) are thought to reflect a system operating close to a critical regime (Beggs and Plenz, 2003;Shew et al., 2009;Chialvo, 2010). Given the importance of neural synchrony and its fluctuation for communication between brain regions, it is of interest to ask whether in human brain signals there is evidence that interareal synchronization exists within or close to a critical regime (Kitzbichler et al., 2009).
In this paper we explore whether there are LRTCs in the moment-to-moment fluctuations of the MEG phase difference between different brain regions and whether, if present, these are disrupted by finger movement. Using resting state human MEG data, we determine first if there is power law scaling and if so we derive the DFA exponent of the fluctuations in the rate of change of phase difference between two simultaneously recorded MEG signals obtained from left and right primary motor cortices. We next explore the effects of bilateral muscle activation (finger tapping) on the MEG phase difference scaling and the associated DFA exponent. Finger movement is known to disrupt sensori-motor cortex oscillations (Pfurtscheller, 1977) and may also act as feedback perturbation of the type that has been shown to reduce LRTCs in the amplitude envelope of the EEG and MEG (Linkenkaer-Hansen et al., 2004). Here we use finger tapping for its potential to disrupt the temporal organization of left and right MEG phase differences in the resting state. We show that in the resting state the momentto-moment fluctuations of phase difference between the left and right motor cortex MEG time series has both power law scaling and an estimated DFA exponent >0.5 in the α/µ and β frequency bands, indicative of long-range temporal correlations within moment-to-moment fluctuations of phase difference. Finger tapping movements, when compared to the resting state, are associated with a significantly lower DFA exponent in the MEG moment-to-moment fluctuations of phase difference in the α/µ and β frequency bands. We validate our results by applying the same analysis to data with Gaussian white noise phase difference, recordings from an empty MEG scanner and phaseshuffled time series. The demonstration of the presence of scaling with exponent values consistent with the presence of longrange temporal correlations in this data provides an important connection between the critical brain and neural synchrony paradigms.

Materials and Methods
A key novel component of this work is our application of a method that allows detection of the presence of longrange temporal correlations (LRTCs) in the moment-to-moment fluctuations of phase difference between two neurophysiological time series. In Sections 2.1-2.3, we first detail the steps involved in the extraction of the rate of change of the phase difference time series from the neurophysiological data and then briefly summarize the two techniques used to rigorously assess the presence of LRTCs, along with statistical testing. In Section 2.4 we give details of the subjects and data recordings used to collect the data. Finally, Section 2.5 describes the control time series used to validate the results, including analysis of synthetic noise and empty scanner and phase-shuffled MEG data.
The overall methodology for the extraction and analysis of the neurophysiological signals used in this manuscript follows that proposed by the authors in Botcharova et al. (2014) with modifications suitable to its application to MEG data. The methodology is illustrated in a step-by-step fashion in Figure 1. Below, we detail each step in turn.

FIGURE 1 |
Step-by-step illustration of the method. Each step corresponds to a pair of plots. The time series shown and analyzed in the panels on the left are taken from resting state human MEG data, and those shown and analyzed in the panels on the right are MEG data recorded while the same subject was tapping both index fingers simultaneously. The red and blue colors in steps 2-4 correspond, respectively to the MEG time series for the left and the right motor cortices. (Step 1) Data from the full set of 275 MEG channels is considered (top panel shows the recording set up, panels below this show sensor level MEG in Tesla). (Step 2) The data is beamformed to obtain two time series corresponding to the left and right motor cortices, respectively, measured in nano-Ampere meters (nAm). (Step 3) Each beamformed time series is bandpass-filtered to a frequency of choice as described later in the text. (Step 4) The analytic signal s a (t) for each time series is found and plotted on a phase space diagram with the axes representing the real and imaginary parts of s a (t). (Step 5) The phase difference between the two time series is found by taking its arctangent and the first time derivative of this phase difference is then taken to produce the time series shown as plotted in blue. (Step 6) DFA (bottom panels) is applied to the time series shown in step 5 and ML-DFA is used to ascertain whether the detrended fluctuations follow a power law; if they do, the scaling exponent α is calculated by linear regression. For the data shown here, both DFA plots were judged to follow a power law. The top panel brain image was generated using SPM8 (Litvak et al., 2011).

Beamforming
The MEG data were initially processed using a scalar linearly constrained minimum variance beamformer (van Veen et al., 1997) in order to estimate timecourses of source strength at locations of interest in the left and right motor cortices. The forward models, covariance matrices and calculation of beamformer weights are described in Brookes et al. (2011). The weights for left and right motor cortex were based on those contralateral voxels showing the largest change in β-band power for right and left hand movement, respectively (Pfurtscheller and da Silva, 1999). This was done individually for each of the 7 subjects.

Bandpass filtering
The use of the Hilbert transform (see Section 2.1.3) requires a narrowband signal for the phase and amplitude components of the signal to be straightforward to separate (Boashash, 1992;Onslow et al., 2011). We filtered the data with a finite impulse response (FIR) filter. The filter order of the FIR filter was set to include three cycles of the lower-frequency component of the frequency band. A similar approach has been used in many recent studies, e.g., (Linkenkaer-Hansen et al., 2001).
The frequency bands were selected as 2 Hz intervals spanning the full range of frequencies to 4.5-39.5 Hz following a technique used in Klimesch et al. (1998). In this paper, the authors identify peaks in the power spectrum, which are then used to define frequency bands . . ., f − 2 to f, f to f + 2, f + 2, f + 4, . . . Here, the resting state power spectra have two peaks in two frequency bands one at 10.5 Hz in the α/µ frequency range and another at 21.5 Hz in the β frequency range (see Section 3.1). For this reason the time series were filtered in bands 4.5-6.5 Hz, 6.5-8.5 Hz, . . ., 14.5-16.5 Hz in the θ , α/µ frequency range, and in bands 15.5-17.5 Hz, . . ., 37.5-39.5 Hz in the β and γ frequency ranges.

Phase Difference
The phase relationship φ 1 (t) − φ 2 (t) between two different time series s 1 (t) and s 2 (t) was calculated using the respective Hilbert transform of the signals H[s 1 (t)] and H[s 2 (t)] (Pikovsky et al., 2003): The time series φ 1 (t) − φ 2 (t) is an unbounded process because φ 1 (t) and φ 2 (t) themselves are unbounded as long as the signals s 1 (t) and s 2 (t) continue to evolve as time increases. As we use detrended fluctuation analysis (DFA), see Section 2.2.1, to assess the presence of LRTCs and DFA in its standard form assumes a bounded signal, in this paper, we characterize phase synchronization in terms of the rate of change (first time derivative) of the phase difference time series φ 1 (t) − φ 2 (t).

Characterization of the LRTCs
Long-range temporal correlations (LRTCs) are present in a time series when its autocorrelation function, R ss (τ ), has a decay that takes the form of a power law, such that R ss (τ ) ∼ Cτ −α where τ is the time lag between observations and α is the exponent of the power law. One approach to determining this exponent is to estimate the Hurst exponent, which directly relates to α through α = 2 -2H (Beran, 1994;Taqqu et al., 1995), using DFA (Peng et al., 1994).

Detrended Fluctuation Analysis (DFA)
To calculate the DFA exponent, the time series (in this case the rate of change of phase difference between the two MEGs) is first detrended and then cumulatively summed. The root mean square error is then calculated when this signal is fitted by a line over different window sizes (or scales). Extensions of the technique can be used to fit any polynomial to each window, however, here we only consider linear detrending. If the time series is selfsimilar, there will be power law scaling between the residuals (or detrended fluctuations) and the window sizes. In log space, this power law scaling yields a linear relationship between residuals and window sizes, the so-called DFA fluctuation plot, and the DFA exponent H is obtained using least squares linear regression. A DFA exponent in the range 0.5 < H < 1 indicates the presence of LRTCs. An exponent of 0 < H < 0.5 is obtained when the time series is anti-correlated, H = 1 represents pink noise. Gaussian white noise has an exponent of H = 0.5.

Data Stitching
In our application of DFA to MEG data, we used 20 window sizes with a logarithmic scaling and a minimum window of 1 s, guaranteeing several cycles of the slowest oscillation. To improve the estimate of the DFA, trials during which a given subject was performing the same task were stitched together to produce a single, longer time series for each subject (see below for more detail). We used a maximum window size of 30 s, which is the length of a single trial, and equivalent to N/10 where N is the length of the stitched time series. Stitching is a commonly used step and does not affect the numerical results obtained by DFA if the window sizes are set appropriately (Chen et al., 2002). Since the windows are chosen on a logarithmic scale, there will be some windows at some window sizes in which distinct trials of data overlap. This is acceptable for longer data sets (as formed by stitching here) because at least half of the windows will always be contained within a continuously recorded segment of data. Empirical evidence for the validity of our stitching procedure is provided in Supplementary Material.
When stitching is applied to time series, it is done between steps 5, 6 of our method (Figure 1). Namely, the dividing and stitching together of trials should be performed after beamforming, filtering and calculating the phase difference. This minimizes the chance of edge artifacts from the Hilbert transform (Kramer et al., 2008;Widmann and Schröger, 2012).

Assessing the Validity of DFA
An interpretable DFA exponent is obtained when the log-log plot of fluctuation size against window size is truly linear, i.e., there is power law scaling (see Section 2.2.1). Since there is no a priori means of confirming that a signal is indeed selfsimilar, an exponent can always be obtained even though the DFA fluctuation plot may not necessarily be linear-the only certainty being that it will be increasing (albeit not necessarily monotonically so) with window size.
We used the technique described in Botcharova et al. (2013) to determine whether each DFA fluctuation plot is wellapproximated by a linear model. This is a heuristic technique, which has been tested extensively on known time series in order to verify the robustness of its conclusions (Botcharova et al., 2013). The technique first fits the DFA fluctuation plot with a number of candidate models. These models (polynomial, root, logarithmic, exponential, and spline with 2-4 sections) were chosen based on published characterizations of DFA fluctuation plots obtained with timeseries of known properties, e.g., with superposed trends or with noise. The number of parameters in these models range from 2 for the linear model to 8 for the foursegment spline model. The technique then maximizes a function that is similar in form to a log-likelihood: where lF scaled are normalized fluctuation magnitudes, i indexes the windows, and where f (lns) is the fitted model. The fits of all models are compared using the Akaike Information Criterion (AIC), which discounts for the number of parameters needed to fit the model. The DFA exponent is accepted as being valid only if the best fitting model (that with the lowest value of AIC) is the linear model. Only those time series that were not rejected for not being linear contributed to our results.

Statistical Tests of MEG Data Analysis
We corrected for multiple comparisons using a permutation test in each of the three bands (α/µ, β, and γ ) (Nichols and Holmes, 2002). In particular, we found the mean of the difference in DFA exponent between the resting state and finger tapping task for each subject and each sub-band of the α/µ, β, and γ bands. We then swapped the DFA exponents for resting state and finger tapping task in some of the subjects. This could be done in 2 7 = 128 ways, and we performed each one. For each permutation we found the mean of the difference in DFA exponent between the tasks defined by the shuffled labels. We then took the maximum of these mean differences across all sub-bands. This gave a null distribution of the differences in DFA exponent between tasks that we would expect by chance. We repeated this procedure to create separate significance thresholds for each of three bands of interest (α/µ, β, and γ ).

Participants
Seven healthy right handed subjects took part in the MEG experiments. The study was approved by the University of Nottingham Medical School ethics committee. Details of the experimental procedure are described in Brookes et al. (2011). In brief, MEG data were recorded using the third order gradiometer configuration of a 275 channel CTF MEG system at a sampling rate of 600 Hz. The scanner was housed inside a magnetically shielded room (MSR) and a 150 Hz low pass anti-aliasing hardware filter was applied. All subjects underwent a single experiment which comprised 10 trials of rest and movement. Subjects' eyes were open during both rest and movement conditions. A single trial comprised 30 s during which both left and right index fingers were moved together, at a rate of approximately 1 finger movement per second and 30 s of rest.
The movement itself comprised abductions and extensions of the index fingers. The motor task was cued visually via projection through a waveguide in the MSR onto a back projection screen located 40 cm in front of the subject. During data acquisition the location of the subject's head within the scanner was measured by energizing coils placed at 3 fiducial points on the head (nasion, left preauricular, and right preauricular). If any subject moved more than 5 mm during the experiment, data from that subject was discarded. Following data acquisition, the positions of the coils were measured relative to the subject's head shape using a 3D digitizer (Polhemus isotrack). An MPRAGE structural MR image was acquired using a Philips Achieva 3T MRI system (1 mm 3 isotropic resolution, 256 × 256 × 160 matrix size, TR = 8.3 ms, TR = 3.9 ms, TI = 960 ms, shot interval = 3 s, FA = 8 • and SENSE factor = 3). The locations of the fiducial markers and MEG sensors with respect to the brain anatomy were determined by matching the digitized head surface to the head surface extracted from the 3T anatomical MRI.

Control Data
In order to search for artifactual effects, we applied our method to three different classes of both synthetic and experimental control time series.

Synthetic Data with Noise Phase Difference
Two time series s 1 (t) and s 2 (t) were constructed such that their phase difference was white Gaussian noise with a known DFA exponent of 0.5 (Hosking, 1981). The signals were calculated by Here, ω = 1 and was included because white noise time series alter sharply at each new innovation and the Hilbert transform can produce artifacts when applied to such data. This situation does not arise in physiological data because the sampling frequency with which physiological data are recorded yields a smooth signal.
Seven pairs of time series were created to match the 7 human subjects. Each time series consisted of 10 segments of X(t) of length 30 s which were stitched together to construct the two signals s 1 (t) and s 2 (t) (see Section 2.2.2 for more detail). These signals were not beamformed because there was no possibility of unwanted noise or crosstalk. Window sizes used for application of DFA were identical to those used for physiological data, namely logarithmically spaced windows with minimum 1 s and maximum 30 s. Although there is no sampling frequency associated with noise simulations, the minimum window was set to 600 time samples for consistency.
The difference between the phases of the two signals s 1 and s 2 is precisely tX(t) yielding a first time derivative of X(t).
On applying DFA, we would therefore expect the prescribed exponent to be recovered if no bandpass filter was applied. This control data therefore makes it possible to quantify and rule out any artifact potentially introduced by filtering.

Scanner Noise
We analyzed 7 data sets recorded from an empty MEG scanner. Time series were obtained from sensors corresponding to the approximate locations of the right and left motor cortices, with a sampling frequency of 600 Hz. These continuous time series of length 300 s were divided into 10 segments of length 30 s each, and reordered randomly to replicate the stitching procedure applied to MEG data. This reordering was performed after step 5 of our method (see Figure 1 and Section 2.2.2 for more detail), which ensured consistency with analysis of the human MEG data. DFA exponents were then obtained for the 7 time series pairs, and checked for validity using ML-DFA. Window sizes used for DFA were identical to those used for the physiological data.

Statistical Tests
The DFA exponents of the synthetic pairs of time series constructed in Sections 2.5.1, 2.5.2 were tested for normality using a one-sample Kolmogorov-Smirnov test, for each frequency band individually as for MEG data. These tests showed that a null hypothesis of normally distributed DFA exponents could not be rejected for any synthetic time series pair or frequency band where more than one DFA exponent was judged to be valid by ML-DFA (p-values not shown). This was confirmed using the Jaque-Bera test for normality (data not shown). As normal distributions could not be rejected, DFA exponents obtained from the analysis of the above time series were tested for statistical difference with exponents obtained for resting state human MEG data using a Student's t-test (unpaired).

Phase Shuffling
A permutation 1 of the time derivative of the phase difference time series (the time series in step 5 of our methodology in Figure 1) was performed 100 times yielding 100 DFA exponents. These 100 exponents were tested for significant differences with the 7 human MEG data using the technique presented in Phipson and Smyth (2010). The p-value was approximated using 1 Since shuffling is applied to the time derivative of the phase difference time series, this operation does not equate with phase randomization techniques typically used in the literature to assess the presence of non-linear structure in a given time series.
where b is the number of DFA exponents obtained from the permuted time series that are greater in size than that of the correctly ordered human MEG data, and m is the number of permutations used m = 100. This is an approximation to the p-value when the number of possible permutations of the time series is very large, as in this case where the number of possible permutations is 180,000! (factorial).

Spectral Analysis
We analyzed data from the MEG oscillations of the left and right motor cortices during rest and during finger tapping tasks. The power spectra of time series for all 7 subjects were calculated (see Figure 2) and using pooled power spectral analysis (Halliday et al., 1995) the most prominent spectral peaks were determined (10.5 and 21.5 Hz). As explained in the Methods, these peaks were used to set the boundaries of the ranges of frequency bands used for bandpass-filtering. Figure 3 shows for 7 subjects the mean ± SD valid DFA exponents calculated from the rate of change of phase difference between left and right hemisphere MEGs. The time series were bandpass-filtered in regular 2 Hz bands prior to analysis (see Section 2.1.2). As recommended by Pfurtscheller and da Silva (1999), we state the boundary values of each frequency band. These data are shown for each 2 Hz frequency band across the frequency ranges (Figure 3). Only experiments associated with fluctuation plots judged to show linear scaling using ML-DFA are included in the data (see Table 1). The results reveal an estimated Hurst exponent for the rate of change of phase difference between Frontiers in Physiology | www.frontiersin.org FIGURE 3 | Mean ± SD DFA exponents for "rest" and "both" conditions against the distributions of DFA exponents for phase-shuffled MEG time series for the α/µ, β, and γ frequency ranges (A-C). Mean DFA exponents are in dark blue for resting state and cyan for the finger tapping task, with vertical bars showing standard deviation. The mean DFA exponents of the phase-shuffled resting state are in magenta, with those of the phase-shuffled finger tapping condition in black. Statistical significance for each band and each condition was assessed using permutation p-values calculated with the methodology described in Section 2.5.4 (* for p < 0.05; ** for p < 0.01). left and right motor cortex MEGs that is >0.5 (i.e., with LRTCs) at rest and during finger tapping in the α/µ and β frequency ranges. The mean ± SD DFA values for phase-shuffled MEG data are also shown for each frequency band (Figure 3). Phase shuffling preserves the spectral power of the signal but disrupts its temporal relationships. The results of shuffled data show DFA exponents close to 0.5 (i.e., uncorrelated noise) in all frequency bands as expected. Importantly, the resting state DFA values are significantly greater (as assessed by permutation p-values, see Section 2.5.4) than the corresponding shuffled values for in the α/µ and β frequency ranges. Significant differences between finger tapping DFA values and the corresponding shuffled values are only observed in the α/µ frequency range.

Comparison of DFA Exponents of Resting MEG with Finger Tapping MEG
The mean ± SD differences between the DFA exponents of rest and finger tapping conditions were calculated for each frequency band (Figure 4, blue bars). The corresponding quantities for the phase-shuffled data were also determined (Figure 4, pink bars). Our results show that the average DFA exponents for the resting state are higher than, or comparable to, those obtained from data recorded during the finger tapping task. When correcting for multiple comparisons within each canonical band (see Section 2.3), the difference between these average DFA exponents is shown to be significant in the 12.5-14.5 Hz (p < 0.05) and 23.5-25.5 Hz (p < 0.01) frequency bands (Figure 4). Only exponents associated with those fluctuation plots judged to be valid by ML-DFA are included in the data. The number of subjects and shuffled time series contributing to the calculation of each average DFA exponent is listed in Table 1. In the two frequency bands where statistically significant differences in DFA between rest and movement were seen, there were 6 (out of 7) and 7 (out of 7) valid subject recordings, respectively. We note that these frequency ranges are close, but not identical, to peak power of the MEG in the α/µ and β frequency ranges, respectively (see for comparison pooled spectra in Figure 2). FIGURE 4 | The mean difference of DFA exponents when contrasting rest and movement in the α/µ, β, and γ frequency ranges (A-C). Summary data for the difference between resting state and finger tapping task are shown in blue (rest-both). Data for the corresponding phase-shuffled data are shown in magenta circles (shuffled rest-shuffled both). Statistical significance (corrected for multiple comparisons) is shown by asterisks (* for p < 0.05; ** for p < 0.01).

Validation of Results
The validity of the >0.5 estimated Hurst exponents detected for the resting state MEG was further assessed for each frequency band through comparing them with (a) synthetic data constructed to yield white Gaussian noise phase difference and (b) empty MEG scanner noise. The number of DFA exponents obtained from analysing phase differences in these control signals was comparable to those of human MEG data in each frequency band ( Table 1). There were significant differences for the resting state data against both MEG scanner noise and Gaussian white noise in the β-band (23.5-25.5 Hz). Statistical significance in the α/µ band was only found against white Gaussian noise in the 10.5-12.5 Hz band and scanner noise in the 12.5-14.5 Hz band. There were no differences in the γ frequency range.

Summary
In this paper, it has been shown that resting state human MEG shows power law scaling in the temporal order of the rate of change of the phase difference between left and right motor cortices. The exponent of this scaling was found to be greater than 0.5, which indicates LRTCs in this measure of synchronization. Furthermore, it was demonstrated that these LRTCs are disrupted by finger movement in the 12.5-14.5 Hz and 23.5-25.5 Hz frequency bands.

Methodological Considerations
We have applied a novel method to MEG data analysis (Botcharova et al., 2014). It combines a technique for determining the rate of change of the phase difference between two time series with detrended fluctuation analysis (DFA). DFA produces an exponent for a power law relationship between (temporal) scale and fluctuation magnitude within that scale. A further aspect of our methodology is the application of ML-DFA (Botcharova et al., 2013) to robustly determine whether a power law exists in the fluctuation plot.
Prior to analysis, the MEG data was beamformed, time series from trials corresponding to the same task for the same subject were stitched together, and the data was then bandpass-filtered. As each of these steps introduces complexities and possible artifacts in the analysis, we discuss these in turn below.

MEG Recording
Cross-talk and noise is a potential problem for MEG recordings (Hari, 2004). The effects of these shortcomings are minimized by identifying the motor cortices individually for each subject using a beamforming technique for β frequencies known to be sensitive to motor tasks Brookes et al., 2011). Here we motivated our choice of weights based on the β band because it has been shown to be sensitive to motor tasks (Pfurtscheller and Aranibar, 1979). This gives us two well-defined regions of interest in the sensori-motor cortex. One could argue that in order to make between-band comparisons a wider bandwidth for weight calculation might have been more appropriate. We do think this could be an area for further investigation, but it is interesting to note that counter-intuitively, the fact that the beamformer was used on the β band makes this band the least likely in which we might find interactions, because all linear covariation between the two sources in this band has been removed.

Filtering
The data was bandpass-filtered in bands of 2 Hz using a technique from Klimesch et al. (1998). This method ensures that the frequency bands are sufficiently small to enable the separation of amplitude and phase components of the signal (Boashash, 1992), yet sufficiently wide to avoid artifactual coupling (Rilling and Flandrin, 2008).

Data Stitching
It was necessary to stitch together data for the robust application of the DFA method. A DFA exponent calculated from a time series stitched from shorter segments with comparable exponents of their own will retain validity and will provide a robust estimate of the underlying Hurst exponent if this exponent is in the LRTC range (>0.5) (Chen et al., 2002). Further, it has been demonstrated that removing 50% of a LRTC signal and stitching together the remaining parts (even with different standard deviations) does not affect the scaling behavior of the signal (Chen et al., 2002).
We note that, since each trial is 30 s in length, any periodic artifactual effects of stitching these trials would occur at a frequency of 1/30 Hz. This frequency is below any frequency band used in our analysis. Furthermore, because DFA operates on detrended data and is insensitive to signal amplitude, differences in the mean between segments or the presence of linear trend are unlikely to affect our method (Hu et al., 2001).

Frequency Band Dependence
We observed that the number of DFA exponents found to be valid by ML-DFA increased with higher frequency bands whereas their value fell as the frequency band increased, a possible effect of filtering (see Figure 3). To eliminate this effect from the interpretation of our results, we presented the DFA exponents results as differences between rest and finger tapping conditions within the same frequency band such that any artifactual inflation of the DFA exponent as a result of the choice of frequency band in which it was measured, was controlled. Importantly, there was no difference, for any frequency band, in the DFA exponents of the rate of change of phase difference between the two phaseshuffled time series (rest and movement). These time series have the same power spectra as the original (not phase-shuffled) data. We also compared our results to those obtained when the analysis was performed on Gaussian white noise and MEG scanner noise filtered identically to the rest and movement MEG data (see Table 2). When correcting for multiple comparisons statistically significant increases in DFA exponent values were detected for resting state data in both α/µ and β frequency ranges.

General Discussion
In a previous paper (Botcharova et al., 2014) we applied this method to phase differences generated by the Kuramoto model with added Gaussian noise. No bandpass filtering was needed as the Kuramoto model produces phase signals. The Kuramoto model has a coupling parameter that can be adjusted. The value of the coupling value at which the Kuramoto system shows a critical transition is known. The critical transition is characterized by a global order parameter which reflects the overall organization of the system. Our methodology makes it possible to make observations at a pair-wise level, i.e., between individual pairs of Kuramoto oscillator. As individual Kuramoto oscillator pairs become fully synchronized, their phase difference no longer contains moment-to-moment fluctuations and thus power law scaling in the DFA measure is lost (Botcharova et al., 2014). This occurs at and above the critical coupling parameter. When coupling is just below this level, power law scaling exists in the DFA of the rate of change of phase difference time series in the pairs approaching full synchrony. The DFA exponent of oscillator pairs that are not yet fully synchronized is above 0.5, indicating LRTCs for fluctuations of their rate of change of phase difference. We note that in the present study similar DFA values (>0.5) were obtained for the resting state MEG recorded between left and right motor cortices and that the maximal DFA values were frequency band specific. In the Kuramato model with noise, as coupling is decreased further, the DFA exponent of the rate of change of pairwise phase differences decreases toward 0.5. The results obtained from analysis of left and right motor cortex MEG during movement showed retained power law scaling but with DFA exponents close to 0.5 and significantly less than those DFA exponent values obtained in the α/µ and β frequency ranges in the resting state. These results suggest that in the resting state, MEG signals behave as weakly coupled oscillators whose phase difference fluctuations scale and have non-trivial temporal order (estimated Hurst exponent > 0.5). This behavior is disrupted by finger movement. The Kuramoto model makes it possible to take a system of oscillators through the transition from an uncoupled to a synchronized state. In the human brain we are only privy to snapshots of the working system of MEG oscillations. However, by focusing on the fluctuations of the rate of change of phase difference between brain areas we can begin to understand the interaction of brain oscillations through comparison of the results of this analysis to those obtained when it is applied to the oscillators pairs within the Kuramoto model.
The MEG phase differences of the two motor cortices at rest showed LRTCs of their rate of change consistent with a state prior to full synchrony. In contrast finger movement was found to disrupt the phase difference fluctuations away from the LRTC regime. The principle of correlations breaking down has been observed previously in the amplitude fluctuations of MEG and EEG (Linkenkaer- Hansen et al., 2004;Athreya et al., 2012). The phase relationship of oscillatory neurophysiological signals is a measure of communication between neuronal pools (Singer, 1999;Pikovsky et al., 2003;Schoffelen et al., 2005;Buzsaki, 2006;Fries, 2009;Akam and Kullmann, 2010). We suggest that LRTCs in the phase relationship of motor cortices show that their interaction at rest is in a state of optimal readiness. During movement, this interaction is disrupted because the cortices become engaged in activity. Our results show that in the active state (finger movement) the cortices display a lower predisposition toward synchrony. This may reflect the brain's need to have independent control over the EMGs of the left and right hands. However, there needs also to be the potential for co-operation at the cortical level between left and right movements which might suggest that the predisposition to synchrony should be higher. In this regard a recent modeling study shows that there are oscillator states in which reduced synchrony results from increased rather than decreased interareal coupling strength (Singh et al., 2013). Application of the techniques described in the present paper to Wilson-Cowan units within different coupling strength regimes (Singh et al., 2013) may clarify the relationship between our MEG results and cortical-cortical coupling strength at rest and during bilateral finger movement.
Our results speak also to the emerging literature on reduction in neuronal temporal variability with task and stimulus (Churchland et al., 2010;He, 2013). He (2011) showed when comparing the resting state to a motor task a reduction in the Hurst exponent calculated from amplitude fluctuations of the fMRI signal. This was interpreted as reflecting an activated system with less temporal redundancy and more suited to online information processing. We see similar reductions in the Hurst exponent for inter-hemispheric MEG phase synchrony fluctuations when comparing the resting state with a motor task. The concept of the brain as a dynamic system whose trajectory is dependent on context and task is supported and this view needs now to include phase synchrony as another variable.
Although our findings are related to previous studies of event related synchronization and desynchronization (ERS and ERD), e.g., (Pfurtscheller and Aranibar, 1979;Pfurtscheller, 1992), our results are distinct in a number of ways. First, we remain agnostic about the synchronization level within any specific neural region, but instead study the phase relationship between two brain regions, which is a property orthogonal to oscillation amplitude measures studied previously (Boashash, 1992). Second, our analysis is applied to time series that last several minutes, during continuous finger tapping rather than addressing the more transient effects of movement previously studied (Pfurtscheller, 1977;Pfurtscheller and da Silva, 1999). Further, it should be noted that time-varying left-right MEG synchrony may be characterized using techniques such as wavelet coherence or optical spectral tracking (Brittain et al., 2009). These detect epochs of consistent frequency and phase between signals albeit for the shortest time scales that can be resolved for a given frequency, with an explicit time-frequency trade-off. These techniques either are not sensitive to, or, in the case of optical spectral tracking (Brittain et al., 2009) (arguably the coherencebased technique with highest temporal resolution), explicitly model out through use of a Kalman filter, moment-to-moment phase difference fluctuations. In contrast our method focuses on the fluctuations of rate of change of phase difference themselves and asks if these fluctuations, when viewed as a time series, scale as a power law with a given Hurst exponent. A direct comparison between these two approaches to MEG-MEG coupling is possible but outside the scope of the present paper.
It is important to note that DFA is dependent on finding power law scaling in the detrended fluctuation magnitude and that power law scaling in itself is not a sufficient condition for criticality (Sornette, 2006). We conclude that the interactions between motor brain regions at rest are indicative of a system in a regime that provides a suitable environment for the brain to react readily to external stimuli. Since we do not have access to the full spectrum of brain activity, we cannot actually say whether it is close to a critical transition. However, the ability of our methodology to identify a system with the potential to change could have important implications in the future for assessing the sensitivity of a human subject in responding to drugs or other treatments.