Qephd: Coupling Analysis

Measure of Synchronization Anurak Thungtong January 26, 2009.

1

Introduction

based on Hilbert and Wavelet transforms, generalized synchronization based on nonlinear interdependency and recurrence analysis, and event synchronization. These techniques are evaluated using coupled H´enon, R¨ossler, and Lorenz systems, of which the coupling strength can be controlled. Synchronization analysis has found numerous applications in physiological signal analysis such as cardiorespiratory interactions [2] and [8], and EEG epilepsy analysis [9], [10]. However, the application of synchronization analysis in sleep state classification of neonatal EEG [11] has not yet been reported. There are two sleep states recognized in neonates: active and quiet. During quiet sleep, we would expect stronger synchronization between EEG channels due to the fact that neurons perform high synchrony over the brain area [12]. In this study, we applied various measures of synchronization to the real neonatal EEG recordings to evaluate this hypothesis. The remainder of this report is organized as follows. The background of techniques to measure the intensity of synchronization are presented in section 2. The methodology for evaluating these measures is outlined in section 3. Results of measures for synthetic coupled systems and real EEG signal are presented in section 4. Finally, discussion and conclusions are given in section 5.

Synchronization phenomenon are abundant in science, nature, engineering and social systems. Examples of systems that can operate in synchrony include clocks, singing crickets, cardiac pacemakers, firing neurons and applauding audiences. These phenomena are universal and can be understood within a common framework based on modern nonlinear dynamics [1]. Synchronization of nonlinear periodic systems was first discovered by Huygens in 1673, who noticed that two pendulum clocks that hang on the same beam can synchronize. Following the discovery of the theory of deterministic chaos, the notion of synchronization was extended and has been studied extensively in the case of interacting chaotic oscillators in the early 1990’s [2]. There are several particular cases of synchrony between two chaotic oscillators. Two coupled chaotic oscillators are said to be in generalized synchrony when there is a functional relationship connecting their dynamical variables [3]. Complete synchronization can be attained if identical systems are coupled sufficiently so that their states coincided after transients have died out [4]. Phase synchronization is the situation where two coupled chaotic oscillators keep their phases in step with each other while their amplitudes remain uncorrelated [5]. One can formulate two main problems in synchronization analysis [6]. The first problem is to reveal whether the systems under investigation are coupled and to quantify the intensity of interaction, while the second one is to characterize the driver-response relationships, or directionality of coupling. A number of methods have been proposed to measure the intensity of synchronization and specify the direction of coupling, and several methods derived from different principles and concepts are reviewed in [7]. In this study, we explored a variety of synchronization measurement techniques that have been proposed in literature; selected techniques include linear cross correlation, mutual information, phase synchronization

2

Background

In the following xn and yn (n = 1, ..., N ) will denote two simultaneously measured discrete time series of length N from two different systems x and y respectively.

2.1

Linear cross correlation

Cross correlation (CC) is a linear measure that quantifies the linear dependence between two time series [13]. The cross correlation at time delay d ∈ [−D, D] is de-

1

fined as ½ C xy (d) =

1 N −d

PN −d

x ˆn+d yˆn , C xy (−d),

n=1

d≥0 d<0

where DTFT is the Discrete-Time Fourier Transform, and DTFT−1 is the inverse Discrete-Time Fourier Transform. Similarly, hyn is defined from yn .

(1)

2.3.2

where x ˆ denotes a normalization of the variable x, i.e., x ˆ = (x − µx )/σx , µx and σx denotes mean and standard diviation of x, similarly for yˆ. When d = 0, cross correlation is precisely the correlation coefficient which has a maximum value equal to ±1 if x and y are are linearly dependent, and minimum value equal to 0. CC

2.2

= max |C

xy

d

(d)|

The second method used to extract the phase is based on the complex Morlet wavelet transform [19] defined by a positive bandwidth fb , a wavelet center frequency fc , on a K-point regular grid in the interval [lb, ub]. ϕk =

(2)

Derived from information theory [14], mutual information (M I) measures the information that two random variables x and y share: it measures how knowledge of one reduces the uncertainty about the other. Mutual information can be expressed as M X M X i=1 j=1

