## Abstract

Activity of many physiological subsystems has a well-expressed rhythmic character. Often, a dependency between physiological rhythms is established due to interaction between the corresponding subsystems. Traditional methods of data analysis allow one to quantify the strength of interaction but not the causal interrelation that is indispensable for understanding the mechanisms of interaction. Here we present a recently developed method for quantification of coupling direction and apply it to an important problem. Namely, we study the mutual influence of respiratory and cardiovascular rhythms in healthy newborns within the first 6 mo of life in quiet and active sleep. We find an age-related change of the coupling direction: the interaction is nearly symmetric during the first days and becomes practically unidirectional (from respiration to heart rhythm) at the age of 6 mo. Next, we show that the direction of interaction is mainly determined by respiratory frequency. If the latter is less than ≈0.6 Hz, the interaction occurs dominantly from respiration to heart. With higher respiratory frequencies that only occur at very young ages, the dominating direction is less pronounced or even abolished. The observed dependencies are not related to sleep stage, suggesting that the coupling direction is determined by system-inherent dynamical processes, rather than by functional modulations. The directional analysis may be applied to other interacting narrow band oscillatory systems, e.g., in the central nervous system. Thus it is an important step forward in revealing and understanding causal mechanisms of interactions.

- direction of coupling
- development
- transmission dynamics

interaction between physiological (sub-) systems can be addressed by the analysis of the signals they generate. Traditional signal-processing techniques, e.g., cross-spectral (coherence) estimates (20) or computation of mutual information (25), provide symmetric measures of interaction strength. There were also several attempts to address the causality in interaction, based on theoretical information approach (11, 34) and methods of nonlinear dynamics (27, 33). We study directional interrelation using a technique that is particularly suited for the case of weakly interacting rhythmical systems. Thus we exclude from our consideration the input-output systems, i.e., we do not treat the case when, say, *signal 1* can be regarded as delayed and/or (nonlinearly) filtered *signal 2*. We also do not consider the situation when an intervention is possible, i.e., when causal relation can be addressed by application of certain stimuli and analysis of the response to them (for description of such techniques, see, e.g., Ref. 31). In our approach, we restrict ourselves to the case when *1*) the system can be only passively observed and *2*) it can be modeled by two (or several) weakly coupled self-sustained oscillators, i.e., subsystems, having their own rhythms. Such coupled-oscillator models describe a variety of biological phenomena (9, 10, 24, 38). In the framework of coupled-oscillators theory, interaction (coupling) means that the systems mutually perturb their rhythms. It is of importance that these rhythms exist even if the systems were completely uncoupled, i.e., the rhythms exist not due to interaction but are inherent properties of systems. In this context, we understand directionality as the quantity that characterizes which of the two systems influences its counterpart more strongly. This quantity, or directionality index, can be obtained by the method (29, 30) based on the coupled-oscillators theory (24). Particularly, the analysis of directionality can provide new insights into mechanisms of interaction between the respiratory and cardiovascular control systems.

Two centuries ago, Carl Ludwig described modulation effects of the respiratory rhythm on the cardiovascular rhythm (17). His explanations of this observation focused on direct mechanical peripheral coupling. Later it became evident that respiratory and cardiovascular rhythms may influence each other due to central nervous mechanisms (14-16, 28) as well as interaction mechanisms in the periphery (23). Nowa-days, the modulation effect (respiratory sinus arrhythmia) is well-studied (1). Another effect is synchronization (coordination) between respiratory and cardiovascular systems; it has been described and quantified (18, 22, 32). This phenomenon can be assigned to the presence of weak interaction that adjusts average frequencies of the heart beat and respiration to a ratio *n*ω_{1} = *m*ω_{2}, where *n* and *m* are integers.

Although synchronization and modulation indicate the presence of cardiorespiratory coupling, there are no data about the direction of this interaction so far. This is mainly due to the limitations of classical tools of time series analysis that are not able to quantify a bidirectional (but not necessarily symmetric) interrelationship from two signals. In our work, we address this problem and investigate the direction of interaction of heart rhythm and respiration in 1,717 artifact-free data segments (1,355 quiet sleep, 351 active sleep) from 22 healthy sleeping babies who exhibited a broad band of spontaneous breathing frequencies. The age and frequency dependence of the new index was explored by analyzing respiratory signal and heart rhythm derived from the ECG in a cohort setup. In a previous study (29) we analyzed different methods of detection of bidirectional relationships in a small group and found that the evolution map algorithm (EMA) is well suited for our purposes because of its robustness and simplicity.

