## Abstract

The onset of ventricular fibrillation (VF) has been associated with steep action potential duration restitution in both clinical and computational studies. Recently, detailed clinical restitution properties in cardiac patients were reported showing a substantial degree of heterogeneity in restitution slopes at the epicardium of the ventricles. The aim of the present study was to investigate the effect of heterogeneous restitution properties in a three-dimensional model of the ventricles using these clinically measured restitution data. We used a realistic model of the human ventricles, including detailed descriptions of cell electrophysiology, ventricular anatomy, and fiber direction anisotropy. We extended this model by mapping the clinically observed epicardial restitution data to our anatomic representation using a diffusion-based algorithm. Restitution properties were then fitted by regionally varying parameters of the electrophysiological model. We studied the effects of restitution heterogeneity on the organization of VF by analyzing filaments and the distributions of excitation periods. We found that the number of filaments and the excitation periods were both dependent on the extent of heterogeneity. An increased level of heterogeneity leads to a greater number of filaments and a broader distribution of excitation periods, thereby increasing the complexity and dynamics of VF. Restitution heterogeneity may play an important role in providing a substrate for cardiac arrhythmias.

- action potential duration-restitution-heterogeneity
- reentrant arrhythmias
- ventricular fibrillation
- human ventricular myocytes
- computer simulation

sudden cardiac death is the most common cause of death in the industrialized world, and in most cases it is due to ventricular fibrillation (VF) (59). During VF, excitation waves are disturbed and chaotic, causing the contraction of the ventricles to become rapid and uncoordinated. It has been shown in clinical and experimental studies that these turbulent wave patterns are underpinned by reentrant sources of excitation (12, 15, 19, 28, 51, 54). If VF is not halted by means of defibrillation, this condition can be lethal within minutes.

Historically, it has been hypothesized that VF is due to turbulent wave propagation induced by heterogeneity in refractoriness of cardiac tissue (23, 27). The important role of heterogeneity in VF is supported by the strong association of VF with diseases that lead to increased heterogeneity, such as coronary artery disease, cardiomyopathies, and congenital heart disease (11, 59). Because of these pathologies, remodeling of the cellular ionic currents, calcium handling, and gap junctional coupling in cardiac tissue results in increased heterogeneity. Furthermore, VF is associated with diseases such as Long QT, Short QT, and Brugada syndrome (11, 59), in which ion channel mutations lead to an increase of action potential duration (APD) dispersion.

In the last decade, it has been proposed that VF is caused by the occurrence of dynamic instabilities (alternans) induced by steep APD restitution (16, 33, 35, 36, 41, 52), which became known as the restitution hypothesis of VF (41, 52). More recently, intracellular calcium cycling has been shown to cause alternans in APD (40, 53), which may lead to dynamic instabilities and VF. According to the original restitution hypothesis, alternans instability occurs when the slope of the APD restitution curve, which relates the APD to its previous diastolic interval (DI), exceeds one (16, 33). Recent experimental and modeling studies show that, although steepness of restitution curve is an important determinant of instability (2, 3, 55), other factors such as electrotonic coupling, cardiac memory, and conduction velocity restitution should also be taken into account (13, 39, 49, 50). The restitution hypothesis and refractoriness hypothesis of VF are related. Alternans instability induces dynamic heterogeneities that can lead to wave break, whereas static anatomically defined heterogeneities can also lead to wave break.

Restitution-induced VF has been the subject of a variety of computational studies. Because the restitution hypothesis does not require preexisting heterogeneity for the onset of VF, most of the computational studies regarding the restitution hypothesis have assumed that cardiac tissue is homogeneous. However, recent clinical studies have demonstrated that APD restitution properties on the ventricular epicardium and endocardium have a complex spatially nonuniform distribution (26, 31, 58). These studies showed that: *1*) APD restitution properties are spatially heterogeneous, containing multiple sites of steep and shallow restitution slopes; *2*) maximum restitution slopes range from almost zero to well over two; and *3*) restitution properties differ between patients and are dependent on the underlying pathology. Some studies on the effects of heterogeneous APD restitution in cardiac tissue have shown that spatial dispersion in APD restitution may result in wave break and initiation of VF (6, 9, 57). However, the influence of restitution heterogeneity on VF dynamics and organization has not yet been investigated.

