Detection of synchronization from univariate data using wavelet transform

Detecting regimes of synchronization between self-sustained oscillators is a typical problem in studying their interaction. Two types of interaction are generally recognized Blekhman I.I. (1971, 1988); Pikovsky A., Rosenblum M., Kurths J. (2001); Boccaletti S., Kurths J., Osipov G., Valladares D.L., Zhou C. (2002). The first one is a unidirectional coupling of oscillators. It can result in synchronization of a self-sustained oscillator by an external force. In this case the dynamics of the oscillator generating the driving signal does not depend on the driven system behavior. The second type is a mutual coupling of oscillators. In this case the interaction can be more effective in one of the directions, approaching in the limit to the first type, or can be equally effective in both directions. In the event of mutual coupling, synchronization is the result of the adjustment of rhythms of interacting systems. To detect synchronization one can analyze the ratio of instantaneous frequencies of interacting oscillators and the dynamics of the generalized phase difference Pikovsky A., Rosenblum M., Kurths J. (2001). As a quantitative characteristic of synchronization one can use the phase synchronization index Rosenblum M., Pikovsky A., Kurths J., Schafer C., Tass P. (2001); Meinecke F.C., Ziehe A., Kurths J., Müller K.-R. (2005) or the measure of synchronization Hramov A.E., Koronovskii A.A. (2004); Hramov A.E., Koronovskii A.A., Kurovskaya M.K., Moskalenko O.I. (2005).

Synchronization of interacting systems including the chaotic ones has been intensively studied in recent years. The main ideas in this area have been introduced using standard models Blekhman I.I. (1971, 1988); Pecora L.M., Carroll T.L. (1990); Pecora L.M., Carroll T.L., Jonson G.A., Mar D.J. (1997); Pikovsky A., Rosenblum M., Kurths J. (2000); Boccaletti S., Pecora L.M., Pelaez A. (2001); Pikovsky A., Rosenblum M., Kurths J. (2001); Boccaletti S., Kurths J., Osipov G., Valladares D.L., Zhou C. (2002); Rulkov N.F., Sushchik M.M., Tsimring L.S., Abarbanel H.D.I. (1995); Pyragas K. (1996); Hramov A.E., Koronovskii A.A. (2004); Hramov A.E., Koronovskii A.A., Kurovskaya M.K., Moskalenko O.I. (2005)). At present, more attention is focused on application of the developed techniques to living systems. In particular, much consideration is being given to investigation of synchronization between different brain areas Tass et al. (1998, 2003); Meinecke F.C., Ziehe A., Kurths J., Müller K.-R. (2005); Chavez M., Adam C., Navarro, Boccaletti S., Martinerie J. (2005) and to studying synchronization in the human cardiorespiratory system Schäfer C., Rosenblum M.G., Abel H.-H., Kurths J. (1999); Bračič-Lotrič M., Stefanovska A. (2000); Rzeczinski S., Janson N.B., Balanov A.G., McClintock P.V.E. (2002); Prokhorov M.D., Ponomarenko V.I., Gridnev V.I., Bodrov M.B., Bespyatov A.B. (2003); Hramov A.E., Koronovskii A.A., Ponomarenko V.I., Prokhorov M.D. (2006). Investigating such systems one usually deals with the analysis of short time series heavily corrupted by noise. In the presence of noise it is often difficult to detect the transitions between synchronous and nonsynchronous regimes. Besides, even in the region of synchronization a $2\pi$-phase jumps in the temporal behavior of the generalized phase difference can take place. Moreover, the interacting systems can have a set of natural rhythms. That is why it is desirable to analyze synchronization and phase locking at different time scales Hramov A.E., Koronovskii A.A. (2004); Hramov A.E., Koronovskii A.A., Levin Yu.I (2005); Hramov A.E., Koronovskii A.A. (2005); Chavez M., Adam C., Navarro, Boccaletti S., Martinerie J. (2005); Hramov A.E., Koronovskii A.A., Popov P.V., Rempen I.S. (2005).

