Vol. 275, Issue 3, H1092-H1102, September 1998
MODELING IN PHYSIOLOGY
One-dimensional, nonlinear determinism characterizes heart rate
pattern during paced respiration
Katrin
Suder1,
Friedhelm R.
Drepper1,2,
Michael
Schiek1,2, and
Hans-Henning
Abel3
1 Institut für
Biologische Informationsverarbeitung and
2 Zentrallabor für
Elektronik, Forschungszentrum Jülich, 52425 Jülich; and
3 Abteilung
Kardioanästhesie, Städtisches Klinikum, 38302 Braunschweig, Germany
 |
ABSTRACT |
This study focuses on the dynamic pattern of
heart rate variability in the frequency range of respiration, the
so-called respiratory sinus arrhythmia. Forty experimental time series
of heart rate data from four healthy adult volunteers undergoing a
paced respiration protocol were used as an empirical basis. For
pacing-cycle lengths >8 s, the heartbeat intervals are shown to obey
a rule that can be expressed by a one-dimensional circle map
(next-angle map). Circle maps are introduced as a new type of model for
time series analyses to characterize the nonlinear dynamic pattern
underlying the respiratory sinus arrhythmia during voluntary paced
respiration. Although these maps are not chaotic, the dynamic pattern
shows typical imprints of nonlinearity. By starting from a piecewise linear model, which describes the different circle maps obtained from
the empirical time series for various pacing frequencies, time
invariant measures can be introduced that characterize the dynamic
pattern of heart rate variability during voluntary slow-paced respiration.
sinus arrhythmia; cardiovascular model; periodicity; system
physiology; nonlinear dynamics
 |
INTRODUCTION |
THE VARIATION OF heart rate in the frequency range of
respiration, known as respiratory sinus arrhythmia (RSA), was already described by Ludwig in 1847 (18). Since then, many studies have been
carried out to explore the underlying physiological mechanisms (17,
19). The most important ones are the modulation of cardiac filling
pressure by respiratory movements (1), the direct respiratory modulation of parasympathetic and sympathetic neural activity in the
brain stem (10), and the respiratory modulation of the baroreceptor
feedback control (9).
There are different models summarizing the comprehensive results of
the RSA analysis. Whereas mechanistic models (6, 22, 23, 26,
28) support the qualitative understanding of the interplay of the
RSA-generating mechanisms, time series analyses (2, 23, 24, 28) result
in quantitative estimates that are phenomenologically related to the
strength of autonomic cardiac chronotropic innervation. Estimates based
on linear time series analyses focus on amplitudes in either the time
or the frequency domain (3, 15, 21). They can be used to quantify the
strength of coupling between the respiratory and the cardiovascular
system and, in particular, to estimate the sympathetic/parasympathetic balance (8).
As has already been pointed out by various authors (5, 16, 20), the
dynamics of cardiorespiratory synchronization features essential
nonlinearities that cannot be described by using these estimates. We
hypothesize that the characteristic pattern of successive heartbeat
intervals within one respiratory cycle (pattern of the RSA) can be used
to quantify additional information about the nonlinear dynamics of the
cardiac chronotropic control.
The theory of nonlinear dynamic systems provides many methods for the
identification and classification of the dynamic pattern in an observed
time series. Because biological and physiological systems are
characterized by stable dynamic structures far from equilibrium, such
methods can be used to analyze these systems (7, 11).
Complexity measures have been introduced to differentiate distinct
geometric structures that are representations of the dynamics in an
abstract (embedding) space (12, 16, 30). This differentiation was
possible even in those cases in which the dynamic structure (the
attractors) cannot be identified. One intrinsic drawback of this
approach is the difficulty of giving a physiological interpretation of
the complexity measures. However, these measures are already of
clinical relevance. With regard to cardiology, these measures are
discussed in the context of risk stratification of sudden cardiac
death. With a significance level of ~90%, Skinner et al. (27) showed
that, for the heartbeat intervals, a reduction of the correlation
dimension (a special complexity measure related to the dimension of the
attractor) precedes ventricular fibrillation in high-risk patients by
hours. However, most applications of nonlinear time series analyses in
cardiology do not pay special attention to the cardiorespiratory
coupling, which manifests itself in the RSA.
The complexity of spontaneous respiratory movements obviously precludes
the identification of finite dimensional attractors in heart rate
fluctuations within the frequency range of breathing (14). One way to
overcome this problem is to introduce paced breathing. This
standardized condition has already been used in physiological and
clinical studies to analyze and quantify the RSA (4, 15, 20). The main
result of the present paper is that heart rate fluctuations during
voluntarily induced slow-paced breathing obey a one-dimensional,
nonlinear law of motion. Neither the structure of the one-dimensional
deterministic process nor the deviation from it depends on the mean R-R
interval length or the RSA amplitude. However, they do depend on the
respiratory cycle length and vary significantly between individuals.
For all subjects, the one-dimensional law of motion breaks down for
cycle lengths close to that of spontaneous breathing.
The one-dimensional structure becomes apparent when the dynamics of the
amplitudes of the RSA is separated from its form. This can be achieved
by introducing circle maps. The circle map results from the
introduction of polar coordinates into the two-dimensional scatter plot
(embedding) of successive heartbeat intervals. Circle maps have been
used previously (13) to study the oscillatory behavior of externally
stimulated heart tissue.
The different circle maps can be approximated by a piecewise linear
model specified by a few parameters. From these parameters, new
measures characterizing respiration-related heart rate variability can
be derived and interpreted from a physiological point of view.
 |
