Evaluation of a Novel Finite Element Model of Active Contraction in the Heart

Finite element (FE) modeling is becoming a widely used approach for the investigation of global heart function. In the present study, a novel model of cellular-level systolic contraction, which includes both length- and velocity-dependence, was implemented into a 3D non-linear FE code. To validate this new FE implementation, an optimization procedure was used to determine the contractile parameters, associated with sarcomeric function, by comparing FE-predicted pressure and strain to experimental measures collected with magnetic resonance imaging and catheterization in the ventricles of five healthy rats. The pressure-volume relationship generated by the FE models matched well with the experimental data. Additionally, the regional distribution of end-systolic strains and circumferential-longitudinal shear angle exhibited good agreement with experimental results overall, with the main deviation occurring in the septal region. Moreover, the FE model predicted a heterogeneous distribution of sarcomere re-lengthening after ventricular ejection, which is consistent with previous in vivo studies. In conclusion, the new FE active contraction model was able to predict the global performance and regional mechanical behaviors of the LV during the entire cardiac cycle. By including more accurate cellular-level mechanisms, this model could provide a better representation of the LV and enhance cardiac research related to both systolic and diastolic dysfunction.

Introduction

Cardiac muscle constitutes the histological foundation of the heart. The contraction of cardiac muscle cells (i.e., myocytes) generates force and propels blood out of the heart chambers into the circulatory system. In general, the active contraction of myocytes involves complex mechanisms. Briefly, upon activation, the influx of Ca²⁺ ions promotes their binding to the troponin C molecules, altering the shape of the troponin complex and exposing the binding sites on the actin monomers to myosin heads. Once attached to the binding sites, the myosin heads undergo conformational changes and generate force, leading to shortening of the sarcomeres, the contractile units of myocytes (Gordon et al., 2000; Kobayashi and Solaro, 2005). In vitro cellular experiments have shown that the force generated during myocyte contraction depends on the intracellular concentration of free Ca²⁺ and the length, as well as the shortening velocity, of sarcomeres (ter Keurs et al., 1980; Daniels et al., 1984; Keurs et al., 1988).

Computational modeling is an important approach for studying the function of the heart, especially the left ventricle (LV). Previously, several finite element (FE) models have been established to investigate the muscle contraction of the LV during cardiac systole (Guccione and McCulloch, 1993; Guccione et al., 1993; Hunter et al., 1998). Among those widely used FE models is the simple but computationally efficient time-varying “elastance” model, in which the active fiber stress is computed from a function of peak intracellular Ca²⁺ concentration, time and sarcomere length (Guccione et al., 1993; Kerckhoffs et al., 2007). This model has been used to successfully assess the effects of apical torsion on heart function (Trumble et al., 2011), as well as predict decreased contractility in the border zone myocardium of ovine and human hearts with myocardial infarction (Wenk et al., 2011, 2012). However, this modeling approach does not include the force-velocity relation of cardiac muscle contraction, and neglects to represent the processes occurring at the cellular level.

In addition to the global-level FE models, there are several cellular-level models of myocyte contraction (Trayanova and Rice, 2011). These models were developed based on the Huxley 2-state cross-bridge model (Huxley, 1957), and have been shown to reproduce the typical mechanical features of myocyte contraction under certain experimental conditions (Schneider et al., 2006; Campbell et al., 2008; Rice et al., 2008). Recently, a more flexible myocyte contraction model (called MyoSim) has been developed by Campbell (2014). MyoSim extends the Huxley model by incorporating Ca²⁺ activation, cooperative effects and the effects of interfilamentary movement, which implicitly depend on the length and velocity of sarcomeres during contraction (Campbell, 2014). Moreover, the use of cross-bridge distribution techniques allows MyoSim to predict the mechanical behavior of myocytes under a wider range of experimental conditions as compared to the deformation-based methods (Razumova et al., 1999; Rice et al., 2008). However, cellular-level models are unable to predict the global function of the heart. In this regard, a deformation-based myocyte contraction model (Rice et al., 2008) has been previously adapted into a FE model of the LV to investigate the mechanical properties of a mouse heart (Land et al., 2012). This suggests that the incorporation of an enhanced cellular-level contraction model, i.e., MyoSim, may further improve the behavior of a global-level FE model of ventricular function.

Along this line, the present work introduces a 3D FE implementation of cardiac muscle contraction that was developed based on MyoSim. This allows for the coupling of cellular-level mechanisms into a ventricle-level model. The new method was then validated using animal-specific FE models of the whole LV in rats. Specifically, the results of the models were fit to experimental measures of myocardial strain and ventricular hemodynamics, in order to determine the sarcomeric parameters that govern contraction.

Materials and Methods

Experimental Measurements

In order to assess regional wall deformation in the LVs of healthy rats, 3D cine displacement encoding with stimulated echoes (DENSE) cardiovascular magnetic resonance (CMR) imaging was performed on 5 female Sprague-Dawley rats (~6 months; Harlan, Indianapolis, IN, USA) using a 7T Bruker ClinScan system (Bruker, Ettlingen, Germany; Zhong et al., 2011; Haggerty et al., 2013; Zhang et al., 2017), followed by LV pressure measurements using a pressure transducer (SPR-903, Millar Instruments, Houston, TX, USA; Pacher et al., 2008). End systolic (ES) strains, relative to end diastole (ED), were calculated using the software DENSEanalysis (Spottiswoode et al., 2007). ES LV torsion, represented as the circumferential-longitudinal (CL) shear angle α_CL, was calculated as follows (Rüssel et al., 2009; Zhang et al., 2017):

