## Abstract

Theoretical and experimental studies have shown that restitution of the cardiac action potential (AP) duration (APD) plays a major role in predisposing ventricular tachycardia to degenerate to ventricular fibrillation, whereas its role in atrial fibrillation is unclear. We used the Courtemanche human atrial cell model and the Luo-Rudy guinea pig ventricular model to compare the roles of electrical restitution in destabilizing spiral wave reentry in simulated two-dimensional homogeneous atrial and ventricular tissue. Because atrial AP morphology is complex, we also validated the usefulness of effective refractory period (ERP) restitution. ERP restitution correlated best with APD restitution at transmembrane potentials greater than or equal to −62 mV, and its steepness was a reliable predictor of spiral wave phenotype (stable, meandering, hypermeandering, and breakup) in both atrial and ventricular tissue. Spiral breakup or single hypermeandering spirals occurred when the slope of ERP restitution exceeded 1 at short diastolic intervals. Thus ERP restitution, which is easier to measure clinically than APD restitution, is a reliable determinant of spiral wave stability in simulated atrial and ventricular tissue.

- spiral wave breakup
- atrial modeling

ventricular fibrillation (VF) accounts for over 220,000 sudden cardiac deaths per year in the United States alone, and atrial fibrillation (AF) is a common arrhythmia afflicting as many as 5% of Americans over the age of 65 (19) and accounting for one-third of strokes in the elderly. Whereas VF episodes in high-risk patients can be aborted by implantable cardiodefibrillators, the ability to maintain sinus rhythm in patients with AF using currently available drugs is limited, such that clinical therapy is often confined to ventricular rate control and anticoagulation to prevent thromboemboli or catheter ablation to attempt cure in a select group of patients (16). Recent experimental and modeling studies in ventricular muscle suggest that dynamically induced instabilities in cardiac wave fronts, in addition to preexisting heterogeneities, can contribute to the wave break that fuels fibrillation (11, 15, 27-29, 32, 35). Electrical restitution of the action potential (AP) duration (APD) and conduction velocity (CV) in ventricular muscle has been shown to be a key factor regulating this dynamical instability. In ventricular cell models, spiral reentry quickly breaks up into a complex multiple spiral VF-like state if the slope of APD restitution curve is steep (>1) over a sufficiently wide range of diastolic intervals (DI). If this range of DIs is narrowed, the spiral wave remains intact but still hypermeanders chaotically. Only when the slope of APD restitution is <1 everywhere does the spiral wave stabilize to quasiperiodic meander or stationary behavior (11, 27, 28). This type of wave break arises solely from the dynamics of cardiac propagation and is related to steep APD restitution causing oscillations in wavelength before localized wave break (28), without any requirement for fixed heterogeneities in the tissue.

In atrial tissue, the effects of electrical restitution on the dynamical instability of atrial wave fronts have not been characterized and should not be assumed to essential for wave break-perpetuating fibrillation for several reasons. First, the atrium has a complex geometry, including the atrial appendages, the pectinate muscle network, fossa ovalis, and specialized tissues such as the crista terminalis and Bachmann's bundle as well as multiple orifices for veins, arteries, and valves (14). Recent modeling using simplified approximations to structures such as the pectinate muscles and crista terminalis have shown that these anatomic structures may support AF initiation and maintenance (8, 12, 34). Second, atrial tissue displays substantial regional electrophysiological heterogeneity in AP morphology, degree of anisotropy (such as the crista terminalis), and CV (with slow conduction at the inferior vena cava/tricuspid isthmus and perisinus nodal and atrioventricular nodal regions). The atrial AP has a complex morphology and, depending on heart rate and other factors, can show a spike-and-dome appearance, a triangular shape with no sustained plateau, or ventricular-like AP plateau (10, 30). These features may create substantial dispersion of refractoriness (10), which, coupled with slow conduction, provides a substrate for wave break and reentry. Over four decades ago, modeling studies (21) showed how such preexisting electrophysiological heterogeneities such as dispersion of refractoriness could promote AF. Finally, the electrophysiological and structural properties of the atrium become significantly remodeled by the process of AF itself (1, 4, 6, 7, 13, 22, 33), including increased effective refractory period (ERP) heterogeneity (9), so that AF is likely to become chronic if it persists for more than several weeks.