METHODS |
Subjects.
Four men, with ages between 29 and 35 yr and without any history of
cardiopulmonary diseases, participated in this study. All subjects were
nonsmokers for >1 yr. The measurements were taken between 9:00 AM and
2:00 PM, 2-3 h after the subject's last meal or caffeinic
beverage. Care was taken to ensure that subjects were in a relaxed
state. All subjects were informed about the purpose of the study and
gave their consent before the examination.
Experimental design.
The subjects were asked to breathe voluntarily in phase with a growing
(inspiration) and vanishing (expiration) light circle at 10 different
cycle lengths (TT = 5, 5.5, 6, 6.5, 7, 8, 9, 10, 11, and 12 s) for 70 cycles at each frequency. To
check for transient phenomena, the pacing frequencies were arranged to
also include large frequency jumps. The ratio of expiration to
inspiration length was kept constant (4:3). The tidal volume was not
controlled because the subjects were allowed to adjust their alveolar
minute ventilation to maintain the acid-base homeostasis. Note that
tidal volume is considered to have a minor influence on the RSA (15). The paced respiration experiment was split into three sessions of
nearly equal duration (Table 1). The
interval between sessions was 45 min to 3 days. During each session,
electrocardiogram (ECG, bipolar lead I) and respiratory flow
(uncalibrated thermistor signal) were recorded with a sample rate of
1,000 Hz with the use of MacLab 8/s (ADInstruments) in connection with
a Macintosh Quadra 650 personal computer. All measurements were done in
a head-up, 45° tilt position. This body position was chosen as a compromise between a relaxing body position and a position that supports a relatively high degree of attention.
Data.
The cardiac cycle length (R-R interval) was determined algorithmically
with an accuracy of ~1 ms as the time interval between two
successive maxima of R waves of the QRS complexes of the original ECG.
The data set contains only one extrasystolic heartbeat. In a previous
study, the same data set was analyzed with the use of various methods
including symbolic dynamics and the analysis of phase shifts
(24).
To focus on the relevant timescale of the RSA [the respiratory
cycle length (TT)], a
high-pass filter was applied to the heartbeat lengths. The filter was
implemented as the difference to a moving average with an averaging
period equal to the respiratory period. This procedure removes
low-frequency nonstationary baseline trends. A typical piece of the
high-pass filtered R-R interval time series, recorded at
TT = 9 s, is shown in Fig.
1A.

View larger version (28K):
[in this window]
[in a new window]
|
Fig. 1.
Epochs (165 s) of R-R interval time series of subject
1 (A) and
subject 2 (C) during paced respiration with a
respiratory cycle length (TT)
of 9 s (data are high-pass filtered with cutoff equal to
TT). Deviations from mean
heartbeat interval ( ; 1,045 ms
for subject 1 and 870 ms for
subject 2) are shown. Compared with
subject 1, who exhibits a clustering
of heartbeats at long R-R intervals, subject
2 shows clusters of heartbeats at short R-R intervals.
B: epoch of R-R interval time series
computed from a realization of the second surrogate process
(surrogate 2; preserving spectrum and
histogram), which is based on R-R interval time series of
subject 1, partially shown in
A. Compared with heartbeats of
physiological time series (A), which cluster at long R-R
intervals, surrogate heartbeat lengths show generally little clustering
and less regularity. Time offsets of recordings have been set to 0.
|
|
Respiratory flow data were used to check the breathing performance of
the subjects. As can be seen from Table 2,
the relative standard deviation of the
TT intervals is generally small.
It is caused mainly by short periods of reduced attention to the optical metronome.
 |
RESULTS |
From time series to circle maps.
Typical features of the RSA dynamics of a single healthy subject are
the pronounced variability of the amplitude of the RSA (defined as the
difference between the maximal and minimal R-R interval within one
respiratory cycle) and the considerable stability of the relative
partitioning of the respiratory cycle by heartbeat intervals (2, 23),
which we call the characteristic pattern of the RSA dynamics. Comparing
the R-R interval time series of two subjects (Fig. 1,
A and
C), this pattern becomes visible: whereas subject 1 has a clustering of
heartbeats at longer R-R intervals, subject
2 has clusters of heartbeats at shorter intervals. Another feature is the asymmetry between the number of heartbeats with
increasing and decreasing R-R intervals, respectively (2, 22). In this
paper, we introduce a more general method for describing and
quantifying the RSA dynamics.
A first attempt at identifying the underlying process generating this
stable and characteristic pattern of the RSA is to plot the same data
as a two-dimensional scatter plot of successive R-R intervals
(embedding of R-R intervals) (Fig. 2).