In this paper, we study the effects of heterogeneous APD restitution data measured in clinical studies (31) on VF organization in the human heart. We use our integrative anatomic model of the human ventricles (48), which was recently developed and verified against available clinical data on human VF. The model describes action potential characteristics of human ventricular myocytes (49, 50), combined with an anatomically detailed geometry of the human ventricles that incorporates fiber direction anisotropy (18). We introduce APD restitution heterogeneity into our model by incorporating clinically measured restitution data from patients who were undergoing surgery for aortic valve replacement as reported previously (31). Clinical restitution curves were mapped and interpolated to our three-dimensional (3D) ventricular geometry and reproduced by our cell model using appropriate regional variations of the parameters of the electrophysiological model. We used this human ventricular model to study the effects of spatial heterogeneity of APD restitution on the dynamics of VF and found that the organization of VF became more complex with increased heterogeneity. We have quantified these effects in terms of the number of filaments and the local excitation periods.

## METHODS

We developed a heterogeneous model of the human ventricles by incorporating clinically measured spatially heterogeneous APD restitution properties (31) into our integrative model of the human ventricles (48, 49, 50). Below we describe in detail the main features of our heterogeneous model.

#### Ventricular model.

Excitation of the ventricles in a monodomain model can be described using the following partial differential equation (22): (1) where *V*_{m} denotes the transmembrane voltage, *C*_{m} the membrane capacitance, *D*_{ij} the diffusion tensor, *I*_{ion} the sum of the ionic transmembrane currents describing the excitable behavior of the individual ventricular cell; and *x*_{i} and *x*_{j} the tensor notation of the *x*, *y*, and *z* coordinates. There are different models available to describe action potential characteristics of the cell membrane (1, 4, 13, 24, 25, 50). In this study, we used the human ventricular ionic model published by ten Tusscher et al. (49, 50) and refer to this as the TNNP model. This model provides a detailed description of voltage, ionic currents, and intracellular ion concentrations and is based on a wide range of human-based electrophysiological data. The most recent version of the TNNP model (49) includes a more extensive description of intracellular calcium dynamics. However, the model does not describe calcium alternans nor spontaneous calcium release. The model reproduces experimentally measured APD (29) and conduction velocity restitution curves (14). The TNNP model contains the following ionic currents: (2) where *I*_{Na}, is the fast sodium current, *I*_{to} is the transient outward current, *I*_{CaL} is the L-type calcium current, *I*_{Kr} is the rapid delayed rectifier current, *I*_{Ks} is the slow delayed rectifier current, *I*_{K1} is the inward rectifier potassium current, *I*_{NaCa} is the sodium/calcium exchanger current, *I*_{NaK} is the sodium/potassium pump current, *I*_{pCa} is the plateau calcium current, I_{pK} is the plateau potassium current, *I*_{Na,b} is the sodium background current, and *I*_{Ca,b} is the calcium background current. A complete list of all equations can be found in Refs. 49 and 50. We used the “default” parameter settings from Ref. 49; all changes from these default settings are detailed in the text below. More information about homogeneous two-dimensional (2D) and 3D spiral wave dynamics and spiral wave breakup can be found in Refs. 48 and 49.

Geometric data describing the 3D ventricular anatomy and fiber direction field were derived from a normal healthy human heart (18). From this heart, a voxel description of ∼13.5 million points with a spatial resolution of 0.25 mm was obtained. We assume that the transverse conductivity is the same in all directions orthogonal to the direction of the muscle fiber axis and that local conductivity tensors *D*_{ij} can be derived from local muscle fiber directions using: (3) where *D*_{L} and *D*_{T} are the longitudinal and transverse conductivity, respectively, α is the muscle fiber direction, δ_{ij} is the Kronecker delta, and α_{i} and α_{j} are muscle direction fibers. For *D*_{L} we used 162 Ω·cm and for *D*_{T} we used 40.5 Ω·cm, which resulted in conduction velocities of 68 cm/s in the longitudinal direction and 32 cm/s in the transverse direction. Thus the resulting anisotropy ratio was ∼2:1, which is consistent with clinical measurements (20).

