Kinematical theory of spiral waves in excitable media: Comparison with numerical simulations

Physica D 52 (1991) 379-397 North-Holland

Kinematical theory of spiral waves in excitable media: Comparison with numerical simulations A.S. M i k h a i l o v Institute for Theoretical Physics and Synergetics, Unicersity of Stuttgart, Pfaffenwaldring 57/IV, W-7000 Stuttgart 80, Germany and Department of Physics, Moscow State Unil,ersity, 117234 Moscow, USSR

and V.S. Z y k o v Institute of Control Sciences, 117806 Moscow, USSR Received 8 October 1990 Revised manuscript received 6 April 1991 Accepted 9 April 1991 Communicated by A.V. Holden

The kinematical theory of spiral waves is compared with the data of numerical simulations using complete partial differential equations of a reaction-diffusion model. We find good agreement for temporal periods of spiral waves which rotate around minimal holes and in the free medium. No adjustable parameters were used in the calculations.

1. Introduction

Rotating spiral waves in excitable (or oscillatory) media are currently the subject of intensive experimental [1-6], numerical [2, 7-11] and analytical [1, 12-17] investigations. Starting from the pioneering work by Wiener and Rosenblueth [18], it has been known that spirals are formed when an excitation wave rotates around a hole in the excitable medium. However, numerous experiments and simulations (see, e.g., the above mentioned references [1-11]) provide evidence that the same wave patterns are observed, as well, in media free of any obstacles. In this case, when the central region of a spiral wave is considered with sufficiently good resolution, one can often see that the center is occupied by a "core". The

elements inside the core stay unexcited, and the wave rotates around it, as if it were an effective hole. Recently, Keener and Tyson [12] performed a detailed analytical study of spiral waves in the Oregonator model, which is thought to describe the kinetics of the Belousov-Zhabotinsky reaction. Although their work is remarkable in many respects and makes a significant step forward in the analysis of spiral waves in chemical excitable media, it still does not provide a consistent description. It assumes that the core radius is known beforehand and does not give any recipes for calculation of this quantity. Since the core radius determines the rotation frequency of the wave, this also means that no definite estimate for this frequency is available. This inconsistency was

0167-2789/91/$03.50 © 1 9 9 1 - Elsevier Science Publishers B.V. (North-Holland)

380

A.S. Mikhailor, KS. Zykot' / Spiral waves in excitabh' media

pointed out recently by Jahnke and Winfree [11] who carried out an extensive numerical simulation of spiral waves in the Oregonator model. The principal difference of the approach [19-23] which has been developed by the present authors and is known under the name of kinematical theory, lies in the fact that no hole is introduced ad hoc into the theory for spiral waves. The kinematical theory incorporates the dynamical equations which allow to describe the entire evolution of the wave from an arbitrary initial configuration into a steadily rotating wave, when such a steady solution is stable [20], or to describe the complicated meandering regimes [20, 24-26] when they are realized. The theory was also applied to study other nonstationary effects, such as drift in the media with spatial gradients [20, 23, 26-29] or resonance in media with periodically modulated excitability [20, 27, 30]. In the kinematical theory, any excitable medium is characterized by a few parameters. Since the same parameters are involved in determination of different properties, such as the rotation frequency of the wave, the time of relaxation to steady rotation, the speed and the direction of drift in a given gradient or the resonance effects in periodically modulated systems, they can be determined independently from one set of experiments and then used to predict the behavior and the properties of spiral waves under other conditions. However, it is also possible to directly calculate the fundamental parameters of the kinematical theory, if the equations of an excitable medium are known. In this paper we present, for the first time, a systematic procedure that allows to carry out such calculations. We also give an example of this calculation for a particular reaction-diffusion system. Once the kinematical parameters are found, they can be used to predict, without any further adjustments, the basic properties and the behavior of spiral waves. Hence, by comparing the predictions with the actual p h e n o m e n a observed in numerical simulations of reaction-diffusion

systems or in real experiments, one can test the validity of this theoretical approach. In the present paper we undertake two different tests. First, we compare the theoretical predictions for rotation periods of free spiral waves with the data of a direct numerical simulation within a certain range of parameters. Second, we carry out the same comparison for spiral waves which rotate around "minimal" holes. A hole is "minimal" if, when its radius is further diminished, the tip of a spiral wave becomes detached from the hole boundary and, after a transient, a steady rotation is established with the core radius which is larger than the hole which is located in its center. This phenomenon, which was first observed in numerical simulations [7], is described quantitatively by the kinematical theory. In section 2 we provide a brief outline of the kinematical theory of spiral waves. The procedures for estimation of the two critical curvatures, for a solid wavefront and for a broken wave with a free tip, are formulated in section 3. The dependences of the rotation frequency of spiral waves on the medium excitability, that are obtained within the kinematical theory and by numerical simulations, are compared in section 4. In the last section we discuss the obtained results. Finally, the appendix deals with a possible generalization of these results to models with nonvanishing diffusion constants of both the activator and the inhibitor.

2. K i n e m a t i c s of spiral waves

In this paper (except the appendix) we consider only two-component excitable systems with diffusion for the activatory component and a slowly varying inhibitory component. Such excitable media are described by the models

=f(u,~,) + O,,Au, ,', = e g ( u , ~ ) ,

(1)

(2)

381

A.S. Mikhailov, V.S. ZykoL /Spiral waves in excitable media -¢ =

= VOO

¥

C~

O°5

o

Fig. 1. The null-clines for the model of an excitable medium given by eqs. (1), (2) and (47)-(50). where e is a small parameter. T h e null-clines = 0 and g(u, v) = 0 have a single intersection point (fig. 1) that corresponds to the stable uniform rest state of the medium. U n d e r appropriate conditions (see refs. [21, 27]) such an excitable m e d i u m may support the u n d a m p e d p r o p a g a t i o n of waves or pulses (in one dimension). T h e p r o p a g a t i o n velocity V0 of a solitary pulse is a unique property of a given medium. W h e n e is small, it d e p e n d s linearly on this parameter, i.e. as (see ref. [21])

f(u, v)

v,,

