All animal experiments were performed in accordance with national and local ethical guidelines, including the Guide for the Care and Use of Laboratory Animals, the Public Health Service Policy on Humane Care and Use of Laboratory Animals, the Animal Welfare Act, and an approved California Medical Innovations Institute IACUC protocol regarding the use of animals in research.

Two porcine subjects were used in this study: one normal and one with HF. The description of these animals and the creation of HF has been detailed previously (Choy et al., 2018). The ischemia resulted in decline of the animal's ejection fraction (EF) from 56% at the time of coronary artery occlusion to 32% when the animal was sacrificed, 16 weeks later. Measurements of in vivo left ventricular pressure and volume for each subject were recorded at the time of sacrifice (the incorporation of this data is discussed in section 2.5). Excised hearts were arrested in diastole with a saturated solution of potassium chloride and were fixed with buffered formalin (Carson-Millonig formulation).

After fixation, the ventricular cavities were filled with a silicone rubber compound (Polyvinylsiloxane, Microsonic Inc., Ambridge, PA) in order to maintain the geometry during imaging. The hearts were then placed in a plastic cylindrical container filled with a susceptibility-matched fluid (Fomblin, Solvay Solexis, West Deptford, NJ) and held in place using open-cell foam. Anatomical MRI was then performed (Magnetom Prisma 3T, Siemens, Erlangen, Germany) with the following parameters: T1-weighted imaging using a 3D Fast Low Angle SHot (FLASH) sequence (0.3 × 0.3 × 0.8 mm spatial resolution, echo time (TE)/repetition time (TR) = 3.15/12 ms, scan time: 1.5 h); and T2-weighted imaging using a 2D multi-slice Turbo Spin Echo (TSE) sequence (0.3 × 0.3 × 0.8 mm spatial resolution, TE/TR = 94/15,460 ms, scan time: 2 h).

DT-MRI was performed using a readout-segmented diffusion-weighted spin-echo sequence (Porter and Heidemann, 2009) with b-value = 1,000 s/mm² along 30 directions and one b-value = 0 s/mm² reference, TE/TR = 62/18,100 ms and 1.0 × 1.0 × 1.0 mm spatial resolution with 4–6 signal averages to improve signal-to-noise ratio (scan time: 8–12 h). Diffusion tensors were reconstructed from the diffusion-weighted images using linear regression and custom MATLAB (The MathWorks, Inc., Natick, Massachusetts, United States) routines.

Ex vivo MRI data sets were imported and processed in Simpleware ScanIP (Synopsys, Mountain View, USA). Detailed geometric segmentations of the biventricular structure were created along with segmentations of infarcted tissue in the HF subject. Segmentation relied on a combination of well-established techniques including region growing, level-set thresholding, and morphological smoothing (Vadakkumpadan et al., 2010; Setarehdan and Singh, 2012). Manual intervention was used only if needed to eliminate spurious features.

The full ventricular structure including all four valve openings of the heart was reconstructed from the T2-weighted MRI data sets. The segmented geometry and the cavity morphology are shown in Figures 1A,B, respectively. The segmented geometry was meshed with quadratic tetrahedral elements using Simpleware's built in FE meshing suite as shown in Figure 1C and the resultant mesh with the cavities enclosed is shown in Figure 1D. These meshes were imported into the Abaqus software environment (version 6.14, Dasssult Systèmes, Providence, RI, USA), which was chosen as the FE solver for this research. Since these geometries are constructed from ex vivo imaging, they provide the geometry in an unloaded state.

DT-MRI provides diffusion tensors for each voxel that were decomposed into eigenvalues and corresponding eigenvectors. Primary eigenvectors associated with the largest eigenvalue were identified as the orientation of the myofiber (Scollan et al., 1998; Kung et al., 2011). For practical purposes of computational modeling, a fully continuous 3D field representation of the local material coordinates, derived from diffusion tensors, is needed. This was achieved using a linear invariant interpolation method (Gahm et al., 2012) whereby the diffusion tensor is decomposed into invariants and orientations, which are each interpolated in turn and reconstructed into an interpolated tensor at the point of interest x ∈ ℝ³.

