## Abstract

The investigation of foetal reaction to internal and external conditions and stimuli is an important tool in the characterization of the developing neural integration of the foetus. An interesting example of this is the study of the interrelationship between the foetal and the maternal heart rate. Recent studies have shown a certain likelihood of occasional heart rate synchronization between mother and foetus. In the case of respiratory-induced heart rate changes, the comparison with maternal surrogates suggests that the evidence for detected synchronization is largely statistical and does not result from physiological interaction. Rather, they simply reflect a stochastic, temporary stability of two independent oscillators with time-variant frequencies. We reanalysed three datasets from that study for a more local consideration. Epochs of assumed synchronization associated with short-term regulation of the foetal heart rate were selected and compared with synchronization resulting from white noise instead of the foetal signal. Using data-driven modelling analysis, it was possible to identify the consistent influence of the heartbeat duration of maternal beats preceding the foetal beats during epochs of synchronization. These maternal beats occurred approximately one maternal respiratory cycle prior to the affected foetal beat. A similar effect could not be found in the epochs without synchronization. Simulations based on the fitted models led to a higher likelihood of synchronization in the data segments with assumed foetal–maternal interaction than in the segment without such assumed interaction. We conclude that the data-driven model-based analysis can be a useful tool for the identification of synchronization.

## 1. Introduction

Neural integration is an essential part of foetal development. Changes in heart rate variability in the foetus during gestation reflect increasingly complex reactions to various internal and external conditions and stimuli (Lange *et al*. 2005; Van Leeuwen *et al*. 2007). Some of these stimuli may be associated with maternal state. For example, changes in foetal heart rate and heart rate variability have been associated with altered maternal arterial oxygen content (Bekedam *et al*. 1991; Jensen 1996), maternal hypothermia (Hankins *et al*. 1997) or maternal exercise (Webb *et al*. 1994). Furthermore, direct associations between foetal and maternal heart rate have been observed within time frames over hours or days (Patrick *et al*. 1982; Lunshof *et al*. 1998). Therefore, it is conceivable that the neural integration of the foetus may be investigated non-invasively via the foetal heart rate response to specific maternal states. These may include postural changes, rest/stress or controlled respiration. The characteristics of the foetal response may be useful in the detection of prenatal disease or deficits.

We have explored the possibility of foetal–maternal interaction by examining the concurrence of their respective heart rates on a beat-to-beat basis (Van Leeuwen *et al*. 2003). We were able to show that there is a certain likelihood of occasional heart rate synchronization between the mother and the foetus. Epochs of entrained heartbeats showed a tendency to specific phase preferences, which could not be found in similarly analysed surrogate data. We have postulated that the mechanism for this entrainment may originate in the foetal perception of and coupling to maternal rhythms via mechanical or auditory stimuli associated with these rhythms.

We have subsequently attempted to create conditions conducive to this kind of foetal–maternal interaction by performing data acquisition in a setting in which maternal respiration was controlled (Van Leeuwen *et al*. 2008). The results suggested that foetal–maternal heart rate entrainment may be induced by high maternal respiratory rates. One drawback of the data analysis used in these studies is that the evidence is largely statistical. In the comparison of the results from the surrogate and the original (true) data, there was no way of identifying a specific epoch in the original data as having resulted from physiological interaction, instead of simply reflecting a stochastic, temporary stability of two independent oscillators with time-variant frequencies. The aim of this work is to determine direct influences of maternal heartbeat durations on the foetal heartbeat train based on data-driven modelling analysis.

## 2. Method

### (a) Recording and processing foetal and maternal cardiac signals

Foetal and maternal cardiac electrophysiological activity was recorded based on magnetocardiography. The method permits the registration of the magnetic fields produced by the electrical activation of the heart and delivers a relatively noise-free signal, similar to the electrocardiogram (ECG), with high temporal resolution. It is particularly useful in the prenatal setting where the ECG has limited applicability (Koch 2001). Large-array biomagnetometers permit the simultaneous registration of maternal and foetal magnetocardiographic signals.

