## Abstract

In this paper, we present an approach to recover the dynamics from recurrences of a system and then generate (multivariate) twin surrogate (TS) trajectories. In contrast to other approaches, such as the linear-like surrogates, this technique produces surrogates which correspond to an independent copy of the underlying system, i.e. they induce a trajectory of the underlying system visiting the attractor in a different way. We show that these surrogates are well suited to test for complex synchronization, which makes it possible to systematically assess the reliability of synchronization analyses. We then apply the TS to study binocular fixational movements and find strong indications that the fixational movements of the left and right eye are phase synchronized. This result indicates that there might be only one centre in the brain that produces the fixational movements in both eyes or a close link between the two centres.

## 1. Introduction

The concepts of complex synchronization and especially phase synchronization (PS) have been intensively studied in recent years (Pikovsky *et al*. 2001). Indications of PS have been found in many laboratory and natural systems, e.g. population dynamics (Blasius *et al*. 1999), neurological (Tass *et al*. 1998), cardiorespiratory systems (Schäfer *et al*. 1998) and optics (Allaria *et al*. 2001), or during epilepsy (Mormann *et al*. 2003; Sowa *et al*. 2005). The corresponding studies are usually based on the computation of a measure which quantifies dependencies of the instantaneous phases of the time-series. However, even though these measures may be normalized between 0 and 1, experimental time-series often yield values which are neither close to 0 nor to 1 and hence are difficult to interpret. This problem can be overcome if the coupling strength between the two systems can be varied systematically and a rather large change in the measure can be observed (‘active experiment’; Pikovsky *et al*. 2001). PS in natural systems, e.g. during epilepsy or between the heart beats of a mother with the ones of her foetus (Van Leeuwen *et al*. 2003), frequently evades such an experimental manipulation (‘passive experiment’). In some cases, this problem has been tackled by interchanging the pairs of oscillators (Van Leeuwen *et al*. 2003), e.g. the heart beats of other pregnant women were used as ‘natural surrogates’. These surrogates are independent and hence not in PS with the original system. Hence, if the synchronization index obtained for the original data is not significantly higher than the index obtained for the natural surrogates, there is no sufficient evidence to claim synchronization. But, even this approach has some drawbacks. The natural variability and also the frequency of the heart beats of the surrogate mothers are usually slightly different from the ones of the biological mother. Furthermore, the data acquisition can be expensive and at least in some cases, problematic (e.g. in some states of the pregnancy). Moreover, frequently, e.g. in geophysics or neuroscience, such natural surrogates are often not available. In these cases, we propose to test hypotheses on the basis of surrogates generated by a mathematical algorithm. Several approaches in this direction have been published (Palus 1997; Palus & Stefanovska 2003). However, the specificity of these tests is not always satisfactory (Thiel *et al*. 2006).

In this paper, we show how to reconstruct an attractor from the system's recurrences and as a consequence, suggest that they contain the necessary information about the dynamics of a system to produce ‘alternative’ trajectories. We then present an algorithm for the generation of surrogates, which is based on the recurrences of a system. These surrogates mimic the dynamical behaviour of the system, i.e. not only the linear properties but also the nonlinear ones are preserved. They correspond to an independent copy of the total system, and hence they are not in PS with the original one. The main idea is to find points of the measured trajectory which are not only neighbours, but also share the same neighbourhood in phase space (twins), i.e. all other points are either neighbours of both or of neither of them. Once the twins of the trajectory have been localized, new surrogate trajectories are generated by substituting randomly the next step in the trajectory by either its own future or the one of its twin.

## 2. Recurrence plots

