International Journal of Bifurcation and Chaos, Vol. 13, No. 3 (2003) 743–754
c World Scientiﬁc Publishing Company
OBSERVATION OF A PERIOD DOUBLING
BIFURCATION DURING ONSET OF HUMAN
VENTRICULAR FIBRILLATION
MICHAEL SMALL∗
Department of Electronic and Information Engineering,
Hong Kong Polytechnic University, Kowloon, Hong Kong
ensmall@polyu.edu.hk
DEJIN YU and ROBERT G. HARRISON
Department of Physics, Heriot-Watt University, Edinburgh, UK
Received August 30, 2001; Revised March 20, 2002
We ﬁnd evidence of a period doubling bifurcation and transition from pseudo-periodic chaos to
a stable limit cycle during onset of ventricular ﬁbrillation (VF) in humans. Novel radial basis
modeling techniques are applied to time series of spontaneous VF to estimate time dependent
changes in the dynamics. Furthermore, we show that the spectral power content may be used
as a surrogate for the underlying bifurcation parameter. With this spectral power measure we
demonstrate a characteristic transition from pseudo-periodic chaos to a second pseudo-periodic
chaotic regime. The methods described in this paper are generic and applicable to any time
series data.
Keywords: Ventricular ﬁbrillation; period doubling; radial basis modeling.
1. Introduction
Ventricular ﬁbrillation (VF) is a rapidly lethal cardiac
arrhythmia and a common cause of death in
the industrialized world. However, the mechanism
underlying onset of VF and successful treatment via
electrical deﬁbrillation is poorly understood.
Experimental observations in humans and
animals have shown that the cardiac system is inherently
nonlinear and that signiﬁcant nonlinear deterministic
structure can be extracted from ECG
time series recordings [Small et al., 2000c; Small
et al., 2000b; Ravelli & Antolini, 1992]. Experimental
preparations of cardiac cells have been shown
to exhibit phase locking, period doubling and onset
of chaos [Guevara et al., 1981; Chialvo & Jalife,
1990; Savino et al., 1989; Karagueuzian et al., 1993].
However, these experiments consist of externally
driven isolated cellular preparations. Theoretical
models of the atrio-ventricular (AV) junction and
periodically forced cardiac cells predict phase locking
and period doubling bifurcation [West et al.,
1985]. Period doubling, from period 1 to period 2
(known in the medical literature as alternans), can
be occasionally observed in the diseased heart (see,
e.g. [Goldberger et al., 1985]) and has recently been
shown to precede ventricular ﬁbrillation in theoretical
models [Hastings et al., 2000]. Cohen and
colleagues have demonstrated period doubling (up
to ﬁve times the base period) in the canine heart
[Ritzenberg et al., 1984] and have more recently
shown that T-wave alternans is a signiﬁcant predictor
of arrhythmia [Klingenheben et al., 2000].
∗
Author for correspondence.
743
744 M. Small et al.
Although alternans is indicative of a period
doubling bifurcation, a period doubling bifurcation
route to chaos may not necessarily manifest itself
as alternans. Alternans is deﬁned as a doubling of
the underlying period of a time series (see Fig. 3),
a period doubling bifurcation occurs when the order
of the periodic behavior doubles — typically to
two periods very close to one another, i.e. Alternans
is suﬃcient but not necessary. Many authors have
speculated that a period doubling bifurcation may
provide a route to chaos during ventricular arrhythmia
in humans [Goldberger et al., 1988; Hoekstra
et al., 1997; Zhang et al., 1999]. However, no conclusive
evidence has been provided. In this paper
we infer the existence of a period doubling bifurcation
route to chaos prior to the onset of ventricular
tachyarrhythmia in humans through the analysis of
nonlinear models of clinical ECG data. Whereas
alternans is an easily observable and measurable
phenomenon (Fig. 3) we apply nonlinear modeling
techniques to search for more subtle features, not
visually apparent in the original data. A period
doubling bifurcation is not apparent from a cursory
examination of the data in Fig. 1. We deduce the
existence of such a bifurcation in a model of this
data.
To search for time dependent nonlinear dynamic
structure we employ powerful nonlinear
radial basis modeling techniques described in [Judd
& Mees, 1995] and generalized in [Small & Judd,
1998]. These modeling methods have the necessary
nonlinearity to accurately describe cardiac
dynamics (unlike linear or polynomial models)
and are amenable to analysis with the techniques
of nonlinear dynamical systems theory (unlike
neural networks or local linear methods).
Furthermore, these models are generic and the
modeling algorithm is (in general) assumption free.
By employing an information theoretic model selection
criteria this algorithm chooses the simplest
model that is consistent with the observed dynamics.
To accurately model bifurcation in cardiac
dynamics these existing methods must be modiﬁed.
By embedding time as a state variable [Small
et al., 2001] we observe a period doubling bifurcation
and chaotic dynamics during sinus rhythm and
ventricular tachycardia (VT) immediately preceding
VF. However, we ﬁnd that successive runs of the
modeling algorithm do not produce quantitatively
identical results. Therefore, we explicitly calculate
an index which we use as a surrogate for the true
(hidden) bifurcation parameter and embed this as a
state variable. This technique produces characteristic
and reproduceable results. We observe a transition
from pseudo-periodic chaos (sinus rhythm prior
to VT and VF) to a stable limit cycle (during VT)
and to (distinct) pseudo-periodic chaos (VF).
By pseudo-periodic chaos we mean the commonly
observed type of chaos in systems such as
the R¨ossler equations. A periodic orbit becomes
bi-periodic and cascades through successive period
doubling bifurcations; in the Poincar´e section; to
chaos. In the chaotic regime the system exhibits
chaotic, approximately cyclic behavior. The modeling
described in this paper characterizes the transition
from the physiological states of sinus rhythm
to VT and VF as a period doubling bifurcation.
We must emphasize that applying this modeling
algorithm to experimental time series data
cannot prove the existence of a bifurcation to chaos
in the dynamical system underlying the observed
data. However, this information theoretic modeling
algorithm estimates the simplest dynamical system
(model) that is consistent with the observed time series.
We may then apply computational dynamical
systems theory to estimate the bifurcation behavior
of that model. From this we are able to show that
this observed behavior is both consistent with the
time series and also the simplest consistent behavior;
from the class of models employed.
1.1. Data
We examine four recordings of spontaneous initiation
and evolution of VF in human subjects. Such
time series data are limited and diﬃcult to obtain.
Using a new data collection facility [Small
et al., 2000a] we have been able to record time series
showing spontaneous evolution from sinus rhythm
to VF, and subsequent treatment. In two of these
recordings VF is preceded by VT. One of the other
two recordings shows a direct transition from sinus
rhythm to VF. The fourth recording shows
bradycardic (slow rhythm) behavior prior to onset
of VF. For each of these four recordings we
selected an approximately 90 second (45000 data
point) sub-sequence covering the transition from sinus
rhythm, through any intermediate rhythms to
VF. We down-sampled the time series by taking
a four point moving average and sampling every
fourth point, yielding 90 second episodes of approximately
11250 data points. Down sampling was done
Period Doubling Bifurcation During Onset of Human Ventricular Fibrillation 745
22 24 26 28 30 32 34 36
2.5
3
3.5
ECG(mV)
38 40 42 44 46 48 50 52
2.5
3
3.5
ECG(mV)
52 54 56 58 60 62 64 66
2.5
3
3.5
ECG(mV)
Time (Seconds)
Fig. 1. Spontaneous evolution of VF in a human. The data shown here has been down-sampled to 125 Hz. On each
plot the horizontal axis is time in seconds and the vertical axis is surface (ECG) electrical potential in milli-volts (mV). The
top plot shows a short section of sinus rhythm, including ectopic beats. The second plot shows initiation and evolution of
VT and the third shows VF. These plots are contiguous and the horizontal axes are equal. The entire data set is shown in
Fig. 4(a).
0 5 10 15
1
2
3
ECG(mV)
15 20 25 30
1
2
3
ECG(mV)
30 35 40 45
1
2
3
ECG(mV)
Time (Seconds)
Fig. 2. Spontaneous evolution of VT and bigeminy in a human. The data shown here has been down-sampled to
125 Hz. On each plot the horizontal axis is time in seconds and the vertical axis is surface (ECG) electrical potential in
milli-volts (mV). The top plot shows a short section of sinus rhythm and spontaneous evolution of VT. The second plot shows
self termination of VT, following this regular QRS complexes (normal heart beats) and ectopic (abnormal) beats alternate.
In the third panel sinus rhythm is restored. These plots are contiguous and the horizontal axes are equal.
746 M. Small et al.
55 60 65 70
1
1.5
2ECG(mV)
70 75 80 85
1
1.5
2
ECG(mV)
85 90 95 100
1
1.5
2
ECG(mV)
Time (Seconds)
Fig. 3. Period doubling in human ECG. The data shown here has been down-sampled to 125 Hz. On each plot the
horizontal axis is time in seconds and the vertical axis is surface (ECG) electrical potential in milli-volts (mV). This data
shows a spontaneous period doubling in the respiratory rate of a human subject (at around 76 seconds). Prior to this several
ectopic (premature) beats are observed (for example, at 57 and 71 seconds). These plots are contiguous and the horizontal
axes are equal.
to reduce the computational load to a manageable
level. By taking a four point moving average before
down sampling we reduced the amount of temporal
information in the time series but increased the
spatial resolution. After down sampling, the ﬁrst
11000 data points were used to construct a cylindrical
basis model.
Figure 1 shows one of the recordings used in
this study. For comparison we have also applied this
analysis to recordings of other spontaneous evolving
physiological rhythms. Figure 2 shows a recording
of bigeminy (a phenomenon where normal and
abnormal beats alternate) and Fig. 3 shows spontaneous
alternans.
1.2. Overview
We apply nonlinear cylindrical basis modeling
techniques to these data in two slightly diﬀerent
incarnations. In Sec. 2 we describe the nonlinear
cylindrical basis modeling methodology we employ.
Section 3 describes the application of this method
with no further assumptions about the underlying
bifurcation. In Sec. 4 we assume that visually apparent
structure in the time series is indicative of
the underlying bifurcation and estimate a scalar
time dependent parameter as a ﬁrst approximation
to the bifurcation parameter.
2. Cylindrical Basis Models and
Minimum Description Length
Cylindrical basis models are a generalization of the
well known radial basis models [Powell, 1992]. Let
{yt}N
t=1 be a scalar time series on N observations.
For an embedding window dw, a cylindrical basis
models is a function F : Rdw
−→ R such that
F(zt) = a0 +
k
i=1
aiyt− i
+
m
j=1
λjφ
Pj(zt − cj)
rj
, (1)
Period Doubling Bifurcation During Onset of Human Ventricular Fibrillation 747
where zt = (yt−1, yt−2, . . . , yt−dw ) ∈ Rdw
, ai, λj ∈
R are the weights, cj ∈ Rdw
are the centers, and
rj ∈ R+ are the radii. The lags i satisfy 0 < i <
i+1 ≤ dw and the functions Pj(·) are orthogonal
projections from Rdw
to some subset of the coordinate
axis. The functions φ are typically some class
of C2 functions satisfying ∞
0 φ(x)dx < ∞. For the
computations in this paper we use a combination of
Gaussian basis functions and Morlet wavelets.
The essential diﬀerence between cylindrical basis
models and standard radial basis models is the
inclusion of the projection functions Pj(·). These
functions allow for the projection onto diﬀerent, signiﬁcant,
subspaces of Rdw
in diﬀerent parts of phase
space. This is both useful and intuitive since the
complexity of most nonlinear dynamical systems
varies with location in phase space. For example,
the Lorenz attractor is basically two dimensional on
the “wings”, but the central separatrix contains a
three-dimensional structure [Judd & Mees, 1998].
The function F is used to approximate the
evolution operator of a dynamical system by
F(zt) = F(yt−1, yt−2, . . . , yt−dw ) = yt . (2)
With the model size (k, m) ﬁxed, our modeling
algorithm selects ai, i, λj, cj, rj and Pj(·) such
that the model prediction error
Error = e2
t =
N
t=dw+1
(F(zt) − yt)2
, (3)
is minimized. The optimal model size is then selected
by ﬁnding the values of k and m such that
the model description length is minimized. The description
length of a model is roughly the amount
of compression of the original data that is achieved
by describing the model and model prediction errors,
over specifying the original data values [Judd
& Mees, 1995; Small & Judd, 1998; Small et al.,
2001]. The idea of minimum description length was
originally described in [Rissanen, 1989].
For a given model (1) with error speciﬁed by
(3) the description length is approximately
N − dw
2
1 + ln
2πeT e
N − dw
+
1
2
+ ln γ (k + m) −
k+m
j=1
ln δj (4)
where γ is the constant number of (binary) orders
of magnitude required to cover the range of the data
and δj are the precision with which one must specify
each of the (k + m) model parameters (the ai and
λj terms in (1)). The δj terms are computed as
the solution of (Qδ)j = 1/δj where Q is the second
derivative of (1) with respect to the model parameters
Ai and λj.
3. Assumption Free Modeling
The model described by Eqs. (1) and (2) is time
independent. The process we are modeling exhibits
time dependent features and these equations are unlikely
to provide an adequate model. One possible
solution is to model each of the parameters ai, i, λj,
cj, rj and Pj(·) as a function of time. However, this
approach is somewhat impractical. A more feasible
approach is to extend the embedding
zt = (yt−1, yt−2, . . . , yt−dw ) (5)
to
ˆzt = (yt−1, yt−2, . . . , yt−dw , ξ(t)) (6)
so that time t is explicitly included as a dependent
variable of F(ˆzt) [Small et al., 2001]. The aﬃne
transformation ξ is included so that the time index
“stretches” the attractor. If ξ is too large, then
the embedded points ˆzt will be too sparse in phase
space, if ξ is too small then there will eﬀectively be
insuﬃcient separation. For this data we have found
that it is suﬃcient to choose ξ so that the maximum
and minimum values of ξ(t) are the same as
the maximum and minimum of yt. In general, a
more robust choice may be to deﬁne ξ in terms of
the standard deviation of yt.
The modeling algorithm we apply in this paper
must minimize a nonlinear function (description
length) of several variables (ai, i, λj, cj, rj
and Pj(·)). To solve this problem completely is beyond
our limited computational resources. For this
reason the algorithm we have implemented involves
a stochastic search for a local minimum of description
length. For many data sets this is suﬃcient to
obtain reproducible and consistent results [Judd &
Mees, 1998, 1996; Judd & Small, 2000; Small et al.,
1999a; Small et al., 1999b]. We have tested this
modeling algorithm with four systems with known
bifurcation structure and obtained good agreement
between actual and predicted results [Small
et al., 2001]. This modeling methodology is able
748 M. Small et al.
to accurately uncover bifurcation structure from
time series simulations of (i) the R¨ossler equations,
(ii) FitzHugh–Nagumo type simulations of ventricular
arrhythmia, and (iii) nonstationary noise
processes as well as from (iv) experimental data of
infant respiration [Small et al., 2001]. The computational
simulations examined in [Small et al., 2001]
have well known and easily veriﬁable dynamics, this
provides a test of the accuracy of this approach. We
were able to correctly extract a nonstationary trend
from an i.i.d. noise simulation [Small et al., 2001]
with a slight nonstationarity. When applied to a
noisy simulation of the R¨ossler system undergoing
a period doubling bifurcation the modeling algorithm
was able to correctly identify the qualitative
features of the bifurcation from the noisy data (5000
observations) [Small et al., 2001].
The model we utilize here is an ordinary
iterative map with a tunable parameter. This bifurcation
parameter has been allowed to vary in the
construction of the model. By keeping this parameter
ﬁxed and observing the iterated behavior of the
model we have an estimate of the asymptotic dynamics
of the system for that (ﬁxed) value of the
bifurcation parameter. For each such trajectory we
estimate the (asymptotic) peak and trough values
observed. For a periodic orbit these will be ﬁxed;
for period 2 dynamics there will be two distinct amplitudes;
and for chaos the sequence of amplitudes
will be irregular. It is these estimates of the asymptotic
amplitudes that are plotted in Figs. 4 and 6 as
a function of the bifurcation parameter. Representative
trajectories for various values of the bifurcation
parameter are depicted in Fig. 5. For example,
1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 11000
2.5
3
3.5
ECG(mV)
1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 11000
2.6
2.8
3
Peak&Trough
1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 11000
0.2
0.4
0.6
Time index (125 Hz)/Bif. param.
Est.Amplitude
Fig. 4. Period doubling chaos prior to onset of VF. The top panel shows the original time series. The second panel
shows the estimated asymptotic peak (red) and trough (green) values for ﬁxed values of the bifurcation parameter. The
bottom plot shows the estimated asymptotic amplitudes (blue). The horizontal axis is time series datum number, the model
bifurcation parameter is an aﬃne transformation of this. Sinus rhythm is characterized as a periodic orbit, prior to initiation
of VT this periodic orbit undergoes a period doubling bifurcation and becomes chaotic. At onset of VT the chaotic behavior
bifurcates to pseudo-periodic chaos and period 8 dynamics. At the transition to VF this behavior changes to a stable limit
cycle and eventually a stable focus.
Period Doubling Bifurcation During Onset of Human Ventricular Fibrillation 749
100 200 300 400 500 600 700 800 900 1000
datum1000
100 200 300 400 500 600 700 800 900 1000
datum2000
100 200 300 400 500 600 700 800 900 1000
datum4000
100 200 300 400 500 600 700 800 900 1000
datum5000
100 200 300 400 500 600 700 800 900 1000
datum6000
100 200 300 400 500 600 700 800 900 1000
datum8000
Fig. 5. Trajectories estimates from a model. For ﬁxed values of the bifurcation parameter, the system is run for 1000
iterations. The dynamic behavior for various values of the bifurcation parameter are shown (from top): periodic, period 2,
period 3, chaos, period 4 and a ﬁxed point. The chosen ﬁxed values of bifurcation parameter correspond to datum 1000, 2000,
4000, 5000, 6000 and 8000 of Fig. 4 (respectively). Trajectories such as these are used to estimate the amplitudes shown in
Figs. 4 and 6. Note that while these trajectories do not capture the exact quantitative features of the system (for example,
the periodic behavior observed at datum 1000 is unlike sinus rhythm ECG) these features are quantitatively appropriate.
from the representative trajectories shown in Fig. 5
we would deduce (from top to bottom): a periodic
orbit with only one observable amplitude, period 2,
period 3, chaos, period 4, and a ﬁxed point (zero
amplitude).
In Figs. 4 and 6 we provide an estimate, derived
from a single application of our modeling scheme to
the data set of Fig. 1. This estimate is representative
of our other results. The dynamics presented
in this ﬁgure contain all the various complex behaviors
observed in separate runs of the modeling
routine. Some iterations of the modeling routine did
not produce the complete period doubling bifurcation
illustrated here. However, they all exhibited
a complex bifurcation from pseudo-periodic regime
(corresponding to sinus rhythm) to chaos (VT) and
a stable periodic orbit or focus (VF). Models from
each of the four recordings of onset of VF are qualitatively
similar.
We compared this technique and the bifurcation
behavior estimated from models with traditional
frequency domain methods applied directly
750 M. Small et al.
Time (Seconds)
Est.Amplitude
0 10 20 30 40 50 60 70 80 90
0
0.05
0.1
0.15
0.2
0.25
0.3
Fig. 6. Period doubling chaos prior to onset of VF. The calculation of Fig. 4 is shown here as a probability density
plot. The ﬁgure depicts the estimated distribution of asymptotic oscillation amplitudes as a function of time (in seconds).
High likelihood is depicted in bright red and low likelihood in dark blue.
Time (Seconds)
Frequency(Hz)
0 10 20 30 40 50 60 70 80 90
0
10
20
30
40
50
60
Time (Seconds)
Period(log2
(samples))
0 10 20 30 40 50 60 70 80 90
2
4
6
8
10
12
(a)
(b)
Fig. 7. Period doubling chaos prior to onset of VF. (a) Fourier spectrogram (256 FFT with an overlap of 128 points)
and (b) wavelet analysis of the data used in Figs. 4 and 6. In contrast to the analysis in Figs. 4 and 6 these calculations provide
little additional information compared to the original time series. Most of the change in frequency content that is observable
from these plots is also observable directly from the original time series. High density is shown in red and low density in blue.
Period Doubling Bifurcation During Onset of Human Ventricular Fibrillation 751
1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 11000
0.2
0.4
0.6
0.8
Bif.param.(arb.units)
1000 2000 3000 4000 5000 6000 7000 8000 9000 10000
2.6
2.8
3
Peak&Trough
1000 2000 3000 4000 5000 6000 7000 8000 9000 10000
0.1
0.2
0.3
0.4
Time index (125 Hz)/Bif. param.
Est.Amplitude
Fig. 8. Bifurcation diagram estimated with surrogate bifurcation parameter. The top plot shows values used for
the bifurcation parameter as a function of the time index (sampled at 125 Hz) for the same data set as Fig. 1. The second
plot shows the asymptotic peak (red) and trough (green) values as a function of the time index and the third plot shows
the estimated amplitude behavior (blue). The complex period doubling bifurcation evident in Fig. 4(b) is absent, but a clear
transition from large scale R¨ossler type chaos (for sinus rhythm preceding VF) through period 6 behavior (VT) to another
pseudo-periodic chaotic regime (VF), is apparent.
to the data. Figure 6 is an estimate of the probability
distribution of the amplitude of system orbits
for various values of the bifurcation parameter (time
index). The data is the same as the lower panel
of Fig. 4, expressed as a probability distribution.
For comparison we computed a windowed spectrogram
and wavelet transform from the original data
(Fig. 7). While the coarse time dependent features
are evident in these plots, we ﬁnd very little information
not evident in the original time series.
The asymptotic time behavior displayed in Fig. 6 is
achieved by extrapolating the model dynamics for
a model ﬁtted to the data.
Qualitative features of period doubling and onset
of chaos are present in each of the models despite
this complex structure not being obvious in the original
time series. The modeling algorithm is ﬁtting
subtle features of the asymptotic dynamics that
may not be observed directly during the comparatively
rapid change in the system dynamics prior to
(and during) onset of ventricular arrhythmia.
Although essential features are consistent between
most models, this variation between models
is somewhat unsatisfactory. The problem is that in
addition to estimating unknown dynamics we are
relying on the modeling algorithm to estimate an
unknown bifurcation parameter at each point in
time. To overcome this problem we estimated the
bifurcation parameter a priori and built a model
based on this bifurcation parameter.
4. A priori Bifurcation Parameter
Estimation
Generically, the time series recordings of evolution
of VF exhibit two or three of four distinct dynamical
752 M. Small et al.
regimes: sinus rhythm, VT, bradycardia and VF.
A bifurcation parameter must reﬂect this. For the
data presented in this paper, we found that the
proportion of the power spectral content between
4 and 8 Hz performed well. Alternative indices
including entropy, complexity and point-wise correlation
dimension were also considered. However, we
found that spectral power content provided the best
separation between diﬀerent physiological rhythms.
In our experience spectral power content also performed
best for real time identiﬁcation of arrhythmia
[Small et al., 2000a]. Spectral power captures
many of the essential features of the time series
(Fig. 7). However, this spectral variation alone is
insuﬃcient to deduce the period doubling bifurcation
and onset of chaos observed in models of this
data. Figure 8 shows estimates of this quantity for
the same time series as in Fig. 4.
We estimated the power content between 4
and 8 Hz for a 2048 point sliding window (before
down sampling the original recordings) and used
this quantity as a surrogate for the true, unknown,
bifurcation parameter. That is, the embedding (5)
or (6) is replaced by
˜zt = (yt−1, yt−2, . . . , yt−dw , ξ(p(yt))) (7)
where p(yt) is the above mentioned statistic, and
ξ(·) is an aﬃne scale transformation. The remainder
of the modeling process is identical.
Figure 8 depicts a representative result of this
calculation for the data illustrated in Fig. 4. Repeated
runs of the modeling algorithm on the
same data produced equivalent results. The results
between data sets were similar. Whereas the
calculation estimating both model and bifurcation
parameter would typically yield smooth bifurcation
from a periodic orbit, to chaos (sinus rhythm prior
to VT and during VT) and back to a noisy periodic
orbit (VF) this is not the case with an a priori
estimate of the bifurcation parameter. For smooth
and gradual change in the values of this bifurcation
parameter the change in the underlying dynamics is
similarly subtle. However, as this parameter is estimated
directly from the time series the value often
changes fairly rapidly (Fig. 8), and this in turn produces
sudden changes in the underlying dynamics.
Closer examination of small changes in the bifurcation
parameter do show a period doubling bifurcation
during the transition between sinus rhythm and
tachycardia. The weakness of this approach is that
we are constraining the time dependent dynamics
of the model to conform to our expectation of the
system, i.e. ξ(p(yt)).
5. Conclusion
In this paper we have described two alternative
techniques to examine the changing dynamics in experimental
time series. Although we consider only
human ECG recordings the method itself is generic
and equally applicable to any experimental time series.
Estimating the model and bifurcation parameter
simultaneously (Sec. 3) provides a method of
examining the changing dynamics with no a priori
information. For short, noisy, or experimental data
this method has its limitations. Primarily, the computational
burden to estimate the optimal model,
model parameters, and bifurcation parameter from
a single scalar time series is immense. This problem
meant that it was diﬃcult to get fully repeatable results
(for a single time series). Rather than obtain
the global minimum of (4) the modeling algorithm
would halt at a locally optimal value.
To reduce the computational eﬀort we provided
an alternative approach — estimating the model
bifurcation parameter a priori (Sec. 4). This approach
meant that the models produced were more
easily repeatable (for a single time series) and also
consistent (between diﬀerent time series). The main
problem with this approach is ﬁnding a “reasonable”
surrogate for the bifurcation parameter, and
the restriction of the dynamics that this implies.
Estimating the bifurcation parameter and
model simultaneously meant that the problem of
reproducing the dynamics was reduced to ﬁnding
the right nonlinear function from (an aﬃne transformation
of) the time index to the underlying bifurcation
parameter. By assuming an appropriate
value of bifurcation parameter ﬁrst, the existence of
an injective mapping from the surrogate bifurcation
parameter ξ(p(yt)) to the underlying bifurcation parameter
is not guaranteed (i.e. p(yt) is not necessarily
a 1-to-1 function of t). However, we hope that a
“reasonable” choice of the surrogate parameter will
provide a trivial transformation from ξ(p(yt)) to the
underlying bifurcation parameter.
When applied to recordings of initiation of VF
these modeling techniques provide a description of
the simplest dynamical system (within the model
class) consistent with the observed data. Observations
derived from these models are possible dynamical
phenomena consistent with the data. This does
Period Doubling Bifurcation During Onset of Human Ventricular Fibrillation 753
not imply that the data must have been produced
by such phenomena, only that such phenomena are
a simple and compact explanation for the observed
data. Prior to onset of ventricular arrhythmia, sinus
rhythm may behave as a noisy dynamical system
exhibiting R¨ossler type chaos (i.e. chaotic pseudoperiodic
orbits). We also ﬁnd evidence that a period
doubling bifurcation may precede imminent ventricular
tachyarrhythmia. In the arrhythmic state this
dynamical system behaves as either R¨ossler type
chaos (Sec. 3) or a stable focus (Sec. 4). We believe
that this possible contradiction is due to limitations
of the modeling process. When modeling
both dynamics and estimating bifurcation parameter
(Sec. 3) the modeling algorithm is unable to
extract suﬃcient information from the VF and VT
states to model them as more than noisy periodic
orbits. When the bifurcation parameter is assumed
a priori (Sec. 4), and p(yt) is a nonmonotonic function
of t (i.e. the value remains the same at diﬀerent
times, leading to a more densely populated phase
space) arrhythmic (VT and VF) behavior can be
better modeled as pseudo-periodic chaos (R¨ossler
type chaos).
These results imply that both sinus rhythm
and VF may be characterized as chaotic dynamical
systems. However the underlying dynamics during
sinus rhythm, VT and VF are all fundamentally
diﬀerent. This result is therefore consistent with
earlier work estimating dynamical invariants from
time series [Small et al., 2000c; Small et al. 2000b;
Ravelli & Antolini, 1992; Small et al., 1999b; Yu
et al., 2000]. The similarity between this result and
various observations of alternans is signiﬁcant but
not conclusive. We observe an increase in the order
of periodic behavior in the amplitude of respiration
whereas alternans typically manifests as a doubling
in the underlying respiratory rate (and sometimes
observed directly as two distinct amplitudes). We
found no evidence for chaotic bifurcations in recordings
of alternans or bigeminy.
We also found that prior to onset of tachyarrhythmia
the pseudo-periodic chaos observed in
models of these time series becomes strictly periodic
and a period doubling bifurcation may be
observed as a mechanism consistent with transition
between these states. These results are preliminary,
and certainly not clinically signiﬁcant. A
much larger sample needs to be considered before
deﬁnitive statements may be made. Nevertheless,
observation of period doubling bifurcations prior to
onset of VF suggest that the human cardiac system
undergoes a fundamental change in its dynamical
state prior to onset. Methods such as Lyapunov
exponent and point-wise correlation dimension [Yu
et al., 2000] estimation may therefore be applied to
predict imminent arrhythmia. However, such algorithms
need to be substantially modiﬁed to be both
sensitive and accurate for extremely short time series
if they are to be applied as clinical indicators
of imminent arrhythmia.
Acknowledgments
This research was funded through a Research Development
Grant by the Scottish Higher Education
Funding Council (SHEFC), No. RDG/078. M.
Small is supported by a Hong Kong Polytechnic
University Postdoctoral Fellowship award (No. GYW55).
We would also like to acknowledge the
assistance of N. Grubb, K. Fox and the staﬀ of
the Coronary Care Unit of the Royal Inﬁrmary of
Edinburgh.
References
Chialvo, D. R., Gilmour, Jr. R. F. & Jalife, J. [1990] “Low
dimensional chaos in cardiac tissue,” Nature 343,
653–657.
Goldberger, A. L., Bhargava, V., West, B. J. & Mandell,
A. J. [1985] “Nonlinear dynamics of the heartbeat. II.
Subharmonic bifurcations of the cardiac interbeat interval
in sinus node disease,” Physica D17, 207–214.
Goldberger, A., Rogney, D., Mietus, J., Antman, E. &
Greenwald, S. [1988] “Nonlinear dynamics in sudden
cardiac death syndrome: Heartrate oscillations and
bifurcations,” Experientia 44, 983–987.
Guevara, M. R., Glass, L. & Shrier, A. [1981] “Phase
locking, period-doubling bifurcations and irregular
dynamics in periodically stimulated cardiac cells,”
Science 214, 1350–1353.
Hastings, H. M., Fenton, F.H., Evans, S.J., Hotomaroglu,
O., Geetha, J., Gittelson, K., Nilson, J. & Garﬁnkel,
A. [2000] “Alternans and the onset of ventricular
ﬁbrillation,” Phys. Rev. 62, 4043–4048.
Hoekstra, B. P., Diks, C. G., Allessie, M. A. &
DeGoode, J. [1997] “Nonlinear analysis of the pharmacological
conversion of sustained atrial ﬁbrillation
in conscious goats by the class Ic drug cibenzoline,”
Chaos 7, 430–446.
Judd, K. & Mees, A. [1995] “On selecting models for
nonlinear time series,” Physica D82, 426–444.
Judd, K. & Mees, A. [1996] “Modeling chaotic motions
of a string from experimental data,” Physica D92,
221–236.
Judd, K. & Mees, A. [1998] “Embedding as a modeling
problem,” Physica D120, 273–286.
754 M. Small et al.
Judd, K. & Small, M. [2000] “Towards long-term prediction,”
Physica D136, 31–44.
Karagueuzian, H., Khan, S., Hong, K., Kobayashi, Y.,
Denton, T., Mandel, W. & Diamond, G. [1993]
“Action-potential alternans and irregular dynamics
in quinidine-intoxicated ventricular muscle-cells —
implications for ventricular proarrhythmia,” Circulation
87, 1661–1672.
Klingenheben, T., Zabel, M., D’Agostino, R., Cohen, R.
& Hohnloser, S. [2000] “Predictive value of t-wave
alternans for arrhythmis events in patients with congestive
heart failure,” Lancet 356, 651–652.
Powell, M. J. D. [1992] “The theory of radial basis
function approximation in 1990,” Advances in Numerical
Analysis. Volume II: Wavelets, Subdivision Algorithms
and Radial Basis Functions, ed. Light, W.
(Oxford Science Publications) Chap. 3, pp. 105–210.
Ravelli, F. & Antolini, R. [1992] “Complex dynamics underlying
the human electrocardiogram,” Biol. Cybern.
67, 57–65.
Rissanen, J. [1989] Stochastic Complexity in Statistical
Inquiry (World Scientiﬁc, Singapore).
Ritzenberg, A. L., Adam, D. R. & Cohen, R. J. [1984]
“Period multupling-evidence for nonlinear behavior of
the canine heart,” Nature 307, 159–161.
Savino, G. V., Romanelli, L., Gonz´alez, D. L., Piro, O.
& Valentinuzzi, M. E. [1989] “Evidence for chaotic behavior
in driven ventricles,” Biophys. J. 56, 273–280.
Small, M. & Judd, K. [1998] “Comparison of new nonlinear
modeling techniques with applications to infant
respiration,” Physica D117, 283–298.
Small, M., Judd, K., Lowe, M. & Stick, S. [1999a] “Is
breathing in infants chaotic? Dimension estimates
for respiratory patterns during quiet sleep,” J. Appl.
Physiol. 86, 359–376.
Small, M., Yu, D., Harrison, R. G., Robertson, C. Clegg,
G., Holzer, M. & Sterz, F. [1999b] “Characterizing
nonlinearity in ventricular ﬁbrillation,” Comput.
Cardiol. 26, 17–20.
Small, M., Yu, D., Grubb, N., Simonotto, J., Fox, K. &
Harrison, R. G. [2000a] “Automatic identiﬁcation and
recording of cardiac arrhythmia,” Comput. Cardiol.
27, 355–358.
Small, M., Yu, D. & Harrison, R. G. [2000b] “Nonlinear
analysis of human ECG rhythm and arrhythmia,”
Comput. Cardiol. 27, 147–150.
Small, M., Yu, D., Harrison, R. G., Robertson, C.,
Clegg, G., Holzer, M. & Sterz, F. [2000c] “Deterministic
nonlinearity in ventricular ﬁbrillation,” Chaos 10,
268–277.
Small, M., Yu, D. & Harrison, R. G. [2001] “Nonstationarity
as an embedding problem,” in Space Time
Chaos: Characterization, Control and Synchronization
(Birkhauser, Boston), pp. 3–18.
West, B. J., Goldberger, A. L., Rovner, G. &
Bhargava, V. [1985] “Nonlinear dynamics of the heartbeat.
I. The AV junction: Passive conduit or active
oscillator?” Physica D17, 198–206.
Yu, D., Small, M., Harrison, R. G., Robertson, C.,
Clegg, G., Holzer, M. & Sterz, F. [2000] “Measuring
temporal complexity of ventricular ﬁbrillation,” Phys.
Lett. A265, 68–75.
Zhang, X.-S., Zhu, Y.-S., Thakor, N. V., Wang,
Z.-M. & Wang, Z.-Z. [1999] “Modeling the relationship
between concurrent epicardial action potentials and
bipolar electrograms,” IEEE Biomed. 46, 365–376.