= v,,,,(l

(3)

where X is some m o d e l - d e p e n d e n t proportionality factor. The d e p e n d e n c e of the propagation velocity V0 for larger values of e is shown in fig. 2. A m e d i u m is considered to be more excitable #~ if, u n d e r fixed other properties including the diffusion constant, the propagation velocity of travelling pulses in such a m e d i u m is higher. F r o m #1This definition of excitability was used in the previous papers of these authors [21-30] and in refs. [2, 5, 11]. It will

be used throughout the paper. Note that sometimes "excitability" is defined in a different way, as a normalized threshold near equilibrium (see refs. [8, 17]). This should be remembered to prevent confusion.

o.a

o.~

~L

Fig. 2. The dependence of V0 on e for the model (1), (2) and (47) (50).

this point of view, the effective excitability of the m e d i u m can be controlled by variation of the p a r a m e t e r e in eq. (2). Indeed, it follows from (3) that, when e is increased, the pulses begin to move slower. This implies that the effective excitability is lower for larger values of e. As is seen from fig. 2, there is a critical value ecl at which the solution in the form of a travelling solitary pulse disappears. Note that at this critical excitability the pulse still propagates with a nonvanishing velocity. T h e typical profile of a travelling solitary pulse is shown in fig. 3. O n e can distinguish here the top (where the activator density u is large) and the recovery tail (where the inhibitor density v slowly goes back to the ifiitial low value). The durations of the top (T,) and of the recovery tail (T,) are two important properties of the medium; their ratio can be controlled by varying the form of the functions f(u,v) and g(u,v) (see section 4). The pulses in a periodic train begin to feel the presence of others when the time interval between them a p p r o a c h e s the total width of a single pulse, i.e. the interval T,, + T~. T h e n the head of a pulse begins to run into the recovery tail of the preceding pulse and, hence, the propagation velocity decreases (fig. 4). A minimal temporal pe-

A.S. Mikhailot', V.S. Zykot, / Spiral waees in excitable media

382

(a)

Io'-.. Z (b) I I i

0

Fig. 3. The profiles of the activator (a) and the inhibitor (b) distributions in a solitary travelling pulse for the model (1), (2) and (47)-(50).

q.o

0

Tmi n

20

40

T

Fig. 4. The dependence of V on the temporal period T for the model (1), (2) and (47)-(50); e = 0.25.

riod Tmin of stable pulse propagation exists, such that any train of waves with a smaller period is unstable. In a two-dimensional medium, flat or curved waves can propagate. Obviously, the properties of flat waves are the same as those of the travelling pulses in the one-dimensional problem. T h e normal propagation velocity V of the curved waves d e p e n d s on their curvature K. T h e wave which is convex in the direction of propagation has a smaller propagation velocity. The typical dependence V ( K ) for the models (1), (2) is shown in fig. 5 (see further discussion in section 3). We can

0.2

i

0.~

]

K

C[

r

0.6

0.8

K

Fig. 5. The dependence of V on the curvature K for lhe model (1), (2) and (47)-(50); e = 0.25.

note that there is a critical curvature K d, such that no waves with higher curvatures can propagate in a given medium. The propagation velocity V(Kcl) of the wave with the critical curvature K,,I is not zero. If we start with the initial condition in the form of a broken flat wave, then its further evolution would d e p e n d cardinally on the excitability of the medium. It turns out that there is an additional critical value Ec2. For ecl > e > ec2 the flat broken wave retracts at its tip. In the media with e < e,.2 the wave tip sprouts and curls back, thus producing a spiral wave (fig. 6). Hence, the second critical value ec2 is defined by the condition that the broken flat wave propagates in thc normal direction, neither sprouting nor retracting at its tip. Spiral waves exist only in media with e < ec2. T h e numerical simulations of spiral waves in the models (1), (2) reveal [21, 33-35] that their properties essentially d e p e n d on the m e d i u m ' s excitability (i.e. on the value of e) and on the ratio T , / T , of the excitation time to the recovery time. W h e n e approaches e~2, the rotation period T 0 of the spiral wave diverges, together with the core radius (i.e. the radius of a circle a r o u n d which the tip moves). In this region T 0 is much larger than both T,, and T,.

A.S. Mikhailoc, KS. Zykoc / Spiral waces in excitable media (a)

(b)

(c)

Fig. 6. T h e time evolution of a flat b r o k e n w a v e in the m e d i a d e s c r i b e d by the eqs. (1), (2) and ( 4 7 ) - ( 5 0 ) w i t h four different v a l u e s of the p a r a m e t e r e: (a) a - 0.4, (b) e = ec2 = 0.388, (c) e = 0.35, (d) e = 0.3. T h e isolines of c o n s t a n t a c t i v a t o r c o n c e n tration u = 0.4 are s h o w n at c o n s e q u e n t time m o m e n t s s e p a r a t e d by the intervals At = 6.

For smaller values of e (i.e. for higher excitabilities) T 0 decreases and becomes comparable, to T~, and T,. Then, as was found in refs. [26, 31] by numerical simulations of model (47)-(50), the properties of the spiral wave may be very sensitive to the ratio T~,/T,. If the excitation duration T~, were much smaller than the recovery duration T,., numerical simulations [26, 31] would show that steady circulation of spiral waves in the media with high enough excitability would be unstable. Under these conditions, meandering would be observed in model (47)-(50) in refs. [26, 31]. On the other hand, if T, were of the same order of magnitude or larger than T,, the steady rotation of spiral waves in model (47)-(50) would be stable, even in media with high excitability. In this case the spiral wave would have special properties. Its core would be very small (with the radius about the width of the tip) and it would be filled by the inhibitor (while in the other regimes the core would be void of the inhibitor). Perhaps

383

it would be the increased concentration of the inhibitor in the core region which stabilized the rotation and suppressed meandering of spiral waves in this model. Meandering regimes were reported as well for other mathematical models and for experiments with the Belousov-Zhabotinsky reaction (see refs. [2, 5, 8-11]). A detailed discussion of meandering regimes in these other models or experimental media lies out of the scope of the present paper. We can only conjecture that in these cases the ratio of the excitation duration to the recovery duration may also play a significant role in transition from meandering regimes (at small values of this ratio) to steady rotation (when this ratio is about one or larger). In the latter situation we would expect an increased concentration of the inhibitor in the core region, as was found in simulations [26, 31]. The kinematical theory of spiral waves is mainly applicable when the excitation duration T,, of rotating waves is much smaller than the rotation period T 0 of a spiral wave. When this assumption holds, the excitation zone of a spiral wave (which represents the region on the plane where the elements of the medium are in the states corresponding to the top of the pulse) is narrow in comparison with the core radius and the typical distance between the consecutive coils of the spiral. In this case, if we consider the spiral wave pattern with a somewhat poorer spatial resolution, its momentary excitation region would look like a curved broken line of vanishingly small thickness. Its end point would correspond to the tip of the excitation region. Note that, with such a resolution, we would not be able to distinguish the fine structure of the tip and to see that actually it is rounded. It would be simply the end point of the curve. It is therefore meaningless to ask, within such a spatial resolution, what is the curvature radius of the tip: it is zero. In the kinematical theory the excitation zone of the spiral wave is geometrically modelled by a broken curve. This curve is oriented, which means that for each its point the normal propagation

384

A.S. Mikhailot ~, V.S. Zykot / Spiral wa;'es in ~citable media

direction is indicated. Note that the only intrinsic property of such a curve is its local curvature K. It is this property which determines the local normal propagation velocity V of the segments of such a curve. U n d e r typical conditions of excitable media, V is a decreasing function of the curvature K. This d e p e n d e n c e is explained by the fact that, for larger K ' s and hence for smaller curvature radii, the flow of the "activator" goes into a wider angle sector and thus the efl'ective "excitability" is decreased. The dynamics of the curve is not restricted to the normal motion of its segments. The curve can also grow or retract at its end. This corresponds to sprouting or retraction of the tip. It means that, to have a complete dynamical description, one should indicate not only the normal propagation velocity V but formulate, as well, the law of motion of the end point. Consider an initial condition in the form of a fiat broken wave. Provided that the excitability of the m e d i u m is sufficiently high, this wave will start to move in the normal direction with some velocity V0 and, at the same time, sprout in the tangent direction with some speed G~. Thus we see that there are two parameters (V 0 and G 0) associated with the wave's motion. Note as well that, when the excitability of the m e d i u m is decreased, the tangent velocity G 0 decreases and, at some critical excitability, it can change its sign (so that sprouting is replaced by retraction). Suppose now that we start with a broken curt,ed wave, which represents a segment of a circle with some constant curvature K. Then, as it was already mentioned, this would influence the normal propagation velocity V. But this would also modi~, the speed o f sprouting at the end point. Thus we expect that the speed of sprouting would be a certain function of the curvature K, i.e. G = G ( K ) . To prevent possible confusion, we emphasize that K here is not the "curvature of the tip" (it is infinite in the present description), but the curvature of the line itself. W h e n the form of the line is not circular, its curvature K is not constant. T h e standard proce-

dure in this case would be to divide the curve into small segments with local curvatures K. Any such segment would have its own normal propagation velocity V(K). The speed of sprouting of the end point G would be d e t e r m i n e d by the properties of the first segment which includes the end point. Namely, G would be some function of the curvature K 0 of this first segment. W h e n the size of the segments diminishes, K 0 approaches the line curvature at the end point (not in the end point itself). W h e n the local curvature increases, the normal propagation velocity V becomes smaller. The same happens (for the same reasons u n d e r the same conditions) with the sprouting speed G as a function of K u. W h e n K 0 increases, the speed G changes its sign at some critical curvature K~e and for higher K0's the end point is retracting. Near K o = K~2 we can approximately write G = y(Kc2-K~,),

(4)

where y is some positive coefficient. We discuss explicit recipes for calculation of K,. 2 and y in section 3. Note again that K o is defined as a limit:

K0= lim K(I),

(5)

where 1 is the arc length of the curve m e a s u r e d from the end point. Consider a steadily rotating spiral wave. As we know from many numerical simulations and experimental data, its tip moves around a circular core, neither diving into the core region nor departing from it. In the light of the above discussion, it suggests that, in the steady rotation regime, K 0 coincides with Kc2. In other words, u n d e r the conditions of stable steady rotation, the curve that represents the excitation zone of the wave is coming to the core with the curvature which is simply equal to the critical curvature at which sprouting is replaced by retraction. If the local properties of the medium (its excitability, etc.) are not uniform, both the normal

A.S. Mikhailot', V.S. Zykoc / Spiral wal'es in excitable media

propagation velocity V and the sprouting speed G become, in addition, explicit functions of coordinates x and y. This should be taken into account when, for instance, there is a gradient of excitability along a certain direction. Such a situation is realized in the experiments where the drift of spiral waves is observed. Even if the original properties of the excitable medium are uniform, the effective nonuniformity can be created by the propagating waves themselves. After the excitation zone has passed through a given point, the residual concentration of the inhibitor is left there. During the recovery period after the passage of the wave, this residual concentration slowly decreases and eventually the initial state of the medium is recovered. However, if the next wave passes through the same point within the time shorter than the recovery duration Tr, it will effectively experience somewhat lower excitability. It means that both V and G will be some functions of the time interval T measured from the last passage of the excitation zone through a given point of the medium. When two dependences V = V ( K , T ) and G = G ( K , , T ) are specified, the motion of the curved line with an end point is completely determined. The motion of the curve can be most conveniently described in terms of the so called "natural equation" which expresses K as a function of the arc length l measured from the end point, i.e. K = K ( l , t). If the "natural equation" is known, it determines the curve up to its translation and rotation in the plane. The evolution equation for K ( l , t) was derived, using the methods of differential geometry, in refs. [14, 22, 32-35]. It has the form OK O--i- +

IVK dl' +

G

-~-

= - K 2V - -Ol 2 .

(6)

Here V = V ( K ) is a known function of K, so that (6) represents a closed equation for K ( l , t). Since the form of a steadily rotating spiral does not change with time, it corresponds to a stationary solution of eq. (6). The boundary conditions

385

for this equation are chosen as K(0) =Kc2

(7)

(so that the tip neither grows nor contracts, and G = 0) and

lim K(/) =0.

(8)

l ---) ~¢