View larger version (30K):
[in this window]
[in a new window]
|
Fig. 2.
Two-dimensional embedding of R-R intervals recorded from
subject 1 during 70 cycles of paced
respiration (TT = 9 s). R-R
intervals are taken as deviations from mean value
= 1,045 ms. Part of data is
presented in Fig. 1A. Points scatter
around an elliptical structure and do not fall on a discernible
one-dimensional curve. This finding is generally observed when
analyzing R-R interval time series recorded during paced respiration.
n, Time step.
|
|
However, no sharply defined pattern is revealed in this plot; the
points scatter around an elliptical structure. All R-R interval time
series contain this feature, but they do not show a simple nonrandom
structure. The elliptical structure mirrors the basic periodicity of
the R-R intervals corresponding to the periodicity of respiration.
As a next step in the analysis, we suppress the fluctuations of the RSA
amplitude and focus instead on the angular motion of the points in the
plane. This is done by introducing polar coordinates and by neglecting
the fluctuations of the radii. The transformation requires the
definition of a center (radius = 0) (cf. Fig.
3B).
After application of the above-mentioned high-pass filter on the R-R
interval time series, the mean R-R interval length
(
) is suited to define the
center.

View larger version (26K):
[in this window]
[in a new window]
|
Fig. 3.
Graphic explanation of embedding of rotation angles. R-R intervals
L1, ..., L12 occurring
during 1 respiratory cycle TT
(A) are embedded two-dimensionally,
leading to 12 pairs of successive R-R intervals
(L1,L2),
..., (L12,L13)
(B). Rotation angles ( ) are
defined around a center point
( , ).
Maximal heartbeats are close to = 1/4, and minimal heartbeats are
close to = 5/4. To get a continuous graph of the circle map
(C), we shifted all ordinates
(angles at time n + 1) larger than 1.5 by 2.0. Because of the periodicity of angles , this does not
change the dynamics, but it eases visual analysis of all circle
maps.
|
|
For every pair of successive R-R intervals, an angle
is calculated.
Each rotation angle
corresponds to a unique phase of the RSA cycle.
However, because this phase is not identical to the respiratory phases
(inspiration and expiration), we reserve the term "phase" for the
physiological phase. The respiration-related angles are defined to vary
between 0 and 2, corresponding to 0
2
|
(1)
|
where
n is the angle at time step
n,
Ln is the R-R
interval at time step n, and
is the mean R-R interval length.
The transformation is illustrated in Fig. 3 for one respiratory cycle.
The function shown in Fig. 3C is
called a circle map or next-angle map.
Identifying a one-dimensional deterministic process.
If we transform the R-R interval data from subject
1, as shown in Fig. 2, to polar angles
(Eq. 1) and plot the data as a
next-angle map, this two-dimensional embedding of the angles leads to a
one-dimensional curve, and the points cluster close to the line (Fig.
4A).
This curve describes the deterministic dynamics of the rotation angles around the origin in the
Ln+1
versus Ln plane.
The law of motion, which is represented by the graph in Fig. 4, allows us to read off the next angle when the preceding one is known. From
comparison of the scattering of points in the two-dimensional embedding
of the angles (Fig. 4A) with the
scatter plot of the R-R intervals (Fig. 2), it is obvious that the
dynamics of the angles are much less affected by noise than are the
dynamics of the R-R intervals. Thus the restriction on the dynamics of
the angles leads to the disclosure of a one-dimensional deterministic law of motion. In Test for nonlinear deterministic
process using surrogate data, we show formally that
this visual impression is legitimate. Furthermore, this
circle map will be shown to have essential nonlinear features.

View larger version (21K):
[in this window]
[in a new window]
|
Fig. 4.
Two-dimensional embedding of rotation angles based on data from Fig. 2.
A: original high-pass filtered data
(subject 1, 70 cycles,
TT = 9 s).
B: surrogate R-R interval data with
complete destruction of autocorrelation but preserved histogram.
C: amplitude-adjusted surrogate R-R
interval data with preserved histogram and spectrum (a corresponding
time series of R-R intervals is shown in Fig.
1B). Both surrogates originate from
data shown in A.
|
|
As Fig. 5 shows, a one-dimensional law of
motion is typical for cardiorespiratory dynamics during slow-paced
breathing. For TT > ~7 s, a
staircaselike structure of the circle map is observed for all subjects
under examination.

