## Abstract

Ventricular fibrillation (VF) is one of the main causes of death in the Western world. According to one hypothesis, the chaotic excitation dynamics during VF are the result of dynamical instabilities in action potential duration (APD) the occurrence of which requires that the slope of the APD restitution curve exceeds 1. Other factors such as electrotonic coupling and cardiac memory also determine whether these instabilities can develop. In this paper we study the conditions for alternans and spiral breakup in human cardiac tissue. Therefore, we develop a new version of our human ventricular cell model, which is based on recent experimental measurements of human APD restitution and includes a more extensive description of intracellular calcium dynamics. We apply this model to study the conditions for electrical instability in single cells, for reentrant waves in a ring of cells, and for reentry in two-dimensional sheets of ventricular tissue. We show that an important determinant for the onset of instability is the recovery dynamics of the fast sodium current. Slower sodium current recovery leads to longer periods of spiral wave rotation and more gradual conduction velocity restitution, both of which suppress restitution-mediated instability. As a result, maximum restitution slopes considerably exceeding 1 (up to 1.5) may be necessary for electrical instability to occur. Although slopes necessary for the onset of instabilities found in our study exceed 1, they are within the range of experimentally measured slopes. Therefore, we conclude that steep APD restitution-mediated instability is a potential mechanism for VF in the human heart.

- reentrant arrhythmias
- human ventricular myocytes
- restitution properties
- spiral waves
- computer simulation

one of the most extensively investigated hypotheses for ventricular fibrillation (VF) is the so-called restitution hypothesis. In its initial form the hypothesis stated that if the action potential duration (APD) restitution curve has a maximum slope steeper than 1, it will lead to APD alternans (16, 41). In tissue this APD alternans can result in the fragmentation of a spiral wave, leading to fibrillation-like excitation patterns (21, 22, 45, 48). Modeling studies have confirmed that a steep restitution curve indeed promotes instability. However, modeling studies have also shown that the criterion of an APD restitution slope >1 is an oversimplification that only holds for very simple models. In more realistic and complex models, it has been shown that other factors such as short-term cardiac memory, electrotonic interactions between cells, conduction velocity (CV) restitution, the range of diastolic intervals (DIs) over which restitution is steeper than 1, and the range of DIs visited during spiral wave rotation all play an important role in determining whether alternans and spiral breakup will occur (5, 7, 10, 11, 42, 48, 64). In an extensive study of restitution-induced instability, Cherry and Fenton (6), for example, showed that because of strong electrotonic interactions and gradual CV restitution, spiral breakup may not occur even if the APD restitution curve has slope as much as 2.

The restitution hypothesis has promoted a large number of experimental studies. An important goal of these studies was to measure restitution curves and their slopes. Experiments in animal hearts have shown that the maximum slope of the restitution curve can indeed be >1; for example, for pig heart, slopes of 1.29 (27); for dog heart, slopes of 1.13 (25); and for rabbit heart, slopes of 1.3 (18) up to 2.3 (2) were found. However, slopes of <1 have also frequently been reported (18, 27). In experiments it was also found that factors other than mere APD restitution slope are important in determining whether alternans or fibrillation occurs. Wu et al. (66) and Banville and Gray (2) showed that whether VF occurred depended both on APD and CV restitution. Banville et al. (1) also demonstrated that cardiac memory effects play a role in whether alternans and VF occur.

Early studies of APD restitution in the human heart suggested that the maximum slope in humans could be close to or slightly larger than 1 (37, 59). More recently, Pak et al. (43) measured human APD restitution in the right ventricular apex and outflow tract. They found that restitution slopes can be considerably steeper than 1, reporting values of up to 1.28 for the apex and of up to 3.78 for the outflow tract. In addition, they found a positive correlation between arrhythmia inducibility and both steepness and dispersion of the restitution slopes. Yue et al. (67) measured restitution on the left and right ventricular endocardium at 32 sites, reporting restitution slopes varying from site to site and ranging from 0 to 2. The most extensive study on human APD restitution has recently been performed by Nash et al. (38–40). They measured restitution at 256 points covering the complete ventricular epicardium of 14 patients, showing that restitution is distributed heterogeneously but is regionally organized over the epicardial surface, with maximum slopes varying from almost 0 to >3, and a mean value of 1.1.

These recent data imply the possibility of steep restitution-induced instabilities in human cardiac tissue. However, as was mentioned above, the question of whether such instabilities actually do occur also depends on other factors, such as CV restitution, electrotonic interactions and cardiac memory, range of DIs over which the curve is steep, and DIs visited during spiral wave rotation. Because of the substantial limitations to experimental studies involving human cardiac tissue, alternative methods are of great interest. In this study we will investigate the conditions for electrical instability in human ventricular tissue by using the method of mathematical modeling.

For this purpose we developed a new version of our human ventricular cell model (61). The new model is based on recent experimental restitution data (39) and incorporates an improved description of intracellular calcium dynamics, by including subspace calcium dynamics that controls L-type calcium current and calcium-induced calcium release (CICR), by modeling CICR with a four-state Markov model for the ryanodine receptor, and by incorporating both fast and slow voltage-gated inactivation of the L-type calcium current. We use our new model to study the onset of alternans in a single cell, for waves propagating in a ring of cardiac tissue (a model for reentry around a fixed path), and for two-dimensional reentry (spiral waves).