A striking example of interaction between various rhythms is the operation of the human cardiovascular system (CVS). The main rhythmic processes governing the cardiovascular dynamics are the main heart rhythm, respiration, and the process of slow regulation of blood pressure and heart rate having in humans the fundamental frequency close to 0.1 Hz Malpas S. (2002). Owing to interaction, these rhythms appear in various signals: electrocardiogram (ECG), blood pressure, blood flow, and heart rate variability (HRV) Stefanovska A., Hožič M. (2000). Recently, it has been found that the main rhythmic processes operating within the CVS can be synchronized Schäfer C., Rosenblum M.G., Abel H.-H., Kurths J. (1999); Bračič-Lotrič M., Stefanovska A. (2000); Rzeczinski S., Janson N.B., Balanov A.G., McClintock P.V.E. (2002); Prokhorov M.D., Ponomarenko V.I., Gridnev V.I., Bodrov M.B., Bespyatov A.B. (2003). It has been shown that the systems generating the main heart rhythm and the rhythm associated with slow oscillations in blood pressure can be regarded as self-sustained oscillators, and that the respiration can be regarded as an external forcing of these systems Prokhorov M.D., Ponomarenko V.I., Gridnev V.I., Bodrov M.B., Bespyatov A.B. (2003); Rzeczinski S., Janson N.B., Balanov A.G., McClintock P.V.E. (2002).

Recently, we have proposed a method for detecting the presence of synchronization of a self-sustained oscillator by external driving with linearly varying frequency Hramov A.E., Koronovskii A.A., Ponomarenko V.I., Prokhorov M.D. (2006). This method was based on a continuous wavelet transform of both the signals of the self-sustained oscillator and external force. However, in many applications the diagnostics of synchronization from the analysis of univariate data is a more attractive problem than the detection of synchronization from multivariate data. For instance, the record of only a univariate signal may be available for the analysis or simultaneous registration of different variables may be rather difficult. In this paper we propose a method for detection of synchronization from univariate data. However, a necessary condition for application of our method is the presence of a driving signal with varying frequency. For the mentioned above cardiovascular system our method gives a possibility to detect synchronization between its main rhythmic processes from the analysis of the single heartbeat time series recorded under paced respiration.

The paper is organized as follows. In Sec. II we describe the method for detecting synchronization from univariate data. In Sec. III the method is tested by applying it to numerical data produced by a driven asymmetric van der Pol oscillator. In Sec. IV the method is used for detecting synchronization from experimental time series gained from a driven electronic oscillator with delayed feedback. Section V presents the results of the method application to studying synchronization between the rhythms of the cardiovascular system from the analysis of the human heart rate variability data. In Sec. VI we summarize our results.

Let us consider a self-sustained oscillator driven by external force ${\cal F}$ with varying frequency

where H is the operator of evolution, $\varepsilon$ is the driving amplitude, and $\Phi(t)$ is the phase of the external force defining the law of the driving frequency $\omega_{d}(t)$ variation:

In the simplest case the external force is described by a harmonic function ${\cal F}(\Phi(t))=\sin\Phi(t)$.

Assume that we have at the disposal a univariate time series $x(t)$ characterizing the response of the oscillator (1) to the driving force $\cal F$. Let us define from this time series the phase $\varphi_{0}(t)$ of oscillations at the system (1) basic frequency $f_{0}$. The main idea of our approach for detecting synchronization from univariate data is to consider the temporal behavior of the difference between the oscillator instantaneous phases at the time moments $t$ and $t+\tau$. We calculate the phase difference

where $\tau$ is the time shift that can be varied in a wide range. Note, that $\varphi_{0}(t)$ and $\varphi_{0}(t+\tau)$ are the phases of the driven self-sustained oscillator corresponding to oscillations at the first harmonic of the oscillator basic frequency $f_{0}$.