#### Clinical human restitution data.

A recent clinical study by Nash et al. (31) has provided global epicardial restitution properties in 14 cardiac patients. In this study, activation-recovery intervals (ARI) were recorded over the entire ventricular epicardial surface using an epicardial sock containing 256 unipolar contact electrodes (interelectrode spacing ∼10 mm). These ARI values can be compared with APD measured at ∼50% repolarization (APD_{50}) (17). They applied a standard S1–S2 protocol to determine restitution properties. A train of nine S1 stimuli at a basic cycle length (BCL) of 600 ms was applied, followed by an S2 stimulus at various intervals ranging from 600 to 240 ms. For each electrode, up to 24 pairs of values for the DI and subsequent APD_{50} were recorded. We fitted the clinical ARI data of Nash et al. (31) using a least-square monoexponential fit, since this was the same approach used by Nash et al. (31): (4) where the steady-state APD_{50(SS)}, *a*, and *b* are parameters of the fit. The maximum restitution slope *S*_{max} was determined for each electrode using (31): (5) where DI_{min} is the minimum measured nonrefractory DI. Note that this method for maximum restitution slope determination differs from the method used in ten Tusscher et al. (49) who used a piecewise linear approximation to determine the maximum restitution slope. We compared both methods and found that the piecewise calculation systematically underestimated the maximum restitution slope by ∼0.2–0.3 compared with the least-squares monoexponential fit method [Nash et al. (31)].

#### Mapping and extrapolating clinical restitution data.

In the present study, we used datasets of two patients undergoing replacement surgery for aortic valve disease, for which the most complete APD restitution recordings were available, with >85% of valid electrode recordings for each dataset. An electrode was defined to be valid if the data contained at least 18 pairs of DI and APD_{50} values and the monoexponential fit had at least a correlation coefficient of 0.6 or higher, with one meaning perfect statistical correlation and zero meaning there is no correlation at all. *S*_{max} with a slope <0.4 or >3 were manually evaluated and corrected if needed. *Dataset 1* (i.e., ⇓⇓⇓Fig. 4 in Ref. 31) had a mean ± SD *S*_{max} of 1.78 ± 1.04, based on *n* = 238 (93%) valid restitution curves. For *dataset 2*, the mean ± SD *S*_{max} was 0.97 ± 0.97 based on *n* = 227 (89%) valid restitution curves.

To map the clinical restitution data to our human ventricles, we first mapped the positions of the 256 sock electrodes to the surface of our ventricular model. Coordinates of both the epicardial sock and the anatomic model were normalized about the center of the apex, taking into account the anterior, posterior, left, and right sides of the ventricles. Electrode sites were then projected on the epicardial surface of our ventricular model by minimizing the distance to the epicardial surface.

For our study, we needed to assign restitution properties of cardiac tissue throughout the entire 3D ventricular mass. However, data on the full 3D organization of APD restitution in the human ventricles are presently not available. Observations have been reported only for surfaces of the heart, and measuring restitution properties on the epi- and endocardium require different clinical methods. To date, either epicardial or endocardial recordings (but not both) have been reported for individual hearts. Nash et al. (31) used an epicardial sock with electrodes to measure restitution data at 256 locations spread across most of the epicardium. Yue et al. (58) used a catheter balloon with electrodes to measure restitution data at 128 locations on the endocardium of the heart. The epi- and endocardial datasets each show regional organization of restitution slopes into three to five different areas. Slopes varied between regions, and the spatial arrangement of these regions varied considerably from patient to patient. Given these observations regarding restitution slope patterns on the epi- and endocardial surfaces, it is likely that this heterogeneity has a 3D structure. In the absence of experimental 3D data, we extrapolated the 2D surface patterns across the 3D mass. Rather than using different datasets from two different hearts (i.e., one set of data for the epicardial surface and another for the endocardial surface), we used only the epicardial data from Nash et al. (31). To extrapolate the data into 3D, we developed a diffusion-based algorithm to extend the epicardial profile smoothly across the ventricular walls (described below). We do not claim that our algorithm will reproduce the quantitatively correct 3D restitution profile for a particular patient. However, we believe that this is a reasonable first step until more detailed data on 3D organization of human restitution become available.