To ensure only voxels containing cardiac tissue (and not fat, air bubbles or voids) were included in the interpolation, the requirements that eigenvalues of each voxel be strictly positive, and that the fractional anisotropy (FA), an invariant of diffusion tensors commonly used for tissue thresholding (FA > 0.12), were imposed prior to analysis. This value for FA was found experimentally to be the lowest that would fully threshold out non-fibrous tissue and was reasonably different from values of FA for cardiac tissue (Helm et al., 2005; Kung et al., 2011). The inclination angle α_h, defined as the angle between the myofiber projected onto the longitudinal-circumferential tangent plane and the circumferential unit vector (Bovendeerd et al., 1992; Scollan et al., 1998; Toussaint et al., 2013), was quantified for each voxel. Results are presented regionally for the left ventricle (LV), partitioned into the 17 regions following the American Heart Association (AHA) guidelines (Cerqueira et al., 2002). This is a commonly reported quantification of myofiber orientation, which allows us to compare our results with literature findings as a source of validation.

For the HF pig, in vivo pressure and volume measurements were set as target values in the model calibration. For the normal pig, an in vivo volume measurement, and pressure derived from the healthy baseline of a larger in vivo data set (n = 5) (Choy et al., 2018), were similarly used. The complete in vivo measurements used for this study are presented in Supplementary Table S1. Measurements of in vivo strains were also recorded but deliberately excluded from the calibration process to serve as an independent metric to validate the model.

We introduced a closed loop circulatory model adapted from simple lumped parameter representations (Hoppensteadt and Peskin, 1992; Pilla et al., 2009) of different compartments in the cardiovascular system. The ventricular chambers are defined as fluid-filled cavities fully enclosed by the combined meshed faces of the tetrahedral elements on the cavity surface and the surface elements described in section 2.3.1 that close off the chamber. The coupling of the lumped circulatory system and mechanical function was performed in Abaqus. Details on the numerical underpinnings for this are provided in the Abaqus Theory Guide (2014); here, we provide a brief overview concerning the relation of pressure, volume, compliance, resistance and fluid exchange within a lumped system. The volume V_i and pressure P_i inside a fluid cavity chamber i are related in the following manner:

where κ_i is the compliance (inverse of stiffness) of the vasculature. For the LV and RV, the compliance is highly nonlinear, depending on both the strain state of the material and time (outlined in sections 2.6.1 and 2.6.2). Following Eq. (12), we define additional dimensionless compliance vasculature representing key components of the circulatory system. As the FE ventricular model (and not the circulatory circuit) is the central focus of this research, we sought to model a working circulatory system with as few components and assumptions as possible. This resulted in a circulatory system with three compliance vessels representing the systemic arteries (SA), systemic veins (SV) and the pulmonary circuit (P). The FE model of the heart is connected to this lumped representation as illustrated in Figure 3.

Unidirectional fluid exchanges governing the flow between compliance vessels are driven by the pressure gradients between these chambers as defined by:

In each simulated cardiac cycle, the ventricles contract, increasing the pressure in their chambers until it exceeds the pressure in the connected outflow chamber, driving the flow of blood in the circuit and simulating the physiological circulatory system. Since the total amount of volume in the circulatory system is constant, we have the following:

This facilitates a stable limit cycle when simulations occur over multiple cardiac cycles.

Resistance and compliance values were based on literature values (Santamore and Burkhoff, 1991; Hoppensteadt and Peskin, 1992; Watanabe et al., 2004; Pilla et al., 2009) and adapted to ensure realistic flow. The values for these parameters are presented in Supplementary Table S3 along with literature values for comparison.

Flows from the pulmonary circuit into the LV and from the systemic circuit into the RV are set to zero during the contractile stages of heart function (i.e., during isovolumetric contraction, ejection, and isovolumetric relation). Outflow from the LV and RV is always permissible but only occurs when the pressures in these cavities exceed the pressures in the chambers they are ejecting into. Combining these restraints with the imposed unidirectional flow of the circuit is sufficient to produce realistic fluid exchanges analogous to the physiological circulatory system.

The physiological heart does not experience rigid constraints on its motion. Boundary conditions are needed in computational models, however, to prevent rigid body motion and ensure the problem is mathematically well posed. To accomplish this without placing overly restrictive constraints on the model motion, we exploited the coupled degrees of freedom introduced by enclosing the ventricular cavities (section 2.3.1). The slave node of the truncated pulmonary trunk was fixed in all directions; this indirectly enforces a weighted average restraint on the nodes it is coupled with; i.e., the valve ring, such that the average deformation is zero. The valve is still able to expand and contract during the heart cycle, but this motion is “centered” on a fixed point in space.