The data were collected with a 61-channel biomagnetometer (Magnes 1300C, 4D Neuroimaging, San Diego) in which the sensing channels are arranged concentrically on a slightly curved surface with 825 cm^{2} area of coverage. The coil diameter is 28 mm with a spacing of approximately 38 mm. The coil distance to the outer surface is 16 mm or less. Environmental noise cancellation is possible based on ambient noise detectors. Intrinsic system noise is specified as averaged over frequency and across all channels. The data were recorded for 5 min at a sampling rate of 1 kHz and a band pass of 1–200 Hz. The measurements were performed in a shielded room (AK3b, Vakuumschmelze, Hanau, Germany) in order to reduce the influence of electromagnetic artefacts.

The three datasets presented in this study originate from foetus–mother pairs studied at the end of pregnancy in whom we obtained six consecutive 5 min simultaneous foetal and maternal magnetocardiograms (MCGs; Van Leeuwen *et al*. 2008). After an adaptation period of 10–15 min, a baseline MCG was recorded with the mother breathing spontaneously. This was followed by four data acquisitions in which the mother was asked to adapt her respiratory rate by coordinating her inspiration to an auditory cue. The repetition rates of the cues for the four acquisitions were 15, 10, 20 and 12 cycles per minute (cpm), in that order. This was followed by a final sixth acquisition in which the mother again breathed spontaneously. The single acquisitions were separated by a 2–3 min pause without regulated breathing. In each of the 5 min MCG datasets, foetal and maternal R-peaks were identified to an accuracy of 1 ms using QRS templates. The corresponding RR interval time series (B: beat-to-beat interval) were constructed on the basis of the timing of the R-peaks.

Previous analysis (Van Leeuwen *et al*. 2008) has shown that, in particular, at a fast maternal respiratory rate of 20 cpm, the probability of occurrence of epochs of synchronization was the highest. For the present work, we considered only those data acquired at this rate (figure 1).

### (b) Determination of synchronization between foetal and maternal heart rates

In order to find periods of short-term interactions between foetal and maternal heart rate, we constructed the synchrograms using the stroboscopic technique (Schäfer *et al*. 1998). Accordingly, the phase of the maternal heartbeats was set to multiples of two times *π*. The momentary maternal phase of the foetal beat was linearly interpolated by calculation of *ϕ*(*t*_{i})=2*π*(*t*_{i}−*t*_{k})/(*t*_{k+1}−*t*_{k})+2*πk* with *t*_{k}≤*t*_{i}<*t*_{k+1}, where *t*_{k} is the start of the *k*th maternal cycle characterized by the *k*th maternal heartbeat and *t*_{i} is the time of the *i*th foetal heartbeat. To obtain the synchrogram, the relative cyclic phase *ψ*_{m}(*t*_{i})=(*ϕ*(*t*_{i})mod 2*πm*)/2*π* is plotted versus the times *t*_{i} with respect to *m* maternal cycles (figure 2). Horizontal parallel lines in this presentation indicate phase synchronization between the rhythmic maternal and foetal heartbeat. The number of these parallel line *n* embedded in *m* maternal *B*'s indicates the ratio of synchronization *n* : *m*. It must be pointed out that the value of *n* depends on the choice of *m* (Porta *et al*. 2004). The synchronization epochs were determined according to the so-called Toledo procedure (Toledo *et al*. 1998). For fixed combinations of *n* and *m*, the absolute difference *δ*_{i}=|*ψ*_{m,i}−*ψ*_{m,i+n}| was calculated and a sliding window of length *l*=2*m* (Cysarz *et al*. 2004) was moved over all *δ*_{i}. If all *δ*_{i} within the window are smaller than a threshold value *α*=0.025 (Cysarz *et al*. 2004), then the window is considered to be an epoch of synchronization with ratio *n* : *m*. We assumed that the synchronization caused by the short-term regulation of foetal heart rate depends on the maternal breathing cycle, which is reflected by the maternal sinus arrhythmia. The mean length of the maternal breath is 3 s and the mean maternal *B*'s are between 0.6 and 0.8 s. Therefore, we assumed about four maternal beats during the maternal breath and limited *m* to four. Because *B*_{maternal} is not greater than two times *B*_{foetal}, *n* is limited to *n*=2*m*, so we considered the combinations [*n*, *m*]=[2, 1], [3, 2], [5, 2], [4, 3], [5, 3], [5, 4] and [7, 4]. Figure 2 shows the corresponding synchrograms and detected epochs for dataset 1.

