Abstract
In the present report, we introduce an integrative threedimensional electromechanical model of the left ventricle of the human heart. Electrical activity is represented by the ionic TP06 model for human cardiac cells, and mechanical activity is represented by the NiedererHunterSmith active contractile tension model and the exponential Guccione passive elasticity model. These models were embedded into an anatomic model of the left ventricle that contains a detailed description of cardiac geometry and the fiber orientation field. We demonstrated that fiber shortening and wall thickening during normal excitation were qualitatively similar to experimental recordings. We used this model to study the effect of mechanoelectrical feedback via stretchactivated channels on the stability of reentrant wave excitation. We found that mechanoelectrical feedback can induce the deterioration of an otherwise stable spiral wave into turbulent wave patterns similar to that of ventricular fibrillation. We identified the mechanisms of this transition and studied the threedimensional organization of this mechanically induced ventricular fibrillation.
 ventricular fibrillation
 computer simulations
 mechanics
 electrophysiology
 mechanoelectrical feedback
 stretchactivated channels
mechanical activity of the heart is initiated by electrical waves of excitation that propagate through the heart and initiate cardiac contraction. Abnormal excitation of the heart may result in cardiac arrhythmias and the loss of mechanical pump function, leading to sudden cardiac death. Sudden cardiac death caused by cardiac arrhythmias is the most common cause of death in the industrialized world, and, in most cases, this is due to ventricular fibrillation (VF) (73). It has been shown in clinical and experimental studies (11, 17, 26, 41, 46, 62, 68, 70) that VF occurs as a result of the onset of turbulent electrical activation patterns of the heart that are underpinned by multiple reentrant sources of excitation. Mechanisms behind the onset of reentrant sources in the heart and processes resulting in breakup of these sources into complex turbulent activation patterns are of great interest, e.g., in the design of therapeutic strategies to prevent or treat cardiac arrhythmias.
One of the important factors that affects electrical excitation of the heart is mechanoelectrical feedback. It has been shown that mechanical deformation alters the electrical properties of myocytes via stretchactivated channels (59), which can change the shape of the action potential in response to stretch (34, 37). Mechanoelectric feedback has been studied in the clinical community for well over a century (for reviews, see Refs. 34, 37, 38) and may have both proarrhythmic and antiarrhythmic consequences. However, the mechanisms underlying these phenomena are not well understood. One of the problems in experimental studies of mechanoelectrical feedback is that it is difficult to control and record both the electrical and mechanical properties of cardiac tissue and to alter them in a systematic way. Therefore, alternative methods, such as mathematical modeling, must be used to address this field.
We (47) recently investigated the effects of mechanoelectrical feedback on wave propagation using a twodimensional (2D) explicitly coupled electromechanical model. We found that it is possible to induce automatic pacemaker activity (52) and spiral wave breakup (53). Both mechanisms are regulated by mechanoelectrical feedback via stretchactivated channels. However, in these studies, we chose to use a simple phenomenological model of cardiac excitation in a 2D isotropic medium. Active tension was generated using an oversimplified isotropic active tension transient. Furthermore, we used an isotropic MooneyRivlin material response to describe the passive mechanical properties of the medium.
In the present study, we extended our analyses to an anatomically and functionally realistic electromechanical model of the human left ventricle (LV). We used this approach to study the interactive effects of wave dynamics on mechanical deformation as well as the effect of mechanoelectrical feedback on wave dynamics. To this end, in the present report, we introduce a new realistic threedimensional (3D) anisotropic electromechanical model of the human LV that combines electrical waves of excitation and mechanical deformation. We validated our model by simulating a control heart beat and compared the deformations with those found experimentally by Ashikaga et al. (2). We used our model to study electrophysiological and mechanical dynamics during reentry in the presence of cardiac contraction. We found that turbulent patterns of excitation during ventricular tachycardia (VT) and VF resulted in desynchronized cardiac contraction compared with a control heart beat. Furthermore, we found that dynamic regions of stretch can cause an otherwise stationary reentrant wave to break up into smaller waves via mechanoelectrical feedback of stretchactivated channels.
METHODS
In the present report, we introduce a coupled 3D electromechanical model of the human LV that combines electrical waves of excitation and mechanical deformation. In this model, electrical excitation waves are described using a detailed ionic cell model of the human ventricular myocyte (64). Mechanical deformations are represented using finite elasticity theory combined with finite element analysis. Based on the local intracellular Ca^{2+} concentration, contractile forces are represented using equations governing the development of tension due to cellular crossbridge cycling (48). This active tension is coupled to nonlinear stress equilibrium equations, which are solved using the finite element method. In turn, the resulting mechanical deformations have a feedback effect on the excitation properties of the membrane via tissue deformation and stretchactivated channels. The electrical and mechanical components each have their own numerical discretization. For the electrophysiological model, we created a detailed anatomic finite difference geometry of the LV that contains fiber direction anisotropy. For the mechanical model, we created an equivalent tricubic Hermite finite element geometry that contains the same fiber direction anisotropy. Numerically, the coupled electromechanical model combines a finite difference method approach to integrate the excitation equations, with a Galerkin finite element method to solve the nonlinear equations governing the tissue mechanics. More details are contained in the following sections.
Cardiac Excitation
Excitable behavior was modeled using a monodomain description of cardiac tissue. For this, we used the following partial differential equation (32):
Diffusion tensor D_{ij} is derived from the fiber direction field. Assuming that the transverse conductivity is the same in all directions orthogonal to the direction of the muscle fiber axis, we described the ventricular conductivity tensor using the following equation:
Cardiac Mechanics
To model deformation of the cardiac tissue, the mechanical analysis was based on finite deformation elasticity theory. The normal sarcomere length (SL) of cardiac muscle cells is 2.0 μm (27, 56). In mammals, SLs operate in vivo normally between 1.7 and 2.15 μm (58) and in vitro maximally between 1.5 and 2.3 μm (4, 21, 65, 66). Thus, cardiac cells change their length by 15–20% during contraction (40), and their maximal stretch relative to resting length is ∼5–10%.
Following standard continuum mechanics, we used two types of coordinates to refer to a material point in the tissue. Here, x = {_{xi}} describes the present (deformed) spatial position in rectangular Cartesian coordinates of a material particle that occupied the location X = {X_{M}} in the reference (undeformed) configuration. The deformation gradient tensor (F) transforms the undeformed line segment (dX) into the deformed line segment (dx) by dx = F × dX with F{M^{i}} = ∂x_{i}/∂X_{M}. The right CauchyGreen deformation (metric) tensor, C = F^{T}F, describes how each component of the undeformed line segment dX contributes to the squared length of the deformed line segment dx and is independent of rigid body motion. The Green Lagrange strain tensor (E) is defined as follows: E = 1/2(C − I), where I is the unitary tensor.
Stress equilibrium.
Due to active tension in the cardiac wall and the pressure load on the endocardial surface, the LV model will deform until a new equilibrium state is reached. To determine the deformed state, the stress equilibrium equations are solved using a finite element method. These equations arise from the conservation of linear momentum following Newton's law of motion. In the absence of body forces, the static equilibrium equations reduce to the following:
We assumed that the passive myocardium can be modeled as an incompressible hyperelastic material, described by a transversely isotropic mechanical response (18), as follows:
To represent active contractile tension T_{a}, we used the NiedererHunterSmith (NHS) model (48). The NHS contraction model uses three input variables: the intracellular Ca^{2+} concentration, the fiber extension ratio, and the velocity of fiber shortening. The intracellular Ca^{2+} concentration is generated by the TP06 model (Eq. 1), whereas the other two parameters are obtained from the 3D mechanics model. The NHS model is based on experimental data from rat and guinea pig hearts at room temperature and was designed to operate under physiological SLs in the range 1.8 < SL < 2.3 μm. We adjusted the parameters of the NHS model to account for the longer action potential and different Ca^{2+} transients of the TP06 model. On the advice of Dr. Niederer, we first sped up the relaxation kinetics of the NHS model by setting α_{r1} = 10 s^{−1} and α_{r2} = 25 s^{−1}, which is consistent with the faster cardiac relaxation times reported for higher body temperatures (12, 55). We also set the parameters for contractile tension to T_{ref} = 125 kPa and pCa_{50} = 6.5, which is again consistent with the temperature increase (20). Such modifications were necessary because when the original parameter values were used, the longer action potential duration and Ca^{2+} dynamics in human ventricular myocytes (represented by the TP06 model) resulted in the accumulation of contractile tension during higher stimulation frequencies in the NHS model, even though intracellular Ca^{2+} was restored to resting values between stimuli. With these new model parameters, the NHS model generated maximal contractile tension values of ∼100 kPa at resting SL. For further model information and parameters, see Ref. 48.
Mechanoelectrical Feedback
Tissue deformation.
The equations describing electrical activity (Eq. 1) can be rewritten by taking local coordinate deformation into account:
Stretchactivated channels.
The direct influence of deformation on the excitation properties is via stretchactivated current (I_{s}), which was added to the TP06 model. There are three main groups of mechanically activated channels in the heart, but only two of them (cation nonselective channels and K^{+}selective channels) are activated by stretch (37). The overall physiological action of these channels is depolarization of the membrane in response to stretch, as shown in the majority of experimental observations from isolated cardiac tissue. Experimental studies (25, 72, 74) of the electrophysiological properties of stretchactivated channels have shown that they are instantly activated by mechanical stimulation and that the currentvoltage relationship for the most important nonspecific cation channels is linear. On basis of these observations, linear ionic models for I_{s} have been proposed (67, 69). These linear models have been used to study the effects of mechanical stretch on heart tissue using detailed ionic models of cardiac myocytes. We followed this approach and used the following nonselective linear timeindependent description of I_{s}:
The value of E_{s} in most biophysical models is assumed to be around −20 mV (25, 67, 69) and describes the experimentally observed depolarizing effect of I_{s}. The value of G_{s} varies up to 100 μS/μF (36, 37, 67). We used G_{s} values between 0 and 75 μS/μF to investigate the effects of I_{s} on the excitation properties.
Human LV Anatomy
The geometric data describing the 3D ventricular anatomy and fiber direction field were derived from a structurally normal human heart (23, 24) and are described in more detail in Ref. 63. This dataset had an isotropic resolution of 0.5 mm. The muscle fiber direction field was constructed based on general knowledge of the fiber architecture in the human heart (60) combined with detailed data on the fiber architecture of the canine heart (50). We removed points corresponding to the right ventricle for the human dataset. The remaining 1,118,053 voxels were used as a finite difference mesh for the electrical part of our framework.
For the mechanical equations, we developed a finite element mesh based on the finite difference data. We extracted the endocardial and epicardial surfaces from the LV dataset. To fit a finite element mesh to these data, we created an initial mesh in the shape of a truncated ellipsoid using tricubic Hermite hexahedral elements. The mesh consisted of one element in the radial direction, five elements in the longitudinal direction, and eight elements in the circumferential direction (40 tricubic Hermite elements containing 2,386 geometrical degrees of freedom). The model was then fitted using a nonlinear geometric facefitting procedure (8). The endocardial and epicardial surfaces of the mesh were fitted to 34,177 and 103,319 data points, respectively. For this, we used an iterative method in which the Euclidean distance between all data points and their orthogonal projections to the model mesh were minimized.
The 3D model fiber angles were fitted using a similar nonlinear optimization method. The fiber angle field contained 656 degrees of freedom, which were fitted to 10,904 fiber angles. Through the LV wall, fiber angles varied from 13 (endocardium) to −55° (epicardium) in the anterior region, 72 to −32° in the posterior wall, 62 to −57° in the lateral free wall, and 14 to −37.3° in the septum (elevation angles are defined with respect to the local shortaxis plane). These fiber orientations were similar to measurements of fiber orientations reported in Ref. 60. The resulting root mean square error after fitting was 10°, similar to that reported in Ref. 15. For more information on how to create a mathematical model of geometry and fibrous structure of the heart, we refer to Refs. 8 and 50.
Numerical Methods
The coupled electromechanical model was solved using a hybrid approach that incorporates a finite difference method to solve the electrical equations and finite element techniques to compute the deformation of the tissue. The electromechanical simulations are typically performed on two distinct time scales: one for the electrophysiology and a longer time scale for the mechanics.
For electrical excitation, we used a finite difference method with an adaptive time step of up to Δt = 0.02 ms on a grid containing 1,118,053 points with a space step of Δx = 0.5 mm. The Laplacian in Eq. 6 was computed using weights calculated at each point based on the local conductivity and metric tensor. To impose noflux boundary conditions, the current flow through the surface of the heart to nonheart points was set to zero by adjusting the weights of the boundary points (63).
Mechanical equations (Eq. 3) were solved using CMISS (www.cmiss.org), which was developed at the University of Auckland. Each mechanical step was followed by 100 electrical steps. Excitationcontraction coupling was achieved by passing spatial distributions of contractile tension from the electrophysiology model to the mechanical component of our framework. The mechanical boundary conditions included fixation of the upper basal ring and application of a constant endocardial cavity pressure load of 10 kPa. After the solution of each mechanical deformation, metric tensors were updated and passed back to the electrophysiology model at 41,380 points, which were further interpolated to all electrical elements. These metric tensors were then used to compute the new weights for evaluation of the Laplacian in Eq. 6 and to calculate the current generated via stretchactivated channels (Eq. 7).
To simulate cardiac excitation during normal sinus rhythm, we stimulated the entire endocardium with a period of 600 ms. To initiate 3D scroll waves, we used an S1S2 protocol, for which an S1 stimulus was applied in the septum. The S2 stimulus was activated during the refractory tail of the S1 stimulus on the posterior LV wall and was extended from the base to ∼50% of the baseapex dimension, thereby creating a single reentrant scroll wave on the posterior LV. All stimuli were applied at twice the diastolic threshold for a duration of 10 ms.
Pseudobody surface ECGs were computed using the infinite volume conductor approach (57) as follows:
Simulations were run on a Dell Precision Workstation (dual Intel Xeon 2.8 GHz, 2 GB RAM). Simulating 1 s of wave propagation and deformation of the LV took ∼18 h of wall clock computation time, which was split approximately equally between the electrical and mechanical components of our model. Ventricular geometry, wave patterns, and scroll wave filaments were visualized using the marching cubes algorithm for isosurface detection in voxel data and OpenGL for isosurface rendering. Mechanical deformation was visualized using cmgui (www.cmiss.org/cmgui).
RESULTS
Model Validation
We compared our results with the experimental deformations reported by Ashikaga et al. (2) (Fig. 3, atrial pacing). In five canine hearts, they recorded strains referred to local cardiac coordinates in the anterior wall of the LV (between the first and second diagonal branches of the left anterior descending coronary artery) during normal atrioventricular conduction. We only compare endsystolic strains, because Ashikaga et al. (2) refer strains to the enddiastolic state.
We computed strains with respect to cardiac and fiber coordinates under control conditions and compared these with experimental values (2). With reference to Fig. 1A, cardiac coordinate were defined by the orthogonal system (c, l, and r), where c is the shortaxis circumferential direction tangent to the wall plane, l is the longitudinal axis tangent to the ventricular surface from the apex to base, and r is the transmural axis directed from the endocardium to epicardium. We also defined fiber coordinates using the orthogonal system (f, t, and r), where f is the myocyte (fiber) axis, t lies transverse to the fiber axis within the ventricular wall plane, and r is the radial coordinate, as defined above (Fig. 1B).
Figure 2 shows the epi, mid, and endocardial cardiac strains at three different transmural points approximately halfway between the base and apex within the anterior LV wall. We showed cardiac normal strains [circumferential (E_{cc}), longitudinal (E_{ll}), and radial (E_{rr})] and cardiac shear strains [circumferentiallongitudinal (E_{cl}), longitudinalradial (E_{lr}), and circumferentialradial (E_{cr})]. Note that for comparison with the experimental recordings from Ashikaga et al. (2), we calculated these cardiac strains with respect to the enddiastolic configuration with a LV pressure of 1 kPa. It should be noted that no tuning of the mechanical constitutive parameters was performed as part of this study.
We found that the normal strains E_{cc}, E_{ll}, and E_{rr} (denoted as E_{11}, E_{22}, and E_{33}, respectively, in Ref. 2) were generally in good agreement, although our model slightly overestimated endsystolic wall thickening. The transmural gradients in the cardiac normal strains also compared well with the experimental data. The epicardial inplane shear strain E_{cl} (E_{12}) matched well; however, endocardial E_{cl} was small. The magnitude and transmural gradient of the longitudinal transverse shear E_{lr} (E_{23}) were consistent; however, E_{lr} was of the opposite sign to the recordings. This may be explained by the lack of orthotropic sheet structure in our model, which has been shown to affect transmural shear (10). The shortaxis transverse shear E_{cr} (E_{13}) was also in good agreement. All transmural gradients in the cardiac shear strains were consistent with the experimental recordings.
We also computed several electrical and mechanical characteristics of the model during normal cardiac excitation. Figure 3 shows the pseudoECG, which resembles a normal ECG during sinus rhythm. The QRS complex corresponds to the depolarization of the ventricle and lasted ∼90 ms. The T wave represents the repolarization (or recovery) of the ventricle, and the QT interval was ∼360 ms, which was similar to clinically recorded ECGs during control heart cycles (19). Figure 3 also shows the volume fraction, which represents the deformed cavity volume as a proportion of the enddiastolic volume. This corresponds with an endsystolic LV ejection fraction of 58%, which is in accordance with that of the healthy heart (19).
For LV endocardial pacing, we found multiple phases of deformation, as exemplified by the fiber extension ratio λ_{f} (see Fig. 3). We found fiber shortening (λ_{f} < 1) during systolic contraction, over which myocytes change their length by up to 12%. After systole, myocyte stretch (λ_{f} > 1) occurred due to the LV cavity pressure of 10 kPa applied to the endocardial surface. Myocyte stretch relative to the resting length was between 5% and 10%.
Simulating a Stable Reentrant Spiral
As described in methods, we initiated a spiral wave on the posterior wall of the LV. For this simulation, we set up the conductivity of the stretchactivated channels to be zero (G_{s} = 0). The spiral wave remained stable for the entire duration of the simulation. Figure 4 shows the ECG, which was periodic (dominant frequency of 4.0 Hz), resembling ECG signals observed during VT. The dominant frequency corresponds to the period of the spiral wave. Figure 4 also shows the relative change in the LV cavity volume expressed as a fraction of the enddiastolic volume. The first two oscillations correspond to the initiation of the spiral wave during the S1S2 protocol. Subsequently, we observed that the average volume of the LV cavity decreased by ∼50% and oscillated with an amplitude of 8–10% and a dominant frequency of 4.0 Hz. The slow rise in volume fraction can be explained via intracellular Ca^{2+}; at higher frequencies, the peak values of intracellular Ca^{2+} slightly decrease, and thus less contractile tension is generated.
In addition, Fig. 4 shows the epi, mid, and endocardial fiber extension ratio (λ_{f}) at the same locations as in the control simulation (Fig. 2). We observed fiber shortening (λ_{f} < 1) and local wall thickening (not shown) during reentry. Fiber extension oscillated with the same frequency as that of the spiral wave. We also found that the 3D distribution of λ_{f} varied throughout the LV during reentry and that multiple regions of the LV wall were subject to fiber stretch (λ_{f} > 1) at various times. These regions of fiber stretch dynamically shifted across the LV myocardium and were dependent on the phase of the spiral wave andlocation of wavefronts and were often located ahead of the wavefront adjacent to contracted tissue behind the wavefront.
Mechanically Induced Wavebreak
To investigate mechanically induced wavebreak, we incorporated stretchactivated channels with a conductivity of G_{s} = 75 μS/μF into the electrophysiology model and induced a spiral wave on the posterior wall of the LV. All other parameters were the same as for the stable spiral wave simulation. After the initiation of the spiral, we observed wavebreaks, which were enhanced over time by the stretchrelated processes. From the wavefronts (Fig. 5), we observed that the initial spiral wave on the posterior wall remained stable throughout the simulation. We also observed wavebreaks on both the posterior and anterior walls, but not close to the core of the initial spiral wave.
Figure 6 shows the pseudoECG for this simulation. The ECG was more irregular and complex compared with that of the stable spiral simulation. The dominant frequency was 3.6 Hz, similar to frequencies reported in clinical (9, 43, 46, 71) and numerical (63) studies of VF. The slightly lower frequency compared with control can be explained by the activation of I_{s}, which leads to a localized decrease in Na^{+} current, resulting in lower conduction velocities ahead of the wavefront. Figure 6 also shows the time dynamics of the number of filaments found in the LV during 8 s of simulation. Over the first 4 s, the number of filaments steadily increased from 1 (the initial spiral) to over 15. Between 4 and 8 s of time, the number of filaments slightly decreased and oscillated around 13. The average number of filaments found during the last 4 s was 13.6. In addition, Fig. 6 shows a typical example of the number of filaments present during developed VF. As shown in Figs. 5 and 6, there were clearly multiple independent filaments, while the initial spiral wave remained intact.
After the initial phase and resulting mechanically induced wavebreaks, we found that the LV remained in a contracted state, similar to that of the stable spiral simulation. The LV cavity volume fraction, mechanical 3D deformations, and fiber strains were similar to those of the stable spiral simulation shown in Fig. 4.
We analyzed how filaments were created and destroyed over time through death, birth, bifurcation, and amalgamation events (see Fig. 7). The horizontal lines in Fig. 7 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 righthand end of each line indicates when a filament disappears due to the death or the amalgamation of a filament with another filament. In Fig. 7, the short horizontal lines correspond to shortliving filaments and long lines correspond to more persistent filaments. During the 8 s of simulation time, we detected 1,276 filaments with 453 births, 496 deaths, 824 bifurcations, and 781 amalgamation events. As shown in Fig. 7, the initial spiral remained intact for the entire 8 s (bottom red line). Most of the filaments existed for a short period of time, although we also observed a small number of persistent independent filaments.
We investigated the effects of different levels of activation of stretchactivated channels (G_{s}) on the VF dynamics. For this, we repeated the above simulations with G_{s} set to 0, 30, 45, 60, and 75 μS/μF and analyzed the number of filaments present during the simulation (see Fig. 7). For G_{s} = 30 μS/μF, we found that the simulation remained stable during 4 s of simulation time, although a small number of wavebreaks were observed toward the end of the simulation. For larger G_{s} values, we observed that the number of filaments steadily increased over time with more filaments created for larger values of G_{s}, thereby also increasing the complexity of the VF dynamics.
The filament history dynamics shown in Fig. 7 resemble those of mother rotor fibrillation (33). However, this may just be a visual coincidence, since based on these results alone we cannot claim that the wavebreaks are driven by the initial spiral wave, nor can we claim that these breaks disappear if the initial spiral wave is removed. On the other hand, we can explain such filament dynamics by the maximal stretch distribution observed during the simulation. Figure 8 shows the maximal fiber extension values observed at each epicardial point on both the anterior and posterior free walls during the first 4 s of the simulation. We observed regions of stretch (λ_{f} > 1) on the anterior free wall, where most breaks are located. In contrast, we observed that stretch did not occur (λ_{f} < 1) close to or around the core of the initial spiral wave on the posterior wall.
The mechanism by which stretch of the tissue induces wavebreak can be explained in the following way. Stretchactivated channels produce an inward current, which depolarizes the myocyte membrane. This depolarization is not sufficient to generate an action potential (i.e., it is subthreshold); however, it leads to a quite opposite effect: a depression of excitability. This occurs as a result of voltagedependent inactivation of the Na^{+} channel. In 1952, Hodgkin and Huxley (22) showed that slow depolarization of a nerve cell membrane results in the progressive inactivation of Na^{+} channels. Such inactivation was also found in cardiac cells (see, e.g., Ref. 3) and is present in all ionic models of cardiac cells. Thus, stretchinduced membrane depolarization inactivates Na^{+} channels, which carry the main ionic current underlying the initiation of the action potential upstroke. As a consequence, a wave that arrives in a region of substantial stretch cannot propagate through that region, and wavebreak is formed. We (53) have previously presented a detailed study of this process using a lowdimensional model of cardiac tissue. We have also performed simulations using a higherdimensional model, which combined the mechanics used in Ref. 53 with the biophysical TP06 ionic model of cardiac excitation. Video 1 in the Supplementary Material shows the results of simulations performed in a simple 2D geometry on a grid of 705 × 705 elements with space step of 0.25 mm.^{1} We found that in this case, I_{s} induces the break up of an otherwise stable single spiral wave into a complex spatiotemporal pattern. Since these 2D simulations were performed using a space step of 0.25 mm, this also indicates that the break up observed in this study is not a result of a coarse spatial discretization. Upon further investigation, we confirmed that the mechanism of this effect is, again, the voltagedependent inactivation of Na^{+} channels.
DISCUSSION
In this study, we developed a 3D strongly coupled electromechanical model of the human LV to investigate the interactive effects of electrical waves and mechanical contraction during reentry. This 3D electrophysiological model contains a detailed description of human LV anatomy and fiber direction anisotropy and represents the electrophysiological properties of human ventricular myocytes using the ionic TP06 model (64). This model contains detailed descriptions of voltage, ionic currents, and intracellular ion concentrations and is based on a wide range of human electrophysiological data. The intracellular Ca^{2+} concentration of the TP06 model was used as input for the biophysical NHS model of cardiac tissue contraction (48), which represents the contractile tension developed due to molecular crossbridge cycling. This active tension is coupled to a 3D finite element model of passive mechanics, which represents the same human LV anatomy and fiber direction anisotropy as described above for the electrophysiology model. The mechanical model was subjected to a transversely isotropic exponential material response to describe the passive properties. The stress equilibrium equations were solved to calculate the deformed state of the LV and the corresponding local metric tensors. Mechanoelectrical feedback was introduced via local metric tensors by updating the local coordinate system and by calculating the local deformation, which modulates I_{s}.
We validated our LV electromechanical model using a control excitation, for which we stimulated the entire LV endocardium. We observed ∼10% fiber shortening (compared with resting length) and ∼15% wall thickening during systolic contraction. Endsystolic cardiac strain distributions were similar to experimental recordings reported from canine hearts during normal atrioventricular conduction (2). The endsystolic ejection fraction was 58%, corresponding to that of a healthy human heart.
The mechanisms of wavebreak studied in this report are substantially different from other mechanisms of breakup, such as negative filament tension (1, 5), fiber rotation (14), and restitution (54). The onset of breakup by our mechanism does not require the presence of any of the above factors. Mechanically induced breakup is strongly associated with the deformation of the heart and the presence of stretchactivated channels and will not occur in the absence of either of these factors.
One obvious experimental validation of the effects of stretchactivated channels on wave breakup might be a study that compares the onset of fibrillation in a normal heart and in a heart in which stretchactivated channels are blocked (e.g., by gadolinium). One study (6) has shown that gadolinium indeed decreases the vulnerability of rabbit hearts to atrial fibrillation; however, to the best of our knowledge, such experimental studies have not yet been performed for ventricles. It should also be noted that the interpretation of the results of such experiments may be not so straightforward since it is well known that, in addition to block of stretchactivated channels, gadolinium has additional important membrane effects (see, e.g., Ref. 75).
Although our data show a new possibility of wavebreak in the heart, we do not claim that all VF episodes occur due to mechanoelectrical feedback. There are, of course, multiple experimental studies showing that VF can be induced during complete block of cardiac contraction caused by the application of excitationcontraction decouplers used in optical mapping studies (see, e.g., Refs. 13 and 17). On the other hand, the mechanisms considered in the present study are more likely to occur in conditions where stretch of the tissue is more substantial or is combined with other pathological factors, such as ischemia or infarction. In fact, the results reviewed by Janse et al. (29) show that the trigger for VF during phase 1B of ischemia is most likely related to stretch [see also a recent modeling study (30)].
During mechanically induced breakup, we did not observe any wavebreaks close to the core of the spiral, resulting in mother rotorlike VF patterns. Especially in the initial phase of the simulation, we found that wavebreaks were small and dependent on the initial spiral wave. The complexity of the simulation increased as time increased, enhancing wavebreaks over time due to the stretchrelated processes. The complexity of VF dynamics is dependent on the amount of stretch generated (e.g., determined by λ_{f}, intracellular Ca^{2+}, T_{a}, passive material properties, etc.) and the magnitude of I_{s} (determined by G_{s} and E_{s}). We observed that the complexity of VF and number of filaments increased when G_{s} increased.
These preliminary results demonstrate that mechanoelectric feedback can cause an otherwise stationary spiral wave to break up, thus providing one possible explanation for the degeneration of VT to VF activity. The effects shown here should be considered as an indication of the importance of mechanics on wave dynamics and not as a firmly established mechanism. Further experimental research and numerical studies for variations of initial conditions, boundary constraints, and model parameters are needed to quantify the effects of mechanoelectrical feedback on wave dynamics. In general, mechanoelectrical feedback has multiple effects on the electrical activity of the heart, which are not limited to the factors studied in this report. Extensive reviews and a dedicated editorial on this subject are contained in Ref. 38.
Limitations
We adopted an isotonic loading regime with a constant LV endocardial pressure of 10 kPa. During the normal cardiac cycle, there are variations in LV pressure and volume. This can be improved using a Windkessel afterload model (7, 35). However, such pressurevolume loops are typically not present during VT or VF since cardiac output is often significantly compromised.
Our geometric model only contained fiber orientations. Thus, the LV model did not account for sheet or imbrication angles, which may affect local mechanics, such as transverse shear deformations, or the strains near the apex and base (15). In addition, the effects of the right ventricle on electromechanical dynamics were not accounted for in this study, but this will be addressed in future research.
The mechanical properties of the LV myocardium were described using a transversely isotropic material relation, with parameters based on experimental data from the canine ventricular myocardium. The parameters of the active tension model were based on experimental data from rats. For the present purposes, we made minor adjustments to some parameters to avoid tension accumulation during highfrequency contractions. A similar active contractile mechanics model needs to be developed based on human experimental contraction data to determine the correct parameters and maximum tension values.
Mechanoelectrical feedback was modeled via the modulation of diffusion pathways and I_{s}. Mechanical deformation may also affect the passive electrical properties and other ionic currents of cardiac myocytes. At present, however, limited experimental data are available upon which represent such feedback mechanisms.
While we have identified some important issues above, the conclusions of the present study do not depend on these limitations. Indeed, these additional features and effects may introduce other interesting dynamics and electromechanical interactions.
GRANTS
This work was funded by The Netherlands Organization for Scientific Research Grants 814.02.014 and 635.100.004. A. V. Panfilov was partially supported by the Institut des Hautes Etudes Scientifiques (BuressurYvette, France) and the James Simons Foundation. M. P. Nash is supported by a James Cook Research Fellowship from the Royal Society of New Zealand.
DISCLOSURES
No conflicts of interest, financial or otherwise, are declared by the author(s).
ACKNOWLEDGMENTS
The authors gratefully acknowledge advice provided by Dr. S. Niederer regarding the parameter adjustment of the NHS contraction model and thank Dr. R. Hren for providing anatomic data for the simulations.
Footnotes