View larger version (32K):
[in this window]
[in a new window]
|
Fig. 5.
Two-dimensional embedding of rotation angles based on R-R intervals
recorded at TT = 10 s
(A), 8 s
(B), and 6 s
(C). Subject
3 (left) is compared
with other subjects (right), who
have similar mean R-R interval length and mean respiratory sinus
arrhythmia (RSA) amplitude. Although
A-C show the same overall
one-dimensional, staircaselike structure, they show individual
differences in scattering and pattern.
|
|
However, a comparison of the circle maps for different subjects (Fig.
5, from left to
right) and for different pacing
cycle lengths (Fig. 5, from A to
C), also establishes distinct
interindividual differences. Nonetheless, the example pairs with the
same TT, shown in Fig. 5, are
similar in mean R-R interval length (deviation <5%) and mean RSA
amplitude (deviation <15%). Different degrees and locations of
clustering are observed: whereas subject
3 exhibits a clustering of successive rotation angles
preferentially near the maximum
(
n = 0.3),
subject 2 shows more clustering near
the minimum (
n = 1.3). This
could already be seen in Fig. 1. A second difference is the varying
degree of scattering (e.g., Fig. 5, B
and C), indicating that the
transition to the one-dimensional dynamics occurs at different
respiratory cycle lengths and to a different extent. Finally,
subject 3 shows a characteristic kink
between
n = 0.0 and
n = 0.4 for all cycle lengths
that is absent in all other subjects. It is remarkable that there are
cases in which the circle maps of different subjects differ
qualitatively, whereas the commonly used descriptive variables of RSA
analyses, the mean R-R interval, and the mean RSA amplitude do not
differ substantially. Thus we conclude that the analysis of the
next-angle map provides additional information.
Test for nonlinear deterministic process using surrogate data.
To demonstrate whether the analysis of successive rotation angles (Fig.
4A) in fact reveals a deterministic
structure and does not simply result from the embedding or the
high-pass filtering, we used surrogate data from Monte Carlo
realizations of stochastic processes, which preserve certain linear
properties of the original time series (25, 29).
To estimate the effect of the embedding, we first tested the hypothesis
that the time series is generated by a white noise stochastic process.
On the basis of this hypothesis, surrogate 1 data are obtained by random shuffling of the
empirical R-R intervals. This leads to a process that preserves the
histogram.
In the recurrence plot of rotation angles for
surrogate 1 data, the originally
observed structure is clearly lost (Fig.
4B). This means that this first
hypothesis must be rejected, i.e., the one-dimensionality seen in Fig.
4A is not an artifact of the transformation used. However, the embedding does have an effect, because
n and
n+1 are
not independent variables (Eq. 1).
Only part of the plane is filled, and the apparently deterministic
points at (0/
0.5) and (1/0.5) are introduced by the embedding.
To show that the rejection of surrogate
1 data is not only due to the complete destruction of
the autocorrelation, we generated a second surrogate process that does
preserve both the spectrum (i.e., the linear autocorrelation) and the
histogram. This second (null) hypothesis is equivalent to the
assumption of an underlying linear autoregressive process of unlimited
order that is observed by a nonlinear measurement function. The
nonlinearity is defined by the non-Gaussian property of the original
process (25). Because of the high-pass filtering of the original data,
the effective dimension of the linear autoregressive model is of the
order of the number of heartbeats per respiratory cycle. This kind of
autocorrelation-adapted surrogate process generates a spectrum that
cannot be distinguished from the original one.
Compared with the scatter plot of the original data (Fig.
4A), the one-dimensional structure
in the recurrence plot of rotation angles for
surrogate 2 data (Fig.
4C) appears less obvious.
Furthermore, the asymmetry between the minimum and the maximum is lost.
Figure 1 shows the comparison between the original time series
(A) and an example of this surrogate
process (B).
To get a quantitative estimate of the distinctness of the two plots, we
performed a
2 test using a sum
of squared residuals as a discriminator. The residuals were chosen as
the deviations of the data from an optimally adapted piecewise linear
model consisting of six line segment (see Piecewise
linear model). This special choice of
six lines is suggested by the data presented in Figs. 4 and 5. Two of
the segments are chosen a priori as fixed, and the eight parameters of
the other four segments are estimated independently for each surrogate
process. Because of the large number of data points, the sampling
distribution of the square root of the
2 was assumed to be normal. The
mean and average of this distribution were estimated from a sample of
10 surrogate data sets. The corresponding
2 test results in the rejection
of the second null hypothesis with a probability of <0.001. The
square root of the
2 of the
original process is 7.75 standard deviations from the average of the
corresponding surrogate cases. This means that the second null
hypothesis must also be rejected. Furthermore, the fluctuations of the
empirical heart rate data are influenced by nonlinear deterministic
dynamics, which can be described significantly better by a piecewise
linear model than by a linear autoregressive moving average model with
roughly the same number of parameters.
In summary, we identified a one-dimensional, nonlinear deterministic
process underlying the dynamics of the R-R intervals during slow-paced
respiration.
Piecewise linear model.
As mentioned in Test for nonlinear deterministic
process using surrogate data, a piecewise
linear model can be used to quantify the structure of the individual
circle maps and their specific changes due to the protocol parameter
(respiratory cycle length). As Fig. 6
shows, the model consists of six straight lines (12 parameters), two of
which are chosen a priori to have a fixed distance to the bisector (
= 0.5) and fixed slopes parallel to it. This reflects the deterministic
transitions due to the embedding (line segments 3 and
6; see Fig. 4B). To fix the
remaining four lines of the model, we must specify eight parameters.
One possible choice comprises the six positions of the corner points on
the abscissa and two vertical distances of those corner points to the
bisector: x1,
x2,
x3,
x4,
x5,
x6,
h2, and
h5. The four
fixed parameters are
h1 = h4 = 0.5 and
h3 = h6 = 0.5 (Fig. 6;
for details see Ref. 28).