The purpose of this simulation study was to investigate the role of electrical restitution on wave front stability, emphasizing simulated atrial tissue using a realistic physiological human atrial cell model. Two detailed mathematical models of the human atrial AP, both incorporating intracellular Ca^{2+} dynamics as well as ionic currents, have been published, by Courtemanche et al. (3) and by Nygren et al. (23). Although based on comparably rigorous experimental data, the two models have different AP shapes and electrophysiological properties, and thus may be relevant to different regions of the human atrium. The detailed comparison of ionic properties in both the Nygren model and Courtemanche model can be found in Ref. 24. In simulated two-dimensional (2-D) homogeneous tissue, a spiral wave initiated using the Nygren model remains intact, whereas in the Courtemanche model the spiral wave breaks up (24). Because one focus of this study was to examine whether modifying atrial electrical restitution can increase wave stability, we choose the Courtemanche model as our single cell model. Also, because the complex morphology of the atrial AP makes APD measurement more equivocal than in the ventricle, we examined whether ERP restitution could be used as a substitute for APD restitution in atrial as well as ventricular tissue.

### Glossary

- 2-D
- Two-dimensional
- 3-D
- Three-dimensional
- Δ
*t*_{max} - Maximum time step
- Δ
*t*_{min} - Minimum time step
- τ
_{y} - Time constant of a given gating variable
- AF
- Atrial fibrillation
- AP
- Action potential
- APD
- Action potential duration
- APD
_{50} - Time from onset of depolarization until the cell repolarizes to 50% repolarization
- APD
_{90} - Time from onset of depolarization until the cell repolarizes to 90% repolarization
- APD
_{−32 mV} - Portion of the cardiac cycle during which
*V*≥ −32 mV - APD
_{−62 mV} - Portion of the cariac cycle during which
*V*≥ −62 mV *C*_{m}- Total capacitance
- CL
- Cycle length
- CV
- Conduction velocity
- [Ca
^{2+}]_{i} - Intracellular Ca
^{2+}concentration *D*- Isotropic diffusion coefficient
- DI
- Diastolic interval
*E*- Equilibrium potential of a given ionic channel
- ECG
- Electrocardiogram
- ERP
- Effective refractory period
*f*- Gating variable
*f*_{Ca}- Gating variable
*F*- Faraday's constant
*g*- Conductance of a given ionic channel
*G*_{si}- Maximum conductance of L-type Ca
^{2+}channel *h*- Gating variable
*I*_{Ca,b}- Ca
^{2+}background current *I*_{Ca,L}- Slow inward L-type Ca
^{2+}current *I*_{ion}- Total ionic current
*I*_{K1}- Time-independent K
^{+}current *I*_{Kr}- Rapid delayed outward K
^{+}current *I*_{Ks}- Slow delayed outward K
^{+}current *I*_{Kur}- Ultrarapid delayed outward K
^{+}current *I*_{Na}- Fast inward Na
^{+}current *I*_{Na,b}- Na
^{+}background current *I*_{NaCa}- Na
^{+}/Ca^{2+}exchanger current *I*_{NaK}- Na
^{+}-K^{+}pump current *I*_{pCa}- Sarcolemmal Ca
^{2+}pump current *I*_{stim}- Stimulus current
*I*_{to}- Transient outward K
^{+}current *j*- Gating variable
- JSR
- Junctional sarcoplasmic reticulum
- [K
^{+}]_{i} - Intracellular K
^{+}concentration - LR1
- Phase 1 of the Luo-Rudy model
*m*- Gating variable
*n*- Number
- NSR
- Nonjunctional sarcoplasmic reticulum
- [Na
^{+}]_{i} - Intracellular Na
^{+}concentration *O*_{a}- Gating variable
*O*_{i}- Gating variable
- PCL
- Pacing cycle length
- PDE
- Partial differential equation
*R*- Gas constant
- SR
- Sarcoplasmic reticulum
*t*- Time
- T
- Temperature
*U*_{a}- Gating variable
*U*_{i}- Gating variable
*V*- Transmembrance potential
- VF
- Ventricular fibrillation
*X*- Na
^{+}, K^{+}, or Ca^{2+} *X*_{r}- Gating variable
*X*_{s}- Gating variable
- [X]
_{i} - Intracellular concentration of Na
^{+}, K^{+}, or Ca^{2+} - [X]
_{o} - Extracellular concentration of Na
^{+}, K^{+}, or Ca^{2+} *y*- A given gating variable
*y*_{∞}- Steady state of a given gating variable
*z*- Equal to 1, 1, or 2 for Na
^{+}, K^{+}, or Ca^{2+}, respectively

## METHODS

#### Human atrial model.