$\begin{array}{ll}{\alpha }_{CL}=si{n}^{-1}\frac{2E_{cl}}{\sqrt{(1+2E_{cc})(1+2E_{ll})}}& (1)\end{array}$

In order to generate animal-specific FE models, at a minimally loaded state, the LV myocardium was contoured from the CMR images at early-diastolic filling. These contours where then converted into 3D geometric surfaces. All animal procedures were approved by the Institutional Animal Care and Use Committee at the University of Kentucky and were in agreement with the guidelines by the National Institutes of Health for the care and use of laboratory animals (NIH Publication 85–23, revised 1996).

Ventricular FE Model

Each animal-specific LV FE model was created based on the geometric surface data and size of the corresponding rat LV at early diastole (Figure 1). The FE mesh was produced by filling the myocardial wall with 8-node hexahedral brick elements incorporating a trilinear interpolation scheme (TrueGrid; XYZ Scientific, Inc., Livermore, CA, USA). The myocardium was evenly divided into 3 layers (i.e., epicardium, mid-myocardium, and endocardium; Figure 1; Zhang et al., 2015). The initial sarcomere lengths in each of the three layers, defined in the unloaded reference state of the model, were epi: 1,910 nm, mid: 1,850 nm, and endo: 1,780 nm (Guccione et al., 1993). The helical fiber angles in the epicardium, mid-myocardium, and endocardium were assigned as −60, 0, and 60°, respectively, relative to the circumferential direction. The boundary conditions of the LV were assigned to allow the base to expand/contract radially within the plane, with movement in the direction normal to the plane fully constrained. To simulate the passive filling and ejection of the LV, the volumetric flowrate into and out of the LV was estimated from CMR images (Figures 2A,B) and used to drive the entire cardiac cycle.

FIGURE 1

Figure 1. Representative animal-specific FE model of a rat LV. (A) Full view of the LV model with the four segments labeled, (B) Short axis view of a mid-ventricular slice, and (C) Long axis view of a longitudinal slice. Stars represent the points in the model where strain was compared to experimental measurements during the optimization.

FIGURE 2

Figure 2. Representative examples from a single rat of (A) the volumetric flowrate into and out of the LV used to drive the cardiac cycle in the FE model, (B) the LV volume (LVV) generated by the FE model, and (C) the calcium transient in the LV model. (D) The LV pressure-volume relationship generated by all of the FE models and experiments. LV pressure (LVP) data are mean ± standard error (SE); n = 5. LVV data shown in (D) were normalized to end diastolic volume (EDV).

The myocardium of the LV was assumed to be nearly incompressible, transversely isotropic, and hyperelastic. The passive stresses were derived from the following strain energy function:

$\begin{array}{l}W=\frac{C}{2}\{exp[b_f{E}_{11}^2+b_t\left(E_{22}^2+E_{33}^2+E_{23}^2+E_{32}^2\right)\\ +b_{fs}(E_{12}^2+E_{21}^2+E_{13}^2+E_{31}^2)]-1\}& (2)\end{array}$

where E₁₁ is fiber strain, E₂₂ is cross-fiber strain, E₃₃ is radial strain, and the remaining terms are shear strains (Guccione et al., 1993; Wenk et al., 2011). Values for the material constants b_f, b_t, and b_fs were chosen as 18.48, 3.58, and 1.627, respectively, based on previous studies (Guccione et al., 2001; Zhang et al., 2015). The material constant C was adjusted until the LV ED pressure matched the experimentally measured value for each rat. The value was found to be C = 0.262 ± 0.082 kPa (mean ± standard error).

The active material properties were derived from the 2-state MyoSim model of striated muscle (Campbell, 2014). Briefly, the myosin heads can switch between detached and attached states, and can form cross-bridges with different lengths ranging from −10 to 10 nm (Figure 3). Particularly, when a myosin head is bound to a directly opposed binding site, the length of the cross-bridge was assigned a value of x_ps (Figure 3). In the current model, the cross-bridges were classified, based on their lengths, into 21 bins with bin-width of 1 nm (i.e., x₁ to x_n, n = 21; Figure 3). 21 bins were used to maintain an appropriate balance between accuracy and computational efficiency. The population distribution of cross-bridges in the ith bin, i.e., the number of myosin heads attached to the thin filament with a cross-bridge length of x_i at time t, which is defined as A(x_i, t), was calculated with the following set of 21 coupled ordinary differential equations (ODE) using an explicit 4th-order Runge-Kutta method with adaptive time-step size:

$\begin{array}{ll}\frac{\partial A(x_i,t)}{\partial t}=k_1(x_i)D(t)-k_{-1}(x_i)A(x_i,t)(i=1,2,3,\dots ,21)& (3)\end{array}$

where k₁ (x_i) and k₋₁ (x_i) are the strain-dependent rate constants for the attachment and detachment transitions, respectively. D(t) is the number of myosin heads in the detached state and is determined by the following equation:

$\begin{array}{ll}D(t)=N(t)-N_{bound}& (4)\end{array}$

