**Introduction **

Complexity of cardiovascular system (CVS) and its regulatory mechanism draws high researcher’s interest to this system. However, many problems are yet to be solved. It is well known that many interacting oscillation processes of different origin are present in CVS [1]. Many approaches have been utilized to study regulatory mechanisms of CVS, such as analysis of heart rate variability (HRV) [2], blood pressure variability (BPV) [3], peripheral blood flow variability [4-7], analysis of interactions between different oscillations in CVS [7-13], etc. Some authors proposed different mathematical models to describe the loops of autonomous regulation of blood circulation [14-17].

Among all periodic processes presented in CVS special place is taken by the low-frequency oscillations with the main frequency close to 0.1 Hz (we will further call them as “0.1 Hz oscillations”). These oscillations can be detected in various signals from CVS, such as HRV [2], BPV [2], photoplethysmographic waveform variability (PPGV) [9, 18].

0.1 Hz oscillations in HRV and BPV are commonly associated with both baroreflex and central neural regulation [19-25]. However, an origin of oscillations of such frequency range in photoplethysmogram (PPG) is still the subject of discussion. The most popular opinion is that the low-frequency oscillations in PPGV characterize sympathetic regulation of peripheral vessels tone [26-28]. However, it is important to consider that PPG signal contains information about both peripheral blood flow (including microcirculatory bloodstream) and distal arterial bed [29]. It could explain showings of central neural regulation in peripheral blood flow detected by some authors [30]. Systolic oscillations of blood filling analyzed using PPG signals are similar, but not identical, to the oscillations in BPV in large arteries and to the HRV oscillations [31]. Although these oscillations can reflect similar mechanism of vessels tone regulation [32].

In our previous experiments with the respiration frequency varying linearly from 0.05 Hz to 0.20 Hz, we proved the functional independence of 0.1 Hz oscillations in HRV and PPGV [9, 33]. It was also shown that main frequency of these oscillations is inconstant and can vary in sufficiently wide range [33]. This can be caused by the ability of cardiovascular regulation to adapt to different operating conditions.

We also found out that 0.1 Hz oscillations in HRV and PPGV sometimes demonstrate 1:1 synchronization, i.e. mutual adjustment of instantaneous phases and frequencies [18]. This ensures functional interaction between autonomic control loops of various CVS parts. However, constant synchronization is not required to ensure adequate coupling. For example, total length of synchronization epochs between 0.1 Hz oscillations in resting healthy young subjects vary in sufficiently wide range from 20% to 60% of total observation time [18]. To estimate the strength of synchronization between 0.1 Hz oscillations in HRV and PPGV we proposed a quantitative measure – the total percent of phase synchronization (*S* index). This index potentially has a significant clinical importance, for example, for estimation of cardiovascular risk in patients with myocardial infarction [34, 35], ensuring the effectiveness and safety of medical therapy for patients with coronary artery disease and hypertension [36-39], estimation of autonomic dysfunction for perimenopause women [40, 41], and post operational observation of patients after coronary artery bypass grafting [42].

This article describes in detail the method for detection of synchronization between 0.1 Hz oscillations in HRV and PPGV and estimation of total percent of phase synchronization (*S* index). Methodological recommendations for *S* index application in clinical studies are also given.

**Material and Methods**

*Patients*

Clinical adaptability of method for detecting synchronization between 0.1 Hz oscillations in HRV and PPGV was studied in the following test groups (1056 records in total):

- 17 healthy subjects (127 records) (50% women), 26±5 years old;
- 41 patients three weeks after myocardial infarction (34.4% women) (167 records), 55±9 years old;
- 105 hypertensive patients (762 records) (62.9% women), 46±7 years old.

These heterogeneous data allowed us to unify the recommendations for use of the proposed method.

*Biological signals registration*

For every patient in supine position and for some patients in upright position we recorded II standard lead electrocardiogram (ECG) (*Figure* 1a) and right index finger PPG measured with the reflective infrared light photoplethysmogram sensor (*Figure* 1b). Signals registration was conducted by multichannel multiregistrator (electroencephalograph analyzer EEGA-21/26 Encefalan-131-03 with standard sensor package, Medicom MTD, Russia) with 250 Hz sampling rate and 14 bit resolution.

Each record was 10 minutes long. During registration the respiration of patients was spontaneous. Fragments of typical signals are depicted on *Figure *1.

Figure 1. ECG (a) and PPG (b) signals typical for healthy subjects. Both signals are depicted in arbitrary units.

*Data preprocessing *

R-R intervals (RRIs) sequence was extracted from ECG in order to analyze HRV signal (*Figure* 1a). Due to heart rate variations the R-R interval (RRI) measurements are non-equidistant in time. To transform non-equidistant signal to equidistant one, we approximated it by cubic splines and then resampled at regular intervals with 4 Hz sampling rate (in consistence with R.M. Baevsky et al. [43]). Resulting equidistant sequence of RRIs was used for further processing. Extraction of 0.1 Hz oscillations from raw HRV and PPG signals (*Figure* 2) was carried out by means of band pass filtering in the range of 0.06–0.14 Hz [43]. The band pass was selected after studying the dependence of statistical significance level on the band of filtering for the cases of various clinical tasks.

Figure 2. Linear scale Fourier power spectra calculated from HRV (a) and PPG (b) signals recorded from healthy subject. Spectrum spikes close to 0.1 Hz are marked with *f _{v} *symbol. Spikes close to 0.3 Hz are marked with

*f*symbol.

_{r}

*Determination of signals instantaneous phases *

Biological systems are the most complex known objects. Big amount of interacting nonlinear elements (i.e. equations with nonlinear functions should be used in order to mathematically describe them) is typical for these systems. Signals of these systems are exposed to noises of different origin and commonly unsteady (both parameters of the system and statistical properties of the signals themselves vary over time). Because of that it is necessary to develop and use specialized technics in order to study individual behavior of the elements of biological systems and features of their interactions [44, 45].

Historically, the development of methods for processing and analysis of complex signals, including those of biological origin, is the sphere of radiophysics and its subpart – nonlinear dynamics – field of knowlege, oriented on studying of oscillations and wave processes. It is assumed that theories, approaches and methods developed within the sphere of nonlinear dynamics are universal. In other words, results obtained for oscillatory system of one nature (commonly of well studied radiophysical systems and their maththical models), can be used for studying the systems of different nature, including biological objects. Such versatility of methods and approaches is discussed within synergetics concept [46, 47]. One important, yet highly unusual, question connected to the analysis of complex signals, is the analysis of instantaneous phases of oscillations. In this section, the questions concerning instantaneous phases definition and their properties are discussed and several examples are given.

Some of the basic definitions used by radiophysics are illustrated in *Figure *3. In the example displayed in *Figure *3, time dependence of the coordinate of pendulum deviation from equilibrium state – *x*(*t*) (*Figure *3a) is taken as studied signal. In terms of nonlinear dynamics, this value is called dynamical variable, and its discrete representation, appropriate for computer analysis is named as time series, realization, or time realization [44, 45].

