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.

1. 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 $\frac{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 contra-lateral 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., 2001, 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 moment-to-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 phase-shuffled time series. The demonstration of the presence of scaling with exponent values consistent with the presence of long-range temporal correlations in this data provides an important connection between the critical brain and neural synchrony paradigms.

2. Materials and Methods

A key novel component of this work is our application of a method that allows detection of the presence of long-range 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

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).

2.1. Extraction of the Phase Difference Time Series

2.1.1. 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.

2.1.2. 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.

2.1.3. Phase Difference

The phase relationship ϕ₁(t) − ϕ₂(t) between two different time series s₁(t) and s₂(t) was calculated using the respective Hilbert transform of the signals H[s₁(t)] and H[s₂(t)] (Pikovsky et al., 2003):

$\begin{array}{cc}{\varphi }_1(t)-{\varphi }_2(t)={\text{tan}}^{-1}\left\{\frac{H[s_1(t)]s_2(t)-s_1(t)H[s_2(t)]}{s_1(t)s_2(t)+H[s_1(t)]H[s_2(t)]}\right\}& (1)\end{array}$

The time series ϕ₁(t) − ϕ₂(t) is an unbounded process because ϕ₁(t) and ϕ₂(t) themselves are unbounded as long as the signals s₁(t) and s₂(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 ϕ₁(t) − ϕ₂(t).

2.2. 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).

2.2.1. 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 self-similar, 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.

2.2.2. 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).

2.2.3. 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 self-similar, 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 well-approximated 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 four-segment 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

$p(lns)=\frac{f(lns)}{{\displaystyle {\sum }_{i=1}^{n}f}(lns)}$

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.

2.3. 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⁷ = 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 γ).

2.4. 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³ 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.

2.5. 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.

2.5.1. Synthetic Data with Noise Phase Difference

Two time series s₁(t) and s₂(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 $s_1=\text{cos}[(\omega +\frac{X(t)}{2})t]$ and $s_2=\text{cos}[(\omega -\frac{X(t)}{2})t]$. 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₁(t) and s₂(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₁ and s₂ 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.

2.5.2. 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.

2.5.3. 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).

2.5.4. Phase Shuffling

A permutation¹ 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

$p=\frac{b+1}{m+1}$

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).

3. Results

3.1. 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 2

Figure 2. Pooled power spectra for all time series from 7 subjects. The power spectra of the right and left motor cortices during the resting state and during tapping of both fingers are displayed on the same axes in the left and right plots, respectively. Prominent peaks in power during rest occur at 10.5 Hz and at 21.5 Hz. The dark blue line is associated with the rest condition, and the cyan line indicates the task of tapping both fingers (condition labeled Both in the legend).

3.2. Average DFA Exponents of Rate of Change of Phase Difference Compared to Shuffled Data

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 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.

FIGURE 3

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).

TABLE 1

Table 1. Number of DFA exponents determined to be valid by ML-DFA for each of type of time series (labeling as in Table 2).

3.3. 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

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).

3.4. 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.

4. Discussion

4.1. 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.

4.2. 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.

4.2.1. 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 (Hillebrand and Barnes, 2005; Hillebrand et al., 2005; 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.

4.2.2. 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).

4.2.3. 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).

4.2.4. 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 phase-shuffled 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.

TABLE 2

Table 2. P-values for one-tailed Student's t-tests between the mean difference in DFA exponents between resting state and movement conditions against four sets of control data (white noise - Gauss, scanner noise - Scanner, phase-shuffled rest data - ShuffledR, phase-shuffled tapping data - ShuffledB).

4.3. 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 inter-areal 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 coherence-based 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.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

MB was funded by the Centre for Mathematics and Physics in the Life Sciences and Experimental Biology (CoMPLEX), University College London. SF was supported by University College London Hospitals Biomedical Research Centre (BRC), “Moger Moves” and the Szeben Peto fund. This work was supported by an MRC UK MEG Partnership Grant, MR/K005464/1. The WTCN is supported by the Wellcome Trust.

Supplementary Material

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fphys.2015.00183/abstract

Footnotes

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.

References

Akam, T., and Kullmann, D. M. (2010). Oscillations and filtering networks support flexible routing of information. Neuron 67, 308–320. doi: 10.1016/j.neuron.2010.06.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Andres, F. G., Mima, T., Schulman, A. E., Dichgans, J., Hallett, M., and Gerloff, C. (1999). Functional coupling of human cortical sensorimotor areas during bimanual skill acquisition. Brain 122(Pt 5), 855–870.