In order to initiate the cardiac cycle, each compliance vessel was loaded with fluid until it experienced the physiological pressure at ED. Where possible, pressure catheterization readings from in vivo subject-specific measurements were used. The set of initial conditions used in this study is given in Supplementary Table S4.

The initial conditions outlined above are sufficient to initiate the dynamic beating of the computational heart model. From the initial conditions (i.e., ED), contraction is initiated, which increases the pressure in the ventricular chambers. These rise until LV and RV pressures exceed the pressures in C_SA and C_p, respectively. At this stage, fluid exchanges occur and the ventricles empty while the pressures and volumes in C_SA and C_p increase. Once LV and RV pressures drop below the pressures in C_SA and C_p respectively, the fluid exchanges stop, and the LV and RV pressures decrease with the decline in active tension. The entire duration of these “active contraction” processes is 480 ms, which compares well with literature values for the timing of isovolumetric contraction 66–90 ms (Sengupta et al., 2006), ejection 270–347 ms (Beyar and Sideman, 1984a; Sengupta et al., 2006) and isovolumetric relaxation 64–93 ms (Hanrath et al., 1980; Sengupta et al., 2006), respectively.

At the end of isovolumetric relaxation, all components in the circulatory system behave purely passively. Pressure differences between the compliance chambers in C_SV and C_p and the RV and LV, respectively, drive the passive filling of the ventricular chambers. The majority of the volume transferred in the passive filling stage occurs in the early portion of the step when the pressure difference is largest. Passive filling occurs in 300 ms, which is sufficient time for the cavities to inflate to the ED state and results in a heart rate of 77 bpm. Multiple steps of active contraction and passive filling can be simulated in a continuous sequence until convergent behavior over the cardiac cycle is reached.

Mass proportional Rayleigh damping (i.e., a viscous term introduced in the FE system of equations proportional to the mass matrix; Hughes, 2000) is introduced to dampen unrealistic oscillatory behavior of the low frequency modes. Physiologically, these would be eliminated by the surrounding soft connective tissue in the chest cavity. Similarly, isotropic time-dependent linear viscoelasticity is defined as part of the material constitutive behavior to damp out the high frequency response during active contraction. Whereas cardiac tissue is generally known to exhibit viscoelastic behavior, suitable experimental data on porcine cardiac viscoelasticity were not available. Hence, the model incorporates a small amount of viscoelasticity to eliminate unrealistic transient behavior, which is achieved using a Prony series formulation (Dill, 2006) within the Abaqus material definition (Abaqus Theory Guide, 2014).

In vivo strain echo data (TomTec 4D LV-Function, Version 4.6, Build 4.6.3.9, Unterschleißheim, Germany) of the endocardial surface was collected and excluded from calibration to serve as an independent data source to perform model validation. These in vivo strains are calculated by partitioning the endocardial surface into 16 segments, and measuring local deformation in longitudinal and circumferential directions (Pedrizzetti et al., 2014). An illustration of this 16-segment partition is shown in Supplementary Figures S1A,B.

These in vivo strain measurements reference ED as the initial configuration, and as such, provide relative change in length through a single cardiac cycle compared to the ED state. For purposes of validation, we select ES as the primary point of comparison for measuring heart deformation.

To provide comparable strain measures from our FE model, the endocardial nodes were partitioned into the same 16 segment division as the TomTec strain data. This resulted in 12 quadrilateral segments for the LV trunk and four triangular segments for the apical region. Control nodes placed at the corners and midpoints of the regions were identified. Similarly to the speckle-tracking imaging technique used to determine strain in the in vivo case, the nodal deformations of these control points were extracted at different time points in the cardiac cycle. By fitting cubic splines through these control points, longitudinal and circumferential measurements are created. The change in longitudinal and circumferential spline length provides a consistent strain measurement (i.e., engineering strain) analogous to the strain measurements provided by the in vivo TomTec strain measurements. The quadrilateral and triangular surfaces, the control nodes, and the fitted cubic splines are illustrated in Supplementary Figures S1C,D.