where N(t) is the total number of sites that are activated and available for myosin heads to attach to, and N_bound is the number of sites that have already been occupied by myosin heads and is equivalent to the total number of myosin heads in the attached state (i.e., $\sum _{i=1}^{21}A(x_i,t)$). To compute the total number of activated sites, N(t), a forward Euler method was used to solve the following equation:

$\begin{array}{l}\frac{dN\left(t\right)}{dt}=a_{on}\left[C{a}^{2+}\right]\left(N_{overlap}-N\left(t\right)\right)-a_{off}\left({N\left(t\right)-N}_{bound}\right)\\ +k_{plus}{N}_{bound}-k_{minus}\left(N_{overlap}-N\left(t\right)\right)& (5)\end{array}$

where a_on and a_off are rate constants, k_plus and k_minus are constants defining cooperative effects, and N_overlap is the maximum number of binding sites that heads could potentially interact with and was defined by the relative positions of the thick and thin filaments within each half of a sarcomere, as previously reported (Campbell, 2009). A calcium transient [Ca²⁺], which was based on experimentally recorded transients from intact rabbit hearts (Laurita and Singal, 2001), was temporally scaled according to each animal-specific heart rate and used in all 5 FE models in the present study (Figure 2C). The beginning of the transient was synchronized to the beginning of the isovolumic phase in each FE model. To include the effects of interfilamentary movement, linear interpolation was used to displace the population distributions A(x_i, t) by (1/2)Δx, yielding A_s(x_i, t), where Δx is the half-sarcomere length change between time steps (Huxley et al., 1994; Tajima et al., 1994; Campbell, 2014). Since Δx is highly influenced by the velocity of sarcomere shortening/lengthening, the incorporation of interfilamentary movement accounts for, at least in part, the velocity-dependence of force generation.

FIGURE 3

Figure 3. Schematic of the 2-state contraction model based on MyoSim. x₁ to x_n are different lengths of cross-bridges (n = 21 in the proposed FE model); x_ps is the cross-bridge length when the myosin head is bound to a directly opposed binding site; k₁ (x_i) and k₋₁ (x_i) are the strain-dependent rate constants for the attachment and detachment transitions, respectively.

Finally, the active stress along the fiber direction (T) produced by cross-bridges was calculated as follows:

$\begin{array}{ll}T(t)={{\displaystyle \sum }}_{i=1}^{21}\rho k_{cb}{A}_s(x_i,t)(x_i+x_{ps})& (6)\end{array}$

where ρ is the number of myosin heads in a hypothetical half-sarcomere with a cross-sectional area of 1 m², and k_cb is the stiffness of the cross-bridges. To ensure that the FE model generates realistic active stress throughout the entire cardiac cycle, the following conditions were applied within the model framework: (i) N(t) ≥ 0 at any time; (ii) N(t) ≥ N_bound at any time point; (iii) D(t) ≥ 0 at any time; (iv) active fiber stress T = 0 when the sarcomere length was shortened to 1.20 μm. These conditions are necessary in order to maintain a physiological response, i.e., not allowing N(t) to attain a negative value. A summary of parameters used for the active stress calculation is shown in Table 1, which are related to the MyoSim model (Campbell, 2014). In addition, stress components equivalent to 25% of the fiber stress were added to the two cross-fiber directions.

TABLE 1

Table 1. Parameters for active stress calculation with 2-state model.

Both the passive and active material laws were implemented as a user defined material subroutine in the explicit non-linear FE solver LS-DYNA (Livermore Software Technology Corporation, Livermore, CA, USA).

Optimization Procedure

In order to determine the active material parameters used for the 2-state contraction model, numerical optimization was performed with the software LS-OPT (Livermore Software Technology Corporation, Livermore, CA) as previously described (Wenk et al., 2012; Wang et al., 2016). In the current study, a hybrid technique was employed, which utilized global and local search algorithms. Specifically, simulated annealing was used as a global optimizer to locate the region of the parameter space with the highest probability of containing the best set of parameters. Once this region was identified, the sequential response surface method, which is a gradient-based method, was employed to find the local minimum by iteratively reducing the size of the parameter space until the optimal set of parameters was achieved. A total of 10 parameters were optimized within the ranges initially set according to the values used for unloaded twitch contraction (Table 2) (Campbell, 2014). The goal of the optimization was to minimize the objective function (Φ), which was taken to be the sum of the squared error between experimentally measured data and FE predicted results, and was defined as follows:

$\begin{array}{ll}\Phi ={{\displaystyle \sum }}_{n=1}^N{{\displaystyle \sum }}_{i,j=1,2,3}{(E_{ij,n}-{\overline{E}}_{ij,n})}^2+{{\displaystyle \sum }}_{m=1}^6{(\frac{P_m-{\overline{P}}_m}{{\overline{P}}_m})}^2& (7)\end{array}$

TABLE 2

Table 2. Optimization results.

The first term of the objective function represents the errors induced by ES strains, where n is the strain point within the myocardium, N is the total number of strain points (N ≥ 250 points evenly distributed throughout the mid-layer of the FE model), and E_ij and ${\overline{E}}_{ij}$ are FE-predicted and experimentally measured ES strains, respectively. The second term represents the errors due to LV pressures. In this study, 6 pressure points were compared, including pressure at the end of isovolumetric contraction (IVC), three points during systolic ejection (including the peak pressure), ES pressure, and the end of isovolumetric relaxation (IVR). P_m and ${\overline{P}}_m$ are the FE-predicted and experimentally measured LV pressures, respectively.