#### Diffusion domain algorithm.

Our diffusion algorithm contains the following three steps.

First, we extended the measured *S*_{max} values (maximum 256 electrode points on the epicardium) to all heart points by solving the diffusion equations for the maximum restitution slopes in the domain given by the heart anatomy. We used the clinically measured maximum restitution slopes as constraints, i.e., the variable *S*_{max} was kept constant at each electrode position during the diffusion process. We used a value for the nondimensional diffusion coefficient of *D* = 400, which provided the desired spatial smoothness of interpolation. The diffusion equations were solved until the average square of residual at the nonconstraint points decreased below 10^{−10}. As a result, we obtained a smooth spatial interpolation of *S*_{max} across all points in the 3D model except for the regions immediately adjacent to the initial electrodes.

The second step was to remove these sharp local gradients by applying a moving average window of 3 × 3 × 3 voxels in the domain given by the heart anatomy. Note that this averaging also changed the values of *S*_{max} at the reference electrodes (see Fig. 1, *column 2*). As a third step, we scaled back *S*_{max} values for all points using a linear transformation function, which was found from the data presented in Fig. 2. More specifically, we found the best least-square linear fit of the original, clinically measured maximum restitution slope values and the *S*_{max} values found after *step 2* of our interpolation procedure. We then scaled back our data using this linear transformation. Thus we obtained a smooth distribution of *S*_{max} across all heart points in our anatomic model, which corresponded well with the clinical data at the epicardial surface.

Figure 1 shows the application of our algorithm for *dataset 1* (row on *top*) and *dataset 2* (row on *bottom*). In the first column, we show a direct linear surface interpolation of the clinical data using Matlab's GRIDDATA function. In the second column, we show how *S*_{max} values are distributed after the diffusion algorithm and the 3 × 3 × 3 moving average window were applied. The third column shows the maximum restitution slopes after the linear transformation step. Comparing the raw dataset shown in Fig. 1, *column 1*, with the interpolated data shown in *column 3*, we see that we were able to obtain a profile of heterogeneity in all heart points that corresponded well with the clinical data, but was spatially smooth. We chose to use this smoothed representation, since it provided a robust general description of heterogeneity in the heart that helped guard against possible spurious recordings. (Note that the *patient 1* dataset we use here is the same as the one used in Fig. 4 in Ref. 31).

For *dataset 1*, a large region of the posterior wall of the left ventricle contained steep restitution slopes, with values close to three. The right ventricle contains more shallow restitution slopes, with values around 0.5–1. Between these regions, intermediate values were present. The mean ± SD restitution slope was 1.78 ± 0.65. *Dataset 2* had a relatively small region on the anterior wall with high slope values, whereas the rest of the bulk tissue was characterized by shallow restitution. The mean ± SD restitution slope was 1.01 ± 0.48. Both profiles are shown in Fig. 3, where red areas correspond to steep restitution slopes (*S*_{max} = 2.5–3), and blue areas correspond to shallow restitution slopes (*S*_{max} = 0.5–1). Figure 11 shows two short-axis slices with maximum restitution slopes that have been interpolated from the epicardium through to the endocardium (for *dataset 1*).

#### Parameters of the cell model.