The algorithm to generate the surrogates is based on the recurrence matrix(2.1)where denotes the Heaviside function, the maximum norm and *δ* is a predefined threshold. The (if necessary reconstructed) vectors and describe the state of the system at times *i* and *j*. Coding the ‘1's’ in the matrix as black dots and the ‘0's’ as white ones, we obtain the recurrence plot (RP) of the trajectory (Eckmann *et al*. 1987). The seemingly easy-to-interpret representation of the matrix, has led to the definition of many ad hoc measures to quantify structures in RPs. The resulting recurrence quantification analysis (RQA) has made it possible to successfully study the measured time-series by means of RPs (Weber & Zbilut 1994; Marwan *et al*. 2002*a*,*b*).1

It has also been reported that RPs can be used to estimate dynamical invariants from time-series in a rather robust way and independent of the embedding parameters (Thiel *et al*. 2004*a*).

## 3. Information content of recurrence plots

The recurrence matrix as given by equation (2.1) can be considered as a representation of a specific aspect of a dynamical system. It encodes information about when the system recurs to formerly visited states in its phase space. An important question is how much information about a system's dynamics is contained in the recurrence matrix. The trajectory from which it is computed ‘lives’ in an *m*-dimensional space . The matrix *R*_{i,j} on the other hand is two-dimensional and binary. This suggests that much of the information about the dynamics might be lost in the process of computing the RP.

However, many dynamical invariants such as entropies, fractal dimensions and the mutual information can be estimated from RPs (Thiel *et al*. 2004*a*). A rather unexpected result is that it is possible to reconstruct the latter and, then by Taken's theorem, the attractor from it, if the RP is computed from a univariate time-series (Thiel *et al*. 2004*b*). The possibility to reconstruct the attractor means, roughly speaking, that one can recover the dynamics of the system from the RP. Note that the algorithm presented (Thiel *et al*. 2004*b*), cannot be used to reconstruct the attractor if the RP was computed from the vectors in the phase space, as it involves a rank ordering of the entries of the time-series—a procedure which is only possible for scalar values.

This result is quite surprising as the RP can be supposed to be more appropriate if it represents recurrences in the phase space, rather than recurrences of a projection of the attractor to one dimension. This is why one would expect that embedding a scalar time-series is a necessary procedure before computing the RP. Indeed, one can easily see that some information about the system is lost if one embeds—at least in the RP. Let us suppose that we compute the RP for the *x*-component of the Rössler system. Then we obtain(3.1)where denotes the absolute value. The superscript ‘u’ indicates that this matrix is computed from the univariate time-series. If we now want to compute the RP for the embedded time-series (let us assume that the embedding parameters are given by the dimension *m* and the delay *τ*), we can first consider the embedded vectors and then use equation (2.1) to compute the recurrence matrix *R*_{i.j} from these vectors. Alternatively, we can use the formula(3.2)This equation has a very intuitive interpretation. To compute the *R*_{i,j} of an attractor based on , one proceeds as follows. First, one takes *m*−1 copies of the RP of a univariate time-series. One then shifts the *m*th copy steps along the main diagonal and multiplies the entries of the original and the entire shifted matrix pointwise. The result is an RP which equals the RP obtained if one embeds first and then uses equation (2.1). This means that one can either first embed and then compute the RP, or compute a ‘univariate’ RP and then apply our shifting scheme to obtain the same result.

Note that it is hence possible to compute the *R*_{i,j} from the univariate , but not vice versa. A ‘0’ in *R*_{i,j} does not tell us which of the *m* components of the vector was/were further away than *δ*. Hence, we cannot use the algorithm presented in Thiel *et al*. (2004*b*) to reconstruct the dynamics from *R*_{i,j}.

Instead, we attempt to reconstruct the attractor based on an algorithm which has been proposed for the reconstruction of protein structures from distance inequalities by Bohr *et al*. (1993) from the RP alone. This means, based only on the recurrence matrix, we reconstruct the attractor. This algorithm has allowed Bohr *et al*. to reconstruct the structure of proteins containing several hundreds of amino acids from their distance inequalities, i.e. the information if two amino acids are close in space or not. If we reinterpret the amino acids as ‘points in phase space’ and the protein as the ‘attractor’, it becomes clear that a slight modification of the algorithm might allow us to reconstruct the attractor from an RP.