The basic assumption of this work is that, under appropriate conditions, the foetal heart rate may be influenced in the short-term by the maternal heart rate. Thus, we looked for synchronization associated with short-term regulation of the foetal heart rate. It was assumed that such synchronization is only explainable by specific short-term patterns in the foetal heart rate dependent on given maternal patterns. In order to test this, we replaced the short-term fluctuations in *B*_{foetal} by noise. The mean and the standard deviation of the normal distributed stochastically independent random numbers of this process were estimated from the whole foetal *B*'s. If there were substantially different foetal heart rate levels in one dataset then the mean and standard deviation of the noise are adapted to include this dynamic. A surrogate of the foetal beat-to-beat time series in figure 1 is shown in figure 3.

The synchronization analysis performed on the original maternal and foetal *B*'s (demonstrated in figure 2) was repeated for the combinations of original maternal and surrogate foetal *B*'s. As the surrogate time series destroyed short-term foetal heart rate patterns, synchronization resulting from foetal short-term regulation should be eliminated. If the analysis identifies synchronization epochs here, these will result from maternal (original) and foetal (surrogate) *B*'s that are temporarily but stochastically stable at a particular *n* : *m* ratio (Van Leeuwen *et al*. 2003). Synchronization epochs found in the original data during those periods in the surrogate analysis which did not display synchronization were considered to be epochs resulting from short-term foetal–maternal interaction. Thus, it is possible, based on the synchronization analysis of the surrogate data, to identify those epochs in the original data which do not result from such chance ratios. Therefore, we compared the synchrograms before and after replacing the foetal signal. An example of such comparison based on dataset 1 is given in figure 4. It shows the synchrograms using the surrogate data for those maternal cycles for which synchronization epochs could be found in the original data (*m*=2, *m*=3; cf. figure 2*b*,*c*). The epochs of synchronization found in figure 2 are framed in figure 4, showing that most of the 3 : 2 synchronization (figure 2*b*) is reproduced by the surrogate data. On the other hand, the 5 : 3 synchronization (figure 2*c*) is only partly explained. We therefore assumed that the 5 : 3 synchronization between the 100th and 150th second could indeed result from short-term foetal heart rate fluctuation. For further analysis, we selected those episodes of the original time series with synchronization which do not display synchronization for the foetal surrogate data. These episodes are declared as periods of synchronization caused by short-term regulation in the foetal cardiovascular system.

### (c) Data-driven modelling of synchronization

This work deals with the description of the short-term mechanism of foetal heart rate regulation, which leads to a synchronization between the foetal and maternal *B*'s. The aim of this investigation was to look for evidence of the existence of coupling between both heartbeat trains permitting a deeper insight into the character of this interaction.

Synchronization requires two self-sustained linear or nonlinear oscillators, which are weakly coupled. In this context, a nonlinear oscillator is the excitation of the foetal and maternal heart, respectively. The nonlinearity of these cardiac systems is well-known (e.g. Babloyantz & Destexhe 1988; Shono *et al*. 1991). The synchronization of both systems has been described but the character of the coupling between them is yet unknown (Van Leeuwen *et al*. 2003). Recently, we have attempted to create conditions conducive to synchronization using a fixed rate of maternal breathing.

In order to look for coupling between foetal and maternal heart rate, we used a model-based approach that follows the idea of the Granger (1969) causality. Granger causality quantifies the causal directional influence from one time series to another and is based on linear predictions of time series by autoregressive models. Consider two simultaneously recorded time series. An autoregressive model as well as a second autoregressive model with the first as external input is fitted to the single time series. 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.

