Andreas Groth

Visualization of coupling in time series by order recurrence plots

Keywords: numbers: 05.45.Tp, 05.45.Xt, 87.19.La

published or not. The documents may come

Visualization of coupling in time series by order recurrence plots

I. INTRODUCTION

In the analysis of coupled systems various techniques have been developed to detect cooperative behavior from observed time series [START_REF] Brockwell | Time series : theory and methods[END_REF][START_REF] Kantz | Nonlinear Time Series Analysis[END_REF][START_REF] Pikovsky | Synchronization[END_REF]. Depending on the nature of the systems, there are different requirements to the methods. While linear methods based on correlations are not sufficient to deal with nonlinear dependencies, most nonlinear methods require sufficiently long stationary time series. For the case that stationarity holds only for short observation time, cross recurrence plots (CRPs) were introduced [START_REF] Zbilut | [END_REF]5]. However, this method is based on taking distances of trajectories, which is conceptually difficult on physically different systems. A general problem in studying multivariate data from natural systems, for instance EEG data, is that measurement conditions change with time. Among others offset and amplitude range can vary differently within the channels (cf. Fig. 4).

To overcome this problem we consider a special symbolic dynamics of the system, where the time series is encoded by order patterns. This yields to new symbol sequences, which are invariant with respect to certain distortions in amplitude. With this symbolic dynamics a complexity measure was already proposed [6] and successfully applied to epileptic seizure detection [7]. Moreover, a distance between time series was introduced to study similarities and dissimilarities between EEG channels [8]. Following the idea of CRPs we introduce a new visualization tool based on the recurrence of order patterns. Similar to the CRPs we obtain different structures such as dots and lines, from which conclusions to the temporal behavior of e.g. coupling strength and delay can be drawn. In this context a sensitivity of the methodology to phase coupling is shown. This paper is organized as follows. First we define the symbolic dynamics and introduce order recurrence plots. On coupled oscillators we demonstrate the ability to visualize and quantify the coupling strength. Finally, we apply the technique to EEG data. * Electronic address: groth@uni-greifswald.de

x

π x = 0 π x = 1 π x = 2 π x = 3 π x = 4 π x = 5

II. ORDER RECURRENCE PLOTS AND THEIR QUANTIFICATIONS

A. Symbolic dynamics Given a dynamical system represented by a onedimensional time series {x(t)} t the original phase space can be reconstructed by time delay embedding x(t) = x(t), x(t + ϑ), . . . , x(t + (D -1)ϑ) for sufficiently high dimension D [START_REF] Takens | Lecture Notes in Mathematics[END_REF]. The time delay ϑ can be chosen in an appropriate way. In this paper we restrict ourselves to D ≤ 3. Higher dimensions were also studied, but D = 3 gives already convincing results.

To discuss the order structure of the time series, let us start with D = 2. Two relations between x(t) and x(t+ϑ) are possible, neglecting equality [26]. We denote these relations as order patterns π and derive a new symbol sequence

π x (t) = 0 : x(t) < x(t + ϑ) 1 : x(t) > x(t + ϑ) . (1)
This encoding of a trajectory yields to a decomposition of the phase space into two equal areas.

In dimension D = 3 there are D! = 6 order patterns possible (Fig. 1), again neglecting equality. The phase space is decomposed by three planes into six equivalent regions (Fig. 2). These planes are defined by

x(t) = x(t + ϑ), x(t) = x(t + 2 ϑ) and x(t + ϑ) = x(t + 2 ϑ).
The new symbol sequences are invariant with respect to an arbitrary, increasing transformation of the amplitude. Moreover, we are able to compare symbol sequences from physically different systems, which is an important advantage in contrast to CRPs. Later on we

-1 0 1 -1 0 1 -1 1 x(t) x(t+ϑ) x(t+2ϑ) -1 0 1 -1 0 1 -1 1 π x =1 x(t) π x =0 π x =3 π x =5
x(t+ϑ) will also demonstrate robustness against low-frequency trends and strong observational noise.

π x =2 π x =4 x(t+2ϑ)

B. Order recurrence plots

To visualize the dynamics of phase space trajectories, recurrence plots were introduced [START_REF] Eckmann | [END_REF]. This idea was further extended to CRPs to study recurrent dynamics between different systems. Two trajectories are considered as neighboring if ε ≤ x(t) -y(t ′ ) is fulfilled, with ε a cut-off distance.