In our earlier study, we found that the recovery dynamics of the fast sodium current (*I*_{Na}) are considerably slower than what is often assumed in other modeling studies. *I*_{Na} dynamics are an important determinant of CV restitution and of spiral wave rotation speed. Cherry and Fenton (6) have demonstrated that CV restitution plays an important role in the occurrence of electrical instability, whereas we (63) have previously shown that the speed of spiral wave rotation and hence the DIs visited during spiral wave rotation influence whether alternans instability occurs.

Therefore, we focus on the effect of *I*_{Na} recovery dynamics in combination with APD restitution on alternans and spiral breakup. In particular, we use our new model to reproduce a representative range of the experimentally measured restitution curves recently reported by Nash et al. (39), with slopes varying from 0.7 to 2. We combine these different restitution slopes with three types of *I*_{Na} recovery dynamics: slow, medium, and fast recovery. We find that slower *I*_{Na} recovery dynamics suppress steep APD restitution-mediated instability and spiral breakup. Under slow *I*_{Na} recovery dynamics, maximum restitution slopes significantly steeper than 1 are required to generate electrical instability and breakup. However, such slopes are within the range of restitution slopes reported by Nash et al. (39). Therefore, we conclude that steep restitution-mediated alternans and spiral breakup may exist in the human heart.

## MATERIALS AND METHODS

### Model Development

We developed a new version of our human ventricular cell model. In this new model we included a more extensive description of intracellular calcium handling, incorporating a subsarcolemmal space, describing CICR with a Markov-state model for the ryanodine receptor, including both fast and slow voltage inactivation for the L-type calcium current, and applying some minor changes to parameter values and to slow delayed rectifier time dynamics. Important changes are discussed in detail below. Detailed equations can be found in the appendix; the new default parameter setting of our model can be found in Table 1.

#### Calcium dynamics.

A schematic representation of the calcium dynamics in our new model is given in Fig. 1. The model now contains a description of calcium dynamics in the subspace (Ca_{ss}), cytoplasm (Ca_{i}), and sarcoplasmic reticulum (Ca_{SR}). Calcium in the subspace is buffered, and a diffusion flux is added to allow calcium released in the subspace to flow to the bulk cytoplasm. L-type calcium current and the ryanodine calcium release current now inject calcium into the subsarcolemmal subspace, and, in turn, their dynamics are influenced by the subspace calcium concentration. Sarcolemmal Ca^{2+} pump current (*I*_{pCa}), background Ca^{2+} current (*I*_{bCa}), and Na^{+}/Ca^{2+} exchanger current (*I*_{NaCa}) still depend on bulk cytoplasmic calcium.

Figure 2 shows a steady-state 1-Hz epicardial action potential for the new version our human ventricular cell model (Fig. 2*A*) together with the calcium transients occurring in the subspace (Fig. 2*C*) and bulk cytoplasm (Fig. 2*B*). It can be seen that the calcium transient in the subspace is much larger (peak amplitude of 1.77 vs. 0.0013 mM), faster (time to peak ∼2 vs. ∼35 ms), and shorter (duration of ∼35 vs. ∼400 ms) than that in the cytoplasm. Maximum upstroke velocity of the action potential is 289 mV/ms, the plateau level is 24.9 mV, and resting membrane potential is −85.8 mV, all very similar to the previous version of our model. However, APD at 90% repolarization (APD_{90}) now is 306 ms and was 276 ms in the previous version of our model. We decided for this change to let our model APD lie more in the midrange of experimentally recorded APDs in single cell [336, 298, 360, and 310 ms (28–31)] and tissue preparations [351 and 378 ms (8, 9, 14)].

#### L-type calcium current.

L-type calcium current (*I*_{CaL}) (see *Eq. 6* of appendix) now injects calcium into the subspace and in turn is inactivated by the Ca_{ss} via the inactivation gate *f*_{cass} (Fig. 3, *A* and *B*). *I*_{CaL} voltage-clamp experiments indicate the presence of both a fast and slow voltage-dependent recovery process (32, 35, 46). Therefore, we incorporate both a slow voltage inactivation gate *f* and a fast voltage inactivation gate *f*_{2}, similar to the approach taken by Hund and Rudy (19). Figure 3*C* shows the steady-state inactivation curves of *f* and *f*_{2}, and Fig. 3*D* shows the time constant curves of *f* and *f*_{2} together with experimental data. Note that the inactivation rate of *f*_{2} is slower than that of *f*_{cass} but much faster than that of *f*. Note also that as opposed to the *f* gate, but similar to the *f*_{cass} gate, the *f*_{2} gate inactivates incompletely. These two properties enable the *f*_{2} gate to take over initial inactivation from *f*_{cass}, which is necessary because the short duration of the subspace calcium transient leads to recovery of *f*_{cass} during the action potential. In Fig. 3*E*, we show the voltage dependence of our new *I*_{CaL} together with experimental data, demonstrating a good agreement between model and experiment (34).