Cardiac cells are resistively connected by gap junctions between cells. Ignoring the detailed structure of real tissue, we consider a homogeneous continuous conduction model in which
Equation 1where *V* is the transmembrane potential,*C*
_{m} is the total capacitance, and *D* is the isotropic diffusion coefficient determined by gap junction resistance, surface-to-volume ratio, and membrane capacitance. Here we use *C*
_{m} = 100 pF and *D *= 0.001 cm^{2}/ms. *I*
_{ion} is the total ionic current from the Courtemanche model (3), which is given by
Equation 2
where *I*
_{Na} =*g*
_{Na}
*m*
^{3}
*hj*(*V*− *E*
_{Na}) and is the fast inward Na^{+}current; *I*
_{Ca,L} =*g*
_{Ca,L}
*dff*
_{Ca}(*V*− 65) and is the slow inward L-type Ca^{2+} current;*I*
_{Kr} =*g*
_{Kr}
*X*
_{r}(*V*− *E*
_{K})/1 + exp[(*V* + 15)/22.4] and is the rapid delayed outward K^{+} current;*I*
_{Ks} =*g*
_{Ks}
*X*
(*V*− *E*
_{K}) and is the slow delayed outward K^{+} current; *I*
_{to} =*g*
_{to}
*O*
*O*
_{i}(*V*− *E*
_{K}) and is the transient outward K^{+} current; *I*
_{Kur} = 0.005 + 0.05/{1 + exp[−(*V* − 15)/13]} and is the ultrarapid delayed outward K^{+} current; and*I*
_{K1} =*g*
_{K1}(*V* −*E*
_{K})/{1 + exp[0.07(*V* + 80)]} and is the time-independent K^{+} current. The ionic gating variables *m*, *h*, *j*,*d*, *f*, *f*
_{Ca},*X*
_{r}, *X*
_{s},*O*
_{a}, *O*
_{i},*U*
_{a}, and *U*
_{i} are all governed by Hodgkin-Huxley-type differential equations
Equation 3where *y* represents the appropriate gating variable. The detailed formulations of the ionic exchange currents, background currents, and their gating constants can be found in Courtemanche et al. (3). The driving force for *E*
_{Na},*E*
_{Ca}, and *E*
_{K} all have the formulation *E _{X}
* = (

*R*T/

*zF*)[log ([

*X*]

_{o}/[

*X*]

_{i})], where [

*X*]

_{o}and [

*X*]

_{i}are the extracellular and intracellular concentrations for

*X*= Na

^{+}, K

^{+}, or Ca

^{2+}and

*z*= 1, 1, or 2 for Na

^{+}, K

^{+}, or Ca

^{2+}, respectively.

The Courtemanche model includes the dynamics of intracellular Na^{+}, K^{+}, and Ca^{2+} and the dynamics of sarcoplasmic reticulum (SR) calcium storage and release processes, including intracellular Ca^{2+} uptake into the nonjunctional SR (NSR), Ca^{2+} leak from the NSR back into the cytosol, Ca^{2+} release from the junctional SR (JSR) into the cytosol, Ca^{2+} transfer from the NSR to JSR, and calcium buffering. The details of these equations can be found in Courtemanche et al. (3). In the Courtemanche model, [Na^{+}]_{i} and [K^{+}]_{i}never reach an equilibrium state for all basic pacing cycle lengths (PCLs). To eliminate this unphysiological behavior, we simply set [Na^{+}]_{i} = 11.2 mM and [K^{+}]_{i} = 139 mM. Clamping intracellular Na^{+} and K^{+} had no significant effect on the AP or APD restitution characteristics.

In the Courtemanche model, the maximum conductance of various ionic currents are *g*
_{Na} = 1.8 nS/pF,*g*
_{Ca,L} = 0.1238 nS/pF,*g*
_{Kr} = 0.0294 nS/pF,*g*
_{Ks} = 0.129 nS/pF,*g*
_{to} = 0.1652 nS/pF, and*g*
_{K1} = 0.09 nS/pF. We refer to these parameter settings as control conditions. To modify electrical restitution and spiral wave behavior, we altered*g*
_{Ca,L}, *g*
_{Kr}, and*g*
_{Ks} as described in the text.

#### Computer simulations.

Numerical simulations were performed in the isolated cell, one-dimensional (1-D) cable, and 2-D tissue. To simulate a single cell, we integrated the following ordinary differential equation
Equation 4where *I*
_{stim} is the external stimulus current. We used a square wave stimulus (29 pA/pF × 2 ms) at a constant frequency. We used a fourth-order Runge-Kutta method to integrate *Eq. 4
* with a fixed time step = 0.01 ms.

We used a 1-D cable of tissue for measuring APD and ERP restitution (see below), in which propagation is governed by the following partial differential equation
Equation 5
*Equation 5
* was solved using a forward Euler method with a time step = 0.01 ms and a space step = 0.025 cm.