A visualization of the 3D heart structure from T2-weighted MRI data is given in Figure 4, along with our ventricular segmentation. Structural details such as wall thickness and trabecular morphology extracted from the MRI data conform to features reported for the porcine heart (Crick et al., 1998). Segmentations of the ventricles were meshed with roughly 85,000 quadratic tetrahedral elements–a refinement determined from a mesh convergence analysis. This analysis identified converged model behavior (i.e., stress, strain, and cavity expansion) over the entire cardiac cycle for mesh resolutions greater than approximately 50,000 elements.

Results for the inclinations angles α_h for both porcine LVs are presented graphically in Figure 5 for each AHA segment. This regional presentation of α_h, partitioned by AHA region and further subdivided by transverse wall depth, presents in general with narrow distributions (i.e., small variance). Inclination angles in the normal subject, excluding the apex, vary from 66.5 ± 16.6° on the endocardium to −37.4 ± 22.4° on the epicardium in a predominantly linear fashion. Similarly, inclination angles in the HF subject vary from 63.0 ± 18.3° on the endocardium to −43.4 ± 19.8° on the epicardium.

A local orthonormal coordinate system aligned with the myofiber, sheet plane, and sheet normal directions are interpolated to the centroid of each element in the FE mesh. Images revealing the geometry with and without myofiber orientations are presented in Figure 6. In addition to the characteristics presented quantitatively above, other qualitative features are as follows. Firstly, myofiber orientations are predominantly tangential with geometric surfaces. This can be seen in Figures 6B,D, by the abundance of myofibers protruding through cut surfaces relative to those seen protruding through natural physiological surfaces. Secondly, myofibers are closely aligned with papillary structure morphology, as can be seen in Figures 6C–E. Finally, the distribution in inclination angle in the LV, varying from positive on the endocardium through to negative on the epicardium, is easily identified from the global myofiber arrangement as shown in Figure 6F.

The calibrated material parameters fit the human shear data (Sommer et al., 2015) with an R² = 0.997. This excellent fit is illustrated in Supplementary Figure S2 wherein model response (solid lines) is plotted against experimental data (circles). The corresponding material parameters that produced this response are a = 1.05 kPa, b = 7.542, a_f = 3.465(kPa), b_f = 14.472, a_s = 0.481(kPa), b_s = 12.548, a_fs = 0.283(kPa), and b_fs = 3.088.

The subject-specific calibration resulted in suitable values for the passive material scaling parameters A and B, and active material parameters T_MAX and n_s. These values are presented in Table 1 along with the initial unloaded LV cavity volumes V₀, as these are fundamental to the resulting material parameters.

To illustrate the efficacy of the optimization routine, the passive filling curve resulting from optimizing is given in Figure 7A. These passive filling curves fit the analytical Klotz curves with an R² = 0.967 for the normal subject and an R² = 0.995 for the HF subject. The final calibrated EDV also matches closely with target values; i.e., 57.3 vs. 57.8 ml (normal) and 102.4 vs. 103.0 ml (HF). Sample long- and short-axis profiles of the ventricular structure in the unloaded (i.e., initial) and the ED configuration are presented in Figures 7B,C for the normal subject.

Cardiac function was simulated for six consecutive cardiac cycles (Figure 8) to ensure converged solutions were achieved. All values and results reported in the following sections correspond to the final and converged solution; i.e., these values will differ slightly from the calibrated targets that were acquired for only the first beat simulated, which is visible in Figure 8.

LV functional outputs compared reasonably between in vivo and in silico results. This is seen in SV, which was 32.2 vs. 30.9 ml in the normal subject for in silico and in vivo respectively, and 34.5 vs. 33.0 ml in the HF subject for in silico and in vivo respectively. Similarly, ejection fraction (EF) was 54.4 vs. 53.4% in the normal subject for in silico and in vivo respectively, and 35.0 vs. 32.0% in the HF subject for in silico and in vivo, respectively. These differences in SV and EF are considered minor when compared to the biological variations that occur from beat to beat.

