Resonant drift of autowave vortices in two dimensions and the effects of boundaries and inhomogeneities

Chaos, Solirons

& Fractab Vol. 5, ~os 314, pp. S-E-~?& 1995 Copyright 0 1995 Elsevicr Science Ltd Printed in Great Britain. All rights rcscrved c960-0779/95$9.50 + .oo

0960-0779(93)EOO44-C

Resonant Drift of Autowave Vortices in Two Dimensions Effects of Boundaries and Inhomogeneities

and the

V. N. BIKTASHEV Institute of Mathematical

Problems of Biology,

Pushchino, Moscow Region 142292, Russia

A. V. HOLDEN Department

of Physiology and Ccntre for Non-Linear

Studies, University

of Leeds, Leeds LS2 9JT,

UK

Abstract-The effect on autowave vortices in R* of a non-localized, periodic external forcing with a near-resonant frequency is investigated, and the effects of boundaries and inhomogeneities in the medium and mutual interaction between vortices is explored. We adopt an asymptotic approach which results in a phenomenological description of the resonant drift in terms of ODE’s for the location of the vortex core and its rotation phase. We compare the results predicted by the theory with numerical simulations. The behaviour of the vortex is more complicated than the resonant drift described by Davydov et al. [V. A. Davydov, A. S. Mikhailov and P. K. Brazhnik, “Kinetica i gorenie”, p. 39, Institute of Chemical Physics, Chernogolovka (1986) (in Russian); V. A. Davydov, V. S. Zykov, A. S. Mikhailov and P. K. Brazhnik, Drift and resonance of spiral waves in distributed media, fzv. Vuzov-Rad@$izika 31, 547-582 (1988)] and Agladze el al. [K. I. Agladze, V. A. Davydov and A. S. Mikhailov, An observation of resonance of spiral waves in distributed excitable medium, Pis’ma v ZETP 45, 601-603 (1987)]. In particular, a resonantly drifting vortex can be repelled from a boundary (phenomenon of ‘resonant repulsion’). Possible applications of resonant drift in the control of cardiac arrhythmias are. studied. In particular, an approach for overcoming the resonant repulsion is proposed, which involves a feed-back from the medium to the stimulation device. We show that with this technique, the voltage needed for extinguishing vortices can be an order of magnitude less than that needed for the traditional single-pulse defibrillation.

1. INTRODUCTION

In this paper we present a method of controlling the location of an autowave vortex in an excitable medium, using resonant drift. We proceed from considering the symmetry group of the idealized mathematical problem towards proposing a new technology for treating cardiac arrhythmias. Cardiac arrhythmias are any disturbances in the normal cardiac rhythm [4]. Paroxysmal atria1 tachycardia and ventricular fibrillation are dangerous and potentially lethal arrhythmias, and are both oscillatory in nature. There are various hypotheses for their physiological origins, and a widely accepted hypothesis, deriving from Mines [5] and Garey [6], is that their mechanism lies in recirculation (re-entry) of excitation, initiated via a temporary propagation block. This mechanism depends on the assumption that propagation of excitation in cardiac tissue may be explained as the propagation of waves in an excitable medium. 575

576

V. N. BIKTASHEV

and A. V. HOLDEN

