|Home | About | Journals | Submit | Contact Us | Français|
The excitation-contraction coupling properties of cardiac myocytes isolated from different regions of the mammalian left ventricular wall have been shown to vary considerably, with uncertain effects on ventricular function. We embedded a cell-level excitation-contraction coupling model with region-dependent parameters within a simple finite element model of left ventricular geometry to study effects of electromechanical heterogeneity on local myocardial mechanics and global hemodynamics. This model was compared with one in which heterogeneous myocyte parameters were assigned randomly throughout the mesh while preserving the total amount of each cell subtype. The two models displayed nearly identical transmural patterns of fibre and cross-fibre strains at end systole, but showed clear differences in fibre strains at earlier points during systole. Hemodynamic function, including peak left ventricular pressure, maximum rate of left ventricular pressure development, and stroke volume were essentially identical in the two models. These results suggest that in the intact ventricle heterogeneously distributed myocyte subtypes primarily impact local deformation of the myocardium, and that these effects are greatest during early systole.
In the most basic sense, electromechanical coupling occurs at the level of the individual myocyte where changes to membrane potential lead to release of Ca2+ and production of force by the myofilaments. This intrinsic activity is modulated by surrounding myocytes. Ionic currents through gap junctions electrically couple adjacent cells, allowing propagation of action potentials through the myocardium. Mechanical coupling among myocytes arises from cytoskeletal structures linking myocardial cells to each other and to the extracellular matrix. Thus, in a broader sense, electromechanics includes electrical and structural properties of the intact myocardium as these influence the local environment in which myocytes operate.
Regional deformation of the ventricular wall during the cardiac cycle is influenced by myocardial tissue structure (Arts et al. 1979; Arts et al. 1994; LeGrice et al. 1995). Early modeling studies suggested that ventricular torsion and transmural patterns of fibre orientation act to minimize transmural differences in sarcomere length and fibre strain, both during ejection (Arts et al. 1979) and filling (Guccione et al. 1995), thereby distributing the work of ejection equally through the heart wall (Arts et al. 1979; Arts et al. 1994). Studies using biplane radiography of radiopaque markers embedded in the ventricular wall (Costa et al. 1999; Mazhari et al. 2001) as well as MRI tagging methods (Tseng et al. 2000) showed patterns of transmural fiber strain that were in general agreement with those predicted in models.
In recent years, computational and experimental studies have revealed a more complex scenario. Kerckhoffs et al. (2003a) showed using a computational model of LV electromechanics that a physiological pattern of electrical activation produced large gradients of transmural fibre strain in comparison with a synchronously activated model. The authors postulated that a transmural gradient of intrinsic myocyte electromechanical delay opposing that of electrical activation would reduce transmural heterogeneity of systolic fibre strain. Shortly thereafter, a study of isolated canine myocytes demonstrated the existence of just such a transmural distribution, with endocardial myocytes (those earliest activated in the intact heart) showing increased latency to onset of shortening in comparison with midmyocardial and epicardial cells (Cordeiro et al. 2004). At the same time, improved spatial and temporal resolution in the measurement of fibre strains have revealed transmural heterogeneity in the time courses of regional fibre strains during systole (Ashikaga et al. 2004; Ashikaga et al. 2007). These new results do not invalidate the overall hypothesis that realistic LV torsion and fibre orientation act to reduce transmural variation of fibre strain. But, they illustrate the need for more detailed examination of the mechanisms underlying observed transmural deformations during systole.
The most advanced models of ventricular electromechanics depict threedimensional action potential propagation, contractile force generation, and myocardial deformation under realistic hemodynamic loading conditions and within anatomicallydetailed finite element meshes (Usyk et al. 2002; Kerckhoffs et al. 2003a; Usyk & McCulloch, 2003a; Nash & Panfilov, 2004; Watanabe et al. 2004; Nickerson et al. 2005; Sermesant et al. 2005; Kerckhoffs et al. 2008; Niederer & Smith, 2008). The degree of coupling between component models varies, but newer studies have made use of algorithmic advances which allow full (two-way) coupling between and quasisimultaneous solution of action potential propagation, ionic currents, myofilament force production, and tissue mechanics (Nash & Panfilov, 2004; Watanabe et al. 2004; Nickerson et al. 2005; Niederer & Smith, 2008). This form of coupling is useful for investigating questions relating to mechano-electric feedback (Kohl et al. 1999; Niederer & Smith, 2007), where length changes on the level of myocytes influence Ca2+ buffering by myofilaments and the activity of stretch activated ion channels. In many other cases, it has been sufficient to obtain a spatiotemporal sequence of depolarization from the solution of electrophysiology equations which can then be used to drive myofilament force production in a separate simulation of ventricular mechanics (Usyk et al. 2002; Kerckhoffs et al. 2003b; Kerckhoffs et al. 2003a; Usyk & McCulloch, 2003a; Usyk & McCulloch, 2003b; Kerckhoffs et al. 2005).
In recent years a number of fully-coupled cellular electromechanics models have been published. These models simulate various cell responses by simultaneously solving the ordinary differential equations (ODEs) governing ionic currents, myofilament force production, and one-dimensional cell mechanics (Niederer & Smith, 2007; Rice et al. 2008). We recently developed one such model to explore potential sources of heterogeneous excitation-contraction (EC) coupling behavior previously observed in isolated canine myocytes (Campbell et al. 2008). Here, we use this model to examine the effects of heterogeneous myocyte properties (Sicouri & Antzelevitch, 1991; Sicouri et al. 1994; Cordeiro et al. 2004) on transmural systolic fibre mechanics and ventricular function. In particular, we investigate the hypothesis of Kerckhoffs et al. (2003a) that cellular heterogeneity of myocyte EC coupling helps maintain transmural uniformity of systolic fibre strain in the presence of a physiological activation delay.
In the present study, we embedded our recently published canine myocyte model (Campbell et al. 2008) within a model of ventricular action potential propagation and mechanics. This fully-coupled electromechanical model was used to predict the integrative effects of myocyte EC coupling heterogeneity on ventricular wall mechanics and global pump function during ejection.
A model of canine myocyte ionic currents and cytosolic Ca2+ cycling (Greenstein et al. 2006) was previously expanded and modified to reproduce observed transmurallyvarying characteristics in action potential duration, early repolarization, and Ca2+ transient time course (Flaim et al. 2006). Three parameter sets for the model were developed corresponding to published measurements in endocardial (endo), midmyocardial (mid), and epicardial (epi) cells.
In previous work (Campbell et al. 2008), we developed parameter sets for the myofilament model of Rice et al. (2008), which quantitatively reproduced observed characteristics of unloaded shortening in the three canine myocyte subtypes driven by experimentally-derived Ca2+ transients (Cordeiro et al. 2004). In that same study, we combined this re-parameterized version of the Rice model with the excitation-Ca2+ coupling model of Flaim et al. (2006) to create a fully-coupled model of canine myocyte electromechanics (Campbell et al. 2008). The Rice model itself uses a system of coupled ODEs to reproduce the dynamic response of cardiac muscle to time-varying inputs of sarcomere length and cytosolic Ca2+. The model is capable of reproducing a variety of classical responses such as steep activation of force by Ca2+ and a realistic force-velocity relation. Length dependent activation is represented by a model of thick and thin filament overlap according to the classic work of Gordon et al. (1966). Velocity dependence is imparted through crossbridge distortion, whereby changes to sarcomere length alter average distortion of a population of cycling crossbridges (Razumova et al. 1999).
In the present work, we have embedded the fully-coupled myocyte electromechanics model with region-dependent parameters within a model of LV electrical activation and mechanics. Instances of the myocyte model located at each collocation point throughout the mesh take as inputs membrane potential, sarcomere length, and shortening velocity which originate from solutions to the corresponding monodomain and finite deformation partial differential equations. The myocyte model provides as outputs its contribution to changes in membrane potential and active tension.
Initial conditions for myocyte electromechanical model were determined by running it independently of the ventricular model under conditions of unloaded shortening (Rice et al. 2008) for 50 beats, at which point model responses had achieved steady-state. Pacing rate was 1 Hz.
Multi-scale models such as that used in the present study present an inherent challenge in that several time scales have to be taken under consideration. Stable numerical methods for coupling the various model components under mixed time- stepping schemes are critical for computational efficiency (Nickerson et al. 2005; Niederer & Smith, 2008). Niederer and coworkers (2008) recently developed a scheme which allows stable and efficient coupling of velocity-dependent contraction models and finite-element mechanics. In this 'update' scheme, the contraction model outputs an activation value each coarse time step rather than the final active tension. The mechanics solver then uses the activation value to re-compute a length- and velocity-dependent tension with each intermediate Newton iteration until the solution converges and the coarse time step is complete.
We have modified the update scheme of Niederer et al. (2008) to function with the Rice contraction model (2008), which is embedded within the myocyte electromechanics model. Instead of a generic activation level, the entire vector of contraction-associated state variables is output from the myocyte model. ODEs upon which tension directly depends are then re-integrated over the coarse time step for each Newton iteration using an explicit Euler scheme. Sarcomere length and velocity are updated with every re-integration. This scheme stably accommodates relatively large mechanics time steps with minimal loss of accuracy. A coarse time step of 0.5 ms was used in all simulations. ODEs were solved on a finer time step with a Radau solver (Hairer & Wanner, 1999), using a maximum time step 0f 0.5 ms.
Left ventricular (LV) geometry was approximated using a truncated ellipse of revolution (Costa et al. 1996). At zero pressure, focal length was 3.75 cm (Streeter & Hanna, 1973a), endocardial radial coordinate was 0.32, epicardial radial coordinate was 0.65, and longitudinal coordinate of the base was 120 degrees. This initial shape was subdivided 16 times radially and 36 times longitudinally to yield a 576 element mesh (Figure 1A). Interpolating functions were linear in the circumferential, longitudinal, and radial directions, giving rise to a total of 1887 degrees of freedom (Costa et al. 1996). The myocardium was modelled as a slightly-compressible, transversely isotropic material with respect to local fibre coordinates (Guccione et al. 1995). Fibre angles were varied transmurally in a piecewise-linear fashion determined by fit to angles measured in a previous study by our group (Mazhari et al. 2001). As these fibre angles were not directly reported in the original publication, they have been included in the Supplemental Data. Sarcomere length (SL) in the unloaded myocardium was assumed to vary transmurally according to detailed measurements in rat myocardium (Rodriguez et al. 1993), and offset to a mean SL of 1.95 in accordance with published data from canine hearts (Streeter & Hanna, 1973b; Rodriguez et al. 1992).
Boundary constraints were identical to those described in the original work of Costa et al. (1996). Nodes along the basal (μ = 120 degrees) and apical (μ = 5 degrees) boundaries were fixed in the longitudinal direction. The most epicardial node on the basal boundary was also constrained in the radial direction. Finally, the most endocardial of the basal nodes was fixed in the circumferential direction.
The monodomain equation for electrical propagation was solved using a collocation-Galerkin method along with an operator-splitting scheme to incorporate the myocyte electromechanics model (Rogers & McCulloch, 1994; Flaim et al. 2007). Electrical activation initiating contraction was assumed to occur simultaneously everywhere on the endocardial surface. Under this assumption, there are no circumferential gradients in membrane voltage. Accordingly, electrical propagation was solved on a two-dimensional mesh in rectangular Cartesian coordinates representing the cross-section of the mechanics mesh, with identical refinement (Figure 1, B and C). The voltage solution was interpolated using cubic Hermite basis functions, giving rise to a total of 5032 degrees of freedom. Transmural activation time and action potential duration were converged to within 1% at this level of refinement, based on comparison with an equivalent problem having twice as many (32) transmural elements. Error in the calculation of electrical propagation arising from the absence of circumferential curvature in the two-dimensional model was also assessed in preliminary calculations. Briefly, a three-dimensional ‘wedge’ was created by extracting a 30-degree circumferential portion of the original ellipsoid mesh. This wedge was refined four times in the circumferential direction, and converted to rectangular Cartesian coordinates. Action potential propagation was simulated in two- and three-dimensional meshes (with inputs from the mechanics model held constant) and compared. Differences in transmural activation times between the two models were small (<8%). The tissue conductivity coefficient was set at 2.4 mS cm−1 to give a transmural activation time of approximately 24 ms at the thickest part of the ventricular wall (Yan et al. 1998). Geometry of the electrical mesh was updated every 4 ms in accordance with mechanical deformation induced by contraction.
In the complete ventricular model, the myocyte model parameter set (epi, mid, or endo) was specified in an element-wise fashion according to myocardial region (Figure 1B). For purposes of comparison, an alternate distribution of cell subtypes was created in which elements were assigned parameters randomly (Figure 1C). This was done in such a way that the overall ventricular volume corresponding to each cell subtype was identical between the two distributions.
A three-element Windkessel model was used to provide pressure boundary conditions at the endocardial surface of the finite element mesh using a previouslyestablished coupling method (Kerckhoffs et al. 2007). A simplified system in which the Windkessel model was coupled to a time-varying elastance model was used to determine circulatory parameters which best fit averaged LV pressure and aortic flow traces measured in a set of dog experiments (experiments described below). The resulting model parameters are shown in Table 1. LV pressure from the ventricular model and arterial pressure computed from the circulatory model were used to determine timing of events in the cardiac cycle such as aortic valve opening and closure.
The unloaded mesh was inflated to an end-diastolic pressure (EDP) of 3.9 mmHg in the absence of active myofilament tension, to correspond with the average enddiastolic pressure reported in the study of Mazhari et al. (2001). LV volume increased from 30.2 ml in the unloaded state to 38.4 ml at EDP. After inflation, a stimulus current was applied to the endocardial surface to initiate action potential propagation and contraction. All simulations made use of this same protocol unless otherwise noted.
End-diastolic geometry served as the reference for calculation of fibre strains throughout the cycle. All strains reported below reflect averages of strains in five longitudinally-consecutive elements at constant wall depth. Except where specifically noted, these groups of five elements corresponded to an apical region of the myocardium, centered about the prolate spheroidal coordinate of μ = 40 degrees.
All animal studies were performed according to the National Institutes of Heath Guide for the Care and Use of Laboratory Animals. Surgical protocols were approved by the Animal Subjects Committee of the University of California, San Diego, which is accredited by the American Association for Accreditation of Laboratory Animal Care.
Six adult mongrel dogs (19–24 kg) were anesthetized with intravenous propofol (6mg/kg), intubated, and mechanically ventilated with isoflurane (0.5% – 2.5%), nitrous oxide (2 l/min), and medical oxygen (2 l/min) to maintain a surgical plane of anesthesia. Once data collection was complete, animals were removed from respiratory support and euthanized with a lethal dose of propofol. Left ventricular pressure was measured with an 8-Fr pigtail micromanometer catheter with lumen (Millar Instruments, Houston, TX). End-diastole was defined as coinciding with peak of the Rwave of the electrocardiogram (ECG). A Doppler flow probe (model #T208; Transonic Systems, Ithaca, NY) was positioned around the ascending aorta to define the onset and conclusion of ejection.
As described previously (Waldman et al. 1985; Ashikaga et al. 2004), implanted radioplaque marker positions were followed throughout the cardiac cycle. Three columns (4–6 beads) of gold (0.9 mm) beads were plunged into anterolateral wall to form a triad of material makers at 1/3 of the apex-base distance. Using a biplane X-ray system, cineradiographic images were acquired at sampling rate of 125 frames/sec with mechanical ventilation suspended at end expiration. LV and systemic pressures, aortic flow, and surface ECG were digitally acquired simultaneously with these images.
Local fibre architecture (Streeter & Hanna, 1973a) in the region of the transmural bead set was determined (Ashikaga et al. 2004). Mean fibre angles, α, were measured at 1 mm increments across the wall and interpolated with leastsquared best-fit line.
Continuous, non-homogenous 3-D systolic strain distributions were determined in a local Cartesian cardiac coordinate system (X1, X2, X3), defining the circumferential, longitudinal, and radial directions, respectively. Strain was computed with reference to the end-diastolic state in all cases. At each transmural depth the mean fibre angle, α, was used to rotate the cardiac strain tensor into fibre, cross-fibre, radial (Xf, Xcf, X3) coordinate system (Costa et al. 1999). Time courses of myofibre strain, Eff, during successive cardiac contractions as well as transmural distributions at end-systole were obtained for each animal.
Results from the baseline model, in which myocyte subtypes were heterogeneously arranged in layers (BASELINE), were compared with those from a model in which subtypes were homogeneously distributed at random throughout the mesh (HOM) (Fig. 1, B and C).
Global LV pump function during ejection against arterial loading generated by the three-element Windkessel model showed only minor differences between BASELINE and HOM cases. The largest difference was in ejection time, which was 6% shorter in the HOM case (104 for BASELINE vs. 98 ms in HOM). Maximum dP/dt, peak LV pressure, and stroke volume were essentially the same between the two cases (see Table 2). Total transmural delay of action potential propagation, measured in the thickest region of the wall, was also nearly the same (24.0 ms for BASELINE vs. 23.9 in HOM). The spatial pattern of electrical activation produced in the BASELINE model is shown in Figure 2.
The two cases displayed nearly identical transmural gradients of fibre and crossfibre strain at the time of aortic valve closure (Figure 3A). For both simulations, endsystolic fibre and cross-fibre strains were in good qualitative agreement with experiments (Mazhari et al. 2001) in terms of strain magnitudes and transmural trends (Figure 3). The largest negative fibre strain measured in the study of Mazhari et al. (2001) was −0.085 at 30% wall depth from the epicardium; the largest in simulations were −0.063 (BASELINE) and −0.066 (HOM), at wall depths of 26% and 15%, respectively.
Time courses of fibre strain were markedly different between BASELINE and HOM simulations during isovolumic contraction and early ejection, in contrast to their similarity at end-systole (Figure 4). The HOM model produced a time course of fibre strain during isovolumic contraction which had a smaller transmural range than the heterogeneous case (0.05 HOM vs. 0.07 BASELINE at aortic valve opening). Also, in the HOM case fibres at both epicardial and endocardial extremes were stretched prior to ejection; in the BASELINE case only endocardial fibres were stretched during the same period (Figure 4 A and B). Fibre shortening in the BASELINE model was greatest in the epicardial region throughout systole, while in the HOM model mid-myocardial fibre shortening tended to be greatest.
Since the only time courses of transmurally distributed fibre strains reported in the literature are for individual animals (Ashikaga et al. 2004; Ashikaga et al. 2007), here we use previously unpublished measurements from our laboratory in six anesthetized dogs for comparison with model-produced time courses (Figure 4). Fibre strain time courses of individual dogs are presented in the Supplemental Data. Large differences in fibre strain magnitude and transmural range were observed across animals (Figure 4C). When maximum and minimum fibre strains across the wall were examined at each time point and averaged across dogs, a general stretching of fibres (without regard to transmural location) can be observed around the time of aortic valve opening (Figure 4C). On an individual basis, maximum fibre strain at the time of aortic valve opening occurred at the epicardium in four of the six dogs, and was located at the endocardium in the other two (see Supplementary Data).
Additional simulations were performed to assess the sensitivity of global LV function and regional fibre strains to factors other than cell subtype distribution (Table 2). Effects of cell subtype composition were examined by using midmyocardial cell parameters throughout the mesh (column ALLMID in Table 2). Activation time was altered by doubling the conductivity, from 2.4 to 4.8 mS cm−1 (HIGHCOND). Fibre angles were altered in two separate simulations to match the largest and smallest total transmural gradients observed among individual dogs used in the study of Mazhari et al. (2001). In these simulations, the total transmural fibre angle ranges were 170 and 112 degrees (STEEPFIB and SHALLOWFIB, respectively), compared with a range of 145 degrees used for all other simulations. The transmural gradient in sarcomere length was removed and instead set to the average value of 1.95 µm (Rodriguez et al. 1992) throughout the mesh (NOSLGRAD). The EDP used to inflate the model prior to activation was increased by 50% (from 3.9 to 5.9 mmHg) to examine the effect of preload (HIGHEDP). This increase in preload resulted in a new end-diastolic LV volume of 42.7 ml. Finally, effects of afterload on modelled results was assessed by reducing initial aortic pressure by 25%, from 67.5 to 50.6 mmHg (LOWAOP).
Time courses of fibre strain were reduced to two moments in the cycle, aortic valve opening and closure, for purposes of comparison between simulations. At these times, minimum and maximum fibre strains and their corresponding wall depths were determined. Average transmural fibre strain and absolute range were also calculated. To account for possible variation of the results due to the longitudinal wall location at which strains were calculated, these fibre strain parameters were calculated at a basal location (μ = 90 degrees) in addition to the standard apical location described in the Methods (μ = 40 degrees). Thresholds for significance were defined relative to the BASELINE simulation for each measure included in the sensitivity analysis and are as follows: >0.02 change for fibre strain values (including strain range); >20 percentage-point change for percent wall depth of minimum/maximum strain; and >20% change in hemodynamic measures and activation time.
Results from these simulations are compared with those from the BASELINE simulation in Table 2. The largest changes to hemodynamic function were seen in the HIGHEDP and LOWAOP simulations, which showed 47% and 35% increases in stroke volume, respectively. HIGHEDP also produced a 14% increase in maximum dP/dt. For both HIGHEDP and LOWAOP, hemodynamic changes were not associated with significant change in fibre strains. The ALLMID simulation produced a 20% increase in ejection time as well as changes to fibre strains at the time of aortic valve opening. Greatest among these was a reduction in strain range from 0.07 to 0.02. Apart from the HOM and ALLMID simulations, none of the simulated conditions produced abovethreshold changes to strain metrics at aortic valve opening. At end-systole, only the SHALLOWFIB simulation produced significant changes to fibre strains.
All strain differences that were above-threshold when calculated at the apical location remained significant when a more basal wall location was used (lower half of Table 2). Several additional differences became significant under these conditions, including location in the wall of minimum and maximum fibre strains at end-systole for the STEEPFIB simulation.
We have studied the role of myocyte heterogeneity in LV function by coupling a model of length- and velocity-dependent EC-coupling dynamics with finite element models of LV action potential propagation and mechanics. EC-coupling model parameters were spatially varied according to previous work (Flaim et al. 2006; Campbell et al. 2008) to reflect measured regional differences in cell behavior. Results from the model in which epi, mid, and endo cells were arranged in their corresponding layers were compared with those obtained by homogeneously distributing cell subtypes throughout the LV wall. Differences in fibre strain time courses between the models were clearly evident during early systole, but essentially vanished by end-systole. There was no substantial difference in hemodynamic performance between the two models. Several additional simulations in which other factors such as fibre angle, activation time, preload, and afterload were perturbed revealed that transmural strain patterns early in systole were still most sensitive to the cell subtype distribution. Additionally, endsystolic fibre strains were seen to be sensitive to changes in fibre angles and the longitudinal location of strain measurement.
End-systolic fibre strains have been used as an index of regional ventricular function (Mazhari et al. 2001). We used such strains, measured under carefullycontrolled conditions in a previous study by our group (Mazhari et al. 2001), to validate end-systolic strains produced by the model (Figure 3). The lack of substantial differences in end-systolic fibre strains between BASELINE and HOM models was unexpected given large differences seen at other times throughout systole (Figure 3 and Figure 4). The similarity between the models indicates that end-systolic fibre strain was not sensitive to heterogeneous distribution of cell subtypes. That end-systolic fibre strain was also insensitive to the overall cell subtype composition is evident from results of the ALLMID simulation. The sensitivity analysis instead suggests that end-systolic fibre strains are highly dependent on a combination of fibre angle and measurement location. The latter implicates wall thickness or curvature as important determinants of end-systolic fibre strain.
Transmurally heterogeneous arrangement of myocyte subtypes resulted in large gradients of fibre strain around the beginning of ejection (Table 2, Figure 4B). These gradients are larger than those from simulations in which cell subtype assignment was perturbed (HOM, ALLMID, Table 2). While this difference is interesting, variability of experimentally-measured time courses prevents stringent validation of those produced by the BASELINE simulation. The fibre strain range measured in dogs for this study was sufficiently large to accommodate strains produced by either of the two models (Fig. 4).
The sensitivity analysis performed here may partially explain the large variability of experimental time courses. For instance, a combination of perturbations to cell subtype composition (ALLMID) and location of strain measurement (basal) reversed the transmural gradient of fibre strain at aortic valve opening in comparison with BASELINE (see Table 2).
Interestingly, a 27% reduction in transmural activation delay had essentially no effect on fibre strain patterns at either of the two cycle times reported (BASELINE vs. HIGHCOND, Table 2). The present results do not therefore support the notion that activation delay and cellular heterogeneity interact to affect fibre strains as previously proposed (Kerckhoffs et al. 2003a; Cordeiro et al. 2004).
Previous modelling studies have examined the interaction of activation delay with heterogeneous electromechanics (Solovyova et al. 2003; Nickerson et al. 2005). Solovyova and co-workers (Solovyova et al. 2003) modelled two muscle elements coupled in series, one (representing endocardial cells) having slower twitch kinetics than the other. They found that in the presence of an activation delay the heterogeneous duplex produced greater twitch force than either type of homogeneous pairing (both slow or both fast). We failed to observe any connection between heterogeneity and overall contractility. This discrepancy may be attributed to the fact that in our model, heterogeneous muscle fibres are essentially parallel to each other in the ventricular wall, not in series as in the work of Solovyova et al. (2003).
Nickerson et al. (2005) used a representation of LV mechanics similar to that used here to examine the effects of heterogeneous action potential duration on transmural fibre strain. In both their study and ours, stretching of endocardial fibres in late isovolumic contraction was observed. Nickerson et al. also noted a significant difference in fibre strain dispersion between homogeneous and heterogeneous cases during action potential repolarization. In this study we have referenced strains to hemodynamic events rather than action potential repolarization, which prevents direct comparison. However, differences in strain time course can be seen between BASELINE and HOM simulations as late as 80% into the systolic interval (Fig. 4), in spite of their close agreement at the moment of aortic valve closure (Fig. 3A). Thus, the results of both studies are fairly consistent.
Other studies have examined sensitivity of transmural strains to fibre angle distribution as we did in the SHALLOWFIB and STEEPFIB simulations. We found that the range of fibre strains was sensitive to the fibre distribution, which has been shown before by Bovendeerd et al. (1992) and Geerts et al. (2003). Our results are in agreement with the latter study in that a shallower transmural fibre angle distribution resulted in an increase in the coefficient of variation of sarcomere length shortening during ejection (Geerts et al. 2003) and range of end-systolic fibre strains in our study. The insensitivity of range of fibre strains at aortic valve opening to the gradient of resting sarcomere length (BASELINE vs. NOSLGRAD, Table 2) was in agreement with the findings of Kerckhoffs et al. (2003a).
The simplified model of LV geometry used in the present study is appropriate for the examination of transmural gradients in LV strain, but does not account for the activity of the right ventricle (RV). Because coupling with the RV might influence stresses and strains at the LV-RV junction and in the septum, our model applies only to the left ventricular free wall. In addition, the activation sequence used here assumes simultaneous activation of the entire endocardium, which is a simplification with regard to in-vivo measurements of myocardial activation (Scher et al. 1953; Durrer et al. 1970). Finally, we have considered only transmural variations in myocyte properties, whereas apex-base and LV-RV variations have also been identified (Litten et al. 1985; Carnes et al. 2004; Kondo et al. 2006; Krenz et al. 2007). The strong effects of heterogeneous myocyte properties and spatial arrangement on fibre strain time course shown in this study suggest that efforts to represent additional heterogeneities would be worthwhile.
Myocardial cell model parameters used here were chosen to reproduce behaviors observed in isolated cells, which can be quite different from that occurring in vivo. The consequences of this are perhaps most evident in the magnitude of peak end-systolic fibre shortening produced by the present models (~6%), which is substantially lower than the experimental mean observed in data presented in this study (~15%).
The cellular electromechanics model used here, with its accompanying parameter sets for the three myocyte sub-types is subject to its own set of limitations (see Campbell et al. 2008 for discussion). Parameters in the original contraction model of Rice et al. (2008) derive mainly from rat data. In previous work (Campbell et al. 2008), we modified a limited number of these parameters to reproduce transmural variations in time to peak unloaded shortening reported among canine myocytes (Cordeiro et al. 2004). The paucity of canine-specific muscle physiology data and its transmural variation currently prevents a more detailed validation. The high sensitivity of early-systolic fibre strains to cellular electromechanical function shown in this study suggests that future studies are needed in which the location of strain measurements are more carefully controlled and local fiber angle distributions are measured individually for each experimental subject.
The present results suggest that heterogeneous myocyte behavior exerts substantial effects on regional myocardial strains during early systole. At the same time, end-systolic fibre strains appear insensitive to heterogeneities. Instead, strains at end-systole were most sensitive to the combined effects of fibre angle and longitudinal measurement location. Finally, global LV function was found to be insensitive to perturbations of spatial myocyte subtype distribution.
Multi-scale computational models of ventricular electromechanics have the potential to help understand the integrative nature of myocardial function in vivo in terms of underlying cellular and EC coupling properties, however, such models must reflect regional variations in these properties with increasing accuracy in order to be effective.
This study has been supported by an American Heart Association Predoctoral Fellowship (to S.G.C.), Medtronic, UC Discovery Grant ITL06-10159 (to A.D.M.), the National Biomedical Computation Resource (NIH grant P41 RR08605) (to A.D.M), the National Science Foundation (BES-0506252) (to A.D.M) and NIH grant HL32583 (to J.H.O.). This investigation was conducted in part using a facility constructed with support from Research Facilities Improvement Program Grant Number C06 RR-017588-01 from the National Center for Research Resources, National Institutes of Health.
Andrew McCulloch and Jeffrey Omens are co-founders of and consultants to Insilicomed Inc., a licensee of UCSD-owned software used in this research. Insilicomed was not involved in this research. Lawrence Mulligan is an employee of Medtronic.