Relationship of Transmural Variations in Myofiber Contractility to Left Ventricular Ejection Fraction: Implications for Modeling Heart Failure Phenotype With Preserved Ejection Fraction

The pathophysiological mechanisms underlying preserved left ventricular (LV) ejection fraction (EF) in patients with heart failure and preserved ejection fraction (HFpEF) remain incompletely understood. We hypothesized that transmural variations in myofiber contractility with existence of subendocardial dysfunction and compensatory increased subepicardial contractility may underlie preservation of LVEF in patients with HFpEF. We quantified alterations in myocardial function in a mathematical model of the human LV that is based on the finite element method. The fiber-reinforced material formulation of the myocardium included passive and active properties. The passive material properties were determined such that the diastolic pressure-volume behavior of the LV was similar to that shown in published clinical studies of pressure-volume curves. To examine changes in active properties, we considered six scenarios: (1) normal properties throughout the LV wall; (2) decreased myocardial contractility in the subendocardium; (3) increased myocardial contractility in the subepicardium; (4) myocardial contractility decreased equally in all layers, (5) myocardial contractility decreased in the midmyocardium and subepicardium, (6) myocardial contractility decreased in the subepicardium. Our results indicate that decreased subendocardial contractility reduced LVEF from 53.2 to 40.5%. Increased contractility in the subepicardium recovered LVEF from 40.5 to 53.2%. Decreased contractility transmurally reduced LVEF and could not be recovered if subepicardial and midmyocardial contractility remained depressed. The computational results simulating the effects of transmural alterations in the ventricular tissue replicate the phenotypic patterns of LV dysfunction observed in clinical practice. In particular, data for LVEF, strain and displacement are consistent with previous clinical observations in patients with HFpEF, and substantiate the hypothesis that increased subepicardial contractility may compensate for subendocardial dysfunction and play a vital role in maintaining LVEF.


