## Abstract

The asymmetry of coupling between complex systems can be studied by conditional probabilities of recurrence, which can be estimated by joint recurrence plots. This approach is applied for the first time on experimental data: time series of the human cardiorespiratory system in order to investigate the couplings between heart rate, mean arterial blood pressure and respiration. We find that the respiratory system couples towards the heart rate, and the heart rate towards the mean arterial blood pressure. However, our analysis could not detect a clear coupling direction between the mean arterial blood pressure and respiration.

## 1. Introduction

The cardiorespiratory system is complex with direct and indirect interactions in its sub-components. It includes not only mechanical components reflecting the changing pressure in the thoracic region, but also the autonomic nervous control of both systems as well as a control of diaphragm and external intercostal muscles by means of the somatic nervous system. Investigating and understanding the couplings can help to identify and characterize different physiological and even pathological states, important for diagnosing and assessment of diseases [1–3].

Different linear and nonlinear approaches have been applied for studying couplings within the cardiorespiratory system, such as spectral analysis [4–7], Granger causality [8,9], phase dynamics [10,11], conditional information [12,13], joint symbolic dynamics [14,15] and model-based linear closed-loop approaches [16,17]. The main findings are a dependence of the couplings from the body position where the interaction between respiration and heart rate is dominant during the supine position. In contrast to that, the connection between heart rate and systolic blood pressure dominates the upright position. In all these approaches, the considered data have been on a beat-to-beat basis. There are only a few works that use continuous signals [18,19]. The most important difference is the use of a continuous blood pressures with its pulsating characteristic. Systolic and diastolic pressures correspond only to events in this cyclic change. Therefore, an extension to a continuous signal leads to a difficult interpretation. For this reason, the mean blood pressure value is considered, which does not include the pulse pressure (systolic pressure–diastolic pressure) variability. This causes a reduced high-frequency component of the mean blood pressure that corresponds to the mechanical influence of respiration.

Recently, a new nonlinear method for studying coupling directions has been proposed, which is based on recurrence plots (RPs) [20–22]. An RP itself is a powerful concept allowing the investigation of a variety of aspects of complex systems, such as transition studies, dynamical regime characterization, synchronization analysis or surrogate constructions [23,24]. Bivariate extensions are cross RPs and joint RPs, which can be used to study complete and generalized synchronization, respectively [23]. A further RP-based approach to study couplings between two systems is based on the probability of recurrence and can be used to detect phase synchronization [23]. This latter approach can be extended to a conditioned version allowing inference of coupling directions [20]. Its potential has been demonstrated on prototypical model systems, but not yet on experimental data.

Here, we will apply the approach of conditioned mean recurrence probabilities for the first time on experimental data that come from a study on cardiovascular variability in pregnancy and their change during preeclampsia. In general, the heart rate and the mean blood pressure increase with pregnancy [25]. A change of the coupling between heart rate and systolic blood pressure may also be observed [26]. By definition, the blood pressure is significantly larger in preeclampsia than in normal pregnancy. But there is also a change in cardiovascular regulation that is indicated by decreasing respiratory sinus arrhythmia [25], as well as an increased respiratory influence on diastolic blood pressure and a higher number of baroreflex events, an influence of systolic blood pressure on heart rate [27]. The changes of the respiratory influence on heart rate and diastolic blood pressure have been confirmed by a model-based approach [28]. However, there was no indication of an interaction between heart rate and systolic blood pressure; only an indirect coupling from heart rate to systolic blood pressure via diastolic blood pressure was found.

## 2. Data