In the stationary case eq. (6) can be once integrated yielding

dV

[ z V K dl'

dl + K j0

= w.

(9)

The integration constant w in (9) is simply the angular rotation velocity of the spiral wave. Indeed, it follows from (9) that oJ = ( d V / d l ) l _ o. But since G = 0 t h e wavefront is orthogonal to the core at its end point and therefore at l = 0 we have d l = d r (r is the polar radius). Because V is a definite function of K, (9) is an integro-differential equation for K. Note that it includes only the first derivative while we have two boundary conditions (7) and (8). Hence, (9) with (7) and (8) set a nonlinear eigenvalue problem which fixes the value of the rotation frequency w. To solve eq. (9), one should employ some dependence of V on K. If the curvatures K ( l ) for the entire solution are sufficiently small, the linear approximation #2 V = Vo - D K

(10)

can be used. It is convenient to introduce the dimensionless variables K = K / K c 2 and y = Kc2l. Using the dependence (10), we can write eq. (9) in terms of ~2Generally, this approximation holds only for small curvatures K, when DK << V~. However, in the special case of the model (47)-(50), which is used in section 4 for numerical simulations, it remains valid in a wide interval of curvatures K (see fig. 5).

A.S. Mikhaikn', ldS. Zvkor / Spiral wat'es in excitable media

386

DKc2(To)/Vo(To). Then (13) and (14) become the

these new variables as

K

dR-

11)

1 -/3K) d y ' =.Q + / 3 a 7 , '

equations which should be solved to determine the rotation period T0 of steadily rotating spiral waves.