Figure 3. Illustration of some basic radiophysics notions by using a pendulum clock: (a) – pendulum oscillations registration with marker placed upon it, leaving a trace on a paper strip dragged underneath the pendulum at fixed pace. (b) – examples of dependence of pendulum deviation *x* on time *t* registered from pendulum clock. Solid line corresponds to the case of pendulum start with the zero initial phase *ϕ(t)*=0, dashed line corresponds to the case of pendulum start with a nonzero initial phase *ϕ(t)*= *ϕ*_{0} (ref. panel (c)); (c) is the illustration of different phases of pendulum oscillations. In this example, the phase has the meaning of the angle between pendulum hanger and vertical line (pendulum equilibrium state).

It was shown that time dependence of clock pendulum coordinate can be precisely described by harmonic function, i.e. trigonometric function of sinus or cosinus, *Equation* (1) [48].

*x*(*t*)* = A *sin(2*π**ft *+ *φ** _{0}*) (1)

Graphic representation of *x*(*t*) signal is depicted in *Figure *3b*. Equation* (1) includes several parameters: *A* is the oscillation amplitude, i.e. pendulum maximal deviation from an equilibrium state, which could be measured, for example, in centimeters; *f=1/T* is the oscillation frequency (linear frequency) (measured in Hz), i.e. amount of oscillation cycles the pendulum performs during 1 second; period *T* (*Figure *3b*) *is the value reciprocal to frequency, it is the time after which the oscillation waveform repeats itself. For convinience the so-called rotational frequency, univalently related to linear frequency accordingly to formula: *ω**=2**π**ft**, is *commonly used.

Argument of the *sin* function in *Equation* (1) is called instantaneous phase and has the meaning of amount of oscillations the system performs by time *t* since the beginning of measurements, *Equation* (2).

*φ*(*t*)* = 2**π**ft + **φ*_{0}* = **ω**t + **φ** _{0}* (2)

Instantaneous phase, as it follows from the definition and *Equation* (2), can never decrease in time, but can increase, generally with variable rate. Instantaneous phase is measured in radians or degrees (which are univalently related). During the full oscillation cycle the instantaneous phase increases by 2*π* radians or 360 degrees that is the same. *φ** _{0}* is the initial phase of oscillations or, in our example, the angle of pendulum initial deviation (

*Figure*3c).

*x*(*t*) signals with different initial phases are compared in *Figure *3b: solid line corresponds to *φ** _{0}*=0 and dashed line corresponds to

*φ*

*≠0. The signals presented in*

_{0}*Figure*3b are phase-shifted one from another.

As follows from *Equation* (2), instantaneous phase of harmonic signal (considered to be standard simplest radiophysical signal) increases linearly with a constant rate, i.e. its dependence on time is linear. Differentiation of *φ*(*t*) in *Equation* (2) gives:

*φ*(*t*)* = **ω* (3)

Instantaneous phase time derivative is called instantaneous frequency (for periodic oscillations it is constant, time independent value: *ω**=*2*π**ft*)*. *Since geometrical meaning of time derivative is the angle of the plot slope with respect to the horizontal axis, then the higher the instantaneous frequency the greater slope angle *φ*(*t*). Instantaneous phases of periodic, constant-frequency oscillations increase with constant angle, proportional to oscillation frequency, *Equation* (3).

Introduction of instantaneous phase can seem unnecessary for the analysis of oscillations with simple waveform, but it appears to be very informative and sensitive value during analysis of complex signals, which allows describing their properties [49]. However, introduction of instantaneous phase itself is a complex problem during analysis of nonperiodic oscillations, since basic notions of period and frequency cannot be strictly introduced. In this case, notions of instantaneous period and instantaneous frequency are commonly used along with their averaged values: characteristic period and characteristic frequency. Generally, in this case it is impossible to strictly define instantaneous phase (similar to *Equation* (2)) and approximate formula and numerical approaches are used. For example, during analysis of ECG signal RRI duration naturally defines instantaneous period of the heartbeat, and averaged value, i.e. heart rate (HR), defines characteristic period.

Instantaneous phase is commonly introduced via cross-section method during analysis of signals with well-defined period or pulse-shaped signals. *Figure *4 illustrates this method of instantaneous phase introduction in the case of ECG signal. Interval between crossings of R-wave front and chosen *x*(*t*)*=s* cross-section (*Figure *4a) is supposed to be instantaneous period. Between two consequent crossings signals instantaneous phase is supposed to linearly increase by 2*π*. As can be seen from *Figure *4b, phase of nonperiodic signal is increasing nonlinearly and provides information about its characteristics. During realization of some methods based on phase dynamics analysis it is convenient to use so-called wrapped instantaneous phase (*Figure *4с). It can be evaluated as remainder of division of unwrapped instantaneous phase (*Figure *4b) by 2*π* (commonly written as *φ*(*t*) mod 2*π*). Unwrapped phase can also be calculated from wrapped phase. In the case of ECG signal, anaysis of instantenous phases can be used for diagnostics and quantitative analysis of arrythmia [50] and for solving others fundamental and practical problems [51-53].

Figure 4. Instantaneous phase for ECG signal: (a) – ECG signal with cross-section. Borders of instantaneous periods are marked with vertical lines; (b) – ECG signal instantaneous phase introduced via cross-section method; (c) – wrapped phase.

Instantenous phases analysis may be a convenient tool for evaluating of spaces between spikes of PQRST-complex in ECG. In this case, their positions, in terms of instantenous phases, became attached not only to absolute time, but also to the phase of cardiac cycle.

As it is seen from *Figure *4a, phase introducted with the above method is uncertain. In particular, postion of cross-section *s *can be chosen variously, and therefore can affect the result. However, such uncertainty is a price one should inevitably pay to apply methods of phase analysis to complex nonperiodic signals. This problem does not have universal solution.

*Figure *5 illustrate some more examples of instantenous phases introduced for various signals with *x*(*t*)=0 cross-section*. *It is seen that phase can linearly increase for the signal with distinctly varying amplitude (*Figure *5b), similar to harmonic signal (*Figure *5a)*. *Along with this, phase introduced for signal with constant amplitude, but with varying frequency, reflects changes in frequency (*Figure *5c). Moreover, in fact the only way to quantitatively describe the properies of singnal from *Figure *5c is to extract its instantaneous phase. Examples represented in *Figure *5(a-c) illustrate the appropriateness and possibility of separated analysis for amplitudes and phases of complex signals.

Figure 5. Examples of instantaneous phases of different signals oscillations: (a) – harmonic signal with phase linearly increasing at the constant rate; (b) – signal with complex waveform and varying amplitude, but linearly increasing phase; (c) – signal amplitude is constant, however its instantaneous frequency is varying in respond to changes in frequency; (d) – example of instantaneous phase introduced for harmonic signal via Hilbert transformation.