The variation of driving frequency is crucial for the proposed method. Varying in time, the frequency of the external force sequentially passes through the regions of synchronization of different orders $1:1$, $2:1$, …, $n:1$, …, $n:m$, …($n,m=1,2,3,\dots$). Within the time intervals corresponding to asynchronous dynamics the external signal practically has no influence on the dynamics of the basic frequency $f_{0}$ in the oscillator (1) spectrum. Thus, the phase of oscillator varies linearly outside the regions of synchronization, $\varphi_{0}(t)=2\pi f_{0}t+\bar{\varphi}$, where $\bar{\varphi}$ is the initial phase. Then, from Eq. (3) it follows

i.e., the phase difference $\Delta\varphi_{0}(t)$ is constant within the regions of asynchronous dynamics.

Another situation is observed in the vicinity of the time moments $t_{ns}$ where the driving frequency $\omega_{d}(t)\approx(2\pi n/m)f_{0}$ and $n:m$ synchronization takes place. For simplicity let us consider the case of $1:1$ synchronization. In the synchronization (Arnold) tongue the frequency of the system (1) nonautonomous oscillations is equal to the frequency (2) of the external force and the phase difference between the phase of the driven oscillator $\varphi_{0}(t)$ and the phase $\Phi(t)$ of the external force, $\Delta\tilde{\phi}(t)=\varphi_{0}(t)-\Phi(t)$, is governed in a first approximation by the Adler equation Adler R. (1947). It follows from the Adler equation that in the region of $1:1$ synchronization the phase difference $\Delta\tilde{\phi}(t)$ varies by $\pi$.

Representing the driven oscillator phase as $\varphi_{0}(t)=\Delta\tilde{\phi}(t)+\Phi(t)$, we obtain from Eq. (3):

where $\gamma=\Delta\tilde{\phi}(t+\tau)-\Delta\tilde{\phi}(t)\approx\rm const$ is the correction of the phase difference that appears due to synchronization of the system by external force. Expanding the phase $\Phi(t+\tau)$ in a Taylor series we obtain

Taking into account Eq. (2) we can rewrite Eq. (6) as

Thus, the behavior of the phase difference (3) is defined by the law of the driving frequency $\omega_{d}(t)$ variation.

For the linear variation of the driving frequency, $\omega_{d}(t)=\alpha+\beta t$, from Eq. (7) it follows

Consequently, in the region of synchronization the phase difference varies linearly in time, $\Delta\varphi_{0}(t)\sim t$. In the case of the nonlinear variation of $\omega_{d}(t)$, the dynamics of $\Delta\varphi_{0}(t)$ is more complicated. However, if $\omega_{d}(t)$ varies in a monotone way and the time of its passing through the synchronization tongue is small, one can neglect the high-order terms of the expansion and consider the law of $\Delta\varphi_{0}(t)$ variation as the linear one. We will show below that this assumption holds true for many applications.

The absolute value of the change in the phase difference $\Delta\varphi_{0}(t)$ within the synchronization region can be estimated using Eq. (7):

where $t_{1}$ and $t_{2}$ are the time moments when the frequency of the external force passes through, respectively, the low-frequency and high-frequency boundaries of the synchronization tongue. Assuming that the rate of $\omega_{d}(t)$ variation is slow, we can neglect the terms containing the derivatives of $\omega_{d}(t)$ and obtain

where $\Delta\omega=\omega_{d}(t_{2})-\omega_{d}(t_{1})$ is the bandwidth of synchronization.

The obtained estimation corresponds to the case of $1:1$ synchronization, characterized by equal values of the driving frequency $f_{d}$ and the oscillator frequency $f_{0}$, $f_{d}/f_{0}=1$. However, the considered approach can be easily extended to a more complicated case of $n:m$ synchronization. In this case the change in $\Delta\varphi_{0}(t)$ within the region of synchronization takes the value