PubMed Abstract | Google Scholar

Athreya, D. N., Orden, G. V., and Riley, M. A. (2012). Feedback about isometric force production yields more random variations. Neurosci. Lett. 513, 37–41. doi: 10.1016/j.neulet.2012.02.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Beggs, J. M., and Plenz, D. (2003). Neuronal avalanches in neocortical circuits. J. Neurosci. 23, 11167–11177. doi: 10.3410/f.1018065.217446

PubMed Abstract | CrossRef Full Text | Google Scholar

Beran, J. (1994). “Statistics for long-memory processes,” in Monographs on Statistics and Applied Probability, Vol. 61 (New York, NY: Chapman & Hall).

Berger, H. (1929). Uber das elektrenkephalogramm des menschen. Arch. Psychiatr. Nervenkr. 87, 527–580.

Google Scholar

Boashash, B. (1992). Estimating and interpreting the instantaneous frequency of a signal. I. Fundamentals. Proc. IEEE 80, 520–538.

Google Scholar

Botcharova, M., Farmer, S. F., and Berthouze, L. (2013). A maximum likelihood based technique for validating detrended fluctuation analysis (ML-DFA). Available online at: http://arxiv.org/abs/1306.5075

Botcharova, M., Farmer, S. F., and Berthouze, L. (2014). Markers of criticality in phase synchronisation. Front. Syst. Neurosci. 8:176. doi: 10.3389/fnsys.2014.00176

PubMed Abstract | CrossRef Full Text | Google Scholar

Brittain, J.-S., Catton, C., Conway, B. A., Nielsen, J. B., Jenkinson, N., and Halliday, D. M. (2009). Optimal spectral tracking – with application to speed dependent neural modulation of tibialis anterior during human treadmill walking. J. Neurosci. Meth. 177, 334–347. doi: 10.1016/j.jneumeth.2008.10.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Brookes, M. J., Hale, J. R., Zumer, J. M., Stevenson, C. M., Francis, S. T., Barnes, G. R., et al. (2011). Measuring functional connectivity using MEG: methodology and comparison with fcMRI. Neuroimage 56, 1082–1104. doi: 10.1016/j.neuroimage.2011.02.054

PubMed Abstract | CrossRef Full Text | Google Scholar

Buzsaki, G. (2006). Rhythms of the Brain, 1st Edn. London: Oxford University Press.

Chen, Z., Ivanov, P., Hu, K., and Stanley, H. E. (2002). Effect of nonstationarities on detrended fluctuation analysis. Phys. Rev. E 65(4 Pt 1):041107. doi: 10.1103/physreve.65.041107

PubMed Abstract | CrossRef Full Text | Google Scholar

Chialvo, D. R. (2010). Emergent complex neural dynamics. Nat. Phys. 6, 744–750. doi: 10.1038/nphys1803

CrossRef Full Text | Google Scholar

Churchland, M. M., Yu, B. M., Cunningham, J. P., Sugrue, L. P., Cohen, M. R., Corrado, G. S., et al. (2010). Stimulus onset quenches neural variability: a widespread cortical phenomenon. Nat. Neurosci. 13, 369–378. doi: 10.1038/nn.2501

PubMed Abstract | CrossRef Full Text | Google Scholar

Cohen, D. (1968). Magnetoencephalography: evidence of magnetic fields produced by alpha-rhythm currents. Science 161, 784–786.

PubMed Abstract | Google Scholar

Farmer, S. F. (1998). Rhythmicity, synchronization and binding in human and primate motor systems. J. Physiol. 509, 3–14.

PubMed Abstract | Google Scholar

Fries, P. (2009). Neuronal gamma-band synchronization as a fundamental process in cortical computation. Annu. Rev. Neurosci. 32, 209–224. doi: 10.1146/annurev.neuro.051508.135603

PubMed Abstract | CrossRef Full Text | Google Scholar

Halliday, D., Rosenberg, J., Amjad, A., Breeze, P., Conway, B., and Farmer, S. (1995). A framework for the analysis of mixed time series/point process data - theory and application to the study of physiological tremor, single motor unit discharges and electromyograms. Prog. Biophys. and Mol. Bio. 64, 237–278.

PubMed Abstract | Google Scholar