Single Element FE Model

In order to confirm that the optimized parameters, which were determined by fitting organ level data, result in reasonable in vivo cellular level function, a single element FE model was employed. Specifically, the model was subjected to boundary conditions and calcium levels that replicate the following cellular experiments (de Tombe and Stienen, 2007):

1. Maximum tension generation (T_max) at a fixed sarcomere length of 2,300 nm and maximal calcium concentration.

2. Force-calcium relationship, in terms of Calcium sensitivity (pCa₅₀) and Hill coefficient, at a fixed sarcomere length of 2,300 nm.

3. Maximum and minimum tension redevelopment (k_tr) at a sarcomere length of 2,300 nm with 20% length release for 20 ms followed by restretch to 2,300 nm.

Results

The optimizations for all 5 animal cases displayed good convergence, and the values of the 10 active parameters optimized for each case are shown in Table 2. The mean of squared errors between experimentally measured strains and FE predicted data, i.e., Φ/element (strain), is also shown for each case with values <0.1. Using optimized values for the active parameters, the FE models generated a LV pressure-volume (PV) loop that showed strong agreement with experimental results over the entire cardiac cycle, especially during systole (Figure 2D).

In order to assess the accuracy of the FE modeling approach for capturing regional variations in LV wall deformation, all 6 components of ES strain (i.e., E_rr, E_cc, E_ll, E_cl, E_rl, and E_cr) were analyzed within 4 wall segments (i.e., Anterior, Lateral, Posterior, and Septal) in the mid-ventricular region. In terms of the axial components of strain, the FE model predicted similar ES circumferential and radial (E_cc and E_rr) strain distributions in the majority of the mid-myocardium, when compared to the experimental measures, with significant (p < 0.05) deviation seen only in the septal and anterior segments, respectively (Figures 4A,B). The FE model, however, significantly (p < 0.05) underestimated systolic longitudinal strain (E_ll) in the lateral, posterior and septal segments of the mid-myocardium (Figure 4C). This could be due to the assumed longitudinal myofiber angle distribution. While there existed some significant (p < 0.05) differences between FE predicted values and experimental measures for ES shear strains E_cr and E_cl (only in the septal segment of the mid-myocardium), the shear strain E_rl predicted by the model was comparable with experimental data throughout the entire mid-myocardium (Figure 5). Moreover, the mid-ventricular torsion at ES, represented as CL shear angle, was similar to that estimated from experimental measures over the entire mid-myocardium (Figure 6).

FIGURE 4

Figure 4. End systolic axial strain components (A: E_cc, B: E_rr, and C: E_ll) averaged at the mid-LV in the mid-myocardial layer. Data are mean ± SE; n = 5; *p < 0.05 for comparisons between model-predicted and experiment-derived strains using paired two-tailed t-test.

FIGURE 5

Figure 5. End systolic shear strain components (A: E_cr, B: E_cl, and C: E_rl) averaged at the mid-LV in the mid-myocardial layer. Data are mean ± SE; n = 5; *p < 0.05 for comparisons between model-predicted and experiment-derived strains using paired two-tailed t-test.

FIGURE 6

Figure 6. End systolic circumferential-longitudinal (CL) shear angles averaged at the mid-LV in the mid-myocardial layer. Data are mean ± SE; n = 5. No significant differences between model-predicted and experiment-derived results using paired two-tailed t-test.

The FE model predicted sequential relaxation throughout the entire LV myocardium. Specifically, in contrast to the prompt sarcomere re-lengthening that occurred at the base, both the mid-ventricle and apex exhibited prolonged post-systolic shortening (PSS) and delayed re-lengthening and relaxation in the epicardial layer (Figure 7, arrows indicate the beginning of re-lengthening). This sequential re-lengthening and relaxation was mitigated in the mid-myocardium, and eventually reversed in the endocardium (data not shown). At the mid-ventricle, the onset of sarcomere re-lengthening was sequential in the 3 transmural layers. Namely, re-lengthening first occurred in the epicardium, followed by the mid-myocardium and then the endocardium (Figure 8, arrows indicate the beginning of re-lengthening).

FIGURE 7

Figure 7. Profiles of sarcomere length (SL) over 2 cardiac cycles obtained from 3 representative elements located at the base, mid-LV, and apex in the lateral epicardium. Data are mean (thick solid lines) ± SE (shaded area); n = 5. Note that the simulations start from the unloaded state. Arrows indicate the beginning of re-lengthening.

FIGURE 8

Figure 8. Profiles of sarcomere length (SL) over 2 cardiac cycles obtained from 3 representative elements located in the lateral epicardium (epi), mid-myocardium (mid), and endocardium (endo) at the mid-LV. Data are mean (thick solid lines) ± SE (shaded area); n = 5. Note that the simulations start from the unloaded state. Arrows indicate the beginning of re-lengthening.