Hence, the analysis of the phase difference (3) behavior allows one to distinguish between the regimes of synchronous and asynchronous dynamics of driven oscillator. The phase difference $\Delta\varphi_{0}(t)$ is constant for the regions of asynchronous dynamics and demonstrates monotone (often almost linear) variation by the value $\Delta\varphi_{s}$ defined by Eq. (11) within the regions of synchronization.

To define the phase $\varphi_{0}(t)$ of oscillations at the basic frequency we use the approach based on the continuous wavelet transform Koronovskii A.A., Hramov A.E. (2004); Hramov A.E., Koronovskii A.A. (2004, 2005); Hramov A.E., Koronovskii A.A., Kurovskaya M.K., Moskalenko O.I. (2005). It is significant, that the wavelet transform Wav (2004); Koronovskii A.A., Hramov A.E. (2003) is the powerful tool for the analysis of nonlinear dynamical system behavior. The continuous wavelet analysis has been applied in the studies of phase synchronization of chaotic neural oscillations in the brain Lachaux:1999 ; Lachaux:2000 ; Lachaux:2001 ; Lachaux:2002_BrainCoherence ; Quyen:2001_HTvsWVT , electroencephalogram signals Quiroga:2002 , R–R intervals and arterial blood pressure oscillations in brain injury Turalska:2005 , chaotic laser array DeShazer:2001_WVT_LaserArray . It has also been used to detect the main frequency of the oscillations in nephron autoregulation Sosnovtseva:2002_Wvt and coherence between blood flow and skin temperature oscillations BANDRIVSKYY:2004 . In these recent studies a continuous wavelet transform with various mother wavelet functions has been used for introducing the instantaneous phases of analyzed signals. In particular, in Refs. Lachaux:2001 ; Quiroga:2002 a comparison of Hilbert transform and wavelet method with the mother Morlet wavelet has been carried out and good conformity between these two methods has been shown for the analysis of neuronal activity. It is important to note, that in all the above mentioned studies the wavelet transform has been used for the analysis of synchronization from bivariate data, when the generalized phase difference $\Delta\varphi(t)$ of both analyzed rhythms was investigated. The proposed method allows one to detect synchronization from the analysis of only the one signal of the oscillator response to the external force with monotonically varying frequency. Taking into account the high efficiency of the analysis of synchronization with the help of the continuous wavelet transform using bivariate data, we will use the continuous wavelet transform for determining the instantaneous phase of the analyzed univariate signal.

The continuous wavelet transform Wav (2004); Koronovskii A.A., Hramov A.E. (2003) of the signal $x(t)$ is defined as

where $\psi_{s,t_{0}}(t)$ is the wavelet function related to the mother wavelet $\psi_{0}(t)$ as $\psi_{s,t_{0}}(t)=\left({1}/{\sqrt{s}}\right)\psi_{0}\left(({t-t_{0}})/{s}\right)$. The time scale $s$ corresponds to the width of the wavelet function, $t_{0}$ is the shift of the wavelet along the time axis, and the asterisk denotes complex conjugation. It should be noted that the wavelet analysis operates usually with the time scale $s$ instead of the frequency $f$, or the corresponding period $T=1/f$, traditional for the Fourier transform.

The wavelet spectrum

describes the system dynamics for every time scale $s$ at any time moment $t_{0}$. The value of $|W(s,t_{0})|$ determines the presence and intensity of the time scale $s$ at the time moment $t_{0}$. We use the complex Morlet wavelet Grossman A. and Morlet J. (1984) $\psi_{0}(\eta)=({1}/{\sqrt[4]{\pi}})\exp[j\sigma\eta]\exp\left[{-\eta^{2}}/{2}\right]$ as the mother wavelet function. The choice of the wavelet parameter $\sigma=2\pi$ provides the simple relation $f\approx 1/s$ between the frequency $f$ of the Fourier transform and the time scale $s$ Koronovskii A.A., Hramov A.E. (2003).

