Myocardial contractility of the left ventricle (LV) plays an essential role in maintaining normal pump function. A recent ex vivo experimental study showed that cardiomyocyte force generation varies across the three myocardial layers of the LV wall. However, the in vivo distribution of myocardial contractile force is still unclear. The current study was designed to investigate the in vivo transmural distribution of myocardial contractility using a noninvasive computational approach. For this purpose, four cases with different transmural distributions of maximum isometric tension (T_{max}) and/or reference sarcomere length (l_{R}) were tested with animal-specific finite element (FE) models, in combination with magnetic resonance imaging (MRI), pressure catheterization, and numerical optimization. Results of the current study showed that the best fit with in vivo MRI-derived deformation was obtained when T_{max} assumed different values in the subendocardium, midmyocardium, and subepicardium with transmurally varying l_{R}. These results are consistent with recent ex vivo experimental studies, which showed that the midmyocardium produces more contractile force than the other transmural layers. The systolic strain calculated from the best-fit FE model was in good agreement with MRI data. Therefore, the proposed noninvasive approach has the capability to predict the transmural distribution of myocardial contractility. Moreover, FE models with a nonuniform distribution of myocardial contractility could provide a better representation of LV function and be used to investigate the effects of transmural changes due to heart disease.

## Introduction

Myocardial contractility depends on several factors, including action potential morphology, Ca^{2+} dynamics, and sarcomere length [1]. It has been shown in previous studies that these factors vary as a function of transmural location in the left ventricular (LV) wall [2,3]. Myocardial force generation, therefore, may also vary transmurally. In this regard, some previous ex vivo experimental studies have examined the difference in isometric force generation in permeabilized cardiomyocytes, but only in samples from the subepicardium and subendocardium of the LV [4,5]. However, a more recent study examined cardiomyocyte samples from all three transmural regions and showed that maximum isometric force generation exhibits transmural heterogeneity in healthy human hearts with the midmyocardium generating greater force than the subendocardium and subepicardium [6].

In addition to cellular experiments, finite element (FE) modeling has been used, in combination with magnetic resonance imaging (MRI) and ventricular pressure catheterization, to assess in vivo myocardial material properties in the LV. In these studies, numerical optimization was used to minimize the difference between the deformation field calculated from MRI and that predicted by the FE model, using the measured pressure as a boundary condition. Passive material parameters have been estimated in the LV for both healthy [7,8] and infarcted myocardium [9,10]. The maximum contractile force that can be generated within the remote and border zone regions of LVs with myocardial infarction was also assessed using this approach [11–13]. However, transmural variations in force generation were not examined in these studies.

The goal of the current study, therefore, was to estimate the in vivo contractile forces generated in the subendocardial, midmyocardial, and subepicardial regions of healthy porcine ventricles using a computational technique. Briefly, this was accomplished by using a combination of MRI, catheterization, and FE modeling. Properties were determined by using an optimization scheme to minimize the difference between in vivo strains and ventricular volume calculated from MRI and FE model predicted strains and volume. This work was motivated by recent ex vivo transmural measurements of cardiomyocytes [6].

## Method

### Experimental Measurements.

The animals used in this work received care in compliance with the protocols approved by the Institutional Animal Care and Use Committee at the University of Pennsylvania in accordance with the guidelines for humane care (National Institutes of Health Publication 85-23, revised 1996). In order to assess regional wall deformation in healthy adult pigs (*n* = 4; male; approximately 40 kg), 3D spatial modulation of magnetization (SPAMM) MRI was performed with simultaneous LV pressure measurements using a pressure transducer (Millar Instruments, Houston, TX) [10]. The endocardium and epicardium of the LV were contoured from the MR images at early-diastolic filling, end-diastole, and end-systole in order to generate the reference geometry and calculate LV cavity volume. In vivo systolic strain was calculated using a custom optical flow plugin for ImageJ to derive 3D displacement flow fields [14].

### FE Model.

Each animal-specific LV FE mesh (*n* = 4) was produced using trilinear hexahedral brick elements, which were fit to 3D geometric surfaces that were generated from the endocardial and epicardial contours (TrueGrid; XYZ Scientific, Inc., Livermore, CA). A myofiber distribution of −37 deg (at epicardium) to +83 deg (at endocardium), with respect to the circumferential direction, was assigned to the mesh [10]. The LV wall was evenly divided into three transmural layers: subepicardium, midmyocardium, and subendocardium (Fig. 1). The measured pressure was used as a boundary condition on the endocardial surface to simulate end-diastole and end-systole.

where $C$, $bf$, $bt$, and $bfs$ are diastolic myocardial material parameters, $E11$ is strain in fiber direction, $E22$ is cross-fiber in-plane strain, $E33$ is radial strain transverse to the fiber direction, and the remaining terms are shear strains. The diastolic material parameters were determined in a previous study for each animal and incorporated into the animal-specific FE models [10].

**T**

_{0}[1,16], which was defined by a time-varying elastance model. In order to represent end-systole,

**T**

_{0}was reduced to

*m*and

*b*are constants,

