HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS		QUICK SEARCH:   [advanced] Author: Keyword(s): Year:  Vol:  Page:

This Article Abstract Full Text (PDF) All Versions of this Article: 291/3/H1088    most recent 00109.2006v1 Alert me when this article is cited Alert me if a correction is posted Citation Map Services Email this article to a friend Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via ISI Web of Science (15) Citing Articles via Google Scholar Google Scholar Articles by ten Tusscher, K. H. W. J. Articles by Panfilov, A. V. Search for Related Content PubMed PubMed Citation Articles by ten Tusscher, K. H. W. J. Articles by Panfilov, A. V.

Alternans and spiral breakup in a human ventricular tissue model

Department of Theoretical Biology, Utrecht University, Utrecht, The Netherlands

Submitted 30 January 2006 ; accepted in final form 5 March 2006

	ABSTRACT

TOP ABSTRACT MATERIALS AND METHODS RESULTS DISCUSSION GRANTS APPENDIX REFERENCES

 
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

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

TOP ABSTRACT MATERIALS AND METHODS RESULTS DISCUSSION GRANTS APPENDIX REFERENCES

 
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.

View larger version (9K): [in this window] [in a new window]   Fig. 1. Calcium compartments involved in excitation-contraction coupling. L-type calcium channels (ICaL) release calcium in the diadic space or subspace (SS), where sarcolemmal membrane and membrane of the sarcoplasmic reticulum (SR) are in close proximity. Ryanodine receptors (Irel) sense this elevation of calcium in the subspace and respond with a release of calcium from the SR. Through diffusion (Ixfer), the calcium released in the subspace travels to the cytoplasm (cyto). Sodium-calcium exchangers (INaCa) and the sarcolemmal calcium pump (IpCa) pump calcium out of the cytoplasm to the exterior of the cell. Calcium pumps in the membrane of the SR (Iup) pump calcium out of the cytoplasm back into the SR. A small leakage current (Ileak) leaks calcium from the SR to the cytoplasm.

View larger version (7K): [in this window] [in a new window]   Fig. 2. A: steady-state epicardial action potential (V) shape at 1-Hz pacing. B: steady-state cytoplasmic calcium (Cai) transient at 1-Hz pacing. C: steady-state subspace calcium (CaSS) transient at 1-Hz pacing.

View larger version (12K): [in this window] [in a new window]   Fig. 3. A: steady-state subspace calcium inactivation (fcass) curve. B: time constant curve of subspace calcium inactivation (fcass). C: steady-state curves of slow (f) and fast (f2) voltage inactivation. D: time constant curves of slow (f) and fast (f2) voltage inactivation. Experimental (exp) time constants are added for comparison and are taken from Beuckelmann et al. (5), Benitah et al. (3), Mewes et al. (36), Sun et al. (58), Li et al. (32), Pelzmann et al. (46), and Magyar et al. (35). Experimental inactivation time constants represent the slowest time scale of calcium current inactivation, thus corresponding to f gate inactivation. For recovery from inactivation, both fast and slow experimental time constants are available, corresponding to f2 and f gate recovery. E: peak current-voltage relationship of the L-type calcium current as determined in simulated voltage-clamp experiments. Experimental data are from Magyar et al. (34).

View larger version (9K): [in this window] [in a new window]   Fig. 4. A: schematic drawing of the 4-state ryanodine receptor Markov model of Shannon et al. (54) and Stern et al. (57). O, open conducting state; R, resting closed state; I, inactivated closed state; RI, resting inactivated closed state. k3 and k4 are constants; the values of k1 and k2 are dependent on SR calcium concentration: opening rate is sped up and closing rate slowed down by higher SR calcium load. In addition, k1 is multiplied by the square of the subspace calcium concentration, and k2 is multiplied by the subspace calcium concentration to give a final rate constant: opening and closing rate increase with triggering subspace calcium levels, but opening more so than closing. B: nonlinear gain function of calcium-induced calcium release (CICR) as a function of free calcium concentration in the sarcoplasmic reticulum (CaSR). Experimental data from Shannon et al. (53) for rabbit ventricular myocytes are added for comparison.

(1)

(2)

(3)

where = 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 4B 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²/ms; in two dimensions, D_ij = 0.00154 cm²/ms for i = j, and D_ij = 0 cm²/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₂, and f_cass for I_CaL).

We measure APD at either 50% or 90% repolarization and will refer to these values as APD₅₀ and APD₉₀, respectively. APD₉₀ values will be used, unless stated otherwise. APD₅₀ 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 dV/dt during the QRS wave to the maximum dV/dt 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₅₀ 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 x 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