Using the key parts of the cardiac cycle as presented in Figure 8, we assessed myofiber stress and strain values for each ventricle. These values were volumetrically averaged (i.e., normalized by element volume) to remove potential mesh artifacts. The mean volumetrically averaged myofiber stress is at its lowest at the end of relaxation (ER); it increases during passive filling, reaching a peak passive stress at end diastole (ED) and then rapidly increases during the systolic phase of the heart. By the start of ejection (SE), the myofiber stress is higher (order of magnitude greater than passive stress), which enables the continued contraction of the heart through ejection. While the myofiber stresses are high at end systole (ES), they continuously decline from this point, reaching the lowest values at ER, when the cycle repeats.

The mean volumetrically-averaged myofiber stresses in the LV and RV are presented in Table 2. Additionally, myofiber stress contours are presented in Figure 9 over long-axis cut planes of the ventricular structure. These contour plots reveal qualitative details about the myofiber stress distributions associated with geometric position; e.g., peak stresses are seen on the endocardial surface of the LV.

Since the HF subject had additional geometric information regarding the position and degree of pathological tissue (via the h field variable), we could further analyze stress by tissue health. Considering healthy tissue as regions whereby h = 1, infarcted tissue as h = 0, and border-zone tissue as values between these, we found mean myofiber stress was significantly different (p < 0.001) between all three regions, with substantially increased stress values within the infarcted tissue (Table 3).

The differences in stress values between normal and border-zone tissue are muted due to averaging. Analysis of these values as a function of geometric proximity to the fully infarcted/fibrotic tissue reveals more substantial differences (Figure 10). Most interestingly, the border-zone tissue experiences peak stresses roughly 1 mm away from the fibrotic tissue, after which they decrease and converge with healthy tissue values (distances >2 mm) for both ED and ES. Overlapping data points were analyzed for statistical significance and are illustrated in (Figure 10).

Similar to the analysis of myofiber stress, strain results were analyzed at the same key points in the cardiac cycle. The mean volumetrically averaged myofiber strains in the LV and RV are presented in Table 4. A detailed analysis of regional strain in the LV with respect to the local longitudinal and circumferential directions is presented in Section 3.5.

Initial comparison of the endocardial strains revealed strong agreement of global strains (i.e., averaged over the entire surface). For the normal subject, the global circumferential strain (GCS) for the FE simulations was −22.0%, which compares very well with the in vivo measurement of −21.9% from TomTec data. Global longitudinal strain (GLS) was −10.3% for the FE model simulation and −14.7% for the in vivo measurement. The HF subject produced similar comparisons: GCS for the FE simulations and the in vivo measurement were −14.4 and −12.7%, respectively, and GLS for the FE simulations and the in vivo measurement were −4.6 and −8.3%, respectively. In both subjects, circumferential strains were in closer agreement than longitudinal strains, which were both 4% lower than in vivo measurements from echocardiograms.

A regional analysis of comparison for FE model results and the in vivo measured values for strain on the endocardium surface reveals a very strong agreement in circumferential strains for both subjects, as seen in Figures 11B,D. Qualitatively, it is clear that the longitudinal strain behavior from the FE model correlates to the recorded in vivo values, with similar regional patterns displayed in both modalities, as seen in Figures 11A,C.

The geometric model construction used in this study is one of the most sophisticated ventricular structures produced for FE modeling of the heart compared to those found in the literature. While biventricular representations of the heart have become popular geometric choices for FE studies in recent years, most of these models truncate the geometry or exclude the papillary structures. The cardiac model from the Dassault Systèmes Living Heart Project (Baillargeon et al., 2014, 2015; Genet et al., 2016; Sack et al., 2016b) is an exception that does include this level of detail (and atrial structures); however, this heart geometry is not entirely patient-specific.

The geometric model of HF used in this study is a novel extension of how infarcted tissue has typically been presented in computational studies. Previously, studies have represented infarcted tissue using discrete concentric zones accounting for infarction and the border-zone as neat self-contained geometric regions (Wenk et al., 2011; Miller et al., 2013; Berberoglu et al., 2014). Our study incorporates subject-specific pathological detail, allowing for the description of infarcted material that ranges from concentrated to diffuse descriptions, as can be seen in Figure 2. Moreover, this level of subject-specific geometric detail is coupled with subject-specific fiber detail accounting for regionally precise myofiber detail throughout the ventricular structure, which has never been included in models with HF.