Use of intuitive cross-section method for analysis of complex signals results in error, originating from assumption about linearity of instantaneous phase increase during characteristic period. For that reason, other approaches based of phase-plane portrait reconstruction [49, 52, 54-56] are more common. The most common method is based on Hilbert transformation, which is ideal broad-band *–**π**/*2 phase shifter [57, 58]. Every component of Fourier spectrum is getting phase shifted by *–**π**/*2 after Hilbert transformation. For example, if original signal is cosines function, than its Hilbert transformation is sinus with the same frequency and amplitude.

Introduction of phase requires construction of specific plane, first axis of which is signal *x*(*t*) itself, while second one is its Hilbert transformation *h*(*t*). For example, image point (point with (*x*(*t*), *h*(*t*)) coordinates) of the harmonic signals will circle in time around the origin of coordinates. Image points of the complex signals have more complicated trajectories. Signal’s instantaneous phase *φ*(*t*) is introduced as the angle between the axis and line connecting Image point to planes origin of coordinates (point with (*0*, *0*) coordinates) (*Figure* 5d) [49].

For harmonic signals instantaneous phases introduced by any methods (cross-section method, Hilbert transformation and others) match with each other and with *Equation* (2). However, for complex signals different methods of instantaneous phase introduction give different results. In general, extraction of instantaneous phase via Hilbert transformation is considered to be more precise than via cross-section method [49, 59].

*Phase dynamics analysis for coupling systems*

Phase dynamics analysis for interacting systems has special significance. Instantaneous phases of oscillating systems appear to be most sensitive to appearance of weak coupling between systems [44, 45, 60-62]. In this case, with coupling increasing the first changes in system dynamics take place in the phases and only afterwards they can be detected in the amplitudes.

A variety of methods for coupling detection was developed based on phase dynamics analysis. These approaches were successfully tested on data of biological origin [63-65]. Of special interest for studying are effects of phase synchronization between oscillating systems, which appear when their phases adjust to each other in case of strong enough coupling, whereas the amplitudes can be chaotic and uncorrelated [66-70]. An example of coupled systems demonstrating such behavior is represented in *Figure *6. Section 1 corresponds to non-synchronous behavior of *x*(*t*) and *y*(*t*) signals (*Figure *6a). The first vertical dashed line marks the time point when coupling between systems becomes stronger, which leads to appearance of phase synchronization between them, section 2 in *Figure *6. The second vertical dashed line marks the time point when coupling becomes weaker, which leads to desynchronization between systems (section 3).

Figure 6. Phase synchronization between oscillating systems: (а) – signals of coupled systems; (b) – *x*(*t*) and *y*(*t*) signal instantaneous phases introduced via Hilbert transformation; (c) – instantaneous phase difference; Vertical dashed lines mark time points of systems coupling strength alterations: left line is the coupling strengthening, which leads to phase synchronization (section 2), right line is the coupling weakening, that causes desynchronization (section 3).

Instantaneous phases *φ** _{x}*(

*t*) and

*φ*

*(*

_{y}*t*) calculated from

*x*(

*t*) and

*y*(

*t*) signals, respectively, via Hilbert transformation are represented in

*Figure*6b

*.*During sections of non-synchronous behavior (sections 1 and 3) phases are increasing independently. Section 2 corresponds to synchronous behavior. It is seen, that at any time slope angels are equal for instantaneous phases during this section (

*Figure*6b), i.e. instantaneous frequencies of the oscillations match each other.

It is convenient to observe phase synchronization by calculating instantaneous phases difference Δ*φ*(*t*). According to definition of phase synchronization [49], the following situation corresponds to synchronous sections:

Δ*φ*(*t*) = |*φ** _{x}*(

*t*)

*-*

*φ*

*(*

_{y}*t*)| <

*C,*(4)

wher *С* is a small constant value, i.e. horizontal plateau on phase difference Δ*φ*(*t*) corresponds to sections of phase synchronization. Δ*φ*(*t*) can perform minor oscillations around nominal horizontal line because of noises inevitably presented in any experimental signal. Such section (time interval between 2 vertical dashed lines) is clearly seen in *Figure *6c. Sections of non-synchronous behavior correspond to increasing curve in graphic Δ*φ*(*t*) (sections 1 and 3).

Therefore, epochs of phase synchronization can be detected as the intervals of experimental signals with difference between their instantaneous phases staying unchanged for several characteristic periods. Commonly, epochs of phase synchronization can be detected visually by putting studied signals on the same plot (see section 2 in *Figure *6a). However, most reliable diagnostics and quantitative analysis is possible only with calculation of instantaneous phase’s difference.

*Analysis of phase synchronization between 0.1 Hz oscillations in HRV and PPGV *

Analysis of phase dynamics and phase synchronization is also widely used during analysis of records of human biopotentials [71-73]. For that purpuse wide variety of methods were developed, including coherence function [74], synchrogram [75, 76], coefficient of phase coherency [49], methods based on wavelet transform [77-80], etc. However, analysis of complex biological objects often requires to consider features of concrete systems, which commonly leads to necessity of making modifications to existing methods or developing new ones. We developed the method for detection of phase synchronization beetwen the loop of baroreflectory regulation of arterial vessels tone and loop of heart rate regulation [9, 33, 61, 81]. According to our research [82], the method specialization on nonstationary data with frequent changes between synchronous and non- synchronous regions, resultsed in its high sensitivity.

The proposed method for detection of phase synchronization includes the following steps of data processing:

- simultaneous registration of ECG and PPG signals;
- extraction of sequence of RRIs from ECG signal;
- calculation of equidistant RRIs using 5 Hz sampling rate and approximation [2, 43];
- extraction of signals produced by studied regulatory loops via filtering of RRIs and PPG signals in 0.05-0.15 Hz band;
- resampling of filtered PPG signal at 5 Hz sampling rate;
- extraction of instantaneous phases from oscillations using Hilbert transformation;
- calculation of instantaneous phases difference Δ
*φ*(*t*).

The data can be registered with any digital device with simultaneous registration of single ECG lead and single lead of finger PPG in reflective or transmitted light. Recording device should provide transmission of registered data to the computer for processing and analysis. Band pass of recording device should be at least 0.05-60 Hz for both leads, sampling rate at least 120 Hz, and quantization bit rate at least 14 bit [83]. Since the characteristic period of studied rhythms is about 10 second, potentially hardware requirements can be lower; however, that demands special investigation. As it was shown in our previous papers, phase synchronization can be detected from single PPG signal by extracting heart rate information from it using original method [84, 85].

Detection of synchronous epochs was carried out via original automated procedure. Due to presence of noises of various origins in experimental signals of studied systems, it is not trivial to detect gently sloping regions of instantaneous phase difference. The method based on Δ*φ*(*t*) linear approximation in sliding windows was proposed to lower the effects of fluctuations caused by noise. For that purpose the line equation was matched to Δ*φ*(*t*) time series in *b*(*с*) width sliding window via least square method [86]. The procedure consists in choosing of initially unknown coefficients *α* and *β* for abstract line equation: *z*(*t*)* = **α**t *+ *β* to fit this line between the points of region with width *b* of noisy experimental time series Δ*φ*(*t*) (thin line in *Figure* 7). Coefficient *α** = *tg*γ *corresponds to the slope angle *γ* of *z*(*t*) line in relation to horizontal axis (horizontal dashed line in *Figure* 7). *α* is the time derivative of instantaneous phase, i.e. instantaneous frequency mismatch at a time (3). If *z* line is horizontal, as it supposed to be in synchronous sections, then *α* = *γ* = 0. Because of the noises and unavoidable error, caused by finite calculations accuracy, approximating line *z* commonly would be slightly not horizontal during processing of experimental data. It was taken into account during in the proposed method for synchronization detection.