For numerical simulation in 2-D tissue, the conventional Euler method to integrate *Eq. 1
* is computationally tedious and costly for a detailed cellular model. Therefore, we solved *Eq.1
* using the well-known operator splitting method. We split the nonlinear operator (*I*
_{ion} term) and the diffusion operator in *Eq. 1
* into two terms and then integrated the two terms separately and alternatively. We use an alternating-direction implicit method to integrate the partial differential equation (PDE) of the diffusion term and a time-adaptive second-order Runge-Kutta method (Δ*t*
_{min} ≤ 0.01 ms and Δ*t*
_{max} ≤ 0.1 ms) to integrate the ordinary differential equation of the reaction term with all the gating variable equations and the equations describing [Ca^{2+}]_{i}. The time step of integration of the PDE was set to Δ*t*
_{max} to keep all cells synchronized. The space steps (i.e., the size of the model cells) were set at d*x* = d*y* = 0.025 cm. With this approach, the integration speed increased more than 10-fold, with the relative error not exceeding 2% (25). The numerical methods and criteria for assuring numerical stability have been provided in detail previously (25). Because of differences in wavelength for different parameter settings, different tissue sizes (from 5 × 5 to 20 × 20 cm^{2}) were required to sustain spiral wave reentry in the 2-D simulations. This was achieved by varying the number of model cells (from 200 × 200 to 800 × 800), keeping individual cell size constant (0.025 cm). The actual array size used for each simulation is stated. All simulations were written in FORTRAN code and were run on DEC Alpha workstations.

#### Electrophysiological measurements.

APD restitution refers to the relationship between APD and the previous DI. In standard terminology, APD_{50} and APD_{90}are the times from the onset of depolarization until the cell repolarizes to 50% and 90% of the full AP amplitude, respectively. In this modeling study, we used a voltage threshold crossing method to measure APD, approximating APD_{90} with APD_{−62 mV}, defined as the portion of the cardiac cycle (in ms) during which *V* ≥ −62 mV, and DI as the portion during which*V* ≤ − 62 mV. APD_{50} was analogously approximated by APD_{−32 mV}, using −32 mV as the voltage threshold. APD_{−62 mV} and APD_{−32 mV} yielded values close to the true APD_{90} and APD_{50}values. APD_{−62 mV} and APD_{−32 mV} were measured in both single cells and 1-D tissue using an S1S2 pacing protocol in which a premature S2 stimulus scanning diastole was delivered after pacing with an S1 stimulus for 200 beats at a given CL. The pacing site for the 1-D cable simulations was located at the left boundary of the cable, and the recording site was located 1 cm away. CV restitution is defined as the relation between the propagated speed of the S2 AP and the corresponding DI.

ERP restitution was measured using an S1S2S3 pacing protocol. For each DI at which the S2 stimulus was delivered to obtain APD restitution, a third S3 stimulus was delivered at twice the diastolic threshold at progressively shorter DIs until the S3 stimulus failed to capture. The longest S2S3 interval that failed to capture was defined as the ERP of the S2 beat. ERP of the S2 beat was plotted against the DI of the S2 beat to obtain the ERP restitution curve.

Spiral wave reentry in 2-D tissue was initiated by cross-field stimulation of two successive perpendicular rectilinear waves. The tip of spiral wave was defined as the intersection point of the two successive contour lines of voltage in the simulated tissue corresponding to −30 mV measured at every 2 ms. Specifically, we first calculated the contour lines at every 2 ms, i.e.,*t _{n}
* = 2

*n*ms, where

*n*= 1, 2, etc., after the spiral wave was formed. We then calculated the intersection points between the two contour lines at time

*t*and

_{n}*t*

_{n}_{ + 1}for all

*n*. The movement of these intersection points over time traced the tip trajectories of the spiral waves.

## RESULTS

#### APD restitution in the single cell.

Under control conditions with parameters set to their original values (3), Fig. 1
*A*shows that the Courtemanche human atrial model had a resting membrane potential near −81 mV and diastolic [Ca^{2+}]_{i} ≈ 0.1 μM. AP shape and APD varied with the PCL, with a more prominent “spike-and-dome” appearance to the plateau at longer PCLs. APD_{−62 mV} was similar (264 ms) at PCLs of 1,000 and 600 ms, but shortened to 242 ms at a PCL of 400 ms. The 400 ms PCL was also associated with mild resting membrane potential depolarization, a modest increase in diastolic [Ca^{2+}]_{i}, and a decrease in peak systolic [Ca^{2+}]_{i}.