In this section we investigate synchronization between the respiration and rhythmic process of slow regulation of blood pressure and heart rate from the analysis of univariate data in the form of the heartbeat time series. This kind of synchronization has been experimentally studied in Prokhorov M.D., Ponomarenko V.I., Gridnev V.I., Bodrov M.B., Bespyatov A.B. (2003); Hramov A.E., Koronovskii A.A., Ponomarenko V.I., Prokhorov M.D. (2006); Janson:2001_PRL ; Janson:2002_PRE . We studied eight healthy volunteers. The signal of ECG was recorded with the sampling frequency 250 Hz and 16-bit resolution. Note, that according to Circulation:1996 the sampling frequency 250 Hz used in our experiments suffices to detect accurately the time moment of R peak appearance. The experiments were carried out under paced respiration with the breathing frequency linearly increasing from 0.05 Hz to 0.3 Hz within 30 min. We specially included the lower frequencies for paced respiration in order to illustrate the presence of the most pronounced regime of 1:1 synchronization between the respiration and slow oscillations in blood pressure. The rate of respiration was set by sound pulses. The detailed description of the experiment is given in Ref. Prokhorov M.D., Ponomarenko V.I., Gridnev V.I., Bodrov M.B., Bespyatov A.B. (2003).

Extracting from the ECG signal the sequence of R–R intervals, i.e., the series of the time intervals between the two successive R peaks, we obtain the information about the heart rate variability. The proposed method of detecting synchronization from uniform data was applied to the sequences of R–R intervals.

A typical time series of R–R intervals for breathing at linearly increasing frequency is shown in Fig. 8a. Since the sequence of R–R intervals is not equidistant, we exploit the technique for applying the continuous wavelet transform to nonequidistant data. The wavelet spectra $|W(s,t_{0})|$ for different parameters $\sigma$ of the Morlet wavelet are shown in Figs. 8b and 8c for the sequence of R–R intervals presented in Fig. 8a. For greater $\sigma$ values the wavelet transform provides higher resolution of frequency Koronovskii A.A., Hramov A.E. (2003) and better identification of the dynamics at the time scales corresponding to the basic frequency of oscillations and the varying respiratory frequency. In the case of $\sigma=2\pi$ the time scale $s$ of the wavelet transform is very close to the period $T$ of the Fourier transform and the values of $s$ are given in seconds in Fig. 8b. Generally, the time scale $s$ is related to the frequency $f$ of the Fourier transform by the following equation:

Because of this, the units on the ordinates are different in Figs. 8b and 8c. The wavelet spectra in these figures demonstrate the high-amplitude component corresponding to the varying respiratory frequency manifesting itself in the HRV data. The self-sustained slow oscillations in blood pressure (Mayer wave) have in humans the basic frequency of about 0.1 Hz, or respectively, the basic period close to 10 s. The power of this rhythm in the HRV data is less than the power of respiratory oscillations. As the result, the time scale $s_{0}$ is weakly pronounced in the spectra.

Fig. 9 presents the phase differences $\Delta\varphi_{0}(t)$ calculated for R–R intervals of four subjects under respiration with linearly increasing frequency. All the curves in the figure exhibit the regions with almost linear in the average variation of $\Delta\varphi_{0}(t)$ indicating the presence of synchronous dynamics. In particular, the region of $1:1$ synchronization is observed within the interval 200–600 s when the frequency of respiration is close to the basic frequency of the Mayer wave. This region is marked by arrow. In this region the frequency of blood pressure slow oscillations is locked by the increasing frequency of respiration and increases from 0.07 Hz to 0.14 Hz. Outside the interval of synchronization, $t<200$ s and $t>600$ s, the phase differences demonstrate fluctuations caused by the high level of noise and nonstationarity of the experimental data. Some of these fluctuations take place around an average value as well as in the case of the driven van der Pol oscillator affected by noise (see Fig. 3). The frequency of blood pressure slow oscillations demonstrates small fluctuations around the mean value of about 0.1 Hz outside the interval of synchronization.

