Markers of criticality in phase synchronization

1. Introduction

The concept of the brain as a dynamical system close to a critical regime is attractive because systems close to criticality are thought to maximize their dynamic range of information processing and communication, show efficiency in transmitting information and a readiness to respond to change (Linkenkaer-Hansen et al., 2001, 2004; Beggs and Plenz, 2003; Stam and de Bruin, 2004; Kinouchi and Copelli, 2006; Sornette, 2006; Shew et al., 2009; Werner, 2009; Chialvo, 2010; Beggs and Timme, 2012; Meisel et al., 2012; Shew and Plenz, 2013).

A number of modeling studies have shed important light on the behavior of neurally inspired systems close to their critical dynamical range (Kitzbichler et al., 2009; Shew et al., 2009; Breakspear et al., 2010; Daffertshofer and van Wijk, 2011; Poil et al., 2012). To date there have been two significant experimental observations suggesting that the brain may operate at, or near, criticality. These are: (i) the discovery that the spatio-temporal distribution of spontaneous neural firing statistics can be characterized as neuronal avalanches with a power law distribution of avalanche size (Beggs and Plenz, 2003) and (ii) the presence of long-range temporal correlations (LRTCs) in the amplitude fluctuations of neural oscillations, typically bandpassed MEG or EEG (Linkenkaer-Hansen et al., 2001; Hardstone et al., 2012). The mechanisms by which avalanches and LRTCs of oscillation amplitude may maximize the dynamic range of information processing and communication are still to be fully understood and experimental and computational neuroscience data linking the two phenomena are only just beginning to emerge (Plenz and Chialvo, 2009; Poil et al., 2012).

Population coding approaches to neuronal information storage and transmission show that both changes in the firing rate and changes in neuronal synchronization and desynchronization of action potentials are required to indicate changes in signal salience (Pfurtscheller, 1977, 1992; Singer, 1999; Baker et al., 2001; Schoffelen et al., 2005). At a coarser spatio-temporal scale, extracellular brain signals (local field potentials, corticography, EEG, and MEG), which depend on recordings within the brain, at the brain surface and at the scalp are observed to be quasi-oscillatory (brain oscillations) and in the resting state contain spectral peaks within distinct frequency bands sitting on a 1/f decrease in power with increasing frequency (Buzsaki, 2006). Brain oscillations both in the resting state and during task conditions show short-range and long-range synchronization when examined both from the phase and amplitude envelope perspectives (Wang, 2010). Primarily neuroscience has focused on the detection of synchronization between areas either at zero phase lag, or with a fixed phase delay. This is in part a consequence of the fact that the averaging necessary to extract evidence of signal correlation requires a consistent phase relationship between the two signals for at least some period of the recording.

Importantly, neural synchronization is weak and it fluctuates spontaneously over time. A number of experiments have shown neural synchronization to be consistently modulated by cognitive, perceptual and motor tasks supporting the idea that synchronization and de-synchronization within and across frequency bands may play an important role in communication within the nervous system (Conway et al., 1995; Farmer, 1998; Baker et al., 1999; Singer, 1999; Pikovsky et al., 2003; Schoffelen et al., 2005; Buzsaki, 2006; Doesburg et al., 2009; Fries, 2009; Akam and Kullmann, 2010). Changing synchronization patterns may indicate an evolution in the relationship and exchange of information (Pikovsky et al., 2003). Neural synchronization can exist between nearby and distant regions, across a range of time scales, and can be characterized using a number of techniques based on time- and frequency-domain techniques as well as mutual information (Halliday et al., 1998; Schoffelen et al., 2005; Buzsaki, 2006; James et al., 2008; Brittain et al., 2009; Siegel et al., 2012).

Neuronal synchronization occurs when the mutual influence of neurons on each other causes them to fire close together in time. It is favored by oscillatory activity. Oscillators can be tipped in and out of weak synchonization through shared noise, a phenomenon first appreciated by Huygens (Pikovsky et al., 2003). Therefore, weak yet variable synchrony between neuronal oscillators may easily emerge within complex and highly interactive neural networks. In this paper the term synchronization will be used to encapsulate both zero and fixed phase lag synchrony but also situations in which any non-trivial phase relationship exists between signals. Importantly, we will introduce a new methodology to demonstrate that non-fixed yet non-random phase relationships between signals are present in models of critical synchronization and we will show that, in principle, the methodogy can be applied to neural data in order to further explore the relationship between neural synchronization and systems operating close to a critical regime.

Recent evidence supporting the idea of criticality in the dynamics of the resting state brain activity and the appreciation that synchronization is an important extractable property of neural spatio-temporal dynamics has led researchers to ask whether neuronal synchrony can have properties consistent with a dynamical system at criticality. These approaches identify power law distributions in neural synchronization where synchronization has been defined as phase consistency between two thresholded time series, e.g., see the phase lock interval (PLI) measure and the lability of global synchronization (GLS) measure in Kitzbichler et al. (2009). These findings are of considerable interest, however, the results supporting power law behavior of PLI have been shown by the present authors to be vulnerable to data pooling and therefore may not provide robust estimates of critical synchronization in neural time series data (Botcharova et al., 2012, see also Shriki et al., 2013).

As discussed above, LRTCs (these will be formally defined in Section 2.3) exist in dynamical systems thought to operate close to a critical regime (Linkenkaer-Hansen et al., 2001). They are typically identified by the autocorrelation function of the time series decaying in the form of a power law (Granger and Joyeux, 1980). The detrended fluctuation analysis (DFA) technique allows a characterization of LRTCs through an exponent similar to the Hurst exponent. DFA has been widely used in order to demonstrate the presence of LRTCs in a number of natural and human phenomena (see Peng et al., 1994, 1995a,b; Stanley et al., 1994; Hausdorff et al., 1995; Bak, 1996; Robinson, 2003; Karmeshu and Krishnamachari, 2004; Wang et al., 2005; Samorodnitsky, 2006; Hardstone et al., 2012, for examples). In neurophysiology, the finding of LRTCs in amplitude fluctuations of the bandpass filtered MEG and EEG (Linkenkaer-Hansen et al., 2001, 2004) has inspired us to develop a methodological framework that can be used to to verify the presence or absence of power law scaling of detrended fluctuations and where power law scaling is present to estimate and ascertain non-trivial DFA exponents in the moment to moment fluctuations of phase synchronization (quantified in terms of the rate of change of phase difference time series) between pairs of neuronal oscillation time series. It should be noted here that our focus on the rate of change of phase difference time series means that our framework is not reliant on the definition of (discrete) phase locking events. It is therefore expected to contribute insights regarding phase synchronization that corroborate or complement those provided by the study of intermittency in phase synchronization (e.g., Gong et al., 2007).

The methodology is tested as follows: (i) on synthetic time series where their phase difference has known temporal properties with a known DFA exponent. Using these simulations we demonstrate the method's ability to recover known DFA exponents in the phase difference, and we test the method's robustness to additive noise in such signals; (ii) the method is tested on two classical models of criticality, Ising and Kuramoto (Ising, 1925; Onsager, 1944; Kuramoto, 1975, 1984), from which time series and their pairwise phase differences can be extracted. The output of these models is examined using our method for those parameter values that determine the sub-critical, critical, and super-critical regimes. The classical Kuramoto model is tuned close to the physiological β frequency range of MEG and EEG and examined with additive noise. We show from this analysis that a rise in DFA exponent associated with robust power law detrended fluctuation scaling occurs close to the critical regimes of both the Ising model and the Kuramoto model with noise.

We next use our methodology to examine a system of Kuramoto oscillators, operating in a range of frequencies close to the physiological γ frequency range of MEG and EEG that are connected through a network constructed based on empirical estimations of brain connectivity parameters with time delays, noise and non-uniform connectivity (Cabral et al., 2011). From these simulations, we determine the parameters at which this system shows evidence of LRTCs in the rate of change of phase differences and we relate the presence of LRTCs to the network's connectivity.

Finally, we demonstrate that in principle this methodology may be applied to neurophysiological data through analysing pairs of human EEG and EMG time series. These preliminary results suggest that LRTCs can be detected in the phase synchronization between oscillations in human neurophysiological recordings.

We present and discuss our methodology in detail and we offer an interpretation of its results in relation to the emerging literature on neural synchrony and criticality within neural systems. We suggest that the existence of a valid DFA exponent in fluctuations of a phase difference measure suggests that the fluctuations are governed by non-local behavior, with all scales contributing to system's behavior.

2. Materials and Methods

We seek to characterize the presence of LRTCs in the (time-varying) phase difference between two time series. These time series may be physiological signals such as EEG, MEG, or EMG, time series extracted from a simulation or physical model, or data recorded from other natural phenomena. Below, we present the detail of the various components of our proposed methodology, including a technique used to calculate phase differences, DFA and the recently introduced ML-DFA method for validating the output of DFA. Figure 1 illustrates the application of our methodology to neurophysiological data using two sample MEG time series. We note that for these signals, we bandpassed filter the data to a frequency band of interest, however, this step will be omitted in model data considered further in the manuscript.

FIGURE 1

Figure 1. Step-by-step illustration of the proposed method. We use two sample MEG signals from the left and right motor cortex, displayed throughout panels 1–4 in red and blue, respectively. Panel 2 shows an optional bandpass filtering step. In panel 3 the instantaneous phases of the two time series are calculated using the Hilbert transform. Panel 4 shows the unwrapped phases leading to a time-varying phase difference displayed in panel 5. In panel 6, the rate of change of this phase difference is calculated. This step is illustrated using two plots, each showing a different time scale in the x-axis. These two time scales correspond to the minimum and maximum window sizes used in the DFA analysis, see Section 2.4. Panel 7 shows the resulting DFA fluctuation plot. The validity of this plot is determined using ML-DFA, see Section 2.5. In this case, the validity of the DFA plot was confirmed, with a DFA exponent of 0.57.

2.1. Signal Phase

The phase of a single time series s(t) is calculated by first finding its analytic signal:

$\begin{array}{cc}{s}_a(t)=s(t)+H\left[s(t)\right]& (1)\end{array}$

where H[s(t)] is the Hilbert transform:

$\begin{array}{cc}H[s(t)]=\text{p}.\text{v}.{\displaystyle {\int }_{-\infty }^{\infty }s}(\tau )\frac{1}{\pi (t-\tau )}d\tau & (2)\end{array}$

and p.v. indicates that the transform is defined using the Cauchy principal value.

2.2. Phase Difference

The signal phase is defined such that it belongs to a range ϕ(t) ∈ [0, 2π] or ϕ(t) ∈ [−π, π]. When a single oscillatory cycle is completed the phase returns to its starting value. A time-varying phase therefore has the properties of a sawtooth function (see panel 3 in Figure 1). In order to turn the phase into a continuous signal, the phase is unwrapped, so that at each discontinuity, a value of 2π is added to the phase (Freeman and Rogers, 2002; Freeman, 2004).

The phase difference ϕ₁(t) − ϕ₂(t) between two different time series s₁(t) and s₂(t) is 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\}& (3)\end{array}$

Full synchronization between the two signals is indicated by a constant difference in phase over some time period (Pikovsky et al., 2003). 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 shall use DFA, see Section 2.4, 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 time derivative of the phase difference time series ϕ₁(t) − ϕ₂(t), i.e., the rate of change of the phase difference.

2.3. Long-Range Temporal Correlations

The autocorrelation function R_ss(τ) of a signal s(t) quantifies the correlation of a signal with itself at different time lags τ (Priemer, 1990), formally:

$\begin{array}{cc}{R}_{ss}(\tau )={\displaystyle {\int }_{\infty }^{-\infty }s}(t+\tau )\overline{s}(t)dt& (4)\end{array}$

where $\overline{s}$(t) is the complex conjugate of s(t) and therefore $\overline{s}$ = s(t) if s(t) is real-valued.