To reduce the slope of APD restitution, we decreased the maximum conductance (*g*
_{Ca,L}) of*I*
_{Ca,L}. At 65% and 90%*I*
_{Ca,L} block, respectively, Fig. 1, *B*and *C*, shows that rate-dependent modulation of the AP was virtually eliminated. Rate-dependent changes in the [Ca^{2+}]_{i} transient were preserved at 65% block but eliminated at 90% *I*
_{Ca,L} block, at which point *I*
_{Ca,L} was too small to trigger SR Ca^{2+} release. Figure 2 shows the corresponding APD_{−62 mV} and APD_{−32 mV}(approximating APD_{90} and APD_{50}, respectively) restitution curves for the control and 90%*I*
_{Ca,L} cases, measured at the three different PCLs (S1S1 = 1,000, 600, or 400 ms) using the S1S2 pacing protocol illustrated in Fig. 2
*A*. The shapes of the APD_{−62 mV} and APD_{−32 mV} curves differed markedly, with the slope of the APD_{−32 mV} restitution curve becoming negative at short DI for both the control (Fig. 2
*B*) and 90% *I*
_{Ca,L} block (Fig. 2
*C*) cases. For the control case, both APD_{−62 mV} and APD_{−32 mV} restitution curves were rate sensitive, whereas with 90%*I*
_{Ca,L} block, restitution curves were rate insensitive.

In the control case, the rate sensitivity of APD had interesting features. Figure3
*A* shows that both APD_{−62 mV} and APD_{−32 mV} of the S1 beats shortened as PCL decreased, but APD_{−62 mV} and APD_{−32 mV} of the S2 beats (S1S2 = 500 ms) were slightly longer at short PCL. The corresponding [Ca^{2+}]_{i}transients provide the explanation. At a shorter PCL, peak systolic [Ca^{2+}]_{i} for the S1 beat was decreased due to partial refractoriness of SR Ca release, which shortened APD. For the S2 beats, which all had the same 500-ms coupling interval, however, SR Ca^{2+} release and systolic [Ca^{2+}]_{i}were similar for all three PCLs. Minimal differences in APD of the S2 beats under these conditions were related to recovery processes of other ionic currents. With 90% *I*
_{Ca,L} block, no [Ca^{2+}]_{i} transients were triggered at any PCL, so that secondary rate-sensitive effects of the [Ca^{2+}]_{i} transient on APD were eliminated (Fig. 2
*C*).

#### APD and CV restitution in 1-D tissue.

In cardiac tissue, coupling of myocytes to neighboring cells generates diffusive currents that affect APD restitution and other electrophysiological properties. In addition, CV restitution is an important determinant of wave stability (26, 31) and cannot be measured in a single cell. To address these issues, we paced a 1-D cable (equivalent to a planar wave in 2-D or 3-D) and scanned diastole with premature S2 beats. APD restitution and CV restitution were measured 1 cm from the pacing site by plotting APD and local CV as functions of DI (see methods).

Figure 4, *left* and*middle*, compares APD_{−62 mV} and APD_{−32 mV} restitution obtained from the 1-D cable at PCLs of 1,000, 600, and 400 ms under four conditions: control,*I*
_{Ca,L} blocked by 46%,*I*
_{Ca,L} blocked by 65%, and*I*
_{Ca,L} blocked by 65% coupled with a ninefold increase in *I*
_{Ks} and *I*
_{Kr}. For control conditions, APD_{−62 mV} and APD_{−32 mV} were rate dependent and had markedly different shapes, as in the isolated cell. However, compared with the isolated cell, the restitution slope was steeper at short DIs (e.g., maximum slope 4.0 vs. 2.7 for APD_{−62 mV} in the single cell at a PCL of 1,000 ms). At a PCL of 400 ms, the slope of APD_{−32 mV}restitution at short DIs was steeper for APD_{−32 mV} than APD_{−62 mV}, but the converse was true for PCLs of 600 and 1,000 ms. When *I*
_{Ca,L} was blocked by 46% or more, APD restitution became independent of PCL. Maximum APD restitution slope decreased, but remained >1 at short DIs when*I*
_{Ca,L} was blocked by 46%. With 65% block, however, maximum APD restitution slope <1 for all DIs could be achieved by concomitantly increasing (ninefold)*I*
_{Ks} and *I*
_{Kr}.

Figure 5 shows that none of these interventions had any significant effect on CV restitution, as expected, because under normal excitability conditions, CV is determined primarily by *I*
_{Na} recovery properties, which were not modified.

#### ERP restitution versus APD restitution.

The APD_{−32 mV} and APD_{−62 mV}restitution curves in Figs. 2 and 4 had markedly different shapes and slopes, making the characterization of the APD restitution slope highly dependent on the repolarization level used to define APD. This ambiguity arose largely because PCL and premature stimuli had complex effects on the spike-and-dome morphology of the atrial AP plateau, which affects APD_{−32 mV} to a greater extent than APD_{−62 mV}. To develop a more consistent criterion for atrial restitution steepness, we examined ERP restitution, because wave break in fibrillation results directly from wave fronts encountering still-refractory tissue. ERP restitution was measured using an S1S2S3 protocol (see methods) and correlated with both APD restitution and spiral wave stability.