We give the same argument here as Bohr *et al*. (1993) but adapt the nomenclature to our problem. From now on, we use the Euclidean distance in equation (2.1). Let us first suppose that instead of the recurrence matrix (2.1), we have the distance matrix(3.3)where denotes the Euclidean distance. We now write out the components of the position vector of the *i*th point in the phase space assuming first a three-dimensional space. One can then find suitable coordinates for each point by minimization of the cost function(3.4)If all coordinates are suitably chosen, *H* vanishes. If only distance inequalities are known, i.e. we have got information only about whether two points are closer than a given *δ* or not, we have got to modify the cost function. Bohr *et al*. have proposed the following modification. They have introduced(3.5)where *R*_{i,j} is the recurrence matrix and the function is given by(3.6)Here *s* is a free parameter, which has to be suitably tuned. *δ*′ would optimally be the threshold *δ* from the RP, which usually is unknown when one attempts to reconstruct the attractor from the RP.2 However, choosing , one only stretches the reconstructed structure by the factor . Note that is a smoothed step function, which allows us to minimize the cost function.

The procedure of the minimization is computationally rather expensive. Suppose, for the sake of the argument that we want to reconstruct an RP of a structure in the phase space that consists of 10 000 points, then we have in a three-dimensional space 30 000 parameters with respect to which we have to minimize the cost function. We can expect to find a local minimum rather than a global minimum. This leads to a faulty reconstruction. One could, in a next step, compute the RP from the reconstructed (multivariate) time-series and would find the points that are not correctly reconstructed. An iterative procedure might help to improve the results.

Figure 1*a* shows a sine function embedded in . The reconstructed trajectory (figure 1*b*) is a closed curve and at least topologically equivalent to a circle in the phase space. Figure 1*c*,*d* show the same for a plane in and figure 1*e*,*f* for the chaotic Lorenz attractor (for equations see Lorenz (1963)). In all the three cases, the corresponding correlation entropies and dimensions are the same. Furthermore, after rank ordering the respective *x*-, *y*- and *z*-components of the original and reconstructed attractors, the rank-ordered components can be mapped onto each other by a continuous function (after a suitable rotation of the reconstructed trajectory in the phase space). This indicates a topological equivalence between the original and the reconstructed trajectory.

The algorithm performs reasonably well, given the large number of variables with respect to which the cost function has to be minimized. This shows that the distance inequalities, i.e. the RP, contain enough information to reconstruct the attractor. If one computes the RP from the reconstructed set, one finds rather strong deviations from the original RP. This indicates that we find only a local minimum of the cost function. The iterative procedure described above might yield better results. Note that the only crucial parameter for the reconstruction is the threshold *δ* used for the computation of the recurrence matrix (equation (2.1)). Our simulations show that there exists a broad interval, from which *δ* can be chosen. As expected, the performance of the algorithm improves with increasing length of the time-series.

*These simulations raise the question as to whether distance inequalities or recurrences are, in general, enough to reconstruct an attractor. Is the RP in some sense equivalent to the ‘knowledge of the dynamics’?* Alternatively, we can ask what two structures have in common if they have the same recurrence matrix. Obviously, two structures or set of points that are equivalent up to translation, rotation and reflection have the same recurrence matrix. Isometric transformations do not change the recurrence matrix. However, even non-isometric transformations can leave the recurrence matrix unchanged.3 A hypothesis would be that two point sets (i.e. points on an attractor) with the same recurrence matrix have to be at least topologically equivalent; i.e. there is a homeomorphism between the two. This would mean that the dynamics is, in some sense, the same in both systems. However, in one dimension this can be shown to be wrong.

Even though this problem is still not solved, the results obtained so far suggest that the RP contains a very large amount of (if not all the relevant) information of the dynamical system. In §4 we will make use of this fact to construct surrogates to test for synchronization.

## 4. Generation of surrogates