The single element FE model results, which utilized the animal-specific parameters determined from the optimization of the five cases, are shown in Table 3. The maximum tension generation (T_max) was found to be 135.3 ± 2.3 kPa. The calcium sensitivity, indicated by the level of calcium at which 50% of the maximum tension is developed (pCa₅₀), was found to be 6.46 ± 0.024 and the Hill coefficient was 2.41 ± 0.057. The force-pCa curve for Case 1 is shown in Figure 9. The maximum tension redevelopment (k_tr−max) at saturated calcium was 96.1 ± 1.91 s⁻¹ and the minimum tension redevelopment (k_tr−min) at low calcium was 16.1 ± 1.80 s⁻¹.

TABLE 3

Table 3. Results of the single element simulations of maximum tension generation, calcium sensitivity, and tension redevelopment.

FIGURE 9

Figure 9. Force-pCa curve for Case 1. The square markers represent the values from the single element simulations, while the curve fit is the Hill equation with a value of pCa₅₀ = 6.53 and a Hill coefficient = 2.51.

Discussion

In the present study, a novel FE model of active contraction was developed by incorporating a 2-state myocyte contraction model (i.e., MyoSim) into animal-specific ventricle models. This new model includes both length- and velocity-dependence of force generation, and, therefore, has the potential to provide a more accurate representation of cardiac function at the global level, since it is embedded into the 3D FE framework.

Among the several improvements of the new FE contraction model is the inclusion of cooperative effects in dynamic coupling of binding sites and myosin heads. In real muscles, the status of the binding sites influences the status of the myosin heads, and vice versa. This dynamic coupling is further influenced by the status of near neighbor molecules, which is referred to as cooperative effects and is determined by the values of k_plus and k_minus. Incorporation of this molecular effect allows the proposed FE model to simulate cardiac muscles with altered Ca²⁺ sensitivity, for example, due to a disease state. This could be assessed by optimizing related parameters (i.e., a_on, a_off, k_plus and k_minus) using experimental data from a diseased animal model, which is beyond the scope of the current study. Moreover, a cross-bridge distribution approach is employed in the proposed FE model to solve the kinetics of binding sites and myosin heads. This approach provides greater control on the strain-dependence of myosin kinetics, enabling the proposed model to mimic the contractile properties of muscles under a wider range of conditions. In addition, the proposed model incorporates the effects of interfilamentary movement. According to Huxley (1957), the length of the cross-bridges changes when the filaments move. The interfilamentary movement, especially at a higher speed relative to the rate of cross-bridge cycling, can perturb the normal cross-bridge distributions, and thus impact the production of active force. This incorporation accounts, at least in part, for the velocity-dependency in the FE model of muscle contraction. Since all of the parameters can be adjusted to fit cellular-level experimental data, this FE contraction model is directly linked to cellular mechanical properties. This is critical, especially when the model is used to relate changes in the mechanical properties of cardiac muscle at the cellular level to variations in the global function of the LV.

In order to validate this new FE contraction model, parameters were optimized to match the experimental data from normal healthy rats. Overall, the proposed model demonstrated strong agreement with the experimental results in terms of LV global function (i.e., PV loop) and regional deformation (i.e., systolic strains and ventricular torsion), although there existed some deviations from experimental measures of strain, mainly in the septal region. Particularly, the notable underestimation of systolic E_ll may be due to the assumption of constant fiber angle distribution along the longitudinal direction and could be improved by using more realistic distributions of fiber angles (from DT-MRI for example). The proposed FE model was able to capture the heterogeneity of ventricular relaxation, which has been observed in experiments of intact ventricles in vivo (Sengupta et al., 2006; Ashikaga et al., 2007). This is important particularly because the observed delay of sarcomere re-lengthening in the epicardial apex has been associated with pressure drop during early diastole (Sengupta et al., 2006), and therefore may impact diastolic filling. As such, the proposed FE contraction model was able to match experimental PV measurements of not only the systolic phase, but also the diastolic filling phase. This capability could be beneficial in facilitating research related to diastolic dysfunction (Sugiura et al., 2012). Moreover, pilot studies using the time-varying elastance model, a commonly used FE model which does not include a force-velocity relation, failed to optimize the material properties using the same optimization procedure to match the experimental pressure measurements. In addition, these models exhibited synchronized relaxation throughout the entire ventricle. These observations lend further support for the improvement of this new FE model implementation in providing a more accurate representation of LV global function.

It should be noted that some of the optimized values for the active parameters varied from animal to animal, particularly for those used to calculate attachment/detachment rate constants. However, the sensitivity studies previously conducted by Campbell showed that the contraction model (MyoSim) was less sensitive to changes in the strain-dependence of myosin kinetics, as compared to changes in calcium sensitivity which is mediated predominantly by a_on (Campbell, 2014). In line with this, the present study showed similar calcium sensitivity in healthy rats with a narrow range of a_on values. Moreover, although there existed some differences between the optimized values in the present study and those used for unloaded twitch contraction, the maximum proportions of total available binding sites (N: 0.15~0.25) and bound myosin heads (N_bound: 0.13~0.23) were comparable to the results reported for unloaded twitch contraction (Campbell, 2014). In terms of the single element FE models, which were used to assess in vivo cellular-level function, the values for pCa₅₀ were found to be greater than those typically measured with in vitro permeabilized cell preparations. However, it has been shown in several studies that intact cells, which are evaluated at in vivo temperatures, exhibit greater calcium sensitivity (de Tombe and Stienen, 2007; Chung et al., 2016). Thus, the values found in the present study provide a reasonable representation of in vivo values. Additionally, the maximum tension redevelopment (k_tr−max) at saturated calcium was found to be greater than those measured with in vitro experiments. But, it has been shown (de Tombe and Stienen, 2007) that as the temperature of these experiments is increased, the value of k_tr−max also increases. Thus, it can be inferred from de Tombe and Stienen (2007) that the values found in the current study are a reasonable estimate of in vivo function. The value of the Hill coefficient was found to be lower than those measured from in vitro studies (Dobesh et al., 2002; de Tombe and Stienen, 2007), but was still in a range that indicates cooperativitly is in effect. Finally, the maximum tension generated in the current study is in agreement with values found in previous studies (Land et al., 2012).