The critical role of myofiber orientation on mechanical and electrical function is well recognized (Beyar and Sideman, 1984b; Bovendeerd et al., 1992; Chen et al., 2005; Bishop et al., 2009; Wang et al., 2013). Fiber orientation, along with a second orthogonal vector within the laminar sheet, provides sufficient information to establish a local orthonormal coordinate system whereby each basis vector corresponds to a principal anisotropic direction. To the best of our knowledge, the present models are the first mechanical models to include subject-specific local coordinates derived from DTMRI.

The values of α_h in the normal subject conform well to other studies of large animals in terms of mean values (Streeter and Bassett, 1966; Streeter et al., 1969; Nielsen et al., 1991; Geerts et al., 2002; Helm et al., 2005; Ennis et al., 2008) and typical standard deviations expected from regional DTMRI myofiber analysis (Scollan et al., 1998; Lombaert et al., 2012). Our data reveals larger angles (i.e., close alignment with the longitudinal direction) in the anterior and inferior endocardial region when compared to other AHA regions circumferentially. These regions typically contain large papillary muscles on the endocardium, which explain this observation. These endocardial values of α_h in the anterior and inferior regions tend to plateau with values before declining linearly as a position of wall depth. The study of fiber orientation in papillary muscles is relatively unreported in the literature, due to difficulty in segmentation from in vivo imaging. In the early experimental work of Streeter et al. (1969), however, the authors also identified a plateau of myofiber angle α_h of roughly 90° in papillary muscles.

Myofiber orientations in the failing heart have striking similarities to the healthy heart in anterior, inferior, and the majority of septal (i.e., excluding basal anteroseptal) regions. Our results show that the largest discrepancies in myofiber angles correspond to regions in the LV free wall (i.e., basal and mid lateral regions), regardless of underlying fibrotic tissue content. This suggests that together with underlying pathology, mechanical factors relating to the position and function of an LV region are linked to its susceptibility to remodel.

Lower values of the myofiber inclination angle α_h in the vulnerable anterolateral regions are consistent with values reported in other studies investigating change in fiber orientation due to infarction (Holmes et al., 1997; Wu et al., 2007, 2009). In contrast, the inferolateral regions in the current HF model present with higher values of α_h, especially near the endocardial surface. Here, α_h plateau with values highly aligned longitudinally, indicative that the inferior papillary muscle was segmented within these regions and may be clouding the results. Even with consistent segmentation techniques, biological variance between subjects is an unavoidable challenge when sources of discrepancies are being determined. This likely explains the disagreement seen in regions with little fibrotic content.

The choice to use human shear experimental data, instead of porcine data from an earlier study (Dokos et al., 2002) was made after we calibrated our model to both and found that the porcine data produced much stiffer material behavior, especially in the fiber direction. As the human study was performed over a decade after the porcine study, more sophisticated methods were utilized in preventing contracture of heart muscles in the specimen extraction process, making the results likely more reliable representations of in vivo material response. Our calibration methods are proficient in capturing experimentally recorded orthotropic material response, as seen in Supplementary Figure S2, and when scaled, match the predicted Klotz curve for passive filling accurately (Figure 7). This two-stage method of cardiac tissue calibration, which utilizes small specimen data and pressure-volume data, is becoming a common method for obtaining realistic (and subject-specific) parameters (Krishnamurthy et al., 2013; Gao et al., 2015; Sack et al., 2016c). It ensures that the resulting material parameters produce a material that conforms to the anisotropy uncovered from mechanical experimentation and realistic in vivo function.

The determination of passive material parameters for the myocardial tissue ultimately depends on three factors: (1) the unloaded volumes V₀ of the ventricles; (2) the target ED volumes; and (3) the assumption regarding infarct stiffness. Our calibration techniques found that the HF subject yielded material parameters roughly an order of magnitude greater (considering the linear scaling coefficient A) (Table 1) than the normal subject, resulting in a much stiffer passive filling curve (Figure 7). This increase in remote stiffness is in line with other studies that have found that the remote tissue experiences changes in material properties (Bogaert et al., 2000; McGarvey et al., 2015), as well as the underlying physiological process of adverse remodeling, which results in collagen content, deposition, and cross-linking increasing in both infarcted and remote regions of the heart (Holmes et al., 2005; van den Borne et al., 2008; Fomovsky and Holmes, 2010). Another noteworthy finding is that the unloaded LV cavity in the subject with HF is almost three times as large as that of the healthy counterpart (Table 1).