pxy i,j log2

pxy i,j pxi pyj

fb 2π

In the same way, 2.3.3

k = 1, ..., K

(5)

y wm

m = 1, ..., N + K − 1

(6)

is defined from yn .

Computing phase synchronization index

A couple methods for computing interaction between phases of two identical systems θnx and θny has been widely used. These methods are mean phase coherence and mutual information between phase. In this study, we used mean phase coherence method which based on the statistic of circular data [18]. Let θnxy = θnx − θny denotes the phase different of two time series, the index of phase synchronization is defined as 1 − CV where CV denotes the Circular Variance, i.e., và !2 à !2 u X X 1u PS = t cos(θnxy ) (7) sin(θnxy ) + N n n

(3)

phase synchronization

Similar to cross correlation between two time series x and y, phase synchronization measures the correlation between the phases of two time series signals. There are two steps in computing phase synchronization: extracting phase and computing synchronization index. We studied two ways to assign a phase to a univariate signal: Hilbert (P SH ) and Wavelet (P SW ) transforms. 2.3.1

−k2 2

e 2fb ej2πfc k ,

x wm = xn ⊗ ϕk ,

where M categories were used to partition the data, pxi is the marginal probability of x, and pxy i,j is the joint probability of x and y. If x and y are independent then px py = pxy so M I = 0. If they are identical then px = py = pxy so M I > 0.

2.3

1 √

The convolution of the time series with the Morlet wavelet yields a complex time series of wavelet coefficients

Mutual information

MI = −

Wavelet phase transform