To model cardiac tissue with the desired regional APD restitution properties, we fitted parameters of the TNNP model to different maximum restitution slopes. For each computational point of the ventricular geometry, the value of *S*_{max} was used to individually assign cell model parameters. To generate restitution curves, we paced a single cell model with a series of 10 S1 stimuli at a BCL of 600 ms, followed by a single S2 stimulus delivered at a specified DI after the last S1 action potential. The restitution curve was generated by decreasing the S1–S2 interval (and hence DI) and plotting the APD_{50} of the S2 action potential vs. the preceding DI value. We reproduced a variety of slopes between 0.4 and 3 by varying the parameters *G*_{Kr} (maximal conductance of the rapid delayed rectifier current), *G*_{Ks} (maximal conductance of the slow delayed rectifier current), *G*_{pCa} (maximal conductance of the plateau calcium current), *G*_{pK} (maximal conductance of the plateau potassium current), and *t*_{f}. (time constant for the inactivation kinetics in L-type calcium current). Figure 4 shows the values of these parameters used to obtain the given maximum restitution slopes. Note, that for slopes *S*_{max} < 0.4 we used *S*_{max} = 0.4, and for slopes *S*_{max} > 3 we used *S*_{max} = 3 (gray lines in Fig. 4). Figure 5 shows an example of some APD and APD restitution curves for different parameter settings, together with experimental data to which these restitution curves were fitted. More information about fitting of experimental restitution curves using the TNNP model can be found in ten Tusscher et al. (49).

#### Numerical approach.

Equations for the gating variables in the TNNP model were integrated using a Rush and Larsen (45) integration scheme. To integrate *Eq. 1*, we used a forward Euler scheme with a time step of *t* = 0.02 ms and a space step of *x* = 0.25 mm, and the following Laplacian was evaluated at each point in the human ventricular geometry: (6) which can be discretized into the following equation: (7) where L is an index running over the 18 neighbors of the point (*i*, *j*, *k*) and the point itself, and *w*_{l} are the weights defined for each neighbor point containing information about the voltage contributed to the Laplacian at each point (*i*, *j*, *k*). Weights are calculated at each point using the local conductivity tensors calculated by *Eq. 3*. To ensure a zero axial current flow from heart points to nonheart points, no flux boundary conditions were imposed by setting weights to zero if a neighbor lay outside the heart geometry.

To initiate 3D scroll waves, we used an S1–S2 protocol in which the S2 stimulus was activated in the refractory tail of the S1 stimulus, thereby creating a scroll wave. Stimulus currents were applied at two times the diastolic threshold value. For each restitution slope profile, we performed simulations of VT/VF by initiating a scroll wave in the free wall of the left ventricle. Simulations lasted for 8 s or were terminated if there was no activity present.

Assuming that the medium is an infinite volume conductor, electrograms can be calculated using the dipole source density of the membrane potential in all voxel points using (37): (8) where *D* is the diffusion tensor, ▿ is the Nabla operator (in this case the gradient), V is the domain of integration (i.e., the ventricular volume), and **r** is the vector from each point to the recording electrode, which was placed 10 cm from the center of the ventricles in the anterior direction of the transverse plane.

Scroll wave filaments were detected using an algorithm proposed by Fenton and Karma (13). If there is a spiral wave in 2D (a scroll wave in 3D) present, the core (filament) can be defined as the point(s) for which the excitation wavefront and waveback meet. This can be calculated as the intersection points of the −60 mV isopotential line and the dV/d*t* = 0 isoline. Voxel data corresponding to these intersection points were then stored. Individual filaments were detected by iteratively joining neighboring voxels that were designated as filament points. All voxels that belonged to the same filament were assigned a unique identifier. Filaments were determined at 10-ms intervals. Of these filaments, we tracked: the time of birth, time of death, time of bifurcation, time of amalgamation, lifespan, filament from which a filament bifurcated, filament into which a filament amalgamated, and ultimate filament to which a filament can be traced back through the bifurcations events.

All simulations were coded in C++ and MPI and were run on 16 processors of a Beowulf cluster consisting of 16 Dell 650 Precision Workstations (dual Intel Xeon 2.66 GHz). Simulation of 1 s of wave propagation in the ventricles took ∼6.5 h of wall clock computation time. Ventricular geometry, wave patterns, and scroll wave filaments were visualized using the marching cubes algorithms for isosurface detection in voxel data and OpenGL for isosurface rendering.

## RESULTS

#### Wave fronts and ECG signal.