↵1 Supplemental Material for this article is available online at the American Journal of PhysiologyHeart and Circulatory Physiology website.
 Copyright © 2010 the American Physiological Society
REFERENCES
 1.↵
 2.↵
 3.↵
 4.↵
 5.↵
 6.↵
 7.↵
 8.↵
 9.↵
 10.↵
 11.↵
 12.↵
 13.↵
 14.↵
 15.↵
 16.↵
 17.↵
 18.↵
 19.↵
 20.↵
 21.↵
 22.↵
 23.↵
 24.↵
 25.↵
 26.↵
 27.↵
 28.↵
 29.↵
 30.↵
 31.↵
 32.↵
 33.↵
 34.↵
 35.↵
 36.↵
 37.↵
 38.↵
 39.↵
 40.↵
 41.↵
 42.↵
 43.↵
 44.↵
 45.↵
 46.↵
 47.↵
 48.↵
 49.↵
 50.↵
 51.↵
 52.↵
 53.↵
 54.↵
 55.↵
 56.↵
 57.↵
 58.↵
 59.↵
 60.↵
 61.↵
 62.↵
 63.↵
 64.↵
 65.↵
 66.↵
 67.↵
 68.↵
 69.↵
 70.↵
 71.↵
 72.↵
 73.↵
 74.↵
 75.↵