View larger version (23K):
[in this window]
[in a new window]
|
Fig. 6.
Piecewise linear circle map with 6 line segments. Two lines are chosen
a priori with fixed slopes
(a3 = a6 = 1.0) and
fixed distances
(h1 = h4 = h3 = h6 = 0.5) to
bisector
[ai, slope
of line i;
xi, coordinate
of corner point i,
hi, vertical
distance of corner point i to line of
identity]. Note that angles vary between 0 and 2, corresponding
to 0 2 , and that parameters chosen to describe the model
are presented in larger type size.
|
|
The free model parameters were estimated by iterative linear regression
for all respiratory cycle lengths and all subjects. For simplicity,
only the error of the ordinate was taken into account. For example, for
subject 1 at
TT = 9 s, the following parameter values were determined:
x1 = 0.066, x2 = 0.248, x3 = 0.814, x4 = 1.030, x5 = 1.388, x6 = 1.948, h2 = 0.051, and h5 = 0.089.
For all cycle lengths and subjects, the estimated parameters correspond
to invertible circle maps. Because maps generating chaotic dynamics are
noninvertible (11), it is clear that we do not encounter deterministic
chaos in our analysis.
There is a second embedding effect that can be used as a consistency
check for the estimation. This second effect can best be understood
from Fig. 2. There are two ways to obtain the number of heartbeats
below
: one counts the angles in
either the range from 0.5 to 1.5, corresponding to the cases in which Ln <
, or the range from 1.0 to 2.0, corresponding to Ln+1 <
. Because the two sets are
related by a time shift, both numbers must be equal. This consistency check was easily passed by 39 of 40 estimations. This additional constraint reduces the number of independent model parameters from
eight to seven.
As a second test for the validity of the model, the root mean square
error of the fit was calculated for all 40 time series (Fig.
7). This error is a measure of how well the
different line segments approximate the empirical data. It
is particularly small for high pacing length. In the range between 8 and 12 s, the angle dynamics (and therefore the pattern of the RSA) are
described well by the one-dimensional model. For respiratory cycle
lengths <8 s, the mean square error increases strongly; the
one-dimensional description starts breaking down (cf. Fig. 5,
B and
C). This result is in agreement with
the identification of two different cardiorespiratory synchronization
regimes (5-7 s and 8-12 s) within the same data set using
symbolic dynamics (24).

View larger version (22K):
[in this window]
[in a new window]
|
Fig. 7.
Root mean square error of piecewise linear model optimally adapted to
each empirical data set of all subjects. Model is a very good
approximation in range between 12 and 8 s. For shorter
TT, error increases strongly. If
a uniform distribution of points in squares (cf. Fig. 4B) is assumed, a
mean square error of 5/12 can be calculated.
|
|
In addition, the piecewise linear model defined herein can be used to
reconstruct a time series on the basis of the estimated model
parameters (28).
Characteristics of heart rate pattern.
To be able to quantify the pattern of the heartbeat intervals within
one respiratory cycle, it is preferable to construct measures that are
sensitive to all possible characteristics of the pattern. Furthermore,
it should be possible to obtain these measures directly from the
original time series as well as from the model introduced in
Piecewise linear model. The invariance of the law of motion of the angles (Eq. 1) implies that the pattern of the RSA dynamics has
invariant characteristics that can be used as descriptors. Because the
model is defined by seven independent parameters, the maximal number of
descriptors should not be greater than that. It is convenient to
introduce the following characteristic measures.
Asymmetry index 1 (A1) describes
the relative asymmetry between the number of heartbeats with interval
lengths longer than
compared
with those shorter than
. As can be seen from Fig. 2 or 3, this corresponds to the two ranges of angles,
0 to 1 and 1 to 2. Asymmetry index 2 (A2) describes
the relative asymmetry between the number of increasing and decreasing heartbeat intervals (asymmetry of changes). Because an adjacent pair of
R-R intervals close to
= 1/4 corresponds to a local maximum of R-R
intervals, and because the pair close to
= 5/4 corresponds to a
local minimum, A2
is defined as the difference between the number of angles in the
interval 1/4
5/4 and the remaining number, divided (normalized)
by the total number of angles. Both of these indexes are useful in
characterizing the asymmetry of the map. The asymmetry is a direct
expression of the nonlinearity of the dynamics.
Two other measures,
C1 and
C2, are used to
quantify the degree of clustering at the accumulation points.
Accumulation points are points near the bisector, because the fixed
points of a one-dimensional map are on the bisector. Empirically, the
degrees of clustering can be obtained as the logarithmic ratio of the
maximal step length of all angles compared with the average minimal
step length for heartbeat intervals either longer than
(0 to 1) or shorter than
(1 to 2), respectively. The main
effect comes from the variation of the minimal step length, because the maximal step length can be approximated by 0.5 (cf. Fig. 6).
C1 is a measure
of the clustering of angles near the maximal R-R interval, whereas
C2 corresponds to
clustering near the minimal R-R interval.
Two other invariants,
x2 and
x5, specify the
degree to which the accumulation points of the circle map coincide with
the maximum or the minimum of the R-R intervals at
= 1/4
(corresponding to a maximum) and
= 5/4 (corresponding to a
minimum).
The winding number (W) is commonly
used to describe the average number of R-R intervals in one respiratory
cycle. The notion of a winding number results from viewing the combined
angular dynamics of the two oscillators, heart and respiration, as a
motion on a torus. W is the ratio of
the average number of faster oscillations (windings) per single cycle
of the slower oscillator.
Note that all these measures describing the coupling between
respiration and heart rate can be determined without measuring respiration. They can be obtained directly from the experimental time
series, as described in METHODS, and
also from the piecewise linear model, e.g., from a reconstruction of
time series based on the model. Another useful feature of the model is
that all the characteristic measures have simple dependence on the
model parameters (see APPENDIX).
Figure 8 shows a comparison of both
approaches, the direct determination from the empirical time series and
the approximation based on the model parameters. The standard deviation of the direct measures was obtained by subdividing the time range into
three intervals. Therefore, it can be taken as a measure to describe
the variance in time.