TOP ABSTRACT MATERIALS AND METHODS RESULTS DISCUSSION GRANTS APPENDIX REFERENCES

 
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₅₀ (see METHODS). Model restitution curves are plotted for both APD₅₀ and APD₉₀ 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. 5A); for setting 2, it is 1.1 (Fig. 5B); for setting 3, it is 1.4 (Fig. 5C); and for setting 4, it is 1.9 (Fig. 5D). Table 1 lists the complete default parameter setting, and Table 2 lists the parameter values varied to obtain these different restitution slopes.

View larger version (19K): [in this window] [in a new window]   Fig. 5. Single-cell action potential duration (APD) restitution. A–D: restitution curves obtained using a S1-S2 restitution protocol for a basic cycle length (BCL) of 600 ms measuring APD at 50% (APD50) and 90% repolarization (APD90) for 4 different parameter settings of our model (Table 2). APD is plotted against diastolic interval (DI). Experimental activation recovery interval (ARI) restitution curves (exp) are from Nash et al. (38–4040). E–H: restitution curves obtained using a dynamic restitution protocol measuring APD90 for the 4 same parameter settings. APD is plotted against DI. Gray shading is used to indicate the region of the restitution curve where the slope exceeds 1. I–L: same dynamic restitution curves as in E–H, but now plotting APD against stimulation period. Splitting of the restitution curve indicates the presence of 2 APDs for a single period: APD alternans.

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₉₀. Maximum restitution slopes for our four setting are now 0.8 (Fig. 5E), 1.1 (Fig. 5F), 1.4 (Fig. 5G), and 1.8 (Fig. 5H), 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. 5I), no alternans occurs. For setting 2 (Fig. 5J), 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. 5F 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. 5K), APD alternans occurs for periods <280 ms, corresponding to DIs < 74 ms, for which restitution slope in Fig. 5G is >1, and reaches a maximum amplitude of 100 ms. For setting 4 (Fig. 5L), alternans occurs for BCL of 320 ms and smaller, corresponding to DIs < 94 ms, for which slope in Fig. 5H 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 x _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 6A 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. 6B 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.

View larger version (6K): [in this window] [in a new window]   Fig. 6. APD restitution in a ring for the 12 different model settings. Restitution curves are obtained by initiating a pulse propagating in one direction that propagates until steady-state APD is reached, after which the ring length and hence BCL is shortened. A: restitution for the model setting with a maximum slope of 0.7 combined with the standard, intermediate, and Luo-Rudy (33) model (LR) fast INa dynamics. B: restitution for the model setting with a maximum slope of 1.1 combined with the standard, intermediate, and LR fast INa dynamics. C: restitution for the model setting with a maximum slope of 1.4 combined with the standard, intermediate, and LR fast INa dynamics. D: restitution for the model setting with a maximum slope of 1.9 combined with the standard, intermediate, and LR fast INa 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.

View larger version (99K): [in this window] [in a new window]   Fig. 7. Snapshots of spiral wave dynamics for the 12 different model settings. Columns (left to right): column 1, standard INa dynamics; column 2, intermediate INa dynamics; column 3, LR INa dynamics. Rows (top to bottom): row 1, slope 0.7 setting; row 2, slope 1.1 setting; row 3, slope 1.4 setting; row 4, slope 1.9 setting. If spiral breakup occurs, the time of first wave break is indicated at top right corner.

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.

View larger version (12K): [in this window] [in a new window]   Fig. 8. A: conduction velocity (CV) restitution for the shallow parameter setting of our model (slope 0.7) combined with the standard, intermediate, and LR fast INa dynamics. Restitution was measured in a cable of 800 cells. B: period (3 top lines) and DI (3 bottom lines) as a function of time, for spiral wave rotation for the shallow parameter setting of our model combined with the standard, intermediate, and LR fast INa dynamics.

Figure 8B 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.

View larger version (29K): [in this window] [in a new window]   Fig. 9. Snapshots of spiral wave dynamics for 4 different model settings. A: slope 1.1, standard INa dynamics GNa x 7. B: slope 1.1, LR INa dynamics GNa x 0.4. C: slope 1.4, standard INa dynamics GNa x 7. D: slope 1.4, LR INa dynamics GNa x 0.4.

	DISCUSSION