To validate this approach, we first tested the predictive value of ERP restitution in simulated 2-D ventricular tissue. Previous studies have established a close relationship between ventricular APD restitution slope and spiral wave behavior (15, 27), but whether ERP restitution slope is similarly predictive has not been reported to our knowledge. With the use of phase 1 of the Luo-Rudy ventricular AP model (18) in simulated homogenous 1-D and 2-D tissue, Fig. 6 shows that APD_{−32 mV} and APD_{−62 mV} curves are monophasic and qualitatively similar to each other. For both measurements, the range of DIs over which the APD restitution slope is >1 determined the spiral wave phenotype in homogenous 2-D tissue (Fig.6, last two columns on *right*), as reported previously. The ranges of DIs with slopes >1 for APD_{−62 mV} and ERP restitution in Fig. 6, *A–D*, are summarized in Table1. ERP restitution (*middle*) paralleled APD restitution and was similarly predictive.

In the human atrial model, Fig. 4, *right*, shows ERP restitution curves corresponding to the APD restitution curves in the*left* and *middle*. As expected intuitively, ERP restitution paralleled APD_{−62 mV} restitution more closely than APD_{−32 mV} restitution.

#### Spiral reentry behavior in 2-D atrial tissue and correlation with electrical restitution.

For the control parameter settings of the Courtemanche model, a spiral wave initiated with cross-field stimulation spontaneously broke up into multiple spiral waves that were self-sustaining provided that the 2-D tissue exceeded a critical size (Fig.7
*C*). For tissue sizes smaller than the critical size, breakup either did not occur, or occurred transiently, with the tissue eventually becoming quiescent (Fig. 7,*A* and *B*). The critical tissue size was ∼16.5 × 16.5 cm^{2}, which was obtained empirically by numerical simulation. For 10 × 10- and 15 × 15-cm^{2} tissues (Fig. 7, *A* and *B*), the spiral wave tip moved irregularly toward the boundary before the first wave break occurred and disappeared from the boundary, terminating reentry. In these two cases, the tissue became quiescent after 5 and 8 s, respectively, with the exact duration depending strongly on the conditions used to initiate the spiral wave. Although no wave break occurred in Fig. 7, *A* and *B*, the contours of both the wave back and wave front were irregular, reflecting their dynamical instability. For 20 × 20-cm^{2} tissue (Fig.7
*C*), the initiated spiral wave broke into complex multiple wavelets after several seconds, which were sustained throughout the 25-s simulation.

Figure 8 shows that the character of spiral wave reentry changed when electrical restitution was altered by modifying *I*
_{CaL}, *I*
_{Kr}, and*I*
_{Ks} in the Courtemanche human atrial model. Four conditions were studied: control, *I*
_{Ca,L} blocked by 46%, *I*
_{Ca,L} blocked by 65%, and*I*
_{Ca,L} blocked by 65% coupled with a ninefold increase in *I*
_{Ks} and *I*
_{Kr}. The slope of ERP restitution was >1 at short DIs in the first two conditions (Fig. 4, *A* and *B*) and <1 at all DIs in the second two cases (Fig. 4, *C* and *D*). Figure8,* A–D*, shows the corresponding spiral wave phenotypes for these different parameter settings of the Courtemanche human atrial model in Fig. 4, *A–D*, using sufficiently large tissue size to permit self-sustaining reentry. In each case, a single spiral wave was initiated using cross-field stimulation, and its subsequent evolution was tracked. As shown in the snapshots of voltage at various times after spiral wave initiation (*left*) and the traces of the spiral wave tip trajectory (*right*), the initiated spiral wave spontaneously broke up into complex multiple spirals after several seconds, and the tip trajectory was fully irregular for the control case (Fig. 8
*A*). In this case, maximum ERP restitution slope was 2.34 and was >1 over a 33-ms range of DIs. The pseudo-ECG showed a fibrillation-like pattern. When the maximum slope of ERP restitution was decreased to 1.68 and the range of DIs for which ERP restitution slope was >1 narrowed to ∼9 ms by blocking*I*
_{CaL} by 46% (Fig. 8
*B*), the spiral wave meandered chaotically (as seen from the irregularities in the tip trajectory) but did not break up. The corresponding pseudo-ECG showed polymorphic tachycardia rather than fibrillation. Further block of*I*
_{Ca,L} block to 65% led to further reductions in the maximum slope of ERP restitution to 0.66 so that the slope was now <1 over all DIs (Fig. 8
*C*). This produced a quasiperiodically meandering spiral wave (as seen by the smooth flower pattern of the tip trajectory), with the pseudo-ECG now becoming nearly monomorphic. To produce a completely stationary spiral wave, however, blockade of *I*
_{Ca,L} alone was not sufficient. Concomitantly increasing *I*
_{Kr} and*I*
_{Ks} at least ninefold shortened APD sufficiently so that the tip of the spiral wave did not encounter incompletely repolarized tissue, and the spiral wave then became completely stationary (Fig. 8
*D*). The maximum slope of ERP restitution decreased further to 0.08, and the pseudo-ECG now showed monomorphic tachycardia (atrial flutter-like). The tip trajectory was a circle.