Figure 7. Illustration of the method for automated detection of phase synchronization epochs from the signal of instantaneous phases difference Δ*ϕ*(*t*). Borders of phase synchronization section are marked with vertical dashed lines.

Because of nonstationarity and complexity of signals produced by studied systems, short gently sloping regions of phase difference can appear as a consequence of random short-term matches between instantaneous frequencies of not even coupled systems. To exclude these sections from analysis, the minimal length of the detected synchronous epochs *l*(*c*) was restricted to about 2 characteristic periods of oscillations.

We carried out a special investigation in order to increase the sensitivity and specificity of the proposed method for detection of phase synchronization beetwen the loop of baroreflectory regulation of arterial vessels tone and loop of heart rate regulation, by defining the values of its free parameters. Defined values of parameters were the following: *b*=13 (c), |*α*|≤0.01, *l*=16 (c) [81].

Thus, the algorithm was proposed to detect epochs of phase synchronisation between the studied systems by analysing the time seriece of instantaneous phases difference:

- To define slope α of straight line
*z*, its equation is matched by least square method to the region of*b*seconds, commonly called sliding window. - Slope
*α*is calculated for the next region of Δ*φ*(*t*) with the same length (sliding window), but shifted by 1 discrete count (minimal possible time step). The procedure is repeated multiple times, i.e. sliding window goes through all Δ*φ*(*t*) signal and for each step new value of*α*is calculated*.* - If the approximated line
*z*stays close to horizontal one, i.e. the absolute value of slope angle |*α*| is smaller than specific value 0.01, in each forthcoming sliding window for a region with duration of at least*l*seconds, then this region is considered to be a region of phase synchronization between the studied systems.

Index *S* is the total percent of phase synchronization that was proposed as quantitative measure to characterize the strength of phase synchronization [9, 18]. It is calculated as the sum of all detected synchronization epochs divided by the total length of a record, and then expressed as a percentage.

*Statistical significance analysis of total percent of phase synchronization *

Random fluctuations of frequency and phase are common for complex nonstationary signals of biological origin. Accidental coincidence between the frequencies and phases can appear even for not coupled signals and can be falsely detected as sections of phase synchronization. Such events will decrease the value of calculated index *S*, decreasing the specificity of results [87]. False detections are much more likely to appear in the case when characteristic frequencies of studied signals are close to each other.

Therefore, during analysis of experimental data it is important to evaluate the possibility of index *S* to take a certain value due to random fluctuations of signals and not as a result of specific coupling dynamics between the studied systems. Such procedure is called statistical significance analysis of the results. Approach based on generation of a group of artificially-synthesized surrogate data was used to evaluate the statistical significance in a majority of our papers [88, 89]. These data reproduce some statistical properties of original signals. However, all couplings that could take place between them are intentionally destroyed.

Most commonly we used Amplitude Adjusted Fourier Transform (AAFT) surrogate data to test the statistical hypothesis about uncoupled linear systems [90, 91].

The method is based on estimation of registered signal periodogram (Fourier power spectrum of time series estimated without time averaging), with further randomization of its harmonics phases, while maintaining their power. *N* pairs of surrogate signals are generated (commonly 100 or 1000 pairs) by randomizing the phases of harmonics in Fourier spectrum and further calculation of inverse Fourier transform. While the length and spectral properties of resulting signals match with the corresponding properties of experimental signals all couplings between surrogate signals are intentionally destroyed. *S* value is calculated for every pair of surrogate time series from the generated set. By doing so, nonzero values of *S* can only appear as a result of random matches between instantaneous phases of surrogate signals. If *S* value calculated from experiment is greater than the value calculated from surrogate data, then statistical hypothesis about uncoupled signals is considered to be invalid. I.e. calculated *S* value is considered to be not accidental, but defined by coupling of the systems. Calculated *S* value is consider to invalidate the statistical hypothesis about uncoupled signals with 0.95 probability, if experimental *S* value is greater than at least 95% of indexes, calculated from surrogate data. It also can be stated that *S* index value is statistically significant at the level *p*>0.05. For example, if set of surrogate signals is *N*=100 pairs, then for experimental *S* value to be statistically significant at the level *p*>0.05, it is needed to be greater than at least 95 values calculated from surrogate data set. Significance level *p*>0.05 means that at most 5% of experimental indexes, are defined by random fluctuations and not by the presence of phase synchronization [92].

Statistical significance level is typically 0.05 or 0.01 (when no more than 1% of random conclusions is acceptable), and its selection is defined by formulation of concrete research problem.

Insignificant result does not state that synchronization is absent, but that chosen analysis method is unable to make reliable conclusion about its presence between concrete pair of time series (which could be caused by high level of noise, signals distortions, insufficient length of the record and other factors).

Other methods of surrogate data generation also can be used to control the statistical significance of the results. For example, the approach based on mixing of instantaneous phases by 2*π* intervals [93] is widely used or the method based on random choice of experimental realizations from single group of patients (i.e. the first signal from one patient and the second from another) [94]. Other methods that test various statistical hypotheses and based on various prior guesses about data properties [87-89] are also employed.

**Results**

With the proposed method the total percent *S* of phase synchronization between 0.1 Hz oscillations in HRV and PPGV was calculated for all subjects. Example of processing and analysis of data from healthy subject is presented in *Figure *8.

Figure 8 (part 1). Results of application of the proposed method of phase synchronization diagnostics to 25 year healthy man. (a) – RRI variability, i.e. sequence of R-R intervals length (dots) and its approximation via cubic splines for obtaining equidistant RRIs. (b) – PPG. (c) and (d) – equidistant RRIs and PPG power spectra, respectively. Vertical dashed line marks 0.1 Hz that is the center of the signal filtering band pass. **(e) and (f) – signals x(t) and y(t) – RRIs and PPG signals filtered in 0.05-0.15 Hz band (solid line), dashed line represents their Hilbert transformations. (g) – instantaneous phases of the signals: x(t) is shown by solid line and y(t) is shown by dashed line. (h) – difference between instantaneous phases of x(t) and y(t) signals. Detected sections of phase synchronization are marked with brackets.**

Distributions of *S* index values for studied groups of healthy subjects and subjects at three weeks after myocardial infarction are presented in *Figure *9a*. *It appears that *S* index in average is higher for healthy people (33.3±16.2% represented as М±σ), than for patients with myocardial infarction (15.7±9.4%), which correlates well with our earlier results [18]. Therefore, the proposed method provides adequate separation of groups with different functional status of CVS (in particular, healthy people and myocardial infarction patients).