View larger version (32K):
[in this window]
[in a new window]
|
Fig. 8.
Comparison of empirically obtained characteristic measures (experiment)
with those of analytic approximation (model) for
subject 1: asymmetry
index 1, asymmetry index
2, degree of clustering
1, degree of clustering
2, deviation of the accumulation points to minima
(x5 1.25), and winding number. Standard deviation of empirically obtained
measures describes variance in time.
|
|
The agreement between both approaches is reasonably good for
W. The disagreement for the other
measures results from imperfections of the model fit as well as from
the analytic approximation described in the
APPENDIX. Both of these are of the
same order of magnitude. However, despite this disagreement, both
measures allow a clear distinction of the frequency dependence of the
next-angle map of the four subjects. The differences are particularly
pronounced for the degrees of clustering (Fig.
9).

View larger version (15K):
[in this window]
[in a new window]
|
Fig. 9.
Comparison of 7 experimentally obtained degrees of clustering
(C1 and
C2).
A: subject
2; B:
subject 3;
C: subject
4. All subjects show a distinct frequency dependence
with different saturation level for long
TT.
|
|
The following characteristic changes in
C1 and
C2 are observed.
For respiratory cycle lengths shorter than 8 s, both measures show
similar changes. They increase monotonically with almost the same
slope. W shows a similar behavior. For
long cycle lengths, all subjects show a saturation effect. However, the
subjects differ in their saturation level: subject
4 shows a high clustering for short heartbeat
intervals, subject 2 has a high degree
of clustering for long heartbeat intervals, and
subjects 1 and
3 have both measures nearly equal.
Because of their discriminatory properties and stabilities, C1 and
C2 might serve as
robust characteristics for a subject, in addition to the
asymmetry index 2 (2).
Together, these seven measures encode the systematic information about
the pattern of heart rate variability during slow-paced respiration. In
the validity range of the model, these measures are guaranteed to be
time invariants. Furthermore, for respiratory cycle lengths shorter
than 8 s, these measures are still invariant in time.
 |
DISCUSSION |
The main result of the present study is that heart rate variability
during paced respiration at longer cycle lengths (8-12 s) obeys a
dynamic rule that can be expressed by a one-dimensional, nonlinear
circle map (next-angle map). This map describes the rotation dynamics
in the two-dimensional embedding space spanned by two successive
heartbeat intervals. These dynamics represent properties of the
characteristic pattern of R-R intervals within one respiratory cycle.
The existence of this law of motion implies that the relative
partitioning of the respiratory cycle by heartbeat intervals is
invariant in time compared with the amplitude of the RSA, for which
such an invariance could not be demonstrated (28).
Data presented by the empirically obtained circle maps for four healthy
adult volunteers breathing at 10 different respiratory cycle lengths
led us to define a piecewise linear, phenomenological model that
captures essential features of heart rate dynamics. The description of
heart rate variability by this nonlinear model was contrasted with its
description by a linear autoregressive model with the same number of
parameters. The nonlinear model provides a significantly better fit to
the empirical data than to the surrogate data. This is further proof
that the fluctuations of the empirical heart rate data have essential
nonlinear dynamic features (5, 16, 20). The parameters of the piecewise
linear model were related to characteristic time invariants. Two of
them, A1 and
A2, are direct
measures of the nonlinearity of the model.
These characteristic invariants, quantifying properties of the pattern
of heartbeat intervals within one respiratory cycle, can also be
obtained directly from the R-R interval time series without using the
model. Their algebraic relationship to the model parameters guarantees
their invariance in time. This time invariance holds not only for the
longer respiratory cycle lengths, for which the one-dimensional
deterministic process is identified, but also for cycle lengths in the
range of spontaneous breathing, for which no similar deterministic
one-dimensional structure can be detected (14).
The observation of one-dimensional dynamics for slow pacing of the
cardiorespiratory system is indeed surprising, because this system is
known to be influenced by physiological processes on many different
timescales, including longer timescales such as circadian rhythm or
hormonal cycles with cycle lengths ranging from minutes to hours. A
transition to one-dimensional dynamics can only be expected for
periodically forced systems with a maximal time constant of all
internal feedback mechanisms that is either shorter than, or of the
same order of magnitude as, the cycle length of the driving system
(entrainment of faster internal degrees of freedom). For nonlinear
coupled oscillators, the frequency range of the entrainment (1:1 Arnold
tongue) can be relatively broad. The fact that such a transition to
one-dimensional dynamics has only been found for pacing cycle lengths
longer than 7 s should be seen as further proof that the
cardiorespiratory system has at least one important degree of
freedom that is not entrained by the respiration for
shorter cycle lengths. Therefore, this additional degree of freedom
(oscillator) must have a cycle length distinctly longer than 7 s. A
similar transition at the same cycle length has been found using
symbolic dynamics (24).
The present analysis leads to the hypothesis that this intrinsic degree
of freedom is identical to the so-called low-frequency component or
Mayer wave of the cardiorespiratory system, which is a well-known
oscillatory phenomenon with a cycle length of ~10 s (8) that can
preferentially be identified in blood pressure data.
There is no doubt that the pattern of the RSA is influenced by dynamic
properties of the baroreflex feedback loop and of the intracentral
cardiorespiratory feedback mechanisms. The pronounced interindividual
variation of the one-dimensional map at fixed cycle length and of the
respiratory frequency-dependent transition to it underlines the
expectation that the analysis of the structure of circle maps might be
a good probe for the characterization of the dynamic properties of the
baroreflex and of its interaction with different hypothesized
cardiorespiratory oscillators.
The one-dimensional nonlinear model is able to reproduce features of
the cardiorespiratory system that cannot be expressed by linear
autoregressive models. Thus we expect the characteristic invariants to
express additional physiological features of the cardiorespiratory
system that are inaccessible using the standard measures. Therefore,
the measures based on the pattern of heartbeat intervals within one
respiratory cycle supplement the standard measures based on amplitudes
in the time or frequency domain. The latter are known to be good
indicators for the strength of the parasympathetic-mediated baroreflex
feedback (8). In addition, the quantitative results from nonlinear time
series analyses may be useful in estimating parameters of mechanistic
models. The detailed physiological interpretation of the newly defined
measures should be explored in further studies.
 |