where .Q w/VoK¢. 2 and [3 = D K c 2 / G I. The boundary conditions then take the form

When the rotation period T 0 is known, the radius R 0 of the core circle (around which the tip goes) can be found as

K(0) = 1,

Ro = K,(To) T , , / 2 w .

lira K ( y ) = 0 .

(12)

Note that, since the linear approximation (10) for the propagation velocity is used, this implies that /3 should be a small parameter. The solution of (11) with the boundary conditions (12) yields ,(2 as a certain function of /3. This function was determined in ref. [19] by numerical integrations. In the interval from 1 >/3 > 0 it is approximated with the accuracy of up to one percent by the analytical expression ~Q

0.685/3 ~/2 - - (I.06/3 - 0.2 9 3/3"9 +

....

(13)

The rotation period T~ = 2"rr/w of spiral waves is then given by

2"rrD

(14)

T,, - V,~/3~(/3) " Note that expression (14) is valid only while the rotation period T0 is much larger than the recovely duration T, (this holds for the media with very small critical curvatures Kc2, s e e section 3). Otherwise we should take into account the dependences of V0 and Kc2 on the time interval T between passage of the consecutive waves. Such dependences Vo(T) and Kc2(T) must be determined beforehand for the particular considered model. When a spiral wave is rigid, the interval of time between passage of the consecutive waves through any point of the medium (except the core) is the same and equal to the wave period T0. Hence, we should put Vo(T o) and K~2(T o) in expressions (13) and (14), and to assume that /3 =

(15)

Eq. (9) can be used as well to construct the solution for a spiral wave that rotates around the hole of the radius R, assuming that the wave touches the boundary of this hole. If the boundary is not permeable for diffusion, the wave front should be orthogonal to it. Then the solution will be formally the same as that found above for a spiral wave in absence of a hole. However, the curvature K o at the hole boundary should not furthermore coincide with Kc:; instead K o should be determined by the solution. First, the rotation period T~¢ of the wave can be found from the equation

R =

K~(TR)T R / 2 V "

(16)

After that, substituting T R a s To into (13) and (14) and replacing K,, 2 in these equations by the unknown quantity K 0, we can solve them and find the front curvature at the boundary of the hole. For large holes (R--+ ac) the explicit expression can be found: K,, = 1.287(Vo/DR2) '/3.

(17)

Generally, the wave curvature at the boundary increases when the hole is diminished. But we know that wave curvature cannot be too l a r g e - otherwise the wave will not propagate. Numerical simulations [8, 21] showed that, if the hole is gradually diminished, the wave suddenly separates from its boundary at some critical hole radius R c. After that the free tip is formed that moves around the core circle which is larger than the hole, and thus do not feel its presence. An

387

A.S. Mikhailot', KS. Zykoc /Spiral waces in excitable media

important question is: At which curvature K 0 does this happen? Suppose that a wave moves touching the boundary of the region which is impermeable for diffusion. Then the solution of eqs. (1) and (2) can be obtained if we take the respective solution for a wave in the medium without such a region and cut this solution along the boundary. This implies that the propagation properties of such a wave would be the same as those for the infinitely extended medium. Particularly it means that the wave will break if its curvature at the boundary will locally exceed the critical curvature K~. When this happens at the boundary (provided that the curvature is maximal there), it would look like detaching of the wave from the boundary and formation of a free tip. The above arguments indicate that the spiral wave rotating around a hole will detach from its boundary when K 0 will approach the first critical value K ~ . The rotation period of the wave just before the detachment can be found from (13) and (14) if we replace there (and in the expression for/3) K~2 by Kcl.

3. Determination of the kinematical parameters The kinematical theory of spiral waves that was briefly outlined in section 2 is phenomenological. It includes several parameters which specify the excitable medium and should be found again for any particular model given by (1) and (2). In the present section we discuss how such parameters can be calculated analytically or extracted from simple numerical simulations. First we consider the situation when the distance between the propagating waves is much larger than the length of the recovery tail. Then the kinematical equations include five parameters, i.e. V0, D, Kcl, Kc2 and y. V0 is the propagation velocity of a solitary pulse. It can be easily found in numerical simulations. There are also many analytical calculations

of this property for various excitable media (see e.g. ref. [27]). To determine the coefficient D in the approximate dependence V = V0 - DK of the propagation velocity on the curvature, we derive the exact expression for the dependence V = V ( K ) in the model (1) and (2). The solitary flat wave propagating along the direction x is given by the solution u = .(~).

,, = ~ ' ( ~ : ) .

~ =x-

v,,t

that satisfies the equations

- V0 d u / d ~ = f ( u , t , ) + D . ( d 2 u / d ~ 2 ) , - V0 dz,/d~ = e g ( u , l , ) ,

(18)

with the boundary condition that for I~l --* 00 the medium approaches its stationary rest state. Consider a piece of the wave that represents a large segment of the circle of radius R = 1/K. Let us go into the polar coordinate system (r, ¢) whose center lies in the center of this circle. Since the wave is circular, there will be no dependence on the polar angle ¢ and, hence, eqs. (1) and (2) will take the form

ti = f ( u , ~,) + ( D , / r ) ( O u / O r ) + D,(OZu/O~2), i, = eg(u,L').

(19)

Suppose further that the curvature radius R is very large, much larger than the width of the wave. Since the activator density u varies significantly only within the width of the wave, we can approximately replace r by R in the second term of the first equation. After this we can go to the local coordinate system which moves with the same velocity V as the wave. In this system (19) takes the form

-( V + D.K ) du/d~ = f( u,v) + Du(d2u/dL2), - V d L ' / d ~ = e g ( u , c).

(20)

388

A.S. Mikhailoc, V..S. Z v k o v / S p i r a l waves in excitabh" media

Here we have already assumed that the wave profile is stationary in the moving coordinate system. Comparing (20) with (18), we see that the two sets of equations are very similar. They can be made identical if we slightly transform (20). Let us introduce the modified p a r a m e t e r e* = e ( 1 + D , , K / V ) .

(21)

Then (20) will be written as

-(V + D,,K)du/d~=f(u,v)

+ D,,(d u / d ~ " ) ,

- ( V + D,,K) d v / d s ~ = e * g ( u , v ) .

(22)

These equations formally describe the propagation of a solitary pulse in the same medium but with a different excitability (remember that the inverse of e specifies the excitability of the medium). Therefore the sum V + D . K will be equal to the propagation velocity Vo(e*) of a solitary pulse in the medium with the p a r a m e t e r e*, i.e.

V + D , , K = Vo(e* ) .

(23)

According to (3), this velocity is

Vo( e* ) = <,,,( 1 - Xe* ).

(24)

Substituting (21) and (23) into (24) we find the equation for V as a function of K:

can be approximated by the linear function V = Vo - DK with D = D, 1 + Xe 1 -Xe

(27)

The determination of the first critical curvature Kc~ involves a certain ambiguity. It is defined by the condition that the propagation of a curved wave (or of a periodic curved wavetrain) becomes unstable at K = Kc~ (therefore the wave breaks if locally its curvature K exceeds K~). The direct use of this definition in simulations would require, however, too much computation. The typical dependence of the propagation velocity V on the wave curvature K has the form shown in fig. 5. No steady propagation is possible above a certain value of K~ at which the derivative d V / d K goes to infinity. Evidently, it might very well happen that this steady propagation becomes unstable even earlier, at some lower value of the curvature. However, since no solution exists for K > K~ at all, the first critical curvature cannot exceed Ks. Hence, Ks represents an upper estimate for Kc~. Below we simply take K~ as the first critical curvature Kcl. Consider eqs. (21) and (23). For a given parameter e, they determine the propagation velocity V and the renormalized p a r a m e t e r e* as some functions of the curvature K. It is convenient to rewrite them in the form 6 ~

V+D,K=

~,o[1 - x e ( 1 + D , K / V ) ] .

(25)

e* - t D ' ' K = G , ( e * ) , t

D,,K.

(28) (29)

It has the solution

V = (~,,,/2) {( 1 - Xe - D,,K/Voo) + [(1 - Xe - D,K/Voo) 2 - 4x~D,K/V,,,,]

'/). (26)

The dependence (26) has two branches; the lower of them corresponds to the unstable solutions. In the region of small curvatures the upper branch

The first of these equations determines the renormalized p a r a m e t e r e* as a function of K. Its substitution into (29) yields V = V(K). Since the dependence of V0 on e* is known (cf. fig. 2), the solutions of (28) can be easily analyzed by a graphic construction. Suppose that we have plotted on the same graph the two functions which are the left and the right sides of eq. (28). Then the intersection points give the

389

A.S. Mikhailoc, IAS. Zykoc /Spiral wal'es in excitable media Vo

Let us construct the approximate analytical solutions of (31) in the two limiting cases. Suppose that the function Vo(e*) can be approximated by the linear dependence (3). Then (31) reduces to

I I

I I K--Kc

(~*/~)(~*

- ~) = (1 - x~*)/x.

(32)

1o0 I

cl

Its solution is I

e* = ( e / x ) '/2

I

(33)

I i

Substitution of this expression into (28) yields o.2

~

o.~

~;

~ gcl ~ (roo/Ou)[1

Fig. 7. The graphical solution of(28).

roots of this equation (fig. 7). Generally, there are two such roots e T and e~ that correspond, according to (29), to the upper and the lower branches of the dependence V = V(K). If we increase K, the two roots move towards each other. They merge at some value of K = K ~ above which no solution is possible. Hence, K~I can be determined from the condition that at K = Kc~ the curves of the two functions in fig. 7 touch in a single point e*. This implies that at such a point both the values of the functions and of their first derivatives coincide. Therefore we have an additional equation

e

(e* - e ) 2 D " K -

dVo

dK"

(30)

- (X6)1/2] 2

(34)

Note that (33) and (34) are valid only for sufficiently small values of the parameter e. Consider further the opposite limiting case where e is close to the value ecj, above which no flat waves can propagate in the given medium. Let us introduce the difference 8e = ec~ - e. Near the point e* = ec~ the upper branch of the dependence Vo(e*) can be then approximated by

Vo=Vc,[l+e~(Se*)l/2],

(35)

where a is some dimensionless coefficient and V~.~ is the value of the propagation velocity at e = ecl. Substituting (35) into (31) and solving the resulting equations we find E* = E c l -- (O~'2/4) ~E 2.

(36)

Dividing (28) by (30) we find

~. --(~* -~) e

v0(~*) dVo/de* "

This yields the estimate for Kcl in the considered region:

(3~) Kc, = ( V c l / O u ) ( S g / e c l ) [ 1

This equation for e* can be solved if we know the function Vo(e* ).

+ (o<2/4)8~'].

(37)

If the deviation e is very small, we can keep only

390

A.S. Mikhailot', V.S. Zvkot' / Spiral waces in t~citahh, media

Re2 is much l a r g e r than the width of the wave. we can a p p r o x i m a t e l y d e s c r i b e the m o t i o n of such b r o k e n circular wave by the e q u a t i o n s

-(V+D,,K~,)(au/a~) = f ( u , c) q.O

+D,,(a~./a~ ~) + D,,(a~./ay~), - v a t , / a ~ = e-g(u, r ) ,

o

o.4

o°a

C

Fig. 8. Dependence of the critical curvatures K~I and Kc2 on e for the model (1), (2) and (47)-(50).

the linear t e r m in (37) which yields Kcl = ( V c , / D u ) ( s c ,

-- e ) / a c , .

(38)

T h u s the first critical c u r v a t u r e vanishes at e = %p T h e typical d e p e n d e n c e o f the first critical curvature K d on e is shown in fig. 8. T h e s e c o n d critical c u r v a t u r e K~2 can be estim a t e d in a similar m a n n e r . W e know that at e = %2 the b r o k e n flat wave p r o p a g a t e s in the m e d i u m with the c o n s t a n t n o r m a l velocity V~.2, n e i t h e r growing n o r c o n t r a c t i n g in the t a n g e n t direction (fig. 6). In the c o o r d i n a t e system w h e r e such a wave is at rest, it is d e s c r i b e d by the equations

- E2 au/a~ = f ( . , t:) + D,,(O2u/O~ 2) + D , ( a 2 u / ~ ) y 2 ) , - V~2 a v / a ~ = e c 2 g ( u , v ) ,

(39)

(40)

w h e r e ~ = r - Vt a n d y = R~.2q~. C o m p a r i n g (39) and (40), we again n o t e that thc p r o b l e m for the curved wave can be effectively r e d u c e d to the p r o b l e m for the p r o p a g a t i o n of a flat wave at a velocity V + D,,Kc2 = V0(e*) in the m e d i u m with a r e n o r m a l i z e d p a r a m e t e r ~* = t(1 + D , , K ~ . J V ) . Since for K = Kc2 the b r o k e n wave d o e s not grow or r e t r a c t at its tip, the r e n o r m a l i z e d par a m e t e r e* should be equal to %2 (then the velocity ~ ( e - * ) will be equal to ~ 2 ) . This condition yields the e q u a t i o n s %2 = e(1 + D , , K ~ 2 / V

),

Vc2 = V + O u K c 2 .

(41) (42)

Solving (42), we find the s e c o n d critical c u r v a t u r e Kc2 as a function of e:

Kc2 = ( V~2/D,,)( 1 - ~/~c2)"

(43)

W e see that Kc2 vanishes at e = %2. T h e last k i n e m a t i c a l p a r a m e t e r 7 can be also c s t i m a t e d . A c c o r d i n g to (4), the t a n g e n t velocity of growth at the tip of a b r o k e n flat wave is G o = yK~2. H e n c e , by m e a s u r i n g this velocity G 0 in a n u m e r i c a l e x p e r i m e n t , we can easily d c t c r mine

w h e r e s~ = x - Vcet. S u p p o s e f u r t h e r that we have a b r o k e n circular wave with the critical c u r v a t u r e K~2 (i.e. with the radius R,. 2 = 1/K~2). By the definition of Kce, such a wave w o u l d not grow or r e t r a c t at its free end. Let us i n t r o d u c e the p o l a r c o o r d i n a t e system ( r , ~#), so that its c e n t e r coincides with the c e n t e r of the circle. A s s u m i n g that t h e c u r v a t u r e r a d i u s

7 = Go/Kc2,

(44)

p r o v i d e d that the s e c o n d critical c u r v a t u r e Kc2 is a l r e a d y known. T h e above p r o c e d u r e s for d e t e r m i n a t i o n of the p a r a m e t e r s Kc~ a n d Kc2 can be g e n e r a l i z e d to the case when we have a train of the curved

391

A.S. Mikhaikw, V.S. Zykoc / Spiral waces in excitable media

waves with the period T comparable to the recovery duration 7],. The propagating train of flat waves is again described by the solution of eqs. (18), but this time with the periodic boundary conditions:

4. Numerical simulations

In numerical simulations we used a particular realization of model (1), (2) with the functions

f(u,c) =f(u)-c, u(~+L)=u((),

v(s~ + L) = t:(s~),

(45)

where L is the spatial period of the train. These equations determine the propagation velocity Vo(e, L) of the train. If we have a train of curved waves, we can again reduce (cf. (22)) the problem to propagation of a train of flat waves with the same spatial period L in the medium with the renormalized p a r a m e t e r e*. Hence we find again the equations (21) and (23) where, however, we should take V0 = V0(F*, L). Solution of these equations would yield a certain dependence V = V(K) for given e and L. The first critical curvature Kcl(e, L) for a curved wavetrain with the spatial period L can be found from the condition that d V / d K = ~ at K=Kcl(e, L). The temporal period T can be expressed in terms of the spatial period as T = L / V , where V is the propagation velocity of the curved wavetrain that should be also determined from (21) and (23). In combination with the known dependence Kcj(e , L) this would allow to determine K d as a function of the temporal period T. Similar arguments can be used to find the temporal period dependence of the second critical curvature Kc2. NOW we should replace %2 and Vc2 in (43) by ec2(L) , the value of e at which the periodic train of flat broken waves with free tips and propagates without growth or retraction at the wave tips, and by Vc2(L) which is the respective propagation velocity of the train with the spatial period L. The temporal period T can be expressed in terms of the spatial period as

T=L/V=

[%2(L)/e][L/Vc2(L)].

(46)

The dependences ec2(L ) and Vc2(L) can be found from the data of numerical simulations.

f(u) = -k,u,

(47)

~O',

= k r (u - a),

¢r
=k2(1 -u),

l --o'
g(u,v) =h(u,t'), =k~h(u,t,), h(u,L') = k ~ u - t ' .

(48)

for h(u,~') >_0, for h ( u , c ) < O,

(49)

(50)

The following values of the parameters are always taken: k f = 1.7, k g = 2, a =0.1, ~r=0.01. The parameters k~ and k 2 are chosen in such a way that the function f(u) is continuous at u = ¢r and u = l - ~ r ; this yields k l = k f ( a - o - ) / o - = 168.3 and k 2 = k f(1 - cr - a)/cr = 151.3 The diffusion constant is D u = 1. The same model (47)-(50) was used in simulations that are described in book [19]. The properties of the excitable medium described by (1) and (2) with the functions (47)-(50) can be controlled by changing the parameters e and k~,. As it was already noted, the inverse of e specifies the medium excitability. The coefficient k~ allows to control the ratio T~,/T,. of the excitation and recovery durations. During the excitation stage (which corresponds to the top of the pulse, see fig. 3) the activator variable t, increases and hence we should take g(u, v) = h(u, v). On the other hand, in the recovery tail v goes down and hence g ( u , v ) = k,,h(u, c). Therefore, the larger the coefficient k~ is, the shorter the recovery durations T,. This means that the ratio T i T , is inversely proportional to k~. Extensive simulations of spiral waves in the model (1), (2) with the functions (47)-(50) showed [33, 35] that the meandering regimes are observed only for sufficiently high excitabilities (e < 0.3) and provided that k~ < 1.5. Since in the present

392

A.S. Mikhailot

VS. Zykm'

T 100

1 /2 \\ T . 0

002

~¢2 0 . #

~ct

k_.;

Fig. 9. The rotation periods TII of free spiral waves found by direct numerical simulations in the model (1), (2) and (47)-(50) for the different values of the parameter e (the crosses). The rotation periods T E of the spiral waves, attached to the minimal holes, found by the direct numerical simulations in the same model (the black circles). Curve I shows the kinematical theory prediction for T0 when the effects of the recovery tail are neglected. Curve 2 shows the respective prediction for T I. Curve 3 shows the kinematical theory prediction for T I that takes into account the recovery tail effects. The dashed curve shows the computed dependence of Train on c.

paper we consider steadily rotating spirals, in the simulations we chose k~ = 1.5, so that no meandering could a p p e a r (at least for e > 0.025). Fig. 9 shows (by cross marks) the d e p e n d e n c e of the rotation period of spiral waves T 0 on the p a r a m e t e r e which was obtained in ref. [35] by direct numerical simulation of t i m e - d e p e n d e n t equations (1), (2) with the functions (47)-(50). In the same figure we show (by black circles) the c o m p u t e d d e p e n d e n c e on e of the rotation period of a spiral wave that rotates attached to the smallest possible hole (see the discussion in section 2). These d e p e n d e n c e s will be further compared with the predictions of the kinematical theory. If we neglect the effects of the recovery tail, the rotation period of spiral waves T 0 in the kinematical theory is given by the expressions (13)

/Spiral wares

in excitable media

and (14) that include only the parameters V,(e), D(e) and Kc2(e). The c o m p u t e d d e p e n d e n c e V0 = Vo(e) for the given values of the model parameters is presented in fig. 2. Note that it is surprisingly well fitted by the linear function (3) in a very large interval of e, even close to ec2. Such a linear approximation has V0o = 1.5 and X = 0.55. Using this value X and applying (27), we can calculate D for any ~- from the considered interval. The d e p e n d e n c e Kc2(e) is given by eq. (43) that includes two parameters, %2 and V,.z. These were d e t e r m i n e d by us from numerical experiments (fig. 6). We started with the initial condition in the form of a broken flat wave and followed its further evolution. By varying the p a r a m e t e r t it was possible to find such its value e - %2 that the broken flat wave did not grow or retract. This value was found to be ec2 = 0.388. The normal propagation velocity of the flat wave which is observed at e = ec2 is Vc2 = 1.24. Substituting the c o m p u t e d values of the kinematical parameters into eqs. (13) and (14), we obtain the d e p e n d e n c e of the rotation period T~j on e which is shown by curve 1 in fig. 9. Note that it is calculated by neglecting the recovery tail effects. Within the kinematical theory it is also possible to calculate the rotation period of the spiral wave that rotates a r o u n d a circular hole. As it was emphasized in section 2, for any given medium there is a ccrtain minimal size of thc hole, such that for the smaller holes the spiral wave tip becomes d e t a c h e d from the hole boundary. The respective rotation period T~ is a characteristic property of the m e d i u m itself, as universal as the rotation period of spiral waves in absence of any holes. At a first glance it might seem that, to find T l, it is sufficient to replace Kc2 by Kc~ in the expressions (13) and (14). However, the situation is more complicated here. Since the first critical curvature is determined by the condition d V / d K = ~c, the linear approximation V = Vo DK cannot hold until this value of K. Instead,

A.S. Mikhailoc, KS. Zykoc / Spiral wal,es in excitable media

the actual curvature dependence V = V(K) should be substituted into eq. (9). Then this equation can be numerically integrated and, by requiring boundary conditions (7) and (8), one can determine the rotation frequency ~o. Above it was noted (see fig. 2) that the dependence of V on e for the considered model can be very well fitted by the linear dependence (3) almost until eel. Therefore, for this model the dependence of V on K can be described by (26) in a wide range of values of e. The numerical integration of (9) with this particular function V(K) was carried out in refs. [20, 21]. It was found that the rotation period T, of the spiral wave with K ( 0 ) = Kc, can be approximately described by

T,

2-rrD v Zu(/3, ) ,

(51)

where v = 0.51/3, - 0.02/3 2 - 0.16/33

(52)

and /31 = DKc,/V o. Curve 2 in fig. 9 shows the rotation period T 1 as a function of e, computed using (51) and (52). The effects of the recovery tail in the circulation of spiral waves can be neglected only while their rotation periods are much larger than Tm~,,, the minimal temporal period of one-dimensional pulse propagation. The dashed curve in fig. 9 shows the actual dependence of Tmin that was computed for our model. We see that the rotation period T0 of the spiral waves is larger than Tram only for sufficiently large values of e (i.e. close to ec2). At the smaller values of e we should take into account the dependence of V0 and Kc2 on the rotation period that can be obtained by the numerical simulations. Then, to find the rotation period T 0 of a spiral wave, one should take the computed dependences Vo(T) and Kcz(T) for a given value of e, substitute them into (13) and (14) and to solve numerically the resulting equations.

393

Kc

L

0

Tmi n

20

40

T

Fig. 10. The dependence of the first critical curvature Kcl on the temporal period T in the model (1), (2) and (47)-(50): e = 0.25.

At present, we have not yet realized this program of calculations. However, we can report the numerical results for a related problem which consists in calculation of the rotation period T, of a spiral wave attached to the minimal hole. To realize the programme, we started with the computation of the dependences V0(T) (see fig. 4) and ecl(T) for our model. Using them, we computed KcI(T) in accordance with (34) for different values of e (fig. 10). Then these functions were substituted into (51) and (52) which thus became the equations for the determination of T,. Curve 3 in fig. 9 shows the dependence of T~ on e which was obtained by their solution. The main purpose of this p a p e r is to compare the predictions of the kinematical theory with the data of numerical simulations for r e a c t i o n diffusion systems. Since we now determined all principal kinematical parameters for the reaction-diffusion model, given by eqs. (1), (2) and (47)-(50), and obtained the quantitative estimates for the rotation period of spiral waves within this theory, such a comparison can be carried out. The values of the rotation period T0 of spiral waves, computed for the original full r e a c t i o n diffusion equations, are shown by crosses in fig. 9. We can see that in the region of small e (which correspond to high excitabilities) they are close to Tm~n, the minimal temporal period of stable one-

394

A.S. Mikhailor, U.S. Zykoc / Spiral waces in excitable media

dimensional pulse propagation. This means that the rotation regime of spiral waves in this region is essentially controlled by recovery tail effects. For the smaller ratios T,/T~ we would probably find meandering in this region of the parameters e. Since fl)r small e the estimate T,,~i,~~e" ~/2 holds (see, e.g. refs. [19, 2(1, 36]), we can roughly expect that the rotation period of spiral waves has in this region the same asymptotical behavior, provided that no meandering is observed. The rotation period T~ diverges as well near the critical value e~2, where the fiat broken wave can propagate without growth or contraction. Here its computed values agree with the estimate of the kinematical theory (curve 1) obtained when we neglect the recovery tail effects. According to (13) and (14), the rotation period To varies as T o ~ Kc23/2 for very small critical curvatures K~2. As follows from (43), we have K~ ~ (e,. 2 - e ) . Therefore, the rotation period of spiral waves should behave as

To_(ec2_e

) ~/2

(53)

for the right branch of the U-shaped dependence, near the critical value ec2. It is interesting to note that both estimates for T0, obtained either by neglecting the curvature effects (and thus putting T0 = Tram) or by neglecting the recovery tail influence (curve 1), predict somewhat lower values of the rotation period than those actually observed in the numerical simulations. It can be expected that, when these two factors are simultaneously taken into account, the agreement would be even better. In addition to the rotation period T0 of free spiral waves, any excitable medium possesses yet another characteristic property which is the rotation period T~ of spiral waves attached to the minimal possible hole (when the hole radius is further decreased, the rotating wave detaches). The values of T~ which were obtained by numerical simulation of the attached spiral waves for the full reaction-diffusion model are shown by black circles in fig. 9. We see that for the small parame-

t e r s e they almost coincide with the values of the rotation period T0 of free spiral waves. This confirms that the curvature effects are not significant in this region. However, the difference between T~ and T0 becomes large for the higher values of parameter e. We can see from fig. 9 that the attached spiral wave exists even in the interval e,.2 < e < e~.t, where no free spiral waves can be observed. The estimate for iV 1 which was obtained within the kinematical theory, by applying (51) and (52) and neglecting the recovery tail effects, is shown by curve 2 in fig. 9. It can be noted that it has the correct asymptotical behavior near e =ec~ but predicts the values of T t which are smaller than those supplied by numerical simulations. For the spiral waves attached to the minimal hole wc were able to perform the complete kinematical calculation of the rotation period T~, taking into account both the curvature and the recovery tail effects. The results of such a calculation are shown by curve 3 in fig. 9. A good fit with the data of the numerical simulations of T~ with full time-dependent reaction-diffusion equations was thus achieved.

5. Discussion In this paper we performed two numerical tests of the kinematical theory. They consisted of calculation of rotation frequences of free spiral waves and of waves, that rotate around "minimal" holes, for different values of the excitability. These calculations, which used no adjustable parameters, produced good agreement with the data of numerical simulations of spiral waves in the respective reaction-diffusion model. The agrecment is especially close in the case of spiral waves attached to minimal holes. It would be premature to consider these results as the final verification of the kinematical theory. The next task might consist of quantitative comparison between various nonstationary effects (such as the process of relaxation to steady rota-

A.S. Mikhailov, V.S. Zykor /Spiral waves in ~<~'citabh"media

tion or meandering, drift and resonance etc.) and numerical simulations of the same effects in reaction-diffusion models. This also puts forward a task of comparing the kinematical theory with the experimental data, for the cardial tissue and the Belousov-Zhabotinsky reaction. Obviously, under the experimental conditions it is not possible to vary the characteristic time ratio e. However, one can easily control the excitability of the medium (by changing the chemical composition of the nutrition solution in the experiments with the isolated cardial tissue specimens; by varying the chemical composition of the BZ solution or by changing the illumination when the photosensitive modification of the BZ reaction is used). Since, under fixed other parameters, the effective excitability is inversely proportional to ~, the U-shaped dependence of the rotation period of spiral waves, which is mirror-reflected in respect to that shown in fig. 9, is formed. The sharp increase of the rotation period T0 in the region of lower excitabilities q is explained by vanishing of the second critical curvature Kc2 at some critical excitability q~2, i.e. by the dependence K~2 ~ ( q - q~2)- Since, according to (13) and (14), we have T o ~ K~2~/2 for very small critical curvatures, the kinematical theory predicts the universal power law

T0~ (q -qc2) 3/2,

(55)

where q is any p a r a m e t e r of the medium which is related to its excitability; the asymptotical dependence (55) should hold only near the critical value qc2 below which the free spiral waves are not observed. The above results were derived under the assumption that the "diffusion constant of the inhibitor" is zero. This is true for the cardial tissue where (1), (2) are actually the F i t z H u g h - N a g u m o equations and the "inhibitor concentration" t' is the local potassium ion conductivity of the cell membrane. However, it is not directly applicable to the BZ reaction where all the reacting compo-

395

nents diffuse. Nonetheless, as we show in the appendix, a similar procedure for the determination of the kinematical parameters can be suggested in the case with two nonvanishing diffusion constants. Provided that the diffusion constant of the activator is higher than that of the inhibitor, one can expect a similar behavior of the rotation period as a function of the medium excitability. The detailed study of this problem will be a subject of a separate publication.

Acknowledgements The authors acknowledge stimulative discussions with S. Muller, H.L. Swinney and A.T. Winfree. One of the authors (A.S.M.) wants to express his gratitude to H. Haken for the hospitality provided during his long stay in Stuttgart.

Appendix In the present p a p e r we discussed only the situation for models where the inhibitor does not diffuse. This assumption holds, for instance, in the cardial tissue. However, in the chemical medium with the Belousov-Zhabotinsky reaction both the activator and the inhibitor diffuse. Their diffusion constants D , and D, are not precisely known. It can be said only that they must be of the same order of magnitude, because they are principally determined by the molecular weight of a reacting component. The above analysis can be repeated ['or the case of arbitrary diffusion constants of the activator and the inhibitor. This work is now in progress and will be reported elsewhere. In this appendix we only give some basic relationships which allow to determine the kinematical parameters in this more general case. Now we consider a model described by the equations gt = f ( u , c )

t"=eg(u,t')

+ D,,Au, + D At'.

(A.1) (A.2)

396

A.S. Mikhailov, K S . Z y k o v / Spiral waves in excitab& media

The steady propagation of a solitary fiat wave in this model is described by - Vo d u / d ~ = f ( u , , , )

+ D,,(d2u/d~2),

- V0 dv/d~: = e g ( u , v ) + D,(d2v/d~2).

(A.3)

The propagation velocity V0 of a fiat wave is a certain function of e, D,, and D,; it can be found for any particular model. Consider a curved wave in the same medium. If its curvature K is sufficiently small, the steady propagation of such a wave will be approximately described by the equations (cf. (20))

-(V+D,,K)du/d~ =f(u,v)

v 2 + D,(d-u/d~ ),

-(V + O,.K)dc/d~ 2

=eg(u,v)+D,(d

v/d~:~).

(A.4)

Introducing

V + D,, K ~-*=eV+D, K,

The second critical curvature K~., can be obtained if we put e* = e~2 and Vo(e*, D,,, D,* ) = V~, in (A.5) and (A.7). This yields

V + D,, K D,*=D, v + D , K

(A.5)

we can cast them in the form

Vc2

K"2

6c2-6

D,,- D,.

e~2

(A.8)

However, (A.8) is not a final result because V~2 and e~_~ here should be taken for a medium with the renormalized inhibitor diffusion constant

D,* = ( G.2/e ) D,.

(A.9)

Provided we know the dependences ec: = ~~2( D, ) and Vc2 = Vc2(D, ), (A.8) and (A.9) give the system of equations which should be solved to determine Kc2. The situation is simpler when v is close to Vc2(D, ). Then, as it follows from (A.9), the renormalization of the inhibitor diffusion constant is negligible and, therefore, V~.2 and e~2 in (A.8) can be taken for the original diffusion constant D,. In this appendix we do not proceed further with the analysis of the derived relationships. Still it can be noted that, as it could be expected, a special behavior is found when the two diffusion constants D, and D, are very close to each other.

- ( V + D.K ) du/d~

=f(u,v)

+D,(d

2

References

v u/d~-),

-(V + D,K)dv/d~ = ~*g(u,v) + D*(dev/d~2).

(A.6)

Comparing (A.6) with (A.4) we find that

V+ D,,K = Vo(e*, D,,, D,*).

(A.7)

If the dependence V0 = Vo(e, D,,, D,) is known, eqs. (A.5) and (A.7) can be solved to obtain the curvature dependence of the propagation velocity V= V(K). Expanding V(K) in powers of K and keeping only the linear term, one can determine the coefficient D in the linear approximation V= V¢~- DK and find the first critical curvature Kcl from the condition that d V / d K = ~.

[1] A.T. Winfree, Spiral waves of chemical activity, Science 175 (1972) 634 686. [2] W. Jahnke, W.E. Skaggs and A.T. Winfree, Chemical vortex dynamics in the Belousov-Zhabotinsky reaction and in the two-variable Oregonator model J. Phys. Chem. 95 (1989) 74(I-749. [3] S.C. Muller, T. Plesser and B. Hess. The structure of the core of the spiral wave in the Belousov Zhabotinsky reagent, Science 2311 (1985) 661. [4] S.C. Muller, T. Plesser and B. Itess, Two-dimensional spectrophotometry and spiral wave propagation in the Belousov-Zhabotinsky reaction. I1. Geometric and kinematic patterns, Physica D 24 (1987) 87-96. [5] G.S. Skinner and H.L. Swinney, Periodic to quasiperiodic transition of chemical spiral rotation, Physica D 4;4 (1990) 1 16. [6] W.Y. Tam, W. Horsthemke, Z. Noszticzius and tl.L. Swinney, Spiral waves in continuously fed unstirred chemical reactor, J. Chem. Phys. 88 (199171) 3395 3439.

A.S. Mikhailoc, V.S. Z)'kot' / Spiral waces in excitable media [7] A.M. Pertsov, A.V. Panfilov and E.A. Ermakova, Numerical simulation of spiral waves in active media. Physica D 14 (1984) 117-125. [8] A. Karma, Meandering transition in two-dimensional excitable media, Phys. Rev. Lett. 65 (1990) 2824-2827. [9] D. Barkley, A model for fast computer simulation of waves in excitable media, Physica D, submitted. [10] E. Lugosi, Analysis of meandering in Zykov kinetics, Physica D 411 (1989) 331-337. [I 1] W. Jahnke and A.T. Winfree, A survey of spiral waves in the Oregonator model, Int. J. Bifurc. Chaos, submitted. [12] J.P. Keener and J.J. Tyson, Spiral waves in the Be[ousov-Zhabotinsky reaction, Physica D 21 (1986) 3110 324. [13] J.P. Keener and J.J. Tyson, Spiral waves in a model of miocardium, Pbysica D 29 (1987) 215-222. [14] E. Meron and P. Pelc6, Model for spiral wave formation in excitable media, Phys. Rev. Lett. 60 (1988) 1880-1883. [I5] E. Meron, Nonlocal effects in spiral waves, Phys. Rev. Lett. 63 (1989) 684-687. [16] D.A. Kessler and H. Levine, Physica D 49 (19911 90 97. [17] P. Pelc~ and J. Sun, Wave fronts interaction in steadily rotating spirals, Physica D 48 (1991) 353-366. [18] N. Wiener and A. Rosenbluetb, The mathematical formulation of the problem of conduction of impulses in a network of connected excitable elements, specifically in cardiac muscle. Arch. Inst. Cardiol. Mexico 16 (19461 pp. 205-265. [19] V.S. Zykov, Simulation of Wave Processes in Excitable Media (Nauka, Moscow, 19841 (English translation: Manchester Univ. Press, Manchester, 1987). [20] A.S. Mikhailov, Foundations of Synergetics. I. Distributed Active Systems (Springer, Berlin, 19901. [21] V.A. Davydov and A.S. Mikhailov, Spiral waves in distributed active media, in: Nonlinear Waves. Structures and Bifurcations, ed. A.V. Gaponov-Grekhov, M.I. Rabinovich (Nauka, Moscow, 1987) pp. 261 279. [22] P.K. Brazhnik, V.A. Davydov and A.S. Mikhailov, Kinematical approach to description of autowave processes in active media. Theor. Math. Phys. 74 (1988) 30/)-306.

397

[23] V.A. Davydov, V.S. Zykov and A.S. Mikhailov, Kinematical theory of autowave patterns in excitable media, in: Nonlinear Waves in Active Media, ed. J. Engelbrecht (Springer, Berlin 1989) p. 38 51. [24] V.S. Zykov, Kinematics of the nonstationary circulation of spiral waves in excitable media, Biofizika 32 (1987) 337-340. [25] V.S. Zykov, Stability conditions for circulation of spiral waves, Biofizika 34 (1989) 1031-1035. [26] V.S. Zykov and O.L. Morozova, Computer simulations and investigation of spiral wave stability, J. Nonlin. Biol., in press. [27] V.A. Davydov, V.S. Zykov, A.S. Mikhailov and P.K. Brazhnik: Drift and resonance of spiral waves in a distributed excitable medium, Radiophys. Q u a n t u m Electron. 31 (1988) 574 582. [28] V.A. Davydov and V.S. Zykov: Spiral waves in anisotropic excitable media, Zh. Eksp. Teor. Fiz. 95 (1989) 139-148. [29] A.Yu. Abramychev, V.A. Davydov and V.S. Zykov, Drift of spiral waves on the nonuniformly curved surfaces, Zh. Eksp. Teor. Fiz. 97 (1990) 1188-1197. [3/I] K.I. Agladze, V.A. Davydov and A.S. Mikbailov, Observation of a spiral-wave resonance in an excitable distributed medium, JETP Lett. 45 (1987) 767-770. [31] V.S. Zykov, Cycloidal circulation of spiral waves in excitable media, Biofizika 31 (1986) 862-865. [32] V.S. Zykov, in: Control of Complex Systems, ed. Z. Tsypkin (Nauka, Moscow, 1975) pp. 59-65. [33] P.K. Brazhnik, V.A. Davydov and A.S. Mikhailov, in: Kinetics and Combustion (Inst. of Chem. Phys. Publications, Chernogolovka. 1986) pp. 39-43. [34] P.K. Brazhnik, V.A. Davydov, V.S. Zykov and A.S. Mikhailov, Vortex rings in distributed excitable media, JETP 66 (1987) 984-990. [35] R.C. Brower, D. Kessler, J. Koplik and H. Levine, Geometrical models of interface evolution, Phys. Rev. A 29 (1984) 1335-1342. [36] A.S. Mikhailov and V.I. Krinsky, Rotating spiral waves in excitable media: the analytical results, Physica D 9 (1983) 346-376.