- Email: [email protected]
Chaotic meander of spiral waves in the FitzHugh-Nagumo system
Chaos, Sohrons & Fmctais Vol. 5, Nos 314, pp. 66-670. 1995 Copyright 0 1995 Elsevier Science Ltd Printed in Great Britain. AU rights reserved c960-0779/95$9.50 + .oo
Pergamon
0960-0779(93)EO@48-G
Chaotic Meander of Spiral Waves in the FitzHugh-Nagumo System H. ZHANG and A. V. HOLDEN Centre for Nonlinear
Studies and The Department
of Physiology,
LJniversity of Leeds, Leeds LS2 9JT, UK
Abstract-We study the dynamics of hypermeander of a spiral wave produced by the FitzHughNagumo model, and use the surrogate test method to identify the irregular time series of the core position as chaotic. The correlation dimension, and maximum Lyapunov exponent of this time series are estimated.
1. INTRODUCTION
Re-entrant spiral waves are often observed in effectively two-dimensional (2D) excitable systems in physics, chemistry and biology [l-3]. It is believed that recirculation of excitation in cardiac tissue underlies some cardiac arrhythmias [4]. Experimental investigations on the dynamic behaviours of spiral waves in 2D, and scroll waves in 3D excitable media have used parabolic partial differential equations [5], coupled arrays of ordinary differential equations [6,7], coupled map lattice [&lo] and cellular automata [ll-131 as models. For an unstable spiral wave, the spiral wave will break down, producing new interactive vortices that evolve into spatial temporal irregularity [4,X]. A stable spiral wave can follow rigid rotation around a circular core of radius r,,.+ the wave front does not invade this core: in models of cardiac muscle the core is depolarized to about - 40 mV [14-161. However, the tip of the spiral need not follow a circle of radius ro: it can meander in a larger circle, or floret pattern of open or closed curves. These different behaviours are found in parametrically different excitable media, and meander has been found in partial differential equation [5], CML [8-lo] and CA [11-131 models of excitable media. Systematic investigation of meander has been carried out using the FHN model [5] and Zykov kinetic model [17]. In these models, with the change of one or more parameters of the systems, the spiral wave will change from a stationary state of rigid rotation around a point or along a circle, to a non-stationary state, with the core of the spiral wave moving quasi-periodically. At some parameters, the movement of the core is irregular: this is hypermeander. The mathematical basis for meander is poorly understood, and requires a better understanding of the stability of the normal and tangential models of the spiral tip [S]. We adopt a phenomenological approach to meander, as we are primairly interested in the control of spiral waves in the context of the control of reentrant cardiac arrhythmias. We use FitzHugh-Nagumo (FHN) equations as a model of cardiac excitable medium. The parameters to produce stable, meander and hypermeander can be found in {5]. 661
H. ZHANG and A. V. HOLDEN
662
Is the irregular movement of the core of a hypermeandering spiral wave chaotic? (See refs [5,19,20].) We use the surrogate test method [21-231 to identify characteristics of the reconstructed attractor from the time series of core position, and estimate correlation dimension [24] and Lyapunov exponents [25] to quantify the dynamics of the core movement in phase space. 2. FHN MODEL OF CARDIAC TISSUE
The FHN equation is a two-variable system [26,27], which takes the form 2
= f(u, v) = ?+L
dv _ = g(u, v) = E(U + p - yv). dt This two-variable system can be taken as a general model for an excitable system, such as a single cardiac cell. A spatially extended system of FHN model is realized by using a Laplacian operator acting on the excitable variable u as a presentation of the diffusion process: al.4 -=
u - u3/3 - v
at
E -
8V
=
E(U
+
/3 -
+ Au yv).
at
In the model, the temporal and spatial scales are simply refer to time unit (tu) and space unit (su), respectively. One tu corresponds to about 10 ms, one su to about 2.5 mm. 3. NUMERICAL EXPERIMENT
The numerical experiment is carried on 420 x 420 grid of size about 105 x 105 su* (about 262.5 x 262.5 mm*), with FHN parameters at each grid point E = 0.03, y = 0.5 and p = 1.2. We use the explicit Euler integration method to solve the partial differential equation with integration time interval dt = 0.006 tu. Numerical approximation of the Laplace operator acting on the variable u takes the form:
Au=&+&=?x ax*
+;x
Ui,j+l
ay* ui+l,j+l
+
Ui,j-1
+
+
ui+l,j
-
4
’
%j
AX2
3 +
Ui-1,j
ui+l,j-l
+
ui-l,j+l
0/2w*
-I-
ui-l,j-l
-
4
x
uiJ
7
and we use no-flux boundary condition: (nA)& = 0, where n is the normal unit vector to the boundary r. A spiral wave is formed from a broken plane wave. We trace the core position of the spiral wave, and record the time series from several points in the medium. The broken plane wave evolves into a spiral wave in a transition time of 300 tu. The total integration takes 2.4 x lo4 tu. To locate a core for a spiral wave, we note that any excitable cell can be one of four different states, i.e. exciting state, excited state, refractory state and resting state. In a two-dimensional medium, all cells undergoing excitation form a wave front curve, and all
Chaotic meander of spiral waves
663
refractory cells form the wave back curve, while all excited cells form the excited region. A cell undergoing excitation is characterized by u1 < u G u2, and du/dt > 0. A refractory cell is characterised by u1 < u G u2 and du/dt < 0. An excited cell is characterized by u > u2, and resting cell is characterized by u 6 u1 and du/dt = 0 (for non-pacemaker excitable cell). ur and u2 are the selected thresholds. For a spiral wave only the core has, as neighbour, cells undergoing excitation, excited cells, refractory cells and resting cells. We trace the wave front and wave back curves for the spiral wave, obtain the crossing point of these two curves, and if it is also neighboured by excited and resting cells, then we take this point as the core position for the spiral wave. This method is economical and practical, and has been used to locate the filament for a scroll wave in a 3D coupled map lattice model of cardiac excitable medium [8]. We use this method to locate the core position and record the time intervals for the spiral wave front to arrive at a specified point of the medium. We sample the core position (x, y) time series at each 50 integration time intervals. The core path, the time series of the y-component of core position and its power spectrum, and the spiral wave period are shown in Fig. 1. The core of the spiral wave meanders irregularly in a large area of the medium, forming a complicated pattern illustrated in Fig. l(a), similar to that seen in Winfree [5]. We note that the core moves about 100 su in the x-direction, and 70 su in the y-direction. The time series from the y-component of the core position is shown in Fig. l(b). It is apparently irregular. Since the core of the spiral wave meanders in a large area of the medium, the source of the spiral wave changes with time. The period of the spiral wave changes as well, as illustrated as in Fig. l(c). We note that the period of the spiral wave oscillates in a range of 30 to 33 tu.
36 35
32
25
50
75
ld0
I
b 150 0
Fig. 1. The movement of the core for a hypermeandering spiral wave in FHN model with E = 0.03, /? = 1.2. y = 0.5. (a) The complicated pattern formed by the core path. Abscissa: x-coordinate for core position, i(iAxsu). Ordinate: y-coordinate for core position, j(jAysu); (b) the time series from y-component of the core position recording. Abscissa: time (tu). Ordinate: y coordinate for core position, j(jAysu); (c) the fluctuating periods for the spiral wave front arriving at the lattice point (210, 210). Abscissa: ith number of periods. Ordinate: periods (tu); (d) the power spectrum of the time series for y-component of core time series. Abcissa: frequency (Hz). Ordinate: log (Ps/Pi).
664
H. ZHANG and A. V. HOLDEN
The power spectrum of the time series from the y-component of core position reveals that this time series has a broad power spectrum illustrated in Fig. l(d). There are four wide peaks. If we take one time unit as 10 ms, then the four peaks are at 0.16 Hz, 1.95 Hz, 3.09 Hz and 6.09 Hz. We note that the third wide peak of the power spectrum is corresponding to the spiral frequency, the second wide peak is about half the frequency of the spiral wave, and the fourth one is about twice the frequency of the sprial wave. The most interesting wide peak is the one at 0.16 Hz. From the time series, we can see clearly the time series of y-components of the core position are modulated by a slow sinusoidal wave. We think this is a result of the boundary influence. In a unbounded medium, or a bounded medium while the core position of the spiral wave is far from each boundary of the medium, a stable stationary spiral wave will have a constant frequency and a fixed core. However, in a bounded medium, when the core position of the spiral wave is near to the boundary, the influence of the boundary is significant. Numerical experiments have demonstrated the influence of the boundary on the motion of the core of a spiral and the filament of a scroll wave [28,29]. A simple model to study this influence has been proposed [30]. For simplicity we consider a stationary spiral wave in a bounded medium. Take the core position of the spiral wave as (x, y), without any influence except the influence from the boundary of the medium, the equation for the core position is
where C, and C,, denote the influence of boundary on the speed of the movement of core. In the case of impenetrable boundary this influence can be simply taken as an exponential decaying function of the distance from the boundary: a exp (-d/R) [30], where R is a characteristic distance, and d is the distance of the core from the boundary, a is the parameter determined by the properties of the boundary. Suppose a medium is with size of L x L with impenetrable boundaries at x = 0, x = L, y = 0, y, = L. At the point (x, y), the equation of the core under the influence of all boundaries is s
= aexp(-x/R)
- aexp(-(L
- x)/R)
$
= aexp(-y/R)
- aexp(-(L
- y)/R).
The trajectory of the core is governed by dx2 -+dx2= dt2 dt*
If the core position is far from all boundaries, i.e. the distance of the core from each boundary is bigger than the characteristic distance R, or the core is at such a position that the influences from all boundaries balance each other, then the boundary will not influence the core. However, if the core position is near to one of the boundaries, say, near to x = 0 boundary, then the equation will be
Chaotic meander of spiral waves
--dx* + dx2 = $exp(-x/R) 5 + ;exp[-(L - y)/R] - ;exp(-y/R) dt2 dt* so the core will move along an ellipsoid periodically with frequency determined by the properties of the boundary. The four wide peaks in the power spectrum of the time series of the core position imply some regularity of the core movement; however, the complicated pattern formed by the core path, and the broader power spectrum show the irregularity of the time series. 4. HYPERMEANDERING
SPIRAL WAVE
The time series of the core position shown in Fig. l(b) is aperiodic, and has a broad continuous power spectrum in a low-frequency range composed of four wide peaks. Here we take a statistical approach, a surrogate test method, to study the irregularity of the time series. The surrogate test method contains two ingredients: a null hypothesis and statistical measures, based on which we test our hypothesis. In our case, the null hypothesis is: the aperiodic time series of the core position for the hypermeandering spiral wave is produced by a linear stochastic process, and the statistical measures are the characteristics of the reconstructed attractor from the time series. To apply a surrogate test method, we need to generate surrogate data which have similar properties to the original data. There is no relationship between phases in the power spectrum of a stochastic process. Thus stochastic surrogate data can be obtained by inverting a power spectrum exactly equal to that of the original time series, and randomizing its Fourier phases uniformly and independently. We generate surrogate data by the following steps [16]: first, we calculate the power spectrum of the original time series of the core position by FFT analysis, and randomize the phase of the power spectrum distribution by generating a random number uniformly distributed in a range of 0 - 2rr; second, we take the inverse FFT on the phase randomized power spectrum to obtain a new time series. This time series has the same power spectrum distribution and auto-correlation function as the original time series; however there is no relationship between phases. We take this data as surrogate data. We calculate the characteristics for the reconstructed attractors from both the original data and the surrogate data. If the original time series is a stochastic process, the results will not be affected by phase randomization; however, if the original time series is a nonlinear process, then the results of the estimated characteristics will be significantly affected. 4.1. Reconstruction of an attractor We use a time-delayed method to reconstruct an attractor from a time series [31]. Given a time series, x(t), x(t + At), . . ., x(t + nht), which is sampled at regular time interval, At, from the value of a state variable of a system. It is one-dimensional, and represents a projection rr: RM + R from the full state vector, X(t) E R”. Using the time delayed method, we create a vector, x E R” x(t) = (x(t), x(t + r), . . ., x(t + (m - 1)r)) where r = kAt is the time delay, k takes integer value, k = 1, 2, . . . m is the embedding dimension, x(t) acts as a mapping, P: RM + R”, from the full state X(t) to the reconstructed state x(t). It is proved that, when m is big enough, m > 20 + 1, where D is the fractal dimension of the attractor, this procedure gives a reliable reconstruction of the attractor from an infinite long time series [32,33].
666
4.2.
H. ZHANG and A. V. HOLDEN
Correlation
dimension
Correlation dimension is defined as the scaling index of correlation integral over correlation distance [22]. It measures the geometric properties of an attractor in the phase space. For a M-dimensional attractor, the correlation integral C,(r) is based on the statistic of the pairwise distance of two-state points and takes the form
where H is the Heaviside step function, N is the number of points for the discretized state vectors on the trajectory of attractor in the phase space. The discretized state vectors are noted by Xi, xi. II* 11indicates the normal of the vector. For a small correlation distance r, if the power law holds, ChfW N-m,r-+O
-
r
D2
where I changes in rmin- rmax,then we call D2 the correlation dimension of the attractor. In a reconstructed attractor, efficient codes to calculate the correlation dimension are available [34,35]. For a m-embedded attractor, we calculate the correlation integral for each discretised correlation distances, and take the slope of the curve of In C,(r) - In (r) as the correlation dimension, Wm,
r> = d(ln G(r))/d(ln
(r)).
We change the embedding dimension m from a low dimension to a high dimension, until the calculated D2 gets saturated. We calculated the correlation dimension for the reconstructed attractors from the time series of core position and surrogate data. The parameters we used here are: r = 200 tu, i.e. 4 times the recording time interval, m changes from 2 to 8 for the time series of the recording of the hypermeandering spiral wave, while m changes from 2 to 10 for the surrogate time series. The results are illustrated in Fig. 2, with left side for the time series of core position and the right side for the surrogate time series. Figure 2(a) shows the curves of the logarithms of the correlation integral versus the logarithms of the correlation distance for each embedding dimension. Figure 2(b) displays the slopes of the curves in Fig. 2(a). We note that, in both cases, there is a large region of the correlation distance where the curves in Fig. 2(a) are linear, i.e. the slopes of the curves are plateaux. By fitting the curves linearly in this region, we get the slopes of the fitted lines and take it as &(m). Figure 2(c) shows the results of Da(m) for different embedding dimensions m. We note that at m = 3, the correlation dimension of the time series saturates at about 1.13, and at m = 5, the correlation dimension of the surrogate time series situates at about 2.79. The difference between them is significant.
4.3.
Lyapunov
exponent
The Lyapunov exponent is a measure of the sensitivity of the system to initial conditions. For an M-dimensional phase space, we consider an infinitesimal M-sphere of initial conditions, with the centre of the M-sphere on the attractor, and trace the evolution of this M-sphere of initial conditions. The averaged growth rate of the norm of the ith principal axis a,(t) of the M-sphere gives the ith Lyapunov exponent [35,36] Ai = lim,,,
1log, -Ilai(t)llbitss-1 t ( Ibi(t 1 ’
Chaotic meander of spiral waves
Fig. 2. The reconstructed attractor from the y-component of the core position recording.
If a system is in a chaotic state characterized by a sensitivity to initial conditions, then with time the M-sphere of initial conditions will evolve into a M-ellipsoid, since the M-sphere of initial conditions expands in at least one direction, and contracts in other directions. The algorithms for calculating the maximal [25], and all of Lyapunov exponents [36] are available for an attractor reconstructed from a time series. Here we calculate the maximal Lyapunov exponent for the reconstructed attractor using Wolf’s method [25]. The parameters we used here are: r = 200 tu and m = 8 for both cases. The results are given in Fig. 3. Figure 3(a) is the estimated averaged maximum Lyapunov exponent with time for the core time series, and 3(b) is the estimated averaged maximum Lyapunov exponent for surrogate data. In both cases, the convergence of the averaged results is good. With the time evolution, the converged maximum Lyapunov exponent is 2.6 bitss-’ for the core time series, while 5.5 bitss-’ for the surrogate time series. So the difference between them is also significant. 5. CONCLUSION
The estimated results of the characteristics of the reconstructed attractor for the core time series and its surrogate time series are significantly different. The characteristics we used here are correlation dimension, and the maximal Lyapunov exponent. The null hypothesis that the irregular time series of core position for a hypermeander spiral wave we study is governed by a linear stochastic process can be rejected. The irregularity of the core time series is produced by a nonlinear process.
H. ZHANG and A. V. HOLDEN
668 a
0.75
0.5
0.25
1.25
1.0
1.5
2.5
0 1.0
1.1
1.2
1.3
l.4
1.5
1.6
17
1.5-
l.O-
0.5-
0.0
I I I I, 2
I 3
/ I,
I I, 4
I 5
I I,
I I I/, 6
I I I I 7
r 6
2
3
4
5
6
7
6
9
10
Fig. 3. The estimated correlation dimension for the time series of core position (left) and surrogate data (right). (a) The logarithms of correlation integral versus the logarithms of correlation distance for different embedding dimension m. m changes from 2 to 8 for the core time series, and m changes from 2 to 10 for surrogate time series. Abscissa: correlation distance. Ordinate: correlation integral; (b) the slopes of the curves in (a). The plateau region shows good linear scaling for the curves of (a); (c) the estimated correlation dimension with the embedding dimension. Abscissa: embedding dimension. Ordinate: correlation dimension.
The estimated correlation dimension, and the maximal Lyapunov exponent confirm that the hypermeandering spiral wave in the FitzHugh-Nagumo model of cardiac excitable medium at the parameter E = 0.03, y = 0.5, p = 1.2 is chaotic. The numerical experiments demonstrate that the period of the spiral wave front arriving a certain point in the medium is irregular. The period of the spiral wave front is more irregular close to the core region, and further from the core it is less irregular. Thus the spatio-temporal chaos associated with hypermeander is localized, and measures of its irregularity decrease with distance from the vortex centre. It is an example of localized spatio-temporal chaos.
Chaotic meander of spiral waves
1
2
3
4
669
5
6
7
*lo3
5-
0
,III,III,/I,,,,,,,,,,,,, 1
2
3
4
5
*lo'
Fig. 4. The estimated averaged maximum Lyapunov exponent with time for the core time series and surrogate data. Abcissa: time (tu). Ordinate: the estimated biggest Lyapunov exponent. (a) The maximum Lyapunov exponent for the core time series: (b) the maximum Lyapunov exponent for surrogate data.
5
10
15
20
25
30
Fig. 5. The fluctuating periods for the hypermeandering spiral wave front arriving at points with different distances from the mean core position over the medium. The mean point of the core region is at (i, j) = (55,125). Abscissa: the ith number of period. Ordinate: the period (in numbers of iterations). .: at the point (i, j) = (70,70). n : at the point (i, j) = (210,210). 0: at point (i, j) = (350,350). Near to the centre point of the core region of the spiral wave the period of wave front is more irregular; far from the centre point of the core region, the period of wave front is less irregular.
670
H. ZHANG and A. V. HOLDEN
REFERENCES 1. A. T. Winfree, The geometry of biological time. Springer, New York (1980); When Time Breaks Down. Princeton University Press, Princeton, NJ (1987). 2. H. L. Swinnev and V. I. Krinskv.I Waves and Patterns in Chemical and Biolonical Media. North-Holland. Amsterdam (1691). 3. A. V. Holden, M. M. Markus and H. G. Othmer, Nonlinear Wave Processes in Excitable Media, NATO ASI Series 244B. Plenum Press, New York (1991). 4. J. Jalife, Mathematical approaches to cardiac arrhythmias, Ann NY Acad. Sci. 591, (1990). 5. A. T. Winfree, Varieties of spiral wave behaviour: An experimentalist’s approach to the theory of excitable media, Chaos 1, 303-334 (1991). 6. F. J. L. van Capelle and D. Durrer, Computer simulation of arrhythmias in a network of coupled excitable elements, Circ. Res. 47, 454-466 (1980). 7. R. L. Winslow, A. L. Kimball, A. Varghese and D. Noble, Simulation of cardiac sinus and atria1 network dynamics on the connect machine, Physica D64,281-298 (1993). 8. A. V. Holden and H. Zhang, Modelling propagation and re-entry in anisotropic and smoothly heterogeneous cardiac tissue, .I. Chem. Sot. Faraday Trans., 89,2833-2837 (1993). 9. D. Barkley, A model for fast computer simulation of waves in excitable media, Physica D49, 61-70 (1991). 10. R. Kapral, Discrete models for chemical reacting systems, J. Math. Chem. 6, 113-163 (1991). 11. M. Markus and B. Hess, Isotropic cellular automata for excitable media, Nature Lond. 347, 56-58, (1990). 12. A. T. Winfree, E. M. Winfree and H. Seifert, Organisation centres in a cellular excitable medium, Physica D17, 109-115 (1985). 13. M. Gerhardt, H. Schuster and J. Tyson, A cellular automaton model of excitable media including curvature and dispersion, Science 247, 1563-1566 (1991). 14. A. Panfilov and A. V. Holden, Spatio-temporal chaos in model of cardiac electrical activity, Int. J. Bifurcation & Chaos 1, 219-225 (1991). 15. H. Zhang and N. Patel, Spiral wave breakdown in an excitable medium model of cardiac tissue, Chaos, Soliton & Fractals 5(3/4), 635-643 (1995). 16. M. Courtemanche and A. T. Winfree, Re-entrant rotating waves in Beeler-Reuter based model of two-dimensional cardiac conduction, Int. J. Bifurcation & Chaos 1, 431-444 (1991). 17. E. Lugosi, Analysis of meandering in Zykov-kinetics, Physica D40, 331-337 (1989). 18. P. Pelce and J. Sun, On the stability of steadily rotating waves in the free boundary formulation, Physica D63, 273 (1993). 19. 0. E. Rossler and C. Kahlert, Winfree meandering in a 2-dimensional 2-variable excitable medium, Z. Naturforsch A34, 565-570 (1979). 20. Z. Nagy-Ungvarai, J. Ungvarai and S. C. Muller, Complexity in spiral wave dynamics Chaos 3, 15-19 (1993). 21. A. Provenzale, L. A. Smith, R. Vio and G. Mutante, Distinguishing between low-dimensional dynamics and randomness in measured time series, Physica D58, 31-49 (1992). 22. L. A. Smith, Identification and prediction of low-dimensional dynamics, Physica D58, 50-76 (1992). 23. J. Theiler, S. Eubank, A. Longtin, B. Galdrikan and J. D. Farmer, Testing for nonlinearity in time series: the method of surrogate data, Physica D58, 77-94 (1992). 24. P. Grassberger and I. Procaccia, Phy. Rev. Lett 50, 346 (1983). 25. A. Wolf, J. B. Swift, H. L. Swinney and J. A. Vastano, Physica D16, 285 (1985). 26. R. FitzHugh, Thresholds and plateaus in the Hodgkin-Huxley nerve equations, J. Gen. Phys. 43, 867-896 (1960). 27. R. FitzHugh, Impulses and physiological states in theoretical models of nerve membrane, Biophys. 1, 445-466 (1961). 28. A. N. Rudenko and A. V. Panfilov, Drift and interaction of vortices in a two-dimensional inhomogeneous active medium, Studia Biophysics 98, 183-188 (1983). 29. P. J. Nandapurkar and A. T. Winfree, Dynamical stability of untwisted scroll rings in excitable media, Physica D35, 277-288 (1989). 30. V. N. Biktashev, Drift in a reverberator in an active medium due to the interaction with boundaries, in Nonlinear Waves Vol 2. edited by A. V. Gaponov-Grekhov, M. I. Rabinovich and J. Engelbrecht. Springer, Berlin (1989). 31. N. H. Packard, J. P. Crutchfied, J. D. Farmer and R. S. Shaw, Phy. Rev. Lett. 45,712 (1980). 32. F. Takens, in Proc. Dynamical Systems and Turbulence, Warwick, 1980, Lecture Notes in Mathematics 898, Springer, Berlin (1981). 33. R. Mane, in Proc. Dynamical Systems and Turbulence, Warwick, 1980, Lecture Notes in Mathematics 898, Springer, Berlin (1981). 34. P. Grassberger, Phy. Lett A148, 63 (1990). 35. A. V. Holden, J..Hyde and H. Zhang, Computing with the unpredictable chaotic dynamics and fractal structure in the brain, in Application of Fractal and Chaos, edited by R. A. Eanshaw. Springer, Berlin (1993). 36. J.-P. Eckmann, S. 0. Kamphorst, D. Rulle and S. Ciliberto, Phy. Rev. A34, 4971 (1986).
Pergamon
0960-0779(93)EO@48-G
Chaotic Meander of Spiral Waves in the FitzHugh-Nagumo System H. ZHANG and A. V. HOLDEN Centre for Nonlinear
Studies and The Department
of Physiology,
LJniversity of Leeds, Leeds LS2 9JT, UK
Abstract-We study the dynamics of hypermeander of a spiral wave produced by the FitzHughNagumo model, and use the surrogate test method to identify the irregular time series of the core position as chaotic. The correlation dimension, and maximum Lyapunov exponent of this time series are estimated.
1. INTRODUCTION
Re-entrant spiral waves are often observed in effectively two-dimensional (2D) excitable systems in physics, chemistry and biology [l-3]. It is believed that recirculation of excitation in cardiac tissue underlies some cardiac arrhythmias [4]. Experimental investigations on the dynamic behaviours of spiral waves in 2D, and scroll waves in 3D excitable media have used parabolic partial differential equations [5], coupled arrays of ordinary differential equations [6,7], coupled map lattice [&lo] and cellular automata [ll-131 as models. For an unstable spiral wave, the spiral wave will break down, producing new interactive vortices that evolve into spatial temporal irregularity [4,X]. A stable spiral wave can follow rigid rotation around a circular core of radius r,,.+ the wave front does not invade this core: in models of cardiac muscle the core is depolarized to about - 40 mV [14-161. However, the tip of the spiral need not follow a circle of radius ro: it can meander in a larger circle, or floret pattern of open or closed curves. These different behaviours are found in parametrically different excitable media, and meander has been found in partial differential equation [5], CML [8-lo] and CA [11-131 models of excitable media. Systematic investigation of meander has been carried out using the FHN model [5] and Zykov kinetic model [17]. In these models, with the change of one or more parameters of the systems, the spiral wave will change from a stationary state of rigid rotation around a point or along a circle, to a non-stationary state, with the core of the spiral wave moving quasi-periodically. At some parameters, the movement of the core is irregular: this is hypermeander. The mathematical basis for meander is poorly understood, and requires a better understanding of the stability of the normal and tangential models of the spiral tip [S]. We adopt a phenomenological approach to meander, as we are primairly interested in the control of spiral waves in the context of the control of reentrant cardiac arrhythmias. We use FitzHugh-Nagumo (FHN) equations as a model of cardiac excitable medium. The parameters to produce stable, meander and hypermeander can be found in {5]. 661
H. ZHANG and A. V. HOLDEN
662
Is the irregular movement of the core of a hypermeandering spiral wave chaotic? (See refs [5,19,20].) We use the surrogate test method [21-231 to identify characteristics of the reconstructed attractor from the time series of core position, and estimate correlation dimension [24] and Lyapunov exponents [25] to quantify the dynamics of the core movement in phase space. 2. FHN MODEL OF CARDIAC TISSUE
The FHN equation is a two-variable system [26,27], which takes the form 2
= f(u, v) = ?+L
dv _ = g(u, v) = E(U + p - yv). dt This two-variable system can be taken as a general model for an excitable system, such as a single cardiac cell. A spatially extended system of FHN model is realized by using a Laplacian operator acting on the excitable variable u as a presentation of the diffusion process: al.4 -=
u - u3/3 - v
at
E -
8V
=
E(U
+
/3 -
+ Au yv).
at
In the model, the temporal and spatial scales are simply refer to time unit (tu) and space unit (su), respectively. One tu corresponds to about 10 ms, one su to about 2.5 mm. 3. NUMERICAL EXPERIMENT
The numerical experiment is carried on 420 x 420 grid of size about 105 x 105 su* (about 262.5 x 262.5 mm*), with FHN parameters at each grid point E = 0.03, y = 0.5 and p = 1.2. We use the explicit Euler integration method to solve the partial differential equation with integration time interval dt = 0.006 tu. Numerical approximation of the Laplace operator acting on the variable u takes the form:
Au=&+&=?x ax*
+;x
Ui,j+l
ay* ui+l,j+l
+
Ui,j-1
+
+
ui+l,j
-
4
’
%j
AX2
3 +
Ui-1,j
ui+l,j-l
+
ui-l,j+l
0/2w*
-I-
ui-l,j-l
-
4
x
uiJ
7
and we use no-flux boundary condition: (nA)& = 0, where n is the normal unit vector to the boundary r. A spiral wave is formed from a broken plane wave. We trace the core position of the spiral wave, and record the time series from several points in the medium. The broken plane wave evolves into a spiral wave in a transition time of 300 tu. The total integration takes 2.4 x lo4 tu. To locate a core for a spiral wave, we note that any excitable cell can be one of four different states, i.e. exciting state, excited state, refractory state and resting state. In a two-dimensional medium, all cells undergoing excitation form a wave front curve, and all
Chaotic meander of spiral waves
663
refractory cells form the wave back curve, while all excited cells form the excited region. A cell undergoing excitation is characterized by u1 < u G u2, and du/dt > 0. A refractory cell is characterised by u1 < u G u2 and du/dt < 0. An excited cell is characterized by u > u2, and resting cell is characterized by u 6 u1 and du/dt = 0 (for non-pacemaker excitable cell). ur and u2 are the selected thresholds. For a spiral wave only the core has, as neighbour, cells undergoing excitation, excited cells, refractory cells and resting cells. We trace the wave front and wave back curves for the spiral wave, obtain the crossing point of these two curves, and if it is also neighboured by excited and resting cells, then we take this point as the core position for the spiral wave. This method is economical and practical, and has been used to locate the filament for a scroll wave in a 3D coupled map lattice model of cardiac excitable medium [8]. We use this method to locate the core position and record the time intervals for the spiral wave front to arrive at a specified point of the medium. We sample the core position (x, y) time series at each 50 integration time intervals. The core path, the time series of the y-component of core position and its power spectrum, and the spiral wave period are shown in Fig. 1. The core of the spiral wave meanders irregularly in a large area of the medium, forming a complicated pattern illustrated in Fig. l(a), similar to that seen in Winfree [5]. We note that the core moves about 100 su in the x-direction, and 70 su in the y-direction. The time series from the y-component of the core position is shown in Fig. l(b). It is apparently irregular. Since the core of the spiral wave meanders in a large area of the medium, the source of the spiral wave changes with time. The period of the spiral wave changes as well, as illustrated as in Fig. l(c). We note that the period of the spiral wave oscillates in a range of 30 to 33 tu.
36 35
32
25
50
75
ld0
I
b 150 0
Fig. 1. The movement of the core for a hypermeandering spiral wave in FHN model with E = 0.03, /? = 1.2. y = 0.5. (a) The complicated pattern formed by the core path. Abscissa: x-coordinate for core position, i(iAxsu). Ordinate: y-coordinate for core position, j(jAysu); (b) the time series from y-component of the core position recording. Abscissa: time (tu). Ordinate: y coordinate for core position, j(jAysu); (c) the fluctuating periods for the spiral wave front arriving at the lattice point (210, 210). Abscissa: ith number of periods. Ordinate: periods (tu); (d) the power spectrum of the time series for y-component of core time series. Abcissa: frequency (Hz). Ordinate: log (Ps/Pi).
664
H. ZHANG and A. V. HOLDEN
The power spectrum of the time series from the y-component of core position reveals that this time series has a broad power spectrum illustrated in Fig. l(d). There are four wide peaks. If we take one time unit as 10 ms, then the four peaks are at 0.16 Hz, 1.95 Hz, 3.09 Hz and 6.09 Hz. We note that the third wide peak of the power spectrum is corresponding to the spiral frequency, the second wide peak is about half the frequency of the spiral wave, and the fourth one is about twice the frequency of the sprial wave. The most interesting wide peak is the one at 0.16 Hz. From the time series, we can see clearly the time series of y-components of the core position are modulated by a slow sinusoidal wave. We think this is a result of the boundary influence. In a unbounded medium, or a bounded medium while the core position of the spiral wave is far from each boundary of the medium, a stable stationary spiral wave will have a constant frequency and a fixed core. However, in a bounded medium, when the core position of the spiral wave is near to the boundary, the influence of the boundary is significant. Numerical experiments have demonstrated the influence of the boundary on the motion of the core of a spiral and the filament of a scroll wave [28,29]. A simple model to study this influence has been proposed [30]. For simplicity we consider a stationary spiral wave in a bounded medium. Take the core position of the spiral wave as (x, y), without any influence except the influence from the boundary of the medium, the equation for the core position is
where C, and C,, denote the influence of boundary on the speed of the movement of core. In the case of impenetrable boundary this influence can be simply taken as an exponential decaying function of the distance from the boundary: a exp (-d/R) [30], where R is a characteristic distance, and d is the distance of the core from the boundary, a is the parameter determined by the properties of the boundary. Suppose a medium is with size of L x L with impenetrable boundaries at x = 0, x = L, y = 0, y, = L. At the point (x, y), the equation of the core under the influence of all boundaries is s
= aexp(-x/R)
- aexp(-(L
- x)/R)
$
= aexp(-y/R)
- aexp(-(L
- y)/R).
The trajectory of the core is governed by dx2 -+dx2= dt2 dt*
If the core position is far from all boundaries, i.e. the distance of the core from each boundary is bigger than the characteristic distance R, or the core is at such a position that the influences from all boundaries balance each other, then the boundary will not influence the core. However, if the core position is near to one of the boundaries, say, near to x = 0 boundary, then the equation will be
Chaotic meander of spiral waves
--dx* + dx2 = $exp(-x/R) 5 + ;exp[-(L - y)/R] - ;exp(-y/R) dt2 dt* so the core will move along an ellipsoid periodically with frequency determined by the properties of the boundary. The four wide peaks in the power spectrum of the time series of the core position imply some regularity of the core movement; however, the complicated pattern formed by the core path, and the broader power spectrum show the irregularity of the time series. 4. HYPERMEANDERING
SPIRAL WAVE
The time series of the core position shown in Fig. l(b) is aperiodic, and has a broad continuous power spectrum in a low-frequency range composed of four wide peaks. Here we take a statistical approach, a surrogate test method, to study the irregularity of the time series. The surrogate test method contains two ingredients: a null hypothesis and statistical measures, based on which we test our hypothesis. In our case, the null hypothesis is: the aperiodic time series of the core position for the hypermeandering spiral wave is produced by a linear stochastic process, and the statistical measures are the characteristics of the reconstructed attractor from the time series. To apply a surrogate test method, we need to generate surrogate data which have similar properties to the original data. There is no relationship between phases in the power spectrum of a stochastic process. Thus stochastic surrogate data can be obtained by inverting a power spectrum exactly equal to that of the original time series, and randomizing its Fourier phases uniformly and independently. We generate surrogate data by the following steps [16]: first, we calculate the power spectrum of the original time series of the core position by FFT analysis, and randomize the phase of the power spectrum distribution by generating a random number uniformly distributed in a range of 0 - 2rr; second, we take the inverse FFT on the phase randomized power spectrum to obtain a new time series. This time series has the same power spectrum distribution and auto-correlation function as the original time series; however there is no relationship between phases. We take this data as surrogate data. We calculate the characteristics for the reconstructed attractors from both the original data and the surrogate data. If the original time series is a stochastic process, the results will not be affected by phase randomization; however, if the original time series is a nonlinear process, then the results of the estimated characteristics will be significantly affected. 4.1. Reconstruction of an attractor We use a time-delayed method to reconstruct an attractor from a time series [31]. Given a time series, x(t), x(t + At), . . ., x(t + nht), which is sampled at regular time interval, At, from the value of a state variable of a system. It is one-dimensional, and represents a projection rr: RM + R from the full state vector, X(t) E R”. Using the time delayed method, we create a vector, x E R” x(t) = (x(t), x(t + r), . . ., x(t + (m - 1)r)) where r = kAt is the time delay, k takes integer value, k = 1, 2, . . . m is the embedding dimension, x(t) acts as a mapping, P: RM + R”, from the full state X(t) to the reconstructed state x(t). It is proved that, when m is big enough, m > 20 + 1, where D is the fractal dimension of the attractor, this procedure gives a reliable reconstruction of the attractor from an infinite long time series [32,33].
666
4.2.
H. ZHANG and A. V. HOLDEN
Correlation
dimension
Correlation dimension is defined as the scaling index of correlation integral over correlation distance [22]. It measures the geometric properties of an attractor in the phase space. For a M-dimensional attractor, the correlation integral C,(r) is based on the statistic of the pairwise distance of two-state points and takes the form
where H is the Heaviside step function, N is the number of points for the discretized state vectors on the trajectory of attractor in the phase space. The discretized state vectors are noted by Xi, xi. II* 11indicates the normal of the vector. For a small correlation distance r, if the power law holds, ChfW N-m,r-+O
-
r
D2
where I changes in rmin- rmax,then we call D2 the correlation dimension of the attractor. In a reconstructed attractor, efficient codes to calculate the correlation dimension are available [34,35]. For a m-embedded attractor, we calculate the correlation integral for each discretised correlation distances, and take the slope of the curve of In C,(r) - In (r) as the correlation dimension, Wm,
r> = d(ln G(r))/d(ln
(r)).
We change the embedding dimension m from a low dimension to a high dimension, until the calculated D2 gets saturated. We calculated the correlation dimension for the reconstructed attractors from the time series of core position and surrogate data. The parameters we used here are: r = 200 tu, i.e. 4 times the recording time interval, m changes from 2 to 8 for the time series of the recording of the hypermeandering spiral wave, while m changes from 2 to 10 for the surrogate time series. The results are illustrated in Fig. 2, with left side for the time series of core position and the right side for the surrogate time series. Figure 2(a) shows the curves of the logarithms of the correlation integral versus the logarithms of the correlation distance for each embedding dimension. Figure 2(b) displays the slopes of the curves in Fig. 2(a). We note that, in both cases, there is a large region of the correlation distance where the curves in Fig. 2(a) are linear, i.e. the slopes of the curves are plateaux. By fitting the curves linearly in this region, we get the slopes of the fitted lines and take it as &(m). Figure 2(c) shows the results of Da(m) for different embedding dimensions m. We note that at m = 3, the correlation dimension of the time series saturates at about 1.13, and at m = 5, the correlation dimension of the surrogate time series situates at about 2.79. The difference between them is significant.
4.3.
Lyapunov
exponent
The Lyapunov exponent is a measure of the sensitivity of the system to initial conditions. For an M-dimensional phase space, we consider an infinitesimal M-sphere of initial conditions, with the centre of the M-sphere on the attractor, and trace the evolution of this M-sphere of initial conditions. The averaged growth rate of the norm of the ith principal axis a,(t) of the M-sphere gives the ith Lyapunov exponent [35,36] Ai = lim,,,
1log, -Ilai(t)llbitss-1 t ( Ibi(t 1 ’
Chaotic meander of spiral waves
Fig. 2. The reconstructed attractor from the y-component of the core position recording.
If a system is in a chaotic state characterized by a sensitivity to initial conditions, then with time the M-sphere of initial conditions will evolve into a M-ellipsoid, since the M-sphere of initial conditions expands in at least one direction, and contracts in other directions. The algorithms for calculating the maximal [25], and all of Lyapunov exponents [36] are available for an attractor reconstructed from a time series. Here we calculate the maximal Lyapunov exponent for the reconstructed attractor using Wolf’s method [25]. The parameters we used here are: r = 200 tu and m = 8 for both cases. The results are given in Fig. 3. Figure 3(a) is the estimated averaged maximum Lyapunov exponent with time for the core time series, and 3(b) is the estimated averaged maximum Lyapunov exponent for surrogate data. In both cases, the convergence of the averaged results is good. With the time evolution, the converged maximum Lyapunov exponent is 2.6 bitss-’ for the core time series, while 5.5 bitss-’ for the surrogate time series. So the difference between them is also significant. 5. CONCLUSION
The estimated results of the characteristics of the reconstructed attractor for the core time series and its surrogate time series are significantly different. The characteristics we used here are correlation dimension, and the maximal Lyapunov exponent. The null hypothesis that the irregular time series of core position for a hypermeander spiral wave we study is governed by a linear stochastic process can be rejected. The irregularity of the core time series is produced by a nonlinear process.
H. ZHANG and A. V. HOLDEN
668 a
0.75
0.5
0.25
1.25
1.0
1.5
2.5
0 1.0
1.1
1.2
1.3
l.4
1.5
1.6
17
1.5-
l.O-
0.5-
0.0
I I I I, 2
I 3
/ I,
I I, 4
I 5
I I,
I I I/, 6
I I I I 7
r 6
2
3
4
5
6
7
6
9
10
Fig. 3. The estimated correlation dimension for the time series of core position (left) and surrogate data (right). (a) The logarithms of correlation integral versus the logarithms of correlation distance for different embedding dimension m. m changes from 2 to 8 for the core time series, and m changes from 2 to 10 for surrogate time series. Abscissa: correlation distance. Ordinate: correlation integral; (b) the slopes of the curves in (a). The plateau region shows good linear scaling for the curves of (a); (c) the estimated correlation dimension with the embedding dimension. Abscissa: embedding dimension. Ordinate: correlation dimension.
The estimated correlation dimension, and the maximal Lyapunov exponent confirm that the hypermeandering spiral wave in the FitzHugh-Nagumo model of cardiac excitable medium at the parameter E = 0.03, y = 0.5, p = 1.2 is chaotic. The numerical experiments demonstrate that the period of the spiral wave front arriving a certain point in the medium is irregular. The period of the spiral wave front is more irregular close to the core region, and further from the core it is less irregular. Thus the spatio-temporal chaos associated with hypermeander is localized, and measures of its irregularity decrease with distance from the vortex centre. It is an example of localized spatio-temporal chaos.
Chaotic meander of spiral waves
1
2
3
4
669
5
6
7
*lo3
5-
0
,III,III,/I,,,,,,,,,,,,, 1
2
3
4
5
*lo'
Fig. 4. The estimated averaged maximum Lyapunov exponent with time for the core time series and surrogate data. Abcissa: time (tu). Ordinate: the estimated biggest Lyapunov exponent. (a) The maximum Lyapunov exponent for the core time series: (b) the maximum Lyapunov exponent for surrogate data.
5
10
15
20
25
30
Fig. 5. The fluctuating periods for the hypermeandering spiral wave front arriving at points with different distances from the mean core position over the medium. The mean point of the core region is at (i, j) = (55,125). Abscissa: the ith number of period. Ordinate: the period (in numbers of iterations). .: at the point (i, j) = (70,70). n : at the point (i, j) = (210,210). 0: at point (i, j) = (350,350). Near to the centre point of the core region of the spiral wave the period of wave front is more irregular; far from the centre point of the core region, the period of wave front is less irregular.
670
H. ZHANG and A. V. HOLDEN
REFERENCES 1. A. T. Winfree, The geometry of biological time. Springer, New York (1980); When Time Breaks Down. Princeton University Press, Princeton, NJ (1987). 2. H. L. Swinnev and V. I. Krinskv.I Waves and Patterns in Chemical and Biolonical Media. North-Holland. Amsterdam (1691). 3. A. V. Holden, M. M. Markus and H. G. Othmer, Nonlinear Wave Processes in Excitable Media, NATO ASI Series 244B. Plenum Press, New York (1991). 4. J. Jalife, Mathematical approaches to cardiac arrhythmias, Ann NY Acad. Sci. 591, (1990). 5. A. T. Winfree, Varieties of spiral wave behaviour: An experimentalist’s approach to the theory of excitable media, Chaos 1, 303-334 (1991). 6. F. J. L. van Capelle and D. Durrer, Computer simulation of arrhythmias in a network of coupled excitable elements, Circ. Res. 47, 454-466 (1980). 7. R. L. Winslow, A. L. Kimball, A. Varghese and D. Noble, Simulation of cardiac sinus and atria1 network dynamics on the connect machine, Physica D64,281-298 (1993). 8. A. V. Holden and H. Zhang, Modelling propagation and re-entry in anisotropic and smoothly heterogeneous cardiac tissue, .I. Chem. Sot. Faraday Trans., 89,2833-2837 (1993). 9. D. Barkley, A model for fast computer simulation of waves in excitable media, Physica D49, 61-70 (1991). 10. R. Kapral, Discrete models for chemical reacting systems, J. Math. Chem. 6, 113-163 (1991). 11. M. Markus and B. Hess, Isotropic cellular automata for excitable media, Nature Lond. 347, 56-58, (1990). 12. A. T. Winfree, E. M. Winfree and H. Seifert, Organisation centres in a cellular excitable medium, Physica D17, 109-115 (1985). 13. M. Gerhardt, H. Schuster and J. Tyson, A cellular automaton model of excitable media including curvature and dispersion, Science 247, 1563-1566 (1991). 14. A. Panfilov and A. V. Holden, Spatio-temporal chaos in model of cardiac electrical activity, Int. J. Bifurcation & Chaos 1, 219-225 (1991). 15. H. Zhang and N. Patel, Spiral wave breakdown in an excitable medium model of cardiac tissue, Chaos, Soliton & Fractals 5(3/4), 635-643 (1995). 16. M. Courtemanche and A. T. Winfree, Re-entrant rotating waves in Beeler-Reuter based model of two-dimensional cardiac conduction, Int. J. Bifurcation & Chaos 1, 431-444 (1991). 17. E. Lugosi, Analysis of meandering in Zykov-kinetics, Physica D40, 331-337 (1989). 18. P. Pelce and J. Sun, On the stability of steadily rotating waves in the free boundary formulation, Physica D63, 273 (1993). 19. 0. E. Rossler and C. Kahlert, Winfree meandering in a 2-dimensional 2-variable excitable medium, Z. Naturforsch A34, 565-570 (1979). 20. Z. Nagy-Ungvarai, J. Ungvarai and S. C. Muller, Complexity in spiral wave dynamics Chaos 3, 15-19 (1993). 21. A. Provenzale, L. A. Smith, R. Vio and G. Mutante, Distinguishing between low-dimensional dynamics and randomness in measured time series, Physica D58, 31-49 (1992). 22. L. A. Smith, Identification and prediction of low-dimensional dynamics, Physica D58, 50-76 (1992). 23. J. Theiler, S. Eubank, A. Longtin, B. Galdrikan and J. D. Farmer, Testing for nonlinearity in time series: the method of surrogate data, Physica D58, 77-94 (1992). 24. P. Grassberger and I. Procaccia, Phy. Rev. Lett 50, 346 (1983). 25. A. Wolf, J. B. Swift, H. L. Swinney and J. A. Vastano, Physica D16, 285 (1985). 26. R. FitzHugh, Thresholds and plateaus in the Hodgkin-Huxley nerve equations, J. Gen. Phys. 43, 867-896 (1960). 27. R. FitzHugh, Impulses and physiological states in theoretical models of nerve membrane, Biophys. 1, 445-466 (1961). 28. A. N. Rudenko and A. V. Panfilov, Drift and interaction of vortices in a two-dimensional inhomogeneous active medium, Studia Biophysics 98, 183-188 (1983). 29. P. J. Nandapurkar and A. T. Winfree, Dynamical stability of untwisted scroll rings in excitable media, Physica D35, 277-288 (1989). 30. V. N. Biktashev, Drift in a reverberator in an active medium due to the interaction with boundaries, in Nonlinear Waves Vol 2. edited by A. V. Gaponov-Grekhov, M. I. Rabinovich and J. Engelbrecht. Springer, Berlin (1989). 31. N. H. Packard, J. P. Crutchfied, J. D. Farmer and R. S. Shaw, Phy. Rev. Lett. 45,712 (1980). 32. F. Takens, in Proc. Dynamical Systems and Turbulence, Warwick, 1980, Lecture Notes in Mathematics 898, Springer, Berlin (1981). 33. R. Mane, in Proc. Dynamical Systems and Turbulence, Warwick, 1980, Lecture Notes in Mathematics 898, Springer, Berlin (1981). 34. P. Grassberger, Phy. Lett A148, 63 (1990). 35. A. V. Holden, J..Hyde and H. Zhang, Computing with the unpredictable chaotic dynamics and fractal structure in the brain, in Application of Fractal and Chaos, edited by R. A. Eanshaw. Springer, Berlin (1993). 36. J.-P. Eckmann, S. 0. Kamphorst, D. Rulle and S. Ciliberto, Phy. Rev. A34, 4971 (1986).