**Figure 9. Distribution histograms for total percent of phase synchronization S between 0.1 Hz oscillations in HRV and PPGV of healthy subjects (solid line) and patients at 3 weeks after the myocardial infarction (dashed line). (a) – estimation for the entire group; (b) – estimation for the group with significant results (p<0.05).**

Tests with surrogate data have shown that only about a half of sunchronization index *S* values are statisticaly significant (*p*<0.05) for both healthy subjects and patients with myocardial infarction. It appears that statistical significance control of synchronization analysis allows us to increase the method sensivity. It is seen from comparison of *Figure *9a and *Figure *9b that selection of significant results improves patients health status clasterisation (healthy people and patients with myocardial infarction): *S* value for the group of healthy subjects is 45.7±12.5% (М±σ) and for the group of myocardial infarction patients is 19.9±12.0%. Therefore, statistical significance control of quantitative measure of synchronization between 0.1 Hz oscillations in PPGV and HRV improves the effectiveness of the proposed approach.

However, sometimes estimation of statistical significance of *S* index via surrogate data can be neglected, since its low level does not prove the incorrectness of calculated index value, but only indicates the lack of statistical proves to declare its correctness. Empirically evaluated the threshold index value *S** _{c}* can be used as alternative.

Choice of this value is based on the assumption that when the value of index of synchronization between 0.1 Hz oscillations in HRV and PPG is smaller than *S** _{c}*, the autonomic regulation is unable to maintain functional integrity of CVS, lowering its total adaptability. Therefore, regardless of the results of statistical significance evaluation,

*S*values smaller than

*S*

*allow one to estimate the significant desynchronization between 0.1 Hz oscillations.*

_{c}As it is seen from the analysis of all 1056 records, this critical level of synchronization index *S *is about 25% (*Figure *10), and values of *S*≤25% indicate significant desynchronization of 0.1 Hz oscillations. Dependence of calculated conditional probability of getting statistically significant result on the value of *S* index (*Figure* 11) is an additional ground to justify such approach.

Figure 10. Dependence of statistical significance level on the value of total percent of phase synchronization *S*. The region with *S*≤25% values is highlighted with grey color. Horizontal line shows *p*=0.05 significance level.

Figure 11. Conditional probability of making statistically significant conclusion about the presence of synchronization for each chosen value of *S* index.

In studied array of records, the proposed method classified as statistically significant only 21.9% of all results with *S*≤25%, and 67.4% of results with *S*>25%. Therefore, regardless from healthy status, the level of statistical significance *p *will be <0.05 in about 68% of cases (i.e. at level of ±σ) with synchronization index *S*>25%, which is considered to be enough to interpret the study results. This approach is appropriate for the researches where *S* is evaluated retrospectively and repeated registration of biological signals from concrete patients is impossible. However, according to *Figure *11 the threshold value for total percent of synchronization (*S** _{c}*) can differ from 25% for other studied problems. For example, setting

*S*

*level to 20% increased the effectiveness of risk evaluation for myocardial infarction patients, during analysis of*

_{c}*S*index prognostic value [34, 35].

However, inability to evaluate significance of the results for each patient individually heavily reduces the application of this alternative results control method for studying synchronization of 0.1 Hz oscillations in HRV and PPGV. Therefore, in the case of synchronization analysis between 0.1 Hz oscillations in individual patient, ECG and PPG signals should be reregistered to achieve the statistically significant result.

We believe that the proposed method is potentially perspective for clinical application in cardiology and require in-depth study of its diagnostic capabilities in further research.

**Conclusion**

This paper describes in details the method of quantitative evaluation of synchronization strength between 0.1 Hz oscillations in HRV and PPGV, which is based on automatized detection of gently sloping regions of oscillations of instantaneous phase differences with further calculation of total percent of phase synchronization *S*. It is shown that the statistical significance control of the method results using surrogate data increases the sensitivity of the approach, which is especially relevant for individual clinical practice.

The described method was granted by invention patent № 2374986 (RF) from 10 December 2008 (priority since 22 July 2008).

**Conflict of interest: **none declared.

**Acknowledgment**

This work was supported by the President of the Russian Federation, Grant No. MD-4368.2015.7 and by Russian Foundation for Basic Research, Grant No 15-02-03061 and No 16-32-00326.