Table 1 summarizes the value of maximum slope of ERP and APD_{−62 mV} restitution, the range of DI with ERP and APD_{−62 mV} restitution slope >1, the range of CLs during reentry, and the range of DIs visited during reentry in simulated ventricular and atrial tissues for the four cases (Figs. 6 and 8). During reentry in the breakup and chaotic hypermeander regimes, DIs visited regions corresponding to restitution slope >1 as well as DIs with restitution slope <1. Spiral wave stability correlated strongly with decreases in both the maximum slope of ERP or APD_{−62 mV} restitution and the range of DIs with ERP or APD_{−62 mV} slope >1.

## DISCUSSION

In this study, we investigated the role of electrical restitution on spiral wave stability in simulated homogeneous atrial and ventricular 2-D tissue using the Courtemanche atrial model and the Luo-Rudy phase 1 ventricular AP model, respectively. Our major conclusions can be summarized as follows. First, in the physiologically detailed (including intracellular Ca^{2+} dynamics) Courtemanche atrial AP model, APD restitution curves are qualitatively as well as quantitatively sensitive to the level of repolarization used to define APD. This arises primarily from the complex rate sensitivity of the spike-and-dome configuration of the atrial AP plateau and effects of the intracellular Ca^{2+} transient on ionic currents affecting the APD. These factors tend to have larger effects on APD_{−32 mV} (∼APD_{50}) than on APD_{−62 mV} (∼APD_{90}), which leads to ambiguities in attempting to correlate APD restitution slope with spiral wave behavior. In addition, considerable regional variations in atrial AP morphology, particularly with respect to the spike-and-dome feature, are likely to add further site-dependent variability to atrial APD restitution measurements. ERP restitution and CV restitution in atrial tissue, however, can be characterized with similar reliability as in ventricular tissue. Second, to circumvent ambiguities in atrial APD restitution, we reasoned that wave breaks that promote fibrillation occur primarily when waves encounter areas of increased refractoriness. Thus ERP restitution slope might be an equivalent or superior measure of dynamic wave stability than APD restitution slope. We validated this concept first in ventricular tissue using the Luo-Rudy guinea pig ventricular AP model and then applied the same test to atrial tissue. We found that ERP restitution slope correlated better with APD_{−62 mV} restitution than with APD_{−32 mV}restitution and was highly predictive of spiral wave behavior. For the tissue with steep ERP restitution slope (>1), an initiated spiral wave either broke up or exhibited chaotically hypermeandering tip motion, depending on the range of DIs with ERP restitution slope >1. When ERP restitution slope was flat (<1 for all DIs), a single spiral wave remained intact with either quasiperiodic meandering or stable tip motion, as we and others have previously reported for APD_{90}restitution in ventricular AP models (15, 27). These findings support the general validity of ERP restitution slope as a measure of dynamic wave stability in both ventricular and atrial tissue. A caveat, however, is that in our study, we modulated ERP slope specifically by altering *I*
_{Ca,L},*I*
_{Kr}, and *I*
_{Ks}. Similar effects on ERP restitution steepness can be achieved by modifying other combinations of different currents, and we did not directly demonstrate their equivalence with respect to spiral wave behavior. However, previous simulation studies have demonstrated that APD restitution slope is a robust global system parameter controlling spiral wave stability, i.e., that the manipulations of specific ionic currents used to alter restitution slope are not critically important, only their net effect on restitution slope (28). This has been demonstrated in a wide variety of models, ranging from simple phenomenological two- and three-variable models, such as the Karma models (15), to intermediate models, such as the Beeler-Reuter model (5), to physiologically detailed models, such as the various Luo-Rudy versions (18).

