Evolution of cardiorespiratory interactions with age

We describe an analysis of cardiac and respiratory time series recorded from 189 subjects of both genders aged 16–90. By application of the synchrosqueezed wavelet transform, we extract the respiratory and cardiac frequencies and phases with better time resolution than is possible with the marked events procedure. By treating the heart and respiration as coupled oscillators, we then apply a method based on Bayesian inference to find the underlying coupling parameters and their time dependence, deriving from them measures such as synchronization, coupling directionality and the relative contributions of different mechanisms. We report a detailed analysis of the reconstructed cardiorespiratory coupling function, its time evolution and age dependence. We show that the direct and indirect respiratory modulations of the heart rate both decrease with age, and that the cardiorespiratory coupling becomes less stable and more time-variable.


Inverness IV2 3UJ, UK
We describe an analysis of cardiac and respiratory time series recorded from 189 subjects of both genders aged . By application of the synchrosqueezed wavelet transform, we extract the respiratory and cardiac frequencies and phases with better time resolution than is possible with the marked events procedure. By treating the heart and respiration as coupled oscillators, we then apply a method based on Bayesian inference to find the underlying coupling parameters and their time dependence, deriving from them measures such as synchronization, coupling directionality and the relative contributions of different mechanisms. We report a detailed analysis of the reconstructed cardiorespiratory coupling function, its time evolution and age dependence. We show that the direct and indirect respiratory modulations of the heart rate both decrease with age, and that the cardiorespiratory coupling becomes less stable and more time-variable.

Introduction
The cardiorespiratory interaction was first observed by Hales [1]. Since then, the associated physiology and physics have been under continuing investigation. The mechanisms are still subject to debate [2], with the two prevailing approaches focused mainly on the heart rate 2013 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/ by/3.0/, which permits unrestricted use, provided the original author and source are credited. variability (HRV) that occurs as a consequence of the interaction. One approach supposes that baroreflex plays a key role [3,4], whereas the other is based on the hypothesis of sympathetic/parasympathetic balance [5].
In this work, we adopt a different perspective. We treat the cardiac and respiratory functions as being oscillatory and interacting [6], and we consider the coupling to be enabled by two distinct mechanisms. The first is via the vasculature and the flow of blood through the system of closed tubes that connects the heart and the lungs. The second coupling occurs via the neuronal control that provides information transmission between all of the systems involved, i.e. the heart, lungs and vasculature. The oscillatory processes associated with the coupling mechanisms have been identified in earlier studies, and are related to endothelial [7][8][9][10], neurogenic [9,11,12] and myogenic [13,14] activity. We measure the oscillatory functions of heart and respiration, and then apply methods drawn from nonlinear science to detect their instantaneous phases and frequencies. Once these have been established, we study the interaction between heart and respiration, treating them as self-sustained nonlinear oscillators and investigating their mutual coordination [15][16][17] or synchronization [18][19][20][21][22][23][24][25][26], coherence [27][28][29][30][31][32] and direction of coupling [33][34][35][36], allowing for the fact that the characteristic frequencies, couplings and other interacting characteristics vary in time [37,38]. We report below an analysis of the results obtained from a study of 189 (99 male, 90 female) healthy subjects aged 16-90, enabling us to characterize the changes with age that occur in the cardiorespiratory interaction.
In §2, we present the methods of phase and frequency extraction used in the study, and we show that the synchrosqueezed wavelet transform (SWT) [39] offers significant advantages. In §3, we apply wavelet phase coherence (WPC) to investigate the components modulating the heart and respiratory rates. We formulate a model to describe cardiorespiratory interactions in §4. In §5, we extract the time-varying parameters of the model using recently developed methods based on Bayesian inference, and we investigate the age dependence of the derived measures: the directionality of coupling, synchronization and the overall coupling function. Finally, in §6, we summarize and draw conclusions. Details of subjects, measurements, signal pre-processing and statistical analyses are given in the appendices.