A first idea for the generation of surrogates is to change the structures in an RP consistently with the ones produced by the underlying dynamical system. One could then reconstruct a new realization of the trajectory from the modified RP. However, the structures in RPs are not fully understood and one cannot arbitrarily interchange columns in an RP because such a modification changes the distribution of diagonal lines and hence the entropy and predictability of the system (Thiel *et al*. 2004*a*).

Therefore, we propose a modified approach. In general, in an RP, there are identical columns, i.e. (Thiel *et al*. 2004*b*). Thus, there are points which are not only neighbours (i.e. ) but also share the same neighbourhood. Reconstructing the attractor from an RP, the respective neighbourhoods of these points cannot help to distinguish them, i.e. from this point of view they are identical. This is why we will call them *twins*.

Twins are special points of the time-series as they are indistinguishable considering their neighbourhoods but in general different and, hence, have different pasts and—more important—different futures. The key idea of how to introduce the randomness needed for the generation of surrogates of a deterministic system is that one can jump randomly to one of the possible futures of the existing twins.

A surrogate trajectory of with is then generated in the following way: (i) identify all pairs of twins, (ii) choose an arbitrary starting point, say , (iii) if has no twin, the next point of the surrogate trajectory is , and (iv) if has a twin, say , then one can go to either or to , i.e. or with equal probability.4 Steps (iii) and (iv) are then iterated until the surrogate time-series has the same length as the original one.

This algorithm creates twin surrogates (TS) which shadow a (typical) trajectory of the system (Ott 1993; Katok & Hasselblatt 1995). In the limit of an infinitely long trajectory, its surrogates are characterized by the same dynamical invariants and the same attractor. If the measure of the attractor can be estimated from the observed finite trajectory reasonably well, its surrogates share the same statistics. Also their power spectra and correlation functions are consistent with the ones of the original system (Thiel *et al*. 2006). TS give reasonable results for deterministic and also for stochastic systems; the TS of, for example, an ARMA process also show the typical behaviour of a linear Gaussian process if a suitable ‘embedding’ is used. Even the parameters of the process can be estimated correctly from the surrogates. A more detailed analysis of the TS can be found in Thiel *et al*. (2006).

## 5. Test for synchronization

Next, we use the TS to test for PS. The idea behind this approach is similar to the one by means of ‘natural surrogates’ in the mother–foetus heartbeat synchronization (Van Leeuwen *et al*. 2003). Suppose that we have two coupled self-sustained oscillators and . Then, we generate *M* TS of the total system, i.e. and , with . These surrogates are independent copies of the total system, i.e. trajectories of the whole system beginning at different initial conditions. Note that the coupling between and is also mimicked by the surrogates. Next, we compute the differences between the phases of the original system applying for example the analytical signal approach (Pikovsky *et al*. 2001) and compare them with (one can also consider ). Then, if does not differ significantly from with respect to some index for PS, the null hypothesis cannot be rejected, and hence we do not have enough evidence to state PS.

As a test case, we consider two non-identical, mutually coupled Rössler oscillators with a frequency mismatch of *ν*=0.015.5 In this ‘active experiment’, we vary the coupling strength *ϵ* from 0 to 0.08 and compute a PS index for the original trajectory for each value of *ϵ*. Next, we generate 200 TS and compute the PS index between the measured first oscillator and the surrogates of the second one. As PS index, we use the mean resultant length RL of complex phase vectors(5.1)It takes on values in the interval from 0 (non PS) to 1 (perfect PS; Rodriguez *et al*. 1999; Allefeld & Kurths 2004). If , where represents the PS index between the first oscillator and the surrogate *i* of the second one, we state that at a specified *α*-level the null hypothesis cannot be rejected, and hence we do not have enough evidence to claim PS between both oscillators.