We emphasize that modulation and synchronization are different (although not independent) aspects of cardiorespiratory interaction. From the standpoint of nonlinear dynamics and data analysis, this difference is as follows. Modulation (e.g., heart rate variability) can be described and analyzed within the input-output systems framework (31). On the contrary, coupling that leads to adjustment of rhythms, and direction of this coupling in particular, can be understood only within the nonlinear coupled-oscillators theory.

## MATERIALS AND METHODS

*Subjects and protocols.* We measured the ECGs using a bipolar limb lead (Biomonitor 501, Meßgerätewerk Zwönitz) and obtained thoracic respiratory effort using the inductive plethysmographic method (Respitrace, Studley Data Systems, Oxford, UK) in 22 newborn infants during sleep. Measurements were performed on each of the first 5 days of life, then every week and later monthly up to the 6th mo of life. Informed parental consent and local ethics committee approval were obtained.

Data acquisition began 30-60 min after feeding, in the evening hours between 8 PM and 11 PM, and took ∼1 h. Data were stored on a DAT multichannel recorder (DTR-1800, biologic) for further analysis. Periods of quiet and active sleep were determined using criteria of motor behavior (26). Irregular respiration, closed eyes, rapid eye movement, and gross movements of the extremities defined periods of active sleep. Quiet sleep was defined by regular respiration, closed eyes, no apneic episodes, only occasional twitches and other motor activities, and no eye movement or rhythmic movements of the mouth.

*Identifying directionality of coupling.* The detailed mathematical procedures are given in the appendix.

*Data processing.* The data were offline-digitized with a computer-based monitoring system. Respiratory signals were sampled at 10 Hz, and R-waves were detected with the precision of 1 ms (hardware device based on slope and amplitude criterion, equivalent to 1 kHz). Artifact-free, 2-min-long segments of quiet and active sleep of each measurement were chosen for the further analysis. The length of the segments was chosen to optimize the tradeoff between the requirement to have more data points for significant estimates and to have the segment sufficiently short to assume its stationarity (Fig. 1). In total, 1,717 data segments have been analyzed in this study.

From the series of instances when R-wave peaks appear on the ECG, the instantaneous phase θ_{h} of the cardiac signal has been taken as where *t*_{k} is the time of *k*th R-wave occurrence in the ECG. The instantaneous phase of the respiratory signal θ_{r} has been obtained using the concept of the analytic signal based on the Hilbert transform (see Ref. 24 for discussion on phase estimation from data). The Hilbert transform has been applied to the whole segment of respiratory signal. The latter technique requires the removal of the low-frequency baseline oscillations. Therefore, the respiratory signal was first high-pass filtered with a 3rd-order Savitzky-Golay filter with the window width of 2.5 s. To reduce high-frequency noise, the signal was then low-pass filtered with a second-order Savitzky-Golay filter with the window width of 0.2 s.

The presence of interaction between cardiac and respiratory systems, a prerequisite of directionality estimation, is indicated by the presence of respiratory sinus arrhythmia as well as by the results of synchronization analysis performed on the same group of subjects (18). To trace the changes in coupling degree and/or directionality imposed by slow drift of the coupling parameters in the system, the index of directionality *d*^{(h,r)} has been computed from the time series of phase data (θ_{h},θ_{r}). The parameter τ was taken as the average cardiac cycle length within each segment; the results were stable with respect to variation of τ in a wide range.

*Statistics.* To categorize according to age, all measurements were grouped into six age groups. The age limits were chosen in a way that each group contained approximately the same number of measurements (286 ± 59 total, 58 ± 14 active sleep, 228 ± 49 quiet sleep). We analyzed whether the directionality index depends on the age, sleep stage, or both factors by means of the nonparametric Brunner-ANOVA (3). The confidence limit α was set to 0.05.

## RESULTS

There is an evolution from an approximately symmetric interaction during the first days of life toward a dominant unidirectional (from respiration to the heart rate) coupling at the age of 6 mo (Fig. 2). There is no sleep stage dependence and no sleep stage-specific age dependence of the directionality index.

The direction of interaction depends strongly on respiratory frequency. For all measurements with a respiratory frequency of <0.6 Hz, the interaction occurs dominantly in a unidirectional fashion, from respiration to heart rhythm. For breathing frequencies >0.6 Hz, no clear direction can be identified (Figs. 3 and 4).

## DISCUSSION