In signals with short-range or no dependence (Beran, 1994), the autocorrelation function shows a rapid decay. Gaussian white noise, for example, is a signal with no temporal dependence because each successive value of the time series is independent and thus its autocorrelation function decays exponentially. In contrast, a slow decay of the autocorrelation function indicates that correlations persist even across large temporal separations, and this is referred to as long-range dependence (Beran, 1994).

If there is power law decay of the autocorrelation function, namely:

$\begin{array}{cc}{R}_{ss}(\tau )~C{\tau }^{-\alpha }& (5)\end{array}$

where C > 0 and α ∈ (0, 1) are constants, and the symbol ~ indicates asymptotic equivalence (Clegg, 2006), then the time series is said to contain LRTCs. LRTCs are a subject of considerable scientific interest. They have been detected in biological data (Peng et al., 1994; Carreras et al., 1998; Willinger et al., 1999; Linkenkaer-Hansen et al., 2001; Samorodnitsky, 2006; Berthouze et al., 2010) and have been discussed within the context of complex systems operating in a critical regime.

Applying a Fourier transformation to Equation (5), a similar formulation exists for the spectral density of the signal (Clegg, 2006), with f representing frequency:

$\begin{array}{cc}{G}_{ss}(f)~B{f}^{-\beta }& (6)\end{array}$

where β = 1 − α and is also related to the level of temporal dependence.

The exponents α and β in Equations (5, 6) are connected to the Hurst Exponent, H, by α = 2 − 2H and β = 2H − 1 (Beran, 1994; Taqqu et al., 1995).

In practice, finding the exponent α and β is not straightforward for an arbitrary signal. In the time-domain, α is best approximated by the slope of the autocorrelation function in the limit of infinite time lags τ where measurement errors are also largest (Clegg, 2006). Similarly, in the frequency domain, β is best approximated by the shape of the spectral density at large frequency shifts f. Determination of the Hurst exponent for non-stationary signals is not straightforward, and therefore, for practical applications, the related property of self-similarity (see below) is considered.

2.4. Detrended Fluctuation Analysis

DFA may be used to determine the self-similarity of a time series (Peng et al., 1994, 1995b). The application of DFA returns the value of an exponent, which is closely related to the Hurst exponent (Beran, 1994; Clegg, 2006). DFA is often considered to be applicable to both stationary and non-stationary data although recent reports, e.g., Bryce and Sprague (2012), have suggested that the ability of DFA to deal with non-stationary signals is overstated. In Section 2.5, we will describe our approach to mitigating this concern.

To calculate the DFA exponent, the time series 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 box sizes). 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 box sizes. In the log space, this power law scaling yields a linear relationship between residuals and box 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, and H = 1.5 is Brownian noise. Gaussian white noise has an exponent of H = 0.5.

When performing DFA on oscillatory signals, the smallest window length should be large enough to avoid errors in local root mean square fluctuations, and it is typically taken to be several times the length of a cycle at the characteristic frequency in the time series (Linkenkaer-Hansen et al., 2001). If the minimum window size is significantly smaller than this value, then the fluctuation plot will typically contain a crossover at the window length of a single period (Hu et al., 2001). However, for non-oscillatory time series for which there is no characteristic temporal scale and there are rapid changes at each innovation, such as Gaussian white noise or FARIMA time series (see Section 2.6.1), a smaller window size may be used.

The maximum window size should encompass a significant proportion of the time series yet contain sufficient estimates to allow for a robust estimate of the average fluctuation magnitude across the time series. It is typically taken to be N/10 where N is the length of the data (Linkenkaer-Hansen et al., 2001).

In our application of DFA to neurophysiological and model data, we use 20 window sizes with a logarithmic scaling and a minimum window of 8 time steps for simulated data, and 1 s for neurophysiological oscillations (sampled at 512 Hz, band-pass filtered 15.5–27.5 Hz) providing for a minimum of 16 cycles per second. Following Linkenkaer-Hansen et al. (2001) we take a maximum window size of N/10 time steps where N is the length of the time series.

2.5. Assessing the Validity of DFA

As mentioned above, a self-similar process will produce a power law relationship between the magnitude of the detrended fluctuations and the box sizes. In DFA, this power law scaling is characterized in terms of the linear scaling between the log detrended fluctuations and the log box sizes (DFA fluctuation plot). It is beyond the scope of this paper to argue the validity of operating in the log domain (but see Clauset et al., 2006 for a reasoned view as to why this may not be appropriate), however, since the object of DFA is to find evidence for or against scaling and because a valid DFA exponent can only be obtained when the DFA fluctuation plot is indeed linear we have introduced a model selection method for establishing the linearity of DFA fluctuation plots (Botcharova et al., 2013).

Our arguments for adopting a more rigourous approach are as follows: (i) there is no a priori means of confirming that a signal is self-similar, (ii) a DFA fluctuation plot will necessarily increase with window size, (iii) an exponent may be too easily obtained through simple regression analysis producing a statistically significant result with a high r² value even though the linear model may not best represent a given DFA fluctuation plot, (iv) the discovery of an exponent >0.5 with a high r² value may lead to the incorrect conclusion that the signal is self-similar with LRTCs.

Instead of a simple regression we use the model selection technique (ML-DFA) introduced in Botcharova et al. (2013) to determine whether a given DFA fluctuation plot is best-approximated by a linear model. This is a heuristic technique, which has been tested extensively and found to perform well in assessing linearity in the fluctuation plots of the following time series: (i) those with known combinations of short and LRTCs, (ii) self-similar time series with varying Hurst exponent, (iii) self-similar time series with added noise and (iv) time series with known oscillatory structure, e.g., sine waves (Botcharova et al., 2013).

The technique fits the DFA fluctuation plot with a number of different models (see below) and compares the fit of each model 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 is linear. We want to stress that this does not equate to stating that the fluctuation plot is linear. Rather, we do not reject the linear model hypothesis. In what follows, only those time series for which the linear model hypothesis is not rejected (i.e., their DFA fluctuation plot is best-fitted by the linear model) contribute to the DFA exponents presented in the present paper and where appropriate we indicate where linear scaling of the fluctuation plot is lost.

The models included in ML-DFA are listed below (see Botcharova et al., 2013 for a justification), with the a_i parameters to be found. The number of parameters ranges between 2 for the linear model, and 8 for the four-segment spline model.

Polynomial - $f(x)={\displaystyle {\sum }_{i=0}^K{a}_i}{x}^i$ for K = {1, …, 5}

Root - f(x) = a₁(x + a₂)^1/K + a₃ for K = {2, 3, 4}

Logarithmic - f(x) = a₁log(x + a₂) + a₃

Exponential - f(x) = a₁e^a₂x + a₃

Spline with 2, 3 and 4 linear sections.

The first step of ML-DFA is to normalize the fluctuation magnitudes with:

$l{F}_{scaled}=100\times \frac{lF-l{F}_{min}}{l{F}_{max}-l{F}_{min}}$

where lF_min and lF_max are the minimum and the maximum values of vector lF, respectively. A function is then defined:

which is a product across all windows i, and which works in a similar way to a likelihood function, where p(lns) represents the function:

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

where f(lns) is the fitted model. Absolute values are used in order to ensure that p(lns) remains in the range [0, 1], so that a function is rejected if it falls below 0.

The next step is to apply a logarithm to to produce a function that is similar in form to a log-likelihood:

This is maximized to find the parameters a_i necessary for f(lns). It is worth mentioning that the application of the logarithm means that the values belonging to lns are not equally weighted for all i. The larger window sizes have a lower weighting, which is beneficial because these estimates are also the least robust since they have fewer samples associated with them.

Akaike's Information Criterion (AIC) is then computed, which is designed to prevent over-fitting—a situation that should in general be avoided—by taking into account the number of parameters used (Akaike, 1974; Mackay, 2003). For a model using k parameters, with likelihood function log, the AIC is calculated using the following expression:

where k is the number of parameters that the model uses (Akaike, 1974). An adapted formula was proposed by Hurvich and Tsai (1989), which accounts for small sample sizes. The model which provides the best fit to the data is that with the lowest value of AIC. It is important to recall that the AIC can only be used to compare models. It does not give any information as to how good the models are at fitting the data, i.e., it is only its relative value, for different models, that is important; and it would not be possible, for instance, to compare AIC values obtained from different data sets to each other.

2.6. Method Validation

2.6.1. FARIMA processes

An Autoregressive Fractionally Integrated Moving Average model (FARIMA) (Hosking, 1981) can be used to create time series with self-similarity. The model provides a process that can easily be manipulated to include a variable level of LRTCs within a signal, from which DFA should return the exponent used to construct the FARIMA process.

To construct a FARIMA process a time sequence of zero-mean white noise is first generated, which is typically taken to be Gaussian, and necessarily so to produce fractional Gaussian noise. The FARIMA process, X(t), is then defined by parameters p, d, and q and given by:

$\begin{array}{cc}\left(1-{\displaystyle \sum _{i=1}^p{\phi }_i{B}^i}\right){\left(1-B\right)}^{d}X(t)=\left(1+{\displaystyle \sum _{i=1}^q{\phi }_i{B}^i}\right)\epsilon (t)& (7)\end{array}$

B is the backshift operator operator, so that BX(t) = X(t − 1) and B²X(t) = X(t − 2). Terms such as (1 − B)² are calculated using ordinary expansion, so that (1 − B)²X(t) = X(t) − 2X(t − 1) + X(t − 2). While the parameter d must be an integer in the ARIMA model, the FARIMA can take fractional values for d. A binomial series expansion is used to calculate the result:

${(1-B)}^d={\displaystyle \sum _{k=0}^{\infty }\left(\begin{array}{c}d\\ k\end{array}\right)}{(-B)}^k$

The left hand sum deals with the autoregressive part of the model where p indicates the number of back-shifted terms of X(t) to be included, φ_i are the coefficients with which these terms are weighted. The right hand sum represents the moving average part of the model. The number of terms of white noise to be included are q, with coefficients φ_i. In the range |d| < $\frac{1}{2}$, FARIMA processes are capable of modeling long-term persistence (Hosking, 1981). As we will only consider p = 1 and q = 1 throughout the manuscript, we will refer to φ₁ as φ and φ₁ as θ. We set |φ| < 1, |θ| < 1 to ensure that the coefficients in Equation (7) decrease with increasing application of the backshift operator, thereby guaranteeing that the series converges, and X(t) is finite (Hosking, 1981).

A FARIMA(0,d,0) is equivalent to fractional Gaussian noise with d = H − $\frac{1}{2}$ (Hosking, 1981). This produces a time series with a DFA fluctuation plot that has been shown to be asymptotically linear with a slope of d + 0.5 (Taqqu et al., 1995; Bardet and Kammoun, 2008). By manipulating the φ and θ parameters, the DFA fluctuation plots can also be distorted.

2.6.2. Surrogate data

Two time series x₁(t) and x₂(t) can be constructed such that the time derivative of their phase difference is a FARIMA time series X(t) with a known DFA exponent (Hosking, 1981). Concretely, we work backwards from the time series X(t) to which DFA is applied. The phase difference of the two time series Δ(ϕ(t)) will be the cumulative sum of X(t), which is discrete in this case:

$\Delta (\varphi (t))={\displaystyle \sum _{s=1}^{t}X}(s)$

The two phases ϕ_i(t) and ϕ₂(t) of x₁(t) and x₂(t), respectively, must be constructed to have a difference of Δ(ϕ(t)), or some multiple of Δ(ϕ(t)) since DFA is unaffected by multiplying a time series by a constant. We therefore set ${\varphi }_1(t)=\frac{{\displaystyle {\sum }_{s=1}^{t}X}(s)}{2f_s}$ and ${\varphi }_2(t)=-\frac{{\displaystyle {\sum }_{s=1}^{t}X}(s)}{2f_s}$ where f_s takes the role of a nominal sampling rate for the surrogate data.

Since the phase of a cosine signal is equal to its argument, the two signals x₁(t) and x₂(t) are defined as:

$x_1=\text{cos}\left(\omega +\frac{{\displaystyle {\sum }_{s=1}^{t}X}(s)}{2f_s}\right)$

and

$x_2=\text{cos}\left(\omega -\frac{{\displaystyle {\sum }_{s=1}^{t}X}(s)}{2f_s}\right)$

where ω is a constant.

In what follows, we used ω = 1 and f_s = 600. These values were chosen in order to produce a smooth enough phase difference. This was necessary to prevent artifacts produced by the Hilbert transform when applied to non-smooth data. When using physiological data, a high enough sampling rate guarantees that the signals will be smooth.

A hundred time series X(t) were generated using the algorithm described in Hosking (1981) for each of the 11 DFA exponents 0.5, 0.55, 0.6, …, 1. Each simulation contains 2²² = 4194304 innovations. The value of the exponent of X(t) is first computed, the two signals x₁(t) and x₂(t) are then constructed, and the phase analysis method is applied. Window sizes used for application of DFA were logarithmically spaced with a minimum of 600 time steps to correspond to f_s and maximum N/10 where N = 2²² is the length of the time series.

A further control analysis was performed in which a Gaussian white noise time series η_i(t) was added to one of the signals, namely,

${x^{\prime }}_1(t)=\text{cos}(\omega +\frac{{\displaystyle {\sum }_{s=1}^{t}X}(s)}{2})+{\eta }_i(t)$

before the phase analysis method was applied in order to recover the DFA exponent of the phase difference X(t). This allowed us to alter the signal-to-noise ratio of x₁(t) in an additive way, which we may suppose to be the case for noise in a neurophysiological time series. By applying the phase analysis method to signals with additive noise, we were able to test the robustness of the method to noisy data. In this analysis, first we will estimate the extent to which the DFA exponent alters when noise is added. Second, we will assess whether ML-DFA rejects those DFA exponents that we know to contain noise, and if so, we will quantify the level of noise at which exponents are no longer valid.

2.7. Model Simulations

2.7.1. The Ising model

The Ising model is a model of ferromagnetism (Ising, 1925). In two dimensions, the model is implemented on a lattice (grid) of elements, or particles which represent a metallic sheet. A temperature parameter controls the collective magnetization (Onsager, 1944). The Ising model has been recently used as a model for a two-dimensional network of connected and interacting neurons (Kitzbichler et al., 2009).

Each element of the grid is assigned a spin p_i, initially at random, which takes a value +1 (spin up) or −1 (spin down). Spins may switch up and down in time in a fashion influenced by both the energy of the full system and by the spin configuration of other neighboring elements. The energy of the system in a given configuration of spins p is given by the Hamiltonian function $H(p)=-J{\Sigma }_{i,j=nn(i)}^N{p}_i{p}_j$, where j is an index for the four elements that are nearest neighbors nn of each element, i of the square grid. The negative sign is included by convention. The average energy of the system E = < H > where the symbol <> indicates taking the expectation value.

The probability ℙ of a given configuration occurring is then proportional to ℙ = e^−H(p)/kT, where T is the temperature parameter and k is Boltzmann's constant. The system may switch into a new configuration if its associated probability is higher or equal to that of the current configuration. The Ising model is implemented using the Metropolis Monte Carlo Algorithm (Metropolis et al., 1953).

At temperature T = 0, the system is highly ordered and corresponds to a magnetic state (see Figure 2 for an example of an Ising model lattice). With increasing temperature values, the probability of a spin changing increases. As the system temperature increases the spins change more rapidly and the system becomes increasingly disordered and corresponds to a non-magnetic state (Figure 2A). The temperature value at which a transition occurs between the magnetized and non-magnetized states is known as the critical temperature T_c. At this temperature (see Figure 2B), the system will have a large dynamic range and infinite correlation length. However, in practice, this means that the system contains spin clusters of all sizes, and correlations between elements of an infinite system remain finite (Onsager, 1944; Daido, 1989). In other words, the Ising model is predicted to have long-range correlations between its elements at T_c.

FIGURE 2

Figure 2. The Ising model lattice at a single time point once steady state has been reached for 3 different values of the temperature parameter. (A) The Ising lattice at a cold temperature of 1.5. Almost all spins are aligned (white) and there is little change across time. (C) The Ising lattice at a high temperature of T = 10⁵. The spins form a more or less random pattern across the lattice. (B) The Ising lattice near critical temperature, T = 2.3. The lattice contains clusters of spins that are both small and large. Note that these are snapshots and that the spin structure of the model is best appreciated when evolving across time.

The value of the critical temperature T_c was calculated for the two-dimensional Ising model in Onsager (1944), and is given by the solution to the equation

$\text{sinh}\left(\frac{2J}{k{T}_c}\right)=1$

In the implementation of the Ising model used here, the lattice consists of 96 × 96 elements. The constants J and k are set to J = 1 and k = 1 without loss of generality, which gives the critical temperature $T_c=\frac{2}{ln(1+\sqrt{2})}\approx 2.269$.

In order to obtain a time series from this spatial model, we follow the procedure introduced by Kitzbichler et al. (2009). Namely, the lattice is divided into a number of smaller square lattices, which we refer to as sub-lattices, and a number of time series are created by taking an average spin value for each sub-lattice. Here, we use a sub-lattice size of 8 × 8 as in Kitzbichler et al. (2009), but we also investigated other sub-lattice sizes (results not shown) in order to verify that this choice of sub-lattice size did not affect the results. Indeed, previous work by Priesemann et al. (2009) suggests that the sub-sampling of a system may cause it to be mis-classified as sub-critical or supercritical when it is in fact in a critical state.

Pairs of time series, for every possible pairing of sub-lattices belonging to the larger grid, were used as input signals for the phase analysis method. For the sub-lattice of size 8 × 8 considered here, 144 time series could be created allowing for 10,296 pairings. Each time series consisted of 64,000 innovations.

2.7.2. The Kuramoto model

The Kuramoto model is a classical model of synchronization (Acebrón et al., 2005; Chopra and Spong, 2005) and has been used to study the oscillatory behavior of neuronal firing (Pikovsky et al., 2003; Kitzbichler et al., 2009; Breakspear et al., 2010) among many other biological systems.

The Kuramoto model describes the phase behavior of a system of mutually coupled oscillators with a set of differential equations. Each of N oscillators in the system rotates at its own natural frequency {ω_i, i = 1, …, N}, drawn from some distribution g(ω). However, it is attracted out of this cycle through coupling K, which is globally applied to the system. Time t is taken to run for T seconds of length dt = 10⁻³. The differential equation to describe the phase of an oscillator is (Kuramoto, 1975, 1984):

$\begin{array}{cc}{\dot{\varphi }}_i(t)={\omega }_i(t)+\frac{K}{N}{\Sigma }_{j=1}^N\text{sin}({\varphi }_j(t)-{\varphi }_i(t))& (8)\end{array}$

Because the Kuramoto model provides an equation governing the phase evolution of each oscillator in the system, there is no need for the Hilbert transform to recover the phase time series and therefore only the latter stages of the phase analysis method are used (see steps 3–6 in Figure 1).

Kuramoto (1975) showed that the evolution of any phase ϕ_i(t) may be re-expressed using two mean field parameters, which result from the combined effect of all oscillators in the system. Namely, we may write:

$\begin{array}{cc}{\dot{\varphi }}_i(t)={\omega }_i+Kr(t)\text{sin}(\psi (t)-{\varphi }_i(t))& (9)\end{array}$

where ψ(t) is the mean phase of the oscillators, and r(t) is their phase coherence, so that:

$\begin{array}{cc}r(t)e^{i\psi (t)}=\frac{1}{N}{\displaystyle \sum _{j=1}^N{e}^{i{\varphi }_j(t)}}& (10)\end{array}$

This crucially indicates that each oscillator is coupled to the others through its relationship with mean field parameters r(t) and ψ(t), so that no single oscillator, or oscillator pair drives the process on their own. The oscillators synchronize at a phase equal to the mean field ψ(t), and r(t) describes the strength of synchronization, sometimes referred to as the extent of order in the system (Strogatz and Mirollo, 1991; Bonilla et al., 1992). When r(t) = 0, no oscillators are synchronized with each other. When r(t) = 1, all oscillators are entrained with each other.

One solution to Equation (9) is r ≡ 0 for all time and coupling, leaving each oscillator to evolve independently at its own natural frequency. Using a limit of N → ∞, some further deductions can be made, including the fact that when the natural frequency distribution g(ω) is unimodal and symmetric, another solution can be found for ω_i, with r(t) not equivalent to 0 (Kuramoto, 1975). A critical bifurcation occurs for sufficiently high coupling, resembling a second-order phase transition (Miritello et al., 2009) in which the order parameter [here, r(t)] leaves zero and grows continuously with coupling (Strogatz and Mirollo, 1991; Dörfler and Bullo, 2011). The coupling at the bifurcation is referred to as the critical coupling K_c (Dörfler and Bullo, 2011).

In an infinite Kuramoto model, criticality is defined through this point of bifurcation. For a finite system, however, the critical point can only be approximated by this theoretical value. One defining characteristic of the critical coupling for the Kuramoto system is that the greatest number of oscillators come into synchronization at this value. In our study, we deal with finite-sized implementations of the Kuramoto model, and we use this characteristic as a marker of the onset of critical regime in addition to the theoretical value K_c. Specifically, we use a measure characterizing the onset of synchronization with increasing coupling introduced by Kitzbichler et al. (2009). This is the change in the “effective mean-field coupling strength,” Δ(Kr). If the value of Kr exceeds the difference between the natural frequency and the mean phase ω_i − ψ (in modulus), i.e., |ω_i − ψ| < Kr, then oscillator i will synchronize to the mean field (Mertens, 2011). Thus, the value of K at which Kr increases maximally is the coupling value at which the greatest number of oscillators are drawn into the mean field.

In this paper, we consider the Kuramoto model with a noise term added to the phase equation, namely, Equation (8) becomes:

$\begin{array}{cc}{\dot{\varphi }}_i(t)={\omega }_i(t)+\frac{K}{N}{\displaystyle \sum _{j=1}^N\text{sin}}({\varphi }_j(t)-{\varphi }_i(t))+{\eta }_i(t)& (11)\end{array}$

where η_i is a noise input taken to be uncorrelated Gaussian noise with zero mean (〈η_i〉 = 0) and covariance σ²_i/T (〈η_i(t)〉〈η_j(s)〉 = δ_ijδ(t − s)σ²_i/T) where δ_ij is the Kronecker delta, δ(t − s) is the Dirac delta function, σ_i is in radians and T = 1 s here.

This creates a richer structure in the oscillator dynamics, which we suggest may better reflect coupling of neurophysiological oscillators. Furthermore, it has been shown that addition of noise increases the critical regime over a wider range of coupling values (Breakspear et al., 2010). This may allow for the fluctuations of phase difference of a given oscillator pair to persist for longer with increasing coupling before full synchronization is achieved.

Strogatz and Mirollo (1991) analytically derived a formula for the critical coupling in an infinite Kuramoto model with added noise K_c,noise. As the number of oscillators is inevitably finite, this value is only an approximation to the true critical coupling in the system, but we find it useful and it is displayed alongside plots of Δ(Kr), which although originally introduced for a noiseless model, remains a helpful marker of the effective critical coupling in the Kuramoto model when noise levels are not too large (Mertens, 2011).

In this study, we generated time series for 200 oscillators of the Kuramoto model described by Equation (11). Each time series was 6100-timestep long. The standard deviation σ_i was set to 0.32. The distribution of natural frequencies was g(ω) ~ (44π, σ_ω), with standard deviation σ_ω = 15. This corresponds to a normal distribution centered around 22 Hz (which is a unimodal distribution). In order to get an idea of the spread of the distribution, the minimum natural frequency selected from this distribution was 16.3 Hz and the maximum was 27.8 Hz. We selected this frequency range because it spans the β-band of EEG, MEG, and EMG oscillations (Farmer, 1998).

For these parameter values, the critical coupling K_c is equal to:

$K_c=\frac{2\sqrt{2}}{\sqrt{\pi }}{\sigma }_{\omega }~23.93$

The integral for K_c,noise is not analytically calculable for a normal distribution g(ω) ~ , but empirical calculation yields:

$K_{c,noise}~23.85$

2.7.3. The Cabral model

The third model that we consider in this paper was developed by Joanna Cabral and her colleagues, referred to as the Cabral model. It is a modification of the Kuramoto model, combining the dynamics of the Kuramoto oscillators with the network properties observed in the human brain (Cabral et al., 2011).

The Cabral model includes a noise input to the Kuramoto oscillators and situates the 66 oscillators on a connectivity matrix with varying connection strengths and time delays based on empirical measurements of 998 brain regions, which have been down-sampled to 66 (Honey et al., 2009). The list of brain regions considered in this model are given in the supplementary material of Cabral et al. (2011) and are reproduced in the Appendix to the present paper. Specifically, Equation (8) is modified to include a connectivity term C_ij between oscillators j and i, namely,

$\begin{array}{cc}{\dot{\varphi }}_i(t)={\omega }_i(t)+\frac{K}{N}{\Sigma }_{j=1}^N{C}_{ij}\text{sin}({\varphi }_j(t-D_{ij})-{\varphi }_i(t))+{\eta }_i(t)& (12)\end{array}$

where η_i is the noise input previously introduced, and D_ij is the time delay associated with the link between oscillators j and i. The matrix of delays D is extracted from a matrix of empirical distances L between regions using:

$D_{ij}=\frac{\langle D\rangle L_{ij}}{\langle L\rangle }$

and is used to encode the length of time taken by neural activity to traverse the connection space. The connectivity and distance matrices (C and D, respectively) are shown in Figure 12. They can also be visualized through the schematic diagram in Figure 3 in which the thickness and color of the lines represent the weights of the connections between the oscillators denoting individual brain regions. These weights are proportional to the number of fibers that were empirically observed to connect the various regions (Cabral et al., 2011, 2012). Brain regions may be identified by their labels, the abbreviations of which are given in Table A1 in the Appendix.

FIGURE 3

Figure 3. Schematic plot (top view) of the Cabral human brain model showing the connections and connection weights between oscillators which correspond to different brain regions. The weight of the connection lines represent the strength of connectivity between the oscillators. The darkest blue lines are the strongest 1% of connections. The node colors represent oscillators, which model different brain regions as detailed in Cabral et al. (2011). Colors are consistent for homologous regions in the left and right hemispheres. Anterior and posterior, left and right are shown.

In Cabral et al. (2011), the model was used to generate time series which were used as input to a hemodynamic model and bandpass filtered. Each time series was 10⁶ timestep-long, corresponding to 1000 s. The resulting time series were compared to recordings of BOLD fMRI signals using Pearson's correlation coefficient and mean squared error to determine the parameter values K and 〈D〉 that generated the time series which most closely approximated the BOLD data.

In this model, there is no theoretically derived value of critical coupling and Δ(Kr) is only a marker of effective change in coupling that may or may not be critical. We interpret a rise in Δ(Kr) as an increase in order of the system similar to that observed by Kitzbichler et al. (2009).

The phase analysis method presented here was applied to the Cabral model for coupling parameters K ranging from 1 to 20. We note that this encompasses K = 18, the value identified by Cabral et al. (2011) as best approximating human brain resting state BOLD fluctuations. Natural frequencies were drawn from a normal distribution with g(ω) ~ (120π, σ_ω) with standard deviation σ_ω = 5, which corresponds to a normal distribution centered around 60 Hz in the γ frequency band. This was selected because γ oscillations have been shown to play a significant part in the BOLD signal fluctuations (see Cabral et al., 2011 for details).

The standard deviation σ_i of the noise input was set to 1.25. It was found that values of σ_i < 3 did not significantly alter the resulting parameter values of K and 〈D〉. The value 〈D〉 = 11 is taken as in Cabral et al. (2011).

2.7.4. Clusters in the Cabral model

Cabral et al. (2011) identified a number of clusters of oscillators, along with a set of 12 oscillators which are not part of a cluster. These clusters are listed below in Table 1. In our analysis, we considered how each of these different clusters contributed to the overall behavior.

TABLE 1

Table 1. Cluster information.

2.7.5. Disruptions to the Cabral model

In order to investigate the role of connectivity in sustaining LRTCs of rate of change of phase difference, we modified the connectivity matrix C in the Cabral model in two ways, as shown in Figure 4. First, beginning with the empirical connectivity matrix we deleted any connection that extended from one hemisphere into the other. We preserved all the other elements of the model's connectivity and oscillator characteristics.

FIGURE 4

Figure 4. Schematic plot (top view) showing the connections and connection weights between oscillators belonging to two modifications to the connectivity of the Cabral human brain model. (A) The left and right hemispheres of the brain have been disconnected, but connections within each hemisphere are left unchanged. (B) The connections and weights of each node are assigned randomly, but the degree distribution and weight distribution at each node is kept constant. The weight of the connection lines represent the strength of connectivity between the oscillators. The darkest blue lines are the strongest 1% of connections. The node colors represent oscillators, which model different brain regions as detailed in Cabral et al. (2011) and are identical to Figure 3. Colors are consistent for homologous regions in the left and right hemispheres. Anterior and posterior, left and right are shown.

The second exploration involved a reconnection of the connectivity matrix in a random arrangement, while preserving the degree distribution and weight distribution of each oscillator by an algorithm described in Gionis et al. (2007), Hanhijärvi et al. (2009). Specifically, a list of the outgoing weights of each oscillator was made alongside the node from which it extends. Two weights were selected from this list. If they did not belong to the same node, then the nodes were connected to each other with the associated outgoing weights that were selected. These weights were then deleted from the list. To continue the algorithm, two further weights were selected. After the first step, it was necessary to check at each iteration that the nodes were not already connected before connecting them. If the nodes were connected, or if they were the same node, new weights were selected from the list.

Analysis of the random connectivity model and comparison of the results obtained from it to those derived from the disconnected hemisphere model and standard appropriately connected model allowed us to determine the extent to which a realistic connectivity matrix of the human brain predisposes the system to LRTCs in the rate of change of the phase difference between the oscillator pairs representing different brain regions.

2.7.5.1. A note on notation. From this point in the text, all instances of oscillator phase ϕ_i(t) and r(t) will be written as ϕ_i and r for ease of notation, unless stated otherwise. Any quantities that are defined using the phases of one or more oscillators are also implicitly functions of time, although the t is omitted for the same reason.

2.8. Neurophysiological Data

Previously collected neurophysiological data were used to illustrate the application of the method (see James et al., 2008 for full details). Briefly, EEG and EMG signals were simultaneously recorded whilst a healthy adult subject performed a 2-min 10% MVC (maximum voluntary contraction) isometric abduction of the index finger of the right hand. The EMG was recorded using bipolar electrodes situated over the first dorsal interosseous muscle (1DI). The EEG was recorded using a modified Maudsley montage from 24 Ag/AgCl electrodes with impedance <5 kΩ. The data were amplified and bandpass filtered 4–256 Hz and sampled at 512 Hz. We analyzed EEG recorded from over the left sensorimotor cortex. The signal processing pathway was set out as in Figure 1, including bandpass filtering in the β frequency range (15.5–27.5 Hz).

3. Results

3.1. Surrogate Data

The signals described in Section 2.6.2 were analyzed. The scatter plot presented in Figure 5 shows the DFA exponents of the rate of change of phase difference expected from the construction of a FARIMA time series with known parameters against those recovered by applying the phase analysis method. The scatter plot shows a strong linear relationship between the expected and recovered exponents with a slope of 0.998. The fact that the slope is slightly <1 indicates that the recovered exponent was slightly under-estimated by our method. This minor tendency will decrease the likelihood of false positive results.

FIGURE 5

Figure 5. Plot of the recovered against the true DFA exponent for FARIMA time series. The relationship between recovered and true DFA values is well-approximated by a linear trend with a slope of almost 1. The error bars increase very slightly with increasing DFA exponent.

As noise is added to a signal with a known DFA exponent in its phase, the exponent of its phase is found to be reduced. Figure 6 shows that as the noise level is progressively increased, the percentage difference between the known DFA exponent and that recovered by the method increases. When the noise level is above one which causes the percentage difference between known and recovered DFA exponent to exceed approximately 5% (note, as shown in Figure 6, that this noise level depends on the exponent, e.g., 0.1 for true DFA exponent of 1, 0.025 for exponent of 0.75), no values are returned for the recovered DFA exponent. This occurs because the recovered DFA exponents are not considered to be valid by ML-DFA because their associated DFA fluctuation plots are not best approximated by a linear model (see Section 2.5).

FIGURE 6

Figure 6. True and recovered DFA exponents for noisy signals with LRTCs. (A) Recovered DFA exponent values as noise is progressively added. For each of the DFA exponents given in the legend (box insert), a signal x′₁(t) was constructed with a noise level σ ∈ [0, 1], shown on the x axis. The phase synchrony analysis method was applied to x′₁(t) and x₂(t). This was performed 100 times. For DFA exponents corresponding to DFA fluctuation plots that were accepted as linear by ML-DFA, the average value for the 100 signal pairs is shown. There are no data points corresponding to the intermediate noise level of ≈0.1 to ≈0.3 because all 100 DFA fluctuation plots for signals with this noise level were determined to be invalid by ML-DFA. (B) The % difference between recovered and known DFA exponents as a function of the noise added to a signal with a known DFA exponent in its phase. The data shown in this plot is the same as that in (A), but it is expressed in terms of the % difference between true and recovered DFA exponents rather than the raw recovered value. Only noise levels of σ ∈ [0, 0.1] are shown. The colors represent different true DFA exponent values, as indicated by the legend within the inserted box. The dashed line indicates a 5% difference between known and recovered exponents. When the difference between the known and recovered exponent exceeded approximately 5% for any value of the true exponent, the DFA fluctuation plot is not accepted as being linear by ML-DFA and therefore the exponent is not shown on the plots.

As the noise level is increased further, and as it passes a level of ≈0.3–0.4, noise dominates the signal and valid exponents are once again obtained. These exponents are at or close to 0.5 regardless of the value of the known DFA exponents, indicating that the phase relationship of the two signals s₁(t) and s₂(t) is dominated by noise only.

3.2. The Ising Model

Figure 7 shows the results for sub-lattices of size 8 × 8. At a high temperature of T = 10⁵, the average DFA exponent across all pairwise comparisons is 0.57 (see magenta shaded bar). This value is in excess of 0.5 expected for Gaussian white noise and indicates that even at high temperatures there is order within the rate of change of phase difference between pairs of lattice time series. As the temperature is lowered the DFA exponent of the rate of change of phase difference increases steadily reaching a maximum of 0.65 at T = 2.55 (see magenta shaded bar) indicating maximal LRTC just before the critical temperature is reached.

FIGURE 7

Figure 7. Average DFA exponents of rate of change of phase difference between pairs of time series generated by 8 × 8 sub-lattices of the 96 × 96 Ising model lattice. The temperature parameter, T, is varied on the x axis. The average of the valid DFA exponents is shown in pink, and the error bars are a single standard deviation from the mean. The proportion of valid exponents, as calculated by ML-DFA, is denoted by the vertical bars. The theoretical critical parameter T_c is indicated by a red asterisk. A horizontal dashed line at DFA exponent 0.5 is plotted to guide the eye. Validity bars that are referred to in the text are highlighted in magenta.