In the first series of simulations, we initiated a spiral in the free wall of the left ventricle using a S1–S2 protocol. In Fig. 6, we show snapshots of consecutive wave fronts and the ECG pattern for both patient datasets (*dataset 1* in the *top* row, and *dataset 2* in the *bottom* row). For *dataset 1*, we observed that the initial spiral broke down into a chaotic fibrillatory pattern and that the ECG signal was irregular and complex (dominant frequency of 4.6 Hz). For *dataset 2*, the initial spiral remained stable in the LV. Note that the spiral was not anchored and did not (hyper)meander nor drift during the entire simulation. The ECG signal was periodic (dominant frequency of 4.5 Hz), which resembles ECG signals during ventricular tachycardia (VT). The frequency spectra of both ECG signals are shown in Fig. 7. Although the simulations had similar dominant frequencies, the frequency spectrum of *dataset 1* was broader compared with that of *dataset 2*, which had a single dominant peak. These observations are characteristic of the differences between VF and VT, respectively, and the frequencies are similar to those reported in clinical (8, 30, 32, 56) and numerical (48) studies.

#### Number of filaments.

A convenient way to quantify the complexity of excitation patterns in the thick-walled ventricles of the heart is by determining the number of scroll wave filaments of the excitation sources (7, 10, 21, 34, 42, 44, 48, 51). A scroll wave filament is a line around which a 3D spiral wave rotates. When such a filament intersects with the surface of the heart, a phase singularity (PS) manifests on the heart surface. For *dataset 2*, we observed just one reentrant source, represented by a single filament.

In Fig. 8, we show a snapshot of the filaments present for *dataset 1* after 5 s of simulation time. Here the excitation pattern was organized by ∼15 filaments. Most of these filaments extended from the epicardium to the endocardium, with only five filaments located entirely within the myocardial walls. We observed that the majority of filaments were located near the base of the heart, whereas very few filaments were present near the apex. Furthermore, most filaments were present in the left ventricular wall, i.e., the region with high restitution slope values. The filament numbers and shapes fluctuated during the course of the simulation. In Fig. 8, we show the time dynamics of the number of filaments during 8 s of simulation. After an initial phase of growth, the number of filaments fluctuated around 12–20. The average number of filaments was 14.9, which is within the range reported by ten Tusscher et al. (48). We also show the number of PS visible on the epicardium. The average number of PS was 10.7; thus, the ratio of the number of filaments to PS was ∼1.4, consistent with (48).

We analyzed how filaments were created and destroyed over time through death, birth, bifurcation, and amalgamation events. This is shown in Fig. 9 for *dataset 1*. Horizontal lines indicate individual filaments and their lifespan. The start of a line is determined by the birth of a filament or the bifurcation from another filament. The end of the line occurs when the filaments disappear, i.e., the death or the amalgamation of a filament with another filament. Short horizontal lines in Fig. 9 correspond to short living filaments, and long lines correspond to long living filaments. During the simulation, we detected 1,463 filaments and 379 births, 547 deaths, 1,085 bifurcations, and 917 amalgamation events. In Fig. 9 we can see that most of the filaments exist only for a short periods of time, although there are a small number of long-lived filaments.

#### Role of tissue heterogeneity.

To determine the effect of tissue heterogeneity on wave patterns during VF, we performed a series of simulations in which we varied the degree of APD restitution heterogeneity in the heart. We did this by using the following linear transformation: (9) where, *S*_{new} is the new value of the maximum restitution slope at each given point, the average maximum slope taken over all points of *dataset 1* is *S*_{avg} = 1.78, ε is a variable that scales the degree of heterogeneity, and *S*_{max} is the original value of the maximum restitution slope at each given point for *dataset 1*. Clearly, ε = 1 corresponds to the fully heterogeneous simulation of *dataset 1* (presented above), ε = 0 corresponds to a spatially homogeneous restitution model with a slope of 1.78, and intermediate values of ε scale the degree of heterogeneity: larger values of ε correspond to more heterogeneous properties.