The question of which of either systems is the dominating driver for the interaction has been debated for many years; however, no clear answer has emerged. Most authors assume that there is a unidirectional direction of interaction, namely, that respiration acts on the heart rate. This is also reflected by the widely used terminology “respiratory sinus arrhythmia.” The origin of this terminology is based on the teleological argument (14) that the respiratory rhythm is essential for life, whereas fluctuations of heart rate are dispensable. The physiological background of the generation of sinus arrhythmia can be found in a number of scenarios (4), such as direct thoracic mechanical coupling (17), pulmonary vagal afferents (5), baroreflex modulation (7), and central vagal outflow modulation (8).

To shed light on the directionality problem, we have estimated direction of interaction of respiration and heart rate during quiet and active sleep stages within the first 6 mo of life in humans. This setup allowed us to investigate the cardiorespiratory system at different levels of development as well as at different sleep stages, representing various physiological conditions.

We have made two interesting observations. First, around the age of 6 mo, there is a tendency toward a unidirectional interaction from respiration to the heart rate, whereas in the first days of life there is no dominating direction. This age dependency may be related to the age dependency of respiratory frequency (19). Second, for all measurements with a respiratory frequency of <0.6 Hz, the interaction occurs dominantly in a unidirectional fashion, from respiration to heart rhythm. At respiratory frequencies >0.6 Hz the dominating direction is abolished. The mechanism for the latter interesting fact may be explained by the low-pass behavior of the vagal-atrial transmission dynamics. The dynamics include the release of acetylcholine from the nerve endings, its diffusion across the synaptic cleft, the action on sinoatrial pacemaker cells, the dynamics of signal transduction in the pacemaker cells, and the degradation and reuptake processes of acetylcholine. Early observations show a frequency dependency of respiratory modulation of heart rate (1). The phenomenology of the resulting low-pass like behavior has been modeled (6) as well as explored experimentally, e.g., by the analysis of respiratory sinus arrhythmia and voluntarily controlled breathing frequency at fixed tidal volume (13). In their work they determine a corner frequency of 7.2 ± 1.5 breaths/min, which corresponds to 0.12 ± 0.025 Hz. (The corner frequency refers to a -3 db damping.) At higher breathing frequencies, the damping increases, because the dynamics of the transduction system is too sluggish to follow. We argue that at 0.6 Hz, a transmission of information is essentially reduced. Hence, if the respiratory frequency is high, the attenuating effect is too strong, and the dominating direction of interaction is abolished.

Many functional parameters, such as spectral power of breathing and heart rate during sleep, depend on sleep stage (18, 19, 21, 35). However, sleep stage does not determine the value, the age (Fig. 2), and the breathing frequency dependence (Fig. 3) of the directionality index. These three findings suggest that sleep stage is not a decisive factor for the directionality index and that the mechanism behind them is not related to the well-known modulation of function of the autonomic nervous system in relation to sleep stages. Thus the observed relationships are most probably caused by the kinetics involved in the vagal-atrial transmission.

In summary, we conclude that there is a clear dominating effect of respiration on heart rhythm at the age of 6 mo. At high breathing rates >0.6 Hz, the interaction is not unidirectional, most probably due to low information transmission to the cardiac oscillator caused by the low-pass behavior of the information transmission from the vagal nerve to the atrial pacemaker cells.

The directional analysis may be applied to other interacting narrow band oscillatory systems, e.g., in the central nervous system. Thus it is an important step forward in revealing and understanding causal mechanisms of interactions.

## APPENDIX

*Identifying directionality of coupling.* Here we briefly present the technique developed in Refs. 29 and 30. Qualitatively, it can be explained in the following way. Suppose *system 1* affects *system 2*, but there is no back-action. Then in the power spectrum of an observable signal from *system 2*, we will find at least two peaks: in addition to the peak at the average frequency of *system 2*, we will see the peak at the frequency of *system 1*. In case of bidirectional interaction, the spectrum obtained from each system will contain also frequency components of the counterpart. Directionality then will quantify the relative power of all these peaks. The approach can be also explained in terms of mutual predictability that goes back to the Granger concept (11). If *system 1* affects *system 2*, then the future of *system 2* can be better predicted if we take into account the history of both *systems 1* and *2* (and vice versa). Below we discuss the mathematical framework that allows one to quantify directionality by one number. We emphasize that before application to physiological data, the method was tested with simulated data; it was shown that directionality can be reliably revealed from rather short noisy records. The technique was also verified in a specially designed physical experiment with electronic circuits, where periodic and chaotic regimes have been tested, as well as the influence of the external noise (2).