*E*

_{11}is strain in the fiber direction, Ca

_{0}is peak intracellular calcium concentration, $Tmax$ is maximum isometric tension achieved under maximal activation, and $ECa502$ is the length-dependent calcium sensitivity given by

where *B* is a constant, $(Ca0)max$ is maximum peak intracellular calcium concentration, $l0$ is the sarcomere length at which no active tension develops, and $lR$ is the reference sarcomere length in an unloaded state. The parameters *T*_{max} and $lR$ are discussed further in Sec. 2.3. However, all other material properties used for active myocardium are defined in Ref. [16]. A cross-fiber in-plane stress component equivalent to 40% of the fiber component was added based on previous studies of LV contraction [12]. Each FE simulation was conducted in two phases, where the first phase represented passive diastolic filling and the second phase represented active systolic contraction to end-systole. Both the passive and active material laws were implemented through a user-defined material subroutine in the FE solver ls-dyna (Livermore Software Technology Corporation, Livermore, CA).

### Optimization Procedure.

*T*

_{max}within the LV wall. The optimization process was performed using the software ls-opt (Livermore Software Technology Corporation, Livermore, CA) as previously described [9]. Briefly, the genetic algorithm (GA) technique was used to minimize the objective function $(\Phi )$, which was taken to be the sum of the squared error between experimentally measured data and FE predicted results, and was defined as

where $n$ is the strain point within the myocardium, $N$ is the total number of strain points, and $Eij$ and $V$ are the FE predicted end-systolic strain and LV cavity volume, respectively. The over bar variables represent in vivo measured values. A total of *N* = 756 points evenly distributed throughout the FE model were compared to the nearest LV points from the MRI-derived strain data. Specifically, the points were selected at the element centroids in a pattern of 3 transmural × 36 circumferential × 7 longitudinal points. Around the LV wall, this yielded 252 points per transmural layer. The search range for the parameter *T*_{max} was 40 kPa to 400 kPa based on previous studies [11,12,17,18]. In order to determine the transmural distribution of *T*_{max}, in combination with different transmural distributions of $lR$, four cases were tested with the aforementioned optimization procedure (Table 1). As the GA converges, the range over which it searches for optimal values was reduced. Once each optimization was fully converged, the final range of the parameter space, for each of the optimal values, was estimated from the parameters used to generate the 64 simulations that were run during the last iteration.

## Results

All of the optimizations converged to the final set of *T*_{max} parameters after 15 generations of the GA. The values of *T*_{max} optimized for each case are shown in Tables 2–5, along with their final parameter ranges. Notably, as compared to case 1 (Table 2) where both $lR$ and *T*_{max} were uniformly distributed, inclusion of transmural variations in either $lR$ (case 2; Table 3) or *T*_{max} (case 4; Table 5) led to a smaller $\Phi $ value in all four animals; the $\Phi $ values decreased by 13% and 15% in case 2 and case 4, respectively, on average. Inclusion of transmural variations in both $lR$ and *T*_{max} as in case 3 further decreased the averaged $\Phi $ value by 22% as compared to case 1 (Table 2 vs. Table 4). Interestingly, in both case 3 and case 4, most of the animals (three out of four) exhibited higher *T*_{max} in the midmyocardium as compared to the values in the subepicardium and subendocardium (Tables 4 and 5, Figs. 2 and 3). On average, *T*_{max} generated in the midmyocardium was more than 1.7 times of that produced in the other layers in these two cases.

The transmural distribution of systolic strain components (circumferential, longitudinal, and circumferential–longitudinal shear) in the mid-LV free wall was investigated more closely for case 1 and case 3, which exhibited the highest and lowest $\Phi $ value, respectively. The FE results in case 1 and case 3 exhibited a spatially varying distribution of systolic strain similar to the experimental data (Fig. 4). Results of case 3, however, showed better agreement with the experimental measurements as compared to those of case 1 for all of the three strain components assessed. The value of the objective function was also assessed within each transmural layer for case 1 and case 3. The value of $\Phi $ in the subepicardium, midmyocardium, and subendocardium was 17.33, 14.45, and 36.83 for case 1, respectively, and 16.5, 13.3, and 23.6 for case 3, respectively. The fit between the experimental and FE calculated results was improved in all three transmural layers for case 3. Most notably, the value of $\Phi $ was decreased by 35% in the subendocardial layer, which showed the most improvement in the fit.

In order to assess the influence of nonuniform contractility and reference sarcomere length on LV torsion, the twist angle at the apex of the FE models was calculated in each of the four cases (Fig. 5). These calculations were performed using the approach outlined in Ref. [20]. The results indicate that case 3 produced the largest twist angle, compared to the other cases. More specifically, when comparing case 2 and case 3, which have the same nonuniform distribution of $lR$ but different distributions of *T*_{max}, it can be seen that a uniform distribution of *T*_{max} produced a twist angle of 11.5 deg (case 2) compared to a nonuniform distribution, which produced a twist angle of 16.3 deg (case 3). Interestingly, when comparing case 4 and case 3, which both have nonuniform distributions of *T*_{max} but different distributions of $lR$, it can be seen that a uniform distribution of $lR$ produced a twist angle of 12.2 deg (case 4). This implies that the distribution of *T*_{max} has more of an effect on the twist angle than the distribution of $lR$. When both *T*_{max} and $lR$ were assumed to be uniform (case 1), LV twist angle was at a minimum with a value of only 9 deg.