The change in mean DFA has to be seen within the context of the validity of the DFA fluctuation plots. As the system cools toward the critical point the validity of DFA exponents across all pairwise phase differences drops abruptly. The first temperature value for which <100% of the DFA plots are valid is T = 2.75 shown as magenta shaded bar. There is a large fall in DFA fluctuation plot validity as the critical temperature is reached (56–34%). This fall in validity reflects the onset of full synchronization between a number of the time series. At the critical point, T = Tc which occurs between T = 2.25 and T = 2.3 (see magenta shaded bars) the validity is 34% of time series pairs with mean DFA exponent of 0.64. As the Ising model cools below the critical point the DFA validity in general falls and there are no valid DFA fluctuation plots below T = 2.15. As discussed above this occurs because of the loss of fluctuations in the rate of change of phase difference due to full synchronization.

Results obtained for sub-lattice sizes of 32 × 32, 16 × 16, 12 × 12, and 6 × 6 were found to be qualitatively consistent with the results shown in Figure 7 (results not shown).

3.3. The Kuramoto Model

The group average results for the Kuramoto model are shown in Figure 8. As can be seen, the peak average DFA exponent occurs on average at K ≈ 22. The value of the average DFA exponent at this coupling value is 0.65 with standard deviation 0.06, consistent with the rate of change of phase difference showing LRTCs. The peak DFA exponent occurs one coupling value later than the peak of the Δ(Kr) measure, at K ≈ 21. Δ(Kr) represents the coupling value at which the order parameter r increases most, and the point of greatest oscillator coupling flux in the system (Kitzbichler et al., 2009). The peak coupling value Δ(Kr) and the maximum DFA values are just less than the theoretical critical coupling of the infinite Kuramoto system with noise K_c ≈ 23.85. Again, these results must be understood in context of DFA fluctuation plot validity which is 42% of the 199,000 oscillator pairs at K ≈ 22. Once full synchonization occurs between an individual pair of oscillators, their phase difference takes a constant value. ML-DFA detects the resulting loss of scaling by indicating that the DFA fluctuation plot is no longer linear.

FIGURE 8

Figure 8. Results of the phase synchrony analysis method when applied to the Kuramoto model. There are 200 oscillators, with a mean natural frequency of 22 Hz, and a standard deviation of natural frequencies of 15. The theoretical critical coupling K_noise when noise is added is marked with a blue asterisk. The average DFA exponent, order parameter r, its difference Δ(Kr) and the proportion of valid DFA fluctuation plots from the full set of 199,000 pairs are shown. Validity bars that are referred to in the text are highlighted in magenta.

After the peak DFA at K ≈ 22, further increase in K eventually causes full synchronization between all individual oscillator pairs. Across the whole system, fewer than 10% of oscillator pairs yield a valid DFA after the critical coupling is exceeded. When all oscillator pairs are synchronized with each other, the order parameter of the system approaches its maximum level of 1 but the DFA fluctuation measure of rate of change of phase difference is no longer valid.

Analysis of the Kuramoto model with noise suggests that LRTCs in the rate of change of phase difference between oscillator pairs occur when the system is in a state of maximal flux just prior to the onset of full synchronization.

3.3.1. Individual oscillators pairs

Further insights into the rate of change of phase difference fluctuation behavior can be obtained from DFA of individual oscillator pairs. Analysis of a set of 5 oscillator pairs is shown in Figure 9. The top panel shows the change in DFA exponent with coupling K for a pair whose initial frequencies are very close (0.001 Hz apart). The bottom panel shows the changes in DFA exponent for an oscillator pair with initial frequencies that differ by ≈7.0 Hz. The middle panels show oscillator pairs with varying amounts of initial frequency difference (increasing top to bottom). Non valid DFA exponents are not plotted in the left hand panel but the right hand panels indicate for each given pair linear DFA validity “yes” or “no” for a given value of K. At low coupling K, the oscillators do not interact with each other and each evolves at its own natural frequency. The order in the system is low and the DFA exponent ≈0.5 reflects the additive noise which dominates the fluctuations in the rate of change of phase difference. A DFA value of ≈0.5 is also evident in the average DFA (Figure 8). There is almost 100% validity across all pairs because white noise time series are scale-free and therefore the DFA fluctuation plot obtained from analysing them is expected to be linear (Figure 8).

FIGURE 9

Figure 9. Representative relationship of DFA exponents to the coupling parameter K for selected oscillator pairs in the Kuramoto system. (A,C,E,G,I) Show the value of valid DFA exponents, while (B,D,F,H,J) indicate whether the exponent is rejected as invalid by the ML-DFA technique (N) or not (Y). The oscillator numbers and the differences between their natural frequencies are recorded in the legend of (A,C,E,G,I). The first number is the difference in natural frequency (in Hz), and the subsequent pair of numbers identifies which oscillators are being analyzed.

As the coupling parameter K is increased, the DFA exponents of each of the oscillator pairs rise until a peak is reached. The value of K at which a maximal valid exponent is retrieved for these peaks is related to the difference in natural frequencies of the two oscillators as well as their interactions with the noise and the mean field. Oscillator pairs which start further apart in frequency terms develop full synchonization later than those whose initial frequencies are close together. As K increases the DFA exponent of the rate of change of phase difference increases. The pairs with the strongest LRTCs on the basis of the highest DFA exponent value prior to onset of full synchronization are those with the greatest inital frequency difference. Increasing temporal order of the rate of change of phase difference prior to full synchonization of these pairs may indicate a state of pre-synchronization in these pairs.

3.4. The Cabral Model

For the Cabral model we present results regarding both the global behavior of the system through average DFA exponents across all possible pairs of oscillators (Figure 10) and the behavior of the system at cluster level through average DFA exponents of intra-cluster pairs of oscillators (Figure 11).

FIGURE 10

Figure 10. The average DFA exponents of phase synchrony as a function of the coupling parameter, K, in the extended Kuramoto model (Cabral et al., 2011). The model includes 66 oscillators at normally distributed natural frequencies with mean 60 Hz and standard deviation σ_i = 1.25. The connectivity and time delay matrices are set from empirical values. The average of the valid DFA exponents is shown in magenta and the proportion of valid exponents, as calculated by ML-DFA, are indicated by bars. The Kuramoto model order parameter r is in blue, and the quantity Δ(Kr) is in cyan. The peak Δ(Kr) has been used as an indicator of the effective critical coupling. A horizontal line at DFA exponent 0.5 is plotted to guide the eye. The proportion of valid DFA bars for K = 0, K = 1, and K = 6 have been shaded in magenta.

FIGURE 11

Figure 11. Average DFA exponent for intra-cluster pairwise phase differences with increasing coupling parameter K. Where no DFA value appears for a particular cluster, this indicates that there are no valid DFA exponents for the pairwise phase difference within that cluster. The final cluster, which is labeled individual oscillators, consists of a set of nodes that do not fit into any of the clusters as determined by the connectivity and distance matrices but are grouped together to demonstrate their relationship with each other.

3.4.1. Global behavior

The model introduced by Cabral et al. (2011) is affected by rich interplay between the connectivity and distance matrices as well as the noise and natural frequency elements of the system. The average valid DFA exponents for all oscillator pairings (n = 2145) are shown in Figure 10 as the coupling in the system is increased. These average exponent values indicate the presence of LRTCs in the rate of change of phase difference. The peak values of mean DFA exponent correspond to peaks in the change in order paramenter (Δ(Kr)) derived for the classical Kuramoto model and the Kuramoto model with noise, see Kitzbichler et al. (2009) and Figure 8. Such peaks occur when the system undergoes the greatest change in synchronization. The peak in Δ(Kr) corresponds closely to the coupling value that shows maximum mean DFA exponent (K = 5 and 6, respectively—see Figure 10).

The number of pairings that yield valid DFA exponents in the rate of change of their phase difference is equal to 100% when there is no coupling in the system (magenta shaded bar at K = 0), but it falls as coupling is introduced (magenta shaded bar at K = 1). At the coupling value of the DFA peak, K = 6, validity is at 20%, which is higher than the neighboring coupling values (magenta shaded bar at K = 6).

3.4.2. Cluster behavior

At coupling value K = 6, the value at which the global behavior shows peak DFA value, the intra-cluster results indicate that only cluster 4, consisting of oscillators 27–40, shows valid non-trivial DFA exponents. These exponents are consistent with the presence of LRTCs. This suggests that cluster 4 acts as an organizing force in the system when the system is in its greatest state of flux, as demonstrated by a large increase in the order parameter. This cluster corresponds to the most connected brain regions listed in Table 1 and Table A1 in the Appendix.

The connectivity and distance matrices for the Cabral model are shown in Figure 12. The linear coupling between oscillators for two values of K is shown in Figure 13. The central cluster of oscillators with high levels of synchronization is evident from the two correlation matrices. At K = 6 (Figure 13A), i.e., the value at which LRTCs are detected in the rate of change of phase difference, the central oscillator cluster shows evidence of synchronization but with Pearson correlation values of <1.0. As K increases to 18, the value identified by Cabral et al. (2011) as best approximating human brain resting state BOLD fluctuations, it can be seen from Figures 10, 11 that the proportion of oscillator pairs with valid DFA fluctuation plots is low (approximately 5%). Those oscillator pairs that remain and show persistently valid DFA fluctuation plots are predominantly individual oscillators with low average weight per node (0.03) and low average degree distribution (8.59). Their associated DFA exponent is on average 0.5 (see Figure 11). At K = 18, the Cabral model shows strong cluster synchronization. In particular, the central cluster 4 (oscillators 27–40) which contains homologous elements connected across the corpus callosum shows Pearson correlation values close to 1.0 indicative of full synchrony (Figure 13B). Therefore, the results we obtained from the Kuramoto model with noise and those derived from the Cabral model are similar. Both show valid DFA fluctuation plots with LRTCs of the rate of change of phase difference at a coupling value where Δ(Kr) is increasing and loss of validity as full synchronization takes over. As discussed earlier, “criticality” is not defined for the Cabral model but with increasing K there is clearly a change in the system's order which is detected through our method.

FIGURE 12

Figure 12. Connectivity and distance matrices for the Cabral model. Each oscillator number represents a brain region, which is defined in Table A1 in the Appendix. An empty (white) element means that the two regions are not connected. Regions are not connected to themselves so that the diagonals are white. (A) Shows the pairwise connection matrix C between the 66 oscillators. (B) Shows the matrix of pairwise distances L between the brain regions that are represented by the 66 oscillators. Matrix L is symmetric, however, matrix C is not because the connection weights are normalized by row. The values associated with the colors of the plots are defined by the color bars. Red colors in (A) represent higher weights. Red colors in (B) represent longer distance connections.

FIGURE 13

Figure 13. Correlation matrices for all pairs of time series generated by the Cabral model for two coupling values K. (A) K = 6 and (B) K = 18, which corresponds to the oscillator correlation matrix in Cabral et al. (2011). The plots show the value of the Pearson correlation coefficient between all pairwise combinations of the 66 oscillators used in the model.

Figure 14 shows the DFA exponents of the rate of change of phase difference between individual pairs of oscillators in the form of a symmetric lattice of size 66 × 66, where each element in the lattice represents a brain region as detailed in Table A1 of the Appendix. Figure 14A of this figure shows the importance of the central cluster in generating LRTCs of phase synchronization. Importantly it shows this cluster's influence over many of the other oscillators in the Cabral model. Cluster group 4 has the greatest sum of weights per oscillator and the greatest number of connections per oscillator (see Table 1). The correlation between the number of connections of a given oscillator and the average DFA exponent of its rate of change of phase difference with all other oscillators is 0.359, suggesting a relationship between oscillators with large connectivity and those with large DFA exponents in their pairwise phase difference.

FIGURE 14

Figure 14. DFA exponent of the rate of change of phase difference between all pairs of oscillators in the Cabral model at coupling K = 6 in three scenarios. (A) For the empirically observed connectivity matrix of the Cabral model. (B) For a connectivity matrix representing disconnected hemispheres. (C) For random connectivity. Empty (white) elements denote pairs for which no valid DFA exponent was found.

3.4.3. Comparison of the three connectivity structures

In the Cabral model, the Δ(Kr) measure has its peak at coupling value K = 6. Here, we compare the effects of the three connectivity matrices introduced in Section 2.7.5 on the DFA exponents of the pairwise phase difference between oscillators at this coupling value in Figure 14.

