Markers of criticality in phase synchronisation
MMarkers of criticality in phase synchronisation
Maria Botcharova , , Simon F. Farmer , & Luc Berthouze , , ∗ Centre for Mathematics and Physics in the Life Sciences and Experimental Biology, University CollegeLondon, UK Institute of Neurology, University College London, UK The National Hospital for Neurology and Neurosurgery, London, UK Centre for Computational Neuroscience and Robotics, University of Sussex, Falmer, UK Institute of Child Health, London, University College London, UK ∗ Corresponding author:
Abstract.
The concept of the brain as a critical dynamical system is very attractive because sys-tems close to criticality are thought to maximise their dynamic range of information processingand communication. To date, there have been two key experimental observations in support of thishypothesis: i) neuronal avalanches with power law distribution of size and ii) long-range temporalcorrelations (LRTCs) in the amplitude of neural oscillations. The case for how these maximisedynamic range of information processing and communication is still being made and because asignificant substrate for information coding and transmission is neural synchrony it is of interest tolink synchronization measures with those of criticality. We propose a framework for characterisingcriticality in synchronisation based on a new metric of phase synchronisation (rate of change ofphase difference) and a set of methods we have developed centred around the detection of LRTCs.We test this framework against two classical models of criticality (Ising and Kuramoto) and recentlydescribed variants of these models aimed to more closely represent human brain dynamics. Fromthese simulations we determine the parameters at which these systems show evidence of LRTCs inphase synchronisation. We demonstrate proof of principle by analysing pairs of human simultane-ous EEG and EMG time series, suggesting that LRTCs of corticomuscular phase synchronisationcan be detected in the resting state and experimentally manipulated. The existence of LRTCs influctuations of phase synchronisation suggests that these fluctuations are governed by non-localbehaviour, with all scales contributing to system behaviour. This has important implications re-garding the conditions under which one should expect to see LRTCs in phase synchronisation.Specifically, brain resting states may exhibit LRTCs reflecting a state of readiness facilitating rapidtask-dependent shifts towards and away from synchronous states that abolish LRTCs.
Keywords:
Criticality, Long-range temporal correlations, Phase synchronisation, Detrended Fluc-tuation Analysis, Oscillations, Kuramoto, Ising
The concept of the brain as a dynamical system close to a critical regime is attractive because systemsclose to criticality are thought to maximise their dynamic range of information processing and commu-nication, show efficiency in transmitting information and a readiness to respond to change (
Chialvo ,2010;
Shew et al. , 2009;
Stam and de Bruin , 2004;
Beggs and Plenz , 2003;
Shew and Plenz ,2013;
Linkenkaer-Hansen et al. , 2001;
Sornette , 2006;
Beggs and Timme , 2012;
Werner , 2009;
Linkenkaer-Hansen et al. , 2004;
Meisel et al. , 2012;
Kinouchi and Copelli , 2006).A number of modeling studies have shed important light on the behaviour of neurally inspired systemsclose to their critical dynamical range (
Kitzbichler et al. , 2009;
Poil et al. , 2012;
Shew et al. ,2009;
Breakspear et al. , 2010;
Daffertshofer and van Wijk , 2011). To date there have been twosignificant experimental observations suggesting that the brain may operate at, or near, criticality. Theseare: i) the discovery that the spatio-temporal distribution of spontaneous neural firing statistics can becharacterised 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 ofneural oscillations, typically bandpassed MEG or EEG (
Linkenkaer-Hansen et al. , 2001;
Hardstoneet al. , 2012). The mechanisms by which avalanches and LRTCs of oscillation amplitude may maximisethe dynamic range of information processing and communication are still to be fully understood andexperimental and computational neuroscience data linking the two phenomena are only just beginningto emerge (
Plenz and Chialvo , 2009;
Poil et al. , 2012)Population coding approaches to neuronal information storage and transmission show that bothchanges in the firing rate and changes in neuronal synchronisation and desynchronisation of action a r X i v : . [ c ond - m a t . d i s - nn ] A p r otentials are required to indicate changes in signal salience ( Baker et al. , 2001;
Schoffelen et al. ,2005;
Pfurtscheller , 1977, 1992;
Singer , 1999). At a coarser spatio-temporal scale, extracellular brainsignals (local field potentials, corticography, EEG and MEG), which depend on recordings within thebrain, at the brain surface and at the scalp are observed to be quasi-oscillatory (brain oscillations) andin the resting state contain spectral peaks within distinct frequency bands sitting on a 1 /f decrease inpower with increasing frequency ( Buzsaki , 2006). Brain oscillations both in the resting state and duringtask conditions show short-range and long-range synchronisation when examined both from the phaseand amplitude envelope perspectives (
Wang , 2010). Primarily neuroscience has focused on the detectionof synchronisation between areas either at zero phase lag, or with a fixed phase delay. This is in part aconsequence of the fact that the averaging necessary to extract evidence of signal correlation requires aconsistent phase relationship between the two signals for at least some period of the recording.Importantly, neural synchronisation is weak and it fluctuates spontaneously over time. A numberof experiments have shown neural synchronisation to be consistently modulated by cognitive, percep-tual and motor tasks supporting the idea that synchronisation and de-synchronisation within and acrossfrequency bands may play an important role in communication within the nervous system (
Buzsaki ,2006;
Singer , 1999;
Fries , 2009;
Schoffelen et al. , 2005;
Akam and Kullmann , 2010;
Pikovskyet al. , 2003;
Doesburg et al. , 2009;
Farmer , 1998;
Conway et al. , 1995;
Baker et al. , 1999).Changing synchronisation patterns may indicate an evolution in the relationship and exchange of infor-mation (
Pikovsky et al. , 2003). Neural synchronisation can exist between nearby and distant regions,across a range of time scales, and can be characterised using a number of techniques based on time- andfrequency-domain techniques as well as mutual information (
Buzsaki , 2006;
Schoffelen et al. , 2005;
Siegel et al. , 2012;
James et al. , 2008;
Halliday et al. , 1998;
Brittain et al. , 2009).Neuronal synchronisation occurs when the mutual influence of neurons on each other causes themto fire close together in time. It is favoured by oscillatory activity. Oscillators can be tipped in and outof weak synchonization through shared noise, a phenomenon first appreciated by Huygens (
Pikovskyet al. , 2003). Therefore weak yet variable synchrony between neuronal oscillators may easily emergewithin complex and highly interactive neural networks. In this paper the term synchronisation will beused 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 todemonstrate that non-fixed yet non-random phase relationships between signals are present in modelsof critical synchronisation and we will show that, in principle, the methodogy can be applied to neuraldata in order to further explore the relationship between neural synchronisation and systems operatingclose to a critical regime.Recent evidence supporting the idea of criticality in the dynamics of the resting state brain activityand the appreciation that synchronisation is an important extractable property of neural spatio-temporaldynamics has led researchers to ask whether neuronal synchrony can have properties consistent with adynamical system at criticality. These approaches identify power law distributions in neural synchro-nisation where synchronisation has been defined as phase consistency between two thresholded timeseries, e.g., see the phase lock interval (PLI) measure and the lability of global synchronisation (GLS)measure in
Kitzbichler et al. (2009). These findings are of considerable interest, however, the resultssupporting power law behaviour of PLI have been shown by the present authors to be vulnerable to datapooling and therefore may not provide robust estimates of critical synchronisation in neural time seriesdata (
Botcharova et al. , 2012) (see also
Shriki et al. , 2013).As discussed above, long-range temporal correlations (LRTCs) exist in dynamical systems thoughtto operate close to a critical regime (
Linkenkaer-Hansen et al. , 2001). They are typically identifiedby the autocorrelation function of the time series decaying in the form of a power law (
Granger andJoyeux , 1980). The detrended fluctuation analysis (DFA) technique allows a characterisation of LRTCsthrough an exponent similar to the Hurst exponent. DFA has been widely used in order to demonstratethe presence of LRTCs in a number of natural and human phenomena (see
Stanley et al. , 1994;
Bak ,1996;
Samorodnitsky , 2006;
Hardstone et al. , 2012;
Peng et al. , 1994, 1995b,a;
Karmeshu andKrishnamachari , 2004;
Hausdorff et al. , 1995;
Wang et al. , 2005;
Robinson , 2003, for examples).In neurophysiology, the finding of LRTCs in amplitude fluctuations of the bandpass filtered MEG andEEG (
Linkenkaer-Hansen et al. , 2001, 2004) has inspired us to develop a methodological frameworkthat can be used to to verify the presence or absence of power law scaling of detrended fluctuations andwhere power law scaling is present to estimate and ascertain non-trivial DFA exponents in the momentto moment fluctuations of phase synchronisation (quantified in terms of the rate of change of phasedifference time series) between pairs of neuronal oscillation time series.he methodology is tested as follows: i) on synthetic time series where their phase difference hasknown temporal properties with a known DFA exponent. Using these simulations we demonstrate themethod’s ability to recover known DFA exponents in the phase difference, and we test the method’srobustness 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 andtheir pairwise phase differences can be extracted. The output of these models is examined using ourmethod 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 EEGand examined with additive noise. We show from this analysis that a rise in DFA exponent associatedwith robust power law detrended fluctuation scaling occurs close to the critical regimes of both the Isingmodel and the Kuramoto model with noise.We next use our methodology to examine a system of Kuramoto oscillators, operating in a range offrequencies close to the physiological γ frequency range of MEG and EEG that are connected through anetwork 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 theparameters at which this system shows evidence of LRTCs in the rate of change of phase differences andwe relate the presence of LRTCs to the network’s connectivity.Finally, we demonstrate that in principle this methodology may be applied to neurophysiological datathrough analysing pairs of human EEG and EMG time series. These preliminary results suggest thatLRTCs can be detected in the phase synchronisation between oscillations in human neurophysiologicalrecordings.We present and discuss our methodology in detail and we offer an interpretation of its results inrelation to the emerging literature on neural synchrony and criticality within neural systems. We suggestthat the existence of a valid DFA exponent in fluctuations of a phase difference measure suggests thatthe fluctuations are governed by non-local behaviour, with all scales contributing to system’s behaviour.
We seek to characterise the presence of long-range temporal correlations in the (time-varying) phasedifference between two time series. These time series may be physiological signals such as EEG, MEGor EMG, time series extracted from a simulation or physical model, or data recorded from other natu-ral phenomena. Below, we present the detail of the various components of our proposed methodology,including a technique used to calculate phase differences, detrended fluctuation analysis (DFA) and therecently introduced ML-DFA method for validating the output of DFA. Figure 1 illustrates the applica-tion of our methodology to neurophysiological data using two sample MEG time series. We note that forthese signals, we bandpassed filter the data to a frequency band of interest, however, this step will beomitted in model data considered further in the manuscript.
The phase of a single time series s ( t ) is calculated by first finding its analytic signal: s a ( t ) = s ( t ) + H (cid:2) s ( t ) (cid:3) (1)where H (cid:2) s ( t ) (cid:3) is the Hilbert transform: H [ s ( t )] = p.v. (cid:90) ∞−∞ s ( τ ) 1 π ( t − τ ) dτ (2)and p.v. indicates that the transform is defined using the Cauchy principal value. The signal phase is defined such that it belongs to a range φ ( t ) ∈ [0 , π ] or φ ( t ) ∈ [ − π, π ]. When a singleoscillatory cycle is completed the phase returns to its starting value. A time-varying phase thereforehas the properties of a sawtooth function (see panel 3 in Figure 1). In order to turn the phase into acontinuous signal, the phase is unwrapped, so that at each discontinuity, a value of 2 π is added to thephase ( Freeman , 2004;
Freeman and Rogers , 2002). L i n e a r F i t: (cid:95) =0 . Log (ns) Log F ( n s ) (cid:239)
505 x 10 (cid:239) Time (secs) A m p ( µ V ) (cid:239)
202 x 10 (cid:239) Time (secs) A m p ( µ V ) (cid:60)(cid:47) (cid:47) Time (secs) (cid:113) ( t ) (cid:113) ( t ) (cid:239) (cid:54) (cid:113) ( t ) (cid:239)
505 Time (secs) (cid:54) (cid:113) ( t ) / (cid:54) t (cid:239)
202 Time (secs) (cid:54) (cid:113) ( t ) / (cid:54) t
1. Two time series are x ( t ) and x ( t ).2. Physiological signalsmay be filtered to adesired frequency rangeat this stage.3. Time varying phase φ ( t ) of each signal isfound by first calculatingthe analytic signal s a ( t )using Hilbert transform H : s a ( t ) = x ( t ) + H ( x ( t )) s a ( t ) = Ae iφ ( t ) φ ( t ) = tan − H ( x ( t )) x ( t ) .4. The phases φ ( t ) and φ ( t ) are unwrappedto become a smooth signal.5. Instantaneous phasedifference ∆ φ ( t ) is found.It is the tan − of H ( x ( t )) x ( t ) − x ( t ) H ( x ( t )) x ( t ) x ( t ) + H ( x ( t )) H ( x ( t ))6. The time derivate ofthe phase difference∆ φ ( t ) / ∆( t )7. Detrended fluctuation analysis(DFA) is performed onthe phase difference toascertain the presence of LRTCs.The validity of the DFA exponentis tested using the ML-DFA method.Here, the fluctuation plot is notrejected because thebest-fitting model is linear. Figure 1. Step-by-step illustration of the proposed method.
We use two sample MEG signals fromthe left and right motor cortex, displayed throughout panels 1-4 in red and blue, respectively. Panel 2 shows anoptional bandpass filtering step. In panel 3 the instantaneous phases of the two time series are calculated usingthe Hilbert transform. Panel 4 shows the unwrapped phases leading to a time-varying phase difference displayedin panel 5. In panel 6, the rate of change of this phase difference is calculated. This step is illustrated using twoplots, each showing a different time scale in the x-axis. These two time scales correspond to the minimum andmaximum window sizes used in the DFA analysis, see Section 2.4. Panel 7 shows the resulting DFA fluctuationplot. The validity of this plot is determined using ML-DFA, see Section 2.5. In this case, the validity of the DFAplot was confirmed, with a DFA exponent of 0.57. he phase difference φ ( t ) − φ ( t ) between two different time series s ( t ) and s ( t ) is calculated usingthe respective Hilbert transform of the signals H [ s ( t )] and H [ s ( t )] ( Pikovsky et al. , 2003): φ ( t ) − φ ( t ) = tan − (cid:26) H [ s ( t )] s ( t ) − s ( t ) H [ s ( t )] s ( t ) s ( t ) + H [ s ( t )] H [ s ( t )] (cid:27) (3)Full synchronisation between the two signals is indicated by a constant difference in phase over sometime 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 timeincreases. As we shall use detrended fluctuation analysis (DFA), see Section 2.4, to assess the presence oflong-range temporal correlations and DFA in its standard form assumes a bounded signal, in this paper,we characterise phase synchronisation in terms of the time derivative of the phase difference time series φ ( t ) − φ ( t ), i.e., the rate of change of the phase difference. The autocorrelation function R ss ( τ ) of a signal s ( t ) quantifies the correlation of a signal with itself atdifferent time lags τ ( Priemer , 1990), formally: R ss ( τ ) = (cid:90) −∞∞ s ( t + τ )¯ s ( t ) dt (4)where ¯ s ( t ) is the complex conjugate of s ( t ) and therefore ¯ s ( t ) = s ( t ) if s ( t ) is real-valued.In signals with short-range or no dependence ( Beran , 1994), the autocorrelation function showsa rapid decay. Gaussian white noise, for example, is a signal with no temporal dependence becauseeach successive value of the time series is independent and thus its autocorrelation function decaysexponentially. In contrast, a slow decay of the autocorrelation function indicates that correlations persisteven 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: R ss ( τ ) ∼ Cτ − α (5)where C > α ∈ (0 ,
1) are constants, and the symbol ∼ indicates asymptotic equivalence ( Clegg ,2006), then the time series is said to contain long-range temporal correlations (LRTCs). LRTCs are asubject of considerable scientific interest. They have been detected in biological data (
Samorodnitsky ,2006;
Willinger et al. , 1999;
Peng et al. , 1994;
Carreras et al. , 1998;
Berthouze et al. , 2010;
Linkenkaer-Hansen et al. , 2001) and have been discussed within the context of complex systemsoperating in a critical regime.Applying a Fourier transformation to Equation 5, a similar formulation exists for the spectral densityof the signal (
Clegg , 2006), with f representing frequency: G ss ( f ) ∼ Bf − β (6)where β = 1 − α and is also related to the level of temporal dependence.The exponents α and β in Equations 5 and 6 are connected to the Hurst Exponent, H , by α = 2 − H and β = 2 H − Taqqu et al. , 1995;
Beran , 1994).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 timelags τ where measurement errors are also largest ( Clegg , 2006). Similarly, in the frequency domain, β isbest approximated by the shape of the spectral density at large frequency shifts f . Determination of theHurst exponent for non-stationary signals is not straightforward, and therefore, for practical applications,the related property of self-similarity (see below) is considered. Detrended fluctuation analysis (DFA) may be used to determine the self-similarity of a time series (
Penget al. , 1994, 1995b). The application of DFA returns the value of an exponent, which is closely relatedto the Hurst exponent (
Clegg , 2006;
Beran , 1994). DFA is often considered to be applicable to bothstationary and non-stationary data although recent reports, e.g.,
Bryce and Sprague (2012), havesuggested that the ability of DFA to deal with non-stationary signals is overstated. In Section 2.5, wewill describe our approach to mitigating this concern.o calculate the DFA exponent, the time series is first detrended and then cumulatively summed. Theroot 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 scalingbetween the residuals (or detrended fluctuations) and the box sizes. In the log space, this power lawscaling 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 range0 . < H < < H < . H = 1 represents pink noise, and H = 1 . H = 0 . Linkenkaer-Hansen et al. ,2001). If the minimum window size is significantly smaller than this value, then the fluctuation plot willtypically contain a crossover at the window length of a single period (
Hu et al. , 2001). However, fornon-oscillatory time series for which there is no characteristic temporal scale and there are rapid changesat each innovation, such as Gaussian white noise or FARIMA time series (see Section 2.6), a smallerwindow size may be used.The maximum window size should encompass a significant proportion of the time series yet containsufficient estimates to allow for a robust estimate of the average fluctuation magnitude across the timeseries. 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 alogarithmic scaling and a minimum window of 8 time steps for simulated data, and 1 second for neuro-physiological oscillations (sampled at 512Hz, band-pass filtered 15 . − . Linkenkaer-Hansen et al. (2001) we take a maximum window sizeof N/
10 time steps where N is the length of the time series. As mentioned above, a self-similar process will produce a power law relationship between the magnitudeof the detrended fluctuations and the box sizes. In DFA, this power law scaling is characterised in termsof 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
Clausetet al. (2006) for a reasoned view as to why this may not be appropriate), however, since the object of DFAis to find evidence for or against scaling and because a valid DFA exponent can only be obtained whenthe DFA fluctuation plot is indeed linear we have introduced a model selection method for establishingthe 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 meansof confirming that a signal is self-similar, ii) a DFA fluctuation plot will necessarily increase with win-dow size, iii) an exponent may be too easily obtained through simple regression analysis producing astatistically significant result with a high r value even though the linear model may not best representa given DFA fluctuation plot, iv) the discovery of an exponent > . r value may lead tothe 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 Botcharovaet 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 assess-ing linearity in the fluctuation plots of the following time series: i) those with known combinations ofshort and long-range temporal correlations, 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., sinewaves (
Botcharova et al. , 2013).The technique fits the DFA fluctuation plot with a number of different models (see below) andcompares the fit of each model using the Akaike Information Criterion (AIC), which discounts for thenumber of parameters needed to fit the model. The DFA exponent is accepted as being valid only if thebest fitting model is linear. We want to stress that this does not equate to stating that the fluctuationplot is linear. Rather, we do not reject the linear model hypothesis. In what follows, only those time seriesfor which the linear model hypothesis is not rejected (i.e., their DFA fluctuation plot is best-fitted by thelinear model) contribute to the DFA exponents presented in the present paper and where appropriatewe indicate where linear scaling of the fluctuation plot is lost.he 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 ) = (cid:80) Ki =0 a i x i for K = { , ..., } Root - f ( x ) = a ( x + a ) /K + a for K = { , , } 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 normalise the fluctuation magnitudes with: lF scaled = 100 × lF − lF min lF max − lF min where lF min and lF max are the minimum and the maximum values of vector lF respectively. A function L is then defined: L = n (cid:89) i =1 p ( lns ( i )) lF scaled ( i ) 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 ) = (cid:12)(cid:12) f ( lns ) (cid:12)(cid:12)(cid:80) ni =1 (cid:12)(cid:12) f ( lns ) (cid:12)(cid:12) where f ( lns ) is the fitted model. Absolute values are used in order to ensure that p ( lns ) remains in therange [0 , L to produce a function that is similar in form to a log-likelihood: log L = n (cid:88) i =1 lF scaled ( i )log p ( lns ( i ))This is maximised to find the parameters a i necessary for f ( lns ). It is worth mentioning that the appli-cation of the logarithm means that the values belonging to lns are not equally weighted for all i . Thelarger window sizes have a lower weighting, which is beneficial because these estimates are also the leastrobust 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 parametersused ( Akaike , 1974;
Mackay , 2003). For a model using k parameters, with likelihood function log L , theAkaike Information Criterion is calculated using the following expression:AIC = 2 k − L + 2 k ( k + 1) n − k − k is the number of parameters that the model uses ( Akaike , 1974). An adapted formula wasproposed by
Hurvich and Tsai (1989), which accounts for small sample sizes. The model which providesthe best fit to the data is that with the lowest value of AIC. It is important to recall that the AIC canonly be used to compare models. It does not give any information as to how good the models are atfitting the data, i.e., it is only its relative value, for different models, that is important; and it would notbe possible, for instance, to compare AIC values obtained from different data sets to each other.
An Autoregressive Fractionally Integrated Moving Average model (FARIMA) (
Hosk-ing , 1981) can be used to create time series with self-similarity. The model provides a process that caneasily be manipulated to include a variable level of LRTCs within a signal, from which DFA shouldreturn the exponent used to construct the FARIMA process.To construct a FARIMA process a time sequence of zero-mean white noise is first generated, whichis typically taken to be Gaussian, and necessarily so to produce fractional Gaussian noise. The FARIMAprocess, X ( t ), is then defined by parameters p , d and q and given by: − p (cid:88) i =1 ϕ i B i (1 − B ) d X ( t ) = q (cid:88) i =1 ϕ i B i ε ( t ) . (7) is the backshift operator operator, so that BX ( t ) = X ( t −
1) and B X ( t ) = X ( t − − B ) are calculated using ordinary expansion, so that (1 − B ) X ( t ) = X ( t ) − X ( t −
1) + X ( t − d must be an integer in the ARIMA model, the FARIMA can take fractional valuesfor d . A binomial series expansion is used to calculate the result:(1 − B ) d = ∞ (cid:88) k =0 (cid:18) dk (cid:19) ( − B ) k . The left hand sum deals with the autoregressive part of the model where p indicates the number ofback-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 whitenoise to be included are q , with coefficients ϕ i . In the range | d | < , FARIMA processes are capable ofmodelling long-term persistence ( Hosking , 1981). As we will only consider p = 1 and q = 1 throughoutthe manuscript, we will refer to ϕ as ϕ and ϕ as θ . We set | ϕ | < | θ | < X ( t ) is finite ( Hosking , 1981).A FARIMA(0, d ,0) is equivalent to fractional Gaussian noise with d = H − ( Hosking , 1981). Thisproduces a time series with a DFA fluctuation plot that has been shown to be asymptotically linear witha slope of d + 0 . Taqqu et al. , 1995;
Bardet and Kammoun , 2008). By manipulating the ϕ and θ parameters, the DFA fluctuation plots can also be distorted. Surrogate Data.
Two time series x ( t ) and x ( t ) can be constructed such that the time derivativeof 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 differenceof the two time series ∆ ( φ ( t ) will be the cumulative sum of X ( t ), which is discrete in this case: ∆ ( φ ( t )) = t (cid:88) s =1 X ( s )The two phases φ i ( t ) and φ ( t ) of x ( t ) and x ( t ) respectively must be constructed to have a differenceof ∆ ( φ ( t )), or some multiple of ∆ ( φ ( t )) since DFA is unaffected by multiplying a time series by a constant.We therefore set φ ( t ) = (cid:80) ts =1 X ( s )2 f s and φ ( t ) = − (cid:80) ts =1 X ( s )2 f s where f s takes the role of a nominal samplingrate 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 definedas: x = cos( ω + (cid:80) ts =1 X ( s )2 f s )and x = cos( ω − (cid:80) ts =1 X ( s )2 f s )where ω is a constant.In what follows, we used ω = 1 and f s = 600. These values were chosen in order to produce a smoothenough phase difference. This was necessary to prevent artefacts produced by the Hilbert transform whenapplied to non-smooth data. When using physiological data, a high enough sampling rate guarantees thatthe signals will be smooth.A hundred time series X ( t ) were generated using the algorithm described in ( Hosking , 1981) foreach of the 11 DFA exponents 0 . , . , . , ...,
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 logarithmicallyspaced with a minimum of 600 time steps to correspond to f s and maximum N/
10 where N = 2 is thelength of the time series.A further control analysis was performed in which a Gaussian white noise time series η i ( t ) was addedto one of the signals, namely, x (cid:48) ( t ) = cos( ω + (cid:80) ts =1 X ( s )2 ) + η 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 supposeo be the case for noise in a neurophysiological time series. By applying the phase analysis method tosignals with additive noise, we were able to test the robustness of the method to noisy data. In thisanalysis, 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. 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. Atemperature parameter controls the collective magnetisation (
Onsager , 1944). The Ising model has beenrecently used as a model for a 2-dimensional network of connected and interacting neurons (
Kitzbichleret al. , 2009).Each element of the grid is assigned a spin p i , initially at random, which takes a value +1 (spin up)or − p is given by the Hamiltonian function H ( p ) = − JΣ Ni,j = nn ( i ) p i p j , where j is an index for the four elements that are nearest neighbours 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 thesymbol <> indicates taking the expectation value.The probability P of a given configuration occurring is then proportional to P = e − H ( p ) /kT , where T is the temperature parameter and k is Boltzmann’s constant. The system may switch into a newconfiguration if its associated probability is higher or equal to that of the current configuration. TheIsing 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 2for an example of an Ising model lattice). With increasing temperature values, the probability of a spinchanging increases. As the system temperature increases the spins change more rapidly and the systembecomes increasingly disordered and corresponds to a non-magnetic state (Figure 2A). The temperaturevalue at which a transition occurs between the magnetised and non-magnetised states is known as thecritical temperature T c . At this temperature (see Figure 2B), the system will have a large dynamic rangeand infinite correlation length. However, in practice, this means that the system contains spin clusters ofall 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 elementsat T c . A . T empe r atu r e: 1 .
10 20 30 40 50 60 70 80 90102030405060708090 B . T empe r atu r e: 2 .
10 20 30 40 50 60 70 80 90102030405060708090 C . T empe r atu r e: 100000
10 20 30 40 50 60 70 80 90102030405060708090
Figure 2. The Ising model lattice at a single time point once steady state has been reached for 3different values of the temperature parameter.
A. The Ising lattice at a cold temperature of 1 .
5. Almostall spins are aligned (white) and there is little change across time. C. The Ising lattice at a high temperatureof T = 10 . The spins form a more or less random pattern across the lattice. B. The Ising lattice near criticaltemperature, T = 2 .
3. The lattice contains clusters of spins that are both small and large. Note that these aresnapshots 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 2-dimensional Ising model in Onsager (1944), and is given by the solution to the equationsinh (cid:18) JkT c (cid:19) = 1n the implementation of the Ising model used here, the lattice consists of 96 ×
96 elements. Theconstants J and k are set to J = 1 and k = 1 without loss of generality, which gives the criticaltemperature T c = ln (1+ √ ≈ . Kitzbich-ler et al. (2009). Namely, the lattice is divided into a number of smaller square lattices, which we referto 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 × Kitzbichler et al. (2009), but we also investigatedother sub-lattice sizes (results not shown) in order to verify that this choice of sub-lattice size did notaffect the results. Indeed, previous work by
Priesemann et al. (2009) suggests that the sub-samplingof a system may cause it to be mis-classified as sub-critical or supercritical when it is in fact in a criticalstate.Pairs of time series, for every possible pairing of sub-lattices belonging to the larger grid, were usedas input signals for the phase analysis method. For the sub-lattice of size 8 × ,
296 pairings.
The Kuramoto model.
The Kuramoto model is a classical model of synchronisation (
Acebr´on et al. ,2005;
Chopra and Spong , 2005) and has been used to study the oscillatory behaviour of neuronalfiring (
Breakspear et al. , 2010;
Kitzbichler et al. , 2009;
Pikovsky et al. , 2003) among many otherbiological systems.The Kuramoto model describes the phase behaviour of a system of mutually coupled oscillators witha 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 throughcoupling 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):˙ φ i ( t ) = ω i ( t ) + KN Σ Nj =1 sin( φ j ( t ) − φ i ( t )) (8)Because the Kuramoto model provides an equation governing the phase evolution of each oscillatorin the system, there is no need for the Hilbert transform to recover the phase time series and thereforeonly 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 meanfield parameters, which result from the combined effect of all oscillators in the system. Namely, we maywrite: ˙ φ i ( t ) = ω i + Kr ( t )sin( ψ ( t ) − φ i ( t )) (9)where ψ ( t ) is the mean phase of the oscillators, and r ( t ) is their phase coherence, so that: r ( t ) e iψ ( t ) = 1 N N (cid:88) j =1 e iφ j ( t ) (10)This crucially indicates that each oscillator is coupled to the others through its relationship withmean field parameters r ( t ) and ψ ( t ), so that no single oscillator, or oscillator pair drives the process ontheir own. The oscillators synchronise at a phase equal to the mean field ψ ( t ), and r ( t ) describes thestrength of synchronisation, sometimes referred to as the extent of order in the system ( Strogatz andMirollo , 1991;
Bonilla et al. , 1992). When r ( t ) = 0, no oscillators are synchronised with each other.When r ( t ) = 1, all oscillators are entrained with each other.One solution to Equation 9 is r ≡ N → ∞ , some further deductions can be made,including the fact that when the natural frequency distribution g ( ω ) is unimodal and symmetric, anothersolution can be found for ω i , with r ( t ) not equivalent to 0 ( Kuramoto , 1975). A critical bifurcationoccurs 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 ( D¨orflerand Bullo , 2011;
Strogatz and Mirollo , 1991). The coupling at the bifurcation is referred to as thecritical coupling K c ( D¨orfler and Bullo , 2011).In an infinite Kuramoto model, criticality is defined through this point of bifurcation. For a finitesystem, however, the critical point can only be approximated by this theoretical value. One definingcharacteristic of the critical coupling for the Kuramoto system is that the greatest number of oscillatorsome into synchronisation at this value. In our study, we deal with finite-sized implementations of theKuramoto model, and we use this characteristic as a marker of the onset of critical regime in addition tothe theoretical value K c . Specifically, we use a measure characterising the onset of synchronisation withincreasing 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 frequencyand the mean phase ω i − ψ (in modulus), i.e., | ω i − ψ | < Kr , then oscillator i will synchronise to themean field ( Mertens , 2011). Thus the value of K at which Kr increases maximally is the couplingvalue 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:˙ φ i ( t ) = ω i ( t ) + KN N (cid:88) j =1 sin( φ j ( t ) − φ i ( t )) + η i ( t ) (11)where η i is a noise input taken to be uncorrelated Gaussian noise with zero mean ( (cid:104) η i (cid:105) = 0) andcovariance σ i /T ( (cid:10) η i ( t ) (cid:11) (cid:10) η j ( s ) (cid:11) = δ ij δ ( t − s ) σ i /T ) where δ ij is the Kronecker delta, δ ( t − s ) is the Diracdelta function, σ i is in radians and T = 1 second here.This creates a richer structure in the oscillator dynamics, which we suggest may better reflect couplingof neurophysiological oscillators. Furthermore, it has been shown that addition of noise increases thecritical regime over a wider range of coupling values ( Breakspear et al. , 2010). This may allow for thefluctuations of phase difference of a given oscillator pair to persist for longer with increasing couplingbefore full synchronisation is achieved.In
Daniels (2005), the critical coupling for the infinite Kuramoto model with added noise K c,noise is calculated to be: K c,noise = 2 (cid:34)(cid:90) ∞−∞ σ ω σ ω + ω g ( ω ) dω (cid:35) − Again, as the number of oscillators is inevitably finite, this value is only an approximation to the truecritical coupling in the system, but we find it useful and it is displayed alongside plots of ∆ ( Kr ), whichalthough originally introduced for a noiseless model, remains a helpful marker of the effective criticalcoupling 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 Equa-tion 11. Each time series was 6100-timestep long. The standard deviation σ i was set to 0 .
32. The distri-bution of natural frequencies was g ( ω ) ∼ N (44 π, σ ω ), with standard deviation σ ω = 15. This correspondsto a normal distribution centred around 22 Hz (which is a unimodal distribution). In order to get anidea of the spread of the distribution, the minimum natural frequency selected from this distribution was16 . . β -bandof EEG, MEG and EMG oscillations ( Farmer , 1998).For these parameter values, the critical coupling K c is equal to: K c = 2 √ √ π σ ω ∼ . K c,noise is not analytically calculable for a normal distribution g ( ω ) ∼ N , but empiricalcalculation yields: K c,noise ∼ . The Cabral model.
The third model that we consider in this paper was developed by Joanna Cabraland her colleagues, referred to as the Cabral model. It is a modification of the Kuramoto model, com-bining the dynamics of the Kuramoto oscillators with the network properties observed in the humanbrain (
Cabral et al. , 2011).The Cabral model includes a noise input to the Kuramoto oscillators and situates the 66 oscillatorson a connectivity matrix with varying connection strengths and time delays based on empirical measure-ments of 998 brain regions, which have been down-sampled to 66 (
Honey et al. , 2009). The list of brainregions considered in this model are given in the supplementary material of
Cabral et al. (2011) andre reproduced in the Appendix to the present paper. Specifically, Equation 8 is modified to include aconnectivity term C ij between oscillators j and i , namely,˙ φ i ( t ) = ω i ( t ) + KN Σ Nj =1 C ij sin( φ j ( t − D ij ) − φ i ( t )) + η i ( t ) (12)where η i is the noise input previously introduced, and D ij is the time delay associated with the linkbetween oscillators j and i . The matrix of delays D is extracted from a matrix of empirical distances L between regions using: D ij = (cid:104) D (cid:105) L ij (cid:104) L (cid:105) and is used to encode the length of time taken by neural activity to traverse the connection space. Theconnectivity and distance matrices ( C and D , respectively) can be visualised through the schematic dia-gram shown in Figure 3. The thickness and colour of the lines in the diagram represent the weights of theconnections between the oscillators representing individual brain regions. These weights are proportionalto the number of fibres 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 inTable 2 in the Appendix.
Figure 3. Schematic plot (top view) of the Cabral human brain model showing the connectionsand connection weights between oscillators which correspond to different brain regions.
The weightof the connection lines represent the strength of connectivity between the oscillators. The darkest blue lines arethe strongest 1% of connections. The node colours represent oscillators, which model different brain regions asdetailed in
Cabral et al. (2011). Colours 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 toa hemodynamic model and bandpass filtered. The resulting time series were compared to recordings ofBOLD fMRI signals using Pearson’s correlation coefficient and mean squared error to determine theparameter values K and (cid:104) D (cid:105) that generated the time series which most closely approximated the BOLDdata.In this model, there is no theoretically derived value of critical coupling and ∆ ( Kr ) is only a markerof effective change in coupling that may or may not be critical. We interpret a rise in ∆ ( Kr ) as anincrease in order of the system similar to that observed by Kitzbichler et al. (2009).he 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 drawnfrom a normal distribution with g ( ω ) ∼ N (120 π, σ ω ) with standard deviation σ ω = 5, which correspondsto a normal distribution centred 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 Cabralet al. (2011) for details).The standard deviation σ i of the noise input was set to 1 .
25. It was found that values of σ i < K and (cid:104) D (cid:105) . The value (cid:104) D (cid:105) = 11 is taken asin Cabral et al. (2011).
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 behaviour.
Table 1. Cluster information. The 66 oscillators of the Cabral model can be separated into 6clusters, based on their mutual connectivity and distance matrix patterns, and a final set of 12oscillators, which are not considered to belong to a cluster, but are grouped together here forconvenience. The table also states the average sum of weights per node belonging to each clusterand the average number of connections per node (both to 2 d.p.).
Clusters Oscillators Average weight per node Average degree distributionCluster 1 7-17 0.29 19.09Cluster 2 18-22 0.16 15.80Cluster 3 23-26,41-44 0.30 21.00Cluster 4 27-40 0.34 21.71Cluster 5 45-49 0.15 15.60Cluster 6 50-60 0.27 18.73Individual Oscillators 1-6,61-66 0.03 08.59
Disruptions to the Cabral model.
In order to investigate the role of connectivity in sustaining LRTCsof rate of change of phase difference, we modified the connectivity matrix C in the Cabral model in twoways. First, beginning with the empirical connectivity matrix we deleted any connection that extendedfrom one hemisphere into the other. We preserved all the other elements of the model’s connectivity andoscillator characteristics.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 de-scribed in Gionis et al. (2007);
Hanhij¨arvi et al. (2009). Specifically, a list of the outgoing weightsof each oscillator was made alongside the node from which it extends. Two weights were selected fromthis list. If they did not belong to the same node, then the nodes were connected to each other withthe associated outgoing weights that were selected. These weights were then deleted from the list. Tocontinue the algorithm, two further weights were selected. After the first step, it was necessary to checkat each iteration that the nodes were not already connected before connecting them. If the nodes wereconnected, 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 thosederived from the disconnected hemisphere model and standard appropriately connected model allowedus to determine the extent to which a realistic connectivity matrix of the human brain predisposes thesystem to LRTCs in the rate of change of the phase difference between the oscillator pairs representingdifferent brain regions.
A note on notation.
From this point in the text, all instances of oscillator phase φ i ( t ) and r ( t ) will bewritten as φ i and r for ease of notation, unless stated otherwise. Any quantities that are defined usingthe phases of one or more oscillators are also implicitly functions of time, although the t is omitted forthe same reason. igure 4. Schematic plot (top view) showing the connections and connection weights betweenoscillators 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 areleft unchanged. B. The connections and weights of each node are assigned randomly, but the degree distributionand weight distribution at each node is kept constant. The weight of the connection lines represent the strengthof connectivity between the oscillators. The darkest blue lines are the strongest 1% of connections. The nodecolours represent oscillators, which model different brain regions as detailed in
Cabral et al. (2011) and areidentical to Figure 3. Colours are consistent for homologous regions in the left and right hemispheres. Anteriorand posterior, left and right are shown. .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 whilsta healthy adult subject performed a 2-minute 10% MVC (maximum voluntary contraction) isometricabduction of the index finger of the right hand. The EMG was recorded using bipolar electrodes situatedover the first dorsal interosseous muscle (1DI). The EEG was recorded using a modified Maudsley montagefrom 24 Ag/AgCl electrodes with impedance < kΩ . The data were amplified and bandpass filtered4 − β frequencyrange (15 . − . The signals described in Section 2.6 were analysed. The scatter plot presented in Figure 5A shows theDFA exponents of the rate of change of phase difference expected from the construction of a FARIMAtime series with known parameters against those recovered by applying the phase analysis method. Thescatter plot shows a strong linear relationship between the expected and recovered exponents with aslope of 0.998. The fact that the slope is slightly < K nown D F A exponen t R ec o v e r e d D F A ex pon e n t Slope = 0.998
Figure 5. Plot of the recovered against the true DFA exponent for FARIMA time series.
Therelationship between recovered and true DFA values is well-approximated by a linear trend with a slope of almost1. 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 isfound to be reduced. Figure 5B shows that as the noise level is progressively increased, the percentagedifference between the known DFA exponent and that recovered by the method increases. When the noiselevel is above one which causes the percentage difference between known and recovered DFA exponentto exceed approximately 5% (note, as shown in Figure 5B, that this noise level depends on the exponent,e.g., 0 . .
025 for exponent of 0 . ≈ . − .
4, noise dominates thesignal and valid exponents are once again obtained. These exponents are at or close to 0.5 regardless ofthe 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. No i s e L eve l D F A ex pon e n t A . R ecove r ed D F A exponen t Undefined region:non − linear by ML − DFA 1.000.950.900.850.800.750.700.650.600.550.50True DFA exponent 0 0.05 0.1012345678910 No i s e L eve l % d i ff e r e n ce B . % D i ff e r ence be t ween r ecove r ed and t r ue D F A exponen t Undefined region:non − linear by ML − DFA 0.500.550.600.650.700.750.800.850.900.951.00True DFA exponent
Figure 6. True and recovered DFA exponents for noisy signals with LRTCs.
A. Recovered DFAexponent values as noise is progressively added. For each of the DFA exponents given in the legend (box insert),a signal x (cid:48) ( t ) was constructed with a noise level σ ∈ [0 , x axis. The phase synchrony analysismethod was applied to x (cid:48) ( t ) and x ( t ). This was performed 100 times. For DFA exponents corresponding toDFA fluctuation plots that were accepted as linear by ML-DFA, the average value for the 100 signal pairs isshown. There are no data points corresponding to the intermediate noise level of ≈ . ≈ . σ ∈ [0 , .
1] are shown. The colours represent different true DFA exponent values, as indicatedby the legend within the inserted box. The dashed line indicates a 5% difference between known and recoveredexponents. When the difference between the known and recovered exponent exceeded approximately 5% for anyvalue of the true exponent, the DFA fluctuation plot is not accepted as being linear by ML-DFA and thereforethe exponent is not shown on the plots.
Figure 7 shows the results for sub-lattices of size 8 ×
8. At a high temperature of T = 10 , the averageDFA exponent across all pairwise comparisons is 0 .
57 (see magenta shaded bar). This value is in excess of0 . . T = 2 .
55 (see magenta shaded bar) indicating maximal LRTC just before the critical temperature isreached.The change in mean DFA has to be seen within the context of the validity of the DFA fluctuationplots. As the system cools towards the critical point the validity of DFA exponents across all pairwisephase differences drops abruptly. The first temperature value for which < T = 2 .
75 shown as magenta shaded bar. There is a large fall in DFA fluctuation plot validity as thecritical temperature is reached (56% to 34%). This fall in validity reflects the onset of full synchronisationbetween a number of the time series. At the critical point, T = T c which occurs between T = 2 .
25 and T = 2 . .
64. As the Ising model cools below the critical point the DFA validity in general falls and there areno valid DFA fluctuation plots below T = 2 .
15. As discussed above this occurs because of the loss offluctuations in the rate of change of phase difference due to full synchronisation.Results obtained for sub-lattice sizes of 32 ×
32, 16 ×
16, 12 ×
12 and 6 × The group average results for the Kuramoto model are shown in Figure 8. As can be seen, the peakaverage DFA exponent occurs on average at K ≈
22. The value of the average DFA exponent at thiscoupling value is 0 .
65 with standard deviation 0 .
06, consistent with the rate of change of phase differenceshowing LRTCs. The peak DFA exponent occurs one coupling value later than the peak of the ∆ ( Kr )measure, at K ≈ ∆ ( Kr ) represents the coupling value at which the order parameter r increases most, igure 7. Average DFA exponents of rate of change of phase difference between pairs of timeseries generated by × sub-lattices of the × 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 singlestandard deviation from the mean. The proportion of valid exponents, as calculated by ML-DFA, is denoted bythe vertical bars. The theoretical critical parameter T c is indicated by a red asterisk. A horizontal dashed lineat DFA exponent 0.5 is plotted to guide the eye. Validity bars that are referred to in the text are highlighted inmagenta. 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 frequenciesof 15. The theoretical critical coupling K noise when noise is added is marked with a blue asterisk. The averageDFA exponent, order parameter r , its difference ∆ ( Kr ) and the proportion of valid DFA fluctuation plots fromthe full set of 199000 pairs are shown. Validity bars that are referred to in the text are highlighted in magenta. nd the point of greatest oscillator coupling flux in the system ( Kitzbichler et al. , 2009). The peakcoupling value ∆ ( Kr ) and the maximum DFA values are just less than the theoretical critical couplingof the infinite Kuramoto system with noise K c ≈ .
85. Again, these results must be understood incontext of DFA fluctuation plot validity which is 42% of the 199000 oscillator pairs at K ≈
22. Once fullsynchonization occurs between an individual pair of oscillators, their phase difference takes a constantvalue. ML-DFA detects the resulting loss of scaling by indicating that the DFA fluctuation plot is nolonger linear.After the peak DFA at K ≈
22, further increase in K eventually causes full synchronisation betweenall individual oscillator pairs. Across the whole system, fewer than 10% of oscillator pairs yield a validDFA after the critical coupling is exceeded. When all oscillator pairs are synchronised with each other,the order parameter of the system approaches its maximum level of 1 but the DFA fluctuation measureof 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 phasedifference between oscillator pairs occur when the system is in a state of maximal flux just prior to theonset of full synchronisation. Individual oscillators pairs.
Further insights into the rate of change of phase difference fluctuationbehaviour can be obtained from DFA of individual oscillator pairs. Analysis of a set of 5 oscillator pairsis shown in Figure 9. The top panel shows the change in DFA exponent with coupling K for a pairwhose initial frequencies are very close (0 . ≈ . K . At low coupling K , the oscillators do notinteract with each other and each evolves at its own natural frequency. The order in the system is lowand the DFA exponent ≈ . ≈ . K is increased, the DFA exponents of each of the oscillator pairs rise untila peak is reached. The value of K at which a maximal valid exponent is retrieved for these peaks isrelated to the difference in natural frequencies of the two oscillators as well as their interactions withthe noise and the mean field. Oscillator pairs which start further apart in frequency terms develop fullsynchonization later than those whose initial frequencies are close together. As K increases the DFAexponent of the rate of change of phase difference increases. The pairs with the strongest LRTCs on thebasis of the highest DFA exponent value prior to onset of full synchronisation are those with the greatestinital frequency difference. Increasing temporal order of the rate of change of phase difference prior tofull synchonization of these pairs may indicate a state of pre-synchronisation in these pairs. For the Cabral model we present results regarding both the global behaviour of the system throughaverage DFA exponents across all possible pairs of oscillators (Figure 10) and the behaviour of thesystem at cluster level through average DFA exponents of intra-cluster pairs of oscillators (Figure 11).
Global behaviour.
The model introduced by
Cabral et al. (2011) is affected by rich interplay betweenthe 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 thecoupling in the system is increased. These average exponent values indicate the presence of LRTCs inthe rate of change of phase difference. The peak values of mean DFA exponent correspond to peaks in thechange in order paramenter ( ∆ ( Kr )) derived for the classical Kuramoto model and the Kuramoto modelwith noise, see Kitzbichler et al. (2009) and Figure 8. Such peaks occur when the system undergoesthe greatest change in synchronisation. The peak in ∆ ( Kr ) corresponds closely to the coupling valuethat 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 differenceis equal to 100% when there is no coupling in the system (magenta shaded bar at K = 0), but it falls ascoupling 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 neighbouring coupling values (magenta shaded bar at K = 6).
10 15 20 25 30 350.50.60.7 A B Coup li ng K D F A Coup li ng K V a li d ? C D
Coup li ng K D F A Coup li ng K V a li d ? E F
Coup li ng K D F A Coup li ng K V a li d ? G H
Coup li ng K D F A Coup li ng K V a li d ? I J Coup li ng K D F A DFA of Phase Synchrony Validity
10 20 30NY
Coup li ng K V a li d ? Figure 9. Representative relationship of DFA exponents to the coupling parameter K for selectedoscillator pairs in the Kuramoto system. Panels A,C,E,G,I and K show the value of valid DFA exponents,while panels B,D,F,H,J and L 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 inthe legend of panels A,C,E,G,I and K. The first number is the difference in natural frequency (in Hz), and thesubsequent pair of numbers identifies which oscillators are being analysed. igure 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 normallydistributed natural frequencies with mean 60 Hz and standard deviation σ i = 1 .
25. The connectivity and timedelay matrices are set from empirical values. The average of the valid DFA exponents is shown in magenta andthe proportion of valid exponents, as calculated by ML-DFA, are indicated by bars. The Kuramoto model orderparameter r is in blue, and the quantity ∆ ( Kr ) is in cyan. The peak ∆ ( Kr ) has been used as an indicator of theeffective critical coupling. A horizontal line at DFA exponent 0.5 is plotted to guide the eye. The proportion ofvalid DFA bars for K = 0, K = 1 and K = 6 have been shaded in magenta. Cluster behaviour.
At coupling value K = 6, the value at which the global behaviour shows peak DFAvalue, the intra-cluster results indicate that only cluster 4, consisting of oscillators 27-40, shows validnon-trivial DFA exponents. These exponents are consistent with the presence of LRTCs. This suggeststhat cluster 4 acts as an organising force in the system when the system is in its greatest state of flux, asdemonstrated by a large increase in the order parameter. This cluster corresponds to the most connectedbrain regions listed in Table 1 and Table 2 in the Appendix.The connectivity and distance matrices for the Cabral model are shown in Figure 12. The linearcoupling between oscillators for two values of K is shown in Figure 13. The central cluster of oscillatorswith high levels of synchronization is evident from the two correlation matrices. At K = 6 (panel A), i.e.,the value at which LRTCs are detected in the rate of change of phase difference, the central oscillatorcluster shows evidence of synchronization but with Pearson correlation values of < .
0. As K increases to18, the value identified by Cabral et al. (2011) as best approximating human brain resting state BOLDfluctuations, it can be seen from Figures 10 and 11 that the proportion of oscillator pairs with valid DFAfluctuation plots is low (approximately 5%). Those oscillator pairs that remain and show persistentlyvalid 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 (seeFigure 11). At K = 18, the Cabral model shows strong cluster synchronisation. In particular, the centralcluster 4 (oscillators 27 −
40) which contains homologous elements connected across the corpus callosumshows Pearson correlation values close to 1.0 indicative of full synchrony (Figure 13B). Therefore theresults we obtained from the Kuramoto model with noise and those derived from the Cabral model aresimilar. Both show valid DFA fluctuation plots with LRTCs of the rate of change of phase difference ata coupling value where ∆ ( Kr ) is increasing and loss of validity as full synchronisation takes over. Asdiscussed earlier, ’criticality’ is not defined for the Cabral model but with increasing K there is clearlya change in the system’s order which is detected through our method.Figure 14 shows the DFA exponents of the rate of change of phase difference between individual pairsof oscillators in the form of a symmetric lattice of size 66 ×
66, where each element in the lattice representsa brain region as detailed in Table 2 of the Appendix. Panel A of this figure shows the importance ofthe central cluster in generating LRTCs of phase synchronisation. Importantly it shows this cluster’sinfluence over many of the other oscillators in the Cabral model. Cluster group 4 has the greatest sum of igure 11. Average DFA exponent for intra-cluster pairwise phase differences with increasingcoupling parameter K . Where no DFA value appears for a particular cluster, this indicates that there are novalid DFA exponents for the pairwise phase difference within that cluster. The final cluster, which is labelledindividual oscillators, consists of a set of nodes that do not fit into any of the clusters as determined by theconnectivity and distance matrices but are grouped together to demonstrate their relationship with each other. weights per oscillator and the greatest number of connections per oscillator (see Table 1). The correlationbetween the number of connections of a given oscillator and the average DFA exponent of its rate ofchange of phase difference with all other oscillators is 0 . Comparison of the three connectivity structures.
In the Cabral model, the ∆ ( Kr ) measure has itspeak at coupling value K = 6. Here, we compare the effects of the three connectivity matrices introducedin Section 2.7 on the DFA exponents of the pairwise phase difference between oscillators at this couplingvalue in Figure 14.The empirical connectivity matrix showed large DFA exponents indicating the presence of LRTCsat this coupling value for a small number of hub oscillators belonging to cluster 4 (see above). Theseoscillators have a high number of connections and large weights associated with these connections (seeTable 1). When the two hemispheres are disconnected, we see no LRTCs in the DFA exponents of thephase difference at this coupling value. When the distance matrix is preserved, but the connectivity andassociated weights are assigned at random, LRTCs are still present in the DFA exponent of the phasedifferences between oscillators, but a lower value of DFA exponent is obtained. There is no apparentcluster formation when connectivity is random. Neurophysiological data.
Figure 15 illustrates the application of our phase synchrony analysis tech-nique to the human neurophysiological data described in Section 2.8. In this example, a valid DFAexponent of ≈ . The aim of this paper is to introduce a new methodology for eliciting a marker of criticality in neuronalsynchronisation. This methodology relies on the rate of change of the phase difference between two signalsas a (time-varying) measure of phase synchronisation. The presence of long-range temporal correlationsin this quantity is proposed as marker of criticality and is assessed using detrended fluctuation analysis igure 12. Connectivity and distance matrices for the Cabral model.
Each oscillator number representsa brain region, which is defined in Table 2 in the Appendix. An empty (white) element means that the two regionsare not connected. Regions are not connected to themselves so that the diagonals are white. Panel A shows thepairwise connection matrix C between the 66 oscillators. Panel 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 normalised by row. The values associated with the colours of the plotsare defined by the colour bars. Red colours in panel A represent higher weights. Red colours in panel B representlonger distance connections. O s c ill a t o r Numbe r O s c ill a t o r N u m b e r A. Correlation matrix for K=6
10 20 30 40 50 60102030405060 − − O s c ill a t o r Numbe r O s c ill a t o r N u m b e r B. Correlation matrix for K=18
10 20 30 40 50 60102030405060 − − Figure 13. Correlation matrices for all pairs of time series generated by the Cabral model for twocoupling values K . A. K = 6 and B. K = 18, which corresponds to the oscillator correlation matrix in Cabralet al. (2011). The plots show the value of the Pearson correlation coefficient between all pairwise combinationsof the 66 oscillators used in the model.
Figure 14. DFA exponent of the rate of change of phase difference between all pairs of oscillatorsin the Cabral model at coupling K = 6 in three scenarios. A: For the empirically observed connectivitymatrix of the Cabral model. B: For a connectivity matrix representing disconnected hemispheres. C: for randomconnectivity. Empty (white) elements denote pairs for which no valid DFA exponent was found. − T i me ( s ec s ) A m p ( µ V ) A . EE G − − T i me ( s ec s ) A m p ( µ V ) B . E M G − T i me ( s ec s ) A m p ( µ V ) C . Fil t e r e d EE G 16 ⌧
24 Hz − T i me ( s ec s ) A m p ( µ V ) D . Fil t e r e d E M G 16 ⌧
24 Hz − T i me ( s ec s ) E . R a t e o f chang e o f pha s e d i ff e r enc e n F ( n ) F. Linear Fit: (cid:97) =0.59502
Figure 15. Illustration of the method with simultaneous EEG/EMG data.
A and B: One second ofsimultaneously recorded EEG and EMG, respectively. C and D: The signals after bandpass filtering in the β range 16 − ≈ .
60, indicating the presence of LRTCs.
DFA) in combination with the recently proposed ML-DFA, a heuristic technique for validating theoutput of DFA. With these methods, we can first determine the presence or absence of power law scalingusing ML-DFA and secondly the presence or absence of long-range temporal correlations (LRTCs) in thephase synchronisation of two time series based on the value of the DFA exponent. If the method returnsan exponent of ≈ .
5, this indicates a phase relationship similar to white Gaussian noise, however, if theDFA exponent is greater than 0.5, this indicates the presence of LRTCs. Importantly, we can attributesignificance to the loss of power law scaling within the fluctuation plot and draw conclusions based onan exponent value only when the exponent has been recovered from plots that are judged to be valid byML-DFA.
Surrogate Data
It was found that the phase synchrony analysis method recovers a known DFA exponent value in therate 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 synchronisation was perturbed with an additive noise source,it was found that a percentage difference between the true and recovered DFA exponent of above ap-proximately 5% noise caused DFA exponents to be judged as invalid by ML-DFA. When the surrogatedata was characterised by a DFA exponent close to 1, the recovery of this exponent using DFA was moreresistant to noise when compared to surrogate data with a lower DFA exponent of 0.6 (Figure 6). Inthese simulations we used additive noise which was included at the amplitude stage of the surrogate timeseries prior to extraction of the phase using the Hilbert transform.
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, inapplying our method to the Ising model, both of these hypotheses were not fully realised. It was foundthat when the temperature was increased to a very high level of T = 10 , the DFA exponent of the rateof change of phase difference did not fall to 0.5, but remained at ≈ .
57. This did not change when thetemperature was set to an even higher value of T = 10 . This was not a finite size effect of the system, asthe result held when larger lattice sizes (up to 1000 × T is based only on our intuition concerningthe operation of the system. As all elements in the Ising lattice interact with their neighbours it ispossible that some temporal correlation in the rate of change of phase difference may persist regardlessof 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 butwas maximal at T = 2 .
55, just in excess of the critical temperature. As can be seen in Figure 7, theconsistent change in the DFA value and the change in power law scaling behaviour indicates that thephase synchrony analysis method is capturing an important behaviour of the system close to its criticalregime. However, it is important to realise that unless an experimental neuroscientific paradigm can bediscovered that produces similar consistent changes in this measure, neurophysiological data will have tobe intepreted with caution, i.e., we may be able to state that for a given pair of neural oscillation timeseries there exists power law scaling with a DFA exponent indicative of LRTCs in the rate of change oftheir phase difference but we may not know whether for this neural state there may exist other higher (orlower) exponent values. In other words, the technique may provide evidence that the system is orderedin ways that are similar to systems nearing their critical regime but whether the technique will pinpointthe most critical regime in a neural system is open to question. We will consider this further in ourdiscussion of the results of analysing a Kuramoto system with noise.Interestingly, the evolution of the DFA exponent with the temperature parameter shares a key charac-teristic 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 temperatureis increased towards the critical T = T c and a gradual descent as the temperature rises significantly. Itwould be of interest to further assess the extent to which the proposed method captures information flowin the system, e.g., through a comparison of both methods when applied to the Kuramoto model. he Kuramoto model
In the Kuramoto model, the critical transition is characterised in terms of a global order parameter whichreflects the overall organisation of the system. However, through our phase synchrony analysis methodwe are able to make observations at a pair-wise level of Kuramoto oscillators always bearing in mind thateven at the pair-wise level the result is influenced by the oscillators’ interactions with all other oscillatorsin the model. As individual Kuramoto oscillator pairs become fully synchronised, their rate of change ofphase difference no longer contains moment-to-moment fluctuations and thus power law scaling in theDFA measure is lost. This is an important consideration because it emphasises the difference between ourmethod and more standard measures of neural synchrony. Methods for detecting neural synchrony relyon phase consistency to allow averaging out of fluctuations so that a measure of coupling (e.g., coherenceand phase coherence) is obtained. In contrast, the method introduced in this paper is dependent on thefluctuations of the two phase signals and their interaction. Therefore our method detects ’order’ acrosstime 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 synchronisation is well illustrated bothfor 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 atvalues of K just below theoretical critical coupling value. At these values of K , a power law scaling existsfor the rate of change of phase difference, and the DFA exponent of oscillator pairs with different initialfrequencies indicates the presence of LRTCs. At the onset of full synchronisation the number of oscillatorpairs for which DFA is valid drops yet those whose phase differences still possess fluctuations continueto show LRTCs. Once the critical regime has been fully crossed and the order parameter r approaches1, the DFA of the rate of change of phase difference is no longer valid for any oscillator pair.The LRTC behaviour is also clearly explained as the coupling value K decreases towards zero. Ascan be seen in Figure 8, the DFA exponent of the pairwise rates of change of phase difference decreasestowards 0.5 and yet scaling remains valid. These changes in DFA exponent are evident both on the globallevel in the average DFA and for individual oscillator pairs. At K = 0 the phases are independent fromone another yet contain noise; thus the rate of change of phase difference time series contains innovationsthat are random across time with a DFA which is valid and returns the expected exponent of 0.5. Order within the Ising and Kuramoto models
In these models, temperature T (Ising) and coupling K (Kuramoto) play a similar role in controllingthe order within the two systems, and the DFA validity and exponent results obtained from analysis ofrate 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 synchronised state with increasing K . Similarly in the Isingmodel, there is a transition from a very ordered to a disordered system with increasing T . In the humanbrain, we are not able to characterise the system by incrementally tuning a parameter and observing theresult, and we are only privy to snapshots of the working system. However, we can begin to understandthe behaviour of the brain within this range of behaviours by comparing the DFA of the rate of changeof phase difference of pairs of neurophysiological signals to the outcomes of these models of criticality. The Cabral model
We found that LRTCs exist in the rate of change of phase difference between oscillator pairs at parametervalues close to those at which the change in order, ∆ ( Kr ), increases sharply. Extrapolating from theKuramoto model with noise, we suggest that there are important changes in the order of the phasesynchronisation of interacting oscillators in the Cabral model that involve the presence of LRTCs whenthe 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 therange of coupling values 0 −
20. It approaches a level of ≈ . K approaches 20 with maximal rateof change at K ≈
6. Further analysis of the Cabral model indicates that r will gradually reach a valuecloser to 1 as K increases above a value of 60, as seen in Figure 4 of Cabral et al. (2011). Cabralfocussed her attention on K = 18 at which point the model, when fed through the Balloon-Windkesselhemodynamic model, produced an output that closely matched the spatio-temporal correlations seen inthe BOLD signals of the resting state fMRI. We find that at this value, there are no LRTCs detectablein the rate of change of phase difference measure. he 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 importantto note that LRTCs exist for a broader range of coupling values K . This finding agrees with a recentstudy by Moretti and Mu˜noz (2013) in which the authors demonstrated that a network with complexconnectivity, such as that of the Cabral model and, indeed, that of the brain, causes the critical point tobecomes 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 whichhas 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 andMarrelec , 2008; van den Heuvel and Sporns , 2011). These are oscillators 33 and 34, which correspondto the left and right posterior cingulate cortices, and oscillators 32 and 35 which represent the left andright precuneus. These central brain regions are known to be important with a higher metabolic activitythan other regions during the resting state.Importantly, we find that LRTC behaviour of this cluster, and its relationship to the other clustersin the network, is dependent on trans-callosal left-right connectivity. Indeed, disruption of the left-righttrans-callosal connections resulted in a loss of LRTCs in the rate of change of phase difference betweentime 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 phasesof a large number of other oscillators. When these oscillators try to synchronise, we suggest that those thatare well connected will be subjected to conflicting phase inputs from their neighbours and thus increasedvariation in their phase fluctuations, yielding a larger DFA exponent. These variations in fluctuation willin turn feed into the neighbouring oscillators and cause them to also have large variations in fluctuationas they attempt to synchronise with their well-connected neighbour. On the other hand, an oscillatorthat is poorly connected or connected to just one other oscillator may have a more straightforward taskof synchronising with just this (albeit changing) oscillator speed.The LRTCs in the rate of change of phase difference were also disrupted by randomisation of con-nectivity, albeit less severely than when the trans-callosal connections were severed. When a randomconnectivity 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 potentialfuture use of the Cabral model in investigating specific pathological modifications of connectivity andtheir effects on the time-varying synchronisation patterns between different brain regions. The methodhas the potential to be used to trace some types of pathological synchronisation such as may arise inepileptic or Parkinsonian conditions to any roots that they may have either in the connectivity, clusteringor noise input elements of the Cabral model and therefore potentially also of the nervous system.
Neurophysiological data
In order to show proof of principle, we have presented an example of our method’s application to neuro-physiological data, in this case EEG and EMG simultaneously recorded during voluntary muscle contrac-tion. It was through this experimental paradigm that corticomuscular coherence (CMC) in the 16 − β ) frequency range was first discovered by Conway et al. (1995);
Halliday et al. (1998) and shownto 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 ≈ . Muthukumaraswamy ,2011). As discussed earlier, the techniques introduced here allow us to focus on the fluctuations withinthe phase coupling rather than on the averaged measure of coupling. These preliminary results indicatethat the fluctuations in the rate of change of phase difference between simultaneously recorded EEG andEMG show power law scaling and LRTCs within the β frequency range. We suggest that the analysis ofinstantaneous phase diffence of neurophysiological data using the methods described in this paper willallow researchers to investigate the coupling between signals in a way that will allow a new appreciation ofthe relationship between neural synchrony and other oscillator systems approaching their critical regime. LRTCs in rate of change of phase difference and the brain
LRTCs have been associated with model dynamical systems that show efficiency in learning, memoryformation, rapid information transfer and network organisation. The broad dynamical range of whichRTCs are a marker acts to support these functions (
Linkenkaer-Hansen et al. , 2001;
Chialvo , 2010;
Sornette , 2006;
Beggs and Timme , 2012;
Stam and de Bruin , 2004;
Werner , 2010;
Linkenkaer-Hansen et al. , 2004;
Meisel et al. , 2012;
Shew et al. , 2009). It has been argued by a number ofresearchers that these properties if present would be of major benefit to the functions that human braindynamics needs to support and there is now a literature that connects the theory of critical systemswith properties of human brain dynamics (
Kitzbichler et al. , 2009;
Chialvo , 2010;
Shew et al. , 2009;
Linkenkaer-Hansen et al. , 2001;
Beggs and Plenz , 2003).In this paper, we focus on LRTCs, and because of the importance in neuroscience of brain oscillationsand the concept of communication through coherence, we make the link between LRTCs and phasesynchrony. We note that in the model systems that we have explored the highest valid DFA exponentswere recovered when the systems were close to their critical point but in a slightly more disordered statethan at exact criticality. We explained this on the basis of full synchronisation within our model systemsbeing a point at which the rate of change of phase difference is lost (observed in Ising at
T < T c and inKuramoto for increasing K ).In neurophysiological systems, it is important to appreciate that full synchronisation of neural oscil-lators is a pathological state (e.g., observed in the EEG and MEG of epileptic seizures and in EMGsshowing pathological tremor). The healthy resting brain state therefore is characterised by weak andvariable neural synchrony which would be expected to show fluctuations (temporal innovations) in ameasure of the change in phase synchrony, i.e., the rate of change of phase difference. From the perspec-tive of brain dynamics (and muscle activation dynamics) the most important constraints are to avoidpathological synchronisation whilst at the same time maintaining the potential for useful synchronisa-tion. We suggest therefore that in the healthy state the instantaneous phase difference between neuraloscillators will show power law fluctuation plots with a DFA exponent that is either 0.5 or that will showLRTCs. If LRTCs are found in the resting state then they may represent an optimum state of readinessto which the system can readily return if increased synchronisation occurs as a result of sensory stimu-lation, motor task or cognitive action. Such temporary changes in synchronisation may occur in order tosupport communication through coherence. The resting state, however, is characterised by fluctuationsof phase synchrony that have LRTCs and represent the behaviour of weakly coupled oscillators whosesynchrony can be modulated. The hypothesis that the LRTCs of rate of change of phase difference ofbrain 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 importantneed to link the criticality paradigm with the paradigm that attaches functional significance to neuralsynchrony. The methodology presented in this paper takes us some way towards this synthesis. Acknowledgement
Maria Botcharova thanks the Centre for Mathematics and Physics in the Life Sciences and ExperimentalBiology (CoMPLEX), University College London for their funding and continuing support. Simon F.Farmer was supported by University College London Hospitals Biomedical Research Centre (BRC). able 2. A list of the 66 brain regions which are represented by 66 oscillators in the Cabral model.
The abbreviations, full names and oscillator numbers corresponding to the left and the right hemispheres aregiven for each brain region. Oscillator numberAbbreviation Region Right LeftENT Entorhinal cortex 1 66PARH Parahippocampal cortex 2 65TP Temporal pole 3 64FP Frontal pole 4 63FUS Fusiform gyrus 5 62TT Transverse temporal cortex 6 61LOCC Lateral occipital cortex 7 60SP Superior parietal cortex 8 59IT Inferior temporal cortex 9 58IP Inferior parietal cortex 10 57SMAR Supramarginal gyrus 11 56BSTS Bank of the superior temporal sulcus 12 55MT Middle temporal cortex 13 54ST Superior temporal cortex 14 53PSTC Postcentral gyrus 15 52PREC precental gyrus 16 51CMF Caudal middle frontal cortex 17 50POPE Pars opercularis 18 49PTRI Pars triangularis 19 48RMF Rostral middle frontal cortex 20 47PORB Pars orbitalis 21 46LOF Lateral orbitofrontal cortex 22 45CAC Caudal anterior frontal cortex 23 44RAC Rostral anterior cingulate cortex 24 43SF Superior frontal cortex 25 42MOF Medial orbitofrontal cortex 26 41LING Lingual gyrus 27 40PCAL Pericalcarine cortex 28 39CUN Cuneus 29 38PARC Paracentral lobule 30 37ISTC Isthmus of the cingulate cortex 31 36PCUN Precuneus 32 35PC Posterior cingulate cortex 33 34The labels, brain regions and oscillator numbers used in the Cabral model. ibliography
Acebr´on, J. A., P´erez Vicente, C. J., Ritort, F., and Spigler, R. (2005), The Kuramoto model: A simpleparadigm for synchronization phenomena,
Rev. Mod. Phys. , 77, 137 – 185, doi:10.1103/RevModPhys.77.137Akaike, H. (1974), A new look at the statistical model identification,
IEEE T. Automat. Contr. , 19, 6,716–723, doi:10.1109/TAC.1974.1100705Akam, T. and Kullmann, D. M. (2010), Oscillations and filtering networks support flexible routing ofinformation,
Neuron , 67, 2, 308 – 320, doi:10.1016/j.neuron.2010.06.019Bak, P. (1996), How Nature Works: The Science of Self-organized Criticality (Copernicus (Springer)),1st editionBaker, S. N., Kilner, J. M., Pinches, E. M., and Lemon, R. N. (1999), The role of synchronyand oscillations in the motor output.,
Experimental brain research. Experimentelle Hirnforschung.Exp´erimentation c´er´ebrale , 128, 1-2, 109–117Baker, S. N., Spinks, R., Jackson, A., and Lemon, R. N. (2001), Synchronization in Monkey MotorCortex During a Precision Grip Task. I. Task-Dependent Modulation in Single-Unit Synchrony,
J.Neurophysiol. , 85, 2, 869–885Bardet, J.-M. and Kammoun, I. (2008), Asymptotic properties of the detrended fluctuation analysis oflong-range-dependent processes.,
IEEE T. Inform. Theory , 54, 5, 2041–2052Barnett, L., Lizier, J. T., Harr´e, M., Seth, A. K., and Bossomaier, T. (2013), Information flow in a kineticising model peaks in the disordered phase,
Phys. Rev. Lett. , 111, 177203, doi:10.1103/PhysRevLett.111.177203Beggs, J. M. and Plenz, D. (2003), Neuronal avalanches in neocortical circuits,
J. Neurosci. , 23, 35,11167–11177, doi:10.3410/f.1018065.217446Beggs, J. M. and Timme, N. (2012), Being critical of criticality in the brain,
Front. Fract. Physiol. , 3,163, doi:10.3389/fphys.2012.00163Beran, J. (1994), Statistics for long-memory processes, Monographs on statistics and applied probability,61 (Chapman & Hall)Berthouze, L., James, L. M., and Farmer, S. F. (2010), Human EEG shows long-range temporal corre-lations of oscillation amplitude in theta, alpha and beta bands across a wide age range,
Clin. Neuro-physiol. , 121, 8, 1187–1197, doi:10.1016/j.clinph.2010.02.163Bonilla, L. L., Neu, J. C., and Spigler, R. (1992), Nonlinear stability of incoherence and collective synchro-nization in a population of coupled oscillators,
J. Stat. Phys. , 67, 313–330, doi:10.1007/BF01049037Botcharova, M., Farmer, S. F., and Berthouze, L. (2012), Power-law distribution of phase-locking intervalsdoes not imply critical interaction,
Phys. Rev. E , 86, 051920, doi:10.1103/PhysRevE.86.051920Botcharova, M., Farmer, S. F., and Berthouze, L. (2013), A maximum likelihood based technique forvalidating detrended fluctuation analysis (ML-DFA), pre-print , arXiv:1306.5075Breakspear, 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.00190Brittain, 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 anteriorduring human treadmill walking,
J. Neurosci. Meth. , 177, 2, 334 – 347, doi:10.1016/j.jneumeth.2008.10.028Bryce, R. M. and Sprague, K. B. (2012), Revisiting detrended fluctuation analysis,
Sci Rep , 2, 315,doi:10.1038/srep00315Buzsaki, G. (2006), Rhythms of the Brain (Oxford University Press, USA), 1 editionCabral, J., Hugues, E., Sporns, O., and Deco, G. (2011), Role of local network oscillations in resting-statefunctional connectivity,
NeuroImage , 57, 1, 130–139, doi:10.1016/j.neuroimage.2011.04.010Cabral, J., Kringelbach, M. L., and Deco, G. (2012), Functional graph alterations in schizophrenia: aresult from a global anatomic decoupling?,
Pharmacopsychiatry , 45, 57–64Carreras, 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.4438Chialvo, D. R. (2010), Emergent complex neural dynamics,
Nat. Phys. , 6, 1, 744–750, doi:10.1038/nphys1803hopra, N. and Spong, M. W. (2005), On synchronization of Kuramoto oscillators, in IEEE Decis. Contr.P., volume 44 (IEEE), volume 44, 3916–3922, doi:10.1109/CDC.2005.1582773Clauset, A., Moore, C., and Newman, M. E. J. (2006), Structural Inference of Hierarchies in Networks, pre-print , doi:10.1007/978-3-540-73133-7 1, arXiv:physics/0610051Clegg, R. G. (2006), A practical guide to measuring the hurst parameter, in 21st UK PerformanceEngineering Workshop, School of Computing Science Technical Report Series, CSTR-916, Universityof Newcastle, 43–55Conway, B., Halliday, D., Farmer, S., Shahani, U., Maas, P., Weir, A., et al. (1995), Synchronizationbetween motor cortex and spinal motoneuronal pool during the performance of a maintained motortask in man.,
J. Physiol. (Lond.) , 489 ( Pt 3), 917–924Daffertshofer, A. and van Wijk, B. C. M. (2011), On the influence of amplitude on the connectivitybetween phases,
Front. Neuroinform. , 5, 6, doi:10.3389/fninf.2011.00006Daido, H. (1989), Intrinsic Fluctuation and Its Critical Scaling in a Class of Populations of Oscillatorswith Distributed Frequencies,
Progr. Theoret. Phys. , 81, 727–731, doi:10.1143/PTP.81.727Daniels, B. C. (2005), Synchronization of globally coupled nonlinear oscillators: the rich behavior of thekuramoto model,
Ohio Wesleyan Physics Dept., Essay , 7, 2, available onlineDoesburg, 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,
PLoSONE , 4, 7, e6142, doi:10.1371/journal.pone.0006142D¨orfler, F. and Bullo, F. (2011), On the Critical Coupling for Kuramoto Oscillators,
SIAM J. AppliedDynamical Systems , 10, 3, 1070–1099Farmer, S. F. (1998), Rhythmicity, synchronization and binding in human and primate motor systems,
J. Physiol. , 509, 1, 3–14, doi:10.1111/j.1469-7793.1998.003bo.xFarmer, S. F., Bremner, F. D., Halliday, D. M., Rosenberg, J. R., and Stephens, J. A. (1993), Thefrequency content of common synaptic inputs to motoneurones studied during voluntary isometriccontraction in man.,
J. Physiol. (Lond.) , 470, 127–155Fransson, P. and Marrelec, G. (2008), The precuneus/posterior cingulate cortex plays a pivotal role inthe default mode network: Evidence from a partial correlation network analysis,
NeuroImage , 42, 3,1178 – 1184, doi:http://dx.doi.org/10.1016/j.neuroimage.2008.05.059Freeman, W. and Rogers, L. (2002), Fine temporal resolution of analytic phase reveals episodic synchro-nization by state transitions in gamma EEGs.,
J Neurophysiol , 87, 2, 937–45Freeman, W. J. (2004), Origin, structure, and role of background EEG activity. Part 1. Analytic ampli-tude.,
Clin. Neurophysiol. , 115, 9, 2077–88Fries, P. (2009), Neuronal gamma-band synchronization as a fundamental process in cortical computation,
Annu. Rev. Neurosci. , 32, 1, 209–224, doi:10.1146/annurev.neuro.051508.135603Gionis, A., Mannila, H., Mielikinen, T., and Tsaparas, P. (2007), Assessing data mining results via swaprandomization.,
ACM Trans. Knowl. Discov. Data , 1, 3Granger, C. W. J. and Joyeux, R. (1980), An introduction to long-memory time series models andfractional differencing,
J. Time Ser. Anal. , 1, 1, 15–29, doi:10.1111/j.1467-9892.1980.tb00297.xHalliday, D. M., Conway, B. A., Farmer, S. F., and Rosenberg, J. R. (1998), Using electroencephalog-raphy to study functional coupling between cortical activity and electromyograms during voluntarycontractions in humans,
Neurosci. Lett. , 241, 1, 5–8Hanhij¨arvi, S., Garriga, G. C., and Puolam¨aki, K. (2009), Randomization techniques for graphs, in SDM,780–791, doi:http://dx.doi.org/10.1137/1.9781611972795.67Hardstone, 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.00450Hausdorff, J. M., Peng, C. K., Ladin, Z., Wei, J. Y., and Goldberger, A. L. (1995), Is walking a randomwalk? Evidence for long-range correlations in stride interval of human gait.,
J. Appl. Physiol. , 78, 1,349–358Honey, C. J., Sporns, O., Cammoun, L., Gigandet, X., Thiran, J. P., Meuli, R., et al. (2009), Predictinghuman resting-state functional connectivity from structural connectivity,
Proc. Natl. Acad. Sci. U SA. , 106, 6, 2035–2040, doi:10.1073/pnas.0811168106Hosking, J. R. M. (1981), Fractional differencing,
Biometrika , 68, 1, 165–176, doi:10.1093/biomet/68.1.165Hu, K., Ivanov, P. C., Chen, Z., Carpena, P., and Eugene Stanley, H. (2001), Effect of trends on detrendedfluctuation analysis,
Phys. Rev. E , 64, 011114, doi:10.1103/PhysRevE.64.011114Hurvich, C. M. and Tsai, C.-L. (1989), Regression and time series model selection in small samples,
Biometrika , 76, 2, 297–307, doi:10.1093/biomet/76.2.297sing, E. (1925), Beitrag zur Theorie des Ferromagnetismus,
Zeitschrift f¨ur Physik A Hadrons and Nuclei ,31, 1, 253–258, doi:10.1007/BF02980577James, L. M., Halliday, D. M., Stephens, J. A., and Farmer, S. F. (2008), On the development of humancorticospinal oscillations: age-related changes in EEG-EMG coherence and cumulant,
Eur. J. Neurosci. ,27, 12, 3369–3379, doi:10.1111/j.1460-9568.2008.06277.xKarmeshu and Krishnamachari, A. (2004), Sequence variability and long-range dependence in DNA: Aninformation theoretic perspective,
Lecture Notes in Computer Science , 3316, 1354–1361Kinouchi, O. and Copelli, M. (2006), Optimal dynamical range of excitable networks at criticality,
Nat.Phys. , 2, 5, 348–351, doi:10.1038/nphys289Kitzbichler, M. G., Smith, M. L., Christensen, S. R., and Bullmore, E. (2009), Broadband criticality ofhuman brain network synchronization,
PLoS Comput. Biol. , 5, 3, e1000314, doi:10.1371/journal.pcbi.1000314Kuramoto, Y. (1975), Self-entrainment of a population of coupled nonlinear oscillators, in H. Araki, ed.,International Symposium on Mathematical Problems in Theoretical Physics, Lecture Notes in Physics,Vol. 39, (Springer, New York, NY, USA), 420–422Kuramoto, Y. (1984), Chemical Oscillations, Waves, and Turbulence (Springer–Verlag, New York)Linkenkaer-Hansen, K., Nikouline, V. V., Palva, J. M., and Ilmoniemi, R. J. (2001), Long-range temporalcorrelations and scaling behavior in human brain oscillations.,
J. Neurosci. , 21, 4, 1370–1377Linkenkaer-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, 1, 203–218, doi:10.1111/j.1460-9568.2004.03116.xMackay, D. J. C. (2003), Information Theory, Inference and Learning Algorithms (Cambridge UniversityPress), first editionMeisel, C., Storch, A., Hallmeyer-Elgner, S., Bullmore, E., and Gross, T. (2012), Failure of adaptiveself-organized criticality during epileptic seizure attacks,
PLoS Comput. Biol. , 8, 1, e1002312, doi:10.1371/journal.pcbi.1002312Mertens, D. C. (2011), Population-specific predictions for the finite Kuramoto model and collectivesynchronization in a system with resonant coupling, Ph.D. thesis, University of Illinois at Urbana-ChampaignMetropolis, N., Rosenbluth, A., Rosenbluth, M., Teller, A., and Teller, E. (1953), Equation of statecalculations by fast computing machines,
J. Chem. Phys. , 21, 1087Miritello, G., Pluchino, A., and Rapisarda, A. (2009), Phase transitions and chaos in long-range modelsof coupled oscillators,
Europhys. Lett. , 85, 1, 10007Moretti, P. and Mu˜noz, M. A. (2013), Griffiths phases and the stretching of criticality in brain networks.,
Nat. Commun. , 4, 2521Muthukumaraswamy, S. D. (2011), Temporal dynamics of primary motor cortex γ oscillation amplitudeand piper corticomuscular coherence changes during motor control., Exp. Brain. Res. , 212, 4, 623–633,doi:10.1007/s00221-011-2775-zOnsager, L. (1944), Crystal statistics. I. A two-dimensional model with an order-disorder transition,
Phys. Rev. , 65, 3-4, 117–149Peng, 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, 2, 1685–1689Peng, 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-4Peng, C. K., Havlin, S., Stanley, H. E., and Goldberger, A. L. (1995b), Quantification of scaling exponentsand crossover phenomena in nonstationary heartbeat time series,
Chaos , 5, 1, 82–87Pfurtscheller, G. (1977), Graphical display and statistical evaluation of event-related desynchronization,
Electroen. Clin. Neuro. , 43, 757 – 760Pfurtscheller, G. (1992), Event-related synchronization (ERS): an electrophysiological correlate of corticalareas at rest,
Electroen. Clin. Neuro. , 83, 1, 62 – 69, doi:10.1016/0013-4694(92)90133-3Pikovsky, A., Rosenblum, M., and Kurths, J. (2003), Synchronization: A Universal Concept in NonlinearScience, Cambridge Nonlinear Science Series (Cambridge University Press), doi:10.1119/1.1475332Plenz, D. and Chialvo, D. R. (2009), Scaling properties of neuronal avalanches are consistent with criticaldynamics, pre-print , arXiv:0912.5369Poil, S.-S., Hardstone, R., Mansvelder, H. D., and Linkenkaer-Hansen, K. (2012), Critical-State Dynam-ics of Avalanches and Oscillations Jointly Emerge from Balanced Excitation/Inhibition in NeuronalNetworks,
J. Neurosci. , 32, 29, 9817–9823, doi:10.1523/jneurosci.5990-11.2012Priemer, R. (1990), Introductory Signal Processing (World Scientific Publishing, River Edge, NJ, USA)riesemann, V., Munk, M., and Wibral, M. (2009), Subsampling effects in neuronal avalanche distribu-tions recorded in vivo,
BMC Neuroscience , 2009, 1 – 20, doi:10.1186/1471-2202-10-40Robinson, P. M., ed. (2003), Time series with long memory (Oxford University Press)Samorodnitsky, G. (2006), Long range dependence,
Found Trends Stoch Syst , 1, 3, 163–257Schoffelen, J.-M., Oostenveld, R., and Fries, P. (2005), Neuronal Coherence as a Mechanism of EffectiveCorticospinal Interaction,
Science , 308, 5718, 111–113, doi:10.1126/science.1107027Shew, W. L. and Plenz, D. (2013), The functional benefits of criticality in the cortex,
The Neuroscientist ,19, 1, 88–100, doi:10.1177/1073858412445487Shew, W. L., Yang, H., Petermann, T., Roy, R., and Plenz, D. (2009), Neuronal avalanches implymaximum dynamic range in cortical networks at criticality,
J. Neurosci. , 29, 49, 15595–15600, doi:10.1523/jneurosci.3864-09.2009Shriki, O., Alstott, J., Carver, F., Holroyd, T., Henson, R. N. A., Smith, M. L., et al. (2013), Neuronalavalanches in the resting MEG of the human brain.,
J. Neurosci. , 33, 16, 7079–7090Siegel, M., Donner, T. H., and Engel, A. K. (2012), Spectral fingerprints of large-scale neuronal interac-tions,
Nat. Rev. Neurosci. , 13, 2, 121–134, doi:10.1038/nrn3137Singer, W. (1999), Neuronal synchrony: a versatile code for the definition of relations?,
Neuron , 24, 1,49–65Sornette, D. (2006), Critical Phenomena in Natural Sciences: Chaos, Fractals, Self-Organization andDisorder: Concepts and Tools (Springer), 2nd editionStam, C. J. and de Bruin, E. A. (2004), Scale-free dynamics of global functional connectivity in thehuman brain,
Hum. Brain. Mapp. , 22, 2, 97–109, doi:10.1002/hbm.20016Stanley, 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, 1, 214–253Strogatz, S. H. and Mirollo, R. E. (1991), Stability of incoherence in a population of coupled oscillators,
J. Stat. Phys. , 63, 3, 613–635, doi:10.1007/BF01029202Taqqu, M. S., Teverovsky, V., and Willinger, W. (1995), Estimators for long-range dependence: Anempirical study,
Fractals , 3, 785–798van den Heuvel, M. P. and Sporns, O. (2011), Rich-club organization of the human connectome.,
J.Neurosci. , 31, 44, 15775–15786, doi:10.1523/JNEUROSCI.3539-11.2011Wang, F., Yamasaki, K., Havlin, S., and Stanley, H. E. (2005), Scaling and memory of intraday volatilityreturn intervals in stock market, pre-print , arXiv:physics/0511101Wang, X.-J. J. (2010), Neurophysiological and computational principles of cortical rhythms in cognition.,
Physiological reviews , 90, 3, 1195–1268, doi:10.1152/physrev.00035.2008Werner, G. (2009), Fractals in the nervous system: conceptual implications for theoretical neuroscience, pre-print , arxiv:0910.2741Werner, G. (2010), Fractals in the nervous system: conceptual implications for theoretical neuroscience.,
Front. Physiol. , 1, 1–15, doi:10.3389/fphys.2010.00015Willinger, W., Taqqu, M. S., and Teverovsky, V. (1999), Stock market prices and long-range dependence,