- Special focus issue on cardiovascular physics edited by N. Wessel, J. Kurths, W. Ditto, and R. Bauernschmitt.
*Chaos*2007; 17(1). - Heart rate variability: standards of measurement, physiological interpretation and clinical use. Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology.
*Circulation*1996; 93(5): 1043–1065. (PMID: 8598068) (doi: 10.1161/01.CIR.93.5.1043) - Parati G, Saul JP, Di Rienzo M, Mancia G. Spectral analysis of blood pressure and heart rate variability in evaluating cardiovascular regulation. A critical appraisal.
*Hypertension*1995; 25: 1276–1286. (PMID: 7768574) (doi: 10.1161/01.HYP.25.6.1276) - Hilsted J. Peripheral blood flow control in diabetes mellitus.
*Acta Physiol Scand Suppl*1991; 603: 47-51. (PMID: 1789129) - Nishihara N, Iwakiri H, Ogiuchi H, Ozaki M. Diagnosis of the autonomic nervous function during sedation by peripheral blood flow and psychogenic perspiration.
*Masui*2003; 52: 128-134. Japanese (PMID: 12649866) - Koryakov VV. Device diagnosis and combined treatment of hyperventilation syndrome.
*Ter Arkh*2012; 84(3): 28-31. Russian (PMID: 22708419) - Zhang Q, Patwardhan AR, Knapp CF, Evans JM. Cardiovascular and cardiorespiratory phase synchronization in normovolemic and hypovolemic humans.
*Eur J Appl Physiol*2015; 115: 417-427. (PMID: 25344797) (doi: 10.1007/s00421-014-3017-4) - Prokhorov MD, Ponomarenko VI, Gridnev VI, et al. Synchronization between main rhythmic processes in the human cardiovascular system.
*Phys Rev E*2003; 68: 041913. (PMID: 14682979) - Karavaev AS, Prokhorov MD, Ponomarenko VI, et al. Synchronization of low-frequency oscillations in the human cardiovascular system.
*Chaos*2009; 19: 033112. (doi: 10.1063/1.3187794) (PMID: 19791992) - Kirilina TV, Krasnikov GV, Tankanag AV, et al. Spatial synchronization of the blood flow oscillations in human skin microcirculation.
*Regional Haemodynamics and Microcirculation*2009; 8(3): 32-36. Russian - Kiselev AR, Gridnev VI. Oscillatory processes in vegetative regulation of cardiovascular system.
*Saratov Journal of Medical Scientific Research*2011; 7: 34-39. Russian (Link) - Orini M, Bailón R, Mainardi LT, Laguna P. Time-frequency phase differences and phase locking to char-acterize dynamic interactions between cardiovascular signals.
*Conf Proc IEEE Eng Med Biol Soc*2011; 2011: 4689-4692. (doi: 10.1109/IEMBS.2011.6091161) (PMID: 22255384) - Liao F, Jan YK. Enhanced phase synchronization of blood flow oscillations between heated and adjacent non-heated sacral skin.
*Med Biol Eng Comput*2012; 50: 1059–1070. (PMID: 22936012) (PMCID: PMC3469726) (doi: 10.1007/s11517-012-0948-y) - Seidel H, Herzel H. Bifurcations in a nonlinear model of the baroreceptor-cardiac reflex.
*Physica D: Nonlinear Phenomena*1998; 115: 145-160. (doi: 10.1016/S0167-2789(97)00229-7) - Ottesen JT. Modelling the dynamical baroreflex-feedback control.
*Mathematical and Computer Modelling*2000; 31: 167-173. (doi: 10.1016/S0895-7177(00)00035-2) - Kotani K, Struzik ZR, Takamasu K, et al. Model for complex heart rate dynamics in health and disease.
*Physical Review E*2005; 72: 041904. (doi: 10.1103/PhysRevE.72.041904) - Silvani A, Magosso E, Bastianini S, et al. Mathematical modeling of cardiovascular coupling: Central autonomic commands and baroreflex control.
*Auton Neurosci*2011; 162: 66-71. (PMID: 21550860) (doi: 10.1016/j.autneu.2011.04.003) - Kiselev AR, Bespyatov AB, Posnenkova OM, et al. Internal synchronization of the main 0.1-Hz rhythms in the autonomic control of the cardiovascular system.
*Human Physiology*2007; 33: 188–193. (doi: 10.1134/S0362119707020089) - deBoer RW, Karemaker JM, Strackee J. Hemodynamic fluctuations and baroreflex sensitivity in humans: a beat-to-beat model.
*Am J Physiol*1987; 253: H680- H689. (PMID: 3631301) - Montano N, Gnecchi-Ruscone T, Porta A, et al. Presence of vasomotor and respiratory rhythms in the discharge of single medullary neurons involved in the regulation of cardiovascular system.
*J Auton Nerv Syst*1996; 57: 116-122. (PMID: 8867094) - Cooley RL, Montano N, Cogliati C, et al. Evidence for a central origin of the low-frequency oscillation in RR-interval variability.
*Circulation*1998; 98: 556-561. (PMID: 9714113) (doi: 10.1161/01.CIR.98.6.556) - Whittam AM, Claytont RH, Lord SW, et al. Heart rate and blood pressure variability in normal subjects com-pared with data from beat-to-beat models developed from de Boer’s model of the cardiovascular system.
*Physiol Meas*2000; 21: 305-318. (PMID: 10847197) (doi: 10.1088/0967-3334/21/2/310) - Cohen MA, Taylor JA. Short-term cardiovascular oscillations in man: Measuring and modeling the physiologies.
*J Physiol (London)*2002; 542: 669-683. (PMID: 12154170) (doi: 10.1113/jphysiol.2002.017483) (PMCID: PMC2290446) - Malpas S. Neural influences on cardiovascular variability: Possibilities and pitfalls.
*Am J Physiol Heart Circ Physiol*2002; 282: H6-H20. (PMID: 11748042) - Gridnev VI, Kiselev AR, Kotel’nikova EV, et al. Influence of external periodic stimuli on heart rate variability in healthy subjects and in coronary heart disease patients.
*Human Physiology*2006; 32: 565–573. (doi: 10.1134/S0362119706050100) - Bernardi L, Radaelli A, Solda PL, et al. Autonomic control of skin microvessels: assessment by power spectrum of photoplethysmographic waves.
*Clin Sci (London)*1996; 90: 345–355. (PMID: 8665771) (doi: 10.1042/cs0900345) - Middleton PM, Chan GS, Steel E, et al. Fingertip photoplethysmographic waveform variability and systemic vascular resistance in intensive care unit patients.
*Med Biol Eng Comput*2011; 49: 859–866. (doi: 10.1007/s11517-011-0749-8) (PMID: 21340639) - Middleton PM, Tang CH, Chan GS, et al. Peripheral photoplethysmography variability analysis of sepsis patients.
*Med Biol Eng Comput*2011; 49: 337–347. (doi: 10.1007/s11517-010-0713-z) (PMID: 21153887) - Rhee S, Yang BH, Asada H. Theoretical evaluation of the influence of displacement on finger photoplethysmography for wearable health monitoring sensors. In: Symposium on Dynamics, Control, and Design of Biomechanical Systems ASME International Mechanical Engineering Congress and Exposition, Nashville, Tennessee, November 14-19, 1999.
- Bernardi L, Hayoz D, Wenzel R, et al. Synchronous and baroceptor-sensitive oscillations in skin microcirculation: evidence for central autonomic control.
*Am J Physiol Heart Circ Physiol*1997; 273: 1867-1878. (PMID: 9362255) - González H, Infante O, Lerma C. Response to active standing of heart beat interval, systolic blood volume and systolic blood pressure: recurrence plot analysis Translational Recurrences.
*Springer Proceedings in Mathematics & Statistics*2014; 103: 109-123. (doi: 10.1007/978-3-319-09531-8_7) - Millasseau SC, Guigui FG, Kelly RP, et al. Noninvasive assessment of the digital volume pulse. Comparison with the peripheral pressure pulse.
*Hypertension*2000; 36: 952–956. (PMID: 11116106) (doi: 10.1161/01.HYP.36.6.952) - Karavaev AS, Kiselev AR, Gridnev VI, et al. Phase and frequency locking of 0.1-Hz oscillations in heart rate and baroreflex control of blood pressure by breathing of linearly varying frequency as determined in healthy subjects.
*Human Physiology*2013; 39(4): 416–425. (doi: 10.1134/S0362119713010040) - Kiselev AR, Gridnev VI, Prokhorov MD, et al. Evaluation of 5-year risk of cardiovascular events in patients after acute myocardial infarction using synchronization of 0.1-Hz rhythms in cardiovascular system.
*Ann Noninvasive Electrocardiol*2012; 17(3): 204–213. (PMID: 22816539) (doi: 10.1111/j.1542-474X.2012.00514.x) - Kiselev AR, Gridnev VI, Karavaev AS, et al. Evaluation of five-year risk of lethal outcome and development of cardiovascular disorders in patients with acute myocardial infarction on basis of 0.1-Hz rhythms synchronization in cardiovascular system.
*Saratov Journal of Medical Scientific Research*2010; 6: 328-338. Russian (Link) - Kiselev AR, Gridnev VI, Posnenkova OM, et al. Assessment of dynamics of the autonomic cardiovascular system regulation based on low-frequency rhythm synchronization in patients with ischemic heart diseases complicated by myocardial infarction treated with metoprolol.
*Ter Arkh*2007; 79(4): 23-31. Russian (PMID: 17564014) - Kiselev AR, Gridnev VI, Karavaev AS, et al. Effect of carvedilol and metoprolol on vegetative regulation of heart and microcirculation in patients with hypertension and high body mass.
*Rational Pharmacother Card*2009; (3): 55-61. Russian (doi: 10.20996/1819-6446-2009-5-3-55-61) - Kiselev AR, Gridnev VI, Prokhorov MD, et al. Selection of optimal dose of beta-blocker treatment in myocardial infarction patients based on changes in synchronization between 0.1 Hz oscillations in heart rate and peripheral microcirculation.
*J Cardiovasc Med*2012; 13: 491-498. (doi: 10.2459/JCM.0b013e3283512199) (PMID: 22343262) - Kiselev AR, Gridnev VI, Prokhorov MD, et al. Effects of antihypertensive treatment on cardiovascular autonomic control.
*Anadolu Kardiyol Derg*2014; 14: 701-710. (PMID: 25188759) (doi: 10.5152/akd.2014.5107) - Neyfel’d IV, Kiselev AR, Karavaev AS, et al. Peculiarities of indexes of autonomic regulation of blood circulation and heart rate variability in perimenopausal women.
*Annaly Aritmologii*2014; 11(2): 98-108. Russian (doi: 10.15275/annaritmol.2014.2.4) - Neufeld IW, Kiselev AR, Karavaev AS, et al. Autonomic control of cardiovascular system in pre- and postmenopausal women: a cross-sectional study.
*J Turk Ger Gynecol Assoc*2015; 16: 11-20. (doi: 10.5152/jtgga.2015.15201) (PMID: 25788843) - Shvartz VA, Kiselev AR, Karavaev AS, et al. Study of the dynamics of nonlinear properties of baroreflex regulation of circulation in patients undergoing surgery coronary artery bypass grafting (study design).
*Bulletin of Medical Internet Conferences*2014; 4(9): 1042-1048. Russian (Link) - Baevskiy RM, Ivanov GG, Chireykin LV, et al. The analysis of heart rate variability using different electrocardiographic systems.
*Vestnik Aritmologii*2002; (24): 65-86. Russian - Bezruchko BP, Smirnov DA. Mathematical modeling and chaotic time series. Saratov, Russia: GosUNTs "Kolledzh", 2005; 299 p. Russian
- Bezruchko BP, Smirnov DA. Extracting knowledge from time series (an introduction to nonlinear empirical modeling). Berlin, Germany: Springer-Verlag Berlin Heidelberg, 2010; 405 p. (doi: 10.1007/978-3-642-12601-7)
- Malinetskiy G, Potapov A, Podlazov A. Nonlinear Dynamics: approaches, results, hope. Moscow, Russia: Librokom, 2011; 280 p. Russia
- Trubetskov DI. Introduction to synergetics. Oscillations and waves. Moscow, Russia: Librokom, 2011; 224 p. Russia
- Trubetskov DI, Rozhnev AG. Linear oscillations and waves. Moscow, Russia: Fizmatlit, 2001; 416 p. Russian
- Pikovsky A., Rosenblum M., Kurths J. Synchronization: а universal concept in nonlinear sciences. USA: Cambridge University Press, 2001; 441 p.
- Umapathy K, Nair K, Masse S, et al. Phase mapping of cardiac fibrillation.
*Circ Arrhythm Electrophysiol*2010; 3: 105-114. (PMID: 20160178) (doi: 10.1161/CIRCEP.110.853804) - Latka M, Turalska M, Glaubic-Latka M, et al. Phase dynamics in cerebral autoregulation. Am J Physiol Heart Circ Physiol 2005; 289: H2272- H2279. (PMID: 16024579) (doi: 10.1152/ajpheart.01307.2004)
- Kantz H, Kurths J, Mayer-Kress G. Nonlinear analysis of physiological data. Berlin, Germany: Springer Berlin Heidelberg, 2012; 344 p. (doi: 10.1007/978-3-642-71949-3)
- Laughner JI, Ng FS, Sulkin MS, et al. Processing and analysis of cardiac optical mapping data obtained with potentiometric dyes. Am J Physiol Heart Circ Physiol 2012; 303: H753–H765. (PMID: 22821993) (doi: 10.1152/ajpheart.00404.2012)
- Cohen L. Time-frequency analysis electrical engineering signal processing. Prentice Hall PTR, 1995; 299 p.
- Zhang Y, Wang S, Ji G, Dong Z. An improved quality guided phase unwrapping method and its applications to MRI.
*Progress In Electromagnetics Research*2014; 145: 273-286. (doi: 10.2528/PIER14021005) - Taner MT, Koehler F, Sheriff RE. Complex seismic trace analysis.
*Geophysics*1979; 44: 1041-1063. (doi: 10.1190/1.1440994) - Ifeachor EC, Jervis BW. Digital signal processing: a practical approach (2nd Edition). Prentice Hall, 2002; 933 p.
- Gabor D. Theory of communication. Part 1: The analysis of information.
*Journal of the Institution of Electrical Engineers - Part III: Radio and Communication Engineering*1946; 93(26): 429-441. (doi: 10.1049/ji-3-2.1946.0074) - Khovanov IA, Khovanova NA. Methods of processing time series. Saratov, Russia: GosUNTs "Kolledzh", 2001; 120 p. Russian
- Tass P, Smirnov D, Karavaev A, et al. The causal relationship between subcortical local field potential oscillations and Parkinsonian resting tremor.
*J Neural Eng*2010; 7(1): 16009. (PMID: 20083863) (doi: 10.1088/1741-2560/7/1/016009) - Kiselev AR, Khorev VS, Gridnev VI, et al. Interaction of 0.1-Hz oscillations in heart rate variability and distal blood flow variability.
*Human Physiology*2012; 38(3): 303–309. (doi: 10.1134/S0362119712020107) - Mase M, Glass L, Disertori M, Ravelli F. The AV synchrogram: a novel approach to quantify atrioventricular coupling during atrial arrhythmias.
*Biomedical Signal Processing and Control*2013; 8(6): 1008–1016. (doi: 10.1016/j.bspc.2013.01.004) - Tallon-Baudry C, Bertrand O, Fischer C. Oscillatory synchrony between human extrastriate areas during visual short-term memory maintenance.
*J Neurosci*2001; 21(20): RC177. (PMID: 11588207) - Rosenblum MG, Pikovsky AS, Kurths J, et al. Phase synchronization: from theory to data analysis. In: Neuro-informatics, edited by F. Moss and S. Gielen. Handbook of Biological Physics. Vol. 4. New York: Elsevier Science, 2000: 279-321.
- Quian Quiroga R, Kreuz T, Grassberger P. Event synchronization: a simple and fast method to measure synchronicity and time delay patterns.
*Phys Rev E*2002; 66: 041904. (doi: 10.1103/PhysRevE.66.041904) - Pikovskii AS. Synchronization and stochastization of array of self-excited oscillators by external noise.
*Radiophys Quantum Electron*1984; 27(5): 390-395. (doi: 10.1007/BF01044784) - Pikovskii A. On the interaction of strange attractors.
*Zeitschrift für Physik B Condensed Matter*1984; 55(2): 149-154. (doi: 10.1007/BF01420567) - Afraimovich VS, Verichev NN, Rabinovich MI. Stochastic synchronization of oscillations in dissipative systems.
*Radiophys Quantum Electron*1986; 29(9): 795-803. (doi: 10.1007/BF01034476) - Aranson IS, Rul'kov NF. Nontrivial structure of synchronization zones in multidimensional systems.
*Phys Lett A*1989; 139(8): 375-378. (doi: 10.1016/0375-9601(89)90581-1) - Pecora LM, Carroll TL. Synchronization of chaotic systems.
*Chaos*2015; 25: 097611. (doi: 10.1063/1.4917383) - Rosenblum MG, Kurths J, Pikovsky A, et al. Synchronization in noisy systems and cardiorespiratory interaction. Engineering in Medicine and Biology 1998; 17(6): 46-53. (doi: 10.1109/51.731320)
- Tass P, Rosenblum MG, Weule J, et al. Detection of n:m phase locking from noisy data: application to magnetoencephalography. Phys Rev Lett 1998; 81: 3291-3294. (doi: 10.1103/PhysRevLett.81.3291)
- Wu MC, Hu CK. Empirical mode decomposition and synchrogram approach to cardiorespiratory synchronization. Physical Review E 2006; 73: 051917. (doi: 10.1103/PhysRevE.73.051917)
- White LB, Boashash B. Cross spectral analysis of nonstationary processes.
*IEEE Transactions on Information Theory*1990; 36: 830-835. (doi: 10.1109/18.53742) - Schafer C, Rosenblum MG, Abel HH, Kurths J. Synchronization in the human cardiorespiratory system.
*Physical Review E*1999; 60: 857–870. (PMID: 11969830) (doi: 10.1103/PhysRevE.60.857) - Bartsch R, Kantelhardt JW, Penzel T, Havlin S. Experimental evidence for phase synchronization transitions in the human cardiorespiratory system. Phys Rev Lett 2007; 98: 054102. (PMID: 17358862) (doi: 10.1103/PhysRevLett.98.054102)
- Hramov AE, Koronovskii AA, Ponomarenko VI, Prokhorov MD. Detection of synchronization from univariate data using wavelet transform.
*Phys Rev E*2007; 75: 056207. (doi: 10.1103/PhysRevE.75.056207) - Hramov AE, Koronovsky AA, Ponomarenko VI, Prokhorov MD. Detecting synchronization of self-sustained oscillators by external driving with varying frequency.
*Phys Rev E*2006; 73: 026208. (doi: 10.1103/PhysRevE.73.026208) - Koronovskii AA, Ponomarenko VI, Prokhorov MD, Hramov AE. Diagnostics of the synchronization of self-oscillatory systems by an external force with varying frequency with the use of wavelet analysis.
*Journal of Communications Technology and Electronics*2007; 52(5): 544-554. (doi: 10.1134/S1064226907050087) - Koronovskii AA, Hramov AE, Ponomarenko VI, Prokhorov MD. Method of studying the synchronization of self-sustained oscillations using continuous wavelet analysis of univariant data.
*Technical Physics. The Russian Journal of Applied Physics*2007; 52(9): 1106-1116. (doi: 10.1134/S1063784207090022) - Bezruchko BP, Gridnev VI, Karavaev AS, et al. Technique of investigation of synchronization between oscillatory processes with the frequency of 0.1 Hz in the human cardiovascular system.
*Izvestiya VUZ. Applied Nonlinear Dynamics*2009; 17(6): 44-56. Russian - Borovkova EI, Karavaev AS, Ponomarenko VI, Prokhorov MD. Comparison of methods for phase synchronization diagnostics from test data modeling nonstationary signals of biological nature.
*Izvestiya of Saratov University. New series. Series Physics*2015; 15(3): 36-42. (doi: 10.18500/1817-3020-2015-15-3-36-42) - Kulminsky DD, Borovkova EI, Khorev VS, Mironov SA. Development of the device for daily monitoring of state of cardiovascular system by analyzing its rhythms synchronization.
*Bulletin of Medical Internet conferences*2014; 4(7): 962-966. Russian - Kulminsky DD, Astakhov OV, Borovkova EI, Kiselev AR. Assessment of cardiovascular system based on the evaluation of its rhythms sinchronization on univariate signal photoplethysmogram. In: Proceeding of All-Russian youth scientific conference "Actual problems of biomedical engineering." Saratov, 2013: 330-335. Russian
- Borovkova EI, Karavaev AS, Kiselev AR, et al. Method for diagnostics of synchronization of 0.1 Hz rhythms of cardiovascular system autonomic regulation in real time.
*Annaly Aritmologii*2014; 11(2): 129-136. Russian (doi: 10.15275/annaritmol.2014.2.7) - Linnik YV. The method of least squares, and the foundations of mathematics and statistical processing of observations of the theory. Moscow, Russia: Fizmatlit Publ., 1958; 336 p. Russian
- Trukhacheva NV. Mathematical statistics in biomedical research using Statistica package. Moscow, Russia: Geotar Media, 2012; 379 p. Russian
- Nait-Ali A. Advanced Biosignal Processing. Springer Science & Business Media, 2009; 378 p. (doi: 10.1007/978-3-540-89506-0)
- Chan K-S, Tong H. Chaos: A Statistical Perspective. Springer Series in Statistics. Springer Science & Business Media, 2013; 300 p. (doi: 10.1007/978-1-4757-3464-5)
- Theiler J, Eubank S, Longtin A, et al. Testing for nonlinearity in time series: the method of surrogate data.
*Physica D*1992; 58: 77-94. (doi: 10.1016/0167-2789(92)90102-S) - Schreiber T, Schmitz A. Improved surrogate data for nonlinearity tests.
*Phys Rev Lett*1996; 77: 635–638. (doi: 10.1103/PhysRevLett.77.635) (PMID: 10062864) - Ayvazyan SA. Applied Statistics. Moscow, Russia: Ripol Klassik Publ., 1983; 384 p. Russian
- Brea J, Russell DF, Neiman AB. Measuring direction in the coupling of biological oscillators: A case study for electroreceptors of paddlefish.
*Chaos*2006; 16: 026111. (doi: 10.1063/1.2201466) (PMID: 16822043) - Toledo E, Rosenblum MG, Kurths J, Akselrod S. Cardiorespiratory synchronization: Is it a real phenomenon? In: Computers in Cardiology. A. Murray and S. Swiryn eds. Hannover: IEEE Computer Society Press, 1999: 237-240. (doi: 10.1109/CIC.1999.825950)

*Received 12 January 2016, Accepted 18 February 2016*

© 2016, Kiselev A.R., Karavaev A.S., Gridnev V.I., Prokhorov M.D., Ponomarenko V.I., Borovkova E.I., Shvartz V.A., Ishbulatov Y.M., Posnenkova O.M., Bezruchko B.P.

© 2016, Russian Open Medical Journal

*Correspondence to* Anton R. Kiselev. Address: Research Institute of Cardiology, 141, Chernyshevsky str., Saratov, 410028, Russia. Phone: +7 (8452) 201 899. E-mail: kiselev@cardio-it.ru