## Discussion

The goal of the current study was to estimate the transmural distribution of in vivo myocardial contractile force using a combination of MRI, catheterization, FE modeling, and numerical optimization. For this purpose, the maximum isometric tension, *T*_{max}, was investigated particularly because it is an important determinant of cardiomyocyte contractility. A recent ex vivo study has shown that maximum isometric force generation exhibits transmural heterogeneity in healthy human hearts [6]. Consistent with these results, the current study showed that when the same level of *T*_{max} is used in all three transmural regions, the objective function values are larger than those with a nonuniform distribution (Tables 2–5). This is reinforced by the fact that case 1 showed larger differences relative to the experimental measurements in the reported systolic strain components. This indicates that the strain field in a model with nonuniform contractility distribution (case 3) is more representative of the deformation that occurs in vivo. In support of this, the model with transmural variations in *T*_{max} exhibited better agreement with the in vivo experimental measures in terms of systolic strains.

Recently, a sensitivity study was performed, which evaluated the influence of different transmural distributions of contractile force, and showed that the transmural distribution affects LV deformation via altered torsion [20]. It is thought that transmurally varying properties are necessary for maintaining healthy pump function, which is directly linked to torsion [2]. This may explain, at least in part, the better fit exhibited by models assuming a nonuniform contractility distribution across the LV wall. In particular, case 3 showed better agreement with experimental circumferential–longitudinal shear strains, which are linked to torsion. Moreover, in the current study, the fit was further improved by incorporating a nonuniform transmural distribution of the reference sarcomere length, which represents a more realistic morphology of the LV at an unloaded state.

In order to more directly assess changes in LV torsion, the apical twist angle was measured in each of the models. It was found that a nonuniform distribution of *T*_{max} produced the largest value of twist angle compared to the cases with a uniform distribution. This implies that transmurally varying contractile properties could be a key factor in maximizing LV torsion. Although residual stress was not included in the model, the nonuniform transmural distribution of the reference sarcomere length incorporates an important factor that influences systolic function. This, together with transmurally varying contractility, may serve to homogenize the force generation during LV systole for an optimal pump function. In addition, changes in the magnitude of the active stress in the cross-fiber direction also impacted the fitness of FE models to experimental data (data not shown). Again, in these cases, models assuming transmurally varying *T*_{max} and $lR$ provided the best fit.

Interestingly, the midmyocardium from most of the animals tested in the current study exhibited the highest contractile force, which was followed by the subepicardium and subendocardium, in terms of *T*_{max} (Figs. 2 and 3). These results show strong agreement with the recent ex vivo experimental results of Haynes et al. [6]. In addition, previous studies have reported that the epicardium contracts slightly more than the endocardium, which is consistent with our study results [4,5]. Finally, consistent with previous studies [21–23], the systolic strain components exhibited transmural variation.

It should be noted that the end-systolic volume calculated in each of the FE models was within 5% of the experimentally measured value from MRI. This implies that all the converged models matched the actual volume very closely. Therefore, when comparing the $\Phi $ values of all the cases, the differences are primarily attributed to the model fit to the experimental strain. For example, comparing the average $\Phi $ for case 1 (68.65) to that of case 3 (53.47) implies that the decrease in $\Phi $ is caused by a better fit to the experimental strain data when using a heterogeneous sarcomere length and isometric tension distribution. Although case 3 exhibited the best fit, there were some deviations from the experimental strains (Fig. 4). This could be explained by the fiber angle distribution, which was based on histology rather than diffusion tensor MRI data (since it was not available). A previous study, however, showed that the transmural fiber orientation had less impact on the deformation of the LV during systole, especially the torsional deformation, compared to the transmural distribution of contractility [20]. Alternatively, this difference may be attributed to the assumption of transverse isotropy and/or the boundary conditions used. These limitations will be addressed in future studies.

## Conclusion

In conclusion, the current study showed that a nonuniform transmural distribution of myocardial contractile force produced the best agreement between in vivo measured strain from MRI and that predicted by the FE models. More experiments are needed to confirm the in vivo results of the current study. Despite the limited statistical power due to the small sample size (*n* = 4), the results support recent experimental ex vivo measurements on cardiomyocytes [6], which showed that the midmyocardium generates the greatest force. Since heart disease has been associated with altered myocardial contractility in specific transmural regions [6], the incorporation of transmural variation of active properties, therefore, may provide a better representation of how disease alters LV function.

## Acknowledgment

This study was supported by the National Institutes of Health Grant Nos. R01 HL063954 (R. Gorman), R01 HL111090 (J. Burdick), R01 HL73021 (J. Gorman), and by a grant from the National Science Foundation CMMI-1538754 (J. Wenk).