Because the respiratory movement of the mother modulates her *B*'s (respiratory sinus arrhythmia), we included this signal in the description of the changes in foetal *B*'s. Additionally, we try to describe the maternal heart rate variability the other way around for comparison. We used a nonlinear additive autoregressive process with external influences (equations (2.1) and (2.2)), a beat-to-beat model which considers linear as well as nonlinear interaction because the form of the possible coupling is unknown,(2.1)and(2.2)The external influences are predecessors of the linearly interpolated maternal values at the time of the foetal heartbeats ( in equation (2.1)) and the foetal values in equation (2.2). The autoregressive part contains previous values of the foetal signal (equation (2.1)) and (equation (2.2)). All other unexplainable disturbances in the foetal as well as maternal beat-to-beat signal were summarized by a stochastic term and , respectively. Since we had no information about the transformations of the predictors *f*_{j}, *g*_{j}, *k*_{j} and *l*_{j}, we used a non-parametric iterative regression (see appendix A*a* and A*b*) because the form of the functions is unknown. In order to take the idea of Granger causality into account, this fitting routine is combined with a model selection algorithm (see appendix A*c*). This procedure selects stepwise single predictors from the *p* (in our case, 10 foetal beats), first predecessors of the foetal and maternal *B*'s, which predict mostly the response variable *B*_{foetal}. At the same time, the required linear or nonlinear transformations *f*_{j}, *g*_{j}, *k*_{j} and *l*_{j} of these predictors are estimated by an enclosed non-parametric regression called backfit (figure 5). The selection is stopped if the prediction does not increase significantly after adding a new predictor (measured by the GCV; see appendix A*c*). The mean value of the response variable 〈*B*_{foetal}〉 and 〈*B*_{maternal}〉 must be included because the backfit is unique except for a constant. The automatic selection of external predictors would indicate the existence of the coupling between foetal and maternal *B*'s. Their lags and transformations would give us information about the time scale and the kind of coupling, respectively.

In order to test the relevance of this approach, we looked at episodes of the time series both with and without assumed synchronization caused by a short-term regulation (i.e. original and surrogate) and analysed both in the same way. We expected that in the surrogate cases no external predictors will be selected for describing the short-term fluctuation of the foetal *B*'s.

We also tested our approach on the basis of the simulation of the foetal *B*'s by the fitted models. We used the original interpolated maternal signal for simulation. The standard deviation of the stochastic term was estimated by the deviation of the difference between the response variable and its one-step prediction of the model. The transformations *f*_{i} of the new foetal values were linearly interpolated. If the model is adequate, this should result in epochs of synchronization with the same ratio as in the original datasets and no synchronization with respect to the selected parts without synchronization (figure 6).

## 3. Results

In all three cases, synchronization was found in the original data. For the first example (figure 1), synchronization of ratios 3 : 2 and 5 : 3 could be detected (figure 2). The second example showed the combinations 2 : 1, 3 : 2, 5 : 3, 4 : 3 and 7 : 4 and in the third dataset only 5 : 3 and 7 : 4 synchronization could be found (synchrograms not shown). In all three cases, the majority of the detected synchronization events had the ratio 5 : 3.

The synchronization analysis of the foetal surrogates and the maternal signals showed that the largest part of the detected epochs resulted from stable runs of specific mean foetal and maternal *B*'s and could be reproduced by white noise. Only a few episodes with a maximum length of 10 s could not be explained by the surrogate results (figure 4). These were assumed to represent synchronization resulting from foetal short-term regulation. In order to produce reliable results from the model fitting, we chose a larger part of the dataset (length of 50 s), which includes these short epochs as well as periods without synchronization. In the first case (figure 1), we chose the interval from 100 to 150 s (figure 2). By contrast, the interval of 1–50 s (in figure 1) was presumed to be periods without synchronization associated with short-term regulation. The considered part is framed in figure 1, for example.

Fitting the proposed model to the foetal *B*'s in the selected periods of the original data, we found that there are external maternal predecessors (lag 5) also selected for the section where synchronization caused by short-term regulation is included.

The selected predictors and their estimated weight function of the framed section in figure 1*b* can be seen in figure 5. On the other hand, for the segment without synchronization coming from foetal short-term regulation only foetal predecessors are included in the model (equation (2.1)).