One limitation of the proposed FE contraction model is the use of a simpler 2-state cross-bridge myosin scheme rather than a more complex 6-state model (Campbell, 2014). However, the 2-state scheme was chosen for a balance between computational efficiency and accuracy. Enhancements to the 2-state model will be pursued in future work. Moreover, the proposed model did not include the effects of cellular shortening on Ca²⁺ transient. But, this effect is relatively small, and the proposed model includes the shortening effects on the total number of binding sites activated by Ca²⁺ as a compensating method. When comparing the FE model strain to the experimental strain, only n = 5 samples were used. Finally, the FE models did not include the right ventricle (RV). This likely caused the deformation in the septal region of the FE models to deviate from the experimental measurements from CMR. To overcome this limitation, the RV will be included in future models.

In conclusion, the proposed FE contraction model successfully predicted both the global function and regional deformation of the LV. The capability of the proposed model to capture an important feature of ventricular relaxation makes it a powerful tool for the future investigation of diastolic dysfunction. Moreover, the incorporation of cellular-level mechanisms may enable the proposed model to assess how pharmaceutical treatments that target cellular function influence global ventricular function.

Author Contributions

XZ wrote the finite element code, performed the computational analysis, collected the experimental data, and wrote most of the manuscript. Z-QL analyzed the experimental data in terms of 3D surface generation and regional strain. KC developed the cellular-level code and helped implement it into the finite element framework. JW helped implement the finite element code, analyze both the experimental/computational results, write the manuscript, and developed the conceptual design of the work.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors wish to thank Hua Wang and Amir Nikou for their help with CMR scans. This study was supported by an award from the American Heart Association 14BGIA18850020 (JW), a grant from the National Science Foundation CMMI-1538754 (JW), grants from the National Institutes of Health R01 HL090749 (KC) and P30 GM110787 (LH), and a postdoctoral fellowship from the University of Kentucky Center for Computational Sciences (XZ).

References

Ashikaga, H., Coppola, B. A., Hopenfeld, B., Leifer, E. S., Mcveigh, E. R., and Omens, J. H. (2007). Transmural dispersion of myofiber mechanics: implications for electrical heterogeneity in vivo. J. Am. Coll. Cardiol. 49, 909–916. doi: 10.1016/j.jacc.2006.07.074

PubMed Abstract | CrossRef Full Text | Google Scholar

Campbell, K. S. (2009). Interactions between connected half-sarcomeres produce emergent mechanical behavior in a mathematical model of muscle. PLoS Comput. Biol. 5:e1000560. doi: 10.1371/journal.pcbi.1000560

PubMed Abstract | CrossRef Full Text | Google Scholar

Campbell, K. S. (2014). Dynamic coupling of regulated binding sites and cycling myosin heads in striated muscle. J. Gen. Physiol. 143, 387–399. doi: 10.1085/jgp.201311078

PubMed Abstract | CrossRef Full Text | Google Scholar

Campbell, S. G., Flaim, S. N., Leem, C. H., and Mcculloch, A. D. (2008). Mechanisms of transmurally varying myocyte electromechanics in an integrated computational model. Philos. Trans. A Math. Phys. Eng. Sci. 366, 3361–3380. doi: 10.1098/rsta.2008.0088

PubMed Abstract | CrossRef Full Text | Google Scholar

Chung, J. H., Biesiadecki, B. J., Ziolo, M. T., Davis, J. P., and Janssen, P. M. (2016). Myofilament calcium sensitivity: role in regulation of in vivo cardiac contraction and relaxation. Front. Physiol. 7:562. doi: 10.3389/fphys.2016.00562

PubMed Abstract | CrossRef Full Text | Google Scholar

Daniels, M., Noble, M. I., Ter Keurs, H. E., and Wohlfart, B. (1984). Velocity of sarcomere shortening in rat cardiac muscle: relationship to force, sarcomere length, calcium and time. J. Physiol. 355, 367–381. doi: 10.1113/jphysiol.1984.sp015424

PubMed Abstract | CrossRef Full Text | Google Scholar

de Tombe, P. P., and Stienen, G. J. (2007). Impact of temperature on cross-bridge cycling kinetics in rat myocardium. J. Physiol. 584, 591–600. doi: 10.1113/jphysiol.2007.138693

PubMed Abstract | CrossRef Full Text | Google Scholar

Dobesh, D. P., Konhilas, J. P., and De Tombe, P. P. (2002). Cooperative activation in cardiac muscle: impact of sarcomere length. Am. J. Physiol. Heart Circ. Physiol. 282, H1055–H1062. doi: 10.1152/ajpheart.00667.2001