Phase and frequency extraction
To obtain the phase and frequency [40] of the respiratory and cardiac signals, we used the recently developed SWT [39]. The SWT T f (ω, t) is constructed from the wavelet transform (WT) W f (a, t) of the signal as where W l = [ω l , ω l+1 ) specifies a frequency bin of central frequency ω * l = (ω l + ω l+1 )/2 and width ω l = ω l+1 − ω l . For each wavelet scale a and time t i , we calculate the corresponding frequency where arg denotes the phase of the complex coefficient, and division by 2π is needed to go from circular frequency to the natural frequency in hertz. The sum in (2.1) is over scales a such that the calculated ω f (a, t) lies in the frequency bin W l at time t.
The WT W f (a, t) is calculated using the Morlet mother wavelet [41, p. 70], and we use scales a j = 2 j/n v (2/f s ), such that each dyadic interval is divided into n v equi-log-spaced bins (n v is chosen such as to provide appropriate numerical resolution; we use n v = 64). In the latter expression, f s denotes the sampling frequency and j is such that a j ranges from 2/f s to T/8.6 (so there should be at least eight oscillations for us to regard the apparent wavelet component as being physically meaningful), with T being the overall time of the signal. For the SWT (2.1), we use similar equilog-spaced bin edges ω l whose spacings are determined by the same n v = 64, although we now investigate only the frequency region within which the corresponding first harmonics reside:  Figure 1. Synchrosqueezed wavelet transforms of (a) the ECG and (b) respiration. The support of the main curve, i.e. the selected region, from which phase and frequency were extracted, is shown by the grey lines. The oscillations in the ECG curve are at the respiration frequency and correspond to respiratory sinus arrhythmia. (Online version in colour.) The advantage of the SWT over the WT is that it is far more compressed than the WT, which enables one to distinguish more accurately between different components in the signal. It becomes possible to select only one component (e.g. the main oscillation) and to filter out the others (e.g. higher harmonics). The signal s(t) can then be reconstructed from its SWT as where C ψ is a normalization constant [39]. We follow the support of the main oscillation (e.g. the first harmonic of the ECG), selecting a region of width n v /2 = 32 in frequency bins around it, adjusted so as to contain only the selected curve and no others. This enables us to extract the phase φ(t) and frequency f (t) as and where the sums in (2.2) correspond only to the selected region. Extracting the frequency from the whole region proved to be more accurate than just taking it to be the centre of the main curve's support, because the SWT curves are often non-uniformly broadened.
The main difficulty lies in automatically tracking the support of the main curve. At time t, we define it to be the region of width n v /2 having the maximum functional M(n, t), withn denoting its central bin: Here, ωn is the central frequency of the current window, ω prev is the frequency (as in (2.2)) of the previous window (the corresponding term is used to suppress jumps) andω is the mean frequency, extracted as in (2.2), but over the whole range of ω (the corresponding term is used to give more stability to the main curve). The Gaussian term is needed to give more weight to the central part of the window, i.e. if there are two or more peaks in one window, to select the support of one of them. Despite its complexity, the method provides for good flexibility and the possibility of accurate tuning. Based on the characteristic variations of the fundamental frequencies, the optimal parameters were determined to be λ = 5, κ = 0 for ECG and λ = 10, κ = 10 for respiration. Figure 1 shows examples of the extracted supports of the SWT curves corresponding to the ECG and respiration first harmonics. The advantage of using SWT for phase and frequency extraction over the more commonly used marked events (MEs) method is that one obtains phase and frequency values at the same sampling frequency as the signal, whereas the ME method gives values only at particular points (e.g. peaks). Thus, the effective sampling frequency of the characteristics extracted by MEs is, in fact, the (time-varying) frequency at which the MEs occur, i.e. the reciprocal of the current interpeak interval (usually around 1 Hz for ECG and 0.25 Hz for respiration); all other values, obtained by interpolation to match the original sampling frequency, are unphysical.
Note that heart activity is sometimes modelled through use of the integral pulse frequency modulation (IPFM) model [42] which, in effect, assumes that only the R-peak positions in the ECG are informative. Within this framework, much effort has been directed towards the reconstruction of instantaneous heart frequency (IHF), treated as the input signal. These methods involve the analysis of an RR interval series or a pulse-train equivalent of the ECG [43], from which the continuous IHF is obtained by low-pass filtering, spline interpolation or other techniques [42,44]. With such an approach, the interbeat information is necessarily neglected, so that the IHF cannot be obtained with a physical sampling frequency higher than that at which the R-peaks occur. In this sense, such methods can be regarded as variations on the ME method, but with different implicit or explicit interpolation schemes. They provide a very good approximation within the assumption that IHF modulation can be represented by a few sinusoids of constant amplitude and frequencies no higher than half of the mean heart frequency [44][45][46]. However, there exist many input signals that can give the same peak positions, so one cannot say anything for sure about the IHF between peaks without also taking account of interbeat behaviour. By contrast, the SWT uses all of the information, in effect yielding the IHF with a high effective sampling frequency. For example, in the case when all RR-intervals are the same, but the T-wave occurs at different times after the R-peak, methods based on the IPFM model assumption will yield a constant heart frequency, whereas the IHF extracted from the SWT will be slightly modulated, thus providing an indication of these small changes. Additionally, SWT frequency extraction does not require a high ECG sampling frequency, in contrast to the above-mentioned methods [47]. It is also more universal in that it is applicable to almost any oscillatory signal and not only to series of point events.
The advantage of the phase extracted from SWT over the commonly used Hilbert phase is that the former can be already regarded as a true phase, whereas the latter is a protophase (see [40] for definitions and a discussion of the differences). This is because the SWT phase is extracted from the fundamental frequency and thus is uncorrupted by higher harmonics. As a result, it is almost unaffected by the protophase-to-phase transformation [40], and it does not change (up to a constant phase shift) if we apply an invertible nonlinear transformation to the oscillatory component. For example, it will be the same for sin φ(t) and exp[sin φ(t)], as well as for the respiration signals as simultaneously measured by belt and by thermistor, whereas the corresponding Hilbert phases will differ. Figure 2 shows typical phase and frequency results extracted by this method in comparison with the values extracted by other methods. The advantage over Hilbert phase extraction is clearly evident in figure 2c,d; in figure 2c, we illustrate the well-known fact that Hilbert phase extraction is inapplicable to the ECG, whose phase is usually obtained using the ME method. For clarity, the ME phase is not shown, but the difference is clearly visible in figure 2e,f, where the frequencies extracted by SWT and ME are compared. It should be noted, however, that one might apply the Hilbert transform to the signal after bandpass filtering in the region where the first harmonic resides, and also obtain a well-behaved phase; but the frequency calculated from it will not, in practice, be accurate. Moreover, when there are additional frequency components present, or where the main frequency varies by more than a factor of two (as is often the case for respiration), the frequency range of the first harmonic will inevitably contain additional components or harmonics. By contrast, SWT phase extraction uses a narrow but timedependent adaptive frequency range, picking up only the main curve and thus avoiding the above problems.
Analysing the frequency extracted by SWT, one recovers the well-known result that heart rate becomes less variable with age [48][49][50]. Figure 3 shows how the standard deviation of the extracted frequency depends on age. Spearman's rank correlation coefficient ρ and its significance p (see appendix C) are also shown. From figure 3a, we see that the standard deviation of heart frequency has a negative correlation with age, i.e. the heart rate becomes less variable. There is no comparable effect in the case of respiration (figure 3b). The mean heart and respiratory frequencies do not change significantly with age (not shown), except for heart rate in males (ρ = −0.27, p = 0.01).  (c-f ) The phases and frequencies of the cardiac (φ h , f h ) and respiratory (φ r , f r ) components, extracted from the signals by different methods. In (c,d), the phases extracted by the SWT are compared with those extracted by the Hilbert transform; the phases extracted by the ME method, apart from a constant phase shift in the ECG case, are almost indistinguishable from phases extracted by SWT within the resolution of the figure, and therefore are not shown. In (e,f ), the frequencies extracted by SWT are compared with those found by the ME method; the frequency calculated from the Hilbert phase is usually inaccurate and seldom used. In both our cases, it deviates considerably from the SWT and ME frequencies.