In the case of the second and third example external predictors, respectively, lags 1, 4 and 5 were selected describing the foetal *B*'s in the periods, including epochs of synchronization. On the other hand, there were no maternal terms in the model of foetal *B*'s without the appropriate synchronization for the third case but the ninth maternal predecessor was selected in the second example. On the other hand, there is no selection of foetal predictors in the model of maternal heart rate (equation (2.2)) except that, for the non-synchronized segment of the third case, the seventh foetal predictor was chosen.

To test the validity of the fitted models, the foetal *B*'s were simulated (figure 6). For the first example (figure 1), the model of foetal *B*'s (figure 5) was able to reproduce the synchronization (figure 6*a*(i)). By contrast, the model of foetal *B*'s without synchronization generates only a small epoch with synchronization (figure 6*b*(i)). A similar result was found in the second example (figure 6(ii)). For the third case (figure 6(iii)), however, there was also synchronization for the modelled epoch without synchronization, although the model did not include any external predictor. This synchronization corresponds to a very small variation in the simulated foetal signal and the natural 5 : 3 relation of foetal–maternal heart rates.

## 4. Discussion

It is well known that, in the prenatal condition, the foetus is sensitive to and can react to its environment. Many of the responses to stimulation can involve changes in cardiac function and heart rate. In particular, heart rate is accessible to measurement and remains one of the central descriptors of foetal physiology. The interaction between the temporal aspects of the maternal and foetal cardiac systems has been noted on various time scales (Patrick *et al*. 1982). Although a short-term association between the foetal and maternal heart rates has been described (Hildebrandt & Klein 1979), it has been a challenge to separate temporary chance, i.e. non-physiological, occurrence of numerical ratios of foetal and maternal cardiac oscillating frequencies from real physiological interactions (Van Leeuwen *et al*. 2003). Indeed, the comparison of the surrogate and original data presented in this paper (figure 4) points to the common occurrence of foetal–maternal heart rate ratios, which cannot be unequivocally attributed to physiological interaction. Thus, caution is required in the interpretation of the occurrence of epochs of synchronization between these two systems. One possible approach is to identify all dubious synchronization epochs based on the surrogate approach applied here and then to analyse and compare data segments including plausible cardiac coupling between mother and foetus as well as segments void of such synchronization.