Contraction is initiated by sarcomere shortening in series, which in turn contracts myofibrils. Whereas the mechanical analysis of this multi-scale phenomenon would seem to only occur in the myofiber direction, contraction has also been recorded in cross-fiber directions (Lin and Yin, 1998), which has been linked to the splay and dispersion of myofibers. The active material calibration resulted in T_max = 114.0 kPa and n_s = 0.07 in the normal subject, and T_max = 140.6 kPa and n_s = 0.14 in the HF subject. The value found for T_max is in line with the values reported in Genet et al. (2014): 130–155 kPa. These values of n_s are typically lower than those found in cardiac models that use generic fiber descriptions, which typically apply 40% of active stress in this direction (i.e. n_s = 0.4) (Lee et al., 2013; Genet et al., 2014; Zhang et al., 2015). Our lower values of n_s are due to the accurate subject-specific fiber orientations incorporated in the model, which naturally reduces the inclusion of unrealistically high (and potentially non-physiological) cross-fiber contraction.

The patient-specific metrics used to calibrate the material model of the heart are preserved in the converged cardiac-cycle simulation with an error <5%. This is hardly a shortcoming as a physiological in vivo heart operates within a range of values, often experiencing variances in SV, pressure and timing from cardiac cycle to cardiac cycle. Without any subject-specific pressure-volume data for the RV, we assumed the same material properties as the LV and loaded the RV using values derived from the literature. While it appears this initial loading state was far from the converged solution, it had comparatively little impact on the LV function. Rather, the RV experienced the majority of the functional change to ensure hemodynamic equilibrium. It is reassuring to see that despite the lack of appropriate data for the RV, the combined ventricular function converges to a state closely in line with LV in vivo targets. This unidirectional ventricular dependence is pronounced in the porcine heart, where the LV is overwhelming the major mechanical component of the heart, and may be less pronounced in human models.

Qualitatively, the strain results conform to expectations. The myofiber elongation is at its greatest at ED, is at its minimum during systole (due to the contractile material behavior), and returns to almost the original length before passive filling starts (Table 4). Furthermore, these strain results reveal functional changes due to pathology. The failing heart experiences diminished myofiber strains within the LV at ED (6.9 ± 4.3 vs. 9.6 ± 6.7) due to its stiffer material composition, and diminished myofiber strains at ES (−5.4 ± 6.5 vs. −9.8 ± 5.0) due to the loss of contractile function in the infarcted region.