The empirical connectivity matrix showed large DFA exponents indicating the presence of LRTCs at this coupling value for a small number of hub oscillators belonging to cluster 4 (see above). These oscillators have a high number of connections and large weights associated with these connections (see Table 1). When the two hemispheres are disconnected, we see no LRTCs in the DFA exponents of the phase difference at this coupling value. When the distance matrix is preserved, but the connectivity and associated weights are assigned at random, LRTCs are still present in the DFA exponent of the phase differences between oscillators, but a lower value of DFA exponent is obtained. There is no apparent cluster formation when connectivity is random.

3.4.4. Neurophysiological data

Figure 15 illustrates the application of our phase synchrony analysis technique to the human neurophysiological data described in Section 2.8. In this example, a valid DFA exponent of ≈0.6 was obtained for the rate of change of phase difference between the simultaneously recorded EEG and EMG data during a steady muscle contraction, indicative of the presence of LRTCs. Analysis of amplifier noise and artificially generated noise time series using processing steps identical to those for the EEG and EMG data (signal processing pathway shown in Figure 1) resulted in a valid DFA fluctuation plot but with exponent of 0.48 consistent with uncorrelated noise.

FIGURE 15

Figure 15. Illustration of the method with simultaneous EEG/EMG data. (A,B) One second of simultaneously recorded EEG and EMG, respectively. (C,D) The signals after bandpass filtering in the β range 16–24 Hz. (E) Rate of change of the phase difference between the two bandpass filtered signals for 100 s. (F) DFA fluctuation plot for the rate of change of phase difference time series in (E). The plot was determined to be valid by ML-DFA with a DFA exponent of ≈0.60, indicating the presence of LRTCs.

4. Discussion

The aim of this paper is to introduce a new methodology for eliciting a marker of criticality in neuronal synchronization. This methodology relies on the rate of change of the phase difference between two signals as a (time-varying) measure of phase synchronization. The presence of LRTCs in this quantity is proposed as marker of criticality and is assessed using DFA in combination with the recently proposed ML-DFA, a heuristic technique for validating the output of DFA. With these methods, we can first determine the presence or absence of power law scaling using ML-DFA and secondly the presence or absence of LRTCs in the phase synchronization of two time series based on the value of the DFA exponent. If the method returns an exponent of ≈0.5, this indicates a phase relationship similar to white Gaussian noise, however, if the DFA exponent is greater than 0.5, this indicates the presence of LRTCs. Importantly, we can attribute significance to the loss of power law scaling within the fluctuation plot and draw conclusions based on an exponent value only when the exponent has been recovered from plots that are judged to be valid by ML-DFA.

4.1. Surrogate Data

It was found that the phase synchrony analysis method recovers a known DFA exponent value in the rate of change of phase difference between two signals of surrogate data with a high degree of accuracy (r = 0.998). When the structure of phase synchronization was perturbed with an additive noise source, it was found that a percentage difference between the true and recovered DFA exponent of above approximately 5% noise caused DFA exponents to be judged as invalid by ML-DFA. When the surrogate data was characterized by a DFA exponent close to 1, the recovery of this exponent using DFA was more resistant to noise when compared to surrogate data with a lower DFA exponent of 0.6 (Figure 6). In these simulations we used additive noise which was included at the amplitude stage of the surrogate time series prior to extraction of the phase using the Hilbert transform.

4.2. The Ising Model

We had initially expected to see LRTCs in the Ising model only in the vicinity of the critical parameter, and a DFA exponent of 0.5 when the energy in the system was large (disordered phase). However, in applying our method to the Ising model, both of these hypotheses were not fully realized. It was found that when the temperature was increased to a very high level of T = 10⁵, the DFA exponent of the rate of change of phase difference did not fall to 0.5, but remained at ≈0.57. This did not change when the temperature was set to an even higher value of T = 10¹². This was not a finite size effect of the system, as the result held when larger lattice sizes (up to 1000 × 1000) were used (results not shown). We noted that when pure phase was analyzed, i.e., an uncoupled system of Kuramoto oscillators, DFA exponents of 0.5 were obtained as expected, and therefore, we cannot exclude the possibility that the Hilbert transform induced artifacts may inject some order into the resulting phase time series. However, within the Ising system, the expectation of a DFA exponent of 0.5 at high T is based only on our intuition concerning the operation of the system. As all elements in the Ising lattice interact with their neighbors it is possible that some temporal correlation in the rate of change of phase difference may persist regardless of temperature value, and this may be the cause of a DFA exponent above 0.5.

Importantly, we found that the DFA exponent was indicative of LRTCs at critical temperature but was maximal at T = 2.55, just in excess of the critical temperature. As can be seen in Figure 7, the consistent change in the DFA value and the change in power law scaling behavior indicates that the phase synchrony analysis method is capturing an important behavior of the system close to its critical regime. However, it is important to realize that unless an experimental neuroscientific paradigm can be discovered that produces similar consistent changes in this measure, neurophysiological data will have to be intepreted with caution, i.e., we may be able to state that for a given pair of neural oscillation time series there exists power law scaling with a DFA exponent indicative of LRTCs in the rate of change of their phase difference but we may not know whether for this neural state there may exist other higher (or lower) exponent values. In other words, the technique may provide evidence that the system is ordered in ways that are similar to systems nearing their critical regime but whether the technique will pinpoint the most critical regime in a neural system is open to question. We will consider this further in our discussion of the results of analysing a Kuramoto system with noise.

Interestingly, the evolution of the DFA exponent with the temperature parameter shares a key characteristic with that of a recently published measure of information flow in the same model (Barnett et al., 2013), specifically, an asymmetry around the critical point, with a sharp rise in the metric as temperature is increased toward the critical T = Tc and a gradual descent as the temperature rises significantly. It would be of interest to further assess the extent to which the proposed method captures information flow in the system, e.g., through a comparison of both methods when applied to the Kuramoto model.

4.3. The Kuramoto Model

In the Kuramoto model, the critical transition is characterized in terms of a global order parameter which reflects the overall organization of the system. However, through our phase synchrony analysis method we are able to make observations at a pair-wise level of Kuramoto oscillators always bearing in mind that even at the pair-wise level the result is influenced by the oscillators' interactions with all other oscillators in the model. As individual Kuramoto oscillator pairs become fully synchronized, their rate of change of phase difference no longer contains moment-to-moment fluctuations and thus power law scaling in the DFA measure is lost. This is an important consideration because it emphasizes the difference between our method and more standard measures of neural synchrony. Methods for detecting neural synchrony rely on phase consistency to allow averaging out of fluctuations so that a measure of coupling (e.g., coherence and phase coherence) is obtained. In contrast, the method introduced in this paper is dependent on the fluctuations of the two phase signals and their interaction. Therefore, our method detects “order” across time in the rate of change of phase difference rather than phase consistency between two processes.

The phenomenon of loss of fluctuations at the onset of full synchronization is well illustrated both for the global Kuramoto model and for individual oscillator pairs extracted from the Kuramoto model. In the global analysis the peak in the DFA exponent occurs close to the observed peak of Δ(Kr) and at values of K just below theoretical critical coupling value. At these values of K, a power law scaling exists for the rate of change of phase difference, and the DFA exponent of oscillator pairs with different initial frequencies indicates the presence of LRTCs. At the onset of full synchronization the number of oscillator pairs for which DFA is valid drops yet those whose phase differences still possess fluctuations continue to show LRTCs. Once the critical regime has been fully crossed and the order parameter r approaches 1, the DFA of the rate of change of phase difference is no longer valid for any oscillator pair.

The LRTC behavior is also clearly explained as the coupling value K decreases toward zero. As can be seen in Figure 8, the DFA exponent of the pairwise rates of change of phase difference decreases toward 0.5 and yet scaling remains valid. These changes in DFA exponent are evident both on the global level in the average DFA and for individual oscillator pairs. At K = 0 the phases are independent from one another yet contain noise; thus the rate of change of phase difference time series contains innovations that are random across time with a DFA which is valid and returns the expected exponent of 0.5.

4.4. Order within the Ising and Kuramoto Models

In these models, temperature T (Ising) and coupling K (Kuramoto) play a similar role in controlling the order within the two systems, and the DFA validity and exponent results obtained from analysis of rate of change of phase difference in both of these models mirror each other. In the Kuramoto model, there is a transition from an uncoupled to a synchronized state with increasing K. Similarly in the Ising model, there is a transition from a very ordered to a disordered system with increasing T. In the human brain, we are not able to characterize the system by incrementally tuning a parameter and observing the result, and we are only privy to snapshots of the working system. However, we can begin to understand the behavior of the brain within this range of behaviors by comparing the DFA of the rate of change of phase difference of pairs of neurophysiological signals to the outcomes of these models of criticality.

4.5. The Cabral Model

We found that LRTCs exist in the rate of change of phase difference between oscillator pairs at parameter values close to those at which the change in order, Δ(Kr), increases sharply. Extrapolating from the Kuramoto model with noise, we suggest that there are important changes in the order of the phase synchronization of interacting oscillators in the Cabral model that involve the presence of LRTCs when the order in that system is at or close to a point of maximal change.

It is important to note that the value of r in the Cabral model does not reach a level of 1 in the range of coupling values 0–20. It approaches a level of ≈0.4 as K approaches 20 with maximal rate of change at K ≈ 6. Further analysis of the Cabral model indicates that r will gradually reach a value closer to 1 as K increases above a value of 60, as seen in Figure 4 of Cabral et al. (2011). Cabral focussed her attention on K = 18 at which point the model, when fed through the Balloon-Windkessel hemodynamic model, produced an output that closely matched the spatio-temporal correlations seen in the BOLD signals of the resting state fMRI. We find that at this value, there are no LRTCs detectable in the rate of change of phase difference measure.

4.6. The Role of Connectivity in the Cabral Model