APPENDIX |
Obtaining characteristic measures from the model.
On the basis of the piecewise linear map with parameters as defined in
Fig. 6, an analytic approximation of the characteristic measures can be
given. Three of the invariants may be defined in terms of the average
number of iterations falling into the different line segments
i. The iteration numbers
ki can be approximated analytically (see
Ref. 28 for further details).
By starting from the assumption that the first injection of a
coordinate happens predominantly at the right end of a line segment, an
upper bound for the number of iterations can be derived. First, one
needs the relationship between two successive coordinates (angles)
yk,
yk
1
where
tan
= (yk
yk
1)/(xk
xk
1)
and equals the slope of the line segment under consideration and
y1 is the first
iteration point. By summing up the differences, the
following relationship between the initial coordinate
(y0) and the final
coordinate (yk) can be
obtained
With
the assumption that
yk coincides
with the left end of the line segment, this equation can be interpreted
as an implicit relationship for k.
Solving for k leads to
Using
tan
= ai this can easily be simplified
to the final equation
where
ki is the
average number of iterations in line segment
i,
hi is the
vertical distance from corner point i
to the bisector, and
ai is the slope
of line segment i. Note that
ki is also well defined for noninteger
values.
This relationship is an approximation, because it does not take into
account the exact distribution of the angles, which is not uniform as
can be seen in the empirical data.
By construction,
a3 and
a6 are equal to
1. The corresponding
ki values depend
on the length of line i and the
distance to the bisector:
k3 = (x4
x3)/h3
and k6 = (x1 + 2
x6)/h6.
W is the sum of the six
ki values
By neglecting the difference between the position of the maximum,
x2, and its
theoretical value
= 1/4 and the difference x5
5/4
for the minimum, the following approximations for the asymmetry indexes
can be derived
The degrees of clustering can be expressed in terms of distances to the
bisector
 |
ACKNOWLEDGEMENTS |
The authors thank R. Engbert (Potsdam, Germany), A. Hilgers (London, UK), N. Stollenwerk (Jülich, Germany), and the
members of the Interdisciplinary Hellersen Meetings, initiated by H. P. Koepchen, for valuable comments.
 |
FOOTNOTES |
Present address of K. Suder: Institut für Physiologie, Abteilung
Neurophysiologie, Ruhr-Universität, 44780 Bochum, Germany (E-mail: suder{at}neurop2.ruhr-uni-bochum.de).
Address for reprint requests: F. Drepper, ZEL, Forschungszentrum
Jülich, 52425 Jülich, Germany.
Received 27 March 1997; accepted in final form 18 May 1998.
 |