PubMed Abstract | CrossRef Full Text | Google Scholar

Gordon, A. M., Homsher, E., and Regnier, M. (2000). Regulation of contraction in striated muscle. Physiol. Rev. 80, 853–924. doi: 10.1152/physrev.2000.80.2.853

PubMed Abstract | CrossRef Full Text | Google Scholar

Guccione, J. M., and McCulloch, A. D. (1993). Mechanics of active contraction in cardiac muscle: Part I–Constitutive relations for fiber stress that describe deactivation. J. Biomech. Eng. 115, 72–81. doi: 10.1115/1.2895473

PubMed Abstract | CrossRef Full Text | Google Scholar

Guccione, J. M., Moonly, S. M., Moustakidis, P., Costa, K. D., Moulton, M. J., Ratcliffe, M. B., et al. (2001). Mechanism underlying mechanical dysfunction in the border zone of left ventricular aneurysm: a finite element model study. Ann. Thorac. Surg. 71, 654–662. doi: 10.1016/S0003-4975(00)02338-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Guccione, J. M., Waldman, L. K., and McCulloch, A. D. (1993). Mechanics of active contraction in cardiac muscle: Part II–Cylindrical models of the systolic left ventricle. J. Biomech. Eng. 115, 82–90. doi: 10.1115/1.2895474

PubMed Abstract | CrossRef Full Text | Google Scholar

Haggerty, C. M., Kramer, S. P., Binkley, C. M., Powell, D. K., Mattingly, A. C., Charnigo, R., et al. (2013). Reproducibility of cine displacement encoding with stimulated echoes (DENSE) cardiovascular magnetic resonance for measuring left ventricular strains, torsion, and synchrony in mice. J. Cardiovasc. Magn. Reson. 15:71. doi: 10.1186/1532-429X-15-71

PubMed Abstract | CrossRef Full Text | Google Scholar

Hunter, P. J., McCulloch, A. D., and Ter Keurs, H. E. (1998). Modelling the mechanical properties of cardiac muscle. Prog. Biophys. Mol. Biol. 69, 289–331. doi: 10.1016/S0079-6107(98)00013-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Huxley, A. F. (1957). Muscle structure and theories of contraction. Prog. Biophys. Biophys. Chem. 7, 255–318.

PubMed Abstract | Google Scholar

Huxley, H. E., Stewart, A., Sosa, H., and Irving, T. (1994). X-ray diffraction measurements of the extensibility of actin and myosin filaments in contracting muscle. Biophys. J. 67, 2411–2421. doi: 10.1016/S0006-3495(94)80728-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Kerckhoffs, R. C., Neal, M. L., Gu, Q., Bassingthwaighte, J. B., Omens, J. H., and Mcculloch, A. D. (2007). Coupling of a 3D finite element model of cardiac ventricular mechanics to lumped systems models of the systemic and pulmonic circulation. Ann. Biomed. Eng. 35, 1–18. doi: 10.1007/s10439-006-9212-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Kobayashi, T., and Solaro, R. J. (2005). Calcium, thin filaments, and the integrative biology of cardiac contractility. Annu. Rev. Physiol. 67, 39–67. doi: 10.1146/annurev.physiol.67.040403.114025

PubMed Abstract | CrossRef Full Text | Google Scholar

Land, S., Niederer, S. A., Aronsen, J. M., Espe, E. K., Zhang, L., Louch, W. E., et al. (2012). An analysis of deformation-dependent electromechanical coupling in the mouse heart. J. Physiol. 590, 4553–4569. doi: 10.1113/jphysiol.2012.231928

PubMed Abstract | CrossRef Full Text | Google Scholar

Laurita, K. R., and Singal, A. (2001). Mapping action potentials and calcium transients simultaneously from the intact heart. Am. J. Physiol. Heart Circ. Physiol. 280, H2053–H2060. doi: 10.1152/ajpheart.2001.280.5.H2053

PubMed Abstract | CrossRef Full Text | Google Scholar

Pacher, P., Nagayama, T., Mukhopadhyay, P., Batkai, S., and Kass, D. A. (2008). Measurement of cardiac function using pressure-volume conductance catheter technique in mice and rats. Nat. Protoc. 3, 1422–1434. doi: 10.1038/nprot.2008.138

PubMed Abstract | CrossRef Full Text | Google Scholar

Razumova, M. V., Bukatina, A. E., and Campbell, K. B. (1999). Stiffness-distortion sarcomere model for muscle simulation. J. Appl. Physiol. 87, 1861–1876. doi: 10.1152/jappl.1999.87.5.1861

PubMed Abstract | CrossRef Full Text | Google Scholar

Rice, J. J., Wang, F., Bers, D. M., and De Tombe, P. P. (2008). Approximate model of cooperative activation and crossbridge cycling in cardiac muscle using ordinary differential equations. Biophys. J. 95, 2368–2390. doi: 10.1529/biophysj.107.119487

PubMed Abstract | CrossRef Full Text | Google Scholar

