Topological Constraints On Spiral Wave Dynamics In Spherical ...
PHYSICAL REVIEW E 70, 056203 (2004)
Topological constraints on spiral wave dynamics in spherical geometries
with inhomogeneous excitability
Jörn Davidsen,1,2,* Leon Glass,3,† and Raymond Kapral1,2,‡
1Max-Planck-Institut für Physik Komplexer Systeme, Nöthnitzer Strasse 38, 01187 Dresden, Germany
2Chemical Physics Theory Group, Department of Chemistry, University of Toronto, Toronto, ON M5S 3H6, Canada
3Department of Physiology, McGill University, Montreal, Quebec H3G 1Y6, Canada
(Received 29 August 2003; published 10 November 2004)
We analyze the way topological constraints and inhomogeneity in the excitability influence the dynamics of
spiral waves on spheres and punctured spheres of excitable media. We generalize the definition of an index
such that it characterizes not only each spiral but also each hole in punctured, oriented, compact, two-
dimensional differentiable manifolds and show that the sum of the indices is conserved and zero. We also show
that heterogeneity and geometry are responsible for the formation of various spiral-wave attractors, in particu-
lar pairs of spirals in which one spiral acts as a source and a second as a sink—the latter similar to an antispiral.
The results provide a basis for the analysis of the propagation of waves in heterogeneous excitable media in
physical and biological systems.
DOI: 10.1103/PhysRevE.70.056203
PACS number(s): 89.75. k, 87.19. j, 05.40. a
I. INTRODUCTION
of atrial arrhythmias because the thin walls of the atria can
be described as two-dimensional (2D) inhomogeneous excit-
Geometry and inhomogeneity influence pattern formation
able media with specific geometrical features [7].
in chemical and biological systems [1,2]. One example
where these two factors play a crucial role is in the experi-
mental observations of distinctive spiral-wave dynamics on
II. INDEX THEOREM FOR PHASE SINGULARITIES
the surfaces of spherical beads, which are excitable inhomo-
geneous chemical media [3,4]. A biological example is the
The mathematical description of spiral waves is based on
origin of abnormal cardiac rhythms in the human heart which
the notion of phase which in turn allows one to characterize
depend on the anatomical substrate. The heart possesses a
spiral waves by an index. From this description, a number of
complex nonplanar geometry with multiple chambers, with
topological results placing restrictions on spiral-wave dy-
holes corresponding to valves and blood vessels. Some seri-
namics can be derived [5,10–13].
ous arrhythmias are associated with circulating spiral waves
With the exception of a finite number of singular points,
similar to those observed in chemical media [5]. Since an
with each point in an orientable and compact two-
abnormal anatomical substrate is a common finding in pa-
dimensional differentiable manifold M we identify a unique
tients with some types of cardiac arrhythmias and interven-
phase lying on the unit circle,
S1. The resulting phase
tions that modify the anatomy are an accepted form of
map or phase field is assumed to be continuously differen-
therapy [6], theoretical analyses of the interplay between ge-
tiable, except at the singular points. The manifold can be
ometry of the substrate and dynamics may help in the
triangulated [14] (subdividing it into a set of polygons),
therapy of cardiac arrhythmias.
where none of the edges or vertices of the polygons pass
In this paper we study spiral-wave dynamics on (punc-
through a singularity. The index I (sometimes also called the
tured) spheres with spatially inhomogeneous excitability. We
topological charge or winding number) of a curve C bound-
show for punctured spheres that the sum of indices which
ing a polygon is found by computing the line integral
characterize each spiral has to be zero. Moreover, we dem-
onstrate that topological constraints imposed by the spherical
2 I =
· dl,
1
geometry and inhomogeneity in excitability lead to the for-
C
mation of pairs of spirals, with distinctive transient dynamics
or as stable attractors, in which one spiral acts as a source
where the polygon is always traversed in a clockwise orien-
and a second as a sink leading to a source-sink pair under a
tation. By continuity of
, I must be an integer. The index
broad range of conditions. Our results explain the experi-
of a singular point is uniquely defined as the index of any
mental observations of spirals on spherical beads [3,4].
curve C provided that C encircles the point but no other
While we do not consider detailed models of cardiac-wave
singular points. The index of a curve that does not enclose
propagation, our results may apply to some generic aspects
any singular points is obviously zero.
If the manifold M has no boundaries, each edge of the
triangulation is an edge of two polygons. Since the line in-
*Electronic address: davidsen@mpipks-dresden.mpg.de
tegral adds up the change in phase along the various edges of
†Electronic address: glass@cnd.mcgill.ca
the polygon, the sum of the indices of the singular points for
‡Electronic address: rkapral@chem.utoronto.ca
a phase field in M is
1539-3755/2004/70(5)/056203(6)/$22.50
70 056203-1
©2004 The American Physical Society
DAVIDSEN, GLASS, AND KAPRAL
PHYSICAL REVIEW E 70, 056203 (2004)
III. FITZHUGH-NAGUMO EQUATION
MI = 0 ,
2
The FitzHugh-Nagumo (FHN) equation [17]
where the sum is over all singular points. This follows since
the contribution of the change in phase of each edge to the
u
u3
total integral is counted twice, but since the edge is traversed
= −1 −
+ u − v + D
2
u
u,
in opposite directions each time, the net contribution of each
t
3
edge is zero. This index theorem is applicable to tori and
other surfaces of genus different from 0. However, unlike the
more familiar Poincaré index theorem (see Ref. [15], p. 74,
v = u − v + + D 2
v
v,
4
for vector fields) the sum of the indices of the singular points
t
does not depend on the genus of the surface.
This index theorem for manifolds without boundaries can
where u r , t and v r , t are two scalar fields, 2 is the ratio
be extended to manifolds with boundaries. In the following,
of the time scales associated with the two fields, and Du and
we will consider the case of structures that arise from punc-
Dv are the constant diffusion coefficients, is a prototypical
turing orientable and compact two-dimensional differentiable
model for excitable media. We choose Du= 2 and Dv= 0. The
manifolds. The index of a hole can be uniquely defined as the
real parameters
and
characterize the local dynamics and,
index of a curve C provided that C encircles the hole but no
hence, the local excitability. Whenever 0
1,
2
1,
other holes or singular points and C is positively oriented
and
H
1 −
2 1/2 1 2 + 2 2 − 1 , the FHN system
with respect to the domain which contains the hole and is
3
is excitable. At
H a Hopf bifurcation occurs such that for
bounded by C. Applying this definition and taking the sum-
H the system exhibits oscillations. In the following we
mation in Eq. (2) over the singularities and the hole, or the
take
= 0.2, = 0.2, and
H = 0.863. . ..
holes if there is more than one hole, the index theorem can
We consider a spherical shell whose outer and inner radii
be proved by the same line of arguments as for the case of
are Re and Ri, respectively, and focus on thin spherical shells
manifolds without boundaries.
where Ri= 40, Re= 42. The radii are large enough to avoid a
This extension is important in the heart, for example,
curvature-dependent loss of excitability [18], and the shell is
where the atrium is punctured by valves and veins. In such
sufficiently thin that the dynamics is effectively 2D corre-
cases one is led to consider manifolds with holes—for ex-
sponding to the dynamics on a sphere [19]. The initial con-
ample, a sphere with a hole. A sphere with a hole is topo-
dition is a domain of an “excited” state, adjacent to a domain
logically equivalent to a disk, and, indeed, the results for
of the “refractory” state. Both domains have the forms of
disks and for spheres with holes are consistent: For the disk
slices of the same size oriented from the north to the south
D2, bounded by a curve C,
pole [20] and yield a pair of counterrotating spirals.
In order to apply the topological results to the FHN equa-
tion, it is necessary to first define the phase. We define a
D2I =
· dl,
3
phase,
r , t
based
on
the
equation
tan
r , t
C
= v r , t / u r , t if v r , t
0 and u r , t
0. Thus, singular
points at given t are points r in the medium for which
so that the sum of the indices of the singular points in the
v r , t = 0 and u r , t = 0. For each t, we obtain a continuously
disk is equal to the index of the curve C bounding the disk. If
differentiable phase map Mt =
: , t Dt that associates to
there is a single singular point on the disk, with an index of
each point in a well-defined domain Dt a phase lying on the
+1, the index of the curve bounding the disk will also be +1.
unit circle,
S1. For our FHN medium, the domain is the
Imagine now the boundary of the disk to be brought together
surface of a (punctured) sphere reduced by a finite number of
(like a drawstring bag) so that the boundary of the disk now
points where the phase is singular at fixed t.
defines a hole in the sphere. In this geometry the curve C will
Rotating spiral waves in the FHN equations are obviously
be traversed in an opposite orientation (the hole is now in-
associated with a singular point which is called the spiral
side C) from the direction it was traversed when it was the
core. In what follows, we assume that there are only single-
boundary of the disk. Now if there is a singular point with an
armed spirals so that a clockwise rotating source has an in-
index of +1 on the sphere, the index of the hole is −1, so that
dex of +1 and a counterclockwise rotating source has index
the sum of the indices is again zero.
−1. A clockwise rotating sink has an index of −1, and a
Since it is necessary to conserve the sum of the indices,
counterclockwise rotating sink has index +1. From Eq. (2), it
singularities of index ±1 usually arise and are destroyed in
is impossible to have a single rotating spiral wave on a
pairs of opposite sign [9]. An exception occurs when singu-
sphere; in addition, there must be at least another singular
larities are destroyed by collision into a boundary, so that the
point or a hole with nonzero index.
index of the singular point and the index of the boundary
For excitable media, a nonzero index of a hole implies
both change simultaneously. Another exception occurs if
that wave fronts travel permanently around the hole such that
there are singularities with index different from 1. In such
the numbers of fronts traveling clockwise and counterclock-
cases interactions between different singularities can lead to
wise are different. This includes the particularly simple case
the destruction or creation of odd numbers of singularities
of a single wave front traveling around the hole which can be
[16].
considered as a spiral wave associated with the hole.
056203-2
TOPOLOGICAL CONSTRAINTS ON SPIRAL WAVE …
PHYSICAL REVIEW E 70, 056203 (2004)
FIG. 1. (Color online) Left: spiral waves of excitation (light
fronts) on the sphere for constant
= ex emanating from spiral
cores close to the poles on an equator projection. The white arrows
show the direction of propagation. One annihilation front along the
equator can be identified. Right: sketch of the constant gradient in
the inhomogeneous case. The dashed lines are the equi-
lines and
FIG. 2. (Color online) The local excitability
r t
at the spiral
we choose
cores versus time. Gradient-induced motion of the two spiral cores
max = 1.0 and
min =
ex. The angle
describes the ori-
entation with respect to the axis from pole to pole. The results
leads to a change in the local excitability at the cores with time. The
described in the text do not depend qualitatively on the choice of
spiral period in the final state is T0= 13.2± 0.1. Four different re-
or
gimes can be identified (see text). Inset: the final bound pair of
min and
max as long as they yield stable spirals.
counterrotating spirals in regime III for
= 51.0° is shown on an
IV. SPIRAL-WAVE DYNAMICS IN SPHERICAL
equator projection such that the point of lowest excitability lies on
GEOMETRIES
the central longitude. The spiral closer to the equator has index −1
while the other one has index +1.
A. Dynamics in excitability gradients
cent of the structure of antispirals—i.e., inward moving spi-
First consider homogeneous FHN media with a constant
rals seen in oscillatory media [26]. However, the origin of
= ex 0.9. The wave fronts emanating from the spiral
this inward spiral motion in oscillatory media differs and is
cores with opposite index “annihilate” along the equator
distinct from that observed here. In oscillatory media, either
such that each spiral determines the dynamics on half of the
spirals or antispirals are stable depending on system param-
sphere (see Fig. 1, left panel)—similar to what has been ob-
eters and the wavelength diverges on the border in the pa-
served in Ref. [21] for a different excitable system. We have
rameter space between these two regimes. Thus, antispirals
shown that this behavior is robust with respect to disorder in
exist independently of spirals. This is not the case here be-
excitability with small amplitude and correlation length. If
cause the generation of an inward-moving spiral relies on the
random, uncorrelated spatial variations in
exist on length
presence of a spiral source and spherical geometry. For ex-
scales much smaller than the diameter of the spiral core me-
ample, consider the FHN system with a disk geometry and a
ander [22], the dynamics is able to average over such small-
scale inhomogeneities. The robustness explains why such
states have been experimentally observed in some chemical
reactions on spherical beads which are intrinsically inhomo-
geneous [4].
Applying a constant gradient in
as sketched in the right
panel of Fig. 1 leads to a different scenario. The time evolu-
tion of the spiral pair may be partitioned into four distinct
regimes as shown in Figs. 2 and 3. Because of the gradient,
the frequencies of the two spirals differ since a higher value
of
corresponds to lower excitability, which generally im-
plies a lower spiral frequency [23]. During a short transient
T, the spiral with the higher frequency assumes control of the
dynamics [24] on the sphere. The location on the sphere,
where wave fronts emanating from the two spiral cores an-
nihilate, moves toward the core of the low-frequency spiral.
Finally, the low-frequency spiral core is pushed farther from
the high-frequency spiral core [24,25] (see Fig. 3). After this
FIG. 3. (Color online) The distance between the two spiral cores
short transient, the wave fronts travel from pole to pole, lead-
d t versus t. Gradient-induced motion of the two spiral cores leads
ing to the creation of a source-sink pair. This (intermediate)
to a change in the distance d between the cores with time. The
state is shown in Fig. 4 and corresponds to regime I in Figs.
lower and upper curves correspond to the distance in R3 and S2,
2 and 3. Viewed from the low-excitability end of the sphere,
respectively. The dashed lines are the respective upper bounds
the waves wind into a small region about the core, reminis-
given by the size of the sphere.
056203-3
DAVIDSEN, GLASS, AND KAPRAL
PHYSICAL REVIEW E 70, 056203 (2004)
Not only does the gradient in the FHN medium change
the local excitability but it also induces a drift of the spiral
cores [29]. For our model, the drift is rather slow compared
to the transition to the source-sink pair which takes place
during the transient regime T. This can be seen in Fig. 2.
[The fluctuations in (r t ) are due to spiral meandering.] In
regime I, the dominating spiral drifts toward lower excitabil-
ity and its wave fronts continuously push the other core in
the opposite direction, thus keeping the distance d between
the cores approximately constant (see Fig. 3). The fluctua-
tions in d t are again due to meandering of the spiral cores.
Because of this drift, the local excitabilities at the two spiral
cores approach each other until they become equal.
At this point, regime II in Figs. 2 and 3 is entered. The
dynamics change drastically: the enslaved spiral reverses its
FIG. 4. (Color online) Waves of excitation on the sphere in
drift direction and regains control over its own dynamics.
regime I of the gradient-induced dynamics shown in Figs. 2 and 3.
Both spirals drift toward lower excitability. Due to the geo-
A source-sink pair has formed. For random spatial variations of
metric constraints imposed by the spherical geometry, the
excitability with a correlation length comparable to the diameter of
spirals approach each other until they form a bound pair (at
the spiral core meander, the final state is very similar to the one
t
7500
shown here [22]. Upper panel from left to right: view centered at
).
the north pole, south pole, and the equator. The source at the south
They finally reach a stable state (at t
9500) correspond-
pole has index −1 and the sink at the north pole index +1. The black
ing to regime III in Figs. 2 and 3. Neither the average dis-
circles show possible choices of C. The white arrow shows the
tance between the spirals nor the average local excitability
direction of wave propagation. Lower panel: dynamics at the north
changes further. Yet on top of the persisting unsynchronized
pole. Time increases from left to right with
t = 2.5 between
meandering of the two spiral cores, the bound pair slowly
snapshots.
moves along a (closed) equi-
curve on the surface of the
sphere. The direction of the motion depends on the initial
radial gradient in excitability such that the maximum value
condition—i.e., whether the spiral with positive index was
of
is located in the center of the disk. In this case, a
initially closer to the point of lowest excitability or to the
source-sink pair cannot occur because the high-frequency
point of highest excitability than the spiral with negative
spiral, acting as a source, would push the other (low-
index. The wave dynamics generated by this bound pair is
frequency) spiral out of the system, excluding the presence
shown in the inset of Fig. 2.
of any strong random inhomogeneities in excitability which
While kinematic theory applies only to spirals with large
may pin the low-frequency spiral and prevent its motion
cores, it is instructive to note that this theory predicts that the
(see, e.g., [27]). The lower panel of Fig. 4 shows that rem-
direction of the drift due to gradients depends on the model
nants of the low-frequency spiral persist in a small area
system and its parameters [30]. Although our simulations
around the core. It has been speculated that antispiral waves
have shown the same direction of drift for a range of param-
might occur in cardiac tissue [28]. While excitable media
eters in the meander region of the FHN phase diagram, it is
cannot support a regime of exclusive antispirals due to their
conceivable that, under different circumstances favoring a
excitable character, our results show that source-sink pairs
drift toward higher excitability, one spiral could act as a per-
with similar characteristics could form in the heart where the
manent source and a source-sink pair could be the final at-
underlying dynamics is excitable, the medium is inhomoge-
tractor. Such a scenario would also be consistent with the
neous, and the topology is similar to a (punctured) sphere.
experimental findings in Ref. [3].
Spiral dynamics of the type described above has been
observed by Maselko and Showalter [3] in experiments on
B. Punctured spheres
the excitable Belouzov-Zhabotinsky reaction on spherical
beads. They attributed the generation of spiral source-sink
Next, we consider a homogeneously excitable sphere with
pairs to inhomogeneities in the medium related to differing
a single hole. Two scenarios can be observed depending on
chemical environments. This is consistent with our findings
the location of the hole with respect to the spiral pair. If a
for systems with gradients and is further confirmed by the
spiral wave is not permanently attached to the hole, the dy-
work reported in Ref. [21]. While the generation of source-
namics is very similar to the case without any hole. If one of
sink pairs due to a gradient in the FHN medium investigated
the spirals is permanently attached to the hole, the frequency
here is only an intermediate state (but one that persists for
of this spiral is lowered. The size of the hole determines the
approximately 240 spiral periods in our simulations), random
frequency of the spiral because the wave front has to travel
spatial variations of the excitability with a correlation length
around the hole. The transient dynamics is similar to that in
comparable to the diameter of the spiral core meander or
regime T for the case with a gradient; however, no drift of
larger can lead to a final state consisting of such a source-
the spiral cores is induced and the final state is a spiral
sink pair [22]. This is due to the fact that the source can be
source-sink pair as shown in Fig. 5. Not only is the net index
trapped in a region of depressed local excitability.
conserved during the transition to a spiral source-sink pair
056203-4
TOPOLOGICAL CONSTRAINTS ON SPIRAL WAVE …
PHYSICAL REVIEW E 70, 056203 (2004)
hole, again conserving the topological charge of the hole.
V. CONCLUDING REMARKS
Inhomogeneities due to spatially varying excitability on
(punctured) spherical shells lead to complex spiral-wave dy-
namics and the formation of source-sink spiral pairs in ex-
citable media. The results presented here are immediately
applicable to excitable media in more complicated geom-
etries such as tori or multiholed tori and to situations in
which multiarmed spirals are found. This includes math-
ematical modeling of cardiac tissue. The approach taken in
this paper stresses constraints and aspects that apply to, and
must be observed in, all realistic models of the heart satisfy-
ing certain criteria of continuity. There are also implications
FIG. 5. (Color online) Waves of excitation on the punctured
for the treatment of cardiac arrhythmias. In cardiology it is
sphere propagating towards the hole. A source-sink spiral wave pair
sometimes possible to develop maps showing the timing of
has formed in the final state. Upper panel from left to right: view
the excitation over limited regions of heart [6]. In this case, a
centered at the north pole, south pole, and the equator. Lower panel:
sink might be confused for a source (of the arrhythmia), and
dynamics at the south pole. Time increases from left to right with
this might have implications for the diagnosis of the mecha-
t = 2.5 between snapshots.
nism and the choice of therapy. The current work shows how
but so is the index of the individual spirals: during the tran-
partial knowledge about what is happening in some regions
sition an outgoing-counterclockwise- (clockwise-) oriented
that could be observed might be helpful in establishing prop-
spiral
is
converted
into
an
ingoing-clockwise-
erties of dynamics that could not be observed. While the
(counterclockwise-) oriented spiral. Thus, the formation of
types of sinks we have described here have only been ob-
spiral source-sink pairs conforms with the topological con-
served in chemical media [3,4] so far, we certainly expect
straints.
their existence in the cardiological domain.
If a gradient as well as a hole is present, the spiral drift
discussed earlier also determines the final state, which de-
ACKNOWLEDGMENTS
pends on the hole’s location in the gradient field. For in-
stance, if the point of lowest excitability in the medium is on
We thank F. Chavez, M. Bär, S. G. Whittington, and D.
the hole boundary, simulations show that the spiral which is
Sumners for helpful discussions and G. Rousseau for provid-
not attached to the hole will stabilize close to this point. In
ing numerical tools. This work was supported in part by a
this final state odd numbers of wave fronts are attached to the
grant from MITACS.
[1] J. D. Murray, Mathematical Biology (Springer-Verlag, New
[11] A. T. Winfree and S. H. Strogatz, Physica D 8, 35 (1983).
York, 1989), and references herein.
[12] D. Sumners, in Graph Theory and Topology in Chemistry, ed-
[2] Chemical Waves and Patterns, edited by R. Kapral and K.
ited by R. B. King and D. Rouvray (Elsevier, Amsterdam,
Showalter (Kluwer Academic, Dordrecht, 1996), and refer-
1987), p. 3.
ences herein.
[13] I. Cruz-White, Ph.D. thesis, Florida State University, 2003.
[3] J. Maselko and K. Showalter, Nature (London) 339, 609
[14] J. G. Hocking and G. S. Young, Topology (Dover, New York,
(1989).
1988).
[4] J. Maselko and K. Showalter, React. Kinet. Catal. Lett. 42,
[15] V. Guillemin and A. Pollack, Differential Topology (Prentice
263 (1990).
Hall, Englewood Cliffs, NJ, 1974).
[5] A. T. Winfree, The Geometry of Biological Time (Springer-
[16] R. M. Zaritski and A. M. Pertsov, Phys. Rev. E 66, 066120
Verlag, New York, 2001).
(2002).
[6] K. M. Stein, S. M. Markowitz, S. Mittal, D. J. Slotwiner, S.
[17] We solve Eq. (4) numerically using an algorithm that automati-
Iwai, and B. B. Lerman, Chaos 12, 740 (2002).
cally adjusts the time step to achieve an effcient simulation
[7] Some aspects of cardiac arrhythmias may arise from scroll
while controlling the error in the solution [31].
wave filament instabilities which depend on the 3D nature of
[18] V. A. Davydov, N. Manz, O. Steinbock, and S. C. Müller,
cardiac tissue [8,9].
Europhys. Lett. 59, 344 (2002).
[8] A. T. Winfree, When Time Breaks Down (Princeton University
[19] It is computationally convenient to carry out the calculations in
Press, Princeton, NJ, 1987).
a thin spherical shell constructed from a cubic rectangular grid.
[9] S. Alonso, F. Sagués, and A. S. Mikhailov, Science 299, 1722
The shell thickness is small compared to characteristic wave-
(2003).
lengths and diffusion lengths so that the small shell thickness
[10] L. Glass, Science 198, 321 (1977).
does not influence the spiral wave dynamics.
056203-5
DAVIDSEN, GLASS, AND KAPRAL
PHYSICAL REVIEW E 70, 056203 (2004)
[20] The initial conditions are as in Ref. 32 with spherical coordi-
[27] V. N. Biktashev, A. V. Holden, S. F. Mironov, A. M. Pertsov,
nates
= / 2 and
= .
and A. V. Zaitsev, Int. J. Bifurcation Chaos Appl. Sci. Eng. 11,
[21] H. Yagisita, M. Mimura, and M. Yamada, Physica D 124, 126
1035 (2001).
(1998).
[28] V. K. Vanag and I. R. Epstein, Science 294, 835 (2001).
[22] For a Gaussian distribution of
in boxes of volume 0.53 with
[29] V. N. Biktashev and A. V. Holden, Chaos, Solitons Fractals 5,
mean 0.95 and standard deviation 0.013, the diameter of the
575 (1995).
spiral core meander is roughly 10 space units.
[
[30] A. S. Mikhailov, V. A. Davydov, and V. S. Zykov, Physica D
23] A. T. Winfree, Chaos 1, 303 (1991).
[
70, 1 (1994); V. A. Davydov, V. S. Zykov, A. S. Mikhailov,
24] V. I. Krinsky and K. I. Agladze, Physica D 8, 50 (1983).
[25] E. A. Ermakova, V. I. Krinsky, A. V. Panfilov, and A. M.
and P. K. Brazhnik, Radiophys. Quantum Electron. 31, 574
(
Pertsov, Biofizika 31, 318 (1986).
1988).
[26] Y. Gong and D. J. Christini, Phys. Rev. Lett. 90, 088302
[31] G. Rousseau and R. Kapral, Chaos 10, 812 (2000).
(2003); L. Brusch, E. Nicola, and M. Bär, ibid. 92, 089801
[32] F. Chavez, R. Kapral, G. Rousseau, and L. Glass, Chaos 11,
(2004).
757 (2001).
056203-6