INTRODUCTION
Heart Failure (HF) is the only cardiovascular disease for which incidence, prevalence, morbidity, mortality, and costs are not decreasing. According to the 2017 Update (Benjamin et al., 2017), the prevalence of HF has increased from 5.7 million (2009 to 2012) to 6.5 million (2011 to 2014) in Americans >20 years of age and projections show prevalence will increase 46% by 2030, resulting in over 8 million adults with HF (Heidenreich et al., 2013). In 2012, the total cost for HF was estimated to be $31 billion and projections show that by 2030, the total cost will increase to $70 billion or roughly ∼$244 for every US adult (Heidenreich et al., 2013). Among patients hospitalized for an HF incident, 47% had HF with preserved ejection fraction (HFpEF) or systolic function, which is the focus of this paper.
The mechanism of the development of HFpEF is not well-understood (Aurigemma and Gaasch, 2004;Shah and Solomon, 2012;Steinberg et al., 2012;Sengupta and Marwick, 2018), and optimal treatment options remain unclear (Vasan et al., 1995;Bhuiyan and Maurer, 2011). Recent studies have suggested that HFpEF is associated with transmural changes in myocardial deformation (Shah and Solomon, 2012;Omar et al., 2016Omar et al., , 2017. Understanding the transmural variations in left ventricular (LV) mechanics associated with HFpEF may offer pathophysiological insights for developing potential therapeutic targets. We therefore explored a physics-based mathematical [finite element (FE)] model of the normal human LV to test the hypothesis that reduced subendocardial contractility combined with compensatory high subepicardial contractility may help in preserving LVEF independent of changes in myocardial geometry and material properties. We used our established computational framework in this paper. To the best of our knowledge, this is the first study that quantifies the development of HFpEF based on transmural variation in contractility, using patient-specific parameters.

Patient Data
In vivo echocardiographic recordings were obtained under a protocol approved by our institutional review board. Individual patients provided informed consent and anonymized data were sent to a core laboratory for analysis.

Geometry Considerations
The ventricle model pertains to a normal human subject. The LV was modeled as a truncated thick-walled ellipsoid (Mercier et al., 1982;LeGrice et al., 2001). Based on echocardiography recordings for end diastolic volume (EDV), LV diameter and wall thicknesses for the posterior and septal wall, we back-calculated ellipsoidal surfaces for the endocardium and epicardium at end diastole (ED).
Using a linearly regressed estimation of the unloaded LV cavity volume V0 (Klotz et al., 2006) we scaled the dimensions of the endocardium surface to match the calculated volume V0. The epicardium dimensions were then scaled to maintain the same myocardial wall volume ascertained at the ED configuration (preservation of mass).
TruGrid (XYZ Scientific Applications Inc, Pleasant Hill, California, USA) was used to mesh LV surfaces. The ventricle was meshed to produce eight layers through the radial direction (Figure 1). Finite element calculations were performed in ABAQUS (SIMULIA, Providence, RI, USA). The FE meshes are shown in Figure 1.
We used a rule-based approach coded in MATLAB 2012b (The MathWorks, Inc., Natick, Massachusetts, United States) to assign myofiber orientations to the centroid of each element in the meshed LV geometry. The aggregated myofiber orientation was assumed to present with an angle of −60 • from the local circumferential direction on the epicardium surface that varies linearly through the LV wall thickness to an angle of +60 • on the endocardial surface. This assumption is well-established in LV modeling studies (Carrick et al., 2012;Lee et al., 2013Lee et al., , 2015Genet et al., 2014), and based on histological studies (Streeter et al., 1969), and diffusion tensor MRI studies (Lombaert et al., 2011).

Constitutive Equation and Material Parameters
The material formulation of the LV tissue includes passive and active properties. The passive behavior of the tissue was described using the model introduced by Holzapfel and Ogden (Holzapfel and Ogden, 2009;Göktepe et al., 2011). Briefly, the strain energy function used to compute passive stresses is composed of FIGURE 1 | In normal conditions, contractility (T max ) was uniform in all layers (scenario 1). To simulate no contraction in the subendocardial region, contractility in three layers in white was set to zero (scenario 2). The three layers in red were used to simulate alterations in subepicardial contractility (scenario 3). The three white layers, the two green layers, and the three red layers comprise subendocardial, midmyocardial, and subepicardial regions, respectively. deviatoric ( dev ) and volumetric ( vol ) parts as follows: where a and b represent isotropic stiffness of the tissue, a f and b f represent tissue stiffness in the fiber direction, and a fs and b fs represent the stiffness resultant from connection between fiber and sheet directions; l 1 , l 4i , and l 8fs are invariants, defined as follows: where C is the right Cauchy-Green tensor, and f 0 and s 0 are vectors specifying the fiber and sheet directions, respectively. J is the deformation gradient invariant, and D is a multiple of the Bulk Modulus K (i.e., D = 2/K). The material constants a, a i , and a fs scale the strain-stress curve, whereas material constants b, b i , and b fs determine the shape of the strain-stress curve. To determine these parameters we used the End Diastolic Pressure Volume (ED PV) curve as described by Klotz et al. who reported an analytical expression for the ED PV curve based on a single PV point that is applicable for multiple species, including humans (Klotz et al., 2006). The LV EDV of 53 ml was recorded using echocardiography and the LV EDP of 14.3 mmHg was approximated from echocardiography data using Nagueh's formula (Nagueh et al., 1997).
The optimized material properties were found using an inhouse Python script that minimized the error between the ED PV curve from the FE model and the analytical expression (Klotz et al., 2006). The sequential least squares (SLSQP) algorithm (Jones et al., 2001) was used in the Python script, and ABAQUS was used for the FE modeling, as the forward solver (Table 1 and Figure 2).
The formulation for the active stress has been described extensively in the literature (Guccione and McCulloch, 1993;Walker et al., 2005;Genet et al., 2014;Sack et al., 2016). In short, the active stress in the myofiber direction was calculated as: TABLE 1 | Passive material properties that produced a pressure-volume curve close to the experimental pressure-volume curve ( Figure 2). FIGURE 2 | The passive material properties were determined such that the end diastolic pressure volume (ED PV) curve from finite element model was close to the experimental ED PV curve determined by Klotz et al. (2006).
where T max is the isometric tension at the largest sarcomere length and highest calcium concentration, Ca 0 is the peak intracellular calcium concentration, and 0 when t ≥ t 0 + t r t r = ml + b m, b = constants that govern the shape of the linear relaxation duration and sarcomere length relaxation. Also, where E f f is the Lagrangian strain in the fiber direction, B is a constant that governs the shape of the peak isometric tensionsarcomere length relation, l 0 is the sarcomere length that does not produce active stress, l R is the sarcomere length with the stressfree condition, and (Ca 0 ) max is the maximum peak intracellular calcium concentration. The active stress was added to the passive stress to compute total stress: where S is the total stress. The boundary and load conditions generally follow the ABAQUS Living Heart Model (Baillargeon et al., 2014(Baillargeon et al., , 2015Sack et al., 2016). In particular, the center of the LV proximal cross-section (base) was fixed. The average rotation and translation of nodes of the endocardial annulus were coupled to the center of the LV base. This boundary condition prevents rigid body rotation, but allows inflations and contractions of the annulus. The nodes of the base were fixed in the longitudinal direction. A pressure load was applied to the LV surface to simulate diastole, whereas the contraction of the LV muscles caused systole. Surface-based fluid cavities and fluid exchanges were used to model blood flow (ABAQUS Analysis User's Guide).
When T max is changed in Equation (2), the total contractile force of the tissue is altered, and other parameters related to the passive and active material formulations (Equations 1, 2) either do not change or change in a consistent way. We can prescribe different values of T max in transmural layers to introduce regionally varying contractility throughout the LV. We considered six scenarios with different contractile properties, as explained in Table 2. Homogenous contractile properties were considered in scenario 1, which also served to establish a baseline value for normal T max . T max was calibrated to produce the echocardiogram-recorded value for end-systolic volume (ESV) for this patient (24.8 ml). To simulate the diseased condition, subendocardial contractility was set to zero by setting T max = 0 (scenario 2). To recover ESV, a scenario was considered in which T max was increased in the subepicardial layers (scenario 3). To further assess the effects of transmural contractility, three more scenarios with different contractility in the transmural layers were created. In scenario 4, T max in all regions was reduced by 50%. In scenario 5, T max was set to zero in subepicardial and midmyocardial regions. In scenario 6, T max was set to zero in the subepicardial region.
Where τ is normalized LV torsion; Ø apex and Ø base are rotations in the apex and base, respectively; ρ apex and ρ base are the radius of the apex and base, respectively; and D is the distance between the apex and base (Figure 3).

RESULTS
The EF decreased from 53.2 to 40.5% when T max was set to zero in the subendocardial layers (Table 2 and Figure 4: scenario 2 vs. 1: 23.9% reduction in EF). The depressed contractility in the subendocardial region was enough to drop EF below 50%, producing HF with reduced EF (HFrEF). The EF normalized FIGURE 3 | The torsion of the LV was computed based on the apical and basal rotations, the apical and basal radius, and the distance between the apex and base. The formula used to compute the LV torsion (Equation 4) makes the LV torsion comparable for hearts of different sizes (Aelen et al., 1997;Rüssel et al., 2009). The positive rotation is counterclockwise when seen from apex. Scenario 1 represents the normal condition; scenario 2 represents zero subendocardial contractility; scenario 3 represents zero subendocardial contractility and increased subepicardial contractility (an HFpEF condition); scenario 4 represents decreased contractility in all regions; scenario 5 represents zero midmyocardial and subepicardial contractility, and scenario 6 represents zero subepicardial contractility. For scenario 1, the computational EF matched the experimental EF. For all scenarios, EDV = 53 ml, EDP = 14.3 mmHg. El, Ec, and Er are, respectively, ES strain in longitudinal, circumferential and radial directions. The circumferential and radial strains were computed using the nodes located at the endocardial annulus (base of the LV).
When subendocardial contractility was zero, LV torsion increased (scenario 2 vs. 1). The torsion further increased after contractility in a remaining region was increased to compensate (scenario 3 vs. 1 and 2). The torsion decreased when contractility in all transmural regions decreased by 50% (scenario 4 vs. 1). The torsion reversed when midmyocardial and subepicardial contractility were decreased to zero (scenario 5 vs. 1). The reversed torsion increased when only subepicardial contractility was zero (scenario 6 vs. 5).
Strains (which are independent of displacement boundary conditions) were altered in diseased conditions. The global longitudinal, circumferential, and radial strains decreased in HFpEF, but recovered after subepicardial contractility increased (Table 2, scenarios 2 and 3 vs. scenario 1). In addition, the global strains decreased when contractility decreased by half in all layers, and when subepicardial and midmyocardial contractility were zero, and also when subepicardial contractility was zero (Table 2, scenarios 4, 5, and 6 vs. scenario 1). The direction of circumferential strain changed when midmyocardial and subepicardial contractility were both zero (Scenario 5 vs. 1, Table 2). With normal homogenous contractility (scenario 1), FIGURE 5 | A long-axis view showing that at end systole, with uniform T max (scenario 1), all layers experienced compressive strain in myofiber directions. When subendocardial contractility was zero, the strain pattern was altered (scenarios 1 and 2), but it partially recovered when subepicardial contractility increased (scenarios 1 and 3).
all layers experienced contractile strains (Figures 5-7). Regional changes in contractility to simulate HFrEF (scenario 2) and HFpEF (scenario 3) both presented with tensile strains in the subendocardial regions where contractility was set to zero (Figures 5-7). However, the increased subepicardial contractility in HFpEF had a global effect on strains throughout all layers, reducing the strains in all regions. Qualitatively, the transmural strain curve of the HFpEF case (scenario 3) replicated the pathological HFrEF curve (scenario 2), albeit with strains that were 23.8% lower on average.
ES stress in the myofiber direction was noticeably reduced when subendocardial contractility decreased (scenarios 1 and 2, Figure 8). A trend to recovery in the stress distribution was observed when subepicardial contractility increased (scenarios 1 and 3, Figure 8).
The ES-shortening longitudinal displacement of the LV was profoundly decreased when subendocardial contractility was zero (scenarios 1 and 2, Figure 9). The longitudinal displacement was partially recovered when subepicardial contractility increased (scenarios 1 and 3, Figure 9).
The ES sphericity index (defined as the ratio between the lengths of the LV long axis and the short axis) was approximated as 1.1, 1.0, and 1.1 for scenarios 1, 2, and 3, respectively. In the HFrEF case (scenario 2), the ES sphericity index decreased compared to scenario 1. However, the ES sphericity index in the HFpEF scenario normalized toward the normal scenario. In other words, when the subendocardial contractility was zero, the LV shape became more spherical, compared to scenario 1. The FIGURE 6 | A short-axis view showing that the ES myofiber strain pattern altered when subendocardial contractility was zero (scenarios 1 and 2), but a partial recovery in strain pattern was observed when subepicardial contractility increased (scenarios 1 and 3).
FIGURE 7 | The ES myofiber strain at various points along LV thickness. In the horizontal axis, 0% represents the endocardium and100% represents the epicardium. The alterations in strains in scenario 2 are noticeable, compared to scenario 1. In scenario 3, the tensile strains decreased compared to scenario 2.
shape of the LV recovered toward the normal scenario when subepicardial contractility increased.

DISCUSSION
In this study, we used a realistic FE model of the human LV to examine the role of altered LV systolic mechanics as a mechanism FIGURE 8 | A long-axis view showing the ES myofiber compressive stress decreased when the subendocardial contractility was zero (scenarios 1 and 2). The stress pattern became partially similar to the normal case when subepicardial contractility increased (scenarios 1 and 3). of HFpEF. Our findings support the hypothesis that HFpEF could be a result of lower subendocardial contractility linked with increased subepicardial contractility (Sengupta and Narula, 2008;Shah and Solomon, 2012;Omar et al., 2016Omar et al., , 2017. When subendocardial contractility was zero, LVEF decreased by 23.9% ( Table 2: 53.2% in scenario 1 vs. 40.5% in scenario 2). The EF normalized when subepicardial contractility increased ( Table 2: 53.2% in scenario 3 vs. 53.2% in scenario 1). The change in subepicardial contractility (less than a 40% increase from normal values) resulted in a 31.4% improvement in EF. Unlike scenario 1, scenarios 2 and 3 experienced abnormal strains within the subendocardial region (Figures 5-7), even though scenario 3 experienced normal EF. The ES sphericity index decreased in scenario 2 (1.0) compared to scenario 1 (1.1), but it recovered in scenario 3 (1.1). The LV torsion increased in scenario 2 (26.7 • ) compared to scenario 1 (24.7 • ), and it further increased in scenario 3 (30.2 • ).
The subendocardial region played an important role in the LV systolic mechanics, as our results showed. In particular, when subendocardial contractility was zero, the EF was reduced below 50%. A scenario with EF below 50% and zero subendocardial contractility corresponds to HFrEF (Vasan et al., 1999;Owan and Redfield, 2005;Yancy et al., 2013). Also, reducing EF below 50% by zeroing subendocardial contractility is in line with previous studies that reported the important role of the subendocardial region in the mechanics of the LV (Sabbah et al., 1981;Algranati et al., 2011). Based on our adopted definition of end-systolic elastance (E ES ) (Chen et al., 2001), our results imply that E ES decreased after subendocardial contractility was zero, but E ES recovered when subepicardial contractility increased ( Table 2, scenarios 2 and 3 vs. scenario 1). Increased ESP in HFpEF could be due to alterations in the ejection period of the LV (scenario 2). After subendocardial contractility was lost, the ejection period shortened and ended with a higher pressure.
The EF decreased by 75% when contractility decreased by 50% in all layers (scenario 4 vs. scenario 1, Table 2). On the other hand, setting subepicardial and midmyocardial contractility to zero affected EF more than subendocardial contractility (0.3% in scenario 5 and 12.7% in scenario 6 vs. 40.5% in scenario 2, Table 2). This result illustrates the important role of subepicardial and midmyocardial regions and confirms previous experimental (Haynes et al., 2014) and computational (Wang et al., 2016) studies that indicated the important roles of the epicardium and midmyocardium in systolic mechanics of the LV. Also, based on our adopted definition of E ES , this parameter decreased when contractility decreased in all layers by 50%, and when subepicardial and midmyocardial contractility were zero, and also when subepicardial contractility was set to zero ( Table 2, scenarios 4, 5, and 6 vs. scenario 1).
Quantifying changes in torsional deformation related to changes in transmural contractility revealed an interesting relationship between the two. Abnormally high torsion could be a useful index of pathology, as we showed when subendocardial contractility was lost (Table 2, scenarios 2 and 3 vs. 1). This result confirms previous reports according to which the LV torsion increases in subendocardial ischemia (Prinzen et al., 1984), which has been related to the counter torque applied by the subendocardial region against the subepicardial region (Aelen et al., 1997). Also, this counter torque effect between subendocardium and subepicarium can be seen in scenarios 5 and 6 (Table 2). In these scenarios a negative torsion was seen after midmyocardium and subepicardium contractility was set to zero. The torsion of the LV is strongly coupled to the LV contractility and the inability to complete ejection properly ( Table 2).
The longitudinal strain has been reported as a criterion to diagnose normal and diseased hearts (Henein and Gibson, 1999;Takeda et al., 2001;Yu et al., 2002;Vinereanu et al., 2005). The decreased longitudinal strain in our results ( Table 2) corresponds to clinical studies that reported longitudinal strains decrease in HFpEF (Mizuguchi et al., 2010). Also, the contour of ES longitudinal displacement, which is directly related to longitudinal strain, was noticeably altered in the diseased scenario compared to the normal scenario (Figure 9, scenarios 1 and 2). However, when subepicardial contractility increased, the pattern of longitudinal displacement became more similar to the normal scenario (Figure 9, scenarios 1 and 3). The ES strain pattern across the regional layers of the LV wall (Figure 7) also supported the hypothesis that increased subepicardial contractility in HFpEF improves function globally (Sengupta and Narula, 2008). Moreover, the alterations in circumferential and radial strains are in line with clinical studies that reported these strains decrease in HFpEF (Wang et al., 2008;Mizuguchi et al., 2010). Yet, our results should be interpreted with caution. The circumferential and radial strains for normal conditions (scenario 1) were in line with clinical data reported in the literature, whereas the longitudinal strain was smaller than reported clinical data (Moore et al., 2000;Yingchoncharoen et al., 2013). The methodology of our study is similar to previous computational models of LV in our group. The longitudinal strain results of these previous models have been validated against experimental strain data (for example, Genet et al., 2014).
It has been well-documented that the shape of the LV changes in HF (Grossman et al., 1975;Carabello, 1995;Gaasch and Zile, 2011). In particular, LV concentric hypertrophy is seen in patients with HFpEF (Melenovsky et al., 2007), and exercise capacity is correlated with the sphericity index of the LV (Tischler et al., 1993). In line with previous studies, in our simulations, the shape of the LV was altered when the subendocardial contractility was zero. The ES sphericity index decreased in scenario 2 (1.0) compared to scenario 1 (1.1). When subepicardial contractility increased, the shape of the LV recovered toward the normal case, as seen in the ES sphericity index in scenario 3 (1.1) compared to scenario 2 (1.0). Thus, the increased subepicardial contractility may prevent LV dilatation in HFpEF, and help preserve the LV shape. This phenomenon will further support the normalization of LVEF due to the direct interplay between the LV shape and function (Grossman et al., 1975;Stokke et al., 2017).
In this study we used tissue-level load-independent properties (T max ) to alter myocardium contractility. This approach is more appropriate than using the LV strains. In fact, the popular notion of equating myocardial contractility with strain measurements (that are load dependent) is "off the mark [and] if contractility means anything, it is as an expression of the ability of a given piece of myocardium to generate tension and shortening under any loading conditions" (Reichek, 2013). Therefore, our approach to alter transmural contractility, which might not be feasible using current experimental methods, could lead to a better understanding of the development of HFpEF.
Novel physics-based mathematical modeling was used in this study to examine a possible mechanism underlying preservation of LVEF in HFpEF. The results from the simulations provide evidence of the potential role of myocardial contractility in the genesis of preserved EF in the HFpEF phenotype. Previous studies on HFpEF were mostly based on experimental data (for example, Phan et al., 2009), where the contribution of a single feature like myocardial contractility could not be varied in isolation of other parameters. However, we used FE modeling to simulate and isolate transmural contractility as a feature and study its effect on LV systolic mechanics. Our results provide important first steps toward eventual development of a computational model of HFpEF.

Study Limitations and Future Directions
The simulation addressed only the relationship between transmural myocardial contractility and LV systolic mechanics. Modeling all aspects of HFpEF was beyond the scope of this paper. In clinical conditions, several factors contribute to the development of HFpEF (Bench et al., 2009;Shah and Solomon, 2012;Sengupta and Marwick, 2018). These factors include abnormalities in both the systolic and the diastolic mechanics of the LV (MacIver and Townsend, 2008;MacIver, 2009;Shah and Solomon, 2012), the LV hypertrophy and geometric changes (Aurigemma et al., 1995;Vasan et al., 1999;Adeniran et al., 2015) and material properties of the LV (stiffness). A recent study employed computational models with heterogeneous transmural distributions of T max (Wang et al., 2016). In our investigation, only in one scenario (scenario 1) did we assume T max to be uniform in the transmural direction. Also, since this study focused on the LV, we assumed timing and activation of contractility are homogenous. A more realistic assumption would be to consider the sequence of electrical stimulations in the tissue (Chabiniok et al., 2012;Villongco et al., 2014;Crozier et al., 2016;Giffard-Roisin et al., 2017), particularly in the septum. However, a heterogeneous distribution of T max would be more important if atria were also included in the model. Moreover, we only modeled LV data from one human subject in our study. Modeling data from multiple subjects is the goal of a subsequent study. Here, our intent was to document our modeling methodology and demonstrate its utility.
Alterations in strain distributions might lead to remodeling in the LV tissue (Figures 5-7). The response of myocardial tissue to an altered mechanical environment will likely lead to changes in tissue properties that will in turn affect the LV inflation, contraction, and relaxation. It is well-documented that diastolic LV tissue stiffness becomes abnormally high in HFpEF (Zile et al., 2004). Our study focused on the systolic mechanics of the LV. As a future direction, integration of tissue response in diastole and systole will provide a more realistic and informative model to understand the mechanisms involved during the onset and development of HFpEF. Integration of cell-based crossbridge cycling and contractility could provide more realistic information about the tissue alterations over the course of HFpEF development (Adeniran et al., 2015;Shavik et al., 2017).
Although our study explored a simplified representation of HFpEF (appropriately so, to isolate mechanical effects), the clinical definition and diagnosis of HFpEF and HFrEF are more complex than just calculations of EF (Borlaug and Paulus, 2011). In fact, HFpEF lacks a clear validated diagnostic guideline (Lam, 2010;Oghlakian et al., 2011). In this hypothesis-generating study, we simply assumed EF < 50% represents HFrEF. This assumption is in line with some definitions used for HFrEF in the literature (Vasan et al., 1999;Paulus et al., 2007). However, an EF = 40.5% (scenario 2) might also be defined as borderline HFpEF (for example, Yancy et al., 2013). These points may be considered semantic because they do not affect the conclusions of our study, which quantified alterations in contractility with changes in EF, torsion and strain. It would be interesting to apply these methods to personalized models derived from patients diagnosed clinically with HFpEF and HFrEF.
Several other scenarios need to be investigated, including more graded loss of subendocardial contractility, and graded decrease of subendocardial contractility, with both coupled to a graded increase in subepicardial contractility. Moreover, the definitions of subendocardium, midmyocardium, and subepicardium regions were arbitrary in this study because exact definitions are not available. A more realistic imaging approach might better delineate transmural layers and their related contractility. Furthermore, exercise intolerance has been reported as a key factor in HFpEF (Roh et al., 2017), and could be implemented in our modeling methodology to better understand the mechanisms of HFpEF development. Despite these limitations, this paper reports instructive quantitative information about development of HFpEF, as we could change one aspect of the model (contractility at a particular location) and determine its effects alone.

CONCLUSIONS
The results of this study support the hypothesis that preservation of LVEF in patients with HFpEF could be explained on the basis of reduced subendocardial contractility with a compensatory increase in subepicardial contractility. These findings underscore the roles of regional LV myocardial contractility in HF syndromes and emphasize the importance of computational models in understanding pathophysiological mechanisms underlying complex phenotypic presentations like HFpEF.

AUTHOR CONTRIBUTIONS
PS and JG designed the study. KS developed the constitutive model. YD and KS created the computational models. YD ran simulations, compiled results, and wrote the initial draft of the paper. SS helped with the modeling process. YD, KS, JG, and PS contributed to analysis of the results, and manuscript writing.

FUNDING
This work was supported by NIH grants R01-HL-077921, R01-HL-118627, and U01-HL-119578. Further financial support was provided by the Oppenheimer Memorial Trust (OMT) and the National Research Foundation (NRF) of South Africa. Opinions expressed and conclusions arrived at are those of the authors and are not necessarily to be attributed to the NRF or OMT.