Blood pressure and respiration have been measured on 11 pregnant women multiple times (in total, 23 datasets) in the course of the second and third trimesters of pregnancy [27]. The continuous blood pressure was measured non-invasively via finger cuff (100 Hz, Portapres device model 2, BMI–TNO, Amsterdam, The Netherlands). The respiration curve *R* was recorded via respiratory effort sensors at the chest (sampling rate 10 Hz). The measurements were performed for subjects in a supine position with relaxed breathing at times between 08.00 h and 12.00 h. Measurements with disturbed respiratory signals or pathological respiratory patterns, e.g. Cheyne–Stokes breathing, have been excluded. Based on an algorithm by Suhrbier *et al.* [29], we have extracted the heart beats and calculated an average heart beat rate *H* (in Hz). The main objective of the analysis of heart rate and blood pressure is to investigate the cardiovascular system during normal sinus rhythm. Therefore, we have removed beats not coming from the sinus node of the heart, the so-called ventricular premature complexes that are not directly controlled by the autonomous nervous system. These features are exchanged for random values by an adaptive filter algorithm preserving the time relation (http://tocsy.pik-potsdam.de; [30]). The ratio of the frequency of the respiration and the heart signal does not necessarily have to be an integer number. Therefore, some important variability of the respiration signal is lost if resampled to the beat-to-beat-based time scale. In order to consider the entire variability of the respiration, we have to use the continuous signals. Systolic *S* and diastolic *D* blood pressures are estimated from the maxima and minima of the blood pressure curve. Instead, using systolic *S* and diastolic *D* blood pressures, we have calculated the mean brachial blood pressure and interpolated it to a ‘continuous’ signal because we are interested in temporal continuous values, but both *D* and *S* can only be used in an analysis on a beat-to-beat basis. All time series have been resampled to 10 Hz (figure 1). For the average heart beat rate *H* and the mean brachial blood pressure *B*, this was carried out by using a cubic interpolation.

## 3. Method

An RP is a representation of recurrent states of a dynamical system *X* in its *m*-dimensional phase space. A phase-space trajectory can be reconstructed from a time series by time delay embedding [31],
3.1where *m* is the embedding dimension, *τ* is the delay and *N*′=*N*−(*m*−1)*τ* is the number of phase-space vectors.

The embedding parameters *τ* and *m* can be estimated by using mutual information and false nearest neighbour method [32].

The RP is then a pair-wise test of all phase-space vectors *x*_{i} () among each other, whether or not they are close,
3.2with *Θ*(⋅) being the Heaviside function, ∥*x*_{i}−*x*_{j}∥ the spatial distance between the phase-space vectors and *ε* a threshold for proximity [23]. The indices *i* and *j* range in the interval [1,…,*N*′] and mark the time points along the phase-space trajectory of length *N*′. The binary recurrence matrix *R*^{X} contains the value one for all close pairs (figure 2*a*).

The average of the recurrence matrix is called the recurrence rate and corresponds to the mean probability that any state recurs. The probability that the system recurs to a certain state *x*_{j} can be estimated by the average of the corresponding column of the recurrence matrix, . For two coupled systems *X* and *Y* , we may ask for joint probabilities of recurrences in both systems. Such joint probabilities can be estimated from the joint RP,
3.3which represents simultaneous recurrences in systems *X* and *Y* . Analogously to the recurrence rate, averaging the matrix **JR**^{X,Y} delivers the joint recurrence rate, i.e. the probability *p*(*x*_{j},*y*_{j}) that we find a recurrence in system *X* and in system *Y* simultaneously. Thus, we can calculate the probability that the trajectory of *Y* recurs to the neighbourhood of *y*_{j} under the condition that the trajectory of *X* recurs to the neighbourhood of *x*_{j} by
3.4Its average is the mean conditional probability of recurrence, MCR [20,22],
3.5

In the presence of the asymmetry of the coupling (e.g. suppose *X* to be the driver and *Y* to be the response without loss of generality), we have the relationship
3.6The interpretation of this inequality is based on the difference of complexity between *X* and *Y* . If *X* drives *Y* , the dimension of *Y* will be larger than the dimension of *X* because the dynamics of *Y* is determined by both the states of *X* and *Y* , while *Y* does not influence *X*. Note that this only holds provided the coupling strength is smaller than the threshold for synchronization, as the coupling direction might be lost if both systems become completely synchronized. Increasing the coupling strength from *X* to *Y* increases the complexity of *Y* . This results in a decrease of the recurrence probability *p*(*y*_{j}) that *Y* recurs to the neighbourhood. However, the complexity of *X* remains constant with increasing coupling strength because *X* is independent of *Y* (not vice versa). Hence, the mean recurrence probability of 〈*p*(*x*_{j})〉>〈*p*(*y*_{j})〉, implying . Therefore, we have MCR(*Y* |*X*)<MCR(*X*|*Y*) if *X* is the driver (figure 2*b*).

Following Romano *et al.* [20], the criterion for selecting the threshold value *ε*_{X} and *ε*_{Y} is such that for coupling strength equal to zero, the recurrence rates (recurrence probabilities) in both systems should be equal. However, in a passive experiment where the coupling strength between both interacting systems cannot be adjusted systematically as for our measurement data, we cannot apply directly such criterion to choose *ε*_{X} and *ε*_{Y} because the value of the coupling strength is not known *per se*. Therefore, we apply another criterion in choosing the threshold [34]: we normalize the data beforehand to have zero mean and unit standard deviation, and then we choose *ε*_{X}=*ε*_{Y}=0.1. Consistent results are obtained for threshold values that are varied in the range [0.05, 0.2].

In this work, we are interested in testing the possible interaction direction between a pair of two time series. An investigation of indirect couplings between the three subsystems requires a systematic study involving all three subsystems simultaneously, which will be future work.

In the case of passive experiments, we often have one scalar measurement time series, as we have for the cardiovascular experiment. We need to statistically assess the significance of the so-calculated (often just one) direction value in order to decide whether the value is obtained by chance or whether it is significant. Therefore, we need an appropriate statistical test in order to test the null hypothesis that the two systems *X* and *Y* have an independent recurrence structure. To test such a null hypothesis, we use the phase randomization surrogate test [35]. Random phases are added to the Fourier transformed time series, which is then inversely transformed to derive the new time series (with different phases). When assessing the MCR(*X*|*Y*) value (where *X* and *Y* represent *R*,*H*, and *B* series, respectively) of one subject, we use the phase-randomized time series of the second time series *Y* as a surrogate *Y* ^{s} (we can also use the first time series *X* to create surrogates *X*^{s}, but it does not change the results). Repeating the phase randomization (in our work 100 times), we get an ensemble of many surrogate series *Y* ^{s} and, hence, a distribution of corresponding *MCR* values. The directionality indices MCR(*X*|*Y*) and MCR(*Y* |*X*) for one subject can now be compared with the distribution of MCR(*X*|*Y* ^{s}) and MCR(*Y* ^{s}|*X*), respectively. If *X* and *Y* are independent, the value MCR(*X*|*Y*) will not differ significantly from the distribution of the values MCR(*X*|*Y* ^{s}). Otherwise, i.e. when exceeding the 0.95 quantile, we can reject the null hypothesis, indicating that the obtained values for the directionality indices are significant with 95% confidence.

Summarizing, the following steps have to be undertaken to assess the coupling direction between two time series for each subject:

— choose the significance level

*α*=0.05;— compute MCR(

*X*|*Y*) and MCR(*Y*|*X*);— create 100 phase-randomized surrogates and compute MCR(

*X*|*Y*^{sj}) and MCR(*Y*^{sj}|*X*) for*j*=1,…,100;— calculate the

*α*quantiles of the distributions of MCR(*X*|*Y*^{s}) and MCR(*Y*^{s}|*X*);— if MCR(

*X*|*Y*) and MCR(*Y*|*X*) are larger than the corresponding*α*quantiles, reject the null hypothesis and consider them as significant; and— if MCR(

*X*|*Y*) and MCR(*Y*|*X*) are significant, we compare MCR(*X*|*Y*) and MCR(*Y*|*X*) regarding equation (3.6) in order to find the directionality of the coupling.

## 4. Results

Using mutual information and the method of false nearest neighbours, we have found optimal embedding parameters for *H* as well as for *R* to be *τ*=2 and *m*=3, which resulted from the average over all cases; for *B*, we have found *τ*=4 and *m*=2. The results of our analysis have not changed much when using different embedding parameters.

We have calculated the MCR measures for all combinations between respiration *R* and heart rate *H*, heart rate *H* and mean blood pressure *B*, and respiration *R* and mean blood pressure *B*. First, we check the significance of MCR in order to limit the subsequent directionality study to the significant results. An MCR value would be significant if it exceeds the 0.95 quantile of the surrogate MCR distribution. Based on this test, we find significant MCR indices between respiration *R* and heart rate *H* for all subjects, but between mean blood pressure and heart rate or respiration only for more than half of the subjects, although still a considerable number (figure 3 and table 1).

Next, we study the coupling direction between the significant couplings. According to equation (3.6), we compare which MCR value is larger. First, we check the coupling direction between respiration *R* and heart rate *H* (figure 4*a*). We find 21 significant cases where the MCR(*R*|*H*) value is clearly larger than MCR(*H*|*R*); thus, we can infer a coupling direction from respiration to heart rate. Preeclampsia and the progression of gestation have not caused the significance of MCR(*H*|*R*).

Then, we check the coupling between heart rate *H* and blood pressure *B* (figure 4*b*). Here, we find 15 significant cases where MCR(*H*|*B*) is larger than MCR(*B*|*H*), i.e. a coupling from heart rate *H* to blood pressure *B*. Including the non-significant MCR(*H*|*B*) values, we would have 18 cases with such a coupling direction. In five cases, we found an opposite coupling direction from blood pressure towards heart rate. However, the difference between the two MCR indices is, in more than 15 cases, small, indicating a potential bidirectional coupling.

Finally, a comparison between the significant MCR values of blood pressure *B* and respiration *R* reveals 13 cases with coupling directions from blood pressure to respiration and five cases from respiration to blood pressure (figure 4*d*).

The coupling directions between heart rate and blood pressure as well as between respiration and blood pressure are not as clear as between respiration and heart rate because the differences of the corresponding two MCR measures are small (figure 4*b*,*d*), and there are also some cases with opposite coupling directions (e.g. where MCR(*B*|*H*) is larger than MCR(*H*|*B*) in figure 4*b*). We might guess that this latter result could be due to preeclampsia. However, this is not the case, as the contradictory results appear for preeclampsia as well as for healthy women (for *H* versus *B* in two healthy and one preeclampsia, for *R* versus *B* in three healthy and one preeclampsia women).

Instead using the mean blood pressure *B*, we have also tested the upper envelope of the blood pressure series *S* (as an analogue for a continuous systolic blood pressure). This upper envelope can be interpreted as a representation of the current total peripheral resistance of the smaller blood vessels. Here, we found larger differences between MCR(*H*|*B*) and MCR(*B*|*H*), and finally 18 significant cases with a coupling from heart rate *H* to blood pressure *B* (figure 4*c*). This might be indicative for the mechanical coupling mechanisms affecting the blood pressure by the heart rate [19].

## 5. Conclusions

The investigation of coupling directions from experimental data is a challenging task [36]. The recently developed nonlinear method based on conditional recurrence probabilities [20] allows for a directionality analysis in coupled complex systems. Here, we have successfully applied this approach for coupling analysis in experimental data.

The application to data from the human cardiorespiratory system has clearly revealed a coupling from the respiratory system towards the heart. These findings support the assumption that the respiratory sinus-arrhythmia results from a direct influence of respiration on heart rate (respiratory gate [37]). It is assumed that the respiratory control centres modulate the vagal outflow in the brainstem.

Cardiovascular couplings are not as clearly detected as the respiratory coupling, which suggests that the respiratory-induced oscillation is the carrier of the couplings detected by the beat-to-beat approaches [6,11,13,14,26,28]. The proposed method has been able to detect that heart rate affects blood pressure (through mechanical coupling mechanisms [7,19]); but for some cases, also that arterial pressure affects heart rate (through the baroflex circuit). The small difference between the MCR measures also supports the potential bidirectional nature of the coupling between heart rate and blood pressure. An even less clear result was found for the coupling between respiration and blood pressure. This might be a hint to indirect coupling mechanisms. Moreover, here we have used continuous cardiorespiratory signals instead of beat-to-beat-based signals, which was the base in previous studies. In contrast to beat-to-beat signals, in the blood pressure series, the high-frequency variation is suppressed. These distinctions and also the fact that we extracted the heart rate from the blood pressure measurements might cause the differences in the coupling structure. Nevertheless, the proposed method could lead to additional information about the cardiorespiratory coupling in comparison with the beat-to-beat approach.

The particular database used in our study might also have some impact on our findings. Nevertheless, during our analysis we did not find any evidence that either preeclampsia or the progression of gestation had a significant impact on the results. However, a detailed analysis of the specific effect of pregnancy and preeclampsia on the cardiorespiratory coupling is out of the scope of this paper, and is the subject of future studies. Moreover, a thorough study about the accuracy of the detection of interaction directions (e.g. how much should MCR(*X*|*Y*) differ from MCR(*Y* |*X*)) is, in general, an open problem and also remains a subject for future work.

## Acknowledgements

This study was supported by the Potsdam Research Cluster for Georisk Analysis, Environmental Change and Sustainability (PROGRESS, support code 03IS2191B), the DFG research group ‘Himalaya: Modern and Past Climates (HIMPAC)’, a Hong Kong Polytechnic University Postdoctoral Fellowship and the National Natural Science Foundation of China (grant no. 11135001).

## Footnotes

One contribution of 13 to a Theme Issue ‘Assessing causality in brain dynamics and cardiovascular control’.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.