An accurate determination of stress within complex mechanical problems is an unrivaled advantage of computational modeling. This is especially relevant for cardiac mechanics as changes in ventricular wall stress are thought to initiate pathological remodeling (Pfeffer and Braunwald, 1990; Sutton and Sharpe, 2000; Matiwala and Margulies, 2004). A unique problem pertaining to biological materials is that in vivo stresses cannot be accurately replicated under in vitro experimental protocols and thus, cannot be measured directly (Dorri et al., 2006). Our results show that chronic HF results in increased mean myofiber stresses in both ventricles at ED, and in the LV at ES (Table 2). This is interesting in the context of our strain results; i.e., the subject with HF experiences increased stress while simultaneously experiencing reduced strain. Determining stress directly from strain, without accounting for the anatomic features of the infarcted tissue or a proper constitutive characterization of the myocardial tissue (e.g., in Laplace's law), could lead to highly erroneous conclusions about the actual stress state within the failing heart. This has been demonstrated by Zhang et al. (2011), who compared Laplace's law and FE methods to evaluate stress in an infarcted LV. Their analysis showed that the average stress from Laplace's law is significantly different to the comprehensive stress analysis produced through anatomically accurate FE techniques.

While these increases in myofiber stress are prevalent throughout the myocardium, the infarcted region experiences myofiber stress roughly twice as high as in remote regions (Table 4). Furthermore, the myocardial wall containing the infarcted region, which is geometrically thinner than other regions, experiences complex stress due to its morphology and stiffer, non-contractile material behavior. Along the endocardium it experiences extremely high tensile stresses and along the epicardium it experiences compression (Figures 9A–C). We also discovered that in border-zone regions peak stresses occurred roughly 1 mm from fully infarcted tissue. Since stress is thought to initiate pathological remodeling (Pfeffer and Braunwald, 1990; Sutton and Sharpe, 2000; Matiwala and Margulies, 2004), this may explain the mechanical propagation of infarct expansion.

Due to the lack of comparable porcine computational studies we compared our stress results to available human studies. We found that the LV myofiber stresses in the normal subject compare well with literature values; e.g., our predicted ED stress of 2.1 ± 4.2 kPa vs. 1.47 ± 20.72 (Sack et al., 2016b) and 2.21 ± 0.58 (Genet et al., 2014).

The comparison between in vivo and model-predicted strains resulted in very similar qualitatively and quantitative results. For global values, the circumferential strains matched very well (<0.1% error), whereas the longitudinal strains were under-predicted in the FE model by merely 4%. A regional analysis reinforced this–circumferential strains matched extremely well for both subjects (Figure 11) and longitudinal strains matched with less accuracy (particularly in HF). The mismatch between longitudinal strain results may be due to simplifications inherent in the derivations of strains from echocardiography or our FE model. One the one hand, the surface resolution of echocardiography is typically low resolution (likely excluding finer morphological features), whereas the FE model has a much higher endocardial surface resolution. Echo-derived strains exclude the papillary muscles, and are often determined “sub-endocardially.” These papillary muscles, which pull downward during systole, would clearly alter the longitudinal strain measurement when included in the determination. On the other hand, our FE model excludes the atria and the pericardium. Ventricular motion is likely coupled with the presence (and weight) of the atria and their function (Fritz et al., 2014), which would specifically affect longitudinal strains. To further complicate the comparison, the difficulty in accurately measuring in vivo strains from echocardiography is well documented. For example, consistent strain-data reproducibility (intra- and inter-observer) is the subject of multiple medical studies, some showing very poor outcomes (Gayat et al., 2011; Badano et al., 2012). Currently, the reliability of global deformation metrics from echo-derived 3D strain is strong, (i.e., GCS and GLS), but the reliability to accurately quantify regional deformation is not (Lang et al., 2015). The uncertainty regarding the reliability of regional strains (from either modality) requires further study either using full-heart models or a better source of in vivo strain data (i.e. MRI-derived).

Reliable computational heart models offer the potential to integrate diverse data, produce otherwise unobtainable metrics relating to function, and quantify complex coupled mechanical behavior in an unparalleled manner. This can be especially advantageous as a research modality investigating the efficacy of treatments pre-clinically. Computational models allow for therapeutic parameters to easily be perturbed, whereby various in silico experiments can be investigated to determine optimal treatment efficacy.

As computing resources become cheaper and more efficient, computational modeling is increasingly being viewed as a viable complementing modality in the clinic. As shown in this study, reliable subject-specific models are achievable with the proper inclusion of high quality imaging data. As imaging technology becomes more advanced, it will become feasible to replicate the quality of models using purely in vivo imaging data. This would enable the ability to investigate potential therapies in silico prior to their application to patients–with medical decisions, treatment plans and interventions being highly patient-specific.

Additional techniques are needed to integrate clinical data into the methods presented here. Firstly, methods that can indirectly assess the unloaded geometry from an in vivo representation, such as inverse-displacement techniques (Bols et al., 2013; Rausch et al., 2017), will likely be required. Secondly, if pressure catheterization data is not available, methods to approximate or infer patient pressures will be needed. Cuff pressures could be used to identify the pressure range in the systemic arteries and Nagueh's formula (Nagueh et al., 1997), can be employed on echocardiogram data to obtain approximate LV filling pressures.

While geometrically detailed, our computational model is still lacking physiological features such as the atria, the pericardial sac and the diaphragm. This may be the source of discrepancies seen in longitudinal strains, especially among basal regions. For the normal subject, an average value of EDP had to be used when determining the in vivo target. A further shortcoming of this work is that it has only been applied to a single subject for each condition, but our current research is focused on expanding this to a larger cohort and including failing hearts treated with biomaterial injection therapy. Our material model of the heart is limited in that it does not include dispersion, micro-structural mechanics or regional heterogeneity of the tissue. Finally, only the mechanical component of heart function was simulated and excluded electro-chemical-mechanical considerations (i.e., the polarization of tissue, excitation and propagation phenomena).

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.