Robustly estimating the flow direction of information in complex physical systems

Robustly estimating the flo w direction of information in complex ph ysical systems Guido Nolte ∗ F r aunhofer FIRST.IDA, Kekul ´ estr asse 7, D-12489 Berlin, Germany Andreas Ziehe † T e chnical University of Berlin, Computer Scienc e, Machine L e arning L ab or atory, F r anklinstr. 28/29, 10587 Berlin, Germany V adim V. Nikulin ‡ Dept of Neur ology, Campus Benjamin F r anklin, Charite University Me dicine Berlin, D-12203, Germany Alois Schl¨ ogl § F r aunhofer FIRST.IDA, Kekul ´ estr asse 7, D-12489 Berlin, Germany Nicole Kr¨ amer ¶ T e chnical University of Berlin, Computer Scienc e, Machine L e arning L abor atory, F ranklinstr. 28/29, 10587 Berlin, Germany T om Brismar ∗∗ Kar olinska Institutet, Clinic al Neur ophysiolo gy, Kar olinska Hospital, S-17176 Sto ckholm, Sweden Klaus-Rob ert M ¨ uller †† T e chnical University of Berlin, Computer Scienc e, Machine L e arning L ab or atory, F r anklinstr. 28/29, 10587 Berlin, Germany and F r aunhofer FIRST.IDA, Germany (Dated: F ebruary 14, 2013) W e propose a new measure to estimate the direction of information flux in m ultiv ariate time series from complex systems. This measure, based on the slop e of the phase spectrum (Phase Slope Index) has in v ariance properties that are important for applications in real physical or biological systems: (a) it is strictly insensitive to mixtures of arbitrary indep enden t sources, (b) it gives meaningful results even if the phase spectrum is not linear, and (c) it prop erly weigh ts con tributions from differen t frequencies. Simulations of a class of coupled multiv ariate random data show that for truly unidirectional information flow without additional noise contamination our measure detects the correct direction as goo d as the standard Granger causalit y . F or random mixtur es of indep endent sources Granger Causalit y err one ously yields highly significan t results whereas our measure c orr e ctly b ecomes non-significan t. An application of our no vel metho d to EEG data (88 sub jects in eyes-closed condition) reveals a strikingly clear fron t-to-bac k information flow in the v ast ma jority of sub jects and thus contributes to a b etter understanding of information pro cessing in the brain. P ACS numbers: 05.45.Xt,05.45.Tp,87.19 Keywords: Causality , Granger Causalit y , Phase Slope, Coherence, Electro encephalography (EEG), Alpha Rhythm, Noise T o understand in teracting systems it is of fundamental imp ortance to distinguish the driv er from the recipient, and hence to b e able to estimate the direction of information flo w. If one cannot in terfere with the system, the direc- tion can b e estimated with a temp oral argument: the driver ∗ Electronic address: nolte@first.fraunhofer.de † Electronic address: ziehe@first.fraunhofer.de ‡ Electronic address: v adim.nikulin@c harite.de § Electronic address: alois.schloegl@first.fraunhofer.de ¶ Electronic address: nkraemer@cs.tu-b erlin.de ∗∗ Electronic address: tom.brismar@ki.se †† Electronic address: krm@cs.tu-b erlin.de is earlier than the recipient from which it follo ws that the driv er contains information ab out the future of the recipient not contained in the past of the recipient while the reverse is not the case. This argument is the conceptual basis of Granger Causalit y [1, 2] whic h is probably the most promi- nen t metho d to estimate the direction of causal influence in time series analysis. Granger Causality was originally devel- op ed in econometry , but is applied to many differen t prob- lems in physics, geosciences (cause of climate change), so cial sciences, and biology with sp ecial emphasis on neural system [3, 4, 5, 6, 7]. The difficulty in realistic measurements in complex sys- tems is that asymmetries in detection pow er may as well arise due to other factors, specifically independent background ac- 2 tivit y ha ving nontrivial sp ectral prop erties and even tually b eing measured in unkno wn sup erp osition in the channels. In this case the in terpretation of the asymmetry as a direc- tion of information flow can lead to significant alb eit false results [8]. The purpose of this pap er is to prop ose a nov el estimate of flux direction whic h is highly robust against false estimates caused b y confounding factors of very general na- ture. More formally , we are interested in statistical dep enden- cies in complex ph ysical systems and esp ecially in causal re- lations betw een a signal of interest consisting of tw o sources with time series x i ( t ) for i = 1 , 2. The measured data y ( t ) are assumed to b e a sup erp osition of these sources of in terest and additive noise η ( t ) in the form y ( t ) = x ( t ) + B η ( t ) (1) where η ( t ) is a set of M indep endent noise sources which are mixed into the measuremen t channels b y an unknown 2 × M mixing matrix B . The prop osed metho d is based on the slop e of the phase of cross-sp ectra b etw een tw o time series. A fixed time de- la y for an in teraction betw een t wo systems will affect differ- en t frequency comp onents in different wa ys. This is most easily seen if we assume that the interaction is merely a dela y by a time τ , i.e. y 2 ( t ) = ay 1 ( t − τ ) with a b eing some constan t. In the F ourier-domain this relation reads ˆ y 2 ( f ) = a exp( − i 2 π f τ ) ˆ y 1 ( f ). F or the cross-spectrum S ij ( f ) b et ween the t w o channels i and j one has S 12 ( f ) = h ˆ y 1 ( f ) ˆ y ∗ 2 ( f ) i ∼ exp( i 2 π f τ ) ≡ exp( i Φ( f )) (2) where h·i denotes exp ectation v alue. The phase-spectrum Φ( f ) = 2 πf τ is linear and proportional to the time dela y τ . The slop e of Φ( f ) can b e estimated, and the causal direction is estimated to go from y 1 to y 2 ( y 2 to y 1 ) if it is p ositiv e (negativ e). The idea here is now to define an av erage phase slop e in suc h a w ay that a) this quantit y prop erly represents rela- tiv e time delays of different signals and esp ecially coincides with the classical definition for linear phase sp ectra , b) it is insensitive to signals whic h do not interact regardless of sp ectral con tent and superp ositions of these signals, and c) it prop erly w eights different frequency regions according to the statistical relev ance. This quantit y is termed ’Phase Slop e Index’ (PSI) and is defined as ˜ Ψ ij = =   X f ∈ F C ∗ ij ( f ) C ij ( f + δ f )   (3) where C ij ( f ) = S ij ( f ) p S ii ( f ) S j j ( f ) (4) is the complex coherency , S is the cross-spectral matrix, δ f is the frequency resolution, and = ( · ) denotes taking the imag- inary part. F is the set of frequencies o v er which the slop e is summed. T o see that the definition of ˜ Ψ ij corresp onds to a meaning- ful estimate of the av erage slop e it is conv enien t to rewrite it as ˜ Ψ ij = X f ∈ F α ij ( f ) α ij ( f + δ f ) sin(Φ( f + δ f ) − Φ( f )) (5) with α ij ( f ) = | C ij ( f ) | b eing frequency dep enden t weigh ts. F or smo oth phase sp ectra, sin(Φ( f + δ f ) − Φ( f )) ≈ Φ( f + δ f ) − Φ( f ) and hence Ψ corresp onds to a weigh ted a v erage of the slop e. W e emphasize that since Ψ v anishes if the imaginary part of coherency v anishes it will b e insensitive to mixtures of non-interacting sources [9, 10]. Finally , it is con venien t to normalize ˜ Ψ by an estimate of its standard deviation Ψ = ˜ Ψ std ( ˜ Ψ) (6) with std ( ˜ Ψ) b eing estimated b y the Jackknife metho d. In the examples b elo w we alwa ys show normalized measures of directionalit y , and w e consider absolute v alues larger than 2 as significant. Estimations of cross -spectra is standard [9, 11] but tec h- nical details may differ. Here, we first divide the whole data in to epo chs containing contin uous data (4 seconds duration), then w e divide each epo ch further into segments of time T , here of 2 seconds duration corresponding to a frequency res- olution of δ f = 0 . 5 Hz, m ultiply the data for each segmen t with a Hanning windo w, F ourier-transform the data, and estimate the cross-sp ectra according to Eq.2 as an av erage o ver all segments. The segmen ts hav e 50% ov erlap within eac h ep o ch but not across ep o chs. T o apply the Jac kknife metho d, for each pair of channels we calculate ˜ Ψ k from data with the k .th epo ch remov ed for all k . The standard devia- tion of ˜ Ψ is finally estimated for K ep o chs as √ K σ where σ is the standard deviation of the set of ˜ Ψ k . Our new method is compared to Granger causalit y using Autoregressiv e (AR) mo dels b oth for wide band and narrow band analysis [12] with analogous normalization by the es- timated standard deviation. T o estimate the parameters of the model w e here use the Levinson-Wiggens-Robinson [13] algorithm av ailable in the op en Biosig to olb ox [14]. Granger Causalit y is defined as the difference b et ween the flux from c hannel 1 to 2 and the flux from c hannel 2 to 1 normalized to unit estimated standard deviation. W e first illustrate typical results for t wo simple cases in Fig.1. The upp er panels show a simulation of a strong inter- action from the second (dashed) to the first channel (solid) generated with a simple AR mo del of order one. The second signal is clearly earlier than the first signal. Both metho ds 3 detect this direction correctly from only 2000 data p oin ts. In the low er panels we show a mixture of pink and white noise. In contrast to PSI, Granger causality erroneously still detects a significant direction. 0 10 20 30 −0.015 −0.01 −0.005 0 0.005 0.01 Amplitude [a.u.] Unidirectional flux −50 0 50 G Ψ 0 10 20 30 −0.015 −0.01 −0.005 0 0.005 0.01 Time [bins] Amplitude [a.u.] Correlated noise −10 −5 0 5 10 G Ψ FIG. 1: Upper panels: strong interaction from second (red in left panels) to first (blue in left panels) signal. Low er panels: mixture of brown and white noise. The error bars in the right panels indicate estimated 95% error margins corresp onding to 2 standard deviations. Time series in the left panels were upsampled to 400 Hz. T o study a more general class of signals we sim ulated data with structure y ( t ) = (1 − γ ) x ( t ) N x + γ B η ( t ) N η (7) Here, the signal x ( t ) contains truly unidirectional informa- tion flux and is generated using AR-mo dels of order P = 5 for tw o channels. In general, an AR-model is defined as z ( t ) = P X p =1 A ( p ) z ( t − p ) + ξ ( t ) (8) where A ( p ) are the AR-matrices up to order P and ξ ( t ) is white Gaussian noise with co v ariance matrix Σ chosen here to be the identit y matrix. F or computing Granger Causalit y , the AR mo del w as fitted with order P = 10. All en tries of AR-matrices w ere selected randomly as inde- p enden t Gaussian random num be rs with A 21 ( p ) = 0 for the signal part x ( t ), corresp onding to unidirectional flow from the second to first signal, and A 12 ( p ) = A 21 ( p ) = 0 for the noise part η ( t ), corresp onding to indep endent sources. Noise w as mixed into c hannels with a random 2 × 2 mixing matrix B . Both the signal part and the mixed noise part ( B η ( t )) are normalized b y the F rob enius norms of the respective data matrices ( N x and N η ) and finally added with a parameter γ controlling for the relativ e strength. The time constant implicit in the AR-mo del w as assumed to b e 10 ms, and w e generated 60000 data p oin ts for each system and c han- nel. This corresponds to a Nyquist frequency of 50 Hz and to 10 minutes measurement. W e analyzed systems for all γ in the range [0 , 1] with step 0 . 1. F or each γ we analyzed 1000 randomly selected stable systems with b oth metho ds and b oth for wide band (using all frequencies) and narrow band analysis. F or the narrow band, we used a band of 5 Hz width, centered this band around the sp ectral p eak of the (kno wn) signal of interest and analyzed only cases where the band includes at least 60% of the total p ow er of the signal of interest. 0 0.5 1 0 20 40 60 80 100 % Significant false, wide band Ψ G 0 0.5 1 0 20 40 60 80 100 Significant false, narrow band % 0 0.5 1 0 20 40 60 80 100 Significant correct, wide band γ % 0 0.5 1 0 20 40 60 80 100 Significant correct, narrow band γ % FIG. 2: F raction of significant detections of Granger Causality and PSI as a function of noise lev el γ . The fractions of significant false and significant correct detections as a function of γ are shown in Fig.2. W e observe that for increasing noise level the fraction of significant false detections for Granger Causality comes close to 50% while PSI rarely mak es significan t false detection at all. F or PSI, the worst case observ ed is at γ = 0 . 8 for the wide band with 6% significant false detections. This level can be reduced to ab out 3 . 5% if w e increase the frequency resolution to 0 . 25Hz. Ho wev er, the price is some loss in statistical p ow er and it is imp ortan t to sho w that also the prop osed metho d migh t fail, ev en if it is unlikely in the sense of the present simulation. W e observ e similar significan t correct detection rates for b oth metho ds for small and moderate noise lev els. F or high noise lev el Granger Causality shows a m uc h larger fraction of significan t correct detections which, how ev er, is meaningless giv en the large fraction of significant false detections. After having illustrated the robustness of our new metho d on simulated data w e now apply the PSI to real data, namely 4 EEG. F or this, 88 healthy sub jects were recruited randomly b y the aid of the Swedish p opulation register. During the exp erimen t, which lasted for 15 min utes, the sub jects w ere instructed to relax and keep their eyes closed. Ev ery minute the s ub jects were asked to open their ey es for 5 seconds. EEG was measured with standard 10-20 system consisting of 19 channels. Data were analysed using link ed mastoids reference. The protocol was approv ed b y the Hospital Ethics Committee. The most prominent feature of this measurement is the alpha p eak at around 10 Hz. This rh ythm is believed to represen t a cortico-cortical or thalamo-cortical in teraction. The direction of this interaction is an open question. While it is mostly b elieved that this rh ythm originates in o ccipital areas and spreads to other (more fron tal) areas [15] this view has also b een c hallenged [16]. FIG. 3: Upp er left: signal p ow er as a function of frequency av er- aged ov er the tw o occiptial channels O1 and O2 showing a clear alpha p eak at f = 9 . 5 Hz. Upper right: net information flux at f = 9 . 5 Hz. Lo wer panels: PSI for all channel pairs and all frequencies pro jected on right-to-left and front-to-bac k direction, resp ectiv ely . F or illustration we show results for PSI for one selected sub ject in Fig.3. The p ow er (upp er left panel), av eraged o ver the tw o o ccipital channels (O1 and O2), shows a very strong p eak at 9 . 5 Hz. PSI v alues were calculated for all c hannel pairs with frequency resolution 0 . 5 Hz using a fre- quency band of 5Hz width centered around frequency f . In the upp er right panel we show the net information flux at f = 9 . 5 Hz defined for the i.th c hannel by Ψ net ( i, f ) = P j Ψ ij ( f ) std ( P j Ψ ij ( f ) (9) W e clearly observe that fron tal channels are net drivers ( FIG. 4: Phase Slop e Index for all pairs of channels av eraged ov er all sub jects each at the p eak of the alpha rhythm. The i.th small circle is located at the i.th electro de position and is a contour plot of the i.th row of the matrix with elements Ψ ij . The red color in fron tal circles indicates that the frontal electro des are estimated as the drivers. Ψ net > 0) and o ccipital c hannels net recipients (Ψ net < 0). T o show prefered direction for all pairs of c hannel and for all frequencies w e calculate the resp ective contribution to a giv en direction in the following wa y: for channels i and j with lo cations r i and r j in the tw o dimensional plane (as sho wn in the upp er right panel) resp ectiv ely , we calculate the normalized difference v ector δ r ij = r j − r i | r i − r j | (10) and pro ject it on to the direction of interest, i.e. onto u = ( − 1 , 0) T for right to left direction and onto u = (0 , − 1) T for fron t to bac k direction. W e finally calculate the con tribution of Ψ ij ( f ) to direction u as Ψ i,j ( f , u ) = Ψ ij ( f ) u · δ r ij (11) Results for all c hannel pairs and for all frequencies are sho wn for right-left information flow (low er left panel) and for front-bac k information flow (low er righ t panel). W e do not observe an y prefered direction in the right-left flo w. In con trast, the information flo w in front-bac k direction shows a clear positive plateau at the alpha-frequency (indicated with letter ’B’) meaning that typically the frontal channels are estimated as the driv ers. W e also observ e a p ositive and negativ e peak (indicated with letters ’A’ and ’C’) at frequencies around 7 Hz and 12 Hz, resp ectively . Note that these p eaks differ by the width of the frequency band. They 5 are clearly artefacts due to inadequate settings of the band. Sp ecifically , the alpha rh ythm has a prefered phase (for given c hannel pair) which is irrelev an t for the estimation of the phase slop e unless the alpha-p eak is right at the edge of the frequency band such that the band cov ers only half of the system. W e found a similar structure in about 60% of the sub jects. An av erage ov er all sub jects now sho wing information flux b et ween all sub jects is shown in Fig.4. W e also found a substan tial inter-sub ject v ariability , b oth with regard to PSI and actual phase at the alpha peak. The origin is interesting but so far unclear and go es b eyond the scop e of this letter. Note that Granger Causality did not yield any consistent spatial pattern, presumably for reasons of high false negative rates similar to the ones observed in Fig.2. Recen t neuroimaging studies hav e c hallenged a simple view on a rest condition by showing a presence of default states in the cortex, which displa y complex patterns of neu- ronal activ ation [17, 18]. W e here sho w that not only sp e- cific areas are co-activ ated during rest state, but they also demonstrate at a gross lev el a preferential ”default” mo de of information flow in the cortex. Imp ortantly , the drivers of suc h flow are mostly situated in the fron tal areas, from where man y top-do wn atten tional influences are though t to b e originated [19]. In agreemen t with the ab ov e men tioned imaging studies, our study suggests that the maintenance of vigilance is a pro cess displa ying a co ordination of neuronal activit y with w ell defined driv ers and recipients of informa- tion flow. T o conclude, w e presen ted a new metho d to estimate the direction of causal relations from time series’ based on the phase slop e of the cross-sp ectra. While it is w ell known that this slop e is an indicator of the direction, the crucial p oin t here is that we defined an av erage of the phase slop e suc h that this av erage is insensitive to arbitrary mixtures of inde- p enden t sources with arbitrary sp ectra. W e verified the claimed prop erties of the PSI for random linear systems also sho wing that the most prominen t metho d to estimate direction of information flo w, Granger Causalit y , is highly sensitive to mixtures of indep enden t noise sources. Additionally , w e show ed that in situations with com bined unidirectional flow and undirected noise our metho d cor- rectly distinguished the t wo phenomena - in sharp con trast to Granger Causalit y . A final application of our method to real EEG data sho ws significant and meaningful results from the neuroph ysiological p oint of view and underlines the ver- satilit y of our new metho d as a universal to ol for estimating causal flow in complex physical systems that consist of mix- tures of sub comp onents. Ac kno wledgements. W e ac kno wledge partial funding from DFG, BMBF and EU. [1] C. Granger, Econometrica 37 , 424 (1969). [2] C. Granger, Journal of Econ. Dynamics Con trol 2 , 329 (1980). [3] R. Kaufmann and D. Stern, Nature 388 , 39 (1997). [4] P . Naray an and R. Smyth, Applied Economics 38 , 563 (2006). [5] D. Marinazzo, M. P ellicoro, and S. Stramaglia, Physical Re- view E 73 , No.066216 (2006). [6] A. Bro velli, M. Ding, A. Ledb erg, Y. Chen, R. Nak amura, and S. Bressler, Pro ceedings of the National Academ y of Sci- ences of the United States of America 101 , 9849 (2004). [7] J. Sato, E. Amaro, D. T ak ahashi, M. F elix, M. Brammer, and P . Morettin, Neuroimage 31 , 187 (2006). [8] Z. Alb o, G. D. Prisco, Y. Chen, G. Rangara jan, W. T ruc- colo, J. F eng, R. V ertes, and M. Ding, Biol. Cybern. 90 , 318 (2004). [9] G. Nolte, O. Bai, L. Wheaton, Z. Mari, S. V orbac h, and M. Hallett, Clin. Neurophysiol. 115 , 2292 (2004). [10] G. Nolte, F. Meineck e, A. Ziehe, and K. M ¨ uller, Phys Rev E 73 , 051913 (2006). [11] P . Nunez, R. Sriniv asan, A. W estdorf, R. Wijesinghe, D. T uc ker, R. Silb erstein, and P . Cadusch, Electro en- cephalogr. Clin. Neurophysiol. 103 , 499 (1997). [12] M. Ding, Y. Chen, and S. Bressler, in Handb ook of Time Series Analysis (Whiley , 2006), pp. 437–459. [13] S. Marple, Digital Sp e ctr al Analysis with Applic ations (Pren- tice Hall, Englewoo d Cliffs, NJ, 1987). [14] A. Schl¨ ogl, BIOSIG - an op en source softw are library for biomedical signal pro cessing, http://BIOSIG.SF.NET. [15] F. L. da Silv a, Electro encephal. and Clin. Neurophys. 79 , 81 (1991). [16] J. Ito, A. Nik olaev, and C. v an Leeuw en, Biol. Cyb ern. 92 , 54 (2005). [17] H. Laufs, K. Krako w, P . Sterzer, E. Eger, A. Beyerle, A. Salek-Haddadi, and A. Kleinsc hmidt, Pro c Natl Acad Sci USA 100 , 11053 (2003). [18] D. Man tini, M. Perucci, C. D. Gratta, G. Romani, and M. Corb etta, Pro c Natl Acad Sci USA 104 , 13170 (2007). [19] C. Gilb ert and M. Sigman, Neuron 54 , 677 (2007).