A qualitatively new step in the medical development is an active implementation of the methods of diagnosing functional disturbances in the systems’ interaction, including the nonlinear analysis methods. This is a promising trend since it allows detecting deviations in the internal system operation even before the appearance of morphological changes. In the papers [1, 2], the researches investigated the synchronization between the main heart rhythm and the breathing in children and adults and showed the fundamental importance of studying the circulatory elements synchronization in understanding how the human cardiovascular system is structured. The cardiorespiratory interaction was studied in detail in the works by Stefanovska and Bračič , McGuinness et al. , and Prokhorov et al. . Kralemann et al.  also proposed the method of obtaining the phase-coupling functions which describe the cardiorespiratory interactions and the phase response curve in healthy humans.
In our previous studies we have shown that the low-frequency (LF) oscillations in the heart rate variability (HRV) and the finger photoplethysmographic waveform variability (PPGV), which characterize the state of systemic autonomic regulation of blood pressure, can synchronize with each other during spontaneous breathing [7, 8]. It is known that the LF oscillations in the HRV are associated with the baroreflex regulation of blood pressure [9-11], and similar oscillations in the PPGV are associated with the sympathetic control of vascular tone [12-15]. The total length of the synchronous intervals is associated with the clinical status of a person. In particular, the LF oscillations are normally synchronous with each other for a considerable part of the time, and the total length of time synchronization intervals significantly decreases in patients with cardiovascular diseases [7, 8, 16-18]. From the physiological viewpoint, this synchronization is the result of the adequate functional interaction of the cardiovascular system elements provided by cardiovascular autonomic control. We proposed an indicator to estimate the quality of LF oscillation synchronization in the HRV and the PPGV – the total percentage of synchronization, called index S .
The results of specialized experiments, carried out in this paper, support the fact that in our studies we observe exactly the phase synchronization of two interacting self-oscillating systems . A typical problem in the study of the interaction of complex systems by their signals is the need to avoid false conclusions about communication and synchronization, which are caused by parasitic linear mixing of the detected signals, and not by their dynamic interaction. However, in our studies such a problem is excluded by the method of data registration. From the photoplethysmogram signal, we extracted 0.1 Hz components that reflect the processes of amplitude modulation. These oscillations are observed at the frequency of about 0.1 Hz in the spectrum of the PPG power. To extract these components, we apply direct filtering of the PPG signal. However, from the electrocardiogram signal, we extract 0.1 Hz components that reflect the frequency modulation of this signal. These oscillations are observed in the frequency domain in the form of satellites of the main heart rhythm. It means that if the frequency of the main heart rhythm is 1 Hz, then its frequency modulation signal with a frequency of 0.1 Hz will be observed at a frequency of 0.9 Hz, which excludes its linear parasitic mixing with the components extracted from the PPG. To extract the heart rate signal from the electrocardiogram (ECG), a sequence of RR-intervals (tachogram), is first identified in a standard way. This ensures the transfer of spectral components responsible for frequency modulation to the low-frequency part of the spectrum. These components are further detected by bandpass filtering of the tachogram in the neighborhood of 0.1 Hz.
However, previous studies of the above-mentioned LF oscillation synchronization were carried out on the records of up to 10 minutes in length, due to the ethical and technical features of performing clinical experiments. At the same time, the question of statistical properties of the lengths of the synchronization intervals between LF oscillations in HRV and PPGV and the characteristics of the proposed index S has a fundamental and applied significance. Previously, attempts to analyze the features of the S index dynamics for long-term records were carried out on little statistics in the only paper , which set a rather narrow research goal – to analyze the effect of a two-hour immobilization on cardiovascular autonomic control. A comprehensive study of the statistical properties of the index S has not been accomplished yet.
The purpose of this paper is to study statistical and dynamic properties of the LF oscillation synchronization in the HRV and the finger PPGV in time, which has a fundamental importance in understanding the features of the systemic interaction of the autonomic mechanisms of the circulatory regulation. Likewise, it allows development of recommendations in order to correctly arrange the experimental measurements, which will provide a reliable estimate of the studied parameters, and which will enable a wide applicability of the method, based on evaluating the quality of the LF oscillations synchronization in the HRV and the PPGV, in an experimental and clinical medicine.
Material and Methods
Our study included 42 healthy men aged from 19 to 21 years. The health of all subjects was confirmed by the results of clinical investigation in the Institute of Cardiological Research of Saratov State Medical University n.a. V.I. Razumovsky (Saratov, Russia).
The electrocardiogram (ECG) and the photoplethysmogram (PPG) from the middle finger of the right hand were simultaneously recorded in all subjects during 120 minutes (7,200 seconds) at rest (horizontal position). All signals were sampled at 250 samples per second (sps) and digitized at 14 bits. All subjects were investigated in the afternoon fasting under spontaneous breathing. The signals were recorded in a quiet, temperature-controlled room using a standard electroencephalograph analyzer EEGA-21/26 “Encephalan-131-03” (Medicom MTD Ltd, Taganrog, Russia). An example of recorded signals is shown in Figure 1.
Figure 1. An example of simultaneous records of ECG and finger PPG of subject A.
The analysis of the HRV was carried out according to the method guidelines . First we extracted sequences of RR-intervals from the ECG signal. Then, we obtained equidistant RR-interval time series (tachogram) from a not equidistant sequence of RR-interval time series by using cubic β-splines and resampling with a frequency of 5 Hz. The PPG signal was also resampled with a frequency of 5 Hz. The further stage of processing includes the filtering of tachogram and PPG in the frequency range from 0.06 to 0.14 Hz (using band-pass filter) in order to extract the LF oscillations in the HRV and the PPGV. The frequency range used for filtering is an important parameter of the method. The narrowing of the band makes it possible to improve the definition of instantaneous phases and to reduce the noise level. On the other hand, excessive narrowing of the band will lead to a loss of important information about the dynamics of the system. We carried out special studies on the choice of free parameters of the method, which allowed us to recommend using the band 0.06-0.14 Hz for filtering.
To estimate the synchronization between the LF oscillations in the HRV and the PPGV we used the method we suggested earlier . To calculate the phases of these oscillations we applied the Hilbert transform  and calculated the instantaneous phase difference:
φ = φp – φh ,
where φp is the phase of LF oscillations of PPGV and φh is the phase of LF oscillations of heart rate. The presence of 1:1 phase synchronization is defined by the condition |φ|<const . In this case the phase difference φ(t) fluctuates around a constant value. So the epochs of phase synchronization were detected as the regions of the almost constant phase difference . For automated detection of phase synchronization epochs, we used a procedure based on a linear approximation of instantaneous phase difference in a moving window. Lastly, index S was calculated as the relative time of synchronization between the considered LF oscillations :
where dk is the duration of the kth epochs of synchronization, N is the number of epochs and T is duration of the signal.
We estimated properties of a sequence of durations (L) of the intervals of the phase synchronization between the LF oscillations in the HRV and the PPGV. A typical sequence L(t) is shown in Figure 2a. The horizontal axis is used to specify the time of the beginning of each synchronization interval and the vertical axis – the number of a subject.
To analyze the correlation properties of L(t), we determined the autocorrelation function (ACF) C(t). The ACF was calculated for each subject with the lag t from -2,000 to 2,000 seconds.
The statistical properties of the sequence of synchronization intervals of the LF oscillations in the HRV and the PPGV were studied using P(L) (the probability distribution of L(t)), while analyzing the whole experimental record (Figure 2c).
To estimate the distribution of the index S we divided the experimental records of every person into non-overlapping 10 minute intervals. Records of this duration have been used in our previous studies . For each subject X we estimated a value of the index in every i-th window, then we subtracted a median value (the value was estimated for every subject X) from the sequence . Thus, we received the set of rating values S, which distributed relative to the zero median: . All values ΔS (deviations of the estimation S from the median) obtained from the experimental sample have been used for estimating the probability distribution. The statistics ΔS were calculated: the median of the distribution 0.0%, the average value was 0.2%, the skewness was 0.1, the kurtosis was 1.3.
The effect of the duration of the experimental signal recording on the evaluation properties of the synchronization of the LF oscillations in the HRV and the PPGV was estimated on the basis of the study of how the value of the calculated index S depends on the duration of the time interval w, which was used to calculate S. We established the dependences S(w) for each subject (the value w changed from 10 to 7,200 seconds).
In the paper, we obtained the generalized analysis results of the dependence of the index S fluctuations on the duration of the time interval used for estimation. For this analysis we divided every experimental record into non-overlapping intervals of length w and calculated the evaluation of the standard deviation for every subject X.
In our previous investigation  we obtained evidence of the existence of the index S trend in time during the analysis of the two-hour records of healthy subjects. In this work, we also carried out a diagnosis of the trend of index S in time. For the analysis, we estimated the deviations of the assessment S from the average . The deviations were averaged over the ensemble and the resulting mean values of the deviations DS were drawn on a graph with standard error.
Figure 2. The analysis of the sequence of lengths L of synchronization intervals between the LF oscillations in the HRV and the finger’s PPGV: (a) the horizontal axis indicates time and the vertical axis indicates the number of a subject (white area corresponds to synchronization interval and black area corresponds to non-synchronous state), (b) the autocorrelation function of the sequence L(t), averaged over the entire experimental record of subjects (n=42), (c) the histogram of the probability distribution of the synchronization intervals length, estimated over the entire experimental record of subjects (n=42).
Figure 3. Statistical properties of the total percentage of phase synchronization S: (a) the histogram of the probability distribution of the values of ΔS, calculated for each subject and constructed through the entire experimental record using 10-minute realization intervals. The line indicates the approximation of the experimental distribution by the Gaussian distribution. (b) the dependence of the estimate S on the length of the increasing analyzed interval w for subject B; (c) the dependence of the standard deviation of the estimate S on the length of the analyzed time series, obtained in the analysis of the entire experimental record. The arrow indicates the value of the length of the analyzed time series w=600 seconds, which is most often used in well-known experimental studies [8, 16-18].
The statistical characteristics of the LF oscillations synchronization in the HRV and the PPGV
As a result of the analysis of two-hour records, we revealed an average of about 120 periods of synchronization of the LF oscillations in the HRV and the finger PPGV for every healthy subject. Likewise, in all records we found the presence of sufficiently long continuous synchronization intervals of over 100 seconds in length. The maximum length of such interval for subject A was 118 seconds (Figure 2a) and the maximum length of a continuous synchronization interval among all analyzed data was 15 seconds. The sequences L(t) are visually characterized by irregular noiselike character (Figure 2a).
The ACF averaged over all records (n=42) is shown in Figure 2b. The ACF is decreasing that is typical for a random process. The first 0 of the ACF is observed at a lag of about t@200 seconds (exactly t= 260 seconds).
The histogram of the distribution L(t) shows a characteristic inverse power-law dependence of the form 1/f (Figure 2c), where the long synchronization intervals are less likely than the short ones. The first, second and third quartiles of this distribution were 14.2, 19.0 and 29.0 seconds, respectively. In other words, the shorter synchronization interval is more likely a longer one.
The distribution histogram of the ΔS is shown in Figure 3a. It is demonstrated that the probability distribution of the index S is close to normal, almost symmetrical and has ‘light tails’, which allows one to characterize the statistical properties of the index S values from different experimental records, in particular, the mean value and the standard deviation.
The effect of duration of the experimental signal records on the properties of the index S
A typical example of the S(w) dependence for subject B is shown in Figure 3b. The fluctuations of the estimation of the index S decrease with the increasing duration of the experimental realization w. At large w, S(w) stabilizes and reaches the saturation level. The values of the index S, at which the S (w) dependence stabilizes, vary in different people. The average saturation value for the ensemble is 59.2±5.8% (data presented as mean with standard deviation), 1th, 2nd and 3nd quartiles are 57.5%, 58.8% and 59.6%, respectively.
Figure 3c shows the dependence representing the result of averaging of over the experimental ensemble. The values of w were moved from 60 to 3,000 seconds in increments of 1 second to construct the Figure 3c. The variance of the estimate of S plummets as the duration of experimental realizations increases. The dispersion of the estimates S decreases with the increase of w and the dependence stabilizes. The value of the standard deviation of the studied estimate for the most often used length (600 seconds) of realizations in experimental studies was (600)=9.8%.
We performed the estimation of the index S trend in time. For this purpose, we evaluated the deviations of the estimates S from the mean . Figure 4 shows the dependence of the DS on time for the value of the evaluation interval w=1,000 s. Such length of the window, on the one hand, provides a low level of fluctuations (Figure 3c) and on the other hand, it provides an acceptable time resolution. As it is seen from Figure 4, the value of S shows a small but statistically significant decrease in the middle of the second hour. At the end of the second hour, S increases almost to the original level.
Figure 4. Averaged values of deviations from the mean S, calculated on window of time series of 1,000 s length. Each horizontal axis value on the plot corresponds to the middle of the interval used for the calculation of S. The whiskers designate the standard error.
The obtained results (Figure 2) indicate that the sequence of the synchronization intervals between the LF oscillations in the HRV and the finger PPGV behaves as a weakly correlated random process, demonstrating the 1/f distribution. However, the distribution of the index S, estimated even for 10-minute records, is close to normal (Figure 3a). This behavior of the index S is explained by the conclusions of the central limit theorem, since S is the sum of lengths of the synchronization intervals, which almost do not correlate and have the 1/f distribution. A normality of the S distribution allows, in particular, to use the sample mean and standard deviations to analyze the statistics of S, which is convenient, especially while analyzing small ensembles of short records.
The results of estimating the S value from realizations of increasing length indicate that S tends to a fixed value with the increase of time series length (Figure 3b).
It is important to note that studies of out-patients with cardiovascular diseases are often conducted on rather short records of the biological signals, which is due to various ethical and technical reasons. It is known that the record length affects the accuracy of generally accepted indicators of the HRV [20, 23]. At the same time, short records continue to be used in various studies. The possibility of using short records for estimating time- and frequency-domain parameters of the HRV is confirmed by some authors [24, 25]. It is noted, however, that the parameters associated with low-frequency oscillations are more sensitive to the record length, especially the ones shorter than 10 minutes . The obtained results indicate that a reasonable compromise between the reduction in length of the record and the tendency to obtain the S estimate with a small error is a choice of the record duration of about 600-1,000 seconds. In this case, the standard deviation falls abruptly with regard to shorter length of record and takes the value less than 10%. In the general case, analyzing the records of healthy subjects we can recommend to increase the record length in order to reduce the fluctuations of S (Figure 3c).
The analysis of the index S behavior in time on two-hour records of healthy subjects confirms the earlier researches  on small statistics about the presence of a downward trend of S in the middle of the second hour, explained in the paper by the development of the immobilization stress. In this paper, for the first time, we found an upward trend of the degree of synchronization between the studied circulatory regulation contours at the end of the second hour with respect to the initial values (Figure 4). This result can be caused by anti-stress reactions of the body, which ensure the "habituation" of healthy subjects to immobilization. This observation requires more careful study, in particular, using an analysis of longer records.
In general, the first detailed analysis of the statistical properties of the LF oscillations in the HRV and the finger PPGV shows that although the sequence of the lengths of the phase synchronization intervals behaves like the 1/f random process, the average value of the previously proposed index, named the total percentage of phase synchronization S is, apparently, a stable characteristic inherent in a particular person in the present physiological state.
The fluctuations in the value of the S estimate near the mean value depend largely on the length of the time series. For practical studies we recommend using the experimental records with the length of 10 minutes or more, which will provide a standard deviation of the S estimate less than 10%.
The index S is a consistent and normally distributed estimate of the average time of synchronous behavior of the investigated circuits of the autonomic circulatory regulation. The average saturated value of S, diagnosed from the experimental record analyzed in this study, was 59.2±5.8% (mean ± standard deviation).
This study was funded by the Russian Science Foundation, Grant No. 14-12-00291.
All procedures performed in studies involving human participants were in accordance with the ethical standards of the institutional and national research committee and with the 1964 Helsinki declaration and its later amendments or comparable ethical standards. Design of this study was approved by the Ethics Committee of the Saratov State Medical University n.a. V.I. Razumovsky (Saratov, Russia) in 2017 and written informed consent was obtained from all subjects who participated in this study.
Conflict of interest
The authors declare that they have no conflict of interest.
- Mrowka R, Cimponeriu L, Patzak A, Rosenblum M. Directionality of coupling of physiological subsystems—age related changes of cardiorespiratory interaction during different sleep stages in babies. Am J Physiol Regul Comp Integr Physiol 2003; 145: R1395–R1401. https://dx.doi.org/10.1152/ajpregu.00373.2003.
- Schäfer C, Rosenblum MG, Kurths J, Abel H-H. Heartbeat synchronized with ventilation. Nature 1998; 392: 239–240. https://dx.doi.org/10.1038/32567.
- Stefanovska A, Bračič M. Physics of the human cardiovascular system. Contemp Phys 1999; 40: 31–55. https://dx.doi.org/10.1080/001075199181693.
- McGuinness M, Hong Y, Galletly DC, Larsen PD. Arnold tongues in human cardiorespiratory systems. Chaos 2004; 14: 1–6. https://dx.doi.org/10.1063/1.1620990.
- Prokhorov MD, Ponomarenko VI, Gridnev VI, Bodrov MB, Bespyatov AB. Synchronization between main rhythmic processes in the human cardiovascular system. Phys Rev E 2003; 68: 041913. https://dx.doi.org/10.1103/PhysRevE.68.041913.
- Kralemann B, Frühwirth M, Pikovsky A, Rosenblum M, Kenner T, Schaefer J, Moser M. In vivo cardiac phase response curve elucidates human respiratory heart rate variability. Nat Commun 2013; 4: 2418. https://dx.doi.org/10.1038/ncomms3418.
- Karavaev AS, Prokhorov MD, Ponomarenko VI, Kiselev AR, Gridnev VI, Ruban EI, Bezruchko BP. Synchronization of low-frequency oscillations in the human cardiovascular system. Chaos 2009; 19: 033112. https://dx.doi.org/10.1063/1.3187794.
- Kiselev AR, Karavaev AS, Gridnev VI, Prokhorov MD, Ponomarenko VI, Borovkova EI, et al. Method of estimation of synchronization strength between low-frequency oscillations in heart rate variability and photoplethysmographic waveform variability. Russ Open Med J 2016; 5: e0101. https://dx.doi.org/10.15275/rusomj.2016.0101.
- 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. https://dx.doi.org/10.1152/ajpheart.1987.253.3.H680.
- Goldstein DS, Bentho O, Park MY, Sharabi Y. Low-frequency power of heart rate variability is not a measure of cardiac sympathetic tone but may be a measure of modulation of cardiac autonomic outflows by baroreflexes. Exp Physiol 2011; 96: 1255-1261. https://dx.doi.org/10.1113/expphysiol.2010.056259.
- Moak JP, Goldstein DS, Eldadah BA, Saleem A, Holmes C, Pechnik S, Sharabi Y. Supine low frequency power of heart rate variability reflects baroreflex function, not cardiac sympathetic innervations. Heart Rhythm 2007; 4: 1523–1529. https://dx.doi.org/10.1007/s10286-010-0098-y.
- Bernardi L, Radaelli A, Solda PL, Coats AJ, Reeder M, Calciati A, et al. Autonomic control of skin microvessels: assessment by power spectrum of photoplethysmographic waves. Clin Sci (Lond) 1996; 90: 345–355. https://dx.doi.org/10.1042/cs0900345.
- Colombo R, Marchi A, Borghi B, Fossali T, Rech R, Castelli A, et al. Pulse photoplethysmographic analysis estimates the sympathetic activity directed to heart and vessels. Anesthesiology 2015; 123: 336–345. https://dx.doi.org/10.1097/ALN.0000000000000712.
- Middleton PM, Chan GS, Steel E, Malouf P, Critoph C, Flynn G, et al. Fingertip photoplethysmographic waveform variability and systemic vascular resistance in intensive care unit patients. Med Biol Eng Comput 2011; 49: 859–866. https://dx.doi.org/10.1007/s11517-011-0749-8.
- Nitzan M, de Boer H, Turivnenko S, Babchenko A, Sapoznikov D. Power spectrum analysis of the spontaneous fluctuations in the photoplethysmographic signal. J Basic Clin Physiol Pharmacol 1994; 5: 269–276. https://dx.doi.org/10.1515/jbcpp.1994.5.3-4.269.
- Kiselev AR, Gridnev VI, Prokhorov MD, Karavaev AS, Posnenkova OM, Ponomarenko VI, Bezruchko BP. Selection of optimal dose of beta-blocker treatment in myocardial infarction patients basing on changes in synchronization between 0.1 Hz oscillations in heart rate and peripheral microcirculation. J Cardiovasc Med (Hagerstown) 2012; 13: 491-498. https://dx.doi.org/10.2459/JCM.0b013e3283512199.
- Kiselev AR, Gridnev VI, Prokhorov MD, Karavaev AS, Posnenkova OM, Ponomarenko VI, et al. Evaluation of five-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: 204-213. https://dx.doi.org/10.1111/j.1542-474X.2012.00514.x.
- Kiselev AR, Gridnev VI, Prokhorov MD, Karavaev AS, Posnenkova OM, Ponomarenko VI, Bezruchko BP. Effects of antihypertensive treatment on cardiovascular autonomic control: a prospective study. Anadolu Kardiyol Derg 2014; 14: 701-710. https://dx.doi.org/10.5152/akd.2014.5107.
- Kiselev AR, Shvartz VA, Karavaev AS, Mironov SA, Ponomarenko VI, Gridnev VI, Prokhorov MD. Correlations between cardiovascular autonomic control indices during the two-hour immobilization test in healthy subjects. Open Cardiovasc Med J 2016; 10: 35-43. https://dx.doi.org/10.2174/1874192401610010035.
- Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology. Heart rate variability: standards of measurement, physiological interpretation, and clinical use. Circulation 1996; 93: 1043–1065. https://dx.doi.org/10.1161/01.CIR.93.5.1043.
- Pikovsky A, Rosenblum M, Kurths J. Synchronization: а universal concept in nonlinear sciences. Cambridge, England: Cambridge University Press, 2001; 432 p.
- Tass P, Rosenblum MG, Weule J, Kurths J, Pikovsky A, Volkmann J, et al. Detection of n:m phase locking from noisy data: Application to magnetoencephalography. Phys Rev Lett 1998; 81: 3291–3294. https://dx.doi.org/10.1103/PhysRevLett.81.3291.
- Haaksma J, Dijk WA, Brouwer J, van den Berg MP, Dassent WRM, Mulder B, Crijns HJGM. The influence of recording length on time and frequency domain analysis of heart rate variability. Comput Cardiol 1998; 25: 377–380. https://dx.doi.org/10.1109/CIC.1998.731818.
- McNames J, Aboy M. Reliability and accuracy of heart rate variability metrics versus ECG segment duration. Med Biol Eng Comput 2006; 44: 747-756. https://dx.doi.org/10.1007/s11517-006-0097-2.
- Munoz ML, van Roon A, Riese H, Thio C, Oostenbroek E, Westrik I, et al. Validity of (ultra-)short recordings for heart rate variability measurements. PLoS ONE 2015; 10: e0138921. https://dx.doi.org/10.1371/journal.pone.0138921.
Received 10 July 2018, Revised 14 September 2018, Accepted 26 September 2018
© 2018, Karavaev A.S., Skazkina V.V., Borovkova E.I., Kiselev A.R., Ponomarenko V.I., Kulminskiy D.D., Gridnev V.I., Prokhorov M.D., Bezruchko B.P.
© 2018, Russian Open Medical Journal
Correspondence to Anton R. Kiselev. Address: Institute of Cardiological Research, Saratov State Medical University, Bolshaya Kazachya str., 112, Saratov, 410012, Russia. E-mail: firstname.lastname@example.org.