Although the mechanisms underlying a particular physiological rhythmic behavior might remain largely unknown, macroscopically observed oscillations such as breathing excursions can be represented by the dynamics of a nonlinear self-sustained oscillator. For the simplest case of periodic, limit cycle oscillations, one can introduce the phase variable θ in a unique way: , where ω = 2π*/T* is the natural (angular) frequency and *T* is the period of the oscillation. A remarkable property of the phase variable is its sensitivity to small external perturbations to the oscillator's dynamics. This leads to an important idea of the coupled oscillators theory: weak coupling merely influences the phase dynamics, without altering each oscillator's intrinsic limit cycle. This property underlies the possibility of reducing the description of weakly interacting oscillators to a phase model. Complex (aperiodic) oscillations admit a similar description if a suitable phase variable can be defined (24).

For the case of two weakly coupled oscillators, of our interest here, the phase model can be written as 1 where θ_{1} and θ_{2} denote the continuous phase variables [defined on the whole real line, not on the (0,2π) circle], and the random terms ζ_{1,2} account for amplitude fluctuations in aperiodic (chaotic) oscillators and/or random perturbations that are always present in the real-world systems. The coupling is described through the 2π-periodic functions *f*_{1,2}, and the strength of interaction is given by the parameters ϵ_{1,2}, with ϵ_{1,2} << ω_{1,2} in the assumption of weak coupling. In case of unidirectional driving configurations (e.g., if *system 1* drives *system 2)*, ϵ_{1} = 0 and ϵ_{2} ≠ 0. It means that the evolution of θ_{2} reflects the bivariate dependence *f*_{2}(θ_{1,}θ_{2}), whereas the derivative of θ_{1} depends solely on θ_{1} (Fig. 5). In bidirectional coupling configurations, derivatives of θ_{1} and θ_{2} depend on both phase variables. Provided that in a certain experimental situation, phase variables can be retrieved from the measured time series data, and the assumption of weak coupling is supported by our knowledge of the system, how can we detect the specific coupling configuration of the two oscillators? Here we briefly present our approach to answer this question. We use the EMA (see Refs. 29 and 30 for details). The idea was to observe the evolution of the phase of each oscillator over a certain temporal window of length τ, during which the phase increments are determined both by the natural frequencies ω_{1,2} and the effect of coupling to the other oscillator (see Fig. 5 for illustration). Such phase increments can be considered as generated by an unknown two-dimensional noisy map that can in principle be derived from the continuous time model *Eq. 1*. We estimate the coupling terms *F̄*_{1,2} of the unknown map from time series θ_{1,2}. As *F̄*_{1,2} are periodic with respect to phases, we approximate them by a finite Fourier series

In our computations we have used *m,n* < 4. The coefficients *A*_{m,n} can be obtained by fitting in the least mean square sense the dependencies of Δ_{1,2} on θ_{1,2}. Note that fitting also filters out the noise.

Further, we analytically calculate the cross-dependency coefficients given by These coefficients are integrative measures of how strongly an oscillator is driven and how sensitive it is to the driving. Finally, a directionality index has been introduced Normalized in this way, the value of *d* ranges from -1 in the case of unidirectional coupling (2 → 1) to 1 in the opposite case (1 → 2), while intermediate values correspond to bidirectional coupling configurations (Fig. 6). If *c*_{1} = *c*_{2}, then two systems affect each other to the same extent, and the interaction is symmetrical. The index of directionality characterizes the asymmetry in coupling terms ϵ_{1,2}·*f*_{1,2}(θ_{2,}θ_{1}) and does not reflect the difference in the frequencies ω_{1,2} of autonomous systems (see analytical treatment in Ref. 29 and illustration in Fig. 7). The algorithm fails in the phase-locked regime of noise-free coupled oscillators, what mathematically corresponds to the appearance of a functional relationship between the two phase variables. In the population analyzed here, synchronous epochs were short and sporadic (18); therefore, the application of the algorithm was not limited by this fact.

Concerning the choice of the parameter τ, we mention that in the weak coupling (ϵ_{1,2} << ω_{1,2}) case, the effect of coupling becomes detectable at time scales large or comparable to the characteristic oscillatory period, while at small time scales the evolution of phase is dominated by noise.

Alternative algorithms for weak coupling directionality identification have been proposed and discussed in Ref. 29. The EMA has been exploited here due to its robustness, simplicity, and low computational demands.

## Footnotes

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2003 the American Physiological Society