## Abstract

Pre-eclampsia (PE), a serious pregnancy-specific disorder, causes significant neonatal and maternal morbidity and mortality. Recent studies showed that cardiovascular variability parameters as well as the baroreflex sensitivity remarkably improve its early diagnosis. For a better understanding of the dynamical changes caused by PE, in this study the coupling between respiration, systolic and diastolic blood pressure, and heart rate is investigated. Thirteen datasets of healthy pregnant women and 10 of subjects suffering from PE are included. Nonlinear additive autoregressive models with external input are used for a model-based coupling analysis following the idea of Granger causality. To overcome the problem of misdetections of standard methods in systems with a dominant driver, a heuristic ensemble approach is used here. A coupling is assumed to be real when existent in more than 80 per cent of the ensemble members, and otherwise denoted as artefacts. As the main result, we found that the coupling structure between heart rate, systolic blood pressure, diastolic blood pressure and respiration for healthy subjects and PE patients is the same and reliable. As a pathological mechanism, however, a significant increased respiratory influence on the diastolic blood pressure could be found for PE patients (*p*=0.003). Moreover, the nonlinear form of the respiratory influence on the heart rate is significantly different between the two groups (*p*=0.002). Interestingly, the influence of systolic blood pressure on the heart rate is not selected, which indicates that the baroreflex sensitivity estimation strongly demands the consideration of causal relationships between heart rate, blood pressure and respiration. Finally, our results point to a potential role of respiration for understanding the pathogenesis of PE.

## 1. Introduction

In recent years, the establishment of a model-based analysis for clinical use has been one of the greatest challenges in medical physics (Fenner *et al.* 2008; Gavaghan *et al.* 2008). Therefore, understanding the complex interactions in human physiology is crucial. But the complexity of this system in terms of structural and functional relationships complicates this task. Especially, the detection of changes caused by dynamical diseases is difficult, for instance in pre-eclampsia (PE), which is a serious disorder of the cardiovascular system during pregnancy. It is characterized by hypertension (mean systolic/diastolic blood pressure greater than 140/90 mmHg) and proteinuria of more than 300 mg in 24 h. The manifestation of PE is the main cause of maternal and neonatal morbidity and mortality, and about 2–5% of all pregnancies are affected.

In a recent study (Malberg *et al.* 2007), we have shown that including heart rate, systolic and diastolic blood pressure, and baroreflex sensitivity (strength of the influence of systolic blood pressure changes on the heart rate) in the analysis improves the prediction of the standard diagnostic tool (Doppler sonography) remarkably from 30 to above 70 per cent. The risk of PE was indicated by a combination of uterine perfusion (Chien *et al.* 2000), an increased normalized power of the very low frequency of the heart rate (Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology 1996), an increased power of the high frequency of the diastolic blood pressure (Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology 1996) and an increased number of baroreflex ramps in the range of 4–6 ms mmHg^{−1} (Malberg *et al.* 2007). The selected parameters of diastolic blood pressure variability and baroreflex characterize the short-term behaviour of the cardiovascular system, and are mainly determined by respiratory influences on the blood pressure and heart rate. In this study, the coupling of respiration, blood pressure and heart rate is investigated in the short-term range in order to understand these changes.

There are many approaches to quantify the direction and the strength of these interactions (Porta *et al.* 2002; Rosenblum *et al.* 2002; Nollo *et al.* 2005; Verdes 2005; Marwan *et al.* 2007; Palus & Vejmelka 2007; Faes *et al.* 2008). In most of these studies, either the relationship between the heart rate and the systolic blood pressure or that between the heart rate and the respiration were investigated. However, a consideration of more than one possible coupling, as given here, is rare. Various typical nonlinear phenomena, e.g. synchronization or amplitude–frequency coupling in the cardiovascular system, have been observed and therefore must be regarded in the modelling process. Cohen & Taylor (2002) proposed a nonlinear transformation of the inputs in order to describe known nonlinear phenomena of heart rate, e.g. saturations and threshold effects or hysteresis. Because such transformations are usually unknown, in this study a coupling analysis based on a non-parametric model is used, which follows the idea of Granger causality (Granger 1969; Wessel *et al.* 2007; Riedl *et al.* 2009).