The model-based analysis by means of the nonlinear additive autoregressive process with external influence has been shown to be a suitable tool for describing complex coupling in human heart rate regulation (Wessel *et al*. 2006, 2007; Penzel *et al*. 2007; Riedl *et al*. 2008). Because it is a non-parametrical approach, it includes different classes of nonlinear discrete models. One of them is the threshold autoregressive model that changes its behaviour after passing a specific state. Therewith, it has been possible to produce synchronization if an oscillating external influence was included (Tong 1990). In these datasets, the rhythmic external force is the maternal *B*'s, i.e. the respiratory modulation of this signal. In cases of synchronization, a maternal predecessor with a lag of approximately 4 or 5 was consistently chosen in order to describe the fluctuation in the foetal *B*'s. This suggests that there is relevant information on these changes within the maternal signal which is not included in the foetal signal itself. The selected lags indicate the important role of the maternal respiration represented by sinus arrhythmia in the generation of synchronization because the time span of five foetal heartbeats corresponds to the length of the maternal respiratory cycle of 3 s. It seems that the regulation of the foetal *B*'s is implemented by adapting to changes during the previous maternal respiratory cycle. To investigate the synchronization of both foetal and maternal heart rate with the maternal respiration, an extended experimental protocol is necessary, which would include the measurement of maternal breathing. If the respiration is additionally considered as an external input in the model of the foetal heart rate, it could be possible to decide whether there is a direct respiratory influence or not. However, examining the transformed external predictor in figure 5, we see that there is no continuous adaptation. *g*_{5} is only significantly in the range around . The parametrization of such discontinuity could be managed by a piece-wise linearization. Based on a threshold model, this would imply a constant maternal influence of the fifth maternal predecessor in the range of approximately 0.71 s and, beyond that, this predecessor would have no influence. This example indicates a nonlinear interaction, but, for more reliable results, larger epochs, more than 120 data points, have to be found for an analysis. In the third case, external predictors were also chosen in the model of foetal and maternal heart rate (lag 9 in equation (2.1) and lag 7 in equation (2.2), respectively) for the non-synchronized segment of the third case. However, the interaction in both directions indicates rather a mutual relationship. Simulations of the fitted model (figure 5) showed that the detected synchronization caused by short-term regulation can be reproduced (figure 6*a*(i)). On the other hand, simulations of the foetal *B*'s without synchronization related to short-term regulation also show synchronization (figure 6*b*(i)). In this case, however, it corresponds to a very small deviation in the simulated foetal *B*'s, and its mean value can be regarded as being responsible for this synchronization. Similar results were found in the second example (figure 6(ii)). Here, the ninth maternal predecessor was selected in the model for the part without short-term synchronization. This results from longer term similarities in the foetal and maternal signals, i.e. the shorter foetal *B*'s resulting from foetal movement coincide randomly with shorter maternal *B*'s. In the case of the third example (figure 6(iii)), the assumed synchronization associated with short-term regulation was also reproduced. Unfortunately, the selected part without assumed short-term synchronization is characterized by low variance. So, the fitted model is dominated by the mean value and the synchronization resulted from the above-mentioned statistical preferences between foetal and maternal mean *B*'s. It shows the limits of our approach. If the short-term fluctuations are very low, then we cannot distinguish between real foetal adaptation and noise-like behaviour. Another problem with our approach is its limited resolution. Although the analysis indicates synchronization repeated epochs with lengths of approximately 7 s (15–20 foetal heartbeats), we have to consider segments larger than 100 foetal heartbeats in order to ensure the reliability of the model fit. If an interrelationship between foetal and maternal heart rate is found we only assume that these short epochs appear significantly more often than in parts without supposed coupling. In further work, the coupling analysis has to be extended to more selective approaches. Nonetheless, our results support the presence of the presumed foetal–maternal interaction. This suggests that foetal short-term regulation of beat duration is influenced by maternal beat duration and represents evidence for true synchronization between both cardiac systems. Thus, data-driven model-based analysis can be a useful tool for the identification of interactions in order to quantify strength and frequency of foetal response on external stimuli dependent on the level of neural integration of the foetus. Such information may indicate neuronal prenatal diseases or deficits but this can only be shown in future studies with extended experimental protocols. Furthermore, the examination of the model predictors may help to identify the mechanisms involved in such synchronization.

## Acknowledgments

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) 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*_{i}(*X*_{j}) on the influencing variable *Y*; *ϵ* is a *N*(0, *σ*^{2}) distributed random variable that is stochastically independent of the predictors. If the assumptions for equation (A 1) 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 (A 2)),(A3)where*S*_{j}is a scatter plot smoother (such as a running-mean, running-line regression (see appendix A*b*) or kernel estimators) for the*j*th predictor. , which is a multivariate time series of the random variables*Y*,*X*_{1}, … and*X*_{p}. 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. 89, in the section ‘Fitting additive models’).

#### (b) Running-line smoother

The running-line smoother is defined by(A4)where *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, which is regulated by the number of neighbours in regression. The variance of the smoothing decreases and the bias tends to increase if the number of considered neighbouring points rise and vice versa. This value controls the smoothness 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 predictors are equidistant. More about this smoother can be found in Hastie & Tibshirani (1990, p. 15, in the 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); equation (A 5)) in order to select the best model,(A5)where *λ*_{i} are the smoothing parameters that are used in the backfit; *f*_{j}, *λ*_{j}(*x*_{ij}) are the estimate functions at the target point *x*_{ij} which depend 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 several smoothing matrices in the backfit. The denominator of GCV contains a cost function that takes the increased number of predictors and roughness (degrees of freedom of the fit) for a better prediction into account.

The minimization of GCV is carried out by one predictor at a time, by applying to 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, in the section ‘Adaptive backfitting: BRUTO’).

## Footnotes

↵† These authors contributed equally to the study.

One contribution of 15 to a Theme Issue ‘Addressing the complexity of cardiovascular regulation’.

- © 2009 The Royal Society