Although most of results were obtained at K = 6, selected because it is the peak of Δ(Kr), it is important to note that LRTCs exist for a broader range of coupling values K. This finding agrees with a recent study by 2013 (Moretti and Muñoz (2013) in which the authors demonstrated that a network with complex connectivity, such as that of the Cabral model and, indeed, that of the brain, causes the critical point to becomes a broader critical “region.”

Our examination of oscillator pairs belonging to a single cluster, as defined in Cabral et al. (2011), indicates that the emergence of LRTCs is determined primarily by oscillators belonging to cluster 4 which has a large number of connections and a large sum of connection weights. This cluster is located centrally, and it contains four brain regions of particular importance to the resting state network (Fransson and Marrelec, 2008; van den Heuvel and Sporns, 2011). These are oscillators 33 and 34, which correspond to the left and right posterior cingulate cortices, and oscillators 32 and 35 which represent the left and right precuneus. These central brain regions are known to be important with a higher metabolic activity than other regions during the resting state.

Importantly, we find that LRTC behavior of this cluster, and its relationship to the other clusters in the network, is dependent on trans-callosal left-right connectivity. Indeed, disruption of the left-right trans-callosal connections resulted in a loss of LRTCs in the rate of change of phase difference between time series extracted from the central cluster 4 and the other oscillators in the Cabral network. Intuitively, those oscillators that are connected to many other oscillators in the network will also influence the phases of a large number of other oscillators. When these oscillators try to synchronize, we suggest that those that are well connected will be subjected to conflicting phase inputs from their neighbors and thus increased variation in their phase fluctuations, yielding a larger DFA exponent. These variations in fluctuation will in turn feed into the neighboring oscillators and cause them to also have large variations in fluctuation as they attempt to synchronize with their well-connected neighbor. On the other hand, an oscillator that is poorly connected or connected to just one other oscillator may have a more straightforward task of synchronizing with just this (albeit changing) oscillator speed.

The LRTCs in the rate of change of phase difference were also disrupted by randomization of connectivity, albeit less severely than when the trans-callosal connections were severed. When a random connectivity is assigned, no clusters exist and DFA exponents are significantly reduced.

The results obtained from the phase synchrony analysis method here may pave the way for potential future use of the Cabral model in investigating specific pathological modifications of connectivity and their effects on the time-varying synchronization patterns between different brain regions. The method has the potential to be used to trace some types of pathological synchronization such as may arise in epileptic or Parkinsonian conditions to any roots that they may have either in the connectivity, clustering or noise input elements of the Cabral model and therefore potentially also of the nervous system.

4.7. Neurophysiological Data

In order to show proof of principle, we have presented an example of our method's application to neurophysiological data, in this case EEG and EMG simultaneously recorded during voluntary muscle contraction. It was through this experimental paradigm that corticomuscular coherence (CMC) in the 16–32 Hz (β) frequency range was first discovered by Conway et al. (1995), Halliday et al. (1998) and shown to be the β frequency common drive to human motoneurons first described by Farmer et al. (1993). These preliminary results indicate power law scaling in the DFA plot with a DFA exponent of ≈0.6.

It has been recognized through application of time-varying coherence measures that CMC coherence fluctuates even when a subject attempts to maintain the same motor output (Muthukumaraswamy, 2011). As discussed earlier, the techniques introduced here allow us to focus on the fluctuations within the phase coupling rather than on the averaged measure of coupling. These preliminary results indicate that the fluctuations in the rate of change of phase difference between simultaneously recorded EEG and EMG show power law scaling and LRTCs within the β frequency range. We suggest that the analysis of instantaneous phase diffence of neurophysiological data using the methods described in this paper will allow researchers to investigate the coupling between signals in a way that will allow a new appreciation of the relationship between neural synchrony and other oscillator systems approaching their critical regime.

4.8. LRTCs in Rate of Change of Phase Difference and the Brain

LRTCs have been associated with model dynamical systems that show efficiency in learning, memory formation, rapid information transfer, and network organization. The broad dynamical range of which LRTCs are a marker acts to support these functions (Linkenkaer-Hansen et al., 2001, 2004; Stam and de Bruin, 2004; Sornette, 2006; Shew et al., 2009; Chialvo, 2010; Werner, 2010; Beggs and Timme, 2012; Meisel et al., 2012). It has been argued by a number of researchers that these properties if present would be of major benefit to the functions that human brain dynamics needs to support and there is now a literature that connects the theory of critical systems with properties of human brain dynamics (Linkenkaer-Hansen et al., 2001; Beggs and Plenz, 2003; Kitzbichler et al., 2009; Shew et al., 2009; Chialvo, 2010).

In this paper, we focus on LRTCs, and because of the importance in neuroscience of brain oscillations and the concept of communication through coherence, we make the link between LRTCs and phase synchrony. We note that in the model systems that we have explored the highest valid DFA exponents were recovered when the systems were close to their critical point but in a slightly more disordered state than at exact criticality. We explained this on the basis of full synchronization within our model systems being a point at which the rate of change of phase difference is lost (observed in Ising at T < T_c and in Kuramoto for increasing K).

In neurophysiological systems, it is important to appreciate that full synchronization of neural oscillators is a pathological state (e.g., observed in the EEG and MEG of epileptic seizures and in EMGs showing pathological tremor). The healthy resting brain state therefore is characterized by weak and variable neural synchrony which would be expected to show fluctuations (temporal innovations) in a measure of the change in phase synchrony, i.e., the rate of change of phase difference. From the perspective of brain dynamics (and muscle activation dynamics) the most important constraints are to avoid pathological synchronization whilst at the same time maintaining the potential for useful synchronization. We suggest therefore that in the healthy state the instantaneous phase difference between neural oscillators will show power law fluctuation plots with a DFA exponent that is either 0.5 or that will show LRTCs. If LRTCs are found in the resting state then they may represent an optimum state of readiness to which the system can readily return if increased synchronization occurs as a result of sensory stimulation, motor task, or cognitive action. Such temporary changes in synchronization may occur in order to support communication through coherence. The resting state, however, is characterized by fluctuations of phase synchrony that have LRTCs and represent the behavior of weakly coupled oscillators whose synchrony can be modulated. The hypothesis that the LRTCs of rate of change of phase difference of brain oscillations may be altered through task is an experimentally tractable question.

To conclude the evidence for the brain as a critical system continues to accrue. There is an important need to link the criticality paradigm with the paradigm that attaches functional significance to neural synchrony. The methodology presented in this paper takes us some way toward this synthesis.

Funding

Engineering and Physical Sciences Research Council (EPSRC).

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.

Acknowledgment

Maria Botcharova thanks the Centre for Mathematics and Physics in the Life Sciences and Experimental Biology (CoMPLEX), University College London for their funding and continuing support. Simon F. Farmer was supported by University College London Hospitals Biomedical Research Centre (BRC).

References

Acebrón, J. A., Pérez Vicente, C. J., Ritort, F., and Spigler, R. (2005). The Kuramoto model: a simple paradigm for synchronization phenomena. Rev. Mod. Phys. 77, 137–185. doi: 10.1103/RevModPhys.77.137

CrossRef Full Text | Google Scholar

Akaike, H. (1974). A new look at the statistical model identification. IEEE Trans. Automat. Control, 19, 716–723. doi: 10.1109/TAC.1974.1100705

CrossRef Full Text | Google Scholar

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 | Pubmed Full Text | CrossRef Full Text | Google Scholar

Bak, P. (1996). How Nature Works: The Science of Self-organized Criticality, 1st Edn. Copernicus: Springer. doi: 10.1007/978-1-4757-5426-1

CrossRef Full Text

Baker, S. N., Kilner, J. M., Pinches, E. M., and Lemon, R. N. (1999). The role of synchrony and oscillations in the motor output. Exp. Brain Res. 128, 109–117. doi: 10.1007/s002210050825

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Baker, S. N., Spinks, R., Jackson, A., and Lemon, R. N. (2001). Synchronization in monkey motor cortex during a precision grip task. I. Task-dependent modulation in single-unit synchrony. J. Neurophysiol. 85, 869–885.

Pubmed Abstract | Pubmed Full Text | Google Scholar

Bardet, J.-M., and Kammoun, I. (2008). Asymptotic properties of the detrended fluctuation analysis of long-range-dependent processes. IEEE Trans. Inform. Theory 54, 2041–2052. doi: 10.1109/TIT.2008.920328

CrossRef Full Text | Google Scholar

Barnett, L., Lizier, J. T., Harré, M., Seth, A. K., and Bossomaier, T. (2013). Information flow in a kinetic ising model peaks in the disordered phase. Phys. Rev. Lett. 111, 177203. doi: 10.1103/PhysRevLett.111.177203

Pubmed Abstract | Pubmed Full Text | 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 | Pubmed Full Text | CrossRef Full Text | Google Scholar

Beggs, J. M., and Timme, N. (2012). Being critical of criticality in the brain. Front. Fract. Physiol. 3:163. doi: 10.3389/fphys.2012.00163

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Beran, J. (1994). Statistics for Long-memory Processes. Monographs on statistics and applied probability, Vol. 61. New York, NY: Chapman & Hall.

Google Scholar

Berthouze, L., James, L. M., and Farmer, S. F. (2010). Human EEG shows long-range temporal correlations of oscillation amplitude in theta, alpha and beta bands across a wide age range. Clin. Neurophysiol. 121, 1187–1197. doi: 10.1016/j.clinph.2010.02.163

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Bonilla, L. L., Neu, J. C., and Spigler, R. (1992). Nonlinear stability of incoherence and collective synchronization in a population of coupled oscillators. J. Stat. Phys. 67, 313–330. doi: 10.1007/BF01049037

CrossRef Full Text | Google Scholar

Botcharova, M., Farmer, S. F., and Berthouze, L. (2012). Power-law distribution of phase-locking intervals does not imply critical interaction. Phys. Rev. E 86, 051920. doi: 10.1103/PhysRevE.86.051920

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Botcharova, M., Farmer, S. F., and Berthouze, L. (2013). A maximum likelihood based technique for validating detrended fluctuation analysis (ML-DFA). arXiv:1306.5075

Google Scholar

Breakspear, M., Heitmann, S., and Daffertshofer, A. (2010). Generative models of cortical oscillations: neurobiological implications of the Kuramoto model. Front. Hum. Neurosci. 4:190. doi: 10.3389/fnhum.2010.00190

Pubmed Abstract | Pubmed Full Text | 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 | Pubmed Full Text | CrossRef Full Text | Google Scholar

Bryce, R. M., and Sprague, K. B. (2012). Revisiting detrended fluctuation analysis. Sci. Rep. 2, 315. doi: 10.1038/srep00315

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Buzsaki, G. (2006). Rhythms of The Brain, 1 Edn. London: Oxford University Press. doi: 10.1093/acprof:oso/9780195301069.001.0001

CrossRef Full Text

Cabral, J., Hugues, E., Sporns, O., and Deco, G. (2011). Role of local network oscillations in resting-state functional connectivity. Neuroimage 57, 130–139. doi: 10.1016/j.neuroimage.2011.04.010

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Cabral, J., Kringelbach, M. L., and Deco, G. (2012). Functional graph alterations in schizophrenia: a result from a global anatomic decoupling? Pharmacopsychiatry 45, 57–64. doi: 10.1055/s-0032-1309001

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Carreras, B. A., van Milligen, B., Pedrosa, M. A., Balbín, R., Hidalgo, C., Newman, D. E., et al. (1998). Long-range time correlations in plasma edge turbulence. Phys. Rev. Lett. 80, 4438–4441. doi: 10.1103/PhysRevLett.80.4438

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

Chopra, N., and Spong, M. W. (2005). “On synchronization of Kuramoto oscillators,” in 44th IEEE Conference on Decision and Control, 2005 and 2005 European Control Conference. CDC-ECC '05, Vol. 44 (Seville), 3916–3922.

Google Scholar

Clauset, A., Moore, C., and Newman, M. E. J. (2006). Structural Inference of Hierarchies in Networks. arXiv:physics/0610051

Google Scholar

Clegg, R. G. (2006). “A practical guide to measuring the hurst parameter,” in 21st UK Performance Engineering Workshop, School of Computing Science Technical Report Series, CSTR-916, University of Newcastle (Newcastle, NSW), 43–55.

Google Scholar

Conway, B., Halliday, D., Farmer, S., Shahani, U., Maas, P., Weir, A., et al. (1995). Synchronization between motor cortex and spinal motoneuronal pool during the performance of a maintained motor task in man. J. Physiol. 489(Pt 3), 917–924.

Pubmed Abstract | Pubmed Full Text | Google Scholar

Daffertshofer, A., and van Wijk, B. C. M. (2011). On the influence of amplitude on the connectivity between phases. Front. Neuroinform. 5:6. doi: 10.3389/fninf.2011.00006

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Daido, H. (1989). Intrinsic fluctuation and its critical scaling in a class of populations of oscillators with distributed frequencies. Prog. Theor. Phys. 81, 727–731. doi: 10.1143/PTP.81.727

CrossRef Full Text | Google Scholar

Doesburg, S. M., Green, J. J., McDonald, J. J., and Ward, L. M. (2009). Rhythms of consciousness: binocular rivalry reveals large-scale oscillatory network dynamics mediating visual perception. PLoS ONE 4:e6142. doi: 10.1371/journal.pone.0006142

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Dörfler, F., and Bullo, F. (2011). On the critical coupling for kuramoto oscillators. SIAM J. Appl. Dyn. Syst. 10, 1070–1099. doi: 10.1137/10081530X

CrossRef Full Text | Google Scholar

Farmer, S. F. (1998). Rhythmicity, synchronization and binding in human and primate motor systems. J. Physiol. 509, 3–14. doi: 10.1111/j.1469-7793.1998.003bo.x

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Farmer, S. F., Bremner, F. D., Halliday, D. M., Rosenberg, J. R., and Stephens, J. A. (1993). The frequency content of common synaptic inputs to motoneurones studied during voluntary isometric contraction in man. J. Physiol. 470, 127–155.

Pubmed Abstract | Pubmed Full Text | Google Scholar

Fransson, P., and Marrelec, G. (2008). The precuneus/posterior cingulate cortex plays a pivotal role in the default mode network: evidence from a partial correlation network analysis. Neuroimage 42, 1178–1184. doi: 10.1016/j.neuroimage.2008.05.059

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Freeman, W., and Rogers, L. (2002). Fine temporal resolution of analytic phase reveals episodic synchronization by state transitions in gamma EEGs. J. Neurophysiol. 87, 937–945. doi: 10.1152/jn.00254.2001

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Freeman, W. J. (2004). Origin, structure, and role of background EEG activity. Part 1. Analytic amplitude. Clin. Neurophysiol. 115, 2077–2088. doi: 10.1016/j.clinph.2004.02.029

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | 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 | Pubmed Full Text | CrossRef Full Text | Google Scholar

Gionis, A., Mannila, H., Mielikinen, T., and Tsaparas, P. (2007). Assessing data mining results via swap randomization. ACM Trans. Knowl. Discov. Data 1, 14–23. doi: 10.1145/1297332.1297338

CrossRef Full Text | Google Scholar

Gong, P., Nikolaev, A. R., and van Leeuwen, C. (2007). Intermittent dynamics underlying the intrinsic fluctuations of the collective synchronization patterns in electrocortical activity. Phys. Rev. E 76, 011904. doi: 10.1103/PhysRevE.76.011904

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Granger, C. W. J., and Joyeux, R. (1980). An introduction to long-memory time series models and fractional differencing. J. Time Ser. Anal. 1, 15–29. doi: 10.1111/j.1467-9892.1980.tb00297.x

CrossRef Full Text | 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. doi: 10.1016/S0304-3940(97)00964-6

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Hanhijärvi, S., Garriga, G. C., and Puolamäki, K. (2009). “Randomization techniques for graphs,” in Proceedings of the 9th SIAM International Conference on Data Mining (SDM '09) (Sparks, NV), 780–791.

Hardstone, R., Poil, S.-S., Schiavone, G., Jansen, R., Nikulin, V. V., Mansvelder, H. D., et al. (2012). Detrended fluctuation analysis: a scale-free view on neuronal oscillations. Front. Physiol. 3:450. doi: 10.3389/fphys.2012.00450

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Hausdorff, J. M., Peng, C. K., Ladin, Z., Wei, J. Y., and Goldberger, A. L. (1995). Is walking a random walk? Evidence for long-range correlations in stride interval of human gait. J. Appl. Physiol. 78, 349–358.

Pubmed Abstract | Pubmed Full Text | Google Scholar

Honey, C. J., Sporns, O., Cammoun, L., Gigandet, X., Thiran, J. P., Meuli, R., et al. (2009). Predicting human resting-state functional connectivity from structural connectivity. Proc. Natl. Acad. Sci. U.S.A. 106, 2035–2040. doi: 10.1073/pnas.0811168106

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Hosking, J. R. M. (1981). Fractional differencing. Biometrika 68, 165–176. doi: 10.1093/biomet/68.1.165

CrossRef Full Text | 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 | Pubmed Full Text | CrossRef Full Text | Google Scholar

Hurvich, C. M., and Tsai, C.-L. (1989). Regression and time series model selection in small samples. Biometrika 76, 297–307. doi: 10.1093/biomet/76.2.297

CrossRef Full Text | Google Scholar

Ising, E. (1925). Beitrag zur theorie des ferromagnetismus. Zs. Physik 31, 253–258. doi: 10.1007/BF02980577

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 | Pubmed Full Text | CrossRef Full Text | Google Scholar

Karmeshu., and Krishnamachari, A. (2004). Sequence variability and long-range dependence in DNA: an information theoretic perspective. Lect. Notes Comput. Sci. 3316, 1354–1361. doi: 10.1007/978-3-540-30499-9_210

CrossRef Full Text | Google Scholar

Kinouchi, O., and Copelli, M. (2006). Optimal dynamical range of excitable networks at criticality. Nat. Phys. 2, 348–351. doi: 10.1038/nphys289

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 | Pubmed Full Text | CrossRef Full Text | Google Scholar

Kuramoto, Y. (1975). “Self-entrainment of a population of coupled nonlinear oscillators,” in International Symposium on Mathematical Problems in Theoretical Physics, Lecture Notes in Physics, Vol. 39, ed H. Araki (New York, NY: Springer), 420–422. doi: 10.1007/BFb0013365

CrossRef Full Text | Google Scholar

Kuramoto, Y. (1984). Chemical Oscillations, Waves, and Turbulence. New York, NY: Springer–Verlag. doi: 10.1007/978-3-642-69689-3

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 | Pubmed Full Text | 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 | Pubmed Full Text | CrossRef Full Text | Google Scholar

Mackay, D. J. C. (2003). Information Theory, Inference and Learning Algorithms, 1st Edn. Cambridge: Cambridge University Press.

Google Scholar

Meisel, C., Storch, A., Hallmeyer-Elgner, S., Bullmore, E., and Gross, T. (2012). Failure of adaptive self-organized criticality during epileptic seizure attacks. PLoS Comput. Biol. 8:e1002312. doi: 10.1371/journal.pcbi.1002312

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Mertens, D. C. (2011). Population-Specific Predictions for The Finite Kuramoto Model and Collective Synchronization in A System With Resonant Coupling. Ph.D. thesis, University of Illinois at Urbana-Champaign.

Google Scholar

Metropolis, N., Rosenbluth, A., Rosenbluth, M., Teller, A., and Teller, E. (1953). Equation of state calculations by fast computing machines. J. Chem. Phys. 21, 1087. doi: 10.1063/1.1699114

CrossRef Full Text | Google Scholar

Miritello, G., Pluchino, A., and Rapisarda, A. (2009). Phase transitions and chaos in long-range models of coupled oscillators. Europhys. Lett. 85, 10007. doi: 10.1209/0295-5075/85/10007

CrossRef Full Text | Google Scholar

Moretti, P., and Muñoz, M. A. (2013). Griffiths phases and the stretching of criticality in brain networks. Nat. Commun. 4, 2521. doi: 10.1038/ncomms3521

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Muthukumaraswamy, S. D. (2011). Temporal dynamics of primary motor cortex γ oscillation amplitude and piper corticomuscular coherence changes during motor control. Exp. Brain. Res. 212, 623–633. doi: 10.1007/s00221-011-2775-z

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Onsager, L. (1944). Crystal statistics. I. A two-dimensional model with an order-disorder transition. Phys. Rev. 65, 117–149. doi: 10.1103/PhysRev.65.117

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 | Pubmed Full Text | CrossRef Full Text | Google Scholar

Peng, C. K., Havlin, S., Hausdorff, J. M., Mietus, J. E., Stanley, H. E., and Goldberger, A. L. (1995a). Fractal mechanisms and heart rate dynamics. J. Electrocardiol. 28, 59–65. doi: 10.1016/S0022-0736(95)80017-4

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Peng, C. K., Havlin, S., Stanley, H. E., and Goldberger, A. L. (1995b). Quantification of scaling exponents and crossover phenomena in nonstationary heartbeat time series. Chaos 5, 82–87. doi: 10.1063/1.166141

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Pfurtscheller, G. (1977). Graphical display and statistical evaluation of event-related desynchronization. Electroen. Clin. Neuro. 43, 757–760. doi: 10.1016/0013-4694(77)90092-X

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Pfurtscheller, G. (1992). Event-related synchronization (ERS): an electrophysiological correlate of cortical areas at rest. Electroen. Clin. Neuro. 83, 62–69. doi: 10.1016/0013-4694(92)90133-3

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

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

Google Scholar

Plenz, D., and Chialvo, D. R. (2009). Scaling properties of neuronal avalanches are consistent with critical dynamics. arXiv:0912.5369

Google Scholar

Poil, S.-S., Hardstone, R., Mansvelder, H. D., and Linkenkaer-Hansen, K. (2012). Critical-state dynamics of avalanches and oscillations jointly emerge from balanced excitation/inhibition in neuronal networks. J. Neurosci. 32, 9817–9823. doi: 10.1523/JNEUROSCI.5990-11.2012

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Priemer, R. (1990). Introductory Signal Processing. River Edge, NJ: World Scientific Publishing. doi: 10.1142/0864

CrossRef Full Text | Google Scholar

Priesemann, V., Munk, M., and Wibral, M. (2009). Subsampling effects in neuronal avalanche distributions recorded in vivo. BMC Neurosci. 10:40. doi: 10.1186/1471-2202-10-40

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Robinson, P. M. (ed.). (2003). Time Series with Long Memory. Oxford: Oxford University Press.

Google Scholar

Samorodnitsky, G. (2006). Long range dependence. Found. Trends Stoch. Syst. 1, 163–257. doi: 10.1561/0900000004

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 | Pubmed Full Text | CrossRef Full Text | Google Scholar

Shew, W. L., and Plenz, D. (2013). The functional benefits of criticality in the cortex. Neuroscientist 19, 88–100. doi: 10.1177/1073858412445487

Pubmed Abstract | Pubmed Full Text | 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 | Pubmed Full Text | CrossRef Full Text | Google Scholar

Shriki, O., Alstott, J., Carver, F., Holroyd, T., Henson, R. N. A., Smith, M. L., et al. (2013). Neuronal avalanches in the resting MEG of the human brain. J. Neurosci. 33, 7079–7090. doi: 10.1523/JNEUROSCI.4286-12.2013

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Siegel, M., Donner, T. H., and Engel, A. K. (2012). Spectral fingerprints of large-scale neuronal interactions. Nat. Rev. Neurosci. 13, 121–134. doi: 10.1038/nrn3137

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Singer, W. (1999). Neuronal synchrony: a versatile code for the definition of relations? Neuron 24, 49–65. doi: 10.1016/S0896-6273(00)80821-1

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

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

Google Scholar

Stam, C. J., and de Bruin, E. A. (2004). Scale-free dynamics of global functional connectivity in the human brain. Hum. Brain. Mapp. 22, 97–109. doi: 10.1002/hbm.20016

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Stanley, H., Buldyrev, S., Goldberger, A., Goldberger, Z., Havlin, S., Mantegna, R. N., et al. (1994). Statistical mechanics in biology: how ubiquitous are long-range correlations? Physica A 205, 214–253. doi: 10.1016/0378-4371(94)90502-9

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Strogatz, S. H., and Mirollo, R. E. (1991). Stability of incoherence in a population of coupled oscillators. J. Stat. Phys. 63, 613–635. doi: 10.1007/BF01029202

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

van den Heuvel, M. P., and Sporns, O. (2011). Rich-club organization of the human connectome. J. Neurosci. 31, 15775–15786. doi: 10.1523/JNEUROSCI.3539-11.2011

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Wang, F., Yamasaki, K., Havlin, S., and Stanley, H. E. (2005). Scaling and memory of intraday volatility return intervals in stock market. arXiv:physics/0511101

Pubmed Abstract | Pubmed Full Text | Google Scholar

Wang, X.-J. J. (2010). Neurophysiological and computational principles of cortical rhythms in cognition. Physiol. Rev. 90, 1195–1268. doi: 10.1152/physrev.00035.2008

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Werner, G. (2009). Fractals in the nervous system: conceptual implications for theoretical neuroscience. arxiv:0910.2741

Pubmed Abstract | Pubmed Full Text | Google Scholar

Werner, G. (2010). Fractals in the nervous system: conceptual implications for theoretical neuroscience. Front. Physiol. 1:15. doi: 10.3389/fphys.2010.00015

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Willinger, W., Taqqu, M. S., and Teverovsky, V. (1999). Stock market prices and long-range dependence. Finance Stoch. 3, 1–13. doi: 10.1007/s007800050049

CrossRef Full Text | Google Scholar

Appendix

TABLE A1

Table A1. A list of the 66 brain regions which are represented by 66 oscillators in the Cabral model.