REFERENCES |
1.
Abel, F. L.,
and
J. A. Waldhausen.
Respiratory and cardiac effects on venous return.
Am. Heart J.
78:
266-275,
1969[Medline].
2.
Abel, H.-H.,
C. Schäfer,
M. Schiek,
and
F. R. Drepper.
Methodik der Analyse der Asymmetrie in den Fluktuationen der Herzfrequenz.
Wien Med. Wochenschr.
145:
474-475,
1995[Medline].
3.
Akselrod, S.,
D. Gordon,
F. A. Ubel,
D. C. Shannon,
A. C. Barger,
and
R. J. Cohen.
Power spectrum analysis of heart rate fluctuations: a quantitative probe of beat-to-beat cardiovascular control.
Science
213:
220-222,
1981[Abstract/Free Full Text].
4.
Angelone, A.,
and
N. A. Coulter.
Respiratory sinus arrhythmia: a frequency-dependent phenomenon.
J. Appl. Physiol.
19:
479-482,
1964[Abstract/Free Full Text].
5.
Chon, K. H.,
T. J. Mullen,
and
R. J. Cohen.
Dual-input nonlinear system analysis of autonomic modulation of heart rate.
IEEE Trans. Biomed. Eng.
43:
530-544,
1996[Medline].
6.
DeBoer, R. W.,
J. M. Karemaker,
and
J. Strackee.
Hemodynamic fluctuations and baroreflex sensitivity in humans: a beat-to-beat model.
Am. J. Physiol.
253 (Heart Circ. Physiol. 22):
H680-H689,
1987[Abstract/Free Full Text].
7.
Drepper, F. R.,
R. Engbert,
and
N. Stollenwerk.
Nonlinear time series analysis of empirical population dynamics.
Ecol. Modell.
75/76:
171-181,
1994.
8.
Eckberg, D. L.
Human sinus arrhythmia as an index of vagal cardiac outflow.
J. Appl. Physiol.
54:
961-966,
1983[Abstract/Free Full Text].
9.
Eckberg, D. L.,
Y. T. Kifle,
and
V. L. Roberts.
Phase relationship between normal human respiration and baroreflex responsiveness.
J. Physiol. (Lond.)
304:
489-502,
1980[Abstract/Free Full Text].
10.
Gilbey, M. P.,
D. Jordan,
D. W. Richter,
and
K. M. Spyer.
Synaptic mechanisms involved in the inspiratory modulation of vagal cardio-inhibitory neurons in the cat.
J. Physiol. (Lond.)
365:
65-78,
1984.
11.
Glass, L.,
and
M. C. Mackey.
From Clocks to Chaos. Princeton, NJ: Princeton University Press, 1988.
12.
Grassberger, P.,
T. Schreiber,
and
C. Schaffrath.
Nonlinear time sequence analysis.
Int. J. Bifur. Chaos
1:
521-547,
1991.
13.
Guevara, M. R.,
and
L. Glass.
Phase locking, period doubling bifurcations and chaos in a mathematical model of a periodically driven oscillator: a theory for the entrainement of biological oscillators and the generation of cardiac dysrhythmias.
J. Math. Biol.
14:
1-23,
1982[Medline].
14.
Herzel, H.,
H. Seidel,
and
H. Warzel.
Heart rate, respiration, and baroreflex: entrainement, bifurcation, and chaos.
Wiss. Z Humboldt-Univ. Berl. Rep. Med.
41:
51-57,
1992.
15.
Hirsch, J. A,
and
B. Bishop.
Respiratory sinus arrhythmia in humans: how breathing pattern modulates heart rate.
Am. J. Physiol.
241 (Heart Circ. Physiol. 10):
H620-H629,
1981[Abstract/Free Full Text].
16.
Kanters, J. K.,
M. V. Hojgaard,
E. Agner,
and
N. H. Holstein-Rathlou.
Influence of forced respiration on nonlinear dynamics of heart rate variability.
Am. J. Physiol.
272 (Regulatory Integrative Comp. Physiol. 41):
R1149-R1154,
1997[Abstract/Free Full Text].
17.
Koepchen, H. P.
History of studies and concepts of blood pressure waves.
In: Mechanisms of Blood Pressure Waves, edited by K. Miyakawa,
H. P. Koepchen,
and C. Polosa. Berlin: Springer, 1984, p. 3-23.
18.
Ludwig, C.
Beiträge zur Kenntniss des Einflusses der Respirationsbewegung auf den Blutlauf im Aortensystem.
Arch. Anat. Physiol.
13:
242-302,
1847.
19.
Melcher, A.
Respiratory sinus arrhythmia in man
a study in heart rate regulation mechanisms.
Acta Physiol. Scand. Suppl.
435:
1-31,
1976[Medline].
20.
Novak, V.,
P. Novak,
J. de Champlain,
R. le Blanc,
R. Martin,
and
R. Nadeau.
Influence of respiration on heart rate and blood pressure fluctuations.
J. Appl. Physiol.
74:
617-626,
1993[Abstract/Free Full Text].
21.
Sayers, B. M.
Signal analysis of heart-rate variability.
In: The Study of Heart-Rate Variability, edited by R. I. Kitney,
and U. Rompelman. Oxford, UK: Clarenden, 1980, p. 27-58.
22.
Schiek, M.
Quantifizierung und Modellierung der respiratorischen Sinusarrhythmie.
In: Ber. Forschungszentrum Jülich, 1994, vol. 2899.
23.
Schiek, M.,
F. R. Drepper,
and
H.-H. Abel.
Mathematisches Modell der respiratorischen Sinusarrhythmie.
Wien Med. Wochenschr.
145:
492-494,
1995[Medline].
24.
Schiek, M.,
F. R. Drepper,
R. Engbert,
H.-H. Abel,
and
K. Suder.
Cardiorespiratory synchronization.
In: Nonlinear Analysis of Physiological Data, edited by K. Kantz,
J. Kurths,
and G. Mayer-Kress. Berlin: Springer, 1998.
25.
Schreiber, T.,
and
A. Schmitz.
Improved surrogate data for nonlinearity tests (Abstract).
Phys. Rev. Lett.
77:
635-638,
1996.[Medline]
26.
Seidel, H.,
and
H. Herzel.
Modeling heart rate variability due to respiration and baroreflex.
In: Modeling the Dynamics of Biological Systems, edited by E. Mosekilde,
and O. G. Mouritsen. Berlin: Springer, 1995, p. 205-229.
27.
Skinner, J. E.,
C. M. Pratt,
and
T. Vybrial.
A reduction in the correlation dimension of heartbeat intervals precedes imminent ventricular fibrillation in human subjects.
Am. Heart J.
125:
731-743,
1993[Medline].
28.
Suder, K.
Nichtlineare Analyse eines Modells der respiratorischen Sinusarrhythmie.
In: Ber. Forschungszentrum Jülich, 1996, vol. 3213.
29.
Theiler, J.,
S. Eubank,
A. Longtin,
B. Galdrikian,
and
J. D. Farmer.
Testing for nonlinearity in time series: the method of surrogate data.
Physica D
58:
77-94,
1992.
30.
Wackerbauer, R.,
A. Witt,
H. Atmanspacher,
J. Kurths,
and
H. Scheingraber.
A comparative classification of complexity measures.
Chaos Solitons Fractals
4:
133-173,
1994.
Am J Physiol Heart Circ Physiol 275(3):H1092-H1102
0002-9513/98 $5.00
Copyright © 1998 the American Physiological Society