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 2+ 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 2+ 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 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 2+ concentration, time and sarcomere length 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(Wenk et al., , 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 2state 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 2+ 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 deformationbased 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 cellularlevel 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.

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): 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 . The helical fiber angles in the epicardium, midmyocardium, 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. 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: where E 11 is fiber strain, E 22 is cross-fiber strain, E 33 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 crossbridge 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 1 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 crossbridge 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: where k 1 (x i ) and k −1 (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: 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 x ps is the cross-bridge length when the myosin head is bound to a directly opposed binding site; k 1 x i and k −1 x i are the strain-dependent rate constants for the attachment and detachment transitions, respectively.
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., 21 i =1 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: 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 2+ , 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).
, 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. Finally, the active stress along the fiber direction (T) produced by cross-bridges was calculated as follows: where ρ is the number of myosin heads in a hypothetical halfsarcomere with a cross-sectional area of 1 m 2 , 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.
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: 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 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 P m are the FE-predicted and experimentally measured LV pressures, respectively.  Figure 1 in nm.

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 50 ) 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, 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 | 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.
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 midmyocardium), the shear strain E rl predicted by the model was comparable with experimental data throughout the entire midmyocardium ( 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 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).
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 50 ), 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 −1 and the minimum tension redevelopment (k tr−min ) at low calcium was 16.1 ± 1.80 s −1 .

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 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 | 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.  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 2+ 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 straindependence 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 cellularlevel function, the values for pCa 50 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 2+ 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 2+ 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.