The phase differences in Fig. 9a are plotted for different $\tau$. As the time shift $\tau$ increases, so does the range of $\Delta\varphi_{0}(t)$ monotone variation in the region of synchronization. This result agrees well with the results presented in Sec. III. Similar behavior of $\Delta\varphi_{0}(t)$ is observed for each of the eight subjects studied. In Fig. 9(b) phase differences $\Delta\varphi_{0}(t)$ computed for R-R intervals of another three subjects are presented. The phase differences demonstrate the wide regions of almost linear variation for all the subjects. Such behavior of the considered phase difference cannot be observed in the absence of synchronization, if only the modulation of blood pressure oscillations by respiration is present. These results allow us to confirm the conclusion that the slow oscillations in blood pressure can be synchronized by respiration. However, to come to this conclusion, the proposed method needs only univariate data in distinction to the methods Prokhorov M.D., Ponomarenko V.I., Gridnev V.I., Bodrov M.B., Bespyatov A.B. (2003); Hramov A.E., Koronovskii A.A., Ponomarenko V.I., Prokhorov M.D. (2006) based on the analysis of bivariate data. Note, that paper Prokhorov M.D., Ponomarenko V.I., Gridnev V.I., Bodrov M.B., Bespyatov A.B. (2003) contains the more detailed investigation of synchronization between the respiration and slow oscillations in blood pressure than the present one. Recent reports (see, for examples, Rosenblum:1998_Nature ; Suder:1998_AJP ; Kotani:2000_MIM ) focused on examining the relationship between respiration and heart rate have shown that there is nonlinear coupling between respiration and heart rate. In particular, such coupling is well studied for the respiratory modulation of heart rate Bishop:1981_AJP ; Kotani:2000_MIM known as respiratory sinus arrhythmia. The presence of coupling between the cardiac and respiratory oscillatory processes has been revealed also using bispectral analysis in Jamsek:2003_PRE ; Jamsek:2004_PMB under both spontaneous and paced respiration. Our results are in agreement with those when synchronization between the oscillating processes occurs as the result of their interaction.

We have proposed the method for detecting synchronization from univariate data. The method allows one to detect the presence of synchronization of the self-sustained oscillator by external force with varying frequency. To implement the method one needs to analyze the difference between the oscillator instantaneous phases calculated at time moments shifted by a certain constant value with respect to each other. The instantaneous phases are defined at the oscillator basic frequency using continuous wavelet transform with the Morlet wavelet as the mother wavelet function. The necessary condition for the method application is the variation of the frequency of the driving signal. The method efficiency is illustrated using both numerical and experimental univariate data under sufficiently high levels of noise and inaccuracy of the basic time scale definition.

We applied the proposed method to studying synchronization between the respiration and slow oscillations in blood pressure from univariate data in the form of R–R intervals. The presence of synchronization between these rhythmic processes is demonstrated within the wide time interval. The knowledge about synchronization between the rhythms of the cardiovascular system under paced respiration is useful for the diagnostics of its state N. Ancona, R. Maestri, D. Marinazzo, L. Nitti, M. Pellicoro, G.D. Pinna, S. Stramaglia (2005). The method allows one to detect the presence of synchronization from the analysis of the data of Holter monitor widely used in cardiology.

The proposed method can be used for the analysis of synchronization even in the case when the law of the driving frequency variation is unknown. If the frequency of the external driving varies in the wide range, the analysis of the oscillator response to the unknown driving force allows one to make a conclusion about the presence or absence of synchronization in the system under investigation.

We thank Dr. Svetlana Eremina for English language support. This work is supported by the Russian Foundation for Basic Research, Grants 05–02–16273, 07–02–00044, 07–02–00747 and 07–02–00589, and the President Program for support of the leading scientific schools in the Russian Federation, Grant No. SCH-4167.2006.2. A.E.H. acknowledges support from CRDF, Grant No. Y2–P–06–06. A.E.H. and A.A.K. thank the “Dynasty” Foundation for the financial support.