Figure 2*a* shows the results for RL of the original system (bold line) and the 1%(solid) significance level. Figure 2*b* displays the difference between RL of the original system and the 1, 2 and 5% significance level. For *ϵ*<0.025, RL of the original system is, as expected, below the significance levels and hence the difference is negative, and for higher values of *ϵ* the curves cross (the difference becomes positive). This is in agreement with the criterion for PS via Lyapunov exponents *λ*_{i} (Pikovsky *et al*. 2001), i.e. *λ*_{4} becomes negative at *ϵ*∼0.028 (figure 2*b*), which approximately coincides with the intersection of the curve of RL for the original system and the significance level (zero-crossing of the curves in figure 2*b*). Therefore, we recognize successfully the PS region by means of the TS.

Also note that the significance limit increases when the transition to PS occurs (figure 2*a*). As the TS mimic both the linear and nonlinear characteristics of the system, the surrogates of the second oscillator have in the PS region the same mean frequency as the first original oscillator. Hence is rather high. However, and do not adapt to each other, as they are independent. Hence, the value of RL for the original system is significantly higher than the . We state in conclusion that even though the obtained value for a normalized PS index is higher than 0.97 (right side of figure 2*a*), this does not offer conclusive evidence for PS. Hence, *the knowledge of the PS index alone does not provide sufficient evidence for PS*.6

Next, we perform an analysis of the power of the test for *ϵ*=0 and *ν*=0. For 100 random initial conditions of the Rössler system and a significance level of *α*=1%, the null hypothesis was erroneously rejected only in 1 out of the 100 cases. This is a rather auspicious result as, due to the identical frequencies, it is extremely difficult to detect PS in this case. In the case of *ϵ*=0.02 (e.g. no PS) and a frequency mismatch *ν*=0.015, there were no erroneous rejections of the null hypothesis. Finally, for PS (*ϵ*=0.045 and *ν*=0.015), in all 100 test runs, the null hypothesis was correctly rejected. These results indicate that the power of the test is rather good. A more detailed analysis of this test can be found in Thiel *et al*. (2006).

## 6. Application to eye movements

Next we apply our algorithm to check fixational movements of the two eyes for PS. During fixation of a stationary target, our eyes perform small involuntary and allegedly erratic movements to counteract retinal adaptation. If these eye movements are experimentally suppressed, retinal adaptation to the constant input induces very rapid perceptual fading (Riggs *et al*. 1953; Coppola & Purves 1996). Moreover, statistical correlations show a time-scale separation from persistence to antipersistence (Engbert & Kliegl 2004). Persistence on the short time-scale counteracts retinal fading, whereas antipersistence on the long time-scale contributes to stability of ocular disparity. According to current textbook knowledge, the fixational movements of the left and right eye are correlated very poorly at best (Ciuffreda & Tannen 1995). Therefore, it is highly desirable to examine these processes from a perspective of PS. We analyse the data of two subjects. Each performed three trials, in which they fixated a small stimulus (black square on a white background, 3× pixels on a computer display) with a spatial extent of 0.12°, or 7.2 arc min. Eye movements were recorded using an EyeLink-II system (SR Research, Osgoode, Ontario, Canada) with a sampling rate of 500 Hz and an instrument spatial resolution less than 0.005°. Figure 3*a*,*b* shows a segment of the horizontal and vertical component of the eye movements, respectively, for one person.