We introduce a similar technique, where the neighborhood is defined by the local order structure. Phase space vectors x(t) and y(t ′ ), of the same or different dynamical systems, are defined as neighboring if their order patterns coincide: π x (t) = π y (t ′ ). With this definition we introduce the order recurrence plot (oREP)

R(t, t ′ ) = 1 : π x (t) = π y (t ′ ) 0 : otherwise , (2)
a plot of a boolean matrix of size M × N , where M, N are the lengths of the order pattern sequences.

Similarly to the CRPs we observe different structures in the oREPs. If there are no dependencies the plot is dominated by single dots (Fig. 3a). Strong dependencies yield to diagonal lines, where R(t + k, t ′ + k) = 1 for 0 ≤ k ≤ l with l the length of the line. Here the order pattern sequences π x and π y coincide at a certain period of time. But in contrast to a CRP, the diagonal lines are also preserved on a time-dependent distortion of the amplitudes (Fig. 3b). On time series with different time scales we observe curved lines (Fig. 3c), which indicate the lines of synchronization [11]. In practical time series we normally obtain a mixture of different structures (Fig. 3d).

What we get from all these plots is a global impression of the dynamics of the underlying system. Dependencies over short as well as long distances are obtained, but often we are only interested in short-time dependencies.

Therefore, we slightly rewrite (2) in analogy to the close returns plot [12]

R(t, τ ) = 1 : π x (t) = π y (t + τ ) 0 : otherwise , (3)
where we focus the study to an area close to the main diagonal. Diagonal lines R(t + k, t ′ + k) = 1 are transformed to horizontal lines R(t + k, τ ) = 1 and it is more convenient to study a longer range in time (Fig. 4).

C. Optimal embedding parameter

To reconstruct the dynamics by time delay embedding x(t) we have to find a proper ϑ. As an idealization of systems with periodic behavior, we study a simple harmonic oscillator, where the trajectory describes an ellipse in phase space. We have x(t) = A sin(2πt/T ), where A is the amplitude, and T is the period length.

At D = 2 the trajectory crosses every half period the main diagonal and we get a piecewise constant π x (t). Independently of ϑ we have prob{x(t) < x(t + ϑ)} = prob{x(t) < x(t + ϑ)} = 1/2 and the oREP shows a checkerboard. With ϑ we control just the time of crossing and with it the initial phase. However for practical reasons, especially on noisy signals, an appropriate ϑ has to be chosen to minimize the number of false neighbors.

For periodic functions we propose ϑ = T /4.

The more interesting case, where we focus on is D = 3. Now the probabilities of the order patterns depend on the time delay and with it the oREP. For small embedding delays we almost obtain upward π x (t) = 0 or downward patterns π x (t) = 5 and the oREP does not significantly differ from D = 2. But with increasing embedding delays the other order patterns become more frequent, where we get uniform distribution at ϑ = T /3. At this value the trajectory of the harmonic oscillator is maximally enfolded in this special phase space decomposition (Fig. 2b) and it describes a circle on a plane, orthogonal to the main diagonal. With this proper time delay we obtain distinct diagonal lines as in Fig. 3b, consisting of small squares of length T /6.

D. Recurrence rate

Since the introduction of the recurrence plots, several measures were defined in the recurrence quantification analysis to describe the underlying system [13,14]. We consider the simplest quantity and count the number of dots in the oREP over time t as a function of τ . By analogy we call it the recurrence rate of order patterns

RR(τ ) = t R(t, τ ) . (4)
It represents a statistical measure of similarities between two dynamical systems. In statistics we find a comparable approach, describing relations between two samples

a ij b ij = 1 over discordance a ij b ij = -1 k = i<j a ij b ij n(n -1)/2 , ( 5
)
where -1 ≤ k ≤ 1. Thus (4) can be interpreted at D = 2 as Kendall's tau for time-dependent samples. Concordance means π x (t) = π y (t + τ ) and discordance π x (t) = π y (t + τ ). In contrast to (5) we have 0 ≤ RR(τ ) ≤ 1. In [START_REF] Ferguson | [END_REF] Kendall's tau on time series was already discussed, but with a completely different intention. To define tests for the fit of models the sum over ϑ with fixed τ was analyzed. Here we consider RR(τ ) as a function of the time lag τ between the time series with fixed but adapted delay embedding parameter ϑ.

III. MEASURE OF PHASE COUPLING

To quantify the strength of coupling we analyze the distribution of the recurrence rate RR(τ ) over τ . For strong coupling we get distinct horizontal lines in the oREP, what yields to distinct maxima in RR(τ ). For weak coupling RR(τ ) tends to an uniform distribution. We normalize the recurrence rate by rr(τ ) = RR(τ )/ τ RR(τ ) and introduce a coupling index by means of Shannon entropy

