Measuring Neural Synchrony by Message Passing Justin Dauwels Amari Research Unit RIKEN Brain Science Institute Wako-shi, Saitama, Japan justin@dauwels.com Francois Vialatte, Tomasz Rutkowski, Andrzej Cichocki ¸ Advanced Brain Signal Processing Lab RIKEN Brain Science Institute Wako-shi, Saitama, Japan {fvialatte,tomek,cia}@brain.riken.jp Abstract A novel approach to measure the interdependence of two time series is proposed, referred to as "stochastic event synchrony" (SES); it quantifies the alignment of two point processes by means of the following parameters: time delay, variance of the timing jitter, fraction of "spurious" events, and average similarity of events. SES may be applied to generic one-dimensional and multi-dimensional point processes, however, the paper mainly focusses on point processes in time-frequency domain. The average event similarity is in that case described by two parameters: the average frequency offset between events in the time-frequency plane, and the variance of the frequency offset ("frequency jitter"); SES then consists of five parameters in total. Those parameters quantify the synchrony of oscillatory events, and hence, they provide an alternative to existing synchrony measures that quantify amplitude or phase synchrony. The pairwise alignment of point processes is cast as a statistical inference problem, which is solved by applying the max-product algorithm on a graphical model. The SES parameters are determined from the resulting pairwise alignment by maximum a posteriori (MAP) estimation. The proposed interdependence measure is applied to the problem of detecting anomalies in EEG synchrony of Mild Cognitive Impairment (MCI) patients; the results indicate that SES significantly improves the sensitivity of EEG in detecting MCI. 1 Introduction Synchrony is an important topic in neuroscience. For instance, it is hotly debated whether the synchronous firing of neurons plays a role in cognition [1] and even in consciousness [2]. The synchronous firing paradigm has also attracted substantial attention in both the experimental (e.g., [3]) and the theoretical neuroscience literature (e.g., [4]). Moreover, medical studies have reported that many neurophysiological diseases (such as Alzheimer's disease) are often associated with abnormalities in neural synchrony [5, 8]. In this paper, we propose a novel measure to quantify the interdependence between point processes, referred to as "stochastic event synchrony" (SES); it consists of the following parameters: time delay, variance of the timing jitter, fraction of "spurious" events, and average similarity of the events. The pairwise alignment of point processes is cast as a statistical inference problem, which is solved by applying the max-product algorithm on a graphical model [6]. In the case of one-dimensional point processes, the graphical model is cycle-free and statistical inference is exact, whereas for multi-dimensional point processes, exact inference becomes intractable; the max-product algorithm is then applied on a cyclic graphical model, which not necessarily yields the optimal alignment [6]. Our experiments, however, indicate that the it finds reasonable alignments in practice. The SES parameters are determined from the resulting pairwise alignments by maximum a posteriori (MAP) estimation. The proposed method may be helpful to detect mental disorders such as Alzheimer's disease, since mental disorders are often associated with abnormal blood and neural activity flows, and changes in the synchrony of brain activity (see, e.g., [5, 8]). In this paper, we will present promising results on the early prediction of Alzheimer's disease from EEG signals based on SES. This paper is organized as follows. In the next section, we introduce SES for the case of one-dimensional point processes. In Section 3, we consider the extension to multi-dimensional point processes. In Section 4, we use our measure to detect abnormalities in the EEG synchrony of Alzheimer's disease patients. 1 2 One-Dimensional Point Processes Let us consider the one-dimensional point processes ("event strings") X and X in Fig. 1(a) (ignore Y and Z for now). We wish to quantify to which extent X and X are synchronized. Intuitively speaking, two event strings can be considered as synchronous (or "locked") if they are identical apart from: (i) a time shift t ; (ii) small deviations in the event occurrence times ("event timing jitter"); (iii) a few event insertions and/or deletions. More precisely, for two event strings to be synchronous, the event timing jitter should be significantly smaller than the average inter-event time, and the number of deletions and insertions should comprise only a small fraction of the total number of events. This intuitive concept of synchrony is illustrated in Fig. 1(a). The event string X is obtained from event string X by successively shifting X over t (resulting in Y ), slightly perturbing the event occurrence times (resulting in Z ), and eventually, by adding (plus sign) and deleting (minus sign) events, resulting in X . Adding and deleting events in Z leads to "spurious" events in X and X (see Fig. 1(a); spurious events are marked in red): a spurious event in X is an event that cannot be paired with an event in X and vice versa. The above intuitive reasoning leads to our novel measure for synchrony between two event strings, i.e., "stochastic event synchrony" (SES); for the one-dimensional case, it is defined as the triplet (t , st , spur ), where st is the variance of the (event) timing jitter, and spur is the percentage of spurious events spur = nspur + n pur s , n + n (1 ) with n and n the total number of events in X and X respectively, and nspur and n pur the total number s of spurious events in X and X respectively. SES is related to the metrics ("distances") proposed in [9]; those metrics are single numbers that quantify the synchrony between event strings. In contrast, we characterize synchrony by means of three parameters, which allows us to distinguish different types of synchrony (see [10]). We compute those three parameters by performing inference in a probabilistic model. In order to describe that model, we consider Fig. 1(b), which shows a symmetric procedure to generate X and X . First, one generates an event string V of length , where the events Vk are mutually independent and uniformly distributed in [0, T0 ]. The strings Z and Z are generated by delaying V over -t /2 and t /2 respectively and by (slightly) perturbing the resulting event occurrence times (variance of timing jitter equals st /2). The sequences X and X are obtained from Z and Z by removing some of the events; more precisely, from each pair (Zk , Zk ), either Zk or Zk is removed with probability ps . This procedure amounts to the statistical model: p(x, x , b, b , v , t , st , ) = p(x|b, v , t , st )p(x |b , v , t , st )p(b, b |)p(v |)p()p(t )p(st ), (2 ) where b and b are binary strings that indicate whether the events in X and X are spurious (Bk = 1 if Xk is spurious, Bk = 0 otherwise; likewise for Bk ); the length has a geometric prior p() = (1 - ) - with (0, 1), and p(v |) = T0 . The prior on the binary strings b and b is given by p(b, b |) = (1 - ps )n+n p2-n-n = (1 - ps )n+n ps spur , s ntot (3 ) with ntsotur = nspur + n pur = 2 - n - n the total number of spurious events in X and X , nspur the p s number of spurious events in X : kn nspur = b k = - n , (4 ) =1 and similarly, n pur s the number of spurious events in X : n pur s = kn b = - n. k (5 ) =1 The conditional distributions in X and X are equal to: kn N x t st p(x|b, v , t , st ) = k - vik ; - , 22 =1 p(x |b , v , t , t ) = kn =1 N where Vik is the event in V that corresponds to Xk (likewise Vi ), and N (x; m, s) is a univariate k Gaussian distribution with mean m and variance s. Since we do not wish/need to encode prior information about t and st , we adopt improper priors p(t ) = 1 = p(st ). 2 x 2 t t , k - vi ; k 22 1-bk 1-bk (6 ) , (7 ) with (xjk , x ) the pairs of non-spurious events, nnon-spur = n + n - the total number of nonjk spurious event pairs, and = ps T0 ; in the example of Fig. 1(b), J = (1, 2, 3, 5, 6, 7, 8), J = (2, 3, 4, 5, 6, 7, 8), and nnon-spur = 7. In the following, we will denote model (8) by p(x, x , j, j , t , st ) instead of p(x, x , b, b , t , st , ), since for given x, x , b, and b (and hence given n, n , and nnon-spur ), the length is fully determined, i.e., = n + n - nnon-spur ; moreover, it is more natural to describe the model in terms of J and J instead of B and B (cf. RHS of (8)). Note that B and B can directly be obtained f r o m J an d J . It also noteworthy that T0 , and ps do not need to be specified individually, since they appear in (8) only through . The latter serves in practice as a knob to control the number of spurious events. I B X X Z V t Z X X B I 1 1 0 2 00 34 0 6 0 00 7 89 Z t 2 Eventually, marginalizing (2) w.r.t. v results in the model: p nnon-spur k t ot p(x, x , b, b , t , st , ) = (x, x , b, b , v , t , st , )dv nspur N (x - xjk ; t , st ), j k (8 ) =1 2 0 34 00 5 1 6 0 78 9 00 0 0 t 2 T0 Y (a) Asymmetric procedure (b) Symmetric procedure Figure 1: One-dimensional stochastic event synchrony. Given event strings X and X , we wish to determine the parameters t and st , and the hidden variables B and B ; the parameter spur (cf. (1)) can obtained from the latter : n n k =1 b k + k =1 b k spur = . (9 ) n+n There are various ways to solve this inference problem, but perhaps the most natural one is cyclic maxi(0) ^(0) mization: first one chooses initial values t and st , then one alternates the following two update rules ^ until convergence (or until the available time has elapsed): ^(i) ^(i) (^(i+1) , ^(i+1) ) = argmax p(x, x , j, j , , s ) j j (1 0 ) b,b t t ^ (t (i+1) , st ^ (i+1) j j ) = argmax p(x, x , ^(i+1) , ^(i+1) , t , st ). t ,st (1 1 ) The update (11) is straightforward: ^(i+1) t = 1 nnon-spur 1 (i+1) nnon-spur (i+1) i+ ) n(on-s1ur n p i+ ) n(on-s1ur n p k k x (i+1) - xj (i+1) , j =1 k k (1 2 ) (i+1) st ^ = ^ (xj (i+1) - xj (i+1) - t k k (i+1) 2 ). (1 3 ) =1 The update (10) can readily be carried out by applying the Viterbi algorithm ("dynamic programming") on an appropriate trellis (with the pairs of non-spurious events (xjk , x ) as states), or equivalently, by jk applying the max-product algorithm on a suitable factor graph [6]; the procedure is similar to dynamic time warping [7]. 3 Multi-Dimensional Point Processes In this section, we will focus on the interdependence of multi-dimensional point processes. As a concrete example, we will consider multi-dimensional point processes in time-frequency domain; the proposed algorithm, however, is not restricted to that particular situation, it is applicable to generic multidimensional point processes. 3 Suppose that we are given a pair of (continuous-time) signals, e.g., EEG signals recorded from two different channels. As a first step, the time-frequency ("wavelet") transform of each signal is approximated as a sum of (half-ellipsoid) basis functions, referred to as "bumps" (see Fig. 2 and [17]); each bump is described by five parameters: time X , frequency F , width X , height F , and amplitude W . The resulting bump models Y = ((X1 , F1 , X1 , F1 , W1 ), . . . , (Xn , Fn , Xn , Fn , Wn )) and Y = ((X1 , F1 , X1 , F1 , W1 ), . . . , (Xn , Fn , Xn , Fn , Wn )), representing the most prominent oscillatory activity, are thus 5-dimensional point processes. Our extension of stochastic event synchrony to multi-dimensional point processes (and bump models in particular) is derived from the following observation (see Fig. 3): bumps in one time-frequency map may not be present in the other map ("spurious" bumps); other bumps are present in both maps ("non-spurious bumps"), but appear at slightly different positions on the maps. The black lines in Fig. 3 connect the centers of non-spurious bumps, and hence, visualize the offset between pairs of non-spurious bumps. We quantify the interdependence between two bump models by five parameters, i.e., the parameters spur , t , and st introduced in Section 2, in addition to f , the average frequency offset between non-spurious bumps, and sf , the variance of the frequency offset between non-spurious bumps. We determine the alignment of two bump models in addition to the 5 above parameters by an inference algorithm similar to the one of Section 2, as we will explain in the following; we will use the notation = (t , st , f , sf ). Model (8) may naturally be extended in time-frequency domain as: p(y , y , j, j , ) ntsotur p · where the offset x - xk in time and offset fk - fk in frequency are normalized by the width and height k respectively of the bumps; we will elaborate on the priors on the parameters later on. In principle, one may determine the sequences J and J and the parameters by cyclic maximization along the lines of (10) and (11). In the multi-dimensional case, however, the update (10) is no longer tractable: one needs to allow permutations of events, the indices jk and jk are no longer necessarily monotonically increasing, and as a consequence, the state space becomes drastically larger. As a result, the Viterbi algorithm (or equivalently, the max-product algorithm applied on cycle-free factor graph of model (14)) becomes impractical. x - x N f - f k k k k ; t , st ; f , sf xk + x fk + fk k =1 sp s, p(t )p t (f )p f k N nnon-spur (1 4 ) We solve this problem by applying the max-product algorithm on a cyclic factor graph of the system at hand, which will amount to a suboptimal but practical procedure to obtain pairwise alignments of multidimensional point processes (and bump models in particular). To this end, we introduce a representation of model (14) that is naturally represented by a cyclic graph: for each pair of events Yk and Yk , we introduce a binary variable Ckk that equals one if Yk and Yk form pair of non-spurious events and is zero otherwise. Since each event in Y associated to at most one event in Y , we have the constraints: kn C1k = S1 {0, 1}, =1 kn C2k = S2 {0, 1}, . . . , =1 kn Cnk = Sn {0, 1}, (1 5 ) =1 and similarly, each event in Y is associated to at most one event in Y , which is expressed by a similar set of constraints. The sequences S and S are related to the sequences B and B (cf. Section 2) as follows: Bk = 1 - Sk an d Bk = 1 - Sk . (1 6 ) In this representation, the global statistical model (14) can be cast as: p(y , y , b, b , c, ) kn kn kn · · =1 =1 N kn ( [bk - 1] + [bk ]) =1 kn ( [b - 1] + [b ]) k k ck k sp s p(t )p t (f )p f (1 7 ) =1 =1 [bk + kn N f - f x - x k k k k ; t , st ; f , sf xk + x fk + fk k ck k - 1 ] =1 kn =1 [b + k kn ck k - 1 ] =1 . Since we do not need to encode prior information about t and f , we choose improper priors p(t ) = 1 = p(f ). On the other hand, we have prior knowledge about st and sf . Indeed, we expect a bump in one time-frequency map to appear in the other map at about the same frequency, but there may be some timing offset between both bumps. For example, bump nr. 1 in Fig. 3(a) (t = 10.7s) should be paired with bump nr. 3 (t = 10.9s) and not with nr. 2 (t = 10.8s), since the former is much closer in frequency than 4 the latter. As a consequence, we a priori expect smaller values for sf than for st . We encode this prior information by means of conjugate priors for st and sf , i.e., scaled inverse chi-square distributions. A factor graph of model (17) is shown in Fig. 4 (each edge represents a variable, each node corresponds to a factor of (17), as indicated by the arrows at the right hand side; we refer to [6] for an introduction to factor graphs). We omitted the edges for the (observed) variables Xk , Xk , Fk , Fk , Xk , Xk , Fk , and Fk in order not to clutter the figure. Time-frequency map Time-frequency map Bump model Bump model Figure 2: Two-dimensional stochastic event synchrony. We determine the alignment C = (C11 , C12 , . . . , Cnn ) and the parameters = (t , st , f , sf ) by maximum a posteriori (MAP) estimation: (c, ) = argmax p(y , y , c, ), ^^ c, (1 8 ) where p(y , y , c, ) is obtained from (17) by marginalizing over b and b . From c, we obtain the estimate spur as: ^ ^ n n n^ n ^ n + n - 2 k =1 k =1 c k k ^ k =1 b k + k =1 b k = . (1 9 ) spur = ^ n+n n+n The MAP estimate (18) is intractable, and we try to obtain (18) by cyclic maximization: first, the para(0) (0) (0) ^(0) meters are initialized: t = 0 = f , st = s0,t , and sf = s0,f , then one alternates the following ^ ^ ^ two update rules until convergence (or until the available time has elapsed): ^ c(i+1) = argmax p(y , y , c, (i) ) ^ c (2 0 ) (2 1 ) ^ (i+1) = argmax p(y , y , c(i+1) , ). ^ ^ The estimates (i+1) is available in closed-form. Update (20), i.e., finding the optimal pairwise alignment ^ C for given values (i) of the parameters , is less straightforward: it involves an intractable combinatorial optimization problem. We attempt to solve that problem by applying the max-product algorithm to the (cyclic) factor graph depicted in Fig. 4 [6]. Let us first point out that, since the alignment C is com^ ^ puted for given = (i) , the (upward) messages along the edges are the point estimate (i) (cf. (20)); equivalently, for the purpose of computing (20), one may remove the edges and the two bottom nodes in Fig. 4; the N -nodes then become leaf nodes. The other messages in the graph are iteratively updated according to the generic max-product update rule (see Table 1, row 1). In particular, the upward messages along the edges Ckk and the downward messages along the edges Bk and Bk are trivial since they are transmitted by leaf nodes: they are equal to the functions represented by those lead nodes; for example, the downward message along edge Bk is equal to µ(bk ) = [bk ] + [bk - 1]. The max-product update rule for the other messages are non-trivial, they can be found in Table 1 (row 2 and 3). The resulting inference algorithm for computing (20) is summarized in Table 2. The messages µ(ckk ) ¯ and µ (ckk ) propagate upward along the edges ckk towards the -nodes connected to the edges 30 25 20 15 10 5 00 30 25 20 15 10 5 00 1 2 3 f [ Hz ] f [ Hz ] 5 10 t [s] 15 20 5 10 t [s] 15 20 (a) Bump models of two EEG channels. (b) Non-spurious bumps (spur = 27%); the black lines connect the centers of non-spurious bumps. Figure 3: Spurious and non-spurious activity. 5 B1 B1 B2 B2 ¯ ... ¯ Bn Bn [bn ] + [bn - 1] µ ¯ ¯ ¯ ¯ [bn + µ µ µ n k =1 cnk - 1] = C11 N = C12 N ... = = C21 N = C22 N ... = ... = Cn1 N N = Cn2 N ... = ttt C1n ... N C2n ... N Cnn ... N cn n ^ (k ) = (t , st , f , sf ) = ^ (k ) -x n n x n + x ; t , s t n x N -fn n fn +fn f ; f , sf p(t , st , f , sf ) = p(t )p(st )p(f )p(sf ) Figure 4: Factor graph of model (17). Bk and Bk respectively (see Fig. 4, left hand side); the messages µ (ckk ) and µ (ckk ) propagate ¯ downward along the edges ckk from the -nodes connected to the edges Bk and Bk respectively. After initialization (22) of the messages µ (ckk ) and µ (ckk ) (k = 1, 2, . . . , n; k = 1, 2, . . . , n ), one alternatively updates (i) the messages µ(ckk ) (23) and µ (ckk ) (24), (ii) the messages µ(ckk ) (25) and µ (ckk ) (26), until convergence; it is noteworthy that, although the max-product algorithm is not guaranteed to converge on cyclic graphs, we observed in our experiments (see Section 4) that alternating the updates (23)­(26) always converged to a fixed point. At last, one computes the marginals p(bk ) (27) and p(b ) (28), and from the latter, one determines the k decisions ^k and ^ (29), which indicate whether the k -th and k -th bump in the first and second sequence b bk respectively are considered to be spurious or not. From the decisions ^, one may obtain decisions c as b ^ follows: If ^k = 1, then ckk = 0 for all k , whereas if ^k = 0, then ckk = 1 for k = argmax µ ^ ^ b b (ck = 1)/µ (ck = 0), and ckk = 0 otherwise. Alternatively, one may obtain decisions c from ^ the decisions ^ : if ^ = 1, then ckk = 0 for all k , whereas if ^ = 0, then ckk = 1 for k = b bk ^ bk ^ argmax µ(ck = 1)/µ(ck = 0), and ckk = 0 otherwise. Both sets of decisions c are not necessarily ^ ^ equal; moreover, the decisions c are not guaranteed to be consistent. In our experiments, however, both ^ ^ solutions c were always identical and consistent. We also applied the sum-product algorithm instead of the max-product algorithm, and observed that the resulting configurations c were in that case not always ^ consistent (consistency in only about 70% of our experiments, cf. Section 4); both the max-product and the sum-product algorithm always converged. X1 X2 . . . Y 1 Xm µ(y ) maxx1 ,x2 ,...,xm f (x1 , . . . , xm , y )µ(x1 ) . . . µ(xm ) f (x1 , x2 , . . . xm , y ) X1 X2 .. . 2 = Y Xm [x1 - x2 ] . . . [xm - y ] µ (y = 0) µ(y = 1) µ (x1 = 0)µ(x2 = 0) . . . µ(xm = 0) µ(x1 = 1)µ(x2 = 1) . . . µ(xm = 1) X 1 , . . . , Xm , Y { 0 , 1 } B B ¯ X1 ... 3 Xm µ (xk = 0) µ(xk = 1) µ (b = 0) µ(b = 1) m µ 1 (b=1) ax µ(b=0) , max=k m axk µ(xk =1) µ(xk =0) µ(x =1) µ(x =0) [b + , X1 , . . . , Xm { 0 , 1 } m k=1 xk - 1] 1 Table 1: Message computation rules for the max-product algorithm. 4 Diagnosis of MCI from EEG We analyzed rest eyes-closed EEG data recorded from 21 sites on the scalp based on the 10­20 system. The sampling frequency was 200 Hz, and the signals were band pass filtered between 4 and 30Hz. The 6 Initialization µ(ckk ) = µ (ckk ) N N f - f x - x k k k k ; t , st ; f , sf xk + x fk + fk k ck k (2 2 ) Iteratively compute messages until convergence A. Downward messages: µ m (ckk = 0) ax ( , max =k µ(ck = 1)/µ(ck = 0)) µ(c = 1) 1 µ kk m (ckk = 0) ax ( , max=k µ (ck = 1)/µ (ck = 0)) µ (ckk = 1) 1 B. Upward messages: x - x N f - f ck k k k k k ; t , st ; f , sf xk + xk fk + fk N x - x N f - f ck k k k k k µ (ckk ) µ(ckk ) ; t , st ; f , sf xk + x fk + fk k µ(ckk ) µ (ckk ) ( 23) 24) ( N (2 5 ) (2 6 ) Marginals (bk p(bk p (bk p(b k Decisions p = 0) = 1) = 0) = 1) bk m axk µ(ckk = 1)/µ(ckk = 0) m axk µ (ckk = 1)/µ (ckk = 0) an d ^ = argmax p(b ) bk k b k ( ( 27) 28) ^k = argmax p(bk ) b (2 9 ) Table 2: Inference algorithm. subjects comprised two study groups: the first consisted of a group of 22 patients diagnosed as suffering from MCI, who subsequently developed mild AD. The other group was a control set of 38 age-matched, healthy subjects who had no memory or other cognitive impairments. Pre-selection was conducted to ensure that the data were of a high quality, as determined by the presence of at least 20s of artifact free data. We computed a large variety of synchrony measures for both data sets; the results are summarized in Table 3. We report results for global synchrony, obtained by averaging the synchrony measures over 5 brain regions (frontal, temporal left and right, central, occipital). For SES, the bump models were clustered by means of the aggregation algorithm described in [17]. The strongest observed effect is a significantly higher degree of background noise (spur ) in MCI patients, more specifically, a high number of spurious, non-synchronous oscillatory events (p = 0.00076). We verified that the SES measures are not correlated (Pearson r) with other synchrony measures (p > 0.10); in contrast to the other measures, SES quantifies the synchrony of oscillatory events (instead of more conventional amplitude or phase synchrony). Combining spur with ffDTF yields good classification of MCI vs. Control patients (see Fig.5(a)). Interestingly, we did not observe a significant effect on the timing jitter st of the non-spurious events (p = 0.069). In other words, AD seems to be associated with a significant increase of spurious background activity, while the non-spurious activity remains well synchronized. Moreover, only the non-spurious activity slows down (p = 0.0012; see Fig.5(c)), the average frequency of the spurious activity is not affected in MCI patients (see Fig.5(c)). In future work, we will verify those observations by means of additional data sets. 0.45 0.4 0.35 MCI CTR 19 18 17 16 15 14 19 18 17 spur fspur 0.3 0.25 0.2 0.15 fnon-spur CTR 16 15 14 13 12 13 0.1 0.045 0.05 Fi2 j 0.055 0.06 12 MCI CTR MCI (a) spur vs. ffDTF (b) Av. frequency of the spurious (c) Av. frequency of the nonactivity (p = 0.87) spurious activity (p = 0.0019) Figure 5: Results. References [1] F. Varela, J. P. Lachaux, E. Rodriguez, and J. Martinerie, "The Brainweb: Phase Synchronization and LargeScale Integration", Nature Reviews Neuroscience, 2(4):229­39, 2001. 7 Measure p-value References Measure p-value References Measure p-value References Measure p-value References Measure p-value References Measure p-value Cross-correlation 0 .0 2 8 Gr a n g e r c o h e r e n c e 0 .1 5 Kullback-Leibler 0 .0 7 2 Nk 0 .0 3 2 Hilbert Phase 0 .1 5 1 5 [24] st 0 .0 6 9 Coherence 0 .0 6 0 [16] Partial Coherence 0 .1 6 Renyi ´ 0 .0 7 6 Sk 0 .2 9 [15] Wavelet Phase 0 .0 8 2 spur 0.00076 Phase Coherence 0 .7 2 PDC 0 .6 [13] Jensen-Shannon 0 .0 8 4 [23] Hk 0 .0 9 0 Evolution Map 0 .0 7 1 5 [19] Corr-entropy 0 .2 7 [18] DT F 0 .3 4 Jensen-Renyi ´ 0 .1 2 S-estimator 0 .3 3 [21] Instantaneous Period 0 .0 2 0 Wave-entropy 0 .0 1 2 [20] ffDTF 0.0012 IW 0 .0 6 0 d DT F 0 .0 3 0 I 0 .0 8 0 [22] Table 3: Sensitivity of synchrony measures for early prediction of AD (p-values for Mann-Whitney test; * and ** indicate p < 0.05 and p < 0.005 respectively). N k , S k , and H k are three measures of nonlinear interdependence [15]. [2] W. Singer, "Consciousness and the Binding Problem," Annals of the New York Academy of Sciences, 929:123­ 146, April 2001. [3] M. Abeles, H. Bergman, E. Margalit, and E. Vaadia, "Spatiotemporal Firing Patterns in the Frontal Cortex of Behaving Monkeys," J. Neurophysiol, 70(4):1629­1638. 1993. [4] S. Amari, H. Nakahara, S. Wu, and Y. Sakai, "Synchronous Firing and Higher-Order Interactions in Neuron Pool," Neural Computation, 15:127­142, 2003. [5] H. Matsuda, "Cerebral Blood Flow and Metabolic Abnormalities in Alzheimer's Disease," Ann. Nucl. Med., vol. 15, pp. 85­92, 2001. [6] H.-A. Loeliger, "An Introduction to Factor Graphs," IEEE Signal Processing Magazine, Jan. 2004, pp. 28­41. [7] C. S. Myers and L. R. Rabiner, "A Comparative Study of Several Dynamic Time-Warping Algorithms for Connected Word Recognition," The Bell System Technical Journal, 60(7):1389­1409, September 1981. [8] J. Jong, "EEG Dynamics in Patients with Alzheimer's Disease," Clinical Neurophysiology, 115:1490­1505 (2004). [9] J. D. Victor and K. P. Purpura, "Metric-space Analysis of Spike Trains: Theory, Algorithms, and Application," Network: Comput. Neural Systems, 8:17, 164, 1997. [10] H. P. C. Robinson, "The Biophysical Basis of Firing Variability in Cortical Neurons," Chapter 6 in Computational Neuroscience: A Comprehensive Approach, Mathematical Biology & Medicine Series, Edited By Jianfeng Feng, Chapman & Hall/CRC, 2003. [11] E. Pereda, R. Q. Quiroga, and J. Bhattacharya, "Nonlinear Multivariate Analysis of Neurophsyiological Signals," Progress in Neurobiology, 77 (2005) 1­37. [12] M. Breakspear, "Dynamic Connectivity in Neural Systems: Theoretical and Empirical Considerations," Neuroinformatics, vol. 2, no. 2, 2004. [13] M. Kaminski and Hualou Liang, "Causal Influence: Advances in Neurosignal Analysis," Critical Review in ´ Biomedical Engineering, 33(4):347­430 (2005). [14] C. J. Stam, "Nonlinear Dynamical Analysis of EEG and MEG: Review of an Emerging Field," Clinical Neurophysiology 116:2266­2301 (2005). [15] R. Q. Quiroga, A. Kraskov, T. Kreuz, and P. Grassberger, "Performance of Different Synchronization Measures in Real Data: A Case Study on EEG Signals," Physical Review E, vol. 65, 2002. [16] P. Nunez and R. Srinivasan, Electric Fields of the Brain: The Neurophysics of EEG, Oxford University Press, 2006. [17] F. Vialatte, C. Martin, R. Dubois, J. Haddad, B. Quenet, R. Gervais, and G. Dreyfus, "A Machine Learning Approach to the Analysis of Time-Frequency Maps, and Its Application to Neural Dynamics," Neural Networks, 2007, 20:194­209. [18] Jian-Wu Xu, H. Bakardjian, A. Cichocki, and J. C. Principe, "EEG Synchronization Measure: a Reproducing Kernel Hilbert Space Approach," submitted to IEEE Transactions on Biomedical Engineering Letters, Sept. 2006. [19] M. G. Rosenblum, L. Cimponeriu, A. Bezerianos, A. Patzak, and R. Mrowka, "Identification of Coupling Direction: Application to Cardiorespiratory Interaction," Physical Review E, 65 041909, 2002. [20] C. S. Herrmann, M. Grigutsch, and N. A. Busch, "EEG Oscillations and Wavelet Analysis," in Todd Handy (ed.) Event-Related Potentials: a Methods Handbook, pp. 229-259, Cambridge, MIT Press, 2005. [21] C. Carmeli, M. G. Knyazeva, G. M. Innocenti, and O. De Feo, "Assessment of EEG Synchronization Based on State-Space Analysis," Neuroimage, 25:339­354 (2005). [22] A. Kraskov, H. Stogbauer, and P. Grassberger, "Estimating Mutual Information," Phys. Rev. E 69 (6) 066138, ¨ 2004. [23] S. Aviyente, "A Measure of Mutual Information on the Time-Frequency Plane," Proc. of ICASSP 2005, vol. 4, pp. 481­484, March 18­23, 2005, Philadelphia, PA, USA. [24] J.-P. Lachaux, E. Rodriguez, J. Martinerie, and F. J. Varela, "Measuring Phase Synchrony in Brain Signals," Human Brain Mapping 8:194208 (1999). 8