- Email: [email protected]
The simplified FitzHugh-Nagumo model with action potential duration restitution: Effects on 2D wave propagation
Physica D 50 (1991) 327-340 North-Holland
The simplified FitzHugh-Nagumo model with action potential duration restitution: Effects on 2D wave propagation Boris Y. K o g a n , W a i t e r J. Karplus, Brian S. Billett, Alex T. P a n g Computer Science Department, Unieersity of California, Los Angeles, CA, USA
H r a y r S. K a r a g u e u z i a n and S t e v e n S. K h a n D&ision of Cardiology, Cedars-Sinai Medical Center, Los Angeles, CA, USA Received 30 July 1990 Revised manuscript received 10 December 1990 Accepted 10 December 1990 Communicated by A.T. Winfree
A modification of the simplified FitzHugh-Nagumo (FN) equations is proposed for introducing a residual component of the slow variable, which determines the restitution of action potential duration (APD) also known as the interval-excitation duration relationship. The three-step-wise approximation of e(E)which is widely used in current publications is replaced in a new model by a four-step approximation. This change is used for studying by computer simulation the effects of APD restitution properties independently of the APD and refractory period on 2D wave propagation in an isotropic matrix (made by 128 × 128 nodes). The method for fitting the model to the given experimental restitution data (obtained from myocardial cells) is presented. The computer simulations implemented on a massively parallel computer (Connection Machine) showed at least three important qualitative distinctions in behavior which demonstrate the effect of APD restitution: changes in the speed and wavelength of propagated waves with the period of stimulation, non-stationary propagation of spiral waves, and site-specific induction of spiral waves with premature stimulation not on the tail of the previous wave. Quantitative effects of differing restitution properties are expressed in the size and location of a window of vulnerability in 2D excitable media. These windows are characterized by the appearance of single and double spiral waves in response to premature stimulation applied inside the window. Thus the APD restitution incorporated in the FN model produces a significant effect on the formation and propagation of spiral waves.
1. Introduction
It has long been known for excitable media that if a second excitation is initiated soon after the first one, the duration of the second will be considerably shorter than the first. For biological excitable media the relationship of the interval between two excitations to the duration of the second excitation is known as the interval-duration relation, or action potential duration restitut i o n - hereafter referred to simply as restitution. For example, a cardiac muscle cell requires a time equal to several action potential durations
(APDs) to recover fuliy. If stimulation is applied at relatively short cycle lengths, i.e. while the cell is in a state of partial return to its pre-excited condition, then the response will be different from that obtained when the cell has fully recovered. In the case of restitution a typical outcome is that APs initiated at relatively short diastolic intervals (i.e. the interval between successive APs) will result in action potentials of shorter duration. McAllister, Noble and Tsien [1] showed that restitution in cardiac Purkinje fibers is due to the time-dependent component of potassium current which decays very slowly after the completion of
0167-2789/91/$03.50 © 1991 -Elsevier Science Publishers B.V. (North-Holland)
B.K Kogan et al. /Restitution properties and 2D propagation
328
F(E)
'E ~
I
].dl
1
E
F
.9Is Ill" E
Fig. 1. The piece-wise approximation of functions F(E), f(E), and s(K). tan a = G,: t a n f i = (;~: tan y = ( ; j -/:'th is the threshold potential,
the repolarization phase of the action potential. Introducing a residual component of the slow variable into the FitzHugh-Nagumo model simulates this slow decay of the potassium current. It is evident that the time course of restitution of APD directly affects the degree of dispersion of repolarization, when stimulation is applied during the process of restitution [2]#1. It is known that the greater the dispersion of repolarization, the easier will be the induction of spiral waves (circus movement reentry) [3]. It is therefore essential to incorporate these properties in computer models designed to elucidate the properties of impulse propagation. Over the last 30 years, many valuable contributions have been made in wave propagation studies using axiomatic models [4, 5], and the simplified FitzHugh-Nagumo (FN) equations (see refs. [6-9]). However, the parameters of these models, particularly the function s ( E ) in the simplified FN model, provide limited, if at all, restitution properties. Here we propose a method for incorporating restitution properties (residual component of slow variable) into the simplified FN model. Then we present the results of a #1For the role of A P D restitution curve features on the a p p e a r a n c e of cardiac dysrhytbmias and possible suppression by appropriate pharmacological manipulation see ref. [26].
computer simulation of wave propagation in a 2D piece of an isotropic biological homogeneous excitable media (128 × 128 nodes) fitted to experi mental physiological restitution properties of cardiac tissue [2]. These results show the significant qualitative and quantitative distinctions from the simplified FN model without restitution properties. All computer simulation studies were carried out using a massively parallel computer system which is particularly suited to this class of problems, the Connection Machine (CM). For a brief description see appendix A, and for more details see ref. [10].
2. The FitzHugh-Nagumo simplified equations and their properties The basic FN simplified equations when reduced to dimensionless form [6] are aE
a-7- = A E + F ( E )
- I + l~tm,,
aat/ = e(E) I f ( E ) - I],
(1)
(2)
where E is the fast variable (membrane potential displacement between the interior and exterior of
B.Y Kogan et al. / Restitution properties" and 2D propagation
the cell); I the slow variable (generalized outward current), A = a2/0x 2 + ~ 2 / ~ y 2 is the two-dimensional Laplacian operator, F(E) the c u r r e n t voltage characteristic of the fast inward current, f(E) the current-voltage characteristic of the slow outward current and e(E) a small parameter (inversely proportional to the time constant of the slow outward current). The initial and boundary conditions used in the simulation are
E(x,y,O) =l(x,y,O)
=0,
The piece-wise approximation of the functions F ( E ) , f ( E ) , and e(E) is shown in fig. 1. The widely used set of parameters (see refs. [6-9, 11-13]) hereafter referred to as the standard set is: G, = 1,
Gf =
e=e= =0.5
1,
G r = 30,
Eth =
0.16
if 0.00 < E < 0.01,
= e 2=0.01
if0.01
= e 3 = 0.5
if E > 0.95.
(3)
We define an action potential as the time course of the membrane potential which is above the threshold potential Etn. This action potential is
329
initiated by a proper stimulation and formed by the activation and inactivation of inward and outward currents. The relationships between the action potential duration (APD), refractory period (R), and the model parameters G~, Gf, and e are presented in ref. [6] for the model without diffusion (point model). These dependencies reflect qualitatively the essential properties of heart muscle cells, but do not correctly express the restitution properties. The normalized restitution curve for the FN simplified model with standard parameters is compared to normalized data obtained in physiological experiments in fig. 2. The comparison shows the marked disagreement of the model and experimental data. We introduced the normalized form of the restitution curve to facilitate the comparison of computer simulation data obtained in dimensionless form and dimensional data from physiological experiments. Usually, the restitution curve is presented as a function: APD = f(DI), where APD designates the duration of the action potential caused by the second stimulus applied after the specified diastolic interval (DI). The normalized representation of the restitution curve (RC) has the form APD2/APD l = h(DI/APD~),
1.1
1.0
0.9
0.8
*
0.7 0
I
I
1
2
b)
I
DI/DI
3
4
Fig. 2. The normalized restitution curve for (a) results of computer simulation performed by the author using the FitzHugh-Nagumo simplified model with set of parameters [6], and (b) canine heart muscle tissue in normal conditions (control), obtained from experimental data and normalized for comparison [2].
B.Y. Kogan et al. ,/Restitution properties and 2D propagation
330
1.2
r
:
i
i
i E
1.0
i
j
....
i
...
0.8,~
0.6-~
/
0.4 0.2
/ .Y
j
i
/
0.0 i
I 0 •2
o
2op
I aoo
4oo
5oo
t = 0.234; APD2/APD
I
600
dine Fig. 3. T h e E(t) a n d l(t) f o r s h o r t d i a s t o l i c i n t e r v a l . D I / A P D
where A P D l is the duration of the A P obtained after the full completion of the recovery processes. H e r e a f t e r A P D 1 = D 1 and A P D 2 = D 2. T h e transients of E and 1 in a cardiac cycle with a small D I are shown in fig. 3 for the standard p a r a m e t e r s of the F N simplified model. O n e can see that the outward current I reaches zero just after the A P returns to the resting potential. Thus, a stimulus applied after a short D I p r o d u c e s a subsequent A P approximately of the same duration as the previous one, because the time course of the restitution is very fast (almost instantaneous). We therefore modified the characteristics of the transient outward current I ( t ) , by slowing down its time course of inactivation during the interval of negative A P values.
3. S o m e a p p r o a c h e s to restitution curve fitting O n e m e t h o d of modifying the transient outward current I ( t ) is to use two values of e dependent on the sign of d l / d t : e(dI/dt)
=e I
ifdI/dt>O,
=ke I
ifdl/dt<0.
(4)
=
1 (in t h e m o d e l o f ref. [6]).
This a p p r o a c h was p r o p o s e d by Zykov [14] to change the refractory period independently of the A P duration in eqs. (1), (2) with the standard set of parameters. The increase in the refractory period to longer than the duration of the A P contradicts the normal properties of heart muscle tissue. C h a n g i n g the p a r a m e t e r k also changes the restitution properties simultaneously with the refractory period. T h e use of only one p a r a m e t e r ( k ) does not permit the required shaping of l ( t ) to fit restitution properties, nor does it permit the separate investigation of the effect of refractory and restitution properties on impulse propagation. O u r a p p r o a c h is to use a more complicated logical function to control the changes of the small p a r a m e t e r e. Namely e = eI
if E < 0.01 and d l / d t
> 11.
= e2
if E > 0.01 and d l / d t
> O,
= ~3
if 1 > Im~,l and d l / d t
= e 4 = ke 2
otherwise.
This mation The satisfy
<_ O,
(5)
leads to a piece-wise exponential approxiof the current function l ( t ) . values for e l, e 2, and e 3 are chosen to the requirements of the given action po-
B. K Kogan et al. / Restitution properties and 2D propagation
the computer model, it is necessary to choose E4 such that
tential duration and its refractory period, e 4 is chosen to meet the required restitution properties. The proper values of e 4 and Imin can be matched to the experimental restitution curve by computer simulation trial and error or by a direct analytical approach. The normal restitution curve is commonly approximated by the formula [15]
D = Do( 1 - e-t~v).
,E"4 = / 3 ( O ~ ) e x p ( Doc)sim .
(6)
D/D~) =/3T.
(8)
This formula is based on the well known physiological fact that the time constant of the restitution curve (/3) is similar to the time constant of the variable x 1 which determines the potassium current ixl at negative potentials [16]. The residual component of the slow variable in the FN model approximates this potassium current. The scale factor rn (Ooc)exp/(O~)sim is used to transform the dimensional /3 to the dimensionless e 4. The results of restitution curve (RC) fitting for heart muscle tissues under normal conditions (control) and with the tissue exposed to quinidine are presented in figs. 5a and 5b, respectively. In both cases the close agreement with experimental data was obtained by making E4 =/3m and choosing an appropriate value for Imin' The maximum and rms error in the region 0.5 < D I / A P D I < 3.0 for the control data are respectively: 0.03 and 0.0066, and for the quinidine are: 0.038 and 0.010. It was found that a small decrease of the parameter G f improved the fitting of the quinidine experimental data. This approach also gives a slight
Here D = A P D 2 and D~ = the plateau A P D after an infinitely long DI; /3 is a constant, and T = APD~ + DI. From (6) follows
-ln(1 -
331
=
(7)
Thus the relationship between - I n ( 1 - D / D ~ ) and T is linear, and the experimental data can be represented by a straight line, which can be obtained by the least-squares method of approximation. The results of this approximation are shown in fig. 4 for two sets of experimental data [2] from normal tissue and tissue intoxicated by quinidine. The slope of the line obtained in this manner directly shows the values of/3. To reproduce the restitution properties of the experimental data in
J
_J
J
..l y
J
,...--
J
J
J
2
.J
f
1
2O0
300 y = 0.0047x y =
0.0055x
400 + 0,14 -
1.08
500
R^2 = 0.973
D-----
R^2 = 0.978
"=
600 700 a) control data b) quinidine data
T
800
Fig. 4. Determination of e 4 from experimental data: (a) control data (D~)ex p = 198 ms; (D~)~i m = 50 units; (e4) c = 0.018, (b) quinidine data (D~)~r, = 242 ms; (D~)~i m = 50 units; (e4) q = 0.024.
B.Y. Kogan et al. /Restitution properties and 2D propagation
332
a)
1.1 1.0 0.9~
0.8-
e
t
a
l 8
0.7 0
b)
n
I
I
1
2
control simulation I
DI/D1
3
1.1
1.0
N
................ O
0.9• 0.8-
........o ...... exl~riment with Quinidine 8 simulation
0.70.6 0
I
I
1
2
I
DI/D1
3
Fig. 5. Comparison of restitution curve experimental data [2] to simulate results: (a) normal tissue (control). (b) tissue under quinidine influence. The model parameters that provide the best coincidence with experimental data are: for control RC: G r 311, G t = 0.7, G~ = 1.0, e I = 0.5, e-: - 0.02, e~ - 0.5, e 4 - ke2 = 0.018, lmin = (I.35; For RC with quinidine: Gr = 3(1, G~ =/).7, (;~ - 1,0, E1 : 0 . 5 , E 2 : 0 . 0 2 , 6 3 : 0 . 5 , 6 4 : ke~ = 0.024, [rain- 0.90. Note that beyond the left portion of the graph is the refractory, period. 1.2 1.0 0.8
..........1
0.6 r~
],/I, "
0.4
II .....
'Ii I]
0.2 0.0
1 I Il:: /,
•
I
-0.2 0
25
I
,~D1
50
75
O0
125
J f
150
I DI ~ - - - - - ~ D 2 time
Fig. 6. The transients of l(t) and E(t) in the F4 model. Model parameters are the same as for the control model in fig. 5a. D I / A P D 1 = 0.27; APD2/APD I = 0.71. i n c r e a s e in t h e r e f r a c t o r y p e r i o d in c o m p a r i s o n with AP duration.
4. The influence of restitution properties on wave propagation in 2D excitable m e d i a
T h e t e m p o r a l c h a n g e s o f E ( t ) a n d I ( t ) in fig. 6 c o r r e s p o n d t o t h e p o i n t in fig. 5a w i t h D I / D
1=
T h e s t u d y o f t h e i n f l u e n c e o f h e a r t m u s c l e cell
0.2. C o m p a r i s o n w i t h fig. 3 d i s t i n c t l y s h o w s t h a t
restitution properties on excitation wave propa-
t h e r e s t i t u t i o n p r o p e r t i e s a r e p r e s e n t in t h e m o d i -
gation
fied m o d e l only.
FitzHugh-Nagumo
was
carried
out
using
the
modified
e q u a t i o n o f t h e f o r m o f (1),
B.Y. Kogan et al. / Restitution properties and 2D propagation
(2), and (5). The equations were implemented in a 128 × 128 grid of nodes on the Connection Machine (CM-2) [10]. It is worthwhile to begin with two simple cases: propagation of waves with rectilinear and circular fronts. These two idealized cases facilitate the separation of the effect of wavefront curvature from the effect of the restitution. For both cases we will consider the propagation of two waves generated successively with a specified time interval and the propagation of repeated waves with a specified period. For rectilinear and circular waves, the initial excitation of the tissue produces a wave which propagates with a velocity and a wavelength strictly corresponding to the cell's parameters, its coupling coefficients and the characteristics of the front curvature. If the next wave is generated
333
after a time interval when all inactivating processes are complete, no changes are observed in the propagation of the second wave. However, when the next wave is generated after a time interval when the inactivating processes of the first wave are not yet complete, the wave will propagate through recovering areas of slow outward current ( I ) left by the first wave. Thus, each cell will be excited with nonzero initial conditions of variable ( I ) values which grow with decreasing diastolic interval. This decreases the duration of the generated AP as well as the speed (0) and wavelength (A) of the propagated waves. The corresponding relationships for 0 obtained in the course of computer simulations of rectilinear and concentric wave propagation are shown in figs. 7a and 7b. They cover the situa-
a) 1.2 1.1 1.0
0
0.9-
quinidine
0.8-
control standard
0.7
O0
b)
I
i
200
300
I
I
400
500
600
1.2 .
8
T
1.1 1.0
~
0.9.
~(" /
0.8 O0
----
m
/
, 200
----
n
----
----
e
f
,
-, 300 T
, 400
con,r,
standard , 500
600
Fig. 7. The relationships of speed 0 to period T between two successive excitations of (a) rectilinear waves and (b) concentric wave.
334
B. Y. Kogan et al. / Restitution properties and 2D propagation
tions with slightly, normally and increasingly pronounced restitution properties, corresponding to the simplified FN model with standard parameters, the e~4 FN model fitted to experimental control data, and the ~4 FN model fitted to quinidine data [2]. In the case of circular wave propagation, the 0 values change with wave radius. Therefore, 0 are averaged over the range from the first appearance of a wave to its disappearance. In fig. 7 one can see that the simplified FN model with standard parameters shows no dependence of wave propagation speed on period of stimulation T, whereas the e 4 model reflects a decrease in speed with T as occurs in real tissue. In the case of periodic excitation with a cycle shorter than the period necessary for complete inactivation of outward current (return to its pre-excited condition), it is possible to observe the p h e n o m e n o n of both wavelength and speed of propagation alternation which resembles the well-known rhythm alternation in cell activity [15]. In computer simulations with the modified FN model (with e 4) we observed that the wavelength and speed of propagation are decreased by a decrease of stimulation period T. There exists some critical value of T = Tmin when only odd waves can propagate. This may be a reflection of an actual 2:1 block, and thus may not necessarily be true A P D alternants, since wavelength is a product of A P D and speed of propagation. The value of ~nin varies with restitution properties, e.g. for certain values of T, a simulation of quinidine intoxicated tissue (pronounced restitution properties) exhibits a 2:1 block, while at the same T , a simulation of tissue with normal restitution properties does not. It is impossible to observe any of these p h e n o m e n a in the simplified FN model with standard parameters. We conclude that restitution properties lead to the appearance of nonuniformities in the excitable characteristics (in time and space) of the otherwise uniform tissue when periodic stimulation occurs and when the inactivation processes of the outward current are not complete. This fact gives a new insight into formation and propa-
gation of waves with more complicated forms, in particular, spiral waves (reentry). Gulko and Petrov [17] and independently Winfree [18] were the first to show by computer simulation that spiral waves can arise and propagate in originally homogeneous excitable media (for heart muscle and chemical reactions respectively). The results in ref. [17] were obtained with the physiologically based simplified point model which possessed the restitution properties. Much later, Pertsov and Panfilov [1 1] obtained qualitatively similar results with a simplified FN model which had only negligible restitution properties. Using the model of ref. [17] it was impossible to separately show the effect of restitution properties on spiral wave propagation. Our modification of the FN equations makes it possible to observe the effect of restitution properties on spiral wave generation and propagation. Plates la and Ib show computer simulation results for the models with instant and slow inactivation of generalized outward current (the simplified FN model with standard parameters, and the ~4 FN model respectively). Here we simultaneously display the potential and outward current distributions in space and time. Fig. 8 shows the color designations used to identify levels of AP and outward current. When a premature stimulation is applied at the same site s2 for both models, spiral waves are generated only in the model with slow inactivation of the outward current (plate Ib). For the model with instant inactivation of outward current, premature stimulation will induce spiral waves only when applied on the tail of the previous wave. We found that premature beats can generate single and double spiral waves rotating in different directions with established restitution properties when applied in much larger regions which are distant from the tail of the previous wave. The restitution properties determine the size and location in space (the corresponding times and the levels of the outward current) of the area in the wake of the propagated wave where the application of appropriate
B. Y. Kogan et al. /Restitution properties and 2D propagation
Action Potential
Action Potential
Outward Current
335
Outward Current
sl
tl = 5 units
t2 85
s2
(a)
85
s2 t3 124
124
t4 145
130
t5 165
162
t6 198
190
(b
Plate I. Example of computer simulation results with premature stimulus application. (a) Simplified FitzHugh-Nagumo model with standard parameters (instant inactivation of outward current). Sl is the initial stimulus, $2 the site of premature stimulus.
B.Y. Kogan et al. /Restitution properties and 2D propagation
336
I
E
I
!
O) .
.
.
I
,I f),I
g)
h)
.
l Eth
,l J__
Jd )I I
Imin
I I
Action Potential - E a) b) c) d)
red/yellow - depolarization green - refractory dark green - relative refractory blue - sub-threshold
!
!
outward current - I e) red - activation f) yellow - fast inactivation : I > Imin g) green - slow inactivadon: e4 h) blue - 1= 0
t
Fig. 8. S c h e m a t i c r e p r e s e n t a t i o n of color d e s i g n a t i o n s for action p o t e n t i a l E and c u r r e n t 1 phases.
excitations causes the spiral waves. This area we will term the window of uulnerability (WV) by analogy with that defined by Q u a n and R u d y in ref. [19], where one-dimensional A P p r o p a g a t i o n a p p e a r e d as a result of unidirectional block in ring-shaped excitation tissue. T h e major distinctions from ref. [19] of the W V in our case are: W V is two-dimensional in space, has a g e o m e t r y which is similar to that of the previous wave (because there is no diffusion of the current L), and the stimulation inside the W V leads to the a p p e a r a n c e of double and single spiral wave circulation (see also refs. [12, 13]). It is also useful to define the W V in terms of the value of the decreasing outward current. C o m p u t e r simulation data for the W V are presented in fig. 9. T h e size and location of the W V for the case of circular wave p r o p a g a t i o n are shown in fig. 9a. T h e W V for rectilinear p r o p a g a t e d waves has a rectangular form (fig. 9b). T h e area between the two vertical edges (one which is close and one which is far from the wake of the wave tail) is the space where double spiral waves are g e n e r a t e d in response to the p r o p e r stimulation. T h e area of single spiral wave generation is located near the horizontal edges of the WV. With the decrease of e 4 (from k > 1 to k < 1) the W V displaces in the direction opposite to the previous wave front. For both circular and rectilinear p r o p a g a t e d waves, when the stimulation is applied inside the
WV, the g e n e r a t e d spiral waves have different properties d e p e n d i n g on the distance from the right edge of the W V and on the rate of current I decay (value of k). That is, if the restitution properties are weak (k >> 1) as in ref. [9], the g e n e r a t e d spiral waves quickly b e c o m e stationa r y - t h e i r tip or tips (in the case of double spiral waves) are rotated along a circle of predetermined radius with a constant angular velocity. In the case of normal restitution properties (control), the stimulation applied close to the right edge of W V causes the a p p e a r a n c e of nonstationary spiral wave p r o p a g a t i o n and when applied near the left edge p r o d u c e s a stationary one. W h e n the tissue of the restricted area is u n d e r the influence of quinidine (which leads to an increase of A P duration and of slope of the RC) only the nonstationary spiral waves p r o p a g a t i o n could be observed within the restricted area. D e p e n d i n g on the restitution properties of the tissue, the tip of nonstationary spiral waves could m e a n d e r along curves forming closed or o p e n loop complex curves resembling cycloids. This was shown earlier in ref. [20] using an approximate analytical a p p r o a c h c o m b i n e d with c o m p u t a t i o n s and further studied in ref. [21] by pure c o m p u t e r simulation with a simplified FN point model where the changes of refractory period were tightly connected with those of restitution properties. In ref. [21] the restitution properties were chosen arbitrarily. By contrast, in our c o m p u t e r simulation
337
B.Y. Kogan et aL /Restitution properties and 2D propagation
a) Concentric stimulation •e
b) Rectilin~r stimulation
ee~l
I:::::::::::::::::l • ...... ,
========================
11 l , e
eel
~@O@4HI,4
I,@O@@~1,0 •p'q
~!i!Dliiiii
IP'~" I
1~@41,~@41,el l,e~4,e,e,~l l,@eeee,~
D2
]~3
,~ ~ ,
\
\
[
l~eeeee~ ~e~eeg.t ~eeee,lz,,~
T o,
/
Premature sum
f
be~,l~oe be4,e~ bee~e¢l
Window of vulnerability
Tail of previous wave
x:::=: .i
I
J
,...,
I,@@~l,e@q eeeeee#ele~¢~
i.
t.i:!ii!ii!~?iiii:iil. ::x:=n 1:1::
_,_
li:!:i:i:!:?ii| i
N
,iiiii, Tail of previous wave
Premature stimulus Window dimensions
Concentricstimulus Rectilinear stimulus
Control simulation
D 1 - 20 D 2 - 67 D 3 - 86
D1- 5 D 2 - 67 D3- 89
Quinidine simulation
D1- 20 D 2 - 72 D3- 90
D1- 5 D2- 78 D 3 - 102
Note: The tail (10xl) section of the premature stimulus must lie completely within the window of vulnerability for a spiral wave to occur. The premature stimulus dimensions are 10x20 units. All dimensions correspond to simulation grid dimentions (128 x i28 grid ). Fig. 9. Geometry of the window of vulnerability for rectilinear and concentric stimulation. the restitution properties could be changed without noticeable changes in the refractory period duration and could be adjusted to fit the given experimental data. T h e display of the residual outward current distribution shows the tendency of spiral waves to p r o p a g a t e in the direction in which the instantaneous m i n i m u m of this current is located. In the light of these observations it is necessary to verify the U - s h a p e d e p e n d e n c i e s o f radius of tip wave circulation on the duration of the A P and the dispersion equation previously obtained in the course of simulation [22] and theoretical studies [20, 23].
5. Conclusion T h e restitution properties of excitable media have an important role in the genesis and on the properties of wave propagation. In o r d e r to investigate the influence of restitution properties independently of the refractory period, we p r o p o s e to modify the function e ( E ) in the simplified F N model. A m e t h o d is developed for fitting the model to given experimental restitution data (from heart muscle). These results provide m o r e realistic simulations without substantial complications of the simplified model. T h e p r o p o s e d display of the fast and slow variable ( E and I )
338
B. E Kogan et al. / Restitution properties and 2D propagation
distributions simultaneously is very useful in understanding wave propagation in the presence of restitution properties. The computer simulation results with our modified FN model show at least three qualitative distinction from the simplified FN model with no restitution properties (while all other characteristics and parameters remain unchanged): (1) The development of alternation of speed and wavelength during periodic stimulation at high rates. (2) Non-stationary propagation of spiral waves in media with appropriate restitution properties. (3) Site specific induction of spiral waves with premature stimulation do not fall on the tail of the previous wave. Quantitative effects of differing restitution properties are also shown. In particular, we have shown that quinidine-induced changes in restitution properties cause a shift in the location of the window of vulnerability, and a transition to nonstationary spiral wave propagation. These qualitative and quantitative results for heart muscle require future validation with a point model which reflects physiological properties better than the simplified FN model (for example B e e l e r - R e u t e r ' s [24]). Previous theoretical and simulation studies of excitable media focused on its refractory properties. These studies lead to the dispersion equation and the U-shaped relationship between the radius of the tip of spiral wave circulation and the duration of the action potential. In light of the results presented here, these previous findings should be verified in relation to the real restitution properties.
Acknowledgements T h e authors gratefully acknowledge Dr. William Stevenson of the UCLA Medical Center for his valuable comments to the original text of this paper.
Research in the use of massive parallelism for system simulation in the UCLA Computer Science Department is supported in part by the N A S A / D r y d e n Research Center under Grant NCC 2-374 and NSF Grant BBS 87-14206. Research in experimental cardiac electrophysiology is supported in part by Cedars-Sinai E C H O Research Foundation.
Appendix A. Connection Machine and problem implementation A. 1. Configuration a n d software
A very recent entrant into the high-performance scientific computing field introduced by Thinking Machines Corp. is called the Connection Machine. The full second version of the Connection Machine, the CM-2, with 64K processing elements is packaged into 4096 chips with 16 processing elements on each chip. Within each chip, the processing elements are interconnected to its nearest neighbor in a two-dimensional mesh. The CM-2 provides two types of interprocessor communications. The NEWS (north, east, west, south) grid is a two-dimensional scheme which permits each processor to communicate with its nearest neighbor directly. Processors may also communicate via a 12-dimensional hypercube topology. Each processing element contains a very simple 1-bit arithmetic and logic unit (ALU) and 64K bits of local memory. In addition, each processing element has four 1-bit flag registers, an I / O interface, and share an optional floating point accelerator (among 32 processing elements). The CM-2 may be broken into four different quadrants with a different front-end computer controlling each quadrant. Alternatively, a single front-end computer may control the entire CM-2. Operation on the CM-2 is in single instruction multiple data (SIMD) stream mode. During execution, the front-end computer(s) send a stream of high-level instructions to sequencers which in turn convert these to microcode instructions
B.Y. Kogan et a L / Restitution properties and 2D propagation Table 1 Breakdown of computational task for RK-2 algorithm. numerical integrations calculations of diffusion interactions graphics display
190 s 185 s 15 s
total
390 s
(nanoinstructions ). These
339
In situations where the number of physical processors is not enough to adequately meet the degree of data level parallelism, the CM-2 can support virtual processors.
A.2. Implementation of excitable media
nanoinstructions
are broadcast to the individual processing elements and control their timing and operations. Program development for the CM-2 is supported by three high-level languages: *Lisp, CM Fortran and C*. Software is written in the frontend computer and compiled into the Paris (parallel instruction set) language, which can be issued to the parallel processing units. It is the lowest level protocol by which the front-end computer directs the actions of the CM-2 processors. Paris instructions are sent to a micro-controller which expands them into a series of nanoinstruction, of which there may be from one to several thousand per Paris instruction, and then broadcasts these sequentially to all the processors. Data level parallelism uses a single control sequence or program, where code is executed sequentially. All programs are stored in the frontend computer. Note that the processor memory contains only data and immediate computations, and does not contain code for execution.
All computer simulation described in the p a p e r was performed with a quarter CM-2 configuration with 16K processors. The p a r a m e t e r s used in the simplified F i t z H u g h - N a g u m o equations (standard set) are presented in the main text on p. 338. The resulting p a r a m e t e r s from fitting the new E4 to experimental data from normal and quinidine intoxicated tissue are presented in figs. 4 and 5. In order to take advantage of the data level parallelism in the CM-2, the modified method of line (see ref. [25]) was used. According to this numerical algorithm, the exchange of data between grid nodes can be performed after several steps of integration. The second-order R u n g e - K u t t a formula is used for numerical integration. To reach the required accuracy and stability the step of integration was chosen to be h = 0.025. The exchange of data between nodes was performed after computing three steps of integration. For the chosen e 2 = 0.02 value and for normal cardiac tissue parameters, the conversion for space and time from dimensionless com-
Table 2 Comparison of other simulation experiments. Ref. [22]
Ref. [7]
Ref. [8]
Ref. [10]
HCS-100 (hybrid)
ES 10/40 (IBM)
Cyber 205
CM-2 ( 1 / 4 )
grid dimension
32 × 32
34 x 80
41 x 41 x 41
128 x 128
integration algorithm
continuous (analog)
Euler
Euler
RK-2
0.05
0.03
0.03
type of computer used
integration time step At exchange between grid nodes (time)
0.83
0.05
0.03
0.09 (3 At)
total simulated time (time units)
250
10
100
500
total computer time (s)
600
2400
80
390
340
B. Y. Kogan et al. /Restitution properties anti 2D propagation
puter simulation units into dimensional cardiac tissue values is: 1.0 space unit = 0.1141 cm, 1.0 time unit = 2.0 ms. The space step was chosen as Ax = Ay = 0.85 units. This corresponds for 128 × 128 grids and the excitable media of 4.46 × 4.46 cm 2. The overall simulation run time, 500 time units, requires about 6.5 rain of computing time. This can be broken down among the principle tasks as shown in table 1. Table 2 presents a rough comparison of the computing time requirements of the implementation described in ref. [10] and those reported by other investigators. All three pure digital implementations e m p l o y e d the simplified F i t z H u g h - N a g u m o equations (1), (2), and (3). Only in the hybrid implementation [22] a different, more complicated form of the equation [17] was used.
References [1] R.E. McAllister, D. Noble and R.W. Tsien, Reconstruction of the electrical activity of cardiac purkinje fibers, J. Physiol. 251 11975) 1-59. [2] H.S. Karagueuzian, S.S, Khan, W.J. Mandel and G.A. Diamond, Nonlinear dynamic analysis of temporarily heterogeneous action potential characteristics, PACE (December, 19901. [3] M.A. Allesie, F.1.M. Bonke and F.J.G. Schopman, Circus movement in the rabbit atrial muscle as a mechanism of tachycardia. II. The role of nonuniform recovery, of excitability in the occurrence of unidirectional block, as studied with multiple microelectrodes, Circ. Res. 39 (1976) 168-177. [4] G.K. Moe, W.C. Rheinboldt and J.A. Abildskov, A computer model for atrial fibrillation, Am. Heart J. 67 (1964) 200-220. [5] M. Gerhardt. H. Schuster and J. Tyson. A cellular automation model of excitable media including curvature and dispersion Science, 247 (1990) 1563-1566. [6] A.M. Pertsov, R.N. Chramov and A.V. Panfilov, Sharp increase in refractory period induced by oxidation suppression in F i t z H u g h - N a g u m o model. New mechanism of anti-arrythmic drug action, Biofizika 6 (1981) 1077-1081 [in Russian].
[7] A.M. Pertsov, E.A. Ermakova and A.V. Panfilov, Rotating spiral waves in a modified FitzHugh Nagumo model, Physica D 14 (1984) 117 124. [8] P.I. Nadapurkar and A.T. Winfree, A computation study of twisted linked scroll waves in excitable media, Physica D 29 (1987) 69-83. [9] E.A. Ermakova, A.M. Pertsov and E.E. Shnol, On the interaction of vortices in two-dimensional active media, Physica D 40 (1989) 185-195. [10] B.Y. Kogan, W.J. Karplus and A.T. Pang, Simulation of nonlinear distributed parameter systems on the connection machine, Simulation 15 (1991/) 271 288. [11] A.M. Pertsov and A.V. Panfilov, Spiral waves in active media. Reverberator in FitzHugh-Nagumo model, in: Autowave Processes in Systems with Diffusion (Institute of Applied Physics of the Academy of Science USSR, 1981) pp. 77-89 [in Russian]. [12] A.T. Winfree, Stable particle-like solutions to the nonlinear wave equations of three-dimensional excitable media, SIAM Rev. 32 119911) 1 53. [13] A.T. Winfree, Electrical instability in cardiac muscle: phase singularities and rotors, J. Theor. Biol. 138 (19891 353-405. [14] V.S. Zykov, Simulation of Wave Processes in Excitable Media (Manchester Univ. Press, Manchester 1987) [English translation edited by A.T. Winfree]. [15] E. Carmeliet, Repolarization and frequency in cardiac cells, J. Physiol. (Paris) 73 (1977) 093 923. [16] D. Noble, The Initiation of Heartbeat, 2nd Ed. (Clarendon Press, Oxford, 1979) pp. 82 83. [17] F.B. Gulko and A.A. Petrov, Mechanism of the formation of closed pathways of conduction in excitable media, Biofizika 17 11972) 271-282. [18] A.T. Winfree, Rotating solutions to reaction/diffusion equations, S1AM-AMS Proc. 8 (19741 13-31. [19] W. Q u a n and Y. Rudy, Unidirectional block and reentry of cardiac excitation: A model study, Circulation Res. 66 (1990) 367-382 [20] V.S. Zykov, The kinematics of spiral wave nonstationary circulation in excitable media, Biofizika 32 (1987) 337-34O. [21] E. Lugosi, Analysis of meandering in Zykov kinetics. Physica D 41/( 19891 331-337. [22] B.Y. Kogan, V.S. Zykov and A.A. Petrov, Hybrid computer simulation of stimulative media, in: Simulation of Systems 79 (North-Holland, Amsterdam, 1980) pp. 693 702. [23] J.J. Tyson and J.P. Keener, Singular perturbation theory of traveling waves in excitable media, Physica D 32 (1988) 327 361. [24] G.W. Beeler and tI. Reutcr, Reconstruction of the action potential of ventricular myocardial fibers, J. Physiol. 268 (1977) 177-210. [25] A.A. Petrov, Hybrid method of solving partial ditl'ercntial equations, Avtomat. Telcmekh. 3 (19751. [26] D.R. Chialvo et al., Nature 343 (19901 633-657: Circulation Research 66 (1990) 525-545.
The simplified FitzHugh-Nagumo model with action potential duration restitution: Effects on 2D wave propagation Boris Y. K o g a n , W a i t e r J. Karplus, Brian S. Billett, Alex T. P a n g Computer Science Department, Unieersity of California, Los Angeles, CA, USA
H r a y r S. K a r a g u e u z i a n and S t e v e n S. K h a n D&ision of Cardiology, Cedars-Sinai Medical Center, Los Angeles, CA, USA Received 30 July 1990 Revised manuscript received 10 December 1990 Accepted 10 December 1990 Communicated by A.T. Winfree
A modification of the simplified FitzHugh-Nagumo (FN) equations is proposed for introducing a residual component of the slow variable, which determines the restitution of action potential duration (APD) also known as the interval-excitation duration relationship. The three-step-wise approximation of e(E)which is widely used in current publications is replaced in a new model by a four-step approximation. This change is used for studying by computer simulation the effects of APD restitution properties independently of the APD and refractory period on 2D wave propagation in an isotropic matrix (made by 128 × 128 nodes). The method for fitting the model to the given experimental restitution data (obtained from myocardial cells) is presented. The computer simulations implemented on a massively parallel computer (Connection Machine) showed at least three important qualitative distinctions in behavior which demonstrate the effect of APD restitution: changes in the speed and wavelength of propagated waves with the period of stimulation, non-stationary propagation of spiral waves, and site-specific induction of spiral waves with premature stimulation not on the tail of the previous wave. Quantitative effects of differing restitution properties are expressed in the size and location of a window of vulnerability in 2D excitable media. These windows are characterized by the appearance of single and double spiral waves in response to premature stimulation applied inside the window. Thus the APD restitution incorporated in the FN model produces a significant effect on the formation and propagation of spiral waves.
1. Introduction
It has long been known for excitable media that if a second excitation is initiated soon after the first one, the duration of the second will be considerably shorter than the first. For biological excitable media the relationship of the interval between two excitations to the duration of the second excitation is known as the interval-duration relation, or action potential duration restitut i o n - hereafter referred to simply as restitution. For example, a cardiac muscle cell requires a time equal to several action potential durations
(APDs) to recover fuliy. If stimulation is applied at relatively short cycle lengths, i.e. while the cell is in a state of partial return to its pre-excited condition, then the response will be different from that obtained when the cell has fully recovered. In the case of restitution a typical outcome is that APs initiated at relatively short diastolic intervals (i.e. the interval between successive APs) will result in action potentials of shorter duration. McAllister, Noble and Tsien [1] showed that restitution in cardiac Purkinje fibers is due to the time-dependent component of potassium current which decays very slowly after the completion of
0167-2789/91/$03.50 © 1991 -Elsevier Science Publishers B.V. (North-Holland)
B.K Kogan et al. /Restitution properties and 2D propagation
328
F(E)
'E ~
I
].dl
1
E
F
.9Is Ill" E
Fig. 1. The piece-wise approximation of functions F(E), f(E), and s(K). tan a = G,: t a n f i = (;~: tan y = ( ; j -/:'th is the threshold potential,
the repolarization phase of the action potential. Introducing a residual component of the slow variable into the FitzHugh-Nagumo model simulates this slow decay of the potassium current. It is evident that the time course of restitution of APD directly affects the degree of dispersion of repolarization, when stimulation is applied during the process of restitution [2]#1. It is known that the greater the dispersion of repolarization, the easier will be the induction of spiral waves (circus movement reentry) [3]. It is therefore essential to incorporate these properties in computer models designed to elucidate the properties of impulse propagation. Over the last 30 years, many valuable contributions have been made in wave propagation studies using axiomatic models [4, 5], and the simplified FitzHugh-Nagumo (FN) equations (see refs. [6-9]). However, the parameters of these models, particularly the function s ( E ) in the simplified FN model, provide limited, if at all, restitution properties. Here we propose a method for incorporating restitution properties (residual component of slow variable) into the simplified FN model. Then we present the results of a #1For the role of A P D restitution curve features on the a p p e a r a n c e of cardiac dysrhytbmias and possible suppression by appropriate pharmacological manipulation see ref. [26].
computer simulation of wave propagation in a 2D piece of an isotropic biological homogeneous excitable media (128 × 128 nodes) fitted to experi mental physiological restitution properties of cardiac tissue [2]. These results show the significant qualitative and quantitative distinctions from the simplified FN model without restitution properties. All computer simulation studies were carried out using a massively parallel computer system which is particularly suited to this class of problems, the Connection Machine (CM). For a brief description see appendix A, and for more details see ref. [10].
2. The FitzHugh-Nagumo simplified equations and their properties The basic FN simplified equations when reduced to dimensionless form [6] are aE
a-7- = A E + F ( E )
- I + l~tm,,
aat/ = e(E) I f ( E ) - I],
(1)
(2)
where E is the fast variable (membrane potential displacement between the interior and exterior of
B.Y Kogan et al. / Restitution properties" and 2D propagation
the cell); I the slow variable (generalized outward current), A = a2/0x 2 + ~ 2 / ~ y 2 is the two-dimensional Laplacian operator, F(E) the c u r r e n t voltage characteristic of the fast inward current, f(E) the current-voltage characteristic of the slow outward current and e(E) a small parameter (inversely proportional to the time constant of the slow outward current). The initial and boundary conditions used in the simulation are
E(x,y,O) =l(x,y,O)
=0,
The piece-wise approximation of the functions F ( E ) , f ( E ) , and e(E) is shown in fig. 1. The widely used set of parameters (see refs. [6-9, 11-13]) hereafter referred to as the standard set is: G, = 1,
Gf =
e=e= =0.5
1,
G r = 30,
Eth =
0.16
if 0.00 < E < 0.01,
= e 2=0.01
if0.01
= e 3 = 0.5
if E > 0.95.
(3)
We define an action potential as the time course of the membrane potential which is above the threshold potential Etn. This action potential is
329
initiated by a proper stimulation and formed by the activation and inactivation of inward and outward currents. The relationships between the action potential duration (APD), refractory period (R), and the model parameters G~, Gf, and e are presented in ref. [6] for the model without diffusion (point model). These dependencies reflect qualitatively the essential properties of heart muscle cells, but do not correctly express the restitution properties. The normalized restitution curve for the FN simplified model with standard parameters is compared to normalized data obtained in physiological experiments in fig. 2. The comparison shows the marked disagreement of the model and experimental data. We introduced the normalized form of the restitution curve to facilitate the comparison of computer simulation data obtained in dimensionless form and dimensional data from physiological experiments. Usually, the restitution curve is presented as a function: APD = f(DI), where APD designates the duration of the action potential caused by the second stimulus applied after the specified diastolic interval (DI). The normalized representation of the restitution curve (RC) has the form APD2/APD l = h(DI/APD~),
1.1
1.0
0.9
0.8
*
0.7 0
I
I
1
2
b)
I
DI/DI
3
4
Fig. 2. The normalized restitution curve for (a) results of computer simulation performed by the author using the FitzHugh-Nagumo simplified model with set of parameters [6], and (b) canine heart muscle tissue in normal conditions (control), obtained from experimental data and normalized for comparison [2].
B.Y. Kogan et al. ,/Restitution properties and 2D propagation
330
1.2
r
:
i
i
i E
1.0
i
j
....
i
...
0.8,~
0.6-~
/
0.4 0.2
/ .Y
j
i
/
0.0 i
I 0 •2
o
2op
I aoo
4oo
5oo
t = 0.234; APD2/APD
I
600
dine Fig. 3. T h e E(t) a n d l(t) f o r s h o r t d i a s t o l i c i n t e r v a l . D I / A P D
where A P D l is the duration of the A P obtained after the full completion of the recovery processes. H e r e a f t e r A P D 1 = D 1 and A P D 2 = D 2. T h e transients of E and 1 in a cardiac cycle with a small D I are shown in fig. 3 for the standard p a r a m e t e r s of the F N simplified model. O n e can see that the outward current I reaches zero just after the A P returns to the resting potential. Thus, a stimulus applied after a short D I p r o d u c e s a subsequent A P approximately of the same duration as the previous one, because the time course of the restitution is very fast (almost instantaneous). We therefore modified the characteristics of the transient outward current I ( t ) , by slowing down its time course of inactivation during the interval of negative A P values.
3. S o m e a p p r o a c h e s to restitution curve fitting O n e m e t h o d of modifying the transient outward current I ( t ) is to use two values of e dependent on the sign of d l / d t : e(dI/dt)
=e I
ifdI/dt>O,
=ke I
ifdl/dt<0.
(4)
=
1 (in t h e m o d e l o f ref. [6]).
This a p p r o a c h was p r o p o s e d by Zykov [14] to change the refractory period independently of the A P duration in eqs. (1), (2) with the standard set of parameters. The increase in the refractory period to longer than the duration of the A P contradicts the normal properties of heart muscle tissue. C h a n g i n g the p a r a m e t e r k also changes the restitution properties simultaneously with the refractory period. T h e use of only one p a r a m e t e r ( k ) does not permit the required shaping of l ( t ) to fit restitution properties, nor does it permit the separate investigation of the effect of refractory and restitution properties on impulse propagation. O u r a p p r o a c h is to use a more complicated logical function to control the changes of the small p a r a m e t e r e. Namely e = eI
if E < 0.01 and d l / d t
> 11.
= e2
if E > 0.01 and d l / d t
> O,
= ~3
if 1 > Im~,l and d l / d t
= e 4 = ke 2
otherwise.
This mation The satisfy
<_ O,
(5)
leads to a piece-wise exponential approxiof the current function l ( t ) . values for e l, e 2, and e 3 are chosen to the requirements of the given action po-
B. K Kogan et al. / Restitution properties and 2D propagation
the computer model, it is necessary to choose E4 such that
tential duration and its refractory period, e 4 is chosen to meet the required restitution properties. The proper values of e 4 and Imin can be matched to the experimental restitution curve by computer simulation trial and error or by a direct analytical approach. The normal restitution curve is commonly approximated by the formula [15]
D = Do( 1 - e-t~v).
,E"4 = / 3 ( O ~ ) e x p ( Doc)sim .
(6)
D/D~) =/3T.
(8)
This formula is based on the well known physiological fact that the time constant of the restitution curve (/3) is similar to the time constant of the variable x 1 which determines the potassium current ixl at negative potentials [16]. The residual component of the slow variable in the FN model approximates this potassium current. The scale factor rn (Ooc)exp/(O~)sim is used to transform the dimensional /3 to the dimensionless e 4. The results of restitution curve (RC) fitting for heart muscle tissues under normal conditions (control) and with the tissue exposed to quinidine are presented in figs. 5a and 5b, respectively. In both cases the close agreement with experimental data was obtained by making E4 =/3m and choosing an appropriate value for Imin' The maximum and rms error in the region 0.5 < D I / A P D I < 3.0 for the control data are respectively: 0.03 and 0.0066, and for the quinidine are: 0.038 and 0.010. It was found that a small decrease of the parameter G f improved the fitting of the quinidine experimental data. This approach also gives a slight
Here D = A P D 2 and D~ = the plateau A P D after an infinitely long DI; /3 is a constant, and T = APD~ + DI. From (6) follows
-ln(1 -
331
=
(7)
Thus the relationship between - I n ( 1 - D / D ~ ) and T is linear, and the experimental data can be represented by a straight line, which can be obtained by the least-squares method of approximation. The results of this approximation are shown in fig. 4 for two sets of experimental data [2] from normal tissue and tissue intoxicated by quinidine. The slope of the line obtained in this manner directly shows the values of/3. To reproduce the restitution properties of the experimental data in
J
_J
J
..l y
J
,...--
J
J
J
2
.J
f
1
2O0
300 y = 0.0047x y =
0.0055x
400 + 0,14 -
1.08
500
R^2 = 0.973
D-----
R^2 = 0.978
"=
600 700 a) control data b) quinidine data
T
800
Fig. 4. Determination of e 4 from experimental data: (a) control data (D~)ex p = 198 ms; (D~)~i m = 50 units; (e4) c = 0.018, (b) quinidine data (D~)~r, = 242 ms; (D~)~i m = 50 units; (e4) q = 0.024.
B.Y. Kogan et al. /Restitution properties and 2D propagation
332
a)
1.1 1.0 0.9~
0.8-
e
t
a
l 8
0.7 0
b)
n
I
I
1
2
control simulation I
DI/D1
3
1.1
1.0
N
................ O
0.9• 0.8-
........o ...... exl~riment with Quinidine 8 simulation
0.70.6 0
I
I
1
2
I
DI/D1
3
Fig. 5. Comparison of restitution curve experimental data [2] to simulate results: (a) normal tissue (control). (b) tissue under quinidine influence. The model parameters that provide the best coincidence with experimental data are: for control RC: G r 311, G t = 0.7, G~ = 1.0, e I = 0.5, e-: - 0.02, e~ - 0.5, e 4 - ke2 = 0.018, lmin = (I.35; For RC with quinidine: Gr = 3(1, G~ =/).7, (;~ - 1,0, E1 : 0 . 5 , E 2 : 0 . 0 2 , 6 3 : 0 . 5 , 6 4 : ke~ = 0.024, [rain- 0.90. Note that beyond the left portion of the graph is the refractory, period. 1.2 1.0 0.8
..........1
0.6 r~
],/I, "
0.4
II .....
'Ii I]
0.2 0.0
1 I Il:: /,
•
I
-0.2 0
25
I
,~D1
50
75
O0
125
J f
150
I DI ~ - - - - - ~ D 2 time
Fig. 6. The transients of l(t) and E(t) in the F4 model. Model parameters are the same as for the control model in fig. 5a. D I / A P D 1 = 0.27; APD2/APD I = 0.71. i n c r e a s e in t h e r e f r a c t o r y p e r i o d in c o m p a r i s o n with AP duration.
4. The influence of restitution properties on wave propagation in 2D excitable m e d i a
T h e t e m p o r a l c h a n g e s o f E ( t ) a n d I ( t ) in fig. 6 c o r r e s p o n d t o t h e p o i n t in fig. 5a w i t h D I / D
1=
T h e s t u d y o f t h e i n f l u e n c e o f h e a r t m u s c l e cell
0.2. C o m p a r i s o n w i t h fig. 3 d i s t i n c t l y s h o w s t h a t
restitution properties on excitation wave propa-
t h e r e s t i t u t i o n p r o p e r t i e s a r e p r e s e n t in t h e m o d i -
gation
fied m o d e l only.
FitzHugh-Nagumo
was
carried
out
using
the
modified
e q u a t i o n o f t h e f o r m o f (1),
B.Y. Kogan et al. / Restitution properties and 2D propagation
(2), and (5). The equations were implemented in a 128 × 128 grid of nodes on the Connection Machine (CM-2) [10]. It is worthwhile to begin with two simple cases: propagation of waves with rectilinear and circular fronts. These two idealized cases facilitate the separation of the effect of wavefront curvature from the effect of the restitution. For both cases we will consider the propagation of two waves generated successively with a specified time interval and the propagation of repeated waves with a specified period. For rectilinear and circular waves, the initial excitation of the tissue produces a wave which propagates with a velocity and a wavelength strictly corresponding to the cell's parameters, its coupling coefficients and the characteristics of the front curvature. If the next wave is generated
333
after a time interval when all inactivating processes are complete, no changes are observed in the propagation of the second wave. However, when the next wave is generated after a time interval when the inactivating processes of the first wave are not yet complete, the wave will propagate through recovering areas of slow outward current ( I ) left by the first wave. Thus, each cell will be excited with nonzero initial conditions of variable ( I ) values which grow with decreasing diastolic interval. This decreases the duration of the generated AP as well as the speed (0) and wavelength (A) of the propagated waves. The corresponding relationships for 0 obtained in the course of computer simulations of rectilinear and concentric wave propagation are shown in figs. 7a and 7b. They cover the situa-
a) 1.2 1.1 1.0
0
0.9-
quinidine
0.8-
control standard
0.7
O0
b)
I
i
200
300
I
I
400
500
600
1.2 .
8
T
1.1 1.0
~
0.9.
~(" /
0.8 O0
----
m
/
, 200
----
n
----
----
e
f
,
-, 300 T
, 400
con,r,
standard , 500
600
Fig. 7. The relationships of speed 0 to period T between two successive excitations of (a) rectilinear waves and (b) concentric wave.
334
B. Y. Kogan et al. / Restitution properties and 2D propagation
tions with slightly, normally and increasingly pronounced restitution properties, corresponding to the simplified FN model with standard parameters, the e~4 FN model fitted to experimental control data, and the ~4 FN model fitted to quinidine data [2]. In the case of circular wave propagation, the 0 values change with wave radius. Therefore, 0 are averaged over the range from the first appearance of a wave to its disappearance. In fig. 7 one can see that the simplified FN model with standard parameters shows no dependence of wave propagation speed on period of stimulation T, whereas the e 4 model reflects a decrease in speed with T as occurs in real tissue. In the case of periodic excitation with a cycle shorter than the period necessary for complete inactivation of outward current (return to its pre-excited condition), it is possible to observe the p h e n o m e n o n of both wavelength and speed of propagation alternation which resembles the well-known rhythm alternation in cell activity [15]. In computer simulations with the modified FN model (with e 4) we observed that the wavelength and speed of propagation are decreased by a decrease of stimulation period T. There exists some critical value of T = Tmin when only odd waves can propagate. This may be a reflection of an actual 2:1 block, and thus may not necessarily be true A P D alternants, since wavelength is a product of A P D and speed of propagation. The value of ~nin varies with restitution properties, e.g. for certain values of T, a simulation of quinidine intoxicated tissue (pronounced restitution properties) exhibits a 2:1 block, while at the same T , a simulation of tissue with normal restitution properties does not. It is impossible to observe any of these p h e n o m e n a in the simplified FN model with standard parameters. We conclude that restitution properties lead to the appearance of nonuniformities in the excitable characteristics (in time and space) of the otherwise uniform tissue when periodic stimulation occurs and when the inactivation processes of the outward current are not complete. This fact gives a new insight into formation and propa-
gation of waves with more complicated forms, in particular, spiral waves (reentry). Gulko and Petrov [17] and independently Winfree [18] were the first to show by computer simulation that spiral waves can arise and propagate in originally homogeneous excitable media (for heart muscle and chemical reactions respectively). The results in ref. [17] were obtained with the physiologically based simplified point model which possessed the restitution properties. Much later, Pertsov and Panfilov [1 1] obtained qualitatively similar results with a simplified FN model which had only negligible restitution properties. Using the model of ref. [17] it was impossible to separately show the effect of restitution properties on spiral wave propagation. Our modification of the FN equations makes it possible to observe the effect of restitution properties on spiral wave generation and propagation. Plates la and Ib show computer simulation results for the models with instant and slow inactivation of generalized outward current (the simplified FN model with standard parameters, and the ~4 FN model respectively). Here we simultaneously display the potential and outward current distributions in space and time. Fig. 8 shows the color designations used to identify levels of AP and outward current. When a premature stimulation is applied at the same site s2 for both models, spiral waves are generated only in the model with slow inactivation of the outward current (plate Ib). For the model with instant inactivation of outward current, premature stimulation will induce spiral waves only when applied on the tail of the previous wave. We found that premature beats can generate single and double spiral waves rotating in different directions with established restitution properties when applied in much larger regions which are distant from the tail of the previous wave. The restitution properties determine the size and location in space (the corresponding times and the levels of the outward current) of the area in the wake of the propagated wave where the application of appropriate
B. Y. Kogan et al. /Restitution properties and 2D propagation
Action Potential
Action Potential
Outward Current
335
Outward Current
sl
tl = 5 units
t2 85
s2
(a)
85
s2 t3 124
124
t4 145
130
t5 165
162
t6 198
190
(b
Plate I. Example of computer simulation results with premature stimulus application. (a) Simplified FitzHugh-Nagumo model with standard parameters (instant inactivation of outward current). Sl is the initial stimulus, $2 the site of premature stimulus.
B.Y. Kogan et al. /Restitution properties and 2D propagation
336
I
E
I
!
O) .
.
.
I
,I f),I
g)
h)
.
l Eth
,l J__
Jd )I I
Imin
I I
Action Potential - E a) b) c) d)
red/yellow - depolarization green - refractory dark green - relative refractory blue - sub-threshold
!
!
outward current - I e) red - activation f) yellow - fast inactivation : I > Imin g) green - slow inactivadon: e4 h) blue - 1= 0
t
Fig. 8. S c h e m a t i c r e p r e s e n t a t i o n of color d e s i g n a t i o n s for action p o t e n t i a l E and c u r r e n t 1 phases.
excitations causes the spiral waves. This area we will term the window of uulnerability (WV) by analogy with that defined by Q u a n and R u d y in ref. [19], where one-dimensional A P p r o p a g a t i o n a p p e a r e d as a result of unidirectional block in ring-shaped excitation tissue. T h e major distinctions from ref. [19] of the W V in our case are: W V is two-dimensional in space, has a g e o m e t r y which is similar to that of the previous wave (because there is no diffusion of the current L), and the stimulation inside the W V leads to the a p p e a r a n c e of double and single spiral wave circulation (see also refs. [12, 13]). It is also useful to define the W V in terms of the value of the decreasing outward current. C o m p u t e r simulation data for the W V are presented in fig. 9. T h e size and location of the W V for the case of circular wave p r o p a g a t i o n are shown in fig. 9a. T h e W V for rectilinear p r o p a g a t e d waves has a rectangular form (fig. 9b). T h e area between the two vertical edges (one which is close and one which is far from the wake of the wave tail) is the space where double spiral waves are g e n e r a t e d in response to the p r o p e r stimulation. T h e area of single spiral wave generation is located near the horizontal edges of the WV. With the decrease of e 4 (from k > 1 to k < 1) the W V displaces in the direction opposite to the previous wave front. For both circular and rectilinear p r o p a g a t e d waves, when the stimulation is applied inside the
WV, the g e n e r a t e d spiral waves have different properties d e p e n d i n g on the distance from the right edge of the W V and on the rate of current I decay (value of k). That is, if the restitution properties are weak (k >> 1) as in ref. [9], the g e n e r a t e d spiral waves quickly b e c o m e stationa r y - t h e i r tip or tips (in the case of double spiral waves) are rotated along a circle of predetermined radius with a constant angular velocity. In the case of normal restitution properties (control), the stimulation applied close to the right edge of W V causes the a p p e a r a n c e of nonstationary spiral wave p r o p a g a t i o n and when applied near the left edge p r o d u c e s a stationary one. W h e n the tissue of the restricted area is u n d e r the influence of quinidine (which leads to an increase of A P duration and of slope of the RC) only the nonstationary spiral waves p r o p a g a t i o n could be observed within the restricted area. D e p e n d i n g on the restitution properties of the tissue, the tip of nonstationary spiral waves could m e a n d e r along curves forming closed or o p e n loop complex curves resembling cycloids. This was shown earlier in ref. [20] using an approximate analytical a p p r o a c h c o m b i n e d with c o m p u t a t i o n s and further studied in ref. [21] by pure c o m p u t e r simulation with a simplified FN point model where the changes of refractory period were tightly connected with those of restitution properties. In ref. [21] the restitution properties were chosen arbitrarily. By contrast, in our c o m p u t e r simulation
337
B.Y. Kogan et aL /Restitution properties and 2D propagation
a) Concentric stimulation •e
b) Rectilin~r stimulation
ee~l
I:::::::::::::::::l • ...... ,
========================
11 l , e
eel
~@O@4HI,4
I,@O@@~1,0 •p'q
~!i!Dliiiii
IP'~" I
1~@41,~@41,el l,e~4,e,e,~l l,@eeee,~
D2
]~3
,~ ~ ,
\
\
[
l~eeeee~ ~e~eeg.t ~eeee,lz,,~
T o,
/
Premature sum
f
be~,l~oe be4,e~ bee~e¢l
Window of vulnerability
Tail of previous wave
x:::=: .i
I
J
,...,
I,@@~l,e@q eeeeee#ele~¢~
i.
t.i:!ii!ii!~?iiii:iil. ::x:=n 1:1::
_,_
li:!:i:i:!:?ii| i
N
,iiiii, Tail of previous wave
Premature stimulus Window dimensions
Concentricstimulus Rectilinear stimulus
Control simulation
D 1 - 20 D 2 - 67 D 3 - 86
D1- 5 D 2 - 67 D3- 89
Quinidine simulation
D1- 20 D 2 - 72 D3- 90
D1- 5 D2- 78 D 3 - 102
Note: The tail (10xl) section of the premature stimulus must lie completely within the window of vulnerability for a spiral wave to occur. The premature stimulus dimensions are 10x20 units. All dimensions correspond to simulation grid dimentions (128 x i28 grid ). Fig. 9. Geometry of the window of vulnerability for rectilinear and concentric stimulation. the restitution properties could be changed without noticeable changes in the refractory period duration and could be adjusted to fit the given experimental data. T h e display of the residual outward current distribution shows the tendency of spiral waves to p r o p a g a t e in the direction in which the instantaneous m i n i m u m of this current is located. In the light of these observations it is necessary to verify the U - s h a p e d e p e n d e n c i e s o f radius of tip wave circulation on the duration of the A P and the dispersion equation previously obtained in the course of simulation [22] and theoretical studies [20, 23].
5. Conclusion T h e restitution properties of excitable media have an important role in the genesis and on the properties of wave propagation. In o r d e r to investigate the influence of restitution properties independently of the refractory period, we p r o p o s e to modify the function e ( E ) in the simplified F N model. A m e t h o d is developed for fitting the model to given experimental restitution data (from heart muscle). These results provide m o r e realistic simulations without substantial complications of the simplified model. T h e p r o p o s e d display of the fast and slow variable ( E and I )
338
B. E Kogan et al. / Restitution properties and 2D propagation
distributions simultaneously is very useful in understanding wave propagation in the presence of restitution properties. The computer simulation results with our modified FN model show at least three qualitative distinction from the simplified FN model with no restitution properties (while all other characteristics and parameters remain unchanged): (1) The development of alternation of speed and wavelength during periodic stimulation at high rates. (2) Non-stationary propagation of spiral waves in media with appropriate restitution properties. (3) Site specific induction of spiral waves with premature stimulation do not fall on the tail of the previous wave. Quantitative effects of differing restitution properties are also shown. In particular, we have shown that quinidine-induced changes in restitution properties cause a shift in the location of the window of vulnerability, and a transition to nonstationary spiral wave propagation. These qualitative and quantitative results for heart muscle require future validation with a point model which reflects physiological properties better than the simplified FN model (for example B e e l e r - R e u t e r ' s [24]). Previous theoretical and simulation studies of excitable media focused on its refractory properties. These studies lead to the dispersion equation and the U-shaped relationship between the radius of the tip of spiral wave circulation and the duration of the action potential. In light of the results presented here, these previous findings should be verified in relation to the real restitution properties.
Acknowledgements T h e authors gratefully acknowledge Dr. William Stevenson of the UCLA Medical Center for his valuable comments to the original text of this paper.
Research in the use of massive parallelism for system simulation in the UCLA Computer Science Department is supported in part by the N A S A / D r y d e n Research Center under Grant NCC 2-374 and NSF Grant BBS 87-14206. Research in experimental cardiac electrophysiology is supported in part by Cedars-Sinai E C H O Research Foundation.
Appendix A. Connection Machine and problem implementation A. 1. Configuration a n d software
A very recent entrant into the high-performance scientific computing field introduced by Thinking Machines Corp. is called the Connection Machine. The full second version of the Connection Machine, the CM-2, with 64K processing elements is packaged into 4096 chips with 16 processing elements on each chip. Within each chip, the processing elements are interconnected to its nearest neighbor in a two-dimensional mesh. The CM-2 provides two types of interprocessor communications. The NEWS (north, east, west, south) grid is a two-dimensional scheme which permits each processor to communicate with its nearest neighbor directly. Processors may also communicate via a 12-dimensional hypercube topology. Each processing element contains a very simple 1-bit arithmetic and logic unit (ALU) and 64K bits of local memory. In addition, each processing element has four 1-bit flag registers, an I / O interface, and share an optional floating point accelerator (among 32 processing elements). The CM-2 may be broken into four different quadrants with a different front-end computer controlling each quadrant. Alternatively, a single front-end computer may control the entire CM-2. Operation on the CM-2 is in single instruction multiple data (SIMD) stream mode. During execution, the front-end computer(s) send a stream of high-level instructions to sequencers which in turn convert these to microcode instructions
B.Y. Kogan et a L / Restitution properties and 2D propagation Table 1 Breakdown of computational task for RK-2 algorithm. numerical integrations calculations of diffusion interactions graphics display
190 s 185 s 15 s
total
390 s
(nanoinstructions ). These
339
In situations where the number of physical processors is not enough to adequately meet the degree of data level parallelism, the CM-2 can support virtual processors.
A.2. Implementation of excitable media
nanoinstructions
are broadcast to the individual processing elements and control their timing and operations. Program development for the CM-2 is supported by three high-level languages: *Lisp, CM Fortran and C*. Software is written in the frontend computer and compiled into the Paris (parallel instruction set) language, which can be issued to the parallel processing units. It is the lowest level protocol by which the front-end computer directs the actions of the CM-2 processors. Paris instructions are sent to a micro-controller which expands them into a series of nanoinstruction, of which there may be from one to several thousand per Paris instruction, and then broadcasts these sequentially to all the processors. Data level parallelism uses a single control sequence or program, where code is executed sequentially. All programs are stored in the frontend computer. Note that the processor memory contains only data and immediate computations, and does not contain code for execution.
All computer simulation described in the p a p e r was performed with a quarter CM-2 configuration with 16K processors. The p a r a m e t e r s used in the simplified F i t z H u g h - N a g u m o equations (standard set) are presented in the main text on p. 338. The resulting p a r a m e t e r s from fitting the new E4 to experimental data from normal and quinidine intoxicated tissue are presented in figs. 4 and 5. In order to take advantage of the data level parallelism in the CM-2, the modified method of line (see ref. [25]) was used. According to this numerical algorithm, the exchange of data between grid nodes can be performed after several steps of integration. The second-order R u n g e - K u t t a formula is used for numerical integration. To reach the required accuracy and stability the step of integration was chosen to be h = 0.025. The exchange of data between nodes was performed after computing three steps of integration. For the chosen e 2 = 0.02 value and for normal cardiac tissue parameters, the conversion for space and time from dimensionless com-
Table 2 Comparison of other simulation experiments. Ref. [22]
Ref. [7]
Ref. [8]
Ref. [10]
HCS-100 (hybrid)
ES 10/40 (IBM)
Cyber 205
CM-2 ( 1 / 4 )
grid dimension
32 × 32
34 x 80
41 x 41 x 41
128 x 128
integration algorithm
continuous (analog)
Euler
Euler
RK-2
0.05
0.03
0.03
type of computer used
integration time step At exchange between grid nodes (time)
0.83
0.05
0.03
0.09 (3 At)
total simulated time (time units)
250
10
100
500
total computer time (s)
600
2400
80
390
340
B. Y. Kogan et al. /Restitution properties anti 2D propagation
puter simulation units into dimensional cardiac tissue values is: 1.0 space unit = 0.1141 cm, 1.0 time unit = 2.0 ms. The space step was chosen as Ax = Ay = 0.85 units. This corresponds for 128 × 128 grids and the excitable media of 4.46 × 4.46 cm 2. The overall simulation run time, 500 time units, requires about 6.5 rain of computing time. This can be broken down among the principle tasks as shown in table 1. Table 2 presents a rough comparison of the computing time requirements of the implementation described in ref. [10] and those reported by other investigators. All three pure digital implementations e m p l o y e d the simplified F i t z H u g h - N a g u m o equations (1), (2), and (3). Only in the hybrid implementation [22] a different, more complicated form of the equation [17] was used.
References [1] R.E. McAllister, D. Noble and R.W. Tsien, Reconstruction of the electrical activity of cardiac purkinje fibers, J. Physiol. 251 11975) 1-59. [2] H.S. Karagueuzian, S.S, Khan, W.J. Mandel and G.A. Diamond, Nonlinear dynamic analysis of temporarily heterogeneous action potential characteristics, PACE (December, 19901. [3] M.A. Allesie, F.1.M. Bonke and F.J.G. Schopman, Circus movement in the rabbit atrial muscle as a mechanism of tachycardia. II. The role of nonuniform recovery, of excitability in the occurrence of unidirectional block, as studied with multiple microelectrodes, Circ. Res. 39 (1976) 168-177. [4] G.K. Moe, W.C. Rheinboldt and J.A. Abildskov, A computer model for atrial fibrillation, Am. Heart J. 67 (1964) 200-220. [5] M. Gerhardt. H. Schuster and J. Tyson. A cellular automation model of excitable media including curvature and dispersion Science, 247 (1990) 1563-1566. [6] A.M. Pertsov, R.N. Chramov and A.V. Panfilov, Sharp increase in refractory period induced by oxidation suppression in F i t z H u g h - N a g u m o model. New mechanism of anti-arrythmic drug action, Biofizika 6 (1981) 1077-1081 [in Russian].
[7] A.M. Pertsov, E.A. Ermakova and A.V. Panfilov, Rotating spiral waves in a modified FitzHugh Nagumo model, Physica D 14 (1984) 117 124. [8] P.I. Nadapurkar and A.T. Winfree, A computation study of twisted linked scroll waves in excitable media, Physica D 29 (1987) 69-83. [9] E.A. Ermakova, A.M. Pertsov and E.E. Shnol, On the interaction of vortices in two-dimensional active media, Physica D 40 (1989) 185-195. [10] B.Y. Kogan, W.J. Karplus and A.T. Pang, Simulation of nonlinear distributed parameter systems on the connection machine, Simulation 15 (1991/) 271 288. [11] A.M. Pertsov and A.V. Panfilov, Spiral waves in active media. Reverberator in FitzHugh-Nagumo model, in: Autowave Processes in Systems with Diffusion (Institute of Applied Physics of the Academy of Science USSR, 1981) pp. 77-89 [in Russian]. [12] A.T. Winfree, Stable particle-like solutions to the nonlinear wave equations of three-dimensional excitable media, SIAM Rev. 32 119911) 1 53. [13] A.T. Winfree, Electrical instability in cardiac muscle: phase singularities and rotors, J. Theor. Biol. 138 (19891 353-405. [14] V.S. Zykov, Simulation of Wave Processes in Excitable Media (Manchester Univ. Press, Manchester 1987) [English translation edited by A.T. Winfree]. [15] E. Carmeliet, Repolarization and frequency in cardiac cells, J. Physiol. (Paris) 73 (1977) 093 923. [16] D. Noble, The Initiation of Heartbeat, 2nd Ed. (Clarendon Press, Oxford, 1979) pp. 82 83. [17] F.B. Gulko and A.A. Petrov, Mechanism of the formation of closed pathways of conduction in excitable media, Biofizika 17 11972) 271-282. [18] A.T. Winfree, Rotating solutions to reaction/diffusion equations, S1AM-AMS Proc. 8 (19741 13-31. [19] W. Q u a n and Y. Rudy, Unidirectional block and reentry of cardiac excitation: A model study, Circulation Res. 66 (1990) 367-382 [20] V.S. Zykov, The kinematics of spiral wave nonstationary circulation in excitable media, Biofizika 32 (1987) 337-34O. [21] E. Lugosi, Analysis of meandering in Zykov kinetics. Physica D 41/( 19891 331-337. [22] B.Y. Kogan, V.S. Zykov and A.A. Petrov, Hybrid computer simulation of stimulative media, in: Simulation of Systems 79 (North-Holland, Amsterdam, 1980) pp. 693 702. [23] J.J. Tyson and J.P. Keener, Singular perturbation theory of traveling waves in excitable media, Physica D 32 (1988) 327 361. [24] G.W. Beeler and tI. Reutcr, Reconstruction of the action potential of ventricular myocardial fibers, J. Physiol. 268 (1977) 177-210. [25] A.A. Petrov, Hybrid method of solving partial ditl'ercntial equations, Avtomat. Telcmekh. 3 (19751. [26] D.R. Chialvo et al., Nature 343 (19901 633-657: Circulation Research 66 (1990) 525-545.