The data were first high-pass filtered applying a difference filter with *τ*=40 ms in order to eliminate the slow drift of the data. After this filtering, we find an oscillatory trajectory (figure 4), which has maximum spectral power in the frequency range between 6 and 8 Hz. However, the trajectories of the eyes are rather noisy and non-phase coherent. Therefore, it is cumbersome to estimate the phase of these data. Hence, we apply another measure of PS which is based on the probability of recurrence of a trajectory in the phase space , where is the recurrence matrix (equation (2.1)). The correlation between the probabilities of recurrence of two interacting oscillators(6.1)(where means that the mean value has been subtracted and *σ*_{1} and *σ*_{2} are the standard deviations of *P*_{1}(*τ*), respectively, *P*_{2}(*τ*)), has been proposed to detect PS in non-phase coherent and noisy oscillators, where the phase cannot be estimated directly (Romano *et al*. 2005). Now, we compute 200 surrogates (figure 4*b*) of the left eye's trajectory and compute the recurrence-based synchronization index between them and the measured right eye's trajectory. In figure 4*c*, the results of the test of one trial are visualized. Table 1 summarizes the results for both subjects and all trials. In all cases, the PS index of the original data is significantly different from the ones of the surrogates, which strongly indicates that the concept of PS can be successfully applied to study the interaction between the trajectories of the left and right eye during fixation. This result also suggests that the physiological mechanism in the brain that produces the fixational eye movements controls both eyes simultaneously; i.e. there might be only one centre in the brain that produces the fixational movements in both eyes or a close link between two centres. Our finding of PS between the left and right eye is in good agreement with current knowledge of the physiology of the oculomotor circuitry. In a single-cell study, 66% of abducens motor neurons fired in relation to the movements of either eye, while premotor neurons in the brainstem encode monocular movements (Zhou & King 1998). Thus, motor neurons—as the final common pathway of neural control of eye movements—are candidates for the synchronization of left and right fixational movements. Furthermore, we are interested in whether the fixational movements in the horizontal and vertical direction of one eye are synchronized. Horizontal and vertical saccadic eye movements are controlled in two spatially distinct brainstem nuclei (Sparks 2002). Therefore, we can expect that, on the level of fixational eye movements, horizontal and vertical components are independent. Applying the method of TS to this case, we find that the synchronization index CPR between the *x*- and *y*-component is not significantly different from the one computed for the surrogates. Hence, we do not have sufficient evidence to claim synchronization between the horizontal and vertical components of eye movements.

## 7. Conclusions

In conclusion, we propose a new method for generating surrogates based on the concept of recurrence after presenting results that suggest that recurrences contain a vast amount of information about a dynamical system. These TS correspond to an independent copy of the underlying system; i.e. to a trajectory of the system starting at different initial conditions. We have shown that the TS can be used to test for PS in the well-studied system of two mutually coupled Rössler oscillators; and by means of the TS, we have detected PS in dependence on the coupling strength at several significance levels. Furthermore, we have tested for PS in experiments of binocular fixational movements and found that the left and right eye are in PS, in agreement with physiological results about the functional role of motor neurons in the final common pathway for the control of eye movements. Contrary to popular belief, fixational eye movements are a necessary condition for vision. Thus, an understanding of their dynamics is fundamental for perception and the associated control of spatial attention (Engbert & Kliegl 2003; Laubrock *et al*. 2005; Rolfs *et al*. 2005). First results indicate that the concept of TS can also be used for testing other complex kinds of synchronization, especially generalized synchronization.

## Acknowledgments

We want to thank Bohr *et al*. for providing the code of their reconstruction algorithm. This work has been supported by the DFG priority programme 1114, DFG-KL 955, and the ‘Internationales Promotionskolleg–Helmholtz Center for the Study of Mind and Brain Dynamics’ at the University of Potsdam.

## Footnotes

↵† Present address: Department of Physics, University of Aberdeen, Aberdeen AB24 3UE, UK.

One contribution of 14 to a Theme Issue ‘Experimental chaos II’.

↵A rather comprehensive bibliography about recurrence plots can be found at http://www.agnld.uni-potsdam.de/∼marwan/literature.php or alternatively at http://www.recurrence-plot.tk/bibliography.php.

↵Note that Bohr

*et al*. use futher normalization factors in their algorithm.↵Take a straight line, stretch it suitably and bend it into a half of a circle.

↵If triplets occur one proceeds analogously.

↵The equations are , .

↵Note that the more phase coherent the oscillators are, the more difficult it is to decide whether they are in PS or not. A certain phase diffusion, which allows one to measure the adaptation of the phases of the interacting oscillators is necessary to detect PS. However, the test based on the TS reveals whether there is enough evidence for PS.

- © 2007 The Royal Society