It can be seen that if θnxy are identical for all n (phase lock), then P S = 1. Defining the index this way, however, the zero value can not be obtained even in the case of completely uncorrelated phases. To account for this bias, the index of phase synchronization is compared to the index obtained for the maximum possible uncorrelated phases. This value can easily be estimated by shuffle the phases of one signal before the calculation of the synchronization index. To yield a reliable statistics, 10 different shuffled signals are generated, and the average over the resulting indices is used. A normalized index of phase synchronization is thus obtained using the following deffinition: ( 0 if P S < P S S (8) P S∗ = P S−P S S else 1−P S S

Hilbert phase transform

The first approach is to use the Hilbert transform to determine the phase of a time series based from the analytic signal [15]. Given a time series xn , the analytic x signal is defined as hxn = xn + jx0n = Ax ejθn , where A and θnx denote amplitude and phase of time series x and x0n is the discrete Hilbert transform [16], defined as · · ¸ ¸ 1 · DTFT [xn ] (4) x0n = DTFT−1 DTFT πxn

where P S S is the mean phase index of shuffled signals. 2

2.4

Generalized synchronization

in generalized synchronization, we expect their respective recurrences to occur simultaneuously, and hence RRxy = RRx = RRy , where RRx (RRy ) denotes the average probability of recurrence over time in the state space of x (y), defined as

Several statistical measures have been introduced for the detection of GS such as recurrence analysis [2], the method of mutual false nearest neighbors [21],[3], synchronization likelyhood [20] and [21], and mutual predictability to detect dynamical interdependence [22]. In this study, we explored two methods which are based on the analysis of recurrence and nonlinear interdependence. Both methods have the same basic idea, that is, any two states of x that are close correspond to two states of y that are close. Both algorithms are based on nonlinear time series analysis where the nonlinear system is described by a phase space representation using a time-delay embedˆ , with ding technique. Here xn and yn , n = 1, ..., N ˆ = N − (m − 1)τ the state vector with time-delay τ N and embedding dimension m obtained from False Nearest Neighbor method [23]. 2.4.1

RRx =

By this construction, through proper normalization 0 < R ≤ 1 and low values indicate independence between x and y and higher values indicate higher levels of (generalized) synchronization. 2.4.2

Recurrence of the trajectory of a dynamical system is the return of the trajectory in state space to a neighbourhood of a point where it has been before [16]. Analysis of recurrence involves computing Euclidean distances between all points in the state space. The reˆ is defined as currence at point xi , i = 1, ..., N (9)

rnx =

ˆ (N ˆ −1) N . 2

RR

K 1 X (xn − xin,k )2 K

(13)

k=1

The criterion for the detection of generalized synchronization by means of recurrence has been proposed in [2], where an index that quantifies the degree of similarity between the respective recurrences of both systems is defined. This index is based on the average probability of joint recurrence over time, and is given by ˆ −1 N ˆ N 1 X X x = 2 r · ry Nij i=1 j=i+1 i,j i,j

Nonlinear interdependency

The algorithm for the detection of generalized synchronization by means of nonlinear interdependency has been proposed in [22]. This index of synchronization is based on the mean Euclidean distance of some nearest neighbour of points xn and yn . The algorithm is summarized as follows: Define in,k and jn,k as the time indices of the kth neighbor of xn and yn , respectively. For each xn the squared mean Euclidean distance from xn to xin,k where in,k is the time indices of the k-th neighbor of xn is defined as

where ² is a pre-defined threshold, Θ(x) is the Heaviside function, i.e., Θ(x) = 0 if x ≤ 0 and Θ(x) = 1 if x > 0. Hence, rij = 1 means xi and xj are neighbors, rij = 0 implies that xi and xj are far away from each other. Note that the total number of possible pairs is Nij =

xy

(11)

Hence, the synchronization index (a minor modification from that defined in [2]) is µ ¶ 1 RRxy RRxy GSR = + (12) 2 RRx RRy

Recurrences

¡ ¢ x ri,j = Θ ²− k xi − xj k22 = 1

ˆ −1 N ˆ N 1 X X x r Nij i=1 j=i+1 i,j

Similarily, the squared mean Euclidean distance from xn to xjn,k where jn,k is the time indices of the k-th neighbor of yn is given by rnxy =

K 1 X (xn − xjn,k )2 K

(14)

k=1

If the systems are strongly synchronized, then rnx ≈ rnxy , while rnx ¿ rnxy is always true for independent systems. Accordingly, the nonlinear interdependence measure is defined as N 1 X rnx (15) sxy = N n=1 rnxy

(10)

If both systems x and y are independent, then the average probability of a joint recurrence is expected to be zero or very close to zero. In contrast, if they are

3

2.6

Similar to analysis of recurrence, the index of synchronization is defined as sxy + syx GSS = (16) 2

Proposed in [26] and discussed in [7], event synchronization (ES) is based on the relative timings of certain events extracted from the time series. This algorithm quantifies the level of synchronicity from the number of simultaneous appearances of these events. In this study, events were defined as local maxima and minima in the time series. The data point xn is a local maxima if xn > xn±1 , and xn is a local minima if xn < xn±1 . The local maxima and minima were determined simultaneously by detecting the change of sign at every points of the time series.

where syx is obtained by changing x to y in equation (13) - (14).

2.5

Correlation dimension

The measure of synchronization based on quantifying the correlation dimension (active degrees of freedom) of the coupled systems is proposed in [24]. The measure is defined as the ratio of the sum of the correlation dimensions of the subsystems to the correlation dimension of the coupled dynamical system determined from a concatenated embedding vector obtained from the two subsystems, i.e., CD =

vx + vy vφ

sxn = sign(xn − xn−1 ) + sign(xn − xn+1 )

(17)

exn,max exn,min

ˆ −1 N ˆ N X X ¡ ¢ 2 Θ ² − kxi − xj k22 (18) ˆ (N ˆ − 1) N

ˆ → In the limit with an infinite amount of data (N ∞) and for small ², C is expected to scale like a power law, C(²) ∝ ²v , and we can define the correlation dimension v x of time series x by

vx

=

ˆ , ²) ∂ ln C x (N , ∂ ln ² ˆ , ²) lim lim d(N

²→0 N ˆ →∞

= =

sign(sxn − 2) + 1 1 − sign(sxn + 2)

(21) (22)

To quantify the relative timing of all K events, the number of times that event k appears in time series y after it appears in time series x within a time lag l is counted by

i=1 j=i+1

ˆ , ²) = d(N

(20)

Where sign(.) is the signum function, i.e., sign(x) is equal to 1 if x > 0, -1 if x < 0, and 0 if x = 0. From the condition of maxima and minima, sxn = 2 if and only if xn is a local maxima, and sxn = −2 if and only if xn is a local minima. Let exn,k = 1 if event k = 1, ..., K occurs at point n of the time series x, exn,k = 0 otherwise. Therefore, the occurrence of maxima and minima at point n of the time series x can be determined using the following equations.

where φi = [xi |yi ]. The algorithm that is widely used to calculate correlation dimension is summarized as follow [23]. The correlation sum for a collection of points xi is the fraction of all possible pairs of points which are closer than a given distance ². C x (²) =

Event synchronization

cxy = √

K N −1 X n+l X X

1 Nex Ney

exn,k eym,k

(23)

k=1 n=1 m=n

where Nex denotes number of occurance of all events found in time series x, i.e.,

(19)

N

It is obvious that the two limits defined in this formula can lead to difficulties in the analysis of practiˆ is limited by the sample size cal time series because N of the data, and the range of meaningful choices for ² is limited by the inevitable lack of near neighbours at small length scales. Several algorithms for estimating correlation dimension have been proposed in the literature [25]. In this paper, we use Grassberger-Proccacia estimator, in which v x is the slope of the plot of log C(²) versus log ². The slope is estimated based on the least-squares linear regression.

ex

=

K N −1 X X

exn,k

(24)

k=1 n=1

Similarly, cyx and Ney are defined accordingly, and we obtain the ES measure cxy + cyx (25) ES = 2 This measure is normalized to 0 ≤ ES ≤ 1 with ES = 1 if and only if all events are synchronous. Similar to index of phase synchronization, we use a renormalization which accounts for simultaneous events expected just due to chance (for example, simultaneous 4

local maxima and minima can be a considerable number in high frequency signal even for independent time series.) The normalized index is obtained by: ( 0 if ES < ES S ES ∗ = (26) ES−ES S else 1−ES S

3

H´enon system is generated by iteration. All of the coupled systems have the data length of 4000 points. The computation associated with the coupled H´enon and R¨ossler systems uses the x1 and y1 time series, while the computations associated with the coupled Lorenz system uses the x3 and y3 time series. For each system, fourty one different sets of data were generated corresponding to a different level of coupling strengh C. The coupling strength C was varied from 0 to 0.8 in steps of 0.02 for the coupled H´enon system, while C was varied from 0 to 2 in steps of 0.05 for the coupled R¨ossler and Lorenz systems. For each data set corresponding to each coupling strength, we measure the index of synchronization using the techniques mentioned in section 2. For the application to neonatal sleep EEG, we used EEG data from a study done in the Children’s Hospital of Pittsburgh in 1980, collected by Dr. Mark Scher, a pediatric neurologist at University Hospitals Case Medical Center. We randomly selected EEG recordings from four full-term and four preterm babies. For this data set, approximately the first three hours of each twelve hour recording were visually scored using oneminute epochs by Dr. Scher. Each data set contains fourteen channels of EEG, based on the standard 10-20 EEG electrode placement scheme (illustrated in figure 7). For this analysis we selected only four EEC channels in the central lobe area: C3-Cz, Cz-C4, Fz-Cz, and Cz-Pz. The selected EEG channels were preprocessed to eliminate artifacts and then segmented into multiple epochs with one minute window length. For each epoch, we measured synchronization across all four channels. The parameters for the measures to the different coupled model systems are summarized in table 1.

Materials and Methods

We evaluate these techniques by three coupled dynamical model systems: the coupled H´enon, R¨ossler, and Lorenz systems which include unidirectional coupling of either two H´enon, two R¨ossler, or two Lorenz systems with subsystems x and y. The coupling strength is controled by a system parameter C. The coupled systems are defined as follows. The coupled H´enon system is given by x1,n+1 x2,n+1 y1,n+1 y2,n+1

= 1.4 − x21,n + 0.3x2,n = x1,n ¡ ¢ 2 = 1.4 − Cx1,n y1,n + (1 − C) y1,n + 0.3y2,n = y1,n (27)

The coupled R¨ossler systems is given by x˙1 x˙2 x˙3 y˙1 y˙2 y˙3

= −0.95x2 − x3 = 0.95x1 + 0.15x2 = 0.2 + x3 (x1 − 10) = −1.05y2 − y3 + C(x1 − y1 ) = 1.05y1 + 0.15y2 = 0.2 + y3 (y1 − 10)

(28)

Parameters D fb fc K [lb ub] m t ²

The coupled Lorenz systems is given by x˙1 x˙2 x˙3 y˙1 y˙2 y˙3

= 10(x2 − x1 ) = x1 (28 − x3 ) − x2 8 = x1 x2 − x3 3 = 10(y2 − y1 ) = y1 (28.001 − y3 ) − y2 8 = y1 y2 − y3 + C(x3 − y3 ) 3

H´enon 1 2 6 16 [-4 4] 4 1 0.02

R¨ossler 0 2 1/2π 16 [-4 4] 4 1 0.02

Lorenz 5 2 1/2π 16 [-4 4] 4 1 0.02

EEG 20 0.002 0.1π 4 [−π π] N/A N/A N/A

Table 1: Parameters for the measures to the different coupled

(29)

model system: maximum time delay D, bandwidth fb , center frequency fc , grid K, lower bound lb, upper bound ub, embedding dimension m, time delay t, and threshold ².

The coupled R¨ossler and Lorenz systems are generated by integrating equations (28) and (29) using a 4th order Runge-Kutta integration scheme with a step size of 0.001 and a sampling interval of 0.01. The coupled 5

4

Results

4.1

show rather steep increases for coupling strength in the range 0.6 − 0.8.

Coupled systems

2

2

2

0

0

0

Rossler

−2

0

10 20 C=0

−2

0

10 20 C = 0.6

−2

5

5

0

0

0

0

−5

50

Lorenz

50

50 100 C = 0.06

0

50 100 C = 0.35

−5

10 20 C = 0.8

0

0.4 0.2

0

0.2

0.4 0.6 Coupled strength

0.8

Figure 2: The coupling strength for eight measures of syn-

0

chronization applied to H´ enon systems: linear cross correlation (CC), mutual information (M I), phase synchronization based on Hilbert (P SH ) and Wavelet (P SW ) transforms, event synchronization (ES), generalized synchronization based on nonlinear interdependency (GSS ) and recurrences (GSR ), and correlation dimension (CD).

50 100 C=2

40

The results for R¨ossler systems are shown in figure 3. All measures consistently show monotonically increasing values for higher coupling strength. Correlation dimension has high values for weak coupling as in the H´enon system. Event synchronization starts to increase at intermediate coupling strength only. Mutual information, recurrence, and nonlinear interdepency gradually increase, cross correlation shows steep increases for coupling strength of 0-0.25 and stays quite unchanged after coupling strength is greater than 0.5. The results for the Lorenz systems are presented in figure 4. All measures also consistently exhibit increasing values for higher coupling strengths and identical coupling at a coupling strength of 1.5. Mutual information, event synchronization, nonlinear interdependency, and recurrence have almost identical results. They slightly increase in the uncoupled to intermediate coupling range, and then abruptly increase at a coupling strength of 1.5. Linear cross correlation and synchronization based on correlation dimension start with slightly higher values and gradually increase for higher coupling strengths. However, larger fluctuations can be seen in linear cross correlation and phase synchronization. To compare the different measures of synchronization in terms of their capability to distinguish between

0 0 200 400 0 200 400 0 C=0 C = 0.85

200 400 C=2

Figure 1: Short segments of H´enon, R¨ossler, and Lorenz time series x (solid), and y (dash) at different coupling strengths C.

4.2

0.6

−0.2

20 0

0.8

0

0

5

−5

CC PSH PSW MI ES CD GSS GSR

1 Degree of synchronization

Henon

Figure 1 illustrates short segments of H´enon, R¨ossler, and Lorenz time series x and y at different coupling strengths. It can be seen that the time series exhibit weak synchronization for a small value of C and complete synchronization at the maximum value of C.

Synchronization in coupled systems

Figure 2 shows eight measures of synchronization and their relationship to the coupling strength between two identical H´enon systems. It can be seen that all measures show increasing degrees of synchronization with increasing coupling strength and most of them reflect identical synchronization at a coupling strength of about 0.7. For most measures, the synchronization measure increases monotonically with coupling strength with a more shallow slope in the range 0 − 0.6 and a sharper slope in the range 0.6 − 0.8. The correlation dimension (CD) and nonlinear interdependency (GSS ) have larger fluctuations with coupling strength in the range 0 − 0.6 than can be seen in some of the other measures. Mutual information (M I) and analysis of recurrence (GSR )

6

1 Degree of synchronization

Measures CC P SH P SW MI GSS GSR CD ES Average

CC PSH PSW MI ES CD GSS GSR

0.8 0.6 0.4 0.2

H´enon 0.7463 0.9232 0.9439 0.9512 0.8463 0.9610 0.8000 0.9256 0.8872

R¨ossler 0.9951 0.9037 0.9268 0.9732 0.9780 0.9927 0.6927 0.7659 0.9035

Lorenz 0.7061 0.7049 0.7573 0.7451 0.6902 0.7012 0.7683 0.6720 0.7174

Average 0.8159 0.8439 0.8760 0.8898 0.8390 0.8821 0.7537 0.7878 0.8360

0

Table 2: Degree of monotony of all measures: linear cross cor−0.2

0

0.5

1 1.5 Coupled strength

relation (CC), mutual information (M I), phase synchronization based on Hilbert (P SH ) and Wavelet (P SW ) transforms, event synchronization (ES), nonlinear interdependency (GSS ), phase synchronization by means of recurrences (GSR ), and synchronization based on correlation dimension (CD). Bold text denotes maximum values for each column.

2

Figure 3: Same as figure 2 but applied for R¨ossler system.

1 Degree of synchronization

4.3

CC PSH PSW MI ES CD GSS GSR

0.8

All techniques are implemented to measure synchronization in the EEG; however, only two methods: linear cross correlation and Wavelet phase synchronization provided acceptable results. Figure 5 presents box plots of the synchronization index measured from each of the selected electrode pairs. The results consistently show that the degree of synchronization during quiet sleep seems to be higher than that during active sleep for both measures. Figure 6 summarizes the results from four full-term and four preterm neonates. Note that we applied Principal Component Analysis (PCA) [29] to reduce the dimension of data from six to one so that it is more convenient to present the results. Again, the statistics roughly show that the (transformed) synchronization index in active and quiet sleep is separable in all subjects. This results suggest that using an index of synchronization obtained using either technique as features for sleep state classification in neonatal EEG is possible. However, notice that the degree of synchronization across subjects may not be consistent. Therefore, synchronization analysis may only be suitable for offline classification rather that online (real-time) classification.

0.6 0.4 0.2 0 −0.2

0

0.5

1 1.5 Coupled strength

Synchronization in neonatal EEG

2

Figure 4: Same as figure 2 but applied for lorenz system. different degrees of coupling, we measured degree of monotonicity as proposed in [7]. The results are presented in table 2. From the table, generalized synchronization based on recurrence, linear cross correlation, and phase synchronization based on Hilbert transform provide the highest values of monotonicity for H´enon, R¨ossler and Lorenz systems respectively. The maximum average values for all systems are obtained by Mutual Information (0.8898), Recurrence (0.8821), and Wavelet phase (0.8760). Correlation dimension provides the lowest average value at 0.7537. Regarding the different model systems, the R¨ossler and H´enon systems have the highest average values (0.9035 and 0.8872), the Lorenz system has the smallest average value at 0.7178 due to the fluctuations of the measures as seen in figure 4.

5

Discussion and Conclusions

Seven techniques of synchronization measurement: linear cross correlation, mutual information, phase synchronization based on Hilbert and Wavelet transforms, generalized synchronization based on nonlinear interdependency and recurrence analysis, and event synchro7

Cross correlation

A−B

A−C

A−D

B−C

B−D

C−D

0 1

0 1

0 1

0 1

0 1

0 1

1.5 1 0.5 0

Wavelet

0.9 0.8 0.7 0.6

Figure 5: Statistical Summary (box plots) of synchronization

Figure 7:

The 10-20 EEG electrode placement system widely used in neonatal EEG. This image was obtained from (http://www.immrama.org/eeg/electrode.html).

in active (0), and quiet (1) sleep measured from each pair of EEG channels: A=C3-Cz, B=Cz-C4, C=Fz-Cz, and D=Cz-Pz. The upper figures are cross correlation and the lower figures are Wavelet phase. These results were obtained from a full-term neonate.

Wavelet

Cross correlation

signal. They are able to distinguish different degrees of synchronization during active and quiet sleep in all EEG recording of selected neonatal subjects. All techniques based on nonlinear analysis (GSR , GSS , and CD) were not easily applied to real neonatal EEG because of the difficulty in selecting parameters: embedding dimension m and time-delay τ for phase space reconstruction. These parameters are relatively difficult to determine because of the complexity of EEG signal. Moreover, these technoques have very high computational complexity. We have discussed about measuring degree of synchronization in the identical systems. In subsequent work we will modify these techniques, as well as explore some new techniques to two different systems such as the electrocardiogram and respiratory signals of a sleeping human suffering from sleep disordered breathing, hypertension, etc. An estimation of directionality and causality in coupling is also another problem of interest in synchronization analysis rather than just measuring the strength of synchronization. The knowledge of directionality of coupling is very helpful in better understand the interaction of between systems. Several algorithms proposed for determining the directionality of coupling can be found in [2], [13], [24], [26], and [28].

0 1

0 1

0 1

0 1

0 1

0 1

0 1

0 1

Figure 6: Statistical Summary (box plots) of transformed synchronization data in active (0), and quiet (1) sleep of four fullterm and four preterm neonates. The upper figures are cross correlation and the lower figures are Wavelet phase. nization are implemented and evaluated with coupled H´enon, R¨ossler and Lorenz systems, as well as EEG signals obtained from neonates during sleep. In the synthetic model systems, all techniques have satisfactory ability to distinguish the intensity of synchronization. They exhibit slightly different performance in the different systems, but all measures show some degree of monotonicity with coupling strength. Among these techniques, mutual information has the highest averaged degree of monotonicity, while correlation dimension exhibits the lowest. Measurement of phase synchronization based on the Wavelet transform shows a higher degree of monotonicity than the Hilbert transform. Quantifying generalized synchronization using recurrence analysis performs better than nonlinear interdependency as related to monotonicity. Linear cross correlation and Wavelet phase synchronization work very well in applications of real EEG

References [1] A. Pikovsky, M. Rosenblum, and J. Kurths, Synchronization : A Universal Concept in Nonlinear Sciences, Cambridge University Press 2003. [2] Maria Carmen Romano, Marco Thiel, Jurgen Kurths, Martin Rolfs, Ralf Engbert, and Reinhold Kliegl. Synchronization Analysis and Recurrence

8

in Complex Systems. Handbook of Time Series Analysis pp. 243-276. Wiley-VCH Verlag Berlin GmbH (2006).

the Presence of Noise and Nonlinearity. Handbook of Time Series Analysis pp. 265-282. Wiley-VCH Verlag Berlin GmbH (2006).

[3] L. Kocarev and U. Parlitz, Generalized Synchronization, Predictability, and Equivalence of Unidirectionally Coupled Dynamical Systems. Phys. Rev. Lett. 76, 1816 - 1819 (1996)

[14] R. Gray. Entropy and Information Theory. Springer Verlag, New York, 1990.

[4] H. Fujisaka and T. Yamada, Stability theory of synchronized motion in coupled-oscillator systems, Progr Theor Phys 69 (1983), pp. 3247.

[16] Marple, S.L., Computing the discrete-time analytic signal via FFT. IEEE Transactions on Signal Processing, Vol. 47, No.9 (September 1999), pp.26002603.

[5] Michael G. Rosenblum, Arkady S. Pikovsky, and Jurgen Kurths, Phase Synchronization of Chaotic Oscillators. Phys. Rev. Lett. 76, 1804 - 1807 (1996)

[17] Sanjit K. Mitra, Digital signal processing: a computer-based approach . McGraw-Hill (2002)

[15] D. Gabor. Proc. IEEE London 93 (1946) 429.

[18] N.I. Fisher, Statistical analysis of circular data. Cambridge University Press, London, 1993.

[6] Michael Rosenblum, Laura Cimponeriu, and Arkady Pikovsky. [?]. Handbook of Time Series Analysis pp. 159-180. Wiley-VCH Verlag Berlin GmbH (2006).

[19] Xiaoli Li, Xin Yao, John Jefferys, and John Fox. Computational Neuronal Oscillation with Morlet Wavelet Transform. 27th Annual International Conference of the IEEE.

[7] Kreuz, T and Mormann, F and Andrzejak, RG and Kraskov, A and Lehenrtz, K and Grassberger, P. Measuring synchronization in coupled model systems: A comparison of different approaches. Physica D, 225 . pp. 29-42 (2007).

[20] C.J. Stam, B.W. van Dijk, Synchronization likelihood: An unbiased measure of generalized synchronization in multivariate data sets. Physica D 163 (2002) 236251.

[8] Ralf Mrowka, Andreas Patzak, Michael Rosenblum, quantitative analysis of cardiorespiratory synchronization in infants. Int.J. of. Bifurcation and Chaos, Vol. 10, 2000, pp.2479-. 2488. DSP 2002 -1144.

[21] R. Quian Quiroga, J. Arnhold, and P. Grassberger. Learning driver-response relationships from synchronization patterns. Physical review E V61, Number 5 May 2000.

[9] Altenburg J, Vermeulen RJ, Strijers RL, Fetter WP, Stam CJ. Seizure detection in the neonatal EEG with synchronization likelihood. Clin Neurophysiol. 2003 Jan;114(1):50-5.

[22] J. Arnhold,P. Grassberger, K. Lehnertz and C. E. Elger, A robust method for detecting interdependences: application to intracranially recorded EEG. Physica D: Nonlinear Phenomena, Volume 134, Issue 4, 10 December 1999, Pages 419-430.

[10] A. J. C. Slooter, E. M. Vriens, F. S. S. Leijten, J. J. Spijkstra, A. R. J. Girbes, A. C. van Huffelen and C. J. Stam Seizure detection in adult ICU patients based on changes in EEG synchronization likelihood. Neurocritical Care, Volume 5, Number 3 / December, 2006.

[23] Holger Kantz and Thomas Schriber. Nonlinear Time Series Analysis. Second Edition, Cambridge University Press, 2003. [24] S. Janjarasjitta, and K.A. Loparo. An approach for characterizing coupling in dynamical systems. Physica D: Nonlinear Phenomena Volume 237, Issue 19, 1 October 2008, Pages 2482-2486

[11] Anurak Thungtong, Real-Time Automated Sleep Scoring of Neonates. MS Thesis. Case Western Reserve University.

[25] Review of correlation dimension estimators, http://dissertations.ub.rug.nl/FILES/faculties/ science/1998/s.a.borovkova/c2.pdf

[12] Buzski, Gyrgy. 2006. Rhythms of the brain. Oxford University press. [13] Theoden I. Neto., Thomas L. Carroll, Louis M. Pecora, and Steven J. Schi. Detecting Coupling in 9

[26] Smirnov, Dmitry A.; Bezruchko, Boris P. Estimation of interaction strength and direction from short and noisy time series. Physical Review E, vol. 68, Issue 4, id. 046209 10/2003. [27] C.W. J. Granger, Econometrica 37,424 (1969). [28] Michael G. Rosenblum, Identification of coupling direction: Application to cardiorespiratory interaction, Physical review E, V65 041990. [29] Smith, Lindsey. 2002. A tutorial on Principal Component Analysis. Febuary 26, http://kybele.psych.cornell.edu/ edelman/Psych465-Spring-2003/PCA-tutorial.pdf.

10