Detecting dynamic spatial correlation patterns with generalized wavelet coherence and non-stationary surrogate data
DDetecting dynamic spatial correlation patternswith generalized wavelet coherence and non-stationary surrogate data
M. Chavez and B. Cazelles
2, 3 CNRS UMR-7225, Hˆopital de la Piti´e-Salpˆetri`ere. 75013 Paris, France IRD-UPMC UMI-209, UMMISCO, 93143 Bondy, France CNRS UMR-8197, IBENS, Ecole Normale Sup´erieure. 75005 Paris, France
Time series measured from real-world systems are generally noisy, complex and display statisticalproperties that evolve continuously over time. Here, we present a method that combines waveletanalysis and non-stationary surrogates to detect short-lived spatial coherent patterns from multivari-ate time-series. In contrast with standard methods, the surrogate data used here are realisations ofa non-stationary stochastic process, preserving both the amplitude and time-frequency distributionsof original data. We evaluate this framework on synthetic and real-world time series, and we showthat it can provide useful insights into the time-resolved structure of spatially extended systems.
PACS numbers: 02.50.Sk, 05.45.Tp, 05.45.Xt, 89.75.Fb
Synchronization is a fundamental phenomenon de-scribed in many biological and physical contexts forwhich there are two or more interacting oscillatory sys-tems [1]. The interactions between coupled oscillatorsin real systems continuously create and destroy synchro-nised states, which can be observed as noisy and tran-sient coherent patterns. The statistical detection of spa-tial synchrony in networks of coupled dynamical systemsis therefore of great interest in disciplines such as geo-physics, physiology and ecology [2, 3].Statistical significance of transient coherent patternscannot be assessed by classical spectral measures andtests, which require signals to be stationary [3, 4]. Syn-chrony estimators based on nonparametric methods havethe advantage of not requiring any assumption on thetime-scale structure of the observed signals. Amongthem, measures of synchrony or coherence based onwavelet transforms have been widely used to detect inter-actions between oscillatory components in different realsystems, i. e. neural oscillations, business cycles, climatevariations or epidemics dynamics [2, 3].In recent years, different significance tests for thewavelet cross-spectrum or wavelet coherence have beendeveloped to detect oscillatory patterns with covaryingdynamics [3–5]. Unfortunately, the statistical assump-tions of these tests are not always compatible with thestructure of the data considered, and significance levelsoften depend on the structure of the wavelets applied [5].A rigorous theoretical framework cannot therefore be de-rived, and Monte Carlo simulations have to be performedto estimate the significance level [5].Surrogate data techniques have been proposed as non-parametric resampling methods for testing general hy-potheses on data without making assumptions on theunderlying generating process [6]. However, time seriesmeasurements from real systems generally display irregu-lar fluctuations, long-term trends, or a time-varying spec-tra. Such properties are incompatible with the main as- sumptions of standard surrogate data based on Fouriertransform [7, 8].Recently, parametric models have been also appliedto test wider classes of null hypothesis, including non-stationary behaviour [8]. Some limitations of these ap-proaches include the relatively large basis dimensionneeded to obtain good optimisation, and the monitor-ing needed to control the instabilities in the estimatedmodel [9]. Recent studies have proposed the use of dis-crete wavelet transforms (DWT) for resampling time se-ries such that the multiscale structure of original data ispreserved [10, 11]. The main advantage of DWT is theirability to concentrate the signal’s variance in a limitednumber of coefficients. Nevertheless, the number of datapoints heavily influences this decomposition (the numberof scales); which may render the scale decompositions dif-ficult to interpret [3]. Although continuous wavelets oftenyield a redundant decomposition across scales, they aremore robust to noise as compared with other decompo-sition schemes [3, 4, 12].In this work, we use a continuous wavelet-basedapproach to detect spatial coherent patterns in non-stationary multivariate observations. We generalise thewavelet coherence to multivariate time series and we ex-tend the classic phase-randomised surrogate data algo-rithm to the time-frequency domain for generating non-stationary surrogates. This procedure preserves both theoriginal amplitude and time-frequency energy (spectro-gram) distributions. Compared with other surrogate al-gorithms, our method better replicate the time-frequencystructure of real data. These non- stationary surrogatesare used to assess the significance of transient coherentpatterns found in multivariate time series. We evaluatethe proposed method in different synthetic and real-worldnon-stationary data, and we show that this approach cansubstantially improve the detection of time-varying spa-tial coherence.We start by considering the time-frequency (TF) dis-tributions obtained by correlating a time series x ( t ) a r X i v : . [ phy s i c s . d a t a - a n ] A p r f r e q u e n c y ( H z ) B) Stationary surrogate data C) Non-stationary surrogate dataA) Original EEG data f r e q u e n c y ( H z ) f r e q u e n c y ( H z ) D) Amplitudes distributions
Non-stationary surrogatedataOriginal dataStationary surrogate data -4 -3 -2 -1 0 1 2 3 40.400.30.20.10.5
FIG. 1. Exemple of the time-frequency (TF) structure for different surrogate algorithms applied to epileptic EEG data. A)The original time series after reconstruction by the wavelet filtering, B) surrogate data generated with the iterative AmplitudeAdjusted Fourier Transform (iAAFT) algorithm, C) surrogate generated with our algorithm and D) distributions of amplitudesfor the different algorithms. The color maps code for | ( W x ( t, f )) | values. Wavelet analysis of EEG were done over the frequencyrange 0 −
128 Hz. For better visualisation spectra are displayed only for f <
30 Hz. Refer to Supplementary Material for acomparison with other algorithm (based on DWT) to replicate the TF structure. with a scaled and translated version of a chosen motherwavelet w s,τ ( t ) = √ s w ( t − τs ). In the following, we alwaysconsider the mother complex Morlet wavelet defined as w ( θ ) = π − / exp( − θ / × exp( iω θ ), where the param-eter ω controls the time scale resolution of each wavelet.To quantify the relationships between two non-stationary signals, x i ( t ) and x j ( t ), the wavelet cross-spectrum is given by W i,j ( t, f ) = W i ( t, f ) W ∗ j ( t, f ),where ∗ denotes the complex conjugate operatorand W k ( t, f ) is the wavelet transform of signal x k ( t ). Let us now consider M zero-mean time series x ( t ) , . . . , x M ( t ), and define the complex coherence spec-trum as C i,j ( t, f ) = (cid:10) W i,j ( t,f ) (cid:11) (cid:13)(cid:13)(cid:13) (cid:10) W i,i ( t,f ) (cid:11) (cid:13)(cid:13)(cid:13) (cid:13)(cid:13)(cid:13) (cid:10) W j,j ( t,f ) (cid:11) (cid:13)(cid:13)(cid:13) for i, j =1 , . . . , M , where (cid:104)·(cid:105) denotes a smoothing operator bothin time and frequency [13].In bivariate data analysis, the wavelet coherence is de-fined as Γ i,j ( t, f ) = | C i,j ( t, f ) | . To extend this idea tothe general case of M (cid:62) Σ ( t, f ) at every point in the time-frequency domaincontaining all the pairwise coherence spectra [14]: Σ ( t, f ) = C , ( t, f ) . . . C ,M ( t, f ) C , ( t, f ) 1 . . . C ,M ( t, f )... ... . . . ... C M, ( t, f ) C M, ( t, f ) . . . , (1)The time-varying spatial coherence (TVSC) can be de-fined by Ψ( t, f ) = 1 M − (cid:0) λ Σmax ( t, f ) − (cid:1) , (2)where λ Σmax ( t, f ) denotes the largest eigenvalue of thespectral matrix Σ ( t, f ).The values of Ψ( t, f ) are bounded between 0 (cid:54) Ψ( t, f ) (cid:54)
1, reaching the maximum when all the M signals are locally -in the time-frequency plane- pair-wise correlated ( Σ ( t, f ) becomes an all-ones matrix with λ Σmax ( t, f ) = M ); and the minimum when all signals are completely uncorrelated ( Σ ( t, f ) = I and λ Σmax ( t, f ) =1) .Interestingly, for the case M = 2, Σ ( t, f ) is given bythe matrix (cid:20) C , ( t, f ) C , ( t, f ) 1 (cid:21) , whose largest eigen-value is λ Σmax ( t, f ) = 1 + | C , ( t, f ) | , which yieldsΨ( t, f ) = ( λ Σmax ( t, f ) −
1) = | C , ( t, f ) | . In the bivari-ate case, this therefore reduces the TVSC to the classicdefinition of the wavelet coherence Ψ ( t, f ) = Γ ( t, f ).In wavelet-based analysis, test statistics are stronglyaffected by data’s structure, the mother wavelet’s prop-erties, and by the smoothing applied [5, 13]. In this work,the statistical properties of Ψ( t, f ) under the null hypoth-esis H of M uncorrelated processes are determined byMonte Carlo simulation. To do this, we generate a num-ber of surrogate data realisations (cid:98) x j ( t ) , j = 1 , . . . , K byrepeating the randomisation procedure K times. Thestatistical significance of Ψ( t, f ) values was assessed by az-test to quantify the statistical deviation from values ob-tained in the ensemble of surrogate data. To correct formultiple testing, the false discovery rate (FDR) methodwas applied [16]. With this approach, the threshold ofsignificance was set such that the expected fraction offalse positives over the time-frequency plane is restrictedto q (cid:54) . (cid:98) x ( t ) can be obtained by ran-domising the phase structure of the original signal x ( t )in the time-frequency domain. As the Morlet waveletis a complex function, the resulting wavelet decom-position W x ( t, f ) has both real and imaginary parts.We can therefore write W x ( t, f ) in terms of its phase φ x ( t, f ) = tan − (cid:61) ( W x ( t,f )) (cid:60) ( W x ( t,f )) and modulus | ( W x ( t, f )) | .The wavelet-based surrogate algorithm first generatesa Gaussian white noise time series to match the origi-nal data length and it then derives the wavelet trans-form to extract the phase φ noise ( t, f ). We use this ran- In case of stationary observations, eigenvalues of the covariancematrix are commonly used in radio communications for detectingspatial correlations between multivariate time series [14]
B) Stationary surrogate data C) Non-stationary surrogate dataA) Original measles data
20 200 400 600 800 1000 1200 1400 1600 1800 time (weeks) sc a l e ( w ee k s ) sc a l e ( w ee k s ) sc a l e ( w ee k s ) D) Amplitudes distributions
Non-stationary surrogate dataOriginal dataStationary surrogate data vaccine era
FIG. 2. Exemple of the TF structure for different surrogate algorithms applied to the squared root transformed measles data.Same stipulations as in the caption of Fig. 1. Gray box in upper plot A indicates the vaccine era. Black transparent maps oftime-frequency plots indicate the cone of influence that delimits the regions not influenced by edge effects [15]. domised phase and the WT modulus of the original sig-nal to obtain a surrogate time-frequency distribution W (cid:98) x ( t, f ) = | ( W x ( t, f )) | exp( iφ noise ( t, f )). A surrogatetime series can be reconstructed by taking the real partof the inverse wavelet transform. Finally, the surrogate (cid:98) x ( t ) is rescaled to the distribution of the original databy sorting the data (after the wavelet filtering in the fre-quency band of interest) according to the ranking of thewavelet-based surrogate [6]. As its Fourier-based coun-terpart, our scheme can be iteratively repeated to betteradjust the time-frequency distribution of the surrogate.Other methods based on DWT have been proposedto resample non-stationary time series. Compared tostandard surrogate data, they preserve better the TFstructure of original data. Nevertheless, our algorithmreplicates the TF distributions more accurately (see theSupplementary Material for additional details and com-parisons).To illustrate our surrogate data method, we consideran electroencephalographic (EEG) recording from a pe-diatric subject with intractable epileptic seizures . Al-though our approach is applicable to any neuroimagingfunctional method (e.g. EEG, fMRI, and MEG signals)here we use the EEG as this modality of acquisitionhas the major feature that collective neural behaviors,i.e., synchronization of cortical assemblies are reflectedas time-varying interactions between EEG signals. Thefile studied here contains 22 EEG signals sampled at 256Hz according to the 10-20 bipolar montage [17].The non-stationarity of epileptic EEG signals is clearlyillustrated in Fig. 1-A. One can notice that the frequencycontent of epileptic oscillations may change rapidly acrosstime over a range of frequencies. The time-frequency plotexhibits a short fast oscillatory behavior ( f ≈
15 Hz)around t = 3 s followed by slow and large oscillationsaccompanying the epileptic seizure after t = 6 s. Asdepicted in Fig. 1-B, classical stationary surrogate data(here we used the iAAFT algorithm [6]) is not able toreplicate the non-stationary oscillations embedded in the The EEG data was obtained from the open repository CHB-MIT Scalp EEG Database ( ) original signal. In contrast, the time-varying spectrumof original signal is clearly conserved with our algorithm,as illustrated in Fig. 1-C. Plot in Fig. 1-D confirms thatthe three surrogate algorithms conserve the amplitudedistributions.Another paradigmatic example of non-stationary spa-tial synchrony is that observed in population dynamics.Here, we considered the weekly measles notifications inseven large English cities studied in Refs. [18]. Measlesepidemics generally exhibit a non-stationary dynamicswith a regular and highly epidemics before nationwidevaccination programs, and an irregular and spatially un-correlated dynamics in the vaccine era. As illustrated inFig. 2, the data display multiannual cycles that dramat-ically varies with time, specially after vaccination. Thisrich behavior can not be encompassed by classical sta-tionary surrogate data (Fig. 2-B). Instead, the wavelet-based method perfectly keeps the variations of epidemicperiods observed in the original time series (Fig. 2-C).Throughout this work, the number of scales of thewavelet decomposition was selected such that an accu-rate reconstruction of the original signal was obtained,and the non-stationary oscillations and transient eventsobserved in the time series were accurately captured.We then test now the performance of our frameworkto detect spatial coherent components on two syntheticdatasets with time-varying structure. In the first bench-mark, the spatial system consists of 5 linear oscillatorsdescribed by the following autoregressive (AR) model: x t = 0 . √ x t − − . x t − + (cid:15) t ,x t = 0 . x t − − . x t − + k x t − + (cid:15) t x t = 0 . x t − − . x t − + 0 . x t − + (cid:15) t ,x t = k x t − + 0 . √ x t − + 0 . √ x t − + (cid:15) t ,x t = − . √ x t − + 0 . √ x t − + (cid:15) t (3)where t denotes a discrete time index, (cid:15) i are independentwhite noise processes with zero means and unit variances,and k i are coupling strengths. Here, we set k = 0 and k = 0 .
15 for t < k = − . k = 0 . t (cid:62) t, f ) still provides a qualitativedescription in case of nonlinear oscillators. Indeed, weconsider a network of i = 1 , . . . ,
10 coupled non-identicalchaotic R¨ossler oscillators. The equations of motion read˙ x i = − ω i y i − z i + λ (cid:104)(cid:80) j ξij ( x j − x i ) (cid:105) + σ i η i , ˙ y i = ω i x i + 0 . y i , ˙ z i = 0 . z i ( x i −
10) (4)where λ ( t ) is the time-varying coupling strength, ξij arethe elements of the coupling matrix (a random graph withan average number of links per node k m = 4); ω i is thenatural frequency of the i th oscillator (randomly assignedfrom a uniform distribution with values between 0 . (cid:54) ω i (cid:54) . η i denotes a Gaussian delta correlated noisewith (cid:104) η i ( t ) (cid:105) = 0 and (cid:104) η i ( t ) η i ( t (cid:48) ) (cid:105) = 2 Dδ ( t − t (cid:48) ), D = 0 . λ varies with time as follows: λ = 0 . < t <
900 and λ = 0 .
001 elsewhere. f r e q u e n c y ( H z ) f r e q u e n c y ( H z ) B) Stationary surrogate data testC) Non-stationary surrogate data testA) Original data f r e q u e n c y ( H z ) f r e q u e n c y ( H z ) AR Rossler
FIG. 3. Ψ( t, f ) values estimated from synthetic time seriesand statistical differences with those values obtained from dif-ferent surrogate data. A) The original time series (gray boxesdelimitate the region where the system is synchronized), B)surrogate data test based on the iAAFT algorithm and C)with our algorithm. The color maps code for Ψ( t, f ) values.Unmasked color regions in panels B-C indicate the significantlevels. Black transparent maps indicate the cones of influ-ence [15].
Here, all time series are first centered and set to havezero mean and unit variance. Then, the TVSC values arecomputed. To estimate the distribution of Ψ( t, f ) under H we generate 100 surrogates for each time series.In Figs. 3-(A-C) we report, respectively, the time se-ries generated by the above models, the significant co-herent components detected by Ψ( t, f ) in combinationwith classical surrogate data and with the non-stationarysurrogate data. Results reveal that stationary random-izations detect several large spurious synchrony patcheson the time-frequency plane, e.g. the large patches be-fore t = 1000 for the coupled AR model, or those out of the synchronous region for the coupled R¨ossler model(500 < t < t, f )with non-stationary surrogate data, constitutes a goodcriterion to assess spatial coherence in the case of non-linear dynamical time series.To illustrate our approach on real-world time series, westudy the two spatial systems described above. The sit-uation with EEG data is illustrated in left panels of Fig-ure 4. The first crucial observation is that, as expected inepilepsy dynamics, spatial coherent patterns are not timeinvariant, but instead they exhibit a rich time-frequencystructure during seizure evolution. Results clearly showthat classical surrogate data test may yield to the detec-tion of large synchronous regions, specially at high fre-quency bands ( f (cid:62)
20) Hz. In contrast, non-stationarysurrogates improves the time-frequency localization ofspatial correlation patterns. A first synchronous patternseem to involve the low-amplitude fast oscillations oftenobserved during the first seconds of epileptic seizures. In-terestingly, the absence of significative values of Ψ( t, f )between t = 4 − t = 8 s. This fully agrees withprevious findings suggesting a neural desynchronizationbefore the propagation of seizures which could facilitatethe development of local pathological recruitments [19].Right panels of Figure 4 show the results for themeasles data. We observe from Ψ( t, f ) values that globalinteractions between major epidemics change relativelysmoothly through time. Classic surrogate analysis cancapture epidemic’s dynamics at different scales, but doesnot allow a proper description when they change withtime. Indeed, standard surrogate data test reveals nosignificant spatial correlation patterns. Conversely, ourapproach clearly detects the main changes in spatial cor-relation structure: a high spatial coherence between themajor epidemic (mainly biennial) component of time se-ries is clearly identified in the pre-vaccine era. The inter-actions between the smaller epidemics with longer peri-ods observed after vaccination are not found to be statis-tically significant. This is a remarkable result as it sup-ports previous findings that during the pre-vaccinationera, measles dynamics is characterized by a high spatialcorrelation of biennal epidemic patterns, while the vacci-nation eliminates large epidemics yielding thus a signifi-cant spatial decorrelation [18].In conclusion, we have addressed a fundamental prob-lem in complex systems: detecting, from scalar observa-tions, the time scales involved in spatial interactions of f r e q u e n c y ( H z ) f r e q u e n c y ( H z ) B) Stationary surrogate data testC) Non-stationary surrogate data testA) Original data sc a l e ( w ee k s ) sc a l e ( w ee k s ) EEG Measles te FIG. 4. Ψ( t, f ) values estimated from real spatial systems andstatistical differences with those values obtained from differ-ent surrogate data. For measles data, missing values in eachoriginal time series were imputed using a local average, i.e.the mean of the two neighboring time points. Same stipula-tions as in the caption of Figure 3. See Supplementary Mate-rial for a comparison with other algorithm (based on DWT)to detect these coherent patterns. oscillators with time-varying spectral components. Clas-sical surrogate data tests require time-series to be sta-tionary. Nevertheless, data recorded from real-world sys-tems are generally noisy and non-stationary. In orderto study their interactions we propose a complementaryapproach based on wavelet analysis. Wavelet coherenceis generalized as a method for detecting transient butsignificant coherence between multivariate nonlinear sig-nals. The classic surrogate algorithm is also generalizedto produce non-stationary surrogates. Several artificialand real non-stationary, linear and nonlinear time seriesare examined, in order to demonstrate the advantages ofour approach [20].Other wavelet-based methods have been used to anal-yse the relationships between multivariate signals [1, 21].Nevertheless, standard significance tests assume station-arity of observations, which strongly affects the signifi-cance of the detected coherent patterns. Our results pro-vide evidence of the constructive role of non-stationarysurrogate data to uncover changes of correlation patternsin multivariate time series. Results confirm that ourmethod performs better than stationary Fourier-basedand non-stationary DWT-based surrogate algorithms.The detection of spatial correlations in other multivariatedata (e.g. financial or climate time series) might providemeaningful insights into the structure of other spatiallyextended systems.BC is partially supported by the French Agence Na-tionale de la Recherche with the PANIC project (ANR-14-CE02-0015-01) [1] A. Pikovsky, M. Rosemblum, and J. Kurths,
Synchroniza-tion. A universal concept in nonlinear systems , CambridgeNonlinear Science Series 12, (Cambridge University Press,UK, 2001).[2] M. Le Van Quyen, et al. J. Neurosci. Methods , 83(2001); B. T. Grenfell, O. N. Bjørnstad, and J. Kappey,Nature , 716 (2001).[3] B. Cazelles, et al. Oecologia , 287 (2008); P. Clemson,G. Lancaster, A. Stefanovska, Proceedings of the IEEE, , 223 (2016).[4] C. Torrence, G. P. Compo, Bull Am Meteorol Soc , 61(1998).[5] D. Maraun, and J. Kurths, Nonlin. Processes Geophys. , 505 (2004); D. Maraun, J. Kurths, and M. Holschnei-der, Phys. Rev. E , 016707 (2007); L. W. Sheppard, A.Stefanovska, and P. V. E. McClintock, Phys. Rev. E, ,046205 (2012); B. Cazelles, K. Cazelles, and M. Chavez,J. R. Soc. Interface , 20130585 (2014)[6] T. Schreiber, and A. Schmitz, Physica D , 346–382(2000); V. Venema, F. Ament, and C. Simmer, Nonlin.Processes Geophys. , 321–328 (2006).[7] T. Nakamura, and M. Small, Phys. Rev. E. , 056216(2005); T. Rouyer, et al. Mar. Ecol. Prog. Ser. , 11–23(2008); J. H. Lucio, R. Vald´es, and L. R. Rodr´ıguez, Phys.Rev. E. , 056202 (2012).[8] D. Kugiumtzis, Phys. Rev. E. , R25 (2000); L. Faes, H.Zhao, K. H. Chon, and G. Nollo, IEEE Trans. Biomed.Eng. , 685 (2009).[9] R. Zou, H. Wang, and K. H. Chon, Ann. Biomed. Eng. , 840 (2003); X. Luo, Master’s Thesis, University ofTennessee, (2005).[10] M. Breakspear, M. Brammer, and P. A. Robinson, Phys-ica D , 1 (2003); C. J. Keylock, Physica D , 219(2007);[11] M. Paluˇs, Phys. Rev. Lett. , 134101 (2008); C. J.Keylock, Phys. Rev. E , 032123 (2017).[12] S. Mallat, A wavelet tour of signal processing (AcademicPress, San Diego, 1998).[13] E. A. Cohen, and A. T. Walden, IEEE Trans. Sig. Pro-cess. , 2964 (2010).[14] R. C. Qiu, Z. Hu, H. Li, and M. C. Wicks, Cognitive radiocommunication and networking: Principles and practice (John Wiley & Sons, UK, 2012)[15] As the wavelet is centered close to the edges of time se-ries, edge effects occur. The area of the TF plane wheresuch effects are relevant, the so-called cone of influence,was chosen as the e -folding time of the Morlet waveletfunction [4].[16] Y. Benjamini and Y. Hochberg, J. R. Stat. Soc. Ser. B(Methodol.) , 289 (1995) .[17] A. Shoeb, PhD Thesis, Massachusetts Institute of Tech-nology, 2009; A. L. Goldberger et al., Circulation ,e215 (2000).[18] B. M. Bolker, and B. T. Grenfell, Proc. Natl. Acad. Sci.USA, , 12648 (1996); P. Rohani, D. J. D. Earn, and B.T. Grendfell, Science, , 968 (1999).[19] F. Amor, et al. Neuroimage , 950 (2009); M. Chavez,et al. IEEE Trans. Biomed. Eng. , 571 (2003).[20] Matlab and R codes implementing the non-stationarysurrogate algorithm are available at the open access reposi-tory https://zenodo.org (uploads 1213657 and 1213696) [21] M. Dhamala, G. Rangarajan, and M. Ding, Phys. Rev.Lett. , 018701 (2008); E. K. W. Ng, and J. C. L. Chan,J. Atmos. Oceanic Technol. , 1845 (2012); L. W. Shep- pard, J. R. Bell, R. Harrington, and D. C. Reuman, Nat.Clim. Change , 610 (2016). Supplementary Material
Surrogate data algorithms based on discretewavelet transforms.
Recent studies have propose theuse of discrete wavelet transforms (DWT) for resamplingtime series such that the multiscale structure of originaldata is preserved. To this end, different randomizationschemes have been used to resample the wavelet coef-ficients: random permutations, block resampling, or aniAAFT algorithm have been applied to wavelet coeffi-cients to preserve the local mean and variance at eachscale [S1–S3]. More refined permutations schemes havebeen proposed to preserve the interactions among scalesreproducing thus the multifractal properties of the orig-inal data [S4, S5].In the surrogate algorithm proposed by Ref. [S4]wavelet coefficients are resampled (permuted) on a dyadictree that, when inverted, yields a data set with mul-tifractal characteristics as the original data. Althoughan amplitude adjustment is used to recover the origi-nal distribution of wavelet coefficients at each scale ofthe wavelet decomposition, the histogram of the newtime series does not exactly replicate the histogram ofthe original data. A more refined method (the itera-tive Amplitude Adjusted Wavelet Transform, iAAWT)is proposed in Ref. [S5] in which the wavelet decomposi-tion is obtained by means of a dual-tree complex DWT.At each scale of the decomposition, phases are extractedand randomized. These phases are then combined withthe original amplitudes and the inverse wavelet transformis applied to obtain a new time series. An amplitudeadjustment is applied to preserve the histogram of theoriginal data. Last steps are iterated until convergence.Compared with other methods based on DWT, this algo-rithm has been shown to be more accurate and precise toapproximate the multiscale structure of time series witha time-varying spectrum.The main advantage of all the methods based on aDWT is their ability to concentrate the signal’s variancein a limited number of scales. Nevertheless, as these de-compositions are generally computed on a dyadic base(the number of scales is always 2 n ), the number of datapoints heavily influences these decompositions and par-ticularly the time-scale localization of the signal’s vari-ance; which may render the scale decompositions difficultto interpret, especially their time evolution [S6]. Replication of the time-frequency structure withother algorithm.
To compare the ability of our methodto replicate the time-frequency (TF) distribution of the original data, we estimate the error: e = (cid:107)| W (cid:98) x ( t, f ) | − | W x ( t, f ) |(cid:107) F (cid:107)| W x ( t, f ) |(cid:107) F (S1)were (cid:107) A (cid:107) F = (cid:112) Tr( AA T ) denotes the Frobenius norm ofreal matrix A , and W (cid:98) x ( t, f ) and W x ( t, f ) are the wavelettransforms of surrogate and original data, respectively. r e l a t i v e e rr o r A) EEG data
Convergence10 iterations1 iterationDWT0 r e l a t i v e e rr o r B) Measles data
FIG. S1. Comparsion to preserve the TF distributions of theEEG and measles data. Violin plots represent distributionsof error values over 100 surrogates (gray points). The error ofiAAWT algorithm is significantly larger than those obtainedby our method at different number of iterations (a right-tailedWilcoxon rank-sum test with alpha set at p = 0 . To compare our method, different surrogating algo-rithms are applied to the two time series studied in Figs. 1and 2 of the main text (the EEG and measles time se-ries). We firstly compare our algorithm with the iAAWTmethod (refer to Ref. [S5] for detailed description of thealgorithm) S1 . We also apply an improved iterative ver-sion of our method, in which an approximation of theoriginal TF distribution is refined with each iteration un-til it is determined to be sufficiently accurate, at whichtime the procedure terminates. S1 We have used the code available at the File Exchange site ofMathWorks here
Non-stationary (CWT) surrogate data testsA) Original data f r e q u e n c y ( H z ) f r e q u e n c y ( H z ) f r e q u e n c y ( H z ) F) 500 surrogates f r e q u e n c y ( H z ) f r e q u e n c y ( H z ) F) 100 surrogatesF) 10 surrogates B) Stationary surrogate data test C) Non-stationary (DWT) surrogate data(100 surrogates)(100 surrogates)
FIG. S2. Ψ( t, f ) values estimated from epileptic EEG data and statistical significances obtained from different surrogate datatests. A) The original time series; B) surrogate data test of the time-varying spatial coherence with the iterative AmplitudeAdjusted Fourier Transform (iAAFT) algorithm; C) non-stationary surrogate data test with the iAAWT algorithm; andsignificance tests D) N = 10, E) N = 100, and F) N = 500 surrogate data obtained with the continuous wavelet transform(CWT). Color maps code for Ψ( t, f ) values. Unmasked color regions in panels B-F indicate the significant levels. Blacktransparent maps indicate the cones of influence that delimit the time-frequency regions not influenced by edge effects. In Fig. S1 we report the normalized error producedby different surrogate replicates of original TF distribu-tions. By construction, the histograms of the originaldata are exactly replicated by all the algorithms. Foreach time series, 100 surrogates are generated and theerror computed. Our algorithm is applied in its simpleversion (a single iteration), after ten iterations, and fi-nally for more than 30 iterations. Although the iAAWTalgorithm replicates the multiscale structure of originaldata, our algorithm does a better job in this respect.Even after a single iteration, our algorithm reproducesbetter the non-stationary oscillations of studied time se-ries. The error differences between the iAAWT and ourmethod are all significant as assessed by a non-parametrictest, the Wilcoxon-Mann-Whitney test, with alpha set at p = 0 . Detection of coherent patterns with different sur-rogates. –
We study more in detail the performances ofour algorithm to detect non-stationary spatial coherentcomponents. In addition to the standard stationary sur-rogate data test, we also apply the iAAWT to the samereal-world multivariate data studied in the main text: theepileptic EEG recordings and the measles dataset. Wealso evaluate the detection of significant coherent pat-terns for different number of surrogates generated by thesimplest version of our algorithm (a single iteration). Sig-nificant coherent patterns were detected as statisticallydifferent from those obtained from surrogates, by a z-testcorrected by a FDR at q (cid:54) . −
25 Hzduring practically the whole recording (these fast oscilla-tions are concentrated mainly around 0 < t <
A) Original data sc a l e ( w ee k s ) D) 10 surrogates Non-stationary (CWT) surrogate data tests sc a l e ( w ee k s ) sc a l e ( w ee k s ) sc a l e ( w ee k s ) sc a l e ( w ee k s ) E) 100 surrogates F) 500 surrogatesB) Stationary surrogate data test(100 surrogates) C) Non-stationary (DWT) surrogate data(100 surrogates)
FIG. S3. Ψ( t, f ) values estimated from measles data and statistical significances obtained from different surrogate data tests.Missing values in each original time series were imputed using a local average, i.e. the mean of the two neighboring time points.Same stipulations as in the caption of Figure S2 tionary bispectrum, dynamic phase or cross-frequency in-teractions, among others [S8]. [S1] M. Breakspear, M. Brammer, and P. A. Robinson, Phys-ica D , 1 (2003).[S2] C. Angelini, D. Cava, G. Katul, and B. Vidakovic, Phys-ica D , 24 (2005).[S3] C. J. Keylock, Phys. Rev. E , 036707 (2006); PhysicaD , 219 (2007); Nonl. Processes in Geophys. , 435(2008).[S4] M. Paluˇs, Phys. Rev. Lett. , 134101 (2008). [S5] C. J. Keylock, Phys. Rev. E , 032123 (2017).[S6] S. Mallat, A wavelet tour of signal processing (AcademicPress, San Diego, 1998); B. Cazelles, M. Chavez, D.Berteaux, F. M´enard, J. O. Vik, S. Jenouvrier, and N.C. Stenseth, Oecologia , 287 (2008).[S7] L. Sun, E. Y. Klein, and R. Laxminarayan, Clin. Infect.Dis. , 687 (2012).[S8] T. Subba Rao, and K. C. Indukumar, J. Franklin Inst. , 425 (1996); J. Jamˇsek, A. Stefanovska, and P. V. E.McClintock, Phys. Rev. E , 046221 (2007); J. Jamˇsek,M. Paluˇs, and A. Stefanovska, Phys. Rev. E , 036207(2010); J. A. Schulte, Nonlin. Processes Geophys.,23