The advantage of using ERP restitution in place of APD restitution is that ambiguities about restitution slope arising from effects of the level of repolarization chosen to define APD are eliminated. This is important not only in the atrium, but also potentially in regions of ventricle exhibiting spike-and-dome AP morphology, such as the epicardium. In addition, ERP restitution can be obtained using extracellular pacing/recording electrodes, without an absolute need for accurate measurement of the transmembrane AP by microelectrode, optical dye, or monophasic AP recording techniques. This should facilitate electrical restitution measurements in clinical studies, in which the monophasic AP catheter is the only option for measuring APD restitution and is more cumbersome than extracellular electrogram recordings. Although the DI cannot be measured directly from extracellular electrograms, it can be readily estimated from the ERP during the basic S1 pacing (i.e., the ERP of the S2 beat). On the other hand, a disadvantage of ERP restitution is that it requires a pacing intervention (the S3 beat) at each site at which ERP is to be measured, making high resolution spatial mapping of electrical restitution impractical and inferior to optical mapping of APD restitution in the experimental setting. We are currently investigating whether the need for the S3 beat can be eliminated by using activation-recovery intervals to estimate local ERP (20).

There are several limitations in this modeling study. First, we modified *I*
_{Ca,L}, *I*
_{Ks}, and*I*
_{Kr} to produce altered ERP restitution slopes and spiral wave behaviors without a primary consideration of whether these interventions had existing physiological counterparts. In principle, however, similar changes could be reproduced experimentally by drugs (which may not exist currently but potentially could be developed), because they all involved modulating simulated real ionic currents. Second, to achieve sustained spiral reentry for the control case, we also had to use a larger tissue size than the realistic atrium. However, this is consistent with the observation that when AF is induced in the normal atrium, it is typically nonsustained. Third, APD, ERP, and CV restitution are not unique relationships, but depend on factors such as wave front curvature and memory. Thus electrical restitution experienced by wave fronts during fibrillation may be different from electrical restitution properties measured by extrastimulus or rapid pacing methods. Although electrical restitution measured by the extrastimulus method is only an approximation of that during fibrillation, our findings substantiate that from a practical standpoint the relationship is close enough to provide useful information for assessing wave stability. This is illustrated in Table1, which demonstrates if the range of DIs visited during spiral wave reentry corresponded to the range of DIs with APD restitution slope >1, the spiral wave reentry was unstable (i.e., in the hypermeander or breakup regimes).

A fourth limitation of this study relates to tissue characteristics. To maximize computational efficiency, we used the minimum tissue size required to support sustained reentry, which varied for different spiral wave phenotypes. In this model, however, we substantiated that increasing tissue size beyond the minimum size did not change the spiral wave phenotype, suggesting that the primary role of tissue borders was to extinguish wave fronts rather than create new ones. Below the critical size, the phenotype was also the same, but reentry was nonsustained, consistent with the critical mass hypothesis validated in experimental studies (17). In addition, to isolate dynamically induced effects of electrical restitution on wave stability from the effects of preexisting heterogeneities, we simulated 2-D homogeneous atrial and ventricular tissue rather than 3-D heterogeneous tissue. It is well known that preexisting electrophysiological/anatomic features in 3-D tissue, such as tissue anisotropy, fiber rotation, and complex anatomic structures, affect the stability of the spiral wave reentry (8, 12, 14, 34). In addition, the atrial AP varies substantially from region to region, producing preexisting dispersion of refractoriness affecting spiral wave stability (10, 30). Although our findings indicate that electrical restitution in the atrium is a potentially important determinant of wave stability, they do not address the relative importance of restitution-based wave break compared with heterogeneity-based wave break in AF. At the present time, the role of restitution-based wave break in the maintenance of AF remains controversial, with some investigators believing that wavebreak is an epiphenomenon (2) rather than the engine of fibrillation, as proposed in the original multiple wavelet hypothesis (21). If wave break is the engine of AF, however, our previous modeling work (35) has documented a synergistic relationship between electrical restitution parameters and preexisting tissue heterogeneity in promoting wave break. Specifically, we found that for a given level of preexisting electrophysiological/anatomic heterogeneity, wave break was more easily induced as restitution-based dynamic instability increased. Together, these observations suggest that reducing dynamic instability by flattening APD or ERP restitution may have therapeutic potential for preventing AF and VF. This strategy will need to be tested in future studies using realistic heterogeneous models of atrial and ventricular tissue and then validated in experimental AF and VF models.

## Acknowledgments

We thank Elizabeth M. Cherry and Flavio Fenton for helpful discussions.

## Footnotes

This work was supported by National Heart, Lung, and Blood Institute Specialized Center of Research in Sudden Cardiac Death P50 HL-52319, American Heart Association Western States Affiliate Beginning Grant-in-Aid 0060083Y, and by Laubisch and Kawata endowments.

Address for reprint requests and other correspondence: F. Xie, Div. of Cardiology, 47-123 CHS, UCLA School of Medicine, Los Angeles, CA 90095-1679 (E-mail:fxie{at}mednet.ucla.edu).

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.First published April 11, 2002;10.1152/ajpheart.00898.2001

- Copyright © 2002 the American Physiological Society