ρ π = 1 - -τmax τ =τmin rr(τ ) log rr(τ ) log(τ max -τ min ) . (6)
This gives 0 ≤ ρ π ≤ 1, where ρ π = 0 corresponds to no coupling. The actual maximum of ρ π depends on [τ min , τ max ], if one or more maxima with a distance of a mean recurrence time are contained. As Fig. 2 already indicates the order patterns can be considered as a discretization of an instantaneous phase. In the appendix it is shown for D = 3 that this phase is given as

tan φ π (t) √ 3 = x(t + 2 ϑ) -x(t) x(t + 2 ϑ) -2 x(t + ϑ) + x(t) (7)
the angle between first and second order differences. A harmonic oscillator describes an ellipse in this phase space and φ π is well defined. On more complex oscillators as the Rössler system this phase space also exhibits a clear oscillation around the origin even in the funnel regime. Thus we expect a sensitivity of ρ π to phase coupling.

A common approach to introduce a phase is based on the analytic signal [17]

z(t) = x(t) + i x(t) = A(t)e iφH(t) . (8)
To define a coupling index the distribution of the instantaneous phase differences ∆φ

H (t) = (φ 1 H (t) -φ 2 H (t)
) mod 2π can be analyzed by means of the Shannon entropy [18]

ρ H = 1 -S H /S max , (9)
where S H is the entropy of the distribution of ∆φ H (t) and S max = ln m the maximum entropy on m bins. However, this classical interpretation of a phase holds only for narrow band signals. x(t+2ϑ) Figure 5: Time delay embedding of the Rössler system with ϑ = 1.9 and its decomposition by order patterns equivalent to Fig. 2 For non-phase-coherent oscillators an alternative method was introduced by the curvature of an arbitrary curve [19,20] with ψ(t) = arctan ẋ(t)/ ẏ(t). But this approach is restricted to systems where we are able to obtain at least two components. In [21] this idea was adopted to introduce a phase by the curvature of an analytic signal (8)

φ c (t) = arctan ẋ(t) ẋ(t) , (10)
where just as in ( 7) only a single component is needed.

Analogously to ρ H the coupling strength is calculated from the phase differences ∆φ c (t) = (φ 1 c (t) -φ 2 c (t)) mod 2π by means of Shannon entropy

ρ c = 1 -S c /S max , (11)
with S c the entropy of the distribution of ∆φ c (t).

A. Detection of coupling in model systems

In a first example we consider two coupled Rössler systems

ẋ1,2 = -ω 1,2 y 1,2 -z 1,2 + k (x 2,1 -x 1,2 ) ẏ1,2 = ω 1,2 x 1,2 + a y 1,2 (12) ż1,2 = b + z 1,2 (x 1,2 -c)
in the simple phase-coherent regime (a = 0.15, b = 0.2, c = 10). They are diffusely coupled with strength k and their frequencies are detuned by ω 1,2 = ω ± ∆ω. The equations are numerically solved with ∆t = 0.05. In the time delay embedding we obtain a well defined oscillation around the main diagonal (Fig. 5). Depending on the coupling strength we obtain the phenomena of phase synchronization (PS), where the phase difference is bounded, even in the case of detuned frequencies [22]. The transition from non-PS to PS can also be seen from oREPs. For weak coupling the phases are not locked and we get drifting lines (Fig. 6, upper panel). With increasing coupling strength the phases get almost locked, but the lines are interrupted by phase jumps (Fig. 6, middle panel). For sufficiently high coupling we obtain PS and the order patterns of both oscillators coincide, where we get more or less straight lines (Fig. 6, lower panel). In contrast to this CRPs visualize recurrent amplitudes. Since the phases are synchronized and not the amplitudes we obtain only short lines in CRPs. In [23,Sec. 2.6] it was demonstrated that CRPs are not appropriate to distinguish between PS and non-PS on Rössler systems. If we study the coupling on a whole range of parameters (∆ω, k) we obtain an area of phase locking, called Arnold tongue (Fig. 7). This coincides also with results from [22]. Although the oscillators are coupled by the x-components, we obtain a similar plot for the y-and z-components (not shown here). If we calculate ρ H on the z-component we have to take care, because of the different origin. However the oREP can be calculated straightforward with the same parameters as on the x-and y-components. If we add observational noise to the time series the area of phase locking is well preserved, even at a signal-to-noise-ratio of σ x /σ Noise = 1/3 (Fig. 8). Here ρ π clearly outperforms ρ H . In a second example we consider the more challenging funnel regime of the Rössler system, where we have no straightforward phase definition. The parameter a in (12) governs the topology of the attractor. Above a critical value the trajectory does no more cycle around an origin in the x-y-plane (Fig. 9, upper panel). However there is an appropriate time delay embedding of the attractor, where the trajectory cycles around the main diagonal (Fig. 9, lower panel). Again we observe the transition to PS by growing horizontal lines (Fig. 10). This coincides with [19], where PS was already shown for this system on the basis of ψ. It is remarkable that in contrast to the phase-coherent regime there are no longer horizontal lines beside τ ≈ 0 and with it no further distinct peaks in the recurrence rate. This results from the rather broad distribution of the time scale of one cycle. The irregular structure of the horizontal lines is also an indicator of different time scales, where a closer examination of Fig. 10 would show smaller and larger blocks. The size depends on the length of the cycles, wherefore all blocks are nearly equal in the coherent case. Moreover a much stronger coupling k is needed to obtain PS in the noncoherent regime.