Rüssel, I. K., Tecelao, S. R., Kuijer, J. P., Heethaar, R. M., and Marcus, J. T. (2009). Comparison of 2D and 3D calculation of left ventricular torsion as circumferential-longitudinal shear angle using cardiovascular magnetic resonance tagging. J. Cardiovasc. Magn. Reson. 11:8. doi: 10.1186/1532-429X-11-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Schneider, N. S., Shimayoshi, T., Amano, A., and Matsuda, T. (2006). Mechanism of the Frank-Starling law–a simulation study with a novel cardiac muscle contraction model that includes titin and troponin I. J. Mol. Cell. Cardiol. 41, 522–536. doi: 10.1016/j.yjmcc.2006.06.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Sengupta, P. P., Khandheria, B. K., Korinek, J., Wang, J., Jahangir, A., Seward, J. B., et al. (2006). Apex-to-base dispersion in regional timing of left ventricular shortening and lengthening. J. Am. Coll. Cardiol. 47, 163–172. doi: 10.1016/j.jacc.2005.08.073

PubMed Abstract | CrossRef Full Text | Google Scholar

Spottiswoode, B. S., Zhong, X., Hess, A. T., Kramer, C. M., Meintjes, E. M., Mayosi, B. M., et al. (2007). Tracking myocardial motion from cine DENSE images using spatiotemporal phase unwrapping and temporal fitting. IEEE Trans. Med. Imaging 26, 15–30. doi: 10.1109/TMI.2006.884215

PubMed Abstract | CrossRef Full Text | Google Scholar

Sugiura, S., Washio, T., Hatano, A., Okada, J., Watanabe, H., and Hisada, T. (2012). Multi-scale simulations of cardiac electrophysiology and mechanics using the University of Tokyo heart simulator. Prog. Biophys. Mol. Biol. 110, 380–389. doi: 10.1016/j.pbiomolbio.2012.07.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Tajima, Y., Makino, K., Hanyuu, T., Wakabayashi, K., and Amemiya, Y. (1994). X-ray evidence for the elongation of thin and thick filaments during isometric contraction of a molluscan smooth muscle. J. Muscle Res. Cell Motil. 15, 659–671. doi: 10.1007/BF00121073

PubMed Abstract | CrossRef Full Text | Google Scholar

Keurs, H. E., Bucx, J. J., De Tombe, P. P., Backx, P., and Iwazumi, T. (1988). The effects of sarcomere length and Ca++ on force and velocity of shortening in cardiac muscle. Adv. Exp. Med. Biol. 226, 581–593.

PubMed Abstract | Google Scholar

ter Keurs, H. E., Rijnsburger, W. H., Van Heuningen, R., and Nagelsmit, M. J. (1980). Tension development and sarcomere length in rat cardiac trabeculae. Evidence of length-dependent activation. Circ. Res. 46, 703–714. doi: 10.1161/01.RES.46.5.703

PubMed Abstract | CrossRef Full Text | Google Scholar

Trayanova, N. A., and Rice, J. J. (2011). Cardiac electromechanical models: from cell to organ. Front. Physiol. 2:43. doi: 10.3389/fphys.2011.00043

PubMed Abstract | CrossRef Full Text | Google Scholar

Trumble, D. R., McGregor, W. E., Kerckhoffs, R. C., and Waldman, L. K. (2011). Cardiac assist with a twist: apical torsion as a means to improve failing heart function. J. Biomech. Eng. 133:101003. doi: 10.1115/1.4005169

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, H., Zhang, X., Dorsey, S. M., McGarvey, J. R., Campbell, K. S., Burdick, J. A., et al. (2016). Computational investigation of transmural differences in left ventricular contractility. J. Biomech. Eng. 138, 114501. doi: 10.1115/1.4034558

PubMed Abstract | CrossRef Full Text | Google Scholar

Wenk, J. F., Klepach, D., Lee, L. C., Zhang, Z., Ge, L., Tseng, E. E., et al. (2012). First evidence of depressed contractility in the border zone of a human myocardial infarction. Ann. Thorac. Surg. 93, 1188–1193. doi: 10.1016/j.athoracsur.2011.12.066

PubMed Abstract | CrossRef Full Text | Google Scholar

Wenk, J. F., Sun, K., Zhang, Z., Soleimani, M., Ge, L., Saloner, D., et al. (2011). Regional left ventricular myocardial contractility and stress in a finite element model of posterobasal myocardial infarction. J. Biomech. Eng. 133:044501. doi: 10.1115/1.4003438

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, X., Haynes, P., Campbell, K. S., and Wenk, J. F. (2015). Numerical evaluation of myofiber orientation and transmural contractile strength on left ventricular function. J. Biomech. Eng. 137:044502. doi: 10.1115/1.4028990

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, X., Liu, Z. Q., Singh, D., Wehner, G. J., Powell, D. K., Campbell, K. S., et al. (2017). Regional quantification of myocardial mechanics in rat using 3D cine DENSE cardiovascular magnetic resonance. NMR Biomed. 30:e3733. doi: 10.1002/nbm.3733

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhong, X., Gibberman, L. B., Spottiswoode, B. S., Gilliam, A. D., Meyer, C. H., French, B. A., et al. (2011). Comprehensive cardiovascular magnetic resonance of myocardial mechanics in mice using three-dimensional cine DENSE. J. Cardiovasc. Magn. Reson. 13:83. doi: 10.1186/1532-429X-13-83

PubMed Abstract | CrossRef Full Text | Google Scholar