Wavelet phase coherence
To investigate which oscillatory components are in common to both the respiration and heart activities, we calculated the WPC [27][28][29][30][31] between the IHF and the instantaneous respiration frequency (IRF), extracted from ECG and respiration SWTs, as well as between IHF and respiration. Note that the IHF is commonly, and somewhat confusingly, referred to in the literature [51] as HRV.
We calculated the WPC between the respiration signal and IHF, and between IRF and IHF. In doing so, we considered the frequency range 0.0095-2 Hz, which contains five physiologically relevant intervals [41]. To avoid bias owing to different lengths of physically relevant WT for different frequencies, we used the time interval of the WT corresponding to the lowest frequency (0.0095 Hz) for all scales. Note that, usually, one cannot properly consider the WT at frequencies greater than or approximately equal to 0.5 Hz for IHF (on account of Nyquist's theorem) because IHF extracted by the ME method has effective sampling frequency equal to the cardiac frequency (see above). In the present case, however, with SWT frequency extraction, we can go beyond this limit. Thus, in our case, the IRF and IHF have the physical sampling frequency of 50 Hz, so we can, in principle, go up to 25 Hz.  V  IV  III  II  I  V  IV  III  II  I   phase  For oscillations with lower frequencies, there are fewer cycles within a given time interval. Their phase coherence is therefore expected to be larger, but without implying a higher level of interdependence. We therefore need to find an appropriate threshold (significance level) for each frequency, above which the phase coherence may be regarded as physically meaningful and indicates genuine interdependence.
To estimate significance levels, we use intersubject surrogates. Signals from different subjects must obviously be independent, while having similar characteristic properties. We therefore calculate surrogate values of, for example, the IHF-respiration phase coherence using respiration signal from one subject and the IHF from another. The significance level is then estimated as the 95th percentile of 189 surrogate values. Note that we use a significance level estimated by direct rank ordering (see [52, p. 352]) rather than based on the number of standard deviations (sigmas) above the mean of the surrogates, as is sometimes done. The latter criterion is not always appropriate because it implicitly assumes normally distributed surrogate values, which is untrue, in general, and for our case, in particular.
Use of other surrogate types, such as the widely used FT [52,53] or IAAFT [52,54] surrogates, is also possible. In our case, they may be used to test against the hypothesis that the apparent phase coherence obtained is entirely attributable to bias related to the frequency content and to the finite duration of the signal (and possibly to a distribution of values in the case of IAAFT). One might argue that our intersubject surrogates can sometimes have a lower significance level than (IAA)FT surrogates owing to the different breathing and cardiac frequencies for each subject. Although it appears not to be the case, it would not represent a problem even if true. We wish to test the presence/absence of dependence, rather than being interested in frequency bias in particular, so that the existence of exact matching of, for example, the mean frequency of respiration to the respiratory component in the IHF should not be neglected as is implicitly done by (IAA)FT. Thus, intersubject surrogates seem to represent the more natural choice for our case, as well as being convenient and computationally cheap. Nevertheless, given the high time-variability of cardiorespiratory signals, the WPC suffers bias only owing to the finite time length of the time series, so that (IAA)FT and intersubject surrogate tests become almost equivalent. We have found that they all give the same result in terms of significant WPC peaks. Figure 4 shows the results for all subjects, but we have found that the same considerations apply to males and females separately. In the respiration-IHF phase coherence (figure 4a), there is a highly significant peak in the respiratory frequency range, corresponding to the well-known respiratory sinus arrhythmia (RSA). It was also observed previously in the cardiorespiratory coherence function, although the corresponding significance levels were not always estimated reliably (see [53] for discussion and references therein). We have observed a decrease with age in the respiratory peak WPC (not shown: ρ = −0.38 for males and ρ = −0.35 for females with p < 0.001 for both). Coupled with the fact that the respiratory frequency and its standard deviation do not depend on age (figure 1b), this suggests that RSA becomes more unstable with age.
In the respiration-IHF WPC, there is also a significant peak at the cardiac frequency, which could not have been observed using the earlier ME-or IPFM-based methods of IHF extraction. It suggests that IHF and respiration share a cardiac component. Nevertheless, it is actually a measurement artefact, because the apparent cardiac component in respiration is attributable to the direct influence of heart beats on the respiratory measuring belt, which is not removed completely by the smoothing procedure (see appendix B). However, the presence of a cardiac component in the IHF is something new and interesting. It is clearly evident in the wavelet and SWTs of the IHF (not shown), but will not be discussed further here.
We find no significant coherence at all between the IHF and IRF ( figure 4b). This implies that IHF and IRF do not, in general, share any frequency component, i.e. there is no common influence that measurably affects both the heart and respiration rhythms.

Model of cardiorespiratory interaction
We model the cardiorespiratory system as a pair of noisy coupled oscillators with phases φ h,r , where 'h' implies heart and 'r' implies respiration: The time-derivatives of the phasesφ h,r /2π correspond to IHF and IRF, ω h,r denote natural frequencies, and ξ h,r (t) is assumed to be white Gaussian noise: ξ i (t)ξ j (τ ) = δ(t − τ )E ij . The coupling functions q h,r are separated into three parts: a self-interaction part s h,r which accelerates/slows phase growth depending on the current phase of the oscillator considered; a direct interaction d h,r which accelerates/slows phase growth depending on the current phase of the other oscillator; and an indirect part i h,r which depends on both phases. Note that all terms are allowed to be time-dependent. The self-interaction s h,r reflects the effect of the oscillator on itself. In the case of using the protophase, for example, extracted by Hilbert transform, it could arise simply owing to nonlinearity, i.e. the contribution from the higher harmonics [40]. However, from the definition of phase, it follows that for the true phase the self-interaction can be only coupling-induced [40], being attributable to the interaction with other systems.
As implied by its name, the direct interaction d h,r describes the direct influence, or driving of one oscillator by the other. We will refer to it as the 'RSA-term' because it is the only term that can be responsible for the stable respiratory modulation of IHF (and the corresponding phase coherence seen in figure 4a). This is because clear respiratory frequency modulation ofφ h can be established only by the part dependent solely on φ r .
The last term, i h,r (φ h , φ r , t), describes a more complicated coupling, dependent on both phases. It is called indirect, because if the heart and respiratory activities both interact with some other system, for example, are each coupled to a third oscillator, this will be included in our model (4.1) as some complicated term(s) dependent on both heart and respiratory phases, and so be included within i h,r . Obviously, respiration and heart always interact through some other subsystems but, in the case of indirect coupling, we mean that the intermediate system is influenced by both respiratory and cardiac activity and does not merely transfer signals from one oscillator to the other. Nevertheless, the indirect interactions can also make a contribution to the direct and selfinteraction parts. Such contributions are called coupling-induced and are mainly responsible for s h,r in our case. Let us illustrate the meanings of s h,r , d h,r , i h,r by consideration of the mechanisms of cardiorespiratory interaction. It is well known that the respiratory centre modulates vagal motoneuron responsiveness ('respiratory gating') and thus heart rate [2,3]. Such 'driving' oḟ φ h depends solely on the respiratory phase φ r , thus being completely included within the direct modulation term d h (φ r , t). Other influence onφ h comes from the baroreflexes, which feel arterial pressure and induce heart rate changes in response to its variations [55]. Thus, heart, arterial pressure and respiration interact in a closed loop [56,57]. Because arterial pressure is modulated by both heart and respiration [58], depending on both φ h,r , this mechanism will be included mainly within the indirect term i h (φ h , φ r , t). As already mentioned, it can contribute to s h,r , d h,r as well.
Our model (4.1) includes explicitly only the phases of the cardiac and respiratory activities. In general, IHF and IRF are modulated by other processes as well, for example, modulation of IHF at low frequencies around 0.1 Hz [5,35,[48][49][50]56]. Nonetheless, such additional mechanisms are effectively accounted for by allowing the terms in (4.1) to be time-dependent. At the same time, consideration of other signals such as arterial pressure can provide additional insights into the underlying dynamics. Thus, in [57], the authors studied in addition arterial pressure, making it possible to separate the influence of respiration on IHF into the effect owing to 'respiratory gating' and that owing to baroreflex coupling with arterial pressure fluctuations (modulated in part by respiration). Here, however, we are mainly interested in the age dependence of the coupling parameters and not so much in revealing the physiological origins of the interactions.

Extraction of nonlinear interactions using Bayesian inference
The phase dynamics of cardiorespiratory interactions can be analysed by application of a technique based on Bayesian inference [59][60][61][62]. To extract the parameters of interaction (ω h,r , q h,r ), we apply the Bayesian method extended to incorporate time variability [37,38] to the phases extracted from the SWT. Thus, we model the right-hand sides of (4.1) as a sum of chosen basis functions multiplied by some coefficientsφ h,r = C h,r;1 f 1 (φ h , φ r ) + C h,r;2 f 2 (φ h , φ r ) + · · · + ξ h,r (t), and find most probable values of these coefficients C h,r;i in each time window.
Because all functions of phase, such as q h,r (φ h , φ r , t), should be 2π -periodic over each phase, the most appropriate basis functions are Fourier series. We used sin(mφ h − nφ r ), cos(mφ h − nφ r ) with n = −2, −1, 0, 1, 2 and m = −2, −1, 0, 1, 2, i.e. Fourier series up to second order. In addition, we used synchronization terms with n : m = 3 : 1, 4 : 1, 5 : 1, 6 : 1, 7 : 1, 7 : 2, 9 : 2, 10 : 3, corresponding to the common cardiorespiratory synchronization ratios [63,64]. Altogether, for each ofφ h,r , we used 41 basis functions. We used non-overlapping windows of time length 50 s, chosen to incorporate at least 10 respiratory cycles and thus provide enough information for accurate inference. We perform inference within each window, so all parameters are determined on this time scale (e.g. the synchronization duration is 'quantized' in units of 50 s). The propagation constant [38] was chosen by trial-and-error basis to be p w = 20. Its significance relates to the propagation of information, i.e. the extent to which inference in the current time window takes account of the result obtained in the previous one. For completely independent inference in each window, p w = ∞, whereas p w = 0 means full propagation. The choice of p w does not significantly influence the results.
With the chosen basis functions, the self-interaction (s h,r ) and direct interaction (d h,r ) parts of the coupling function are modelled using sin(φ h,r ), cos(φ h,r ), sin(2φ h,r ), cos(2φ h,r ) (φ h for s h , the IHF and IRF; all other terms depend on both phases and thus generally cannot have purely respiratory or cardiac frequency. It should be noted that incorporation of higher terms, such as Fourier series up to fourth order, does not lead to results qualitatively different from those presented below.

(a) Synchronization
Having determined for some window the most probable coefficients C h,r;i multiplying the basis functions, we can find out whether there is n : m synchronization between φ h and φ r within that window. To do so, we first fix ψ n,m ≡ φ h /n − φ r /m, and for each value integrate (4.1) with inferred right-hand sides, so that φ h /n + φ r /m go from 0 to 2π . Let us denote the value of φ h /n − φ r /m at the end of this procedure as M(ψ n,m ), because it depends on the initial choice of ψ n,m . Then, most probably there is synchronization if there exists a ψ n,m such that M(ψ n,m ) = ψ n,m and, for synchronization to be stable, |dM(ψ n,m )/dψ n,m | < 1. If there is such a root, then we consider there to be synchronization throughout the duration of the selected window, and otherwise no synchronization within that window. A major advantage of the Bayesian approach to the investigation of synchronization is that it is more noise robust [38] than methods based on the synchrogram [64] or on different synchronization indices [65][66][67]. This is because through the inference of noise parameters, one can explore the uncorrupted underlying dynamics. A detailed comparison of approaches is a separate topic, but we comment that we obtained qualitatively the same results by the application of the phase coherence method [65]. Figure 5a summarizes all the information about synchronization obtained in this way, in the form of the relative synchronization duration, i.e. the proportion of the time during which the signals were synchronized at a particular n : m ratio. Note that, because synchronization is quantized in 50 s intervals, the results do not take into account possible shorter episodes of synchronization. Therefore, in general, the relative synchronization duration will actually be higher. Second, the calculated overall synchronization duration relates only to the particular synchronization ratios considered, taking no account of, for example, 9 : 1 synchronization. As can be seen, however, even the 6 : 1 and 7 : 1 ratios are quite rare, whereas the most prevalent ratio is usually 4 : 1 (figure 5c).
No monotonic age dependence of the overall degree of synchronization has been found (figure 5b), and nor has it for synchronization at any particular ratio (not shown). Nevertheless, it seems that certain synchronization ratios may be age-and gender-dependent, for example, 10 : 3 synchronization arises mainly in males aged around 40 (figure 5a,d). However, a more detailed investigation, involving larger numbers of subjects of given age and gender, would be needed to draw strong conclusions about particular age/gender groups (see appendix C).

(b) Directionality
We calculate the directionality of influence as D(t) = (c r→h (t) − c h→r (t))/(c r→h (t) + c h→r (t)), where c r→h (t) is calculated as the norm of the coefficients ( C 2 1 (t) + C 2 2 (t) + · · ·) of the basis functions that depend on φ r (i.e. correspond to d h and i h ) and represents the influence of respiration on the heart; correspondingly, c h→r (t) represents the influence of the heart on respiration. Figure 6 shows the dependence on age of the time-averaged c h→r , c r→h and the resultant coupling directionality D together with its standard deviation. It is clear from figure 6a that the influence of respiration on the heart has significant negative correlation with age, whereas the influence of the heart on respiration is not age-dependent as shown in figure 6b. As a result, the time-average of coupling directionality significantly decreases with age (figure 6c). At the same time, its standard deviation increases with age (figure 6d), implying that the cardiorespiratory interaction also becomes less stable.   Figure 7.
(a-f ) Dependence on age of the different contributions to coupling directionality (see text).
To ascertain whether the decrease of influence of respiration on heart (and thus directionality) reflects only a decrease of RSA, or whether it implies something more, we now consider the relative contributions made by different terms in the directionality. Let us denote by c   Figure 7 shows time-averages of the different contributions to the cardiorespiratory interaction, as well as their mutual relationships. It is evident that both the direct and indirect influences of respiration on the heart are negatively correlated with age ( figure 7a,b), with the direct contribution decreasing faster than the indirect one (figure 7c). At the same time, contributions to the influence of the heart on respiration and the corresponding indirect/direct ratio exhibit no significant age dependence (figure 7d-f ). Thus, the decreases with age of c r→h and coupling directionality are completely attributable to the combined effects of a decrease in RSA and the indirect effect of respiration on the heart.
We also analysed in detail the direct and self-interaction terms (5.1), characterizing their shape variations by the amplitude ratio B(t)/A(t) between two harmonics in (5.1). We found a significant increase with age in the time-averaged B/A ratio (ρ = 0.38, p = 0.00 for males and ρ = 0.29, p = 0.01 for females), implying that the form of RSA becomes less harmonic with age. We also found a significant increase in the standard deviation of this ratio (males: ρ = 0.34, p = 0.00; females: ρ = 0.28, p = 0.01) and a decrease of the phase coherence γ = | e i(2α(t)−β(t)) | between harmonics in (5.1) (males: ρ = −0.29, p = 0.00; females: ρ = −0.33, p = 0.01). This indicates an increase in RSA time-variability with age, consistent with the result of figure 6d, and a decrease in respiratory WPC with age (see §3). No significant correlations with age were observed in the cases of s h , s r , d r .

(c) Coupling function
To gain further insights into the nature of the cardiorespiratory interaction, we now analyse the overall form of the reconstructed coupling functions q h,r (φ h , φ r , t) (4.1). Figure 8 shows the time-averaged versions of the coupling functions q h,r typical of a younger and an older subject. Decrease of the RSA amplitude with age is clearly seen in figure 8a,b. It can also be concluded that the main stable contribution to q h , surviving after time-averaging, remains the RSA irrespective of age. The respiratory coupling q r shown in figure 8c,d seems to be quite irregular and not age-dependent.  However, the time-averaged coupling functions provide only limited information. The full time evolution of the same coupling functions are illustrated through videos. 1 From them, one can see that, in older people, the heart coupling function becomes also less stable in time, dominated by the highly time-variable indirect contributions. At the same time, the dynamics of q r does not seem to change with age, being irregular and unstable.

Summary and conclusions
Our investigations have yielded results in agreement with those reported earlier [41,[48][49][50], while also revealing many new features that invite a reconsideration of some of the phenomena observed previously. In particular are the following.
-Frequencies. Applying the SWT [39] to measured time series, we have obtained the physically relevant instantaneous frequencies of the heart (IHF, commonly referred to as HRV in the literature) and respiration (IRF) with a sampling frequency of 50 Hz, as well as their corresponding phases. In agreement with earlier works [41,[48][49][50], we observed a decrease of IHF variability with age (figure 3a), but no significant correlation with age in IRF variability (figure 3b) or mean respiratory frequency (not shown). -Coherence. Computing the WPC ( §3) and estimating its significance threshold by the use of intersubject surrogates, we found statistically significant peaks at both the respiratory and heart frequencies in the respiration-IHF WPC (figure 4a). The former is attributable to RSA and is in agreement with earlier findings [53], whereas the latter represents a measurement artefact (see appendix B). There is also a decrease with age in the respiration-IHF phase coherence at respiratory frequencies, indicating that RSA becomes less stable. There is no significance in the IRF-IHF WPC (figure 4b), indicating that the heart and respiratory rates are modulated mainly by separate mechanisms. -Bayesian inference. We model the cardiorespiratory system as a pair of noisy coupled oscillators with explicitly time-dependent coupling parameters. Heart and respiratory coupling functions were divided into three parts, representing: coupling-induced selfinteraction; direct driving by the other system; and indirect interactions. The direct influence of respiration on the heart was attributed to 'respiratory gating' [3], i.e. RSA, whereas the indirect influence was related to baroreflex coupling [55] [38,[59][60][61][62], we have obtained the forms and time evolutions of the underlying heart and respiratory coupling functions (figure 8, see also footnote 1), confirming and extending the preliminary results of [38]. They were then used to analyse synchronization, directionality of coupling and the contributions of different mechanisms to the cardiorespiratory interaction. -Synchronization. There is no significant correlation of the overall synchronization duration with age (figure 5b), consistent with recent findings [64] based on the synchrogram method. The most common cardiorespiratory synchronization ratio was found to be 4 : 1 ( figure 5c,d). It seems also that certain synchronization ratios may be characteristic of particular ages and genders (figure 5a,d). -Coupling directionality. The overall influence of respiration on the heart decreases with age (figure 6a), whereas influence in the opposite direction stays constant (figure 6b), leading to a net decrease of coupling directionality with age ( figure 6c). Directionality also becomes more time-variable with age ( figure 6d), corresponding to a progressive impairment of the underlying mechanisms. -Cardiorespiratory interactions. Consistent with earlier findings [41,[48][49][50] the RSA amplitude, i.e. the strength of the direct respiratory modulation of heart rate, decreases with age (figure 7a), as well as changing form and becoming more time-variable (see last paragraph of §5b). Because the breathing rate does not change with age, these effects cannot be attributable to the influence of the breathing frequency on RSA [68,69]. The indirect modulation of the heart by respiration also decreases with age (figure 7b), though not so rapidly as RSA (figure 7c). This effect may correspond to a decrease in baroreflex sensitivity with age [70,71]. The strengths of the direct and indirect influences of heart on respiration (figure 7d,e) do not change with age. -Coupling functions. Underlying all these separate effects, the heart coupling function changes markedly with age, both in its average form (figure 8a,b) and in its timevariability (see footnote 1), whereas the respiratory coupling function seems to be irregular and unaffected by age (figure 8c,d, and see footnote 1).

Appendix B. Filtering and pre-processing
After signals with significant movement artefacts or frequent ectopic beats (more than 50 during the 30 min recording) had been discarded, we were left with 189 sets of data (99 from males, 90 from females). The ECG and respiration signals from each subject were pre-processed as follows. The first and last 1 per cent of each signal were discarded to eliminate the starting and ending transients present in some of them. Next, those signals sampled at 1200 Hz were resampled to our standard 400 Hz. All signals were then detrended using a 200 s window. The belt used to measure respiration 'feels' the heart beats directly, and consequently introduces a spurious heart component into the respiration signal. This effect can be clearly observed during a breath hold. To partly eliminate it, the respiratory signals were smoothed by use of a 1 s moving average. All signals were then resampled to 50 Hz, matched in length, and de-meaned. These procedures were applied carefully, in such a way as to preserve the exact time relationship between each pair of signals.

Appendix C. Statistical methods
To quantify the significance of the observed monotonic age dependence, we have used Spearman's rank correlation coefficient ρ (see [41, §C.2.2]), and (two-tailed) statistical significance was assumed for values of p < 0.05. The correlation coefficient for all males, all females and the merged group, as well as its significance, is presented in each figure rounded to two decimal places. For example, 'p = 0.00' implies p < 0.005.
Because the Spearman's correlation coefficient only takes account of monotonic dependence, the subjects' datasets were divided into age bins, and medians over each bin are shown in the figures in order to better eliminate possibly nonlinear age dependence. The age bins are centred on 20, 30, 40, 50, 60 and 70 years and, except for the last group, are each of 10 years duration (e.g. the 50-year group includes subjects aged 46-55). The final bin includes all subjects older than 65 (up to 90). Age distributions for each gender are presented in figure 9. Note, however, that comparisons between age bins of particular gender may have not enough subjects for good statistics, given the complexity of the cardiorespiratory system and large individual influences. Therefore, our main conclusions are all based on the Spearman's correlation and its significance, which take into account all 99 male and/or 90 female subjects, whereas the information over age bins is provided mainly as a guide to the eye.