#### CICR current.

In our new model, we use a reduced version of the Markov-state ryanodine receptor model developed by Shannon et al. (54) and Stern et al. (57) to describe CICR. Figure 4*A* shows the original four-state Markov model of the ryanodine receptor developed by Shannon et al. (54) and Stern et al. (57). It can be seen that the model incorporates both the influence of subspace calcium concentration (the trigger) and sarcoplasmic reticulum calcium concentration (the load) on receptor opening and closing dynamics by making transition rates depend on these concentrations.

For computational efficiency, we reduced the dimensionality of the Markov model by using quasi-steady-state assumptions to the following set of equations (1) (2) (3) where R̄ = R + O, with R the resting closed state and O the open conducting state; and where *I*_{rel} is the CICR current.

We determined the gain function for the dependence of calcium release on SR calcium load. (Gain is defined as the CICR flux normalized by the L-type calcium current flux triggering it.) Figure 4*B* shows gain in our model as a function of free calcium in the sarcoplasmic reticulum. For qualitative comparison, experimental data from Shannon et al. (53) in rabbit ventricular myocytes are added; unfortunately we could not find data on the gain function in human ventricular myocytes. We can see that our model generates a nonlinear response function, as is the case for the experimental data.

### Numerical Methods

Action potential generation in single cells was described using the following differential equation (23) (4) where *C*_{m} is the membrane capacitance, *V* is the transmembrane potential, *I*_{stim} is the externally applied transmembrane current, and *I*_{ion} is the sum of all transmembrane ionic currents. For *I*_{ion}, we use our new version of the human ventricular cell model.

Similarly, action potential generation and propagation in one-dimensional cables and rings and two-dimensional sheets of cardiac tissue were described by using the following parabolic reaction diffusion equation (23) (5) where *D* is a tensor describing the conductivity of the tissue and ∇ is the one- or two-dimensional gradient operator. In one dimension, *D = 0.00154* cm^{2}/ms; in two dimensions, *D*_{ij} = 0.00154 cm^{2}/ms for *i = j*, and *D*_{ij} = 0 cm^{2}/ms for *i* ≠ *j*. With these values we obtain a maximum planar CV of 68 cm/s, which is the velocity found for conductance along the fiber direction in human myocardium (60).

Physical units used in our model are as follows, time (*t*) in milliseconds, voltage (*V*) in millivolts, current densities (*I*_{X}) in picoamperes per picofarad, and ionic concentrations (*X*_{i}, *X*_{o}) in millimoles per liter.

For single cell simulations, forward Euler integration with a time step of Δ*t* = 0.02 ms was used to integrate *Eq. 4*. For one- and two-dimensional computations, the forward Euler method was used to integrate *Eq. 5* with a space step of Δ*x* = 0.25 mm and a time step of Δ*t* = 0.02 ms . In all cases, the Rush and Larsen integration scheme (52) was used to integrate the Hodgkin-Huxley type equations for the gating variables of the various time-dependent currents (*m*, *h*, and *j* for *I*_{Na}; *r* and *s* for *I*_{to}; *x*_{r1} and *x*_{r2} for *I*_{Kr}; *x*_{s} for *I*_{Ks}; and *d*, *f*, *f*_{2}, and *f*_{cass} for *I*_{CaL}).