B. An application to EEG signals

Several neurological diseases such as epilepsy or Parkinson's disease manifest in a synchronization of neuronal groups. The specific type of phase synchronization has been discussed in [24,25] to detect seizure activity. On intracranial EEG data the phenomenon of phaselocking during seizures were shown. Here, we present the concept of oREP on scalp EEG data, which are more susceptible to artifacts.

To get time-dependent coupling indices ρ π (t), ρ H (t) and ρ c (t) we study the time series in a sliding window analysis. In the first case we start with the calculation of R(t, τ ) over the whole time t. Next, the recurrence rate is calculated in consecutive overlapping frames of length L RR(t, τ ) =

t+L t ′ =t R(t ′ , τ ) , (13)
and from (6) we get ρ π (t). In the second case a meaningful phase is only given on narrow-band signals. For this reason the EEG signals are band-pass filtered. Performing (9) on the filtered signals in a sliding window of length L gives ρ H (t). In the third case the analytic signal is determined over the whole time series without any preprocessing and the derivatives are estimated by a Savitzky-Golay FIR filter. Again ρ c (t) is calculated in a sliding window. We consider EEG signals from a 14 years old child with epileptic disorder. The 19-channel scalp EEG data (international 10-20-system) were sampled with 256 Hz sampling rate and band-pass filtered (0.3 -70 Hz). In Fig. 11 the epileptic seizure onset is shown on two representative EEG channels. Before the seizure there are no dominant oscillations and the oREP shows no significant longer lines of synchronization. This behavior changes during the seizure, where the time series show clear oscillations in the alpha range (≈ 8 -13 Hz, Fig. 11h). These oscillations are synchronized, what is well determined by distinct horizontal lines in the oREP. The resulting coupling strength ρ π clearly reveals the period of the seizure. The time delay is adapted to a mean frequency in the alpha range, where further tests turn out that the coupling index is not sensitive to time delays ϑ = T /3 with corresponding frequency 1/T in the alpha range. With coupling index ρ H of the band-pass filtered signals (IIR Butterworth filter, 8 -13 Hz) we get qualitatively the same shape as with ρ π (Fig. 11f), what suggests a synchronization within the alpha range. Without filtering the period of seizure activity is not reliable revealed (not shown here). Considering ρ c pre-processing is not needed and the seizure is well detected (Fig. 11g). In the overall statistics we get an impression of the power to distinguish between seizure and no seizure (Fig. 12). Obviously ρ π is most suitable to separate both states. On ρ H and ρ c the overlaps between both states are clearly larger, what yields to larger classification errors.

IV. SUMMARY

We have introduced a method to visualize cooperative behavior between two coupled dynamical systems. The trajectories are encoded by its local ordinal structure. With this special phase space decomposition we introduce the order recurrence plots analogously to the idea of cross recurrence plots. The distribution of recurrence points provides an evaluation of similarities between both systems and we have introduced a measure of coupling strength. The potential of this approach were demonstrated on different model systems, where we have shown connections to results from phase synchronization analysis. An important advantage is the conceptually simple and computationally fast algorithm. Invariance with respect to an distortion in amplitude and robustness against trend and noise reveals the method as a helpful tool in applications. In a first short application EEG data were successfully analyzed.