The studied groups and their measurements are described in §2. In §3, the model-based approach that is used for the coupling analysis is described. Our results are summarized in §4 and interpreted in §5. Finally, §6 includes the conclusions of the results and their meaning for the diagnosis of PE.

## 2. Data and preprocessing

Multiple longitudinal measurements of seven healthy pregnant women (13 datasets) and four suffering from PE (10 datasets) in the course of pregnancy are considered (Malberg *et al.* 2007). All women underwent Doppler sonography at the Department of Obstetrics and Gynaecology of the University of Leipzig. Immediately after the Doppler examination, continuous blood pressure was measured non-invasively via finger cuff (100 Hz, Portapres device model 2, BMI−TNO, Amsterdam, The Netherlands). Additionally, the respiration curve (via respiratory effort sensors at the chest; sampling rate 10 Hz) was recorded. The measurements were performed for subjects in a supine position with relaxed breathing at times between 8.00 and 12.00.

The study was approved by the local ethics committee and obtained the informed consent of all subjects.

Measurements with disturbed respiratory signals or pathological respiratory patterns, e.g. Cheyne–Stokes breathing, are not included in the analysis. From the electrocardiogram, the detection of the heart beats is done by an algorithm (Suhrbier *et al.* 2006). Intervals between successive heart beats (*B*_{i}, the beat-to-beat interval) are calculated and assigned to the latter heart beat. The maximum values of the blood pressure curve in each beat-to-beat interval are extracted as the time series of systolic blood pressure (*S*_{i}). The corresponding diastolic blood pressure value is estimated by the next following blood pressure minimum (*D*_{i}). The values of the respiration signal are determined at the times of the *B*_{i} (*R*_{i}, respiration on beat-to-beat basis). The main objective of the analysis of heart rate and blood pressure is to investigate the cardiovascular system during normal sinus rhythm. Therefore, it is necessary to exclude not only artefacts (e.g. double recognition, i.e. R-peak and T-wave recognized as two beats) but also beats not coming from the sinus node of the heart, so-called ventricular premature complexes that are not directly controlled by the autonomous nervous system. These features are excluded from the time series of the *B*_{i} by an adaptive filter algorithm (Matlab implementation is available from http://tocsy.agnld.uni-potsdam.de; Wessel *et al.* 2000). A representative example of extracted time series is shown in figure 1.

## 3. Methods

The aim of this study is the investigation of the coupling between various cardiovascular measurements (*B*_{i}, *S*_{i}, *D*_{i} and *R*_{i}) in healthy pregnant women and in subjects suffering from PE. Causality plays an important role in the explanation of coupling. Therefore, the principle of Granger causality (Granger 1969) is used: given two simultaneously recorded time series, then Granger causality quantifies the causal directional influence from one time series to another. For both time series, autoregressive models with and without the other time series as external input are fitted. If the additional external input leads to a significant reduction in the variance of the predicting error, then the external input is said to have a causal influence on the response variable.

In this study, we consider the time series of *B*_{i}, *S*_{i}, *D*_{i} and *R*_{i}. Because there are possibly nonlinear functional interactions (Cohen & Taylor 2002; Wessel *et al.* 2006; Riedl *et al.* 2008), the following nonlinear additive autoregressive (NAARX) models with external inputs are fitted to the time series of *B*_{i}, *S*_{i}, *D*_{i} and *R*_{i}:
3.1
3.2
3.3and
3.4
The first term of each equation includes the mean value of the considered time series because the non-parametric fit is unique except for a constant. The second and the third terms include the autoregressive part and the external inputs, respectively. The last part of these equations is white noise that describes all other external inputs that are not explicitly considered in the model; and *p* is the maximal order of the models considered. Because there is no information about the weight functions *l*_{i}, *f*_{i}, *g*_{i} and *k*_{i}, they are estimated by a non-parametric iterative least-squares routine (see appendix A) fitting the model. This will ensure the highest possible flexibility of the additive model. Therefore, not only are linear functions possible, but also polynomials, jump functions or transcendental ones, for instance. The weighted predictors are selected by a model selection approach called adaptive backfitting (see appendix A), in order to increase the predictability of the response variable. The models are fitted separately to the time series, i.e. we do not consider closed-loop models. During this algorithm, the different autoregressive predictors as well as external ones are compared with each other with regard to the best prediction. The improvement of the prediction is measured by a cross-validation criterion (see appendix A), which is weighted by the cost of adding a new predictor or increasing the roughness of the estimated curve. This criterion is minimized during the iterations. In order to avoid over-fitting, this loop is stopped when the value of the criterion does not decrease by about more than 10 per cent after the last iteration. Hastie & Tibshirani (1990) indicate that the comparison of two successive values of the cross-validation criterion is equivalent to an approximated *F*-test. Our simulations show that the threshold of 10 per cent corresponds to the significance level of *α*=0.01. The selection of an external input indicates its causal influence on the response variable in relation to Granger causality (Riedl *et al.* 2009). That is, the additional values lead to a significantly better prediction of the response variable than the considered autoregressive part only.

In the case of the considered cardiovascular measurements, the problem of a coupling analysis is that *B*_{i}, *S*_{i} and *D*_{i} contain respiratory-induced oscillations with equal frequency. Simulations show that the proposed method finds not only the real couplings, but also additional false detections that occur randomly. To overcome this problem, a heuristic approach is chosen. We assume that the structure of coupling between *B*_{i}, *S*_{i}, *D*_{i} and *R*_{i} is equal in each subject of a group. That is, datasets of the group are considered as independent runnings of one system (ensemble) that differ by a random effect, initial conditions and boundary conditions. The built ensembles enclose the first 300 data points of each considered dataset in order to minimize the influence of long-term effects, e.g. trends, on the estimation process. We expect that almost all real couplings can be detected, but there are also false detections. Such failures could be caused by corrupt signals and overfitting. In order to distinguish between existing and artificial couplings, the relative frequency of the selection of an external input (table 1) is divided by a chosen threshold. If this relative frequency is greater than or equal to 80 per cent, then it is assumed that the response variable is influenced by the considered external input (table 1). Based on these results, the system of directed interactions between the four measurements is reconstructed (see figure 2). In a second step, the strength and the morphology of the selected couplings are analysed. Therefore, the explained variance of the separate one-step prediction of each external input is calculated and normalized by dividing by the variance of the response variable (table 2). These resulting values are indicators for the strength of coupling between the several external inputs and the response. If it increases, then the strength of the coupling also rises. After this, the frequently selected delayed predictors in each selected external input are searched and their parametrized weight functions compared in the group and between them in order to find out group-specific similarities and differences between healthy and pathological coupling (figures 3 and 4).

In order to test the reliability of this approach, the described analysis is repeated for two additional non-overlapping parts of the datasets (data point intervals [301,600] and [601,900]).

## 4. Results

We obtained the following results analysing the coupling of the cardiovascular measurements in pregnant women. First, we found a respiratory influence on *B*_{i}, *S*_{i} and *D*_{i}, which was detected in the whole ensemble, i.e. in more than 80 per cent of the considered cases. There is also an influence of *D*_{i} on *S*_{i} as well as an influence of *B*_{i} on *D*_{i} (table 1). Additionally, the influence of *B*_{i} on *S*_{i} is detected at the beginning of the measurement (first part) in healthy pregnant women as well as patients suffering from PE. In this part, there is also a detected influence of *D*_{i} on *R*_{i} in the group of patients. Second, we find a significantly (*p*<0.01; Wilcoxon’s sign test) increased value of the relative explained variance in the case of the respiratory influence on *D*_{i} in patients in comparison with healthy subjects (cf. table 2). The values of the relative explained variance are quite similar in all three considered parts. The strongest effect is between *D*_{i} and *S*_{i}, where the diastolic values seem to drive *S*_{i} for both healthy persons and patients. It is followed by the respiratory influence on *B*_{i}. The third largest effect is the respiratory influence on *D*_{i}. Selected couplings of a smaller effect are the influence of *B*_{i} on *D*_{i} and the respiratory influence on *S*_{i} (table 2). It is remarkable that the relative explained variance of the selected couplings is mostly larger than that in the opposite direction. With this finding, the structure of the couplings is reconstructed (figure 2). In the first part, there is the same structure with the mentioned additional selection. Finally, the morphology of the weight function of the frequently selected external predictors is investigated. The weight function (transformation of the diastolic predictor of lag 0 in the NAARX model of *S*_{i}, equation (3.2)) shows an increased linear influence above the mean *D*_{i} in healthy subjects (figure 3*c*). The horizontal shift resulted from the different mean values of *D*_{i} in the considered datasets. In the case of patients, the weight function shows a similar slope for lower values as but a decreasing linear influence above the mean *D*_{i} (figure 3*d*). It is important to emphasize that the morphology of (equation (3.1)) is very different in the groups of healthy subjects and patients (figure 3*a*,*b*). In the case of healthy women, the influence is linear in the range of mean values of respiratory movement, but decreases at the boundaries. For patients, a plateau is found in the range of mean value of the predictor and an increased linear influence at the boundaries. The fit of a piecewise linear function on the estimations of in figure 3*a*,*b* (cut-offs −0.002 and 0.002) shows that there is a significantly different slope in the interval [−0.002, 0.002] (*p*<0.01; Wilcoxon’s sign test). For healthy persons, there is a mean slope of −6800±1000, whereas −1700±800 is the value in the case of PE patients. The weight function (equation (3.3)) has a similar form in both groups, with a linear slope in the range of the mean value of the predictor and lower slopes at the boundaries (figure 4). There is a slightly increased middle range in the case of patients.

## 5. Discussion

In the first step of this analysis, the structure of the couplings between *R*_{i}, blood pressure and *B*_{i} is reconstructed. A main point is to answer the question whether the respiratory fluctuation in *B*_{i} is caused directly by respiratory modulation of the neuronal activity (Eckberg 2003) or indirectly via respiratory-induced *S*_{i} oscillations that are transmitted to *B*_{i} by the baroreflex (deBoer *et al.* 1987). The reconstruction shows that there is no significant influence of *S*_{i} on *B*_{i} in the short-term range, but a respiratory influence on both measurements. This result supports the model of a respiratory modulation of the neuronal activity (Eckberg 2003), which leads to the found respiratory fluctuation of *B*_{i}. The selected respiratory influences on *D*_{i} could be explained by a change of the intrathoracic pressure filling the lungs, which leads to higher pressure on the aortic tree because *D*_{i} is mainly determined by the ‘Windkessel’ function of the aortic tree. *D*_{i} also depends on the cardiac stroke volume, which depends on *B*_{i} (Starling’s law). Instead, the selection of the respiratory influence on *S*_{i} indicates the dependence on the peripheral resistance, which is controlled by the respiratory-modulated neuronal activity. The strong influence of *D*_{i} on *S*_{i} (see figure 1) could be explained by the function of the ‘Windkessel’ system of the aortic tree, which is influenced by stroke volume and the compliance of the vascular system. The ‘Windkessel’ function means that about half of the cardiac stroke volume is buffered by stretching the vessels near the heart and release during the filling time in order to decrease the differences between *S*_{i} (ejection time of the heart) and *D*_{i} (filling time of the heart). This buffer is mainly determined by the elastic property of the aortic tree (compliance rather than stiffness). It is known that there is a decreased compliance caused by PE that results in increased *S*_{i} (Dart & Kingwell 2001).

The explained structure (figure 2) is found in healthy subjects and patients during the second and third parts of the measurement. The additional selected interaction from *B*_{i} to *S*_{i} during the first part seems to be the end of a transient phase, which is caused by a change of posture from sitting or standing to the supine position.

The morphology of the weight functions of the frequently selected external predictors has mostly a nonlinear shape that could be quantified by fitting piecewise linear functions (threshold autoregressive models with external input). The weight function (equation (3.1)) is substantially different in the groups of healthy subjects and patients (figure 3*a*,*b*). In the case of healthy women, there is a linear dependence in the range of mean values of respiratory movement and saturations at the boundaries. These saturations can be explained by the limit of the respiratory gating (Eckberg 2003). For patients, there is no influence in the range of the mean value of the predictor and a linear dependence at the boundaries. It seems that the setting points for gating is dilated to more extreme values of nerve activity, e.g. baroreflex activity caused by extreme mean blood pressure values. This indicates an adaptive neuronal change that is caused by larger blood pressure fluctuations during developing PE. The weight functions (equation (3.3)) and (equation (3.2)) show a linear dependence of the blood pressure value on *R*_{i} in the range of the mean value of the predictor. This respiratory influence of blood pressure can be explained by a direct respiratory modulation of the stroke volume by changing the pleural pressure that caused greater changes in *D*_{i} than in *S*_{i} (Guz *et al.* 1987). At the boundaries of this range, the curves are saturated, which indicates that the respiratory movement at the end of inspiration and expiration leads to no significant additional change in the pleural pressure. The significantly different effect sizes of *R*_{i} on *D*_{i} that are found are explained by a delayed saturation phase in patients.

The weight function (equation (3.3)) shows rather linear dependence of *D*_{i} on *B*_{i}, which describes the second main influence on the stroke volume.

With the analysis of three successive parts of the measurement, it is shown that the coupling analysis by means of an adaptive backfitted NAARX model leads to reliable results. The comparison of the most commonly selected predictors shows that there are group-specific shapes of the estimated weight functions. For values of the relative explained variance of several external inputs smaller than 0.1, this property strongly decreases, characterized by a changed ratio of random effects and deterministic behaviour. There are, however, some limitations to this approach. The most important one is the analysis of short-term couplings. Expected interaction, for instance, the baroreflex or the renin–angiotensin system, seems to work beyond the order of five.

## 6. Conclusions

The investigation of the couplings between *B*_{i}, *S*_{i}, *D*_{i} and *R*_{i} shows that the resulting reliable structure of interaction is equal in healthy pregnant women and patients suffering from PE (figure 2). For both groups, there is a significant respiratory influence on *B*_{i} and blood pressure, an influence of *B*_{i} on *D*_{i} and a dependence of *S*_{i} on *D*_{i} for pregnant women in a supine position with relaxed non-pathological breathing.

Although the structure of the interactions is equal, there are clear differences in effect size and morphology of several interactions, which indicates a neuronal change by the resetting of the respiratory gating to more extreme respiratory values in patients. As expected, there is an increased respiratory influence on *D*_{i} in patients, which corresponds with a higher power of the high frequencies in diastolic blood pressure variability (Malberg *et al.* 2007). Additionally, the nonlinear form of the respiratory influence on the heart rate is significantly (*p*=0.002) different between the two groups. This result supports the model of a respiratory modulation of the neuronal activity (Eckberg 2003), which leads to the found respiratory fluctuation of *B*_{i}.

Interestingly, no influence of *S*_{i} on *B*_{i} is found (figure 2). Therefore, the baroreflex events (Malberg *et al.* 2007) result from the simultaneous short-term influence of respiration on heart rate and systolic blood pressure rather than from influences of systolic blood pressure on heart rate. This finding supports the conclusions of Porta *et al.* (2000), where it is clearly shown that a physiologically meaningful estimation of the baroreflex sensitivity demands the consideration of the causal relationship between heart rate, systolic blood pressure and respiration.

In further studies, this approach should be expanded to larger scales in order to include the baroreflex as well as interaction resulting from pathological respiratory patterns that frequently appear in pregnant women. Another point for future studies is to add the sympathetic activity to the considered cardiovascular measurements because it probably contains information that is masked by the interaction of *S*_{i} and *D*_{i}. Finally, one needs to analyse whether the change of the respiratory influence is an adaptive change as a response of the development of PE or an important factor by itself in the pathogenesis of this disease.

## Acknowledgements

We are grateful for financial support by the Deutsche Forschungsgemeinschaft, grant nos. KU-837/20-2 and KU-837/23-1, as well as by the EU Network of Excellence, grant no. NoE 005137 Biosim, and the EU project BRACCIA.

## Appendix A

#### (a) Fitting additive models

The additive model is defined by
A1
where *C* is a constant value, *X*_{j} is the *j*th predictor having the effect *f*_{j}(*X*_{j}) on the influencing variable *Y* , and *ϵ* is a *N*(0,*σ*^{2})-distributed random variable that is stochastically independent of the predictors. If the assumptions for equation (A1) are correct, then, for any *k*,
A2

By backfitting, the transformations *f*_{k}(*X*_{k}) are estimated as follows:

Set

*C*to the mean value of*Y*and choose the starting values of the transformations .For each

*j*, a non-parametric regression of the partial residue is made over the values of the*j*th predictor (estimation of equation (A2)): A3 Here,*S*_{j}is a scatter plot smoother (like a running mean, running line regression or kernel estimators) for the*j*th predictor; (,*y**x*_{1},…,*x*_{p}) is a multivariate time series of the random variables*Y*,*X*_{1},…,*X*_{p}; and*f*_{k}are the transformed values of*k*th predictor.Continue with step 2 until each estimated transformation no longer changes.

Details can be found in Hastie & Tibshirani (1990, p. 82 in section ‘Fitting additive models’).

#### (b) Running line smoother

The running line smoother is defined by
A4
Here, *s*(*x*_{0}) is the new value of *Y* at the point *x*_{0}; and *a* and *b* are the coefficients of the linear regression over the neighbours of *x*_{0}. Because the estimation is a kind of linear procedure, it can be expressed by the vector product of a smoothing matrix ** S** and the values of the response variable. There is a trade-off between the bias and the variance of the smoothing that is regulated by the number of neighbours in the regression. The variance of the smoothing decreases and the bias tends to increase if the number of considered neighbouring points rises and vice versa. This value controls the smoothing and is called the smoothing parameter

*λ*. Asymptotic analysis of the nearest-neighbour smoother has shown that

*λ*=

*N*

^{4/5}(

*N*, sample size) is the optimal value of the smoothing parameter if the values of the predictors are equidistant. More about this smoother can be found in Hastie & Tibshirani (1990, p. 15 in section ‘Running-mean and running-line smoothers’).

#### (c) Adaptive backfitting

The adaptive backfitting uses the minimization of the approximated cross-validation criterion (general cross-validation (GCV)) in order to select the best model:
A5
Here *λ*_{i} are the smoothing parameters that are used in the backfit; *f*_{j} and *λ*_{j}(*x*_{ij}) are the estimate functions at the target point *x*_{ij}, which are dependent on the chosen smoothing parameter; *n* is the number of considered data points and *p* the number of predictors; and tr *S*_{j} is the trace of the smoothing matrices in the backfit. The criterion is constructed in order to estimate the mean prediction error of the considered additive model. It is based on the estimation of the transformations *f*_{i} at the target point *x*_{i} without the value at this point (Hastie & Tibshirani 1990). The comparison of the estimation and the real value at the target point mimic the training and test set of a cross-validation approach.

The minimization of GCV is carried out for one predictor at a time, by applying the corresponding partial residuals looking for the smoothing parameter to minimize this criterion. Having obtained the best smoothing for each of the *p* available predictors, only the minimizing smooth is updated. This is continued until the criterion converges. Details can be found in Hastie & Tibshirani (1990, p. 262).

## Footnotes

One contribution of 10 to a Theme Issue ‘Experiments in complex and excitable dynamical systems’.

- © 2010 The Royal Society