We performed a series of computations similar to that shown in Fig. 6. The spiral was initiated in the LV wall using a S1–S2 protocol. We used ε values of 0, 0.25, 0.5, 0.75, and 1. The duration of all simulations was 8 s. We investigated the maximum number of filaments, average number of filaments, maximum number of epicardial PS, and average number of epicardial PS for each simulation. The results are shown in Fig. 10, and we see that, for ε = 0 (homogeneous setting), full-blown breakup did not occur for the given slope (1.78). We found that, to induce breakup for our heart model, we needed to increase the restitution slope to >1.8 (results not shown). This is consistent with the slope value reported by ten Tusscher et al. (49), given the difference in the methods used for slope determination (see methods). We see that the maximum number of filaments and the average number of filaments increased for increasing heterogeneity. Similar results were found for the epicardial PS.

To understand the mechanism of increased filament numbers with increasing heterogeneity, we investigated filament locations in the heart. Figure 11 shows all the locations of filament voxels across two slices of the heart together with the maximum restitution slopes (for the 8-s run with ε = 1). The majority of the filaments was located in regions containing steep restitution slopes, and they were evenly distributed across these steep slope regions. Indeed, the mean ± SD restitution slope taken across all voxels traversed by the filaments was significantly steeper (2.17 ± 0.55) than that for the voxels spanning the remainder of the model (1.76 ± 0.65, *P* < 0.005, 2-tailed Student *t*-test).

In our model, the slope of the restitution curve is related to APD (see Fig. 5), and thus gradients in restitution slopes are associated with gradients in APD. Such APD heterogeneities could potentially induce new wave breaks by themselves (even without steep restitution), and this would be expected to occur in areas of abrupt APD variation. However, we did not observe this in our simulations. We observed comparatively fewer numbers of filaments in regions containing restitution slope gradients. Thus the increase in complexity of VF is more likely to be directly due to the increase in size of regions with a steep enough restitution slope (>1.80) rather than an increased occurrence of new breaks at larger gradients of heterogeneities in APD.

An additional fact in favor of this hypothesis is that the electrical turbulence induced by steep restitution would be expected to be homogeneously distributed across entire regions of steep restitution, and this is what we have observed (Fig. 11).

We finally investigated the spatial-temporal distribution of excitation periods for all voxels in the 3D ventricular anatomy for different extents of heterogeneity (ε). For this, we determined excitation periods in every voxel after 2 s of transient processes using a 2-s window. Figure 12 shows the distribution of excitation periods on the epicardium, averaged in a time interval between 2 and 4 s for ε = 1. We see that there are several distinct regions with different excitation frequencies. We followed this distribution until the end of the simulation and found that these regions with different periods are not anatomically predefined and that their location changes in the course of time (data not shown).

We compared the spatial distribution of excitation periods in the interval 2–4 s for all values of ε. For the homogeneous simulation (ε = 0), we found that the average period values were distributed homogeneously and did not change over time. For ε values between 0 and 1, we found that the number of regions with different periods increased as the heterogeneity (ε) in APD slopes increased.

To quantify this increased heterogeneity in excitation period regions, we determined the average period value (dots) and the 95% confidence interval period values (bars) for different values of ε (see Fig. 13). We did this both for excitation period distributions across the entire ventricular mass (gray) and for distributions across the ventricular epicardium (red). Comparing the gray and red bars, we see that the period distributions on the epicardial surface were similar to those measured over the entire ventricular mass. We see that the average periods (gray and red dots) remained approximately equal for all values of ε but that the width of the period distributions increased with increased ε (gray and red bars).

Given the observed temporal variations in excitation period domains and our results on filament locations, we conclude that the increase in complexity of excitation period patterns is primarily caused by the increased number of filaments under increased heterogeneity rather than the increased heterogeneity in anatomically defined APD restitution regions itself.

## DISCUSSION

We have used a detailed human ventricular model to study wave organization during VF in the presence of APD restitution heterogeneity. The model contains a detailed description of human ventricular anatomy, fiber direction anisotropy, and the electrophysiological behavior of human ventricular cells. Recent clinical studies by Nash et al. (31) provided APD restitution properties at 256 locations spanning the entire ventricular epicardium for different patients. This study showed that APD restitution slopes are organized into different regions with slopes ranging from zero to well over three.