TOP ABSTRACT MATERIALS AND METHODS RESULTS DISCUSSION GRANTS APPENDIX REFERENCES

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₂, 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. 8B) and a more gradual CV restitution curve (Fig. 8A). 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² (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² (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 x 20 cm only the alternans-suppressing effect. This suggests that a larger tissue size is needed for discordant alternans to develop. Given that 20 x 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

TOP ABSTRACT MATERIALS AND METHODS RESULTS DISCUSSION GRANTS APPENDIX REFERENCES

 
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

TOP ABSTRACT MATERIALS AND METHODS RESULTS DISCUSSION GRANTS APPENDIX REFERENCES

 
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²⁺ 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²⁺ concentration; Ca_SRtotal is total SR Ca²⁺ concentration; Ca_SStotal is total diadic subspace Ca²⁺ concentration; Ca_i is free cytoplasmic Ca²⁺ concentration; Ca_SR is free SR Ca²⁺ concentration; Ca_SS is free diadic subspace Ca²⁺ concentration; I_rel is CICR current; I_up is SR Ca²⁺ pump current; I_leak is SR Ca²⁺ leak current; I_xfer is diffusive Ca²⁺ current between diadic Ca²⁺ subspace and bulk cytoplasm; O is proportion of open I_rel channels; and 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

 

Address for reprint requests and other correspondence: K. H. W. J. ten Tusscher, Utrecht Univ., Dept. of Theoretical Biology, Padualaan 8, 3584 CH Utrecht, The Netherlands (e-mail: khwjtuss{at}hotmail.com)

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.

	REFERENCES

TOP ABSTRACT MATERIALS AND METHODS RESULTS DISCUSSION GRANTS APPENDIX REFERENCES

 

1. Banville I, Chattipakorn N, and Gray RA. Restitution dynamics during pacing and arrhythmias in isolated pig hearts. J Cardiovasc Electrophysiol 15: 455–463, 2004.[CrossRef][ISI][Medline]
2. Banville I and Gray RA. Effect of action potential duration and conduction velocity restitution and their spatial dispersion on alternans and the stability of arrhythmias. J Cardiovasc Electrophysiol 13: 1141–1149, 2002.[CrossRef][ISI][Medline]
3. Benitah J, Bailly P, D'Agrosa M, Da Ponte J, Delgado C, and Lorente P. Slow inward current in single cells isolated from adult human ventricles. Pflügers Arch 421: 176–187, 1992.[CrossRef][ISI][Medline]
4. Beuckelmann DJ, Nabauer M, and Erdmann E. Characteristics of calcium-current in isolated human ventricular myocytes from patients with terminal heart failure. J Mol Cell Cardiol 23: 929–937, 1991.[CrossRef][ISI][Medline]
5. Cao J, Qu Z, Kim Y, Wu T, Garfinkel A, Weiss JN, Karagueuzian HS, and Chen P. Spatiotemporal heterogeneity in the induction of ventricular fibrillation by rapid pacing, importance of cardiac restitution properties. Circ Res 84: 1318–1331, 1999.[Abstract/Free Full Text]
6. Cherry EM and Fenton FH. Suppression of alternans and conduction blocks despite steep APD restitution: electrotonic, memory and conduction velocity effects. Am J Physiol Heart Circ Physiol 286: H2332–H2341, 2004.[Abstract/Free Full Text]
7. Cytrynbaum E and Keener JP. Stability conditions for the traveling pulse: modifying the restitution hypothesis. Chaos 12: 788–799, 2002.[CrossRef][ISI][Medline]
8. Drouin E, Charpentier F, Gauthier C, Laurent K, and Le Marec H. Electrophysiologic characteristics of cells spanning the left ventricular wall of human heart: evidence for the presence of M cells. J Am Coll Cardiol 26: 185–192, 1995.[Abstract]
9. Drouin E, Lande G, and Charpentier F. Amiodarone reduces transmural heterogeneity of repolarization in the human heart. J Am Coll Cardiol 32: 1063–1067, 1998.[Abstract/Free Full Text]
10. Echebarria B and Karma A. Instability and spatiotemporal dynamics of alternans in paced cardiac tissue. Phys Rev Lett 88: 208101, 2002.[CrossRef][Medline]
11. Fox JJ, Bodenschatz E, and Gilmour RF Jr. Period-doubling instability and memory in cardiac tissue. Phys Rev Lett 89: 138101, 2002.[CrossRef][Medline]
12. Foz JJ, McHarg JL, and Gilmour RF. Ionic mechanism of electrical alternans. Am J Physiol Heart Circ Physiol 282: H516–H530, 2002.[Abstract/Free Full Text]
13. Fox JJ, Riccio ML, Drury P, Werthman A, and Gilmour RF Jr. Dynamic mechanism for conduction block in heart tissue. New J Phys 5: 101.1–101.14, 2003.
14. Gilmour RF, Heger JJ, Prystowsky EN, and Zipes DZ. Cellular electrophysiologic abnormalities of diseased human ventricular myocardium. Am J Cardiol 51: 137–144, 1983.[CrossRef][ISI][Medline]
15. Greenstein JL and Winslow RL. An integrative model of the cardiac ventricular myocyte incorporating local control of Ca²⁺ release. Biophys J 83: 2918–2945, 2002.[ISI][Medline]
16. Guevara MR, Ward A, Shrier A, and Glass L. Electrical alternans and period doubling bifurcations. IEEE Comp Cardiol 562: 167–170, 1984.