Although low dimensions work fine in the future work higher dimensions are of interest. Furthermore other quantities from the recurrence quantification analysis have to be considered, such as average line length or determinism, to describe different nonlinear phenomena.

Appendix

For dimension D = 3 we show connections between an instantaneous phase and order patterns. Therefore we introduce a cylindrical coordinate system (r, φ, z) in terms of the delay embedding coordinates (x(t), x(t + ϑ), x(t + 2 ϑ)). We start with a cartesian coordinate system where one axis is the main diagonal. This gives the first unit vector e ′ 3 = (e 1 + e 2 + e 3 )/ √ 3. We complete this system taking e ′ 2 = (e 3e 1 )/ √ 2 and e ′ 1 = (e 1 -2 e 2 + e 3 )/ √ 6 with the transformation

  x ′ y ′ z ′   =   1/ √ 6 -2/ √ 6 1/ √ 6 -1/ √ 2 0 1/ √ 2 1/ √ 3 1/ √ 3 1/ √ 3     x(t) x(t + ϑ) z(t + 2 ϑ)   .
From this we simply derive the cylindric coordinates by The cylindric coordinate system is oriented with z along the main diagonal and the order pattern of the delay coordinates is completely determined by φ. Therefore the order patterns can be considered as a discretization of the instantaneous phase φ into six values.

r cos φ = x ′ = x(

Figure 1 :

 1 Figure 1: Order patterns at embedding dimension D = 3

Figure 2 :

 2 Figure 2: left: Delay embedding of a sine function and decomposition of the phase space by order patterns at dimension D = 3. right: Same plot with viewing angle in direction of the main diagonal

Figure 3 :Figure 4 :

 34 Figure 3: oREP on different time series. (a) Bivariate Gaussian noise, (b) Periodic functions with time-varying amplitudes and offsets, but constant period, (c) Periodic functions with decreasing and increasing period, respectively and (d) Onset of an epileptic seizure on EEG-channels F7 and Fz. All plots with embedding dimension D = 3. The embedding delays are appropriately chosen as discussed in Sec. II C

Figure 6 :Figure 7 :

 67 Figure 6: oREP of the x-components of two coupled Rössler systems with detuned frequencies ω1,2 = 1 ± 0.015 at different coupling strength k and corresponding recurrence rate rr(τ )

Figure 8 :

 8 Figure 8: Coupling indices on two coupled Rössler systems with observational noise σx/σNoise = 1/3. (a) ρπ and (b) ρH as a function of coupling strength k and detuning ∆ω. Coupling indices are gray-coded, where black corresponds to no coupling. Time series length L = 2 000

Figure 9 :

 9 Figure 9: Upper panel: Projection of the Funnel system onto the x-y-plane. Lower panel: Time delay embedding with D = 3 and ϑ = 0.9. Parameters are ω = 1, b = 0.1, and c = 8.5

Figure 10 :

 10 Figure 10: oREP of two coupled Funnel attractors at different coupling strength k and corresponding return times rr(τ ). Parameters are ω1,2 = 1 ± 0.015, a = 0.29, b = 0.1, and c = 8.5

Figure 11 :Figure 12 :

 1112 Figure 11: Coupling of two EEG channels (a-c) during an onset of an epileptic seizure. (d) oREP with ϑ = 2.7ms and (e) resulting coupling index ρπ. (f) ρH and (g) ρc. All coupling indices with L = 2s. (h) Spectrogram |F(t, ω)| of C3, where black indicates high values

  t) -2x(t + ϑ) + x(t + 2 ϑ) / √ 6 r sin φ = y ′ = x(t + 2 ϑ) -x(t) / √ 2 , z = z ′ . + 2 ϑ) -2 x(t + ϑ) + x(t) ,depending on first and second order differences. If the function x(t) is sufficiently smooth compared to ϑ, we can approximate x(t + 2 ϑ) -x(t) by ẋ(t + ϑ) 2 ϑ and x(t + 2 ϑ) -2 x(t + ϑ) + x(t) by ẍ(t + ϑ) ϑ 2 . Thus for

	With this φ can be written in terms of delay coordinates tan φ = y ′ x ′ = √ 3 x(t + 2 ϑ) -x(t) x(t ϑ = 1 tan φ(t) ≈ 2 √ 3 ẋ(t + 1) ẍ(t + 1) .

Acknowledgments

The author would like to thank C. Bandt for fruitful discussions. This work was supported by the priority programme SPP1114 of the German Research Foundation. The EEG data were provided by H. Lauffer from the Department of Pediatric Medicine of the University of Greifswald.