An excitable medium responds to a small local perturbation by a local, decremental response, and to a sufficiently large and extensive perturbation by a non-decremental propagating strongly non-linear wave that differs from a soliton in that it is annihilated by collision. Excitable media may be modelled by reaction diffusion equations, and are realized in a variety of physical, chemical and biological systems [7, 81. The spread of electrical excitation through heart muscle may be modelled by non-linear wave activity in an excitable medium [9]. An example of a chemical excitable medium is the BelousovZhabotinsky reaction and its variations [lo]. Wiener and Rosenblueth [ll], using a simple, axiomatic model of an excitable medium, have shown that re-entry may take place in the absence of anatomical “holes”, i.e. recirculation is around a “functional” block which is formed by the refractory tail of the same waves of excitation; and such recirculation can provide an arrhythmia with maximal frequency. This kind of re-entry has been demonstrated in cardiac tissue by Alessie et al. [12]; a recent review on this topic can be found in [13]. It is important to note that re-entry around a “functional” block (also called a leading circle, or reverberator) may take place in a completely homogeneous medium. In non-linear science, this kind of re-entry is called an autowave vortex, rotor, or spiral wave, as in a large enough medium it is a source of a rotating wave having a spiral form. Such spiral waves have been intensively studied in chemical systems, e.g. in the Belousov-Zhabotinsky reaction medium [lo, 14, 151. The existence of spiral wave solutions of several particular reaction-diffusion models has been shown both analytically and numerically (see reviews in [16-181). Recently, the existence of vortex solutions was demonstrated for the Beeler-Reuter model [19] for the electrical activity of ventricular muscle [20-221. Autowave vortices in simple mathematical models have some symmetry properties (see Fig. 1). If the medium is unbounded and has properties constant in time and homogeneous

Fig. 1. The autowave vortex of the PDE model (see Appendix). (a-e) Successive positions of the excited region every 45 time steps. The period of the vortex is 280 time steps, the wavelength is about l/2 of the medium size, the core size is of the order of l/10 of the medium size. The vortex rotation is approximately rigid, except close to boundaries. (f) Phase plane of the diffusionless system, with null-clines. The shades of the phase plane code the states of points of the medium in (a-e), e.g. not shaded at (a-e) are points where u < 0.5 and v < 0.5.

Resonant drift of autowave vortices

577

in space, the vortex solution can be periodic in time, and, moreover, be “rigidly rotating”, i.e. appear stationary in a properly rotating frame of reference. Though this property is only approximately valid for real vortices in the BZ reaction and, because of anisotropy and inhomogenities in cardiac muscle, even less valid for cardiac tissue, it allows us to propose a practical method.for the control of re-entrant vortices in heart muscle. As re-entry provides the mechanism for some dangerous cardiac arrhythmias, it is important to prevent the onset of re-entry or, once it is established, to extinguish it rapidly. The first problem is approached by pharmacological means, and the second one needs more urgent action, as death may be imminent. Electrostimulation is an appropriate method. There are various techniques of electrostimulation, including local stimulation to stop re-entrant activity, which may be effective in case of a single re-entry path, or defibrillation (cardioversion), which is the only effective tool to annihilate most highfrequency tachycardias and fibrillation. In defibrillation a high amplitude electric shock is applied throughout heart tissue, i.e. it is not localized. This method forces the heart into a uniform resting state, out of which normal pacemaker driven activity can emerge, and so does not depend on the mechanisms responsible for the pathological re-entrant activity. However the current densities needed can result in tissue damage. There is a real need for less destructive methods of cardiac defibrillation that, by exploiting specific features of re-entrant mechanisms, use low voltage stimulation. The above explanation of the action of the defibrillating shock is not exhaustive, as it does not explain some known experimental phenomena. For instance, if we apply two appropriately timed pulses instead of one, and the amplitude of the first pulse is below the defibrillation threshold, then the amplitude of the second pulse required for defibrillation is lower than that for single pulse [23]. The effect of low-voltage non-localized stimulation on autowave vortices was analysed in [24, 251. It was shown in numerical experiments for simple mathematical models of excitable medium, that a generic result of such a stimulus is the space shift of the functional block (core, rotation centre), whose magnitude is determined by the amplitude of the stimulus and whose direction by the time of its delivery (and so the direction of the virtual shift) rotates synchronously with the vortex. This suggests another interpretation of defibrillation: if the shift of the vortex throws it out of the excitable medium, re-entry ceases. This interpretation means that if the vortex is located asymmetrically in the medium, then the defibrillation threshold can vary several fold, and is dependent on the time of delivery of the defibrillating shock [24, 251. Defibrillation by a single shock does not take into account the oscillatory nature of re-entry and a possible way to decrease the defibrillation voltage would be to use an oscillatory signal. This was attempted by Krinsky et al., and is related to the phenomenon of induced drift [26, 271. In brief, if an excitable medium supporting an autowave vortex is stimulated with a spatially uniform periodic signal with a period shorter than the vortex period but longer than the refractory period, then this stimulation will synchronize the vortex and induce it to drift. In fact the stimulation transforms the vortex into a wave break in the field of waves generated by the spatially uniform stimulation. A problem is the choice of the range of periods, as it is very narrow in cardiac tissue (and some physiologists believe it is absent); in addition, this range is not constant and will vary with time and location in the medium. These problems can be overcome in principle by scanning the period of the stimulation. This has been demonstrated in the BZ reaction medium [26], but we are not aware of any experimental application of this method for cardiac tissue. A resonant influence on a vortex is clearly described by Davydov et al. [l, 21 as the phenomenon of parametric drift. A spatially homogeneous, time varying influence on the medium is mathematically equivalent to a time variation in the medium parameters, hence the term “parametric drift”. This phenomenon was first found in a kinematic model,

578

V. N. BIKTASHEV

and A. V. HOLDEN

derived asymptotically for wavefront motion, with the assumption that the period between waves is long, so that one wave is not influenced by the propagation of the preceding wave [28, 291. Such autowaves can occur in weakly excitable medium, where the spiral has a large core [30-321. For special parameter values, the motion of the vortex tip in such a model can be described by a closed system of ordinary differential equations, and this motion in a time invariant, spatially homogeneous medium is along a circle, the boundary of the core [2, 16, 291. If the medium parameters are changing in time and/or space, the equation has more complex solutions. Davydov et al. considered the cases of a linear parameter gradient in space and periodic sinusoidal oscillations in time. The first case led to a theory of drift due to inhomogeneity, and the second revealed that the vortex tip will undergo a bi-periodic motion, with one frequency Q being that of the stimulation and the other frequency o the natural frequency of the vortex. The resulting trajectories are of the form of epicycloids, and are results of two circular motions, one with the core radius and the other with a radius R = c&2 - co)

(1) where the velocity c is determined by the perturbation amplitude (in their approximation, directly proportional to the amplitude). If Q = o (resonant case), the epicycloid degenerates into a cycloid, and the vortex motion averaged over its own period is along a straight line with the velocity c, the direction of the straight line being determined by the difference between the perturbation phase and vortex rotation phase. The kinematic approach used in [2, 16, 291 enabled the authors to consider analytically some other problems, such as evolution of vortices on curved 2D surfaces, in anisotropic 2D media and in 3D media [33-391; reviews of these results can be found in [16, 29, 401. It is important to note that the features characteristic of this phenomenon are related to the rotational symmetries of the processes involved and do not depend on specific properties of the autowave medium and the nature of the forcing. Though the assumptions used in the kinematic approach are not really applicable to real excitable media as the BZ reaction medium and heart tissue, the phenomenon of parametric drift has been observed in BZ reaction [3]. A natural idea is to apply this effect to eliminate vortices in the heart. Resonant drift may be used to drive the re-entrant vortex to the absorbing, inexcitable epi- or endocardial borders, where the vortex will annihilate. Since resonant drift can be produced with small amplitude forcing, this way may provide an important alternative to cardioversion. Here, as well as in the previous “resonant” effect of induced drift, the problem of frequency choice arises. If the stimulation frequency is not properly chosen, then the vortex will just meander (in theory a bi-periodic motion) within a small area, instead of a directed drift. In the experiments of Agladze et al. [3], this problem was solved by measurement of the natural frequency of the vortex. An alternative approach is by scanning through the range of possible frequencies. However, the inhomogeneity of cardiac tissue will make the resonant frequency depend on the location of the vortex. This complicates the vortex drift and aggravates the problem of frequency choice. Another way of choosing the frequency for “resonant” methods of vortex control is to use feedback. Garfinkel et al. [41], proposed using feedback in an approach that treated the heart as an abstract dynamic system, or a “black box”. This dynamical system was controlled using signals derived from recordings of the heart’s electric activity. Control of chaotic cardiac dynamics was demonstrated. Unfortunately, both the practical applicability, and mechanism of this technique are not clear. In this paper we analyze resonant drift as a possible means of eliminating re-entrant cardiac arrhythmias. First, we propose an asymptotic theory which shows that resonant drift is a generic property of autowave

Resonant drift of autowave vortices

579

vortices, and which describes this phenomenon in an inhomogeneous medium. Then we analyze some of the consequences of the theory, and illustrate the theoretical predictions using a simple FitzHugh-Nagumo excitable medium. The phenomenon of resonant repulsion of drifting vortex from medium boundaries, inhomogeneities and from each other is described and explained. We demonstrated that feed-back stimulation overcomes both the problem of frequency choice and resonant repulsion.

2. ASYMPTOTIC

THEORY

AND PHENOMENOLOGY STIMULATION

OF PERIODIC,

SPATIALLY

UNIFORM

The main ideas of the asymptotic theory we use here can be found e.g. in [42]; the specific form follows the approach used in [43] for another problem. The small parameters are the amplitude of the external influence, the inhomogeneity and the difference between the stimulation frequency and the natural frequency of the vortex; we assume these three small parameters to be of the same order. 2.1. Main assumptions We only consider generic reaction-diffusion can be written in the form: a,u = f(u)

systems on the plane r = (x, y) E R* that

+ DV2u + E~(u, r, t).

(2)

Here u = u(r, t) E I@ denotes an l-component vector (membrane potential, and gating variables), f E IF!’ is an l-component vector (transmembrane currents and gating variable rates) and D E R?’ is a constant 1 x 1 matrix of diffusion coefficients (only one of the diffusion coefficients, that corresponds to membrane potential, is non-zero in cardiac excitation systems). The term h E R’ denotes a perturbation, assumed to be uniformly bounded; and E is a small parameter. The symbol a always denotes partial derivative by the subscript variable, V is the space gradient 3, = (a,, a,). First, we assume that the system (2) with no stimulation (Eh = 0) supports rigidly rotating vortices, which are solutions of the form u = U(r, t) = U(p, 6 + cot),

(3) where p = p(r), 6 = 8(r) are polar coordinates in the (x, y) plane, and U is 2n-periodic with respect to the second argument. In addition, we assume this solution to be stable. Before formulating this assumption, let us note that every function Ug,6r(r, t) = U(r - br, t - f3t)

(4) for arbitrary constants 6r, at, is also a solution to (2), i.e. there is a three-dimensional manifold (family) of solutions, generated by the symmetry group of the system (2). The manifold is not four-dimensional as the symmetry group, because the solution (3) is invariant under simultaneous time translation by dt and space rotation by 66 = -WC?.Thus, the strongest stability requirement possible in this situation, is: Stability postulate: Any solution to (2) at E = 0 with initial conditions sufficiently close to (3) tends to one of the solutions (4) with some 6r, 6t as t + m, i.e. the family (4) is stable as a whole against small perturbations in initial conditions. The precise meaning of this postulate depends on the interpretation of “sufficiently

V. N. BIKTASHEV

580

and A. V. HOLDEN

close”. We believe that real autowave vortices are stable against perturbations small not only in an integral norm, say, L1, but also in some uniform norm, say, C or Ci. We now analyze the consequences of this postulate in the linear approximation for deviation from the family (4).

2.2. Stability to initial condition perturbations Let us consider the behaviour of the solutions to (2) at E = 0 from initial conditions u(r, 0) = U(r, 0) + uO(r)

(5) with small uo. Taking a linear approximation in uo, we can obtain solutions to (2) in the form u = U(p, 6 + wt) + u(r, t),

(6)

where the linear approximation in u requires 3,~ = Lv = DV*u + F(U(p(r),

8(r) + wt))u

(7)

with initial conditions u(r, 0) = uo(r).

(8) Here F(U) = af(u)l,=cl is the matrix of partial derivatives of the reaction terms taken at the solution (3). The analysis of stability of (3) then reduces to the analysis of the properties of the solutions to (7), (8). This analysis is performed more conveniently in a rotating coordinate system,

t-+i=t,

r 9 f:

p(r) = p(i’),

19(?)= 8((r) + wt

(9)

where the linearized equation becomes aiu = ZU = DP20 - O~&U+ F(U(I’))u

(10) where ?* is the Laplacian operator in the rotating frame of reference. Note that the operator 1 in (10) does not depend on time i, and the problem of linear stability is reduced to studying the spectrum of this operator. The Stability postulate of Section 2.1 requires that z has no eigenvalues with positive real parts. Due to the symmetry, there are at least three eigenvalues at the imaginary axis, related to neutral stability to space and time shifts. Indeed, any solution of the family (4) differs from the solution (3) at small &, 6t by a function linearly composed of the three linearly independent functions, which in the rotating frame of references have the form Q. = -a,u

= -oaaU(i;),

VX = -a,U

= -(E3,U(T)cos(~ - 03) - p-l+U(T)sin(B

VY = 4,U

= -(a,U(r”)sin(B

- 2) + p-l~~U(?)cos(B

- wi)),

(11)

- ml)),

which, consequently, are the solutions to the linear equation (10) @his can easily be checked by direct calculations). This means that the linear operator L has at least the following three non-decaying eigenfunctions, with corresponding eigenvalues: vo = -wUa(f), &I = 0, Q, = -t( up - ijj-’ us) epiB,

Al = +iw,

V-1 = -i(UP + ip-‘U,) e+i9,

A-1 = -iw,

02)

581

Resonant drift of autowave vortices

and the solutions vX’,,vY are expressed as vx

= ~‘leioi

vy

= iv,

eioi

+ Qele-ioi, _ iv’_,

(13)

e-iwi.

According to the Stability postulate, there should be no other eigenvalues on the imaginary axis. In the linear approximation, the constants 6r, 6t defining the limit solution (4) correspond to some initial conditions uo, should linearly depend on uo. This can be written as u(i, T) = coVo(r”) + cleiwiVl(r”) + c-lenio’V’_l(?) + o(1)

(14) as ?+ ~0, where the coefficients are determined by some linear functionals taken on the initial conditions: CO,tl

=

(wO,*ico,

O(?,

O)).

(15)

It is easy to_check that the functionals (WA, .) in (15) are the eigenvectors to the adjoint operator to L: (i&,

Ll)

= (iv*, a,u) = a,(w$l)

= n
(16)

(The second and third identities follow from equations (14), (15) due to their linearity and invariance under change of fiducial time.) We choose these eigenvectors so that <%, qJ = 4%,.

(17)

If the functionals are written in an integral form via some regular functions pi’,, (t&,

a) = j-(i;,(2),

. ) d2F,

(18)

where <,> denotes the inner product in C’, then the functions F* obey the equations z+Yk = A*?~‘,,

(19)

where z’ is differential operator formally adjoint to z: z+ = DQ2 + cd8 + F+(U(I’)),

(20)

and F+ denotes transposed matrix to F. The conditions at p + 03 for the problem (19), (20) should be determined from the properties of the adjoint space. If the Stability postulate holds for perturbations that are bounded but not vanishing as p + 03, then the adjoint space must contain only functions integrable over the whole plane. So, the functions FL should be integrable, i.e. rapidly vanishing as p + 03. We shall refer to the characteristic distance of these functions decaying as the sensitivity distance.

2.3.

Stability to perturbation

of the right-hand sides; regular perturbation

technique

Now let us consider the problem (2) with non-zero but small E, and initial conditions equal to the unperturbed vortex with rotation center at the origin u(r, 0) = U(r, 0). The linearization on this vortex u(r, t) = U(r, t) + .w(r, t)

(21)

yields, in the rotating frame of reference, the problem +u = zu + h(U(T),

r”, 2).

(22)

V. N. BIKTASHEV

582

and A. V. HOLDEN

Let (GA, - ) be the eigenvector to the adjoint operator E+ for some eigenvalue A. Then the evolution of the corresponding component of u is determined by equation a;(&

u) = A(WA, u) + (tik, h)

(23) and this equation is explicitly solved. If the sets of @*, PA for all eigenvalues A constitute bases in the corresponding spaces, this yields a general solution to the linear perturbation problem (22): .

i

u(i) = xvAe*’ e-*‘(mk, h(t))dr. I0 A

(24)

However, independent of the assumption for the existence of these bases, we know from the previous section that there are eigenvalues A = 0, kio, for both z and z+. For the corresponding components, the equation (24) can yield unbounded growth for some h, even if h remains uniformly bounded in time e.g. (@o, u(i)) = j(@o, h(r))dr,

(25)

0

which tends to infinity, e.g. if (W,, h(t)) is a non-zero constant. This means that the solution to (2) drifts along the family (4), generically to large distances (“secular growth”), where the regular perturbation technique is no longer valid. 2.4. Stability to perturbation of the right-hand sides; singular perturbation technique To obtain solutions uniformly valid for large ts at any uniformly bounded perturbations h, let us apply the singular perturbation technique (see [43]; an analogous technique was used in [44]). Expand the solutions to (2) in the vicinity of the manifold (4) in the sum u(r, t) = U(r - R(t), = W(r

t - #(t)/o)

+ &U(T, t)

t9(r - R(t)) + cot - a(t))

- R(t)),

+ m(r, t)

(26)

where R(t) and Q(t) are the coordinates of the projection of solution u onto the manifold (4), R = (X, Y) is current rotation centre position, @ is current initial rotation phase, and u is orthogonal to the manifold (4): (wo,+l> u) = 0

(27)

where IV,,,, are the functionals @o,+l defined by (16), carried over to the non-rotating frame of reference with origin at-R(t) and fiducial time cP(t)/w. If the functionals @A are represented by regular functions Yn, then the explicit form of WAmay be written as (WA, -) = I( pA(p(r - R), 8(r - R) -t wt - a)), - ) d2r.

(28)

Naturally, since the projectors % were constant in the rotating frame of reference, they become rotating in the original one, as can be seen from (28). Supposing that the drifting vortex always remains in the c-vicinity of the manifold (4), putting (26) into (2) and leaving terms 0(c), we get d,u

= E[LU + h(U(r + a,U(r

- R, t - @p/o),r, t) + a@U(r - R, t - @/o)@‘(t)

- R, t - @/w)X’(t)

+ k3,,U(r - R, t - @/m)Y’(t))] + O(e2),

(29)

where the prime denotes time derivative, or, in the frame of reference associated with the vortex, t” = t, p(f) = P(r - R), q(T) = 8(r - R) + of -. Q,:

Resonant drift of autowave vortices

583

&E+V= &[ZV + h(U(f), r, t) - co-'V@)W(t) _ (fYl(i;) ei(+Q) + Q-,(t) e-f(w'-Q))xf(t) _ (ij,Yl(qei(+Q)

- iv-,(qe-i(@'-Q))yf(t)]

+ O(E2).

(30)

Now, applying the orthogonality conditions (27)-(30), and using the biorthogonality of eigenfunctions and adjoint eigenfunctions corresponding to different eigenvalues, we get the “motion equations” that determine the slow evolution of the vortex along the manifold (4): a,@ = &O(miO, h) + 0(&2), 3,X = cRe(e-‘(“‘-Q)(@:,,

h)) + S(E*),

(31)

3,X = eIm(e- i(wt-Q)(wl, h)) + 6(&2).

Passing over to the original frame of reference and using the integral expression for the adjoint eigenvectors, we get a,@ = E (p,-,(p(r - R), fir

I

3,~ = ewe

e-i(wt-Q)

( a,y = elm

( yl((p(r

- R) + ot - Q), h) d*r + S(E*), - R), f?(r - R) + tot - a), h) d*r

1 e-i(@t-Q)

(

( pl(p(r

+ O(E*), 1

- R), 2p(r - R) + wt - a), h) d*r

(32)

+ S(E*) 1

1

These expressions describe the evolution of the three vortex coordinates @‘, X, Y under the generic perturbation h(u, r, t). 2.5. Motion equations in the simplest case Consider the special case of h(u, r, t) = H(Qt - cp>H(q, + 2~) = H(q),

(33)

i.e. of a spatially uniform, periodic stimulation with a frequency S2,which we assume to be close to the natural frequency w of the vortex. In this case, averaging the motion of vortex over the period of stimulation, we get the following equation for the average coordinates: a,G = &HO + q&2), 3,X = ~IH~lcos((ot

- S) - (S2t - (p) + argH,)

+ O(E*),

a,P = clHIJsin((wt

- G) - (Qt - q) + arg H,) + a(&*).

(34)

Here Ho, H1 are the Fourier-components of the projections of the stimulation on the slow manifold, H,, =

p,, d*r, H(q) e-‘“‘7 dy, n = 0, 1. (35) i i(i If we denote the current phase difference between the vortex rotation and the stimulation by $9 $=(wts)-(Qt-q?)+argH,, (36)

then (34) can be rewritten as

V. N. BIKTASHEV

584

and A. V. HOLDEN

a,@= w + &HO - 52 + S(2), 3,x = &pzllCOSt$+ 0(&q,

(37)

El,F = EJH1lsin f$ +. 0(S). These motion equations are obtained for the case when the stimulation is described by a space-independent additive term in r.h.s. of (2). They describe a circular motion of the vortex centre (X, Y) with a velocity ~1~~1 and angular frequency o - B + MO; in the resonant case SZ= w + &HO the circular motion degenerates to a straight drift, whose direction is determined by the initial phase difference @ between vortex rotation and external influence. Note here that the resonant frequency does not coincide with the natural frequency of the vortex, but differs from it by a value E& proportional to the amplitude of the external influence. In a more general case, the external influence is described by variations of some of the parameter(s) a of, say, the function f, which does not depend explicitly on r: i3,u = f(u, a(t)) + DV’u;

a! = eA(Qt

- q),

A(q, + 24 = A(q).

(38)

Linearizing in E, we obtain an equation of the form (2), with the additive term being h(u, r, t) = a,f(u, O)A(Qt - ~4

(39) which leads again to the averaged motion equations (37) with the parameters Ho, Hi defined by EM

= /$%d~, EL Lf(W9 5>, 0)) dWwWWinvdrl. (40) f This case corresponds to the parametric drift used in [l-3], and we see that this generic case presents no new features that are not found in the special case (33), which is relevant for applications in cardiology.

2.6.

Inhomogeneous

medium and influence of boundaries

The asymptotic technique described above can be applied to an inhomogeneous medium, i.e. when the medium parameters explicitly depend on the spatial coordinates, e.g. 3,~ = f(u,

a(t), P(r)) + DV*u.

(41) Here the parameter a! is, as before, the stimulation, and p represents inhomogeneity of the medium. We assume that both these parameters are small, i.e. &(I) = eA(Qt - q) and P(r) = d(r), and the stimulation is periodic, A(q, + 2n) = A( Q$. This problem reduces to the one considered above, via linearizing in E, to give (2) with h(u, r, t) = a,f(u,

0, O)A(t)

+ 'Spf(u, 0, O)B(r).

(42) The resulting “forces” acting on the vortex then can be calculated by the general formula (32), which yields the expression for the “forces” in the form n = 0, 1, fL{~Y PI = K{a) + HdBl~ where H,{ a} are given by (40) and Z!&(p) are defined as K{Pl

= #pR(P?

5-L +?f(wb

(43)

699 0))

B(X + pcos(E - q), Y + psin(E - q))dEpdpe-‘“qdq.

(44

Resonant drift of autowave vortices

585

Note, that the effect of inhomogeneity is always resonant, since it acts via a perturbation of the vortex which is naturally synchronized to the vortex rotation. This means that the values of the terms HO,,(~) are completely determined by tire distribution of the inhomogeneity parameter B(r), i.e. given B(r), the formulae (44) yield the functions Ho,IW@). The convergence of the integrals in (39, (40), (44) is an immediate consequence of the Stability postulate, i.e. if the autowave vortex is stable to small perturbations of initial conditions in location and phase shift then the right-hand sides of the equations (37) can be considered as linear operators on the parameters A, B. The origin of the terms H,{#I}(R) and Ho{/3}(R) as integrals of the product of variations caused by inhomogeneities and the sensitivity functions means that if the formulae are applied to smooth inhomogeneities, then C,(R) and C,(R) are proportional to the parameter gradients at the point R, which is in agreement with the results of phenomenological considerations and numerical simulations of [45]. For non-smooth, e.g. stepwise inhomogeneities, (44) results in a smoothed dependence of these functions on the spatial coordinates. For instance, a parameter step will be smeared over the width of the sensitivity distance. This is in qualitative agreement with numerical experiments of [46], where the vortex approaching a stepwise inhomogeneity is influenced by this step over a finite distance. The resulting equations of motion can be written in the form W =

Q,(R)

- Q + CW*V)arg(Hl{4),

a,X = C,(R) + c(R) cos $,

(45)

i3,Y = C,,(R) + c(R) sin 4. (Here and below, we omit dash signs over R, X, Y and the terms O(c2).) The term B,(R) is the local resonant frequency, Q,(R) = o@(R) + &H,(a)

= w + EH~{(u} + cH,{/3}(R),

(46)

c(R) is resonant drift velocity, c(R) = +WNR)I and C,, C,, are the velocities of the drift caused by inhomogeneities C,(R) + I’C,(R) = EHIUXR).

(47)

(48) The “anisotropic” term (S,R * 0) arg ( HI{ (u}) in the first of equations (45) does not appear to be essential for our modelling of PDE simulations, and is omitted below. System (45) presents the superposition principle for the “forces” of different origins producing drift of a vortex. The validity of a superposition principle for the effects of homogeneous stimulation and spatial inhomogeneities encourages us to propose that the same principle applies to the influence of boundaries. Formally, boundaries can be considered as a special, very strong form of inhomogeneity; and so there are no reasons to treat boundary effects with the same formalism as that for small inhomogeneities. However, numerical experiments show that the dynamics of vortices close to boundaries are analogous to those in inhomogeneous media [47]. Both the direction and velocity of drift, as well as deviation of its frequency, are determined by position of the vortex with respect to boundaries, i.e. can be described by the same terms Q,(R), C,.(R), C,(R) in (45). So, for the rest of the paper, we use the system (45) to describe resonant drift for both an inhomogeneous and a bounded medium. Note, that the system (45) (without the

V. N. BIKTASHEV

586

and A. V. HOLDEN

@,R*V)arg(f&W)) was obtained from purely phenomenologic considerations [48, 49’). The above formalism, however, presents the principal way of calculating the explicit form of the terms in (45). Though the problem of finding the projectors WA can be solved only numerically, the formulae (35), (40), (44) can give a way of establishing qualitative relationships between various phenomena. For instance, from the form of (44) we can conclude that the distance at which the vortex responds to the inhomogeneity is the sensitivity distance, i.e. a parameter defined for an unperturbed problem, and is the same for most types of inhomogeneities. 3. PROPERTIES

OF THE CONSTANT FREQUENCY RESONANT WITH NUMERIC EXPERIMENTS

DRIFT AND COMPARISON

In this section we consider some partial solutions of the system (45) and present numerical results that are in qualitative agreement. The numerical PDE model is presented in an Appendix. To get qualitative and, possibly, quantitative agreement with numerical results, we have chosen special values of parameters in the phenomenological system (45), which are also present in the Appendix. This system will be referred to as the ODE model. 3.1.

“Larmor” drift

First, let us analyze the system (45) if C, and C, can be neglected with respect to c and if c(x, y) = const. Excluding #, we get the system d*X/dt* = -&2(X,

Y)dY/dt,

d*Y/dt* = 652(X, Y)dX/dt, 651(R) = t&,(R)

(49)

- SJ,

which coincides with equations of motion for a classical electrically charged particle in a magnetic field perpendicular to its plane of motion; here the frequency deviation SQ is analogous to the magnetic field magnitude. At &2(R) = const. i.e. in the absence of inhomogeneity and interaction with boundaries, we get a drift motion analogous to the Larmor rotation of electrodynamics [50], i.e. the circular drift as described in Davydov et al. [2]. The behaviour of the vortex at a near-resonant frequency of forcing is illustrated in Fig. 2. If the trajectory is always far from medium boundaries, so their influence is weak, then the trajectory of the vortex core is close to a circle.(Fig. 2a), and we observe the analogue of Larmor drift. The ODE model (Fig. 2b) confirms that the deviation of the trajectory from a circle can be explained by a weak influence of the boundaries and slight asymmetry of the initial position of the vortex. The dependence of the resonant frequency Q, and the velocity c of resonant drift on the stimulation amplitude for our PDE model are shown in Fig. 3. It can be seen, that for a large range of biophysically reasonable amplitudes both dependencies are approximately linear. 3.2. Exponential interaction with boundary: resonant repulsion Since we are interested in applications to the control of re-entry in cardiac tissue, we now suppose that 6Q(X, Y) # const., and is caused by interaction with the boundaries. Here we show that, with some assumptions on the form of this interaction, the drifting vortex will be repelled from boundaries, and this repulsion is primarily due to the

Resonant drift of autowave vortices

587

dependence of the natural frequency of the vortex on the distance to the boundary. A numerical illustration of this phenomenon is presented on Fig. 4. For modelling this effect with our ODE model, assume that the boundary is the line X = 0, and the vortex is located in the left-hand part of the plane, for X < 0. Let us assume that the influence of the boundary decays exponentially with distance (this is shown for a special case of autowave medium in [51]). Suppose, that far from the boundary an exact resonance takes place. To estimate the influence of the boundary, put &2(X, Y) = aexp(-(X(/R)

(50) where R is the sensitivity distance of the vortex. With these assumptions system (45) may be solved explicitly. Eliminating X from (45), (50), we get R d2@/dt2 = c cos @d@/dt,

(51)

from where - sin $1~0). ct = R d#&sin C$J (52) I Here we suppose that the vortex approaches the boundary from infinity at the angle &, = c$(t = - ~0). We can see that if sin (c#J)+ sin (&,), then t + f ~0. This means the vortex is repelled from the boundary, the angle of reflection being equal to the angle of incidence: $$ = #(t

= +m)

= ?T - q$).

Fig. 2. (a).

(53)

V. N. BIKTASHEV

588

and A. V. HOLDEN

Fig. 2. Near-resonant drift of the vortex in the cylindrical medium. Trajectory of the vortex rotation center is shown; the direction of drift is signed by arrows. (a) PDE simulation, forcing amplitude 0.025 and period 272 time steps while resonant period at this amplitude is 269 time steps. (b) ODE model, c = -0.0005. The trajectory does not reach the vicinity of boundaries and is close to a circle.

Now let us take into account also the term C,. For the special case considered in reference [51] it also decays with the distance, and the characteristic space lengths for $2, and C, coincide, so we put C,(X, Y) = -bexp(-1X1/R).

(54)

We again observe reflection of the vortex from boundary, but not with the angle of reflection equal to the angle of incidence. Now we have a coupled system for X, @: dX/dt = ccos($) - bexp(X/R); (55)

d$/dt = u exp (X/R), with initial conditions X(t =

-m)

=

f#(t = -m) = A.

-co,

(56)

The corresponding $I~(&,) for these systems can be obtained from the finite equation SiIl(C#+

+

0)

-

Cexp(-&tan(@)

= 0,

i = 0, 1,

tan (0) = b/(aR),

(57)

i.e. the reflection angle @rcan be found by the following procedure: first find a constant C such that @,,obeys (57), then find the next root of (57), to the right if a > 0, and it will provide $r (see Fig. 5, the details of the derivation can be found in [48]). According to this rule, if 8 = 0 (b = C, =‘O, Fig. 5(b)) we obtain the “mirror” reflection, as in (53). For

Resonant drift of autowave vortices 1.8 1.7

-

1.6

-

1.5

-

1.4

-

1.3

-

1.1

-

1.0

-

0.9

-

0.8

-

0.7

-

0.6

-

0.5

-

0

589

I

I

I

I

I

I

I

I

I

1

2

3

4

5

6

7

8

9

10

Fig. 3. The parameters of resonant drift at various forcing amplitudes. They were chosen to achieve a straight drift in the cylindrical medium parallel to the y-axis far from boundaries at some initial phase. The boundary conditions at upper and lower boundaries are periodic. Abscissa: stimulation amplitude A * 100. Circles: (reciprocal resonant frequency minus that of free vortex) * 100, (time steps)-‘. Triangles: resonant drift velocity * 100, space steps/time step. The single pulse defibrillation threshold - which physically is not small - is about 0.2, which corresponds to 20 in this display.

8 # 0 the character of the reflection-incidence mapping can be seen in Figs 5(a, c).%In particular, (i) for 8 < 0 (the free vortex is attracted to the boundary, Fig. 5(c)) there is a range of angles of incidence, that depends on medium properties (b, a, R) and not on c, over which there is no repulsion and the trajectory of the system (45) goes to x + +w in a finite time; (ii) if 8 > 0 is sufficiently non-zero (Fig. 5(a)), then the reflection angles for a wide range of angles of incidence are located in the vicinity of the angle n - 0. We use this qualitative observation below. By varying the initial phase of stimulation in numerical simulations, we measured the dependence of the reflection angle on the incidence angle. The result is shown in Fig. 6(a). In spite of experimental errors related to the difficulty of measuring the angles for the experimentally obtained vortex trajectories, it can be seen that the main qualitative properties of this dependence (see Fig. 5) take place. Namely, the dependence is monotonic and there is a definite plateau over a wide range of incidence angles. The height of this plateau (1.0 rad) enables the estimation of the dimensionless parameter 8 defined in (57) as about 0.54 radians. The dependence for the phenomenological model is presented in Fig. 6(b). It can be seen that within the limits of experimental errors, there is a satisfactory agreement between Fig. 6(a) and (b). 3.3.

Cycloidal drift

Intuitively, the occurrence of resonant repulsion may be explained as follows. The closer the vortex is to the boundary, the faster the changes in direction of the drift. If this is so,

V. N. BIKTASHEV

and A. V. HOLDEN

Fig. 4. A PDE solution illustrating the mechanisms of resonant drift and resonant repulsion. Stimulation amplitude A = 0.055, period 251 time steps. Successive positions of vortex in time intervals exactly 3 periods of stimulation. The trajectory of the core is shown by black dots at each picture. (a-c) The stimuli occur in the same rotation phase, and the trajectory is straight. (c, d) The natural rotation frequency of the vortex increases near the boundary, each successive stimulus occurring at later rotation phase, and the direction of the drift turns. (d-f) The vortex goes away from the boundary, it resumes its original natural frequency, and the trajectory is again straight.

the direction will change unless the vortex starts drifting away from the boundary. Sufficiently far from the boundary the direction stops changing. and the final direction of drift is away from the boundary, i.e. the vortex is repelled. This explanation suggests that the requirement for the strong resonance (X2 + 0,

591

Resonant drift of autowave vortices

1

I

Fig. 5. (a) and (b).

592

V. N. BIKTASHEV

and A. V. HOLDEN

Fig. 5. Maps of the angle of incidence and reflection (see equation (57)). The arcs of the exponential curves between intersections with the sinusoid join corresponding angles of incidence (#so) and reflection (4~~). (a) 13= 0.5; the interval (-x/2, $2) maps onto the interval (n/2, &), a wide range of angles of incidence maps onto the vicinity of the point $ = rr - 0. (b) 0 = 0; there is the “mirror” reflection; r#i = tr - A. (c) 0 = -0.5; not all the angles of incidence have corresponding angles of reflection. The angles of incidence in the region (-x/2, c#J*) appear to lead to attachment or annihilation of the vortex at the boundary, i.e. the angle $J accelerates exponentially to infinity; the angles of incidence from the interval (&, x/2) map into the interval ($2, 3x/2).

X + -w) is not necessary, and repulsion will also take place close to the resonance. In this case, the vortex drifts along circles with a large radius when far from boundary and with smaller radii when close to the boundary. When C, is negligible, and C,,, &2 depend only on x, then the ODE system (45) has a first integral J(tp, x) = c sin Q,- A(x)

(58)

dA = &2(x). dx

(59)

where

In this case, all the loops of the trajectory are similar, and the whole trajectory is of the cycloidal type. This cycloid can easily be obtained for ODE models at C, = 0. The quasi-cycloidal drift obtained both for PDE and ODE models is illustrated on Fig. 7. Because our PDE model has C, # 0, the cycloid is not perfect, and its radii differ for different loops.

Resonant drift of autowave vortices (4

80 1

I

I

I

.

20

I

.

120

I

I

I

I

60

80

100

I

I

20

40

. .

.

100

60

I

I

.

80

40

0

80

I

.

60 -

(b)

I

593

140

160

I

I

140

160

180

60

40

20

0 L 0

I;!0

180

Fig. 6. The dependence of angle of reflection p1 = $1 - rr on the angle of incidence pc = rr/2 - & (degrees). (a) same as in Fig. 2, but with a stimulation period 269 time steps. (b) ODE model. In spite PDE model, parameters I of great numerical fluctuations in (a), the qualitative accordance can be seen. In particular, there is a plateau in the range (80°, 180°) of the height about 60° for both dependencies, which is consistent with the analytical theory (see Fig. S(a)).

3.4.

Limit cycle in a cylindrical

medium

The long-time dynamics of the phenomenological system (45) in a cylindrical medium with exponential interaction with boundaries (50), (54) has a set of limit cycles that are neutrally stable to shifts along Y due to the symmetry of this medium. For distances between the left and right boundaries sufficiently larger than R, the inclination of straight pieces of trajectories is determined by the solutions of the equations c#I~(c#I,,) = & + rr and are close to the plateau height on the graph Fig. 6. On Fig. 8(b) a typical trajectory approaching a limit cycle is shown. The long-time behaviour in PDE simulations also shows something like a limit cycle (see Fig. 8(a)). But here the trajectory is not closed. An obvious explanation of this imperfection is the neutral stability against vertical shifts caused by the symmetry of the problem. Due to this neutrality, slight computational errors caused by the discreteness of

V. N. BIKTASHEV

594

and A. V. HOLDEN

the model medium are not damped, which causes a “stochastic” shift of the trajectory up or down at each next cycle. 3.5.

Stochastic drift

It is easy to see that since round-off errors occur at every time step of computations in a numerical experiment, their influence on the drift trajectory shape is greater with slower drift, i.e. with smaller amplitude of stimulation. This statement is in agreement with the Fig. 9(a), where a typical trajectory of the drift is shown at a lower stimulation amplitude. The trajectory sometimes even curves far from the boundaries. The hypothesis of the stochastic character of the drift at low stimulation amplitudes was checked with the help of the phenomenological model, supplemented with additive Gaussian white noise in the equation for 4 in (45). In Fig. 9(b) one realization of the stochastic process described by such a model is shown. This “stochastic” drift may have a practical interest, as it might simulate the cellular inhomogeneity of cardiac tissue. It is important to note that the repulsion from the boundary still takes place in this case. 3.6.

Resonant repulsion from an inhomogeneity

Variability of the natural frequency of the vortex can be caused not only by interaction with boundary, but also by an inhomogeneity of medium properties or of the stimulation

Fig. 7. (a).

Resonant drift of autowave vortices

595

(b)

Fig. 7. Cycloidal drift of the vortex near the resonance. Parameters are the same as in Fig. 2, with another initial position of the vortex. (a) phenomenological ODE model. (b) PDE simulation. The trajectory is of a large radius when far from boundary and of a small radius when close; it resembles a cycloid. Note that this is the trajectory of the averaged vortex rotation center; the trajectory of the vortex tip is wound around this one.

amplitude. Both these cases are of interest for cardiological applications, because heart tissue is inhomogeneous and it is impossible to provide perfectly homogeneous stimulation. Figure 10 shows typical trajectories for resonant drift in a medium with stepwise inhomogeneity in the parameter & and in a homogeneous medium with stepwise inhomogeneous stimulation. These pictures clearly illustrate the resonant repulsion. The resonantly drifting vortex moves along a circle of a large radius in one part of the medium, and then reaches another part where the resonance no longer occurs. It now moves along a circle of a small radius until it returns to the region of the medium where it is resonant. This effect is in addition to the resonant repulsion from boundaries, which can produce exotic trajectories. Note, that in this case, the vortex may not return to the “resonant” part of the medium, but will be attracted into the non-resonant part. This can be explained by the competition between the resonant repulsion and “spatial” attraction represented by the terms C,,,(R) in the phenomenological theory. The effect of resonant repulsion in an inhomogeneously stimulated medium is provided directly by the dependence of the resonant frequency on the stimulation amplitude, since the medium itself was homogeneous in this numerical simulation. This example shows that it is incorrect to consider the resonant frequency as equal to the frequency of the free vortex rotation.

596

V. N. BIKTASHEV

and A. V. HOLDEN

Fig. 8. The limit cycle of the vortex drifting between two boundaries of cylindrical medium. (a) PDE simulation. (b) the phenomenological model, with all parameters as in Fig. 6.

Resonant drift of autowave vortices

3.7.

597

Resonant repulsion of vortices from each other

When there are two or more vortices in a medium, their rotation frequencies depend on their mutual location, if they are sufficiently close to each other. This means that the resonant repulsion should take place also between the vortices. The properties of this repulsion should be analogous to that from boundary: it is easily seen, for instance, that the motion of a vortex near a straight impermeable boundary is equivalent to the motion of two oppositely rotating and symmetrically located vortices. The resonant repulsion of vortices from each other is illustrated on Fig. 11, where evolution of 8 vortices is shown during 10 pulses of stimulation. It can be seen that some of the vortices annihilate on boundaries or with each other, while others repel and survive. The figure only presents casesof repulsion of oppositely charged vortices; further evolution of the same pattern also shows cases of repulsion of equally charged vortices.

3.8.

General comments

The numerical examples illustrate that resonant repulsion is not restricted to the case of exponential interaction with the boundary, as presented in Section 3.2. In [48], we have shown in the phenomenological theory repulsion of the vortex from the inhomogeneity/ boundary takes place in rather generic conditions:

(4 + .

.

.’

:

. c. ‘+mm~mmmm.. . . . . :

.

. .

n . . .

.

.

.

.

8 l

.

=. .

.

-5 .

.m. =..‘E : .

: I

. . -m.

. mm..: +

. f .=

Fig. 9. (a).

. : . ;

; : .

598

V. N. BIKTASHEV

and A. V. HOLDEN

Fig. 9. Stochastic drift of the vortex with low amplitude of external forcing. (a) From the direct experiments, forcing amplitude 0.004 and period 280 time steps, equal to the period of the free vortex. (b) A realization of the stochastic process, generated by the phenomenological model supplemented with random noise. Calculations made by Euler method with step 200, the noise is simulated by adding Gaussian random value of mean amplitude 0.002 in the equation for a,~$at each integration step.

Theorem. Assume that the functions w(X, Y), C,(X, Y) (1) are independent of Y; (2) are continuous and w > 0, C, < 0 for any X; and (3) may tend to zero as X+ -03, but not as X+ +m. Then any solution X(t), G(t) of (45) that exists for t E (-m, + @J)obeys: X(t + +m) + -m; r$(t + +m) + const. The proof presented in [48] is based on the analysis of the evolution of the function J defined by (58), (59). This formal statement needs some remarks, about cases when the repulsion may not be seen because of being overcome by other factors. Remark 1. The theorem conditions are violated if the trajectory does not exist for t E

(-00, + m). This may take place for some initial conditions if the boundary is attracting, as at 6 < 0 for the exponential interaction (50), (54) when the solution of (45) is defined

it is

Resonant drift of autowave vortices

599

for t < t,, < torn. The physical interpretation of this restriction is that if the boundary is attracting, then the attraction may be greater than the resonant repulsion. Examples of this sort can be seen for the trajectories on Fig. 10, which did not return to the resonant area, but approached a limit cycle within the non-resonant one. Remark 2. It is often the case that a boundary which is attracting at large distances is repulsive at small ones, and the vortex “attaches” to the boundary [47] in the absence of stimulation. Remark 1 is not relevant in such a case. Nevertheless, resonant repulsion may be exceeded by attraction - the conditions under which the theorem holds are violated, as C, changes sign. Recall that the applicability of the phenomenological system (45) is restricted, in particular, by the assumption that interaction with the boundary or inhomogeneity is sufficiently small for the superposition principle to hold. Therefore the system becomes invalid in a sufficiently close vicinity of the boundary. However, the resonant repulsion requires detectable effect of the interaction between the vortex and the boundary. The question arises under what conditions the assumptions for (45) will be applicable. It can be demonstrated that if the resonant repulsion in (45) takes place and the amplitude of the forcing (the coefficient c in (45)) is small enough, then the vortex will not reach any prescribed vicinity of boundary, and so the conditions of applicability of

-m-m-...

‘.

.

‘a

‘.

‘9

‘a \

-0.

‘0 \m+-.&

l\

l\

i 7 : : .' ,*/

t

B ,.$ /’ 4 .* 1

: I

i L ,i ./=

l\

l\

t I

600

V. N. BIKTASHEV

and A. V. HOLDEN

Fig. 10. Resonant drift in inhomogeneous media. Dotted lines show the boundaries of homogeneous parts of the medium, each medium consists of three homogeneous pieces. The direction of drift is shown by arrows; each dot on trajectory corresponds to one period of rotation. (a) Inhomogeneous stimulation amplitude: 0.03125 in the middle part and 0.055 in the lateral parts; the stimulation period 256 time steps is close to the resonant one for the lateral parts. Two trajectories corresponding to different stimulation phases are shown. (b) Medium inhomogeneous in the parameter gNa, its value is 1.0 in the lateral parts and 0.9 in the middle part, stimulation amplitude 0.055, stimulation period 256 time steps. Each time a trajectory approaches the middle part, it either draws into it and enters into a quasi-cycloidal motion, or is repelled back into the near-resonant region. In all cases, the vortex is always repelled from boundaries, in either part of the medium, and does not annihilate.

phenomenological theory will not be violated. As an illustrative example, note that the system (45), (50), (54) for an exponentially decaying interaction is invariant under the transformation x = x - Rlogp, c” = c/p,

(60)

i = t/p, at any ,M> 0. Hence, if at some z( the closest position of vortex to boundary is -/X0/, then at large p and corresponding small c the closest position will be - lzOl = - IX01 - R log p and this can be made as far from the boundary as desired by choosing an appropriate p.

Resonant drift of autowave vortices

Fig. 11. Evolution of multiple vortices under resonant stimulation. Stimulation amplitude is 0.05, stimulation period 254 time steps. Shown are excitation patterns in the same phase of stimulation (“stroboscopic view”) every 2 stimulation periods. Black dots show the “traces” of the tips of the vortices. The vortices either reflect from each other or from boundaries (“r”, vortices 3 on panel (b), 5 and 6 on (c), 3 and 4 on (e)) or annihilate (“x”: 1, 2 (a), 7, 8, (b)). It can be seen that 50% of vortices die out during first 10 stimulation pulses. Later, the rate of vortex annihilation drops (see Fig. 19).

4. SPATIALLY UNIFORM FEED-BACK STIMULATION

4.1.

The phenomenology

of the feedback resonant drift

Resonant repulsion is closely related to the change in the resonant frequency as the vortex drifts through different regions of the medium. It is this change in frequency that violates the resonance condition and produces a curved vortex trajectory. To avoid the resonant repulsion, let us measure the resonant frequency during the process of drift, in situ. One possible way to do this is to monitor activity at a point in the medium by a recording electrode. The resonant frequency is in fact the frequency of real rotation of the drifting vortex, in the reference frame moving with its core. Therefore, the frequency measured by electrode will differ from the resonant one only due to the Doppler effect and is a good approximation to the latter as the drift velocity is small compared to the velocity of wave propagation. This leads to the technique of the feed-back stimulation: the stimulation frequency is always equal to the monitored frequency when the stimulation is synchronized with the recorded activity. The position of the recording electrode, the criterion for detecting the arrival of a wave at the electrode, and the time delay between recording wave arrival and stimulation are all arbitrary. To estimate the importance of these factors we modify the asymptotic theory developed in Section 2 for the case of the feed-back stimulation. In the case of variable-frequency stimulation, the equation for + in (45) becomes Wt9

= Q@(4)

- Q(t)

(61)

602

V. N. BIKTASHEV

and A. V. HOLDEN

and we need one more equation for Q(t), describing formally the functioning of the feedback. Let the delay between monitoring the arrival of the wave and the stimulation be fixed, so that each stimulus is delivered at the same phase of vortex rotation, up to the phase distance between the current position of the vortex core and the recording. Thus the feedback determines the resonant drift direction 4 as a function of current vortex position and location of the registering point: G(t) = *&w)~ Y(t), x,3 Yr). (62) Here x,, yr are the coordinates of the recording electrode point. For example, if we approximate the form of the spiral wavefront by an Archimedean spiral, we get, instead of (45), the following second order ODE system: @,(X7 y, 4, Yr) = Jxx - 4, y - Y,) + 2&X

- 4,

3,X = C,(R) + c(R) COS@~(R- r,)

y - Y,)/A

+ 4

(63)

a,Y = C,,(R) + ~(R)sin@~(R - r,) where p, 6 are the polar coordinates, A is the asymptotic wavelength of the vortex and the symbol 6 denotes the sum of the stimulation delay and the constant determining the direction of the resonant drift at zero stimulation delay. The mechanism of the feedback driven resonant drift is illustrated by Fig. 12. The stimuli were delivered when the wave arrived at the top left corner of the medium. In contrast to the case of constant frequency, the trajectory of the vortex core far from boundaries is a curve, not a straight line. Close to the boundary, there is no resonant repulsion. The

Fig. 12. PDE stimulation illustrating the mechanism of the feed-back driven resonant drift. Stimulation amplitude A = 0.055, each stimulus is delivered at the moment the excitation wave arrived at the left top corner of the medium (shown by a cross). Due to the feed-back, each stimulus occurs at the same rotation phase, up to the phase distance from registration point to the vortex core. Therefore, the trajectory is affected only by “spatial” forces, i.e. usual attraction/repulsion of boundary, without resonant repulsion taking place. As a result, the vortex annihilates on the boundary. The excitation patterns are shown at the moments of stimuli beginnings, after each 2 stimuli.

Resonant drift of autowave vortices

603

trajectory deviates from what it would be in the absence of boundary, but this seems to be due to the C,,, terms - “spatial” interaction with boundary. As a result, the vortex reaches the boundary and annihilates, at a stimulation amplitude which was never sufficient for annihilation using constant frequency stimulation. 4.2.

Trajectories of the feedback driven resonant drift in an homogeneous medium

Typical trajectories of the phenomenological system (63) at C, = C, = 0 are shown in Fig. 13. It can be seen that the recording point is always a singular point in phase portraits, but its particular type is not of interest since the approximation of @ in (63) is not valid as (X, Y) is close to (x,, y,). It is only important to note that all the trajectories are not closed and leave any finite region independent of the specific value of the constant 6. The terms C,, C, in (63) can make the behaviour of trajectories more complicated. In particular, if their values are large compared to the resonant drift velocity c and their sign corresponds to repulsion from the boundary, then the trajectories will not reach the boundary. This is confirmed by the computations with the phenomenological ODE model and direct solution of the PDE model: in Fig. 14 we can see that for large values of c all trajectories reach the boundaries, for smaller values of c all do not. 4.3.

Feedback driven resonant drift in an inhomogeneous

media

Figure 15 shows trajectories of resonant drift under feedback, in inhomogeneous media (the parameter &, varies in the X direction) and in an homogeneous medium under

Fig. 13. (a).

604

V. N. BIKTASHEV

and A. V. HOLDEN

Fig. 13. (b) and (c).

Resonant drift of autowave vortices

Fig. 13. (d) and (e).

605

V. N. BIKTASHEV

and A. V. HOLDEN

Fig. 13. Feed-back resonant drift for the ODE system at C, = Cy = 0. (a)-(c) The phase portraits of the system at 6 equal to 0, 7r/3, 2rr/3 respectively. Note that at 6 = 71,4~/3, 5n/3 the portraits are the same with reversed directions of trajectories. The position of the recording electrode (in centre) is seen as the singular point at the portraits. (d)-(f) The bundles of trajectories emerging from the same initial point at different phase delays (0, n/3, 27r/3etc), shown for three different initial points and recording points. The recording point is shown by a cross in each panel. On (f) it coincides with the initial point.

inhomogeneous amplitude stimulation-this represents the case of a medium with inhomogeneous sensitivity to the stimulation. The phase portraits and bundles of trajectories show that the qualitative. character of the resonant drift is not changed significantly by the inhomogeneity-for comparison, see Fig. 14. In particular, the trajectories do not stop at the parameter values used within the medium, but always reach the boundary. The numerical experiments show that for annihilation of the vortex in an inhomogeneous medium, it is usually sufficient that close to the boundary, the same conditions hold as for a homogeneous medium with the same properties, i.e. the local stimulation amplitude should be greater than the threshold determined by local medium properties. It is interesting to compare qualitative predictions from the phenomenological theory with the numerical results presented on Fig. 15. Let us compare the “phase portraits” for the vortex drifting in these two types of inhomogeneities, and in corresponding homogeneous media. The “portraits” are obtained by performing simulations with various initial positions of the vortex. The simplest qualitative prediction following from the phenomenological theory is that trajectories in such “portraits” should behave as trajectories on the portrait of the ODE system of two variables (63). In particular, the trajectories should not intersect. From Fig. 15 we see that this is approximately valid, especially at low drift velocities. Next, if the sensitivity distance is sufficiently small, then the phase portraits of inhomogeneous media should be compositions of corresponding pieces of phase portraits of homogeneous media. From Fig. 15 we can see that this is indeed valid far from the

Resonant drift of autowave vortices

607

inhomogeneity borders; and visual comparison gives the sensitivity distance as about 5 space steps. This is of the same order of magnitude as the sensitivity distance chosen in Appendix for the ODE model to simulate the interaction of vortex with boundary. 4.4. Comparison of the effectiveness of constant frequency and feedback stimulation for elimination of vortices

In the case of feedback stimulation, there is a threshold value for the stimulation amplitude (which in the phenomenological system (63) is c) that is necessary for extinguishing the vortex by resonant drift to a boundary. There is also a threshold stimulation amplitude below which extinguishing the vortex is impossible using resonant drift produced by constant frequency stimulation. We now compare these thresholds; we use for these estimations the exponential form for interaction with boundary (50), (54). Here we assume that the medium is located at x 3 0; R is the characteristic length of the interaction with the boundary, and a, b, c are some constants, a > 0, b > 0. For the case of constant frequency stimulation, i.e. for trajectories of system (45), we find xCf = min{X(t)lX(-m)

= +m, 4(-m)

E (ST/~,31~/2), t E (-03, +m)}.

(64)

Using the parametric form for the solution of (45), (50), (54), obtained in [48], x = -R log (a-’ dq/dt),

Fig. 14. (a).

(65)

608

V. N. BIKTASHEV

and A. V. HOLDEN

0’)

Fig. 14. (b) and (c).

Resonant drift of autowave vortices

Fig. 14. (d) and (e).

609

V. N. BIKTASHEV

610

and A. V. HOLDEN

Fig. 14. Feed-back resonant drift for the ODE system at nonzero CX, Cr (b)=0.005, c=O.O03) and the trajectories of core in PDE simulation (a) PDE trajectories, A = 0.025, delays from 0 to 200 time steps equidistant with 40 time steps. The vortex does not annihilate at either delay. (b) ODE trajectories, u = 0.0025, 6 equidistant in 27r/7. (c) Phase portrait of ODE system at v = 0.0025 and S = 0. (d) Same as (a), A = 0.050. The vortex annihilates for each delay. (e) Same as (b), o = 0.005. The trajectory computation was stopped if it reached less than 1 space step from a boundary; this occurred for all trajectories. (f) Phase portrait of ODE system at u = 0.005 and 6 = 0.

drp/dt = -uR-‘sin(g,

+ e)cos(e) + Cexp(-vq)

(66)

Y = tan(e) = b/(aR); 8 E (-n/2, a-/2)

(67)

the minimum (64) can be written as -Gf =

-R log (u(aR)-‘(1

+ ~I+“*E(v))

(68)

where the function E() has a simple geometrical interpretation (see Fig. 16). Note that this function is both bounded 1 < E(V) < 2,

v E (0, +w),

E(0) = 2,

E(+w)

= 1,

(69) and monotonically decreasing, so xcf in (68) decreases monotonically with b. The condition that the minimal value of x should lie beyond the boundary, which is a necessary condition for annihilation, is therefore 0 > Ucf

= aR(1 + ti)“2E(v)-1.

(70)

We assumed that the critical position for the vortex is x = 0; choosing any other position does not influence the comparison. For the case with feedback, the necessary condition for the existence of trajectories that

Resonant drift of autowave vortices

611

cross the line x = 0 from right to left is simply v > vfb = b.

(71)

Comparison of (68) and (71) shows that Vfi/Vcf = v(l + y2)+!Z(V)

< 1;

Y = b/(d),

(72)

i.e. that feedback stimulation always needs smaller (for b/(d) << 1 very much smaller) amplitudes for extinguishing a vortex by resonant drift than constant frequency stimulation. The simple reason for validity of inequality (72) is that in case of feedback stimulation we need to overcome only the effect of spatial interaction with boundary, characterized by the coefficient b, while with a constant frequency we deal with both spatial (b) and temporal (a) effects, and it is natural that the relative effectiveness of these stimulations depends on the dimensionless ratio b/(aR). The “defibrillation thresholds” for extinguishing a vortex in the PDE model are illustrated in Fig. 17. The probability is obtained from averaging over a number of experiments (from 10 to 30) over a randomly chosen phase, where phase is relative to the phase of the vortex, and is obtained from the occurrence time for a single defibrillating pulse, the initial time for beginning of periodic series of stimulation, and the delay time for feedback stimulation. In the last case the delay should not be more than the period of rotation minus the duration of the stimulus. For the constant frequency curve the period of stimulation was chosen as close to resonance as possible (within the constraints of the discrete time lattice) from the values given in Fig. 3. In fact this curve (b) is impractical, as II,,,,T

I I I1 I1 I I I I lTTl, I I/I I I I1 I I! I I I I I r! I

612

Fig. 15. (b) and (c).

Resonant drift of autowave vortices

613

614

V. N. BIKTASHEV

and A. V. HOLDEN

Fig. 15. (f) and (g).

Resonant drift of autowave vortices

04

615

X

Fig. 15. Feed-back driven resonant drift in media inhomogeneous (a-d) in stimulation amplitudeand(e-h) in the parameter g&. (a-c) and (e-g) show “phase portraits” of the PDE model, i.e. vortex trajectories obtained in PDE simulations at its various initial positions, at zero stimulation delay, and (d), (h) are bundles of trajectories with the same initial position and properties of the middle part of the medium varying. The constant parameters are: g& = 1 time steps on (a-d) and A = 0.00044 (time steps)-* and stimulation delay 200 time stepson (e-h). The varying parameters are: (a) A = 0.0044 (time steps)-‘, (b) A = 0.0044in the lateral partsandA = 0.0025in the middle part, (c) A = 0.0025, (d) stimulation delay 100 time steps A = O.CWl in the lateral parts, and A = 0.0020, 0.0032, 0.0044, 0.0056 and 0.0064 for the trajectories l-5 respectively; (e) gNa= 1.0, (f) gNa= 1.0 in the lateral parts and 0.9 in the middle part, (g) gNa = 0.9, (h) stimulation delay 200 time steps,gNa= 1.0in the lateral parts and gNa = 0.84, 0.88, 0.92, 0.96, 1.0 for the trajectories 1-5 respectively. The phaseportraitsof inhomogeneous media appear to be constructed of the pieces of the corresponding phase portraits of homogeneous media, with the exception of vicinity of internal borders.

in practice it is impossible to obtain an exact resonant frequency, as the resonant frequency varies with different positions of the vortex in the medium. Feed-back stimulation can also extinguish multiple vortices that represent fibrillation. Since colliding wavefronts annihilate, the recording electrode monitors the activity of only one vortex at a time. A typical process is shown in Fig. 18. First the vortex that influences the recording electrode is driven by resonant drift under feedback, which causes the directed drift of at least one vortex. These “leading” vortices are circled on the figure. Other vortices may or may not annihilate during this stage. Upon reaching the boundary, the leading vortex is extinguished and the electrode monitors the wavefront from another vortex, which in turn is extinguished etc. As a result, all the vortices are progressively extinguished in a time not much longer than that needed for extinguishing one vortex: for the results of Fig. 18 it took 3 times longer to extinguish 9 vortices. Figure 19 compares the efficacy of the continuous frequency and feedback stimulation applied to multiple vortices. The graph shows how the number of vortices in the medium decreases with time, measured as the number of stimuli. We can see that, at these parameter values, the continuous frequency stimulation eliminates the vortices very slowly,

616

V. N. BIKTASHEV

and A. V. HOLDEN

Fig. 16. The geometrical definition for the function E(k) in (16). The graphs of the functions cos(x) and Ce-” are drawn, with constant C chosen so these curves are tangent at x = aarctan(k) (shown by arrow). The function E(k) for this value of k is the maximal distance between these curves to the right of the touching point (shown by arrows). It can be seen that 1 < E < 2 and E diminishes with k.

while feedback stimulation extinguishes all the vortices in two orders of magnitude less time. If these results are extrapolated to cardiac tissue, we would expect feedback stimulation to be effective in a few seconds, at a voltage of about 20% the single defibrillating shock threshold. 5. DISCUSSION

We presented here an asymptotic theory of the resonant drift phenomenon that accounts for the interaction of the vortex with boundaries, inhomogeneities and with each other. This theory leads to a simple description of the drifting vortex, in terms of its core location and rotation phase. The analysis of the phenomenological equations shows that boundaries and other influences change the behaviour of the drifting vortex e.g. “reflection” of the vortex from an impermeable boundary or from other vortices takes place. This phenomenon of reflection is of the vortex, not of its wavefront. The asymptotic theory is based on the symmetry and stability properties of the rotating vortex and therefore is applicable in general to effectively two dimensional autowave media. The phenomenology shows that the principal role in the repulsion is played by the change in the vortex frequency. The specific physical reasons are not important, and the mechanism of repulsion is related to the phase shift between the vortex and stimulation, caused by changes in vortex frequency. To illustrate the theory, we presented numerical simulations of resonant drift in a PDE

Resonant drift of autowave vortices

617

1 .o 0.9 0.8 0.7 0.6 0.5

0.4 0.3 0.2 0.1 0.0

-

v

5

1 6

I

7

8

9

10

F

20

3c

Fig. 17. Comparison of efficacy of different methods of stimulation. Abscissa: the stimulation amplitude * 100. Ordinate: the probability of annihilation of the vortex with occasionally chosen phase parameter. Initial position of vortex for all experiments is close to the center of the medium. Circles: single-pulse stimulation (classical defibrillation). Diamonds: constant frequency stimulation, with resonant frequency. Squares: feed-back stimulation.

Fig. 18. Evolution of multiple vortices under feed-back driven resonant stimulation. Stimulation amplitude is 0.05, stimulation delay is zero. The “leading” vortex is circled on each panel, traces of annihilated vortices are marked by digits in parentheses. (a) Vortex 1 leads, vortex 2 repulses from boundary. (b) Vortices 1, 2, 5, 6 have annihilated, vortex 3 leads. (c) Vortices 3, 8, 9 have annihilated, the only alive vortex 4 leads. Further evolution results in the annihilation of vortex 4.

model of autowave medium. Many qualitative predictions of the theory are observed in the PDE simulations. In particular, the phenomenon of reflection of resonantly drifting vortex is observed. In our model, the influence of the impermeable boundaries or inhomogeneities on the vortex frequency is well known from earlier computations (see e.g. [45, 47, 481). Other phenomena associated with heterogeneities include the movement of a vortex from one part of the medium to another one with different properties, or the vortex may get “attached” to the boundary. Note that such types of behaviour can occur even in the absence of a resonant external forcing [47].

V. N. BIKTASHEV

618

8

0

and A. V HOLDEN

v i

3t

1

10

100

1000

Fig. 19. Comparative efficacy of feed-back stimulation of multiple vortices. Shown is the number of vortices in the medium on the number of stimuli, obtained for different types of stimulation of the same amplitude and initial conditions (the same as on Fig. 11). Circles: constant frequency stimulation. Triangles: feed-back driven stimulation. Apparently, the last vortex which has survived after 1500 stimuli of constant frequency stimulation, will never annihilate; while feed-back extinguishes all the vortices by 30 stimuli.

Resonant drift also provides a means of estimation of parameters characterizing the medium and the vortex, which are difficult to measure directly. For instance, in the case of resonant repulsion, the reflection angle is very sensitive to the parameters of the vortex-boundary interaction, and we could obtain information about these interactions by measuring the dependence of the angle of reflection upon the angle of incidence, boundary properties and the stimulation amplitude. This method will yield new information about the interaction, which is difficult to obtain. In fact, we have already used this idea in this paper: we determined the dimensionless combination of spatial and temporal sensitivity of vortex and characteristic length b/(aR) by analyzing the experimental dependence @i(@,J, and this measurement is more precise than direct measurements. Application of the resonant drift effect to the control of vortices in cardiac tissue presents some serious difficulties. The frequencies of vortices may be not known a priori. Moreover, they may depend on the location of vortices and be different for different vortices, and in ventricular fibrillation we may be dealing with a system of interacting vortices. In principle, this difficulty may be overcome by scanning some reasonable frequency interval. The second difficulty is less obvious and is due to the repulsion effect described in this paper, which directly prevents the vortex from annihilating. To avoid this effect it is not sufficient just to measure accurately the vortex frequency and to select the frequency of stimulation, one should also increase the amplitude. We have shown for both the phenomenologic ODE model and the 2D PDE model that reflection of a resonantly drifting vortex by a boundary can be overcome by using feedback controlled stimulation, and that the threshold for extinguishing a resonantly drifting vortex using feedback is less than the thresholds using a single defibrillating shock. In the

Resonantdrift of autowavevortices

619

example considered the threshold is an order of magnitude less. The method is based on the symmetry of a rigidly rotating vortex, and so is general for all vortices with these properties. However, vortices in real cardiac tissue are asymmetric, mainly due to anisotropy in conduction pathways, are 3dimensional rather than 2-dimensional, and (especially in the ventricle) are multiple, as vortex breakdown kindles new vortices. Although the method has only been illustrated for symmetric vortices, we expect it to be robust for asymmetric vortices associated with linear, elliptical and Y-shaped cores or areas of functional conduction block [52, 531. Although such vortices may be obtained numerically in an anisotropic, homogeneous excitable medium [54] the effectiveness of the method on cardiac tissue requires electrophysiological rather than numerical investigations. Preliminary experiments with the PDE model in 3D have shown that resonant drift with feedback is effective in extinguishing vortices in an isotropic, homogeneous medium with absorbing boundaries. However, given the richness of vortex behaviour in 3D, with vortex filament twist, tension, and drift, the effects of boundaries and resonant drift in 3D media requires further exploration. In the case of multiple vortices, the approach also appears promising, but we used a homogeneous medium, so no new vortices were being formed. One possible scenario of fibrillation [6, 551 is that due to inhomogeneities of the medium, vortices break and give birth to other vortices. If this is the case, we would have to extinguish vortices faster than they were born. Numerical simulations and testing the feed-back stimulation technique, using quantitatively accurate models of cardiac muscle with inhomogeneities is required. Summarizing, the symmetry and stability properties of 2D autowave vortices enables their description as “quasi-particles”, using ordinary rather than partial differential equations, and analysis of the rather complex behaviour of these vortices under perturbations of various types. Though the “quasi-particle” description is formally valid only asymptotically, the main qualitative predictions of this description are valid and allow us to propose a technique for eliminating autowave vortices from a medium, and hence defibrillating cardiac tissue. Acknowledgements-This work was supported in part by NATO grant CRG 920656 (Leeds), and by a grant from the Wellcome Trust and from the Russian Fund for Fundamental Research No. 93-011-16080(Pushchino). REFERENCES 1. V. A. Davvdov, A. S. Mikhailov and P. K. Brazhnik, “Kinetika i gorenie”, p. 39, Institute of Chemical Physics, Chkrnogolovka (1986) (in Russian). 2. V. A. Davvdov. V. S. Zvkov. A. S. Mikhailov and P. K. Brazhnik. Drift and resonance of sniral waves in distributed Media, Zzv. Vizov 1 Radiofizika 31, 574-582 (1988) (Engl: trans. in Sov. Phys.-Radiophysics). 3. K. I. Agladze, V. A. Davydov and A. S. Mikhailov, An observation of resonance of spiral waves in distributed excitable medium, Pis’ma v ZETP 45, 601-603 (1987) (Engl. trans. in Sov. Phys.-JETP Letters). 4. M. J. A. Walker et al. The Lambeth Conventions: guidelines for the study of arrhythmias in ischaemia, infarction, and reperfusion, Cardiovascular Research, 22, 446-455 (1988). 5. G. R. Mines, On dynamic equilibrium in the heart, J. Physiol. (Land.) 46, 349-382 (1913). 6. W. E. Garey, The nature of fibrillary contraction of the heart - Its relation to tissue mass and form, Am. J. Physiol. 33, 397-414 (1914). 7. A. V. Holden, M. Markus and H. G. Othmer (eds), Nonlinear Wave Processes in Excitable Media. Plenum, New York (1991). 8. H. L. Swinney and V. I. Krinsky (eds). Waves in patterns in chemical and biological media, Physica D 49 (special issue) l-256 (1991). 9. J. Jalife (ed). Mathematical approaches to cardiac arrhythmias. Annals N. Y. Acad Sci. 591: 1-416 (1990). 10. A. M. Zhabotinsky, Concentrational Auto-Oscillations. Nauka, Moscow (1974) (in Russian). .l. N. Wiener and A. Rosenblueth, The mathematical formulation of the problem of conduction of impulses in a network of connected excitable elements, specifically in cardiac muscle, Arch. Inst. Cardiol. Mex. 16, 517-525 (1946).

V. N. BIKTASHEV

620

and A. V. HOLDEN

12. M. A. Alessie, F. I. M. Bonke and F. J. G. Schopman, Circus movement in rabbit atria1 muscle as a mechanism of tachycardia, Circ. Res. 33, 54-62 (1973). 13. M. J. Janse, Re-entrant arrhythmias, in Heart and Cardiovascular System, Second Edition, edited by H. A. Fozzard et al. DO.2055-2094. Raven Press, New York (1991). 14. A. M. Zhabotinsky’ and A. N. Zaikin, “Spatial phenomena in an autooscillatory chemical system, in Oscillatory Processes in Biological and Chemical Systems, Vol. 2, pp. 279-283. The Center of Biological Research, Pushchino (1971). 15. A. T. Winfree, Spiral waves of chemical activity, Science 175, 634-636 (1972). 16. V. S. Zykov, Simulation of Wave Processes in Excitable Media. Nauka, Moscow (1984) (in Russian); Manchester University Press (1988) (in English). 17. J. J. Tyson and J. P. Keener, Singular perturbation theory of travelling waves in excitable media, Physica D 32, 327-361 (1988). 18. A. T. Winfree, Stable particle-like solutions to the nonlinear wave equations of three-dimensional active media, SIAM Review 32, l-53 (1990). 19. G. W. Beeler and H. J. Reuter, Reconstruction of the action potential of ventricular myocardium fibers, J. Physiol. 268, 177-210 (1977).

20. M. Courtemanche and A. T. Winfree, Re-entrant rotating waves in a Beeler-Reuter based model of two dimensional cardiac electrical activity, Int. .I Bifurcation & Chaos 1, 431-444 (1991). 21. I. R. Efimov and V. I. Krinsky, A Study of Cycloidal Behaviour of Vortices in the Beeler-Reuter Model of Myocardium. The Center of Biological Research, Pushchino (1991). 22. I. R. Efimov, V. I. Krinsky and J. Jalife, Dynamics of rotating vortices in the Beeler-Reuter model of cardiac tissue, J. Nonlinear Sci., to appear. 23. I. L. Gurvich. Basic Princioles of Heart Defibrillation. Meditsina. Moscow (1975) (in Russian) 24. A. M. Pert&v, V. N. Biktashev, E. A. Ermakova and V. ‘I. Krinsky, On‘the possibility of lowering considerably defibrillation current by determining the right time for application of a defibrillating pulse: mathematical model Biofizika 35 500-503 (1990) (in Russian). 25. V. I. Krinsky, A. M. Pertsov and V. N. Biktashev, Autowave approaches to cessation of re-entrant arrhythmias, Ann. N. Y. Acad. Sci. 591, 232-246 (1990). 26. V. I. Krinsky, K. I. Agladze, Interaction of rotating waves in an active chemical medium, Physica D 8, 50-56 (1983). 27. E. A. Ermakova, V. I. Krinsky, A. V. Panfilov and A. M. Pertsov, Interaction between spiral and flat periodic waves in an active medium, Biofizika 31, 318-323 (1986) (in Russian). 28. P. K. Brazhnik, V. A. Davydov and A. S. Mikhailov, Kinematical approach to description of autowave processes in active media, Teoreticheskaya i Matematicheskaya Fizika 74, 440-447 (1988) (in Russian). 29. A. S. Mikhailov, Kinematics of wave patterns in excitable media, in Nonlinear Wave Processes in Excitable Media, edited by A. V. Holden, pp. 127-144. M. Markus, H. G. Othmer. Plenum, New York (1991). 30. A. M. Pertsov and A. V. Panfilov, Spiral waves in active media. Reverberator in the Fitz-Hugh-Nagumo model, in Auto-wave Processes in Systems with Diffusion, pp. 77-84. IPF, Gorky (1981) (in Russian). 31. A. V. Panfilov and A. M. Pertsov, Mechanism of spiral waves initiation in active media related to the critical curvature phenomenon, Biofizika 27,886-889 (1982). 32. A. V. Panfilov, Reverberators with anomalous core in a relaxational membrane system, Ph.D. thesis, Moscow Institute of Physics and Technology (1983). 33. P. K. Brazhnik, V. A. Davydov, V. S. Zykov and A. S. Mikhailov, Vortex rings in excitable media, ZETP 93, 1725-1736 (1987). 34. P. K. Brazhnik, V. A. Davydov, V. S. Zykov and A. S. Mikhailov, Dynamics of autowave structures, in II All-Union Conference on Mathematical and Computational Methods in Biology, pp. 118-119. Pushchino (1987). 35. V. A: Davydov and V. S. Zykov, Spiral waves in anisotropic excitable media, ZETP 95, 139-148 (1989). 36. P. K. Brazhnik, V. A. Davvdov, V. S. Zvkov and A. S. Mikhailov, Twisted vortex in excitable medium. Izv. VCJZov - Radiofizika

32; (1989).

.

37. A. Yu. Abramychev, V. A. Davydov and A. S. Mikhailov, On the theory of resonance in excitable media, Biofizika 35, 504 (1990). 38. A. Yu. Abramychev, V. A. Davydov and V. S. Zykov, Drift of spiral waves on non-uniformly curved surfaces, ZETP 97, 118881197 (1990). 39. V. A. Davydov and V. S. Zykov, Spiral waves in a circular exitable medium, ZETP 103, 844-856 (1993). 40. A. Yu. Loskutov and A. S. Mikhailov, Introduction into Synergetics. Nauka, Moscow (1990) (in Russian). 41. A. Garfinkel, M. L. Spano, W. L. Ditto and J. N. Weiss, Controlling cardiac chaos, Science 257, 1230-1235 (1992). 42. V. I. Arnold, Additional Chapters of the Theory of Differential Equations. Nauka, Moscow (1978) (in Russian). 43. V. N. Biktashev, On the geometrical optics for nonlinear waves, in Nonlinear Dispersive Waves, edited by L. Debnath. World Scientific, Singapore, (1992). 44. V. N. Biktashev and A. V. Holden, Tension of organising filaments of scroll waves, Phil. Trans. Roy Sot (London) (in press). 45. A. M. Pertsov, E. A. Ermakova, Mechanism of the drift of a spiral wave in an inhomogeneous medium, Biofizika

33, 338-342 (1988).

46. B. N. Vasiev and A. V. Panfilov, A Study of the Process Iof Initiation

of the Reverberator

in the Medium

Resonant drift of autowave vortices

621

Inhomogeneous in the Refractoryness. The Center of Biological Research, Pushchino (1989) (in Russian). 47. E. A. Ermakova and A. M. Pertsov, Interaction of rotating spiral waves with a boundary, Biofizika 31, 855-861 (1986) (in Russian). 48. V. N. Biktashev and A. V. Holden, Resonant drift of the autowave vortex in a bounded medium, Physics Lett. A 181, 216-224 (1993). 49. V. N. Biktashev and A. V. Holden, Design principles of a low-voltage cardiac defibrillator based on the effect of feed-back resonant drift. J. Theoret. Biol. (in press). 50. L. D. Landau and E. M. Lifshitz, Course of Theoretical Physics 2: Field Theory. Pergamon Press, Oxford (1960). 51. V. N. Biktashev, Drift of a reverberator in an active medium due to interaction with boundaries, in Nonlinear Waves ZZ. Dynamics and Evolution, Edited by A. V. Gaponov-Grekhov, M. I. Rabinovich and E. J. Engelbrecht, pp. 87-96. Springer, Berlin (1989). 52. S. Dillon, P. S. Ursell and A. L. Wit, Pseudo-block caused by anisotropic condition: a new mechanism for sustained re-entry, Circulation 72, 1116 (1985). 53. A. V. Panfilov and A. V. Holden, Computer simulation of re-entry sources in myocardium in two- and three-dimensions, J. Theor. Biol. 161; 271-285 (1993). 54. V. I. Krinsky, I. R. Efimov and J. Jalife, Vortices with linear cores in excitable media, Proc Roy Sot Lond. A437: 645-655 (1992). 55. V. I. Krinsky, Fibrillation in excitable media, Problemy Kibernetiki 20, 59-68 (1968) (in Russian). 56. V. Levitin, TraX - Simulation and Analysis of Dynamical Systems” available from: Exeter Software, 100 North Country Road, Setaket, New York 11733.

APPENDIX:

THE NUMERICAL

MODELS

The PDE model

Numerical experiments were performed using a piece-wise linear version of the Fitz-Hugh-Nagumo similar to that proposed in [30], in R*:

model,

ut = gNa(f(u) - u) + V*u + F(r) 1 Ut = E(U)(gU - u) where f(u) - is an N-shaped piecewise-linear continuous function,

f(u) = 1

-g,u,

u < u,

-g1u1 + &i" - Ul)

u, < u < u*

g*(l

u >

-

u)

U2r

E(U) is a piecewise-constant function, E(U)

=

I

&IT

u < lil

.52,

u,

&3>

u >

<

u <

u*

U2,

and the default parameter values are: g = 1.0,

g1 = 4.0,

u, = 0.08,

gf = 0.98,

El = 1.68,

g* := 15,

E* = 0.06,

&3 =

g‘wl = 1.0, 1.50.

Calculations were performed by the explicit Euler method with a 5-node approximation of the Laplacian on a rectangular grid of 60 x 60 internal nodes with the time step of 0.08 and space step of 0.8. The left and right boundaries were assumed impermeable -

alA ax

= 0. x=x nun:x "mr

The upper and lower boundaries were either assumed impermeable, or we used the periodicity conditions: u, ~yuly=y,,” = u, WY=&, to obtain a cylindrical medium, convenient for analyzing interaction with a single straight boundary. The external forcing current was simulated as F(t) added to the right-hand side of the equation for the variable u, F(r) was a series of rectangular pulses of duration of 50 time steps, with amplitude and period varying in different experiments. The default parameters were selected so that the free vortex rotates rigidly with a constant frequency, as illustrated in Fig. 1. The position of the vortex was determined by two methods: (1) “Averaged core”. At first, the maximum of u over a period was calculated at every grid point: u&f, Y, t) = max {u(X, Y, t), t E (t - At/2, t + At/2)}

V. N. BIKTASHEV

622

and A. V. HOLDEN

for the time moments t = nAt, n = 1, 2, . . These are equidistant with a period At, chosen to be 280-300 time steps, that is a little longer than the period of the vortex (which was within the range 260-280 time steps in different situations). Next, the amount u, below some threshold value uI, chosen to be 0.58 in our experiments, was calculated: if u,(x, Y, t) < 4, Ut - 4&, Y, tL 0 otherwise. i > At last, the “center of mass” of the function nd at the moment t was calculated: ud(x,

y,

t)

=

xc(t) = ud(x, Y, t)XdXdY/M(t), I Y,(t) = &,(x, Y, t)YdXdY/M(t), I where M(t) = ud(x, Y, t)dXdY. I The coordinates XC, Y, provide efficient approximants for the rotation centre of the vortex, and since they need not be at a grid node enhance the spatial resolution. (2) “Instantaneous tip”. Sometimes, the procedure of “averaged core” was ineffective -in inhomogeneous media, where it was difficult to select an averaging period, or in the case of multiple vortices. In these cases, we divided the phase plane of the diffusionless system into three parts (first u < 0.5 and v < 0.5, second u > 0.5 and u < 0.5 and third v > 0.5), which divides the spatial excitation pattern into three regions. Then we identified points in contact with nodes of all three types. For each vortex, this procedure detects at least one point, which provides an approximation for the instantaneous tip of the wavefront. This procedure has a lower spatial resolution, but can be applied to any spatial pattern of excitation in R*. The ODE model

Numerical solutions of the PDE model were accompanied by simulations of the phenomenological system, in the form (45) for constant frequency resonant drift and (63) for feedback resonant drift. The functions in (45) were chosen c(X, Y) = c = const.; C,(X, CJX, 0(X,

Y) = a(G(x)

Y) = b(G(X)

- G(L - X));

Y) = b(G(Y)

- G(L - Y));

+ G(L

- X) + G(Y)

G(x)

+ G(L

- Y)) + 6~;

= e-x/R.

and the default parameter values are c = 0.005;

R = 4.0;

a = 0.005;

b = 0.005;

L = 60;

A = 30;

and 60 and 6 varied-these are the deviation of the forcing frequency from the resonant frequency and the stimulation delay respectively. For comparison with PDE simulation in cylindrical medium, we modified the ODE systems by removing the terms G(Y) and G(L - Y) from the equations for Y and @. The phenomenological systems were investigated with the program TraX [55], using Runge-Kutta and Euler methods.