We incorporated the clinical data of Nash et al. (31) into our anatomically realistic model of the human ventricles by mapping from the 256 electrode locations to our ventricular geometry and interpolating the projected data using a diffusion-based algorithm. Restitution slopes were then fitted by tuning selected parameters of the TNNP model, which were individually specified at each voxel based on the interpolated restitution slopes. For this study, we used datasets of 2 patients: *dataset 1* contained regions of shallow and steep restitution slopes, whereas *dataset 2* contained mostly shallow restitution slopes. Thus the heterogeneity profiles of both datasets were different.

Heterogeneities are not necessary to induce VF in the human heart. As reported by ten Tusscher et al. (48), VF can occur in a homogeneous heart if the slope of the restitution curve is >1.5. In our model, a slope value of 1.8 or higher is necessary to induce VF. This difference in slope value is due to the difference in the slope determination method used between ten Tusscher et al. (48) and this paper (see methods). Note that the minimum slope value necessary to induce VF in a homogeneous heart is substantially higher than the value of one predicted from the original restitution hypothesis. As shown by Cherry and Fenton (5) in a Fenton-Karma model (13) and by ten Tusscher et al. (49) in the TNNP model, this higher slope value can be explained by electrotonic interactions and CV restitution properties influencing APD alternans.

We analyzed the dynamics of filaments during VF, the formation via birth and bifurcations, the disappearance via death and amalgamation, and the duration of their lifespan. We found that most rotors were short lived, whereas only a small number of long-lived rotors were present during VF. This is in agreement with the findings of ten Tusscher et al. (48) and experimental observations on epicardially manifested rotors in animal hearts (21, 43, 43).

We found that APD restitution heterogeneity is not only important for the initiation of wave break and reentry (6, 9), but also affects the number of filaments during VF. We found that an increased level of APD restitution slope heterogeneity leads to more filaments and a broader and more dynamic distribution of excitation periods. We conclude that there is an interplay between the dynamic organization and complexity of VF and the underlying APD restitution heterogeneity. Furthermore, our results are consistent with the notion that the increase in wave breaks caused by increasing heterogeneity (ε) is predominantly due to the presence of larger areas with steep restitution slopes rather than increased gradients in restitution slope or APD. However, we do not exclude that other factors, such as geometric effects, anisotropy, or more pronounced APD gradients, potentially play a role in wave break formation.

Finally, we found that the number of epicardial PSs correlated well with the total number of filaments and that distribution of epicardial excitation periods provides a good estimate of the period distribution within the entire ventricular mass.

### Limitations

Clinical data were only available for the epicardial surface of the ventricles. There was no information on restitution properties in the midmyocardium nor endocardium. To overcome this, the epicardial restitution slope data were extrapolated across the entire ventricular walls using a diffusion-based algorithm. Note that a variety of interpolation schemes could have been used to achieve a similar degree of variation. The use of different parameters in our diffusion algorithm may lead to different gradients in slope values for the same set of clinical data. It remains to be investigated how these different heterogeneity profiles with different slope gradients affect the results presented in this study.

Other heterogeneities, such as epi-, endocardial, and M-cells, Purkinje fibers, laminar sheets, and disease conditions (such as fibrosis and gap junction remodeling), were not taken into account in our model and may also contribute to the complexity of excitation patterns during VF.

It should be noted that the mapping of the TNNP model parameters to the restitution slope values that we used is not unique and that the same slopes may likely be obtained using other combinations of parameter settings. However, we believe that this nonuniqueness should not qualitatively affect the results of our study, since it is the slope of the restitution curve that is the main determinant of the instabilities in our model.

Another limitation is that the TNNP model used in our study does not reproduce calcium-driven alternans (38, 46, 47), which is an additional source of instability that can lead to VF.

## GRANTS

This research was funded by the Netherlands Organization for Scientific Research (NWO Grant no. 814.02.014). M. P. Nash is supported by the Marsden Fund Council from New Zealand government funding, administered by the Royal Society of New Zealand.

## Acknowledgments

We thank Dr. R. Clayton for valuable comments.

## 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 © 2008 by the American Physiological Society