Halliday, D. M., Conway, B. A., Farmer, S. F., and Rosenberg, J. R. (1998). Using electroencephalography to study functional coupling between cortical activity and electromyograms during voluntary contractions in humans. Neurosci. Lett. 241, 5–8.

PubMed Abstract | Google Scholar

Hansen, P., Kringelbach, M., and Salmelin, R. (2010). MEG: An Introduction to Methods, 1st Edn. New York, NY: Oxford University Press.

Google Scholar

Hari, R. (2004). “Magnetoencephalography in clinical neurophysiological assessment of human cortical functions,” in Electroencephalography: Basic Principles, Clinical Applications, and Related Fields, Chapter 57 (Philadelphia, PA: Lippincott Williams & Wilkins), 1165–1197.

He, B. J. (2011). Scale-free properties of the functional magnetic resonance imaging signal during rest and task. J. Neurosci. 31, 13786–13795. doi: 10.1523/JNEUROSCI.2111-11.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

He, B. J. (2013). Spontaneous and task-evoked brain activity negatively interact. J. Neurosci. 33, 4672–4682. doi: 10.1523/JNEUROSCI.2922-12.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

Hillebrand, A., and Barnes, G. R. (2005). “Beamformer analysis of MEG data,” in Magnetoencephalography, Vol. 68 of International Review of Neurobiology, ed H. Preissl (Amsterdam: Academic Press), 149–171.

Hillebrand, A., Singh, K. D., Holliday, I. E., Furlong, P. L., and Barnes, G. R. (2005). A new approach to neuroimaging with magnetoencephalography. Hum. Brain. Mapp. 25, 199–211. doi: 10.1002/hbm.20102

PubMed Abstract | CrossRef Full Text | Google Scholar

Hosking, J. R. M. (1981). Fractional differencing. Biometrika 68, 165–176.

Google Scholar

Hu, K., Ivanov, P. C., Chen, Z., Carpena, P., and Eugene Stanley, H. (2001). Effect of trends on detrended fluctuation analysis. Phys. Rev. E 64:011114. doi: 10.1103/physreve.64.011114

PubMed Abstract | CrossRef Full Text | Google Scholar

James, L. M., Halliday, D. M., Stephens, J. A., and Farmer, S. F. (2008). On the development of human corticospinal oscillations: age-related changes in EEG-EMG coherence and cumulant. Eur. J. Neurosci. 27, 3369–3379. doi: 10.1111/j.1460-9568.2008.06277.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Kitzbichler, M. G., Smith, M. L., Christensen, S. R., and Bullmore, E. (2009). Broadband criticality of human brain network synchronization. PLoS Comput. Biol. 5:e1000314. doi: 10.1371/journal.pcbi.1000314

PubMed Abstract | CrossRef Full Text | Google Scholar

Klimesch, W., Russegger, H., Doppelmayr, M., and Pachinger, T. (1998). A method for the calculation of induced band power: implications for the significance of brain oscillations. Electroen. Clin. Neuro 108, 123–130.

PubMed Abstract | Google Scholar

Kramer, M. A., Tort, A. B., and Kopell, N. J. (2008). Sharp edge artifacts and spurious coupling in EEG frequency comodulation measures. J. Neurosci. Meth. 170, 352–357. doi: 10.1016/j.jneumeth.2008.01.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Linkenkaer-Hansen, K., Nikouline, V. V., Palva, J. M., and Ilmoniemi, R. J. (2001). Long-range temporal correlations and scaling behavior in human brain oscillations. J. Neurosci. 21, 1370–1377.

PubMed Abstract | Google Scholar

Linkenkaer-Hansen, K., Nikulin, V. V., Palva, J. M., Kaila, K., and Ilmoniemi, R. J. (2004). Stimulus-induced change in long-range temporal correlations and scaling behaviour of sensorimotor oscillations. Eur. J. Neurosci. 19, 203–218. doi: 10.1111/j.1460-9568.2004.03116.x

PubMed Abstract | CrossRef Full Text

Litvak, V., Mattout, J., Kiebel, S., Phillips, C., Henson, R., Kilner, J., et al. (2011). EEG and MEG data analysis in SPM8. Comput. Intell. Neurosci. 2011:852961. doi: 10.1155/2011/852961

PubMed Abstract | CrossRef Full Text | Google Scholar

Nichols, T., and Holmes, A. (2002). Nonparametric permutation tests for functional neuroimaging: a primer with examples. Hum. Brain Mapp. 15, 1–25. doi: 10.1002/hbm.1058