We measure APD at either 50% or 90% repolarization and will refer to these values as APD_{50} and APD_{90}, respectively. APD_{90} values will be used, unless stated otherwise. APD_{50} values are used to allow comparison between experimental activation recovery interval (ARI) data obtained by Nash et al. (38–40). It has been shown that ARI (measured as the interval from the minimum d*V*/d*t* during the QRS wave to the maximum d*V/*d*t* of the T-wave) typically corresponds with APD of ∼50% repolarization (17).

In single cells, we apply both the standard S1-S2 and the dynamic protocols to determine APD restitution. For the S1-S2 protocol, 10 S1 stimuli are applied at a specified basic cycle length (BCL) followed by a single S2 extrastimulus delivered at some DI after the action potential generated by the last S1 stimulus. An APD restitution curve is generated by decreasing DI and plotting the APDs generated by the S2 stimuli against the preceding DIs. In this study, we use APD_{50} to allow comparison with the experimental studies of Nash et al. (39). For the dynamic restitution protocol, a series of 50 stimuli is applied at a specified BCL, after which cycle length is decreased. The APD restitution curve is obtained by plotting the final APD for each BCL against the final DI.

In cables we apply the dynamic restitution protocol to determine both dynamic APD restitution and CV restitution. We do so by pacing one end of the 800-cell-long cable at a certain BCL until a steady-state APD and CV are reached, after which the cycle length is decreased. In rings, APD and CV restitution can be obtained in the absence of external pacing. Here the restitution protocol is started by applying a single external stimulus generating a wave that is allowed to propagate only in one direction and that will result in the repeated rotation of a wave along the ring. The cycle length of excitation is determined by the ring length. APD and CV can be determined as a function of ring length and hence cycle length and DI. We start with a ring of 1,400 cells (35 cm) and stepwise reduce its length to decrease cycle length and hence obtain a restitution curve.

We use two-dimensional tissue sheets of 1,000 × 1,000 cells (space step Δ*x* = 0.25 mm). In two dimensions, spiral waves are generated by first applying an S1 stimulus producing a planar wavefront propagating in one direction; then, when the refractory tail of this wave crosses the middle of the medium, an S2 stimulus is applied, generating a second wavefront perpendicular to the first. This produces a second wavefront with a free end around which it curls, forming a spiral wave. Stimulus currents lasted for 2 (S1) and 5 (S2) ms and were twice the diastolic threshold.

Single-cell, cable, and ring simulations were coded in C++ and run on a single processor of a Dell 650 Precision Workstation (dual Intel xeon 2.66 GHz). Two-dimensional simulations were coded in C++ and MPI and were run on 20 processors of a Beowulf cluster consisting of 14 Dell 650 Precision Workstations (dual Intel xeon 2.66 GHz).

## RESULTS

### APD Restitution Properties of Human Ventricular Cell Model

Our new model reproduces a representative range of restitution curves recently reported by Nash et al. (38–40). Figure 5, *A–D*, shows restitution curves for four different parameter settings of our model together with four different experimental restitution curves of Nash et al. Experimental and model restitution curves were obtained using an S1-S2 protocol at a BCL of 600 ms. In experiments, ARIs were measured, which correspond to APD_{50} (see methods). Model restitution curves are plotted for both APD_{50} and APD_{90} levels. Note that there is good agreement between experimental and model curves for DIs of 200 ms and smaller, i.e., for DIs that determine the maximum slope of the restitution curve and that determine the onset of instability. For this study we use four parameter settings corresponding to four different slopes of the restitution curves. For *setting 1*, the maximum restitution slope (*S*_{max}) is 0.7 (Fig. 5*A*); for *setting 2*, it is 1.1 (Fig. 5*B*); for *setting 3*, it is 1.4 (Fig. 5*C*); and for *setting 4*, it is 1.9 (Fig. 5*D*). Table 1 lists the complete default parameter setting, and Table 2 lists the parameter values varied to obtain these different restitution slopes.

### Conditions for Alternans in Human Ventricular Tissue

#### Restitution and alternans in single cells.

It has been suggested by Koller et al. (25) that the slope of the dynamic restitution curve is a better predictor for alternans instability than that of the S1-S2 restitution curve. Therefore we also determined dynamic restitution for our four different model settings. Figure 5, *E–H*, shows restitution curves obtained using a dynamic restitution protocol and measuring APD_{90}. Maximum restitution slopes for our four setting are now 0.8 (Fig. 5*E*), 1.1 (Fig. 5*F*), 1.4 (Fig. 5*G*), and 1.8 (Fig. 5*H*), respectively. This agrees with experimental data from Pak et al. (43), who demonstrated similar maximum slopes for S1-S2 and dynamic restitution in the human ventricles.

In Fig. 5, *I–L*, we show APD restitution for the same dynamic restitution protocol but now with APD plotted vs. period. We can see that for *setting 1* (Fig. 5*I*), no alternans occurs. For *setting 2* (Fig. 5*J*), APD alternans occurs, resulting in the splitting of the restitution curve in an upper and lower arm representing the long and short action potentials. Alternans occurs for periods between 236 and 182 ms, corresponding to DIs between 51 and 24 ms for which the slope of the restitution curve in Fig. 5*F* exceeds 1. Alternans reaches a maximum amplitude of 48 ms at a period of 226 ms and then decreases in amplitude and disappears for periods <182 ms, qualitatively similar to what occurs in the model of Fox et al. (12). For *setting 3* (Fig. 5*K*), APD alternans occurs for periods <280 ms, corresponding to DIs < 74 ms, for which restitution slope in Fig. 5*G* is >1, and reaches a maximum amplitude of 100 ms. For *setting 4* (Fig. 5*L*), alternans occurs for BCL of 320 ms and smaller, corresponding to DIs < 94 ms, for which slope in Fig. 5*H* is >1, and reaches a maximum amplitude of 162 ms. Our finding that alternans amplitude increases with restitution slopes is consistent with experimental results (24, 26).

#### Restitution and alternans in a ring.

The occurrence of alternans in a single cell does not necessarily imply that alternans and spiral breakup occurs in tissue. In tissue, additional factors besides APD restitution, such as electrotonic interactions and CV restitution, the range of DIs for which restitution is steep, and the range of DIs visited during spiral wave rotation, all play a role in determining whether alternans occurs (5–7, 10, 11, 42, 48, 64).

In our previous study (61), we found that recovery of inactivation of human cardiac fast sodium channels is considerably slower than that of the phase-one Luo-Rudy model (33), which is used in many cardiac tissue models. *I*_{Na} dynamics may effect period of spiral wave rotation, CV restitution, and spiral wave meander pattern, all of which are known to affect the occurrence of alternans instability. Therefore, in tissue simulations, we investigate the combined effect of APD restitution and *I*_{Na} recovery dynamics on dynamical instability. To do so, we considered three different descriptions for the fast *I*_{Na} dynamics: the default formulation of our model (61) (referred to as standard *I*_{Na} setting); an *I*_{Na} formulation using *m*, *h*, and *j* gate dynamics as in the phase-one Luo-Rudy model (33) (referred to as the LR *I*_{Na} setting); and an intermediate setting. The intermediate setting is similar to the standard *I*_{Na} setting, except that τ_{j,intermediate} = 0.5 × τ_{j,standard}. Note that changes in our model's fast *I*_{Na} dynamics have virtually no effect on the single-cell APD restitution properties shown in Fig. 5.

To examine the combined effect of APD restitution properties and *I*_{Na} recovery dynamics on the occurrence of electrical instability in tissue, we combined the four different APD restitution parameter settings with the three different *I*_{Na} dynamics discussed above. This resulted in a total of 12 different model settings. Using these settings, we studied the stability of wave propagation in a ring, which is a one-dimensional model for reentry along a fixed path. Note that it was shown in Ref. 6 that the conditions for instability in a ring are the same as the conditions for spiral breakup in two dimensions.

In Fig. 6, we show APD as a function of period in the 12 different model settings for a wave circulating on a ring. Figure 6*A* shows restitution for parameter *setting 1* (*S*_{max} = 0.7) combined with either the standard, intermediate, or LR *I*_{Na} dynamics (curves are superimposed). We can see that for none of these settings alternans occurs. The *I*_{Na} dynamics, however, do affect the period for which block occurs: for standard *I*_{Na} dynamics, block occurs for periods <216 ms; for intermediate dynamics, block occurs for periods <202 ms; and for LR dynamics, block occurs for periods <170 ms. In Fig. 6*B* we can see that for *setting 2* (*S*_{max} = 1.1) for standard and intermediate *I*_{Na} dynamics, no alternans occurs, but for LR *I*_{Na} dynamics, alternans occurs for periods between 255 and 170 ms, with a maximum amplitude of 25 ms.

In Fig. 6*C*, we see that for *setting 3* (*S*_{max} = 1.4), no alternans occurs for standard *I*_{Na} dynamics. Both the intermediate and LR dynamics result in alternans. For intermediate *I*_{Na} dynamics, alternans occurs for periods <280 ms, with a maximum amplitude of ∼50 ms; for LR *I*_{Na} dynamics, alternans occurs for periods <300 ms, with a maximum amplitude of ∼90 ms. Finally, in Fig. 6*D* we see that for parameter *setting 4* (*S*_{max} = 1.9), alternans occurs in all three cases. For the standard dynamics, alternans occurs for periods <321 ms, and for the intermediate and LR settings for periods <330 ms. Maximum alternans amplitude increased from around 75 to 125 ms and to 150 ms when going from standard to intermediate to LR I_{Na} dynamics.

Our results show that for fast (LR) *I*_{Na} dynamics, alternans instability occurs for the same APD restitution slopes in a single cell as in a ring of tissue. For slow (standard) and intermediate *I*_{Na} dynamics, this is not the case. For standard and intermediate *I*_{Na} dynamics, a steeper APD restitution slope is required in a ring than in a single cell to generate alternans. These results suggest that faster sodium recovery dynamics lead to earlier occurrence and larger amplitude of alternans and that slower sodium recovery suppresses alternans.

#### Spiral wave dynamics and alternans in a plane.

In Fig. 7 we show snapshots of spiral wave dynamics for the 12 different model settings that we studied in the previous section. The three columns represent the three different *I*_{Na} dynamics (standard, intermediate, LR). The four rows represent the four different parameter settings resulting in the four different APD restitution slopes. If spiral breakup occurs, the time of first wave break is indicated in the top right corner. Spiral wave dynamics were simulated for 4 s. If no spiral breakup occurred in this time frame, spiral wave dynamics were assumed to be stable.

We can see that for parameter *setting 1* (*S*_{max} = 0.7, first row), stable spiral wave dynamics occur for all three *I*_{Na} dynamics. For *setting 2* (*S*_{max} = 1.1, second row), spiral wave dynamics are stable for the standard and intermediate sodium dynamics, but for the LR dynamics, spiral breakup occurs after 1.17 s. For parameter *setting 3* (*S*_{max} = 1.4, third row), spiral wave dynamics are only stable for the standard *I*_{Na} dynamics. For the intermediate dynamics, spiral breakup occurs after 2 s, while for the LR dynamics, breakup occurs after 0.72 s. For *setting 4* (*S*_{max} = 1.9, fourth row), spiral breakup occurs for all three *I*_{Na} dynamics. For the standard sodium dynamics, breakup occurs after 2.16 s; for the other two sodium dynamics, breakup occurs after 1.44 s.

The conditions we find for spiral breakup in two dimensions correspond nicely with the conditions we find for alternans in a ring. Using some extra parameter settings resulting in restitution slopes in between the 1.4 and 1.9 we used here, we found that the occurrence of alternans/spiral breakup under standard sodium recovery required a minimum APD restitution slope of 1.5 (data not shown).

### How Do I_{Na} Dynamics Affect Instability?

From the previous sections we conclude that for LR *I*_{Na} dynamics in our model a slope just over 1 (1.1) suffices to generate alternans and spiral breakup, whereas for the standard *I*_{Na} dynamics in our model, a slope of 1.5 is required to generate alternans and spiral breakup. This implies that slow *I*_{Na} recovery dynamics suppress alternans instability and breakup. Figure 8 shows the effect of *I*_{Na} dynamics on two factors that have been shown to be important for alternans instability: CV restitution and the DIs visited during spiral wave rotation.

Figure 8*A* shows three CV restitution curves generated with our model for the three different *I*_{Na} dynamics. We can see that the standard *I*_{Na} dynamics results in a CV restitution curve that starts gradually declining for DIs of 250 ms and smaller. The LR *I*_{Na} dynamics result in a CV restitution curve that is nearly flat for a broad range of DIs and then declines fast for DIs <65 ms. We can see that the intermediate *I*_{Na} dynamics indeed result in an intermediate CV restitution curve. It has been shown in both analytical and modeling studies that a shallow CV restitution curve suppresses alternans and spiral breakup (6, 7, 10): in tissue, through electrotonic interactions, gradual CV restitution increases the threshold value for APD restitution slopes above which instability occurs above 1.

Figure 8*B* shows periods and DIs visited during spiral wave rotation for the three different *I*_{Na} dynamics. For illustration purposes, we did this for parameter *setting 1* (restitution slope <1) to ensure stable spiral wave rotation for all sodium dynamics. We can see that when going from standard to intermediate to LR sodium dynamics, both period and DI decrease. It is a well-known effect that slowing of spiral wave rotation leads to stabilization. We found in our earlier studies that the mechanism behind this stabilization is that for slower sodium recovery and hence longer periods and DIs, a less steep part of the restitution curve is visited, thus suppressing alternans and spiral breakup (44, 62, 63).

Is it the shallower CV restitution or the longer DI during spiral wave rotation, or both, that leads to suppression of instability? We tried to resolve this issue by changing DI while keeping the shape of the CV restitution curve the same. This can be achieved by changing maximal conductance (*G*_{Na}) of the sodium current. Table 3 shows period and DIs of spiral wave rotation for parameter *setting 1*, either our standard or LR *I*_{Na} dynamics and different *G*_{Na}. From Table 3 we can see that DI can be reduced from 32 ms, occurring for standard *I*_{Na} dynamics, to 22 ms, close to the DI for LR *I*_{Na} dynamics, by increasing *G*_{Na} by a factor of 7. Similarly, we can increase DI from 19 ms, occurring for the LR *I*_{Na} dynamics, to 29 ms, close to the DI for standard *I*_{Na} dynamics, by multiplying *G*_{Na} with a factor of 0.4. Note that although the CV restitution curve will shift as a result of these modifications, its shape will be maintained (not shown). Furthermore, we assume that multiplication of *G*_{Na} by these factors will also work for the other parameter settings of our model.

We simulated spiral wave dynamics for parameter *settings 2* and *3* (slope 1.1 and 1.4, respectively) either using our standard *I*_{Na} dynamics combined with *G*_{Na} × 7 or using LR *I*_{Na} dynamics combined with *G*_{Na} × 0.4. We see (Fig. 9*B*) that for *setting 2* only increasing DI for LR *I*_{Na} dynamics suffices to remove spiral breakup; thus we obtain the same spiral wave dynamics as for our standard *I*_{Na} dynamics. Similarly, Fig. 9*C* shows that for *setting 3*, decreasing DI for standard *I*_{Na} dynamics induces spiral breakup, thus resulting in the same spiral wave dynamics as for LR *I*_{Na} dynamics. These two simulations thus suggest that changes in DI alone can explain the effects of the different *I*_{Na} dynamics on alternans and spiral wave stability observed in this study. However, our other simulation results do not support this conclusion. For *setting 2*, decreasing DI for standard *I*_{Na} dynamics does not result in the onset of spiral breakup, as observed for LR *I*_{Na} dynamics; however, it does result in transient initial wave breaks (Fig. 9*A*). For *setting 3*, increasing DI for LR *I*_{Na} does not prevent spiral breakup (Fig. 9*D*), as expected from computations for standard *I*_{Na} dynamics. Therefore, we conclude that although the DIs visited during spiral wave rotation have a pronounced effect on the onset of alternans and spiral breakup, CV restitution also plays an important role in determining whether instability will occur.

## DISCUSSION

### New Version of Human Ventricular Cell Model

In this study we developed a new version of our human ventricular cell model. In the new model we incorporated a subspace calcium variable that controls the dynamics of the *I*_{CaL} and the ryanodine receptor current. Similar approaches have been taken in many recent models (19, 20, 54). In addition, we changed the gating dynamics of *I*_{CaL}, which now has a fast subspace calcium inactivation gate *f*_{cass}, and a slow and fast voltage inactivation gate *f* and *f*_{2}, resulting in *I*_{CaL} dynamics that agree better with experimental data. We have also replaced the phenomenological description of CICR in our previous model with a reduced version of the ryanodine receptor Markov model developed by Stern et al. (57) and Shannon et al. (54). We demonstrated that with this new description our model is able to reproduce a nonlinear gain function for calcium release as a function of SR calcium load.

Although we improved the description of calcium handling in this new version of our model, our formulation is still a simplification and cannot describe all details of calcium dynamics. For example, to describe gradedness of CICR, much more complex models incorporating the description of individual L-type channels, ryanodine receptors, etc., are necessary (15, 51). However, such models are computationally very demanding and therefore not feasible for modeling cardiac tissue sheets or complete ventricles, for which our model is intended.

### Conditions for Alternans and Spiral Breakup

We use our new model to study the conditions under which electrical instability and spiral breakup occur to establish whether the spiral breakup hypothesis of VF can be valid for human ventricular tissue. We focus both on the role of APD restitution, reproducing a representative range of recently measured human ventricular restitution curves from Nash et al. (39), and on the role of the recovery dynamics of the fast *I*_{Na}, of which we found in our previous study that it is much slower than was previously assumed in a lot of cardiac tissue modeling studies (61).

We find that for fast (LR) *I*_{Na} recovery dynamics, an APD restitution slope just over 1 (1.1) is enough to result in alternans in single cells and rings of cells and to generate spiral breakup in two-dimensional sheets of tissue. For slower (standard and intermediate) *I*_{Na} recovery dynamics, steeper APD slopes are required to get electrical instability in rings and sheets of tissue.

These findings imply that slow *I*_{Na} recovery dynamics have a protective effect against steep APD restitution-mediated instability and spiral breakup. This may be an important finding, given that alternans and breakup are most often studied in models with fast (LR) *I*_{Na} dynamics. As we demonstrated earlier (61), voltage-clamp data for human cardiac fast sodium channels show slow *I*_{Na} recovery dynamics. Here we find that for the standard *I*_{Na} dynamics of our model, which are based on these experimental data, an APD restitution slope significantly steeper than 1, for our particular settings 1.5, is needed to generate spiral breakup. This 1.5 is well within the range of experimentally found values (38–40, 43). Therefore, our results do support that steep APD restitution-mediated instability is a potential mechanism for VF in the human heart.

### Mechanism of Breakup Suppression

Slow *I*_{Na} dynamics lead both to longer periods and DIs of spiral wave rotation (Fig. 8*B*) and a more gradual CV restitution curve (Fig. 8*A*). We performed a qualitative study of the relative importance of DI during spiral wave rotation and more gradual CV restitution for the suppression of instability, by isolating the CV effect and compensating for the DI effect of the different *I*_{Na} dynamics using conductance. We found that in some cases the DI effect can explain differences in spiral breakup or no spiral breakup for the different *I*_{Na} dynamics. However, for other parameter settings, differences in spiral wave stability could not be explained merely by the DI effect. This suggests that spiral wave period is not the only important effect of slower *I*_{Na} recovery that leads to alternans suppression and that the CV effect also plays a role in suppressing instability. Unfortunately, we found no way to change CV restitution shape while keeping the DI of spiral wave rotation constant to further establish the precise importance of CV restitution vs. DI.

### Relation to Work of Others

In analytical studies, it was shown that for waves propagating in tissue, the condition for steep restitution-mediated alternans instability was not simply slope >1 but is slope > 1 + ξ*c*′/*c*^{2} (7, 10). Here ξ represents the strength of action potential morphology (especially repolarization wave back)-dependent electrotonic interactions, *c*′ is the derivative of the CV curve for the relevant DI (i.e., that of rotation on a ring or spiral wave rotation), and *c* is the maximum planar CV. From the formula it follows that the larger the value of ξ (stronger electrotonic interactions) and the larger the value of *c*′/*c*^{2} (large change of CV over range of DIs, or low CV), the steeper is the slope required to generate instability. Here we indeed found that the slower the *I*_{Na} recovery dynamics, and hence the more gradual the CV restitution curve, the steeper was the APD restitution slope needed to generate instability. Our results agree with a previous study from Qu et al. (48), who showed that a CV restitution curve that stays flat over a broad range of DIs, similar to our steep CV curve, promotes spiral breakup, and to results from Fox et al. (13), who showed that chances of conduction block decrease for lower CVs at short cycle lengths, corresponding to our more gradual CV restitution. Our results also agree with results from an extensive simulation study on the dependence of alternans and spiral breakup on APD restitution, CV restitution, electrotonic interactions, and cardiac memory performed by Cherry and Fenton (6). They found that for action potentials with a fast repolarizing wave back, leading to strong electrotonic interactions, a more gradual CV restitution can suppress alternans in a ring and spiral breakup in two dimensions. Note that in our human ventricular cell model, we have a strong repolarizing wave back caused by the large inward rectifier K^{+} current (*I*_{K1}).

Our results may seem dissimilar from results such as those of Cao et al. (5) and Watanabe et al. (65), who showed that a more gradual CV restitution is proarrhythmic because it promotes the spatial discordance of alternans, causing spiral breakup to require a smaller tissue size. However, as noted by Cherry and Fenton (6), gradual CV restitution has a dual effect: on the one hand it can suppress alternans provided that electrotonic interactions are strong, and on the other hand it promotes the spatial discordance of any alternans still present. Whether alternans arises thus depends on the strength of electrotonic interactions, the shape of CV restitution, and tissue size. In our model, we found in the two-dimensional sheets of tissue of 20 × 20 cm only the alternans-suppressing effect. This suggests that a larger tissue size is needed for discordant alternans to develop. Given that 20 × 20 cm is large relative to human heart size, we conclude that for the human ventricles and standard *I*_{Na} dynamics, the alternans-suppressing effect will dominate.

The results we found here also agree with previous work from our group (44, 62, 63). In these studies we found that the presence of inexcitable obstacles or the removal of gap junctional connections between cells leads to suppression of spiral breakup. This suppression was caused by the slowing down of spiral wave rotation, leading to longer DIs, similar to the effect of slower fast *I*_{Na} recovery. Longer DIs shift the position of the spiral wave on the (unchanged) restitution curve to the right and upward, away from the steepest part.

### Limitations

The ionic mechanisms responsible for the variation in APD restitution curves are currently unknown. In this study, to obtain different restitution slopes, we used more general knowledge of the influence of certain parameters on APD and APD restitution. We varied the dynamics of *I*_{CaL}, as it has been shown in both modeling studies (12, 49) and experiments (50) that *I*_{CaL} is a major determinant of restitution slope. APD was then readjusted by varying conductances of *I*_{Kr}, *I*_{Ks}, *I*_{pCa}, and *I*_{pK} (Table 2). If new data on the ionic current differences underlying restitution variation become available, new parameter settings based on these data need to be constructed. However, we expect that restitution slope rather than the parameter setting is important for the observed effects, and hence that results will stay similar.

In this study we find that a slow *I*_{Na} recovery dynamics protect against steep APD restitution-mediated alternans. However, we investigated this in homogeneous rings and sheets of cardiac tissue, whereas real cardiac tissue is heterogeneous in numerous properties, including APD and APD restitution (39, 43), is anisotropic, and has a complex anatomy. What the precise conditions for spiral breakup are under more realistic conditions therefore remains to be investigated. In addition, we did not investigate the influence of the amount of electrotonic coupling and cardiac memory on the conditions for alternans and spiral breakup.

Furthermore, the steep APD restitution we studied here is not the only mechanism that can lead to alternans instability. Other instabilities, for example in the intracellular calcium dynamics (47, 55, 56), can lead to action potential alternans. Finally, other mechanisms than alternans leading to spiral breakup have been suggested to underlie VF, for example, the mother rotor hypothesis.

In conclusion, we propose a novel model for human ventricular cells. We report that *I*_{Na} recovery dynamics play an important role in the occurrence of steep restitution-mediated dynamical electrical instability. The slow recovery dynamics found in experiments on human cardiac *I*_{Na} suppress electrical instability. We conclude that steep restitution-mediated fibrillation can occur in human ventricular tissue. However, an APD restitution slope considerably steeper than 1 is required.

## GRANTS

This work was supported by the Netherlands Organization for Scientific Research through Grant 635100004 of the Research Council for Physical Sciences (K. H. W. J. ten Tusscher).

## APPENDIX

No changes were made to formulations of the following currents: *I*_{Na}, *I*_{to}, *I*_{Kr}, *I*_{K1}, *I*_{NaCa}, *I*_{NaK}, *I*_{pCa}, *I*_{pK}, *I*_{bNa}, and *I*_{bCa}. For these formulations, we refer to their description in the previous version of our model (61).

### L-Type Ca^{2+} Current

(6) (7) (8) (9) (10) (11) (12) (13) (14) (15) (16) (17) (18) (19) (20) (21) (22) (23)

### Slow Delayed Rectifier Current

(24) (25) (26) (27) (28)

### Calcium Dynamics

In the following equations, Ca_{itotal} is total (free + buffered) cytoplasmic Ca^{2+} concentration; Ca_{SRtotal} is total SR Ca^{2+} concentration; Ca_{SStotal} is total diadic subspace Ca^{2+} concentration; Ca_{i} is free cytoplasmic Ca^{2+} concentration; Ca_{SR} is free SR Ca^{2+} concentration; Ca_{SS} is free diadic subspace Ca^{2+} concentration; *I*_{rel} is CICR current; *I*_{up} is SR Ca^{2+} pump current; *I*_{leak} is SR Ca^{2+} leak current; *I*_{xfer} is diffusive Ca^{2+} current between diadic Ca^{2+} subspace and bulk cytoplasm; O is proportion of open *I*_{rel} channels; and R̄ is proportion of closed *I*_{rel} channels (for parameters, see Table 1). (29) (30) (31) (32) (33) (34) (35) (36) (37) (38) (39) (40) (41) (42) (43)

## Acknowledgments

We thank Dr. J. Sneyd, Dr. M. Nash, and Dr. P. Taggart for valuable discussions.

## Footnotes

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.

- Copyright © 2006 by the American Physiological Society