PubMed Abstract | CrossRef Full Text | Google Scholar

Onslow, A., Bogacz, R., and Jones, M. W. (2011). Quantifying phase-amplitude coupling in neuronal network oscillations. Prog. Biophys. Mol. Bio. 105, 49–57. doi: 10.1016/j.pbiomolbio.2010.09.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Peng, C.-K., Buldyrev, S. V., Havlin, S., Simons, M., Stanley, H. E., and Goldberger, A. L. (1994). Mosaic organization of DNA nucleotides. Phys. Rev. E 49, 1685–1689. doi: 10.1103/PhysRevE.49.1685

PubMed Abstract | CrossRef Full Text | Google Scholar

Pfurtscheller, G. (1977). Graphical display and statistical evaluation of event-related desynchronization. Electroen. Clin. Neuro. 43, 757–760.

PubMed Abstract | Google Scholar

Pfurtscheller, G. (1992). Event-related synchronization (ERS): an electrophysiological correlate of cortical areas at rest. Electroen. Clin. Neuro. 83, 62–69.

PubMed Abstract | Google Scholar

Pfurtscheller, G., and Aranibar, A. (1979). Evaluation of event-related desynchronization (ERD) preceding and following voluntary self-paced movement. Electroen. Clin. Neuro. 46, 138–146.

PubMed Abstract | Google Scholar

Pfurtscheller, G., and da Silva, F. H. L. (1999). Event-related EEG/MEG synchronization and desynchronization: basic principles. Clin. Neurophysiol. 110, 1842–1857.

PubMed Abstract | Google Scholar

Phipson, B., and Smyth, G. K. (2010). Permutation p-values should never be zero: calculating exact p-values when permutations are randomly drawn. Stat. Appl. Genet. Molec. 9, 1–16. doi: 10.2202/1544-6115.1585

PubMed Abstract | CrossRef Full Text | Google Scholar

Pikovsky, A., Rosenblum, M., and Kurths, J. (2003). Synchronization: A Universal Concept in Nonlinear Science. Cambridge: Cambridge Nonlinear Science Series; Cambridge University Press.

Google Scholar

Riehle, A., and Vaadia, E. (2004). Motor Cortex in Voluntary Movements: A Distributed System for Distributed Functions. Boca Raton, FL: Taylor & Francis.

Google Scholar

Rilling, G., and Flandrin, P. (2008). One or two frequencies? the empirical mode decomposition answers. Trans. Sig. Proc. 56, 85–95. doi: 10.1109/TSP.2007.906771

CrossRef Full Text | Google Scholar

Schoffelen, J.-M., Oostenveld, R., and Fries, P. (2005). Neuronal coherence as a mechanism of effective corticospinal interaction. Science 308, 111–113. doi: 10.1126/science.1107027

PubMed Abstract | CrossRef Full Text | Google Scholar

Shew, W. L., Yang, H., Petermann, T., Roy, R., and Plenz, D. (2009). Neuronal avalanches imply maximum dynamic range in cortical networks at criticality. J. Neurosci. 29, 15595–15600. doi: 10.1523/JNEUROSCI.3864-09.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Singer, W. (1999). Neuronal synchrony: a versatile code for the definition of relations? Neuron 24, 49–65.

PubMed Abstract | Google Scholar

Singh, R., Menon, S. N., and Sinha, S. (2013). Complex patterns arise through spontaneous symmetry breaking in dense homogeneous networks of neural oscillators. Available online at: http://arxiv.org/abs/1305.7093

Sornette, D. (2006). Critical Phenomena in Natural Sciences: Chaos, Fractals, Self-Organization and Disorder: Concepts and Tools, 2nd Edn. Berlin: Springer.

Taqqu, M. S., Teverovsky, V., and Willinger, W. (1995). Estimators for long-range dependence: an empirical study. Fractals 3, 785–798.

Google Scholar

van Veen, B., van Drongelen, W., Yuchtman, M., and Suzuki, A. (1997). Localization of brain electrical activity via linearly constrained minimum variance spatial filtering. IEEE Trans. Biomed. Eng. 44, 867–880.

PubMed Abstract | Google Scholar

Widmann, A., and Schröger, E. (2012). Filter effects and filter artifacts in the analysis of electrophysiological data. Front. Psychol. 3:233. doi: 10.3389/fpsyg.2012.00233

PubMed Abstract | CrossRef Full Text | Google Scholar