We sought to calibrate mechanical properties of left ventricle (LV) based on three-dimensional (3D) speckle tracking echocardiographic imaging data recorded from 16 segments defined by American Heart Association (AHA). The in vivo data were used to create finite element (FE) LV and biventricular (BV) models. The orientation of the fibers in the LV model was rule based, but diffusion tensor magnetic resonance imaging (MRI) data were used for the fiber directions in the BV model. A nonlinear fiber-reinforced constitutive equation was used to describe the passive behavior of the myocardium, whereas the active tension was described by a model based on tissue contraction (Tmax). isight was used for optimization, which used abaqus as the forward solver (Simulia, Providence, RI). The calibration of passive properties based on the end diastolic pressure volume relation (EDPVR) curve resulted in relatively good agreement (mean error = −0.04 ml). The difference between the experimental and computational strains decreased after segmental strain metrics, rather than global metrics, were used for calibration: for the LV model, the mean difference reduced from 0.129 to 0.046 (circumferential) and from 0.076 to 0.059 (longitudinal); for the BV model, the mean difference nearly did not change in the circumferential direction (0.061) but reduced in the longitudinal direction from 0.076 to 0.055. The calibration of mechanical properties for myocardium can be improved using segmental strain metrics. The importance of realistic fiber orientation and geometry for modeling of the LV was shown.
Heart disease is the leading cause of death in the US and worldwide . Nearly 6.2 million adults had heart failure (HF) between 2013 and 2016 . The cost and prevalence of heart disease is on the rise. Between 2014 and 2015 (average annual), the estimated costs (direct and indirect) of heart diseases were $218.7 billion . The cost of HF was estimated to be $31 billion in 2012 in US alone and is projected to be $70 billion by 2030 . The number of people with HF in the US is projected to increase from more than 5 million in 2012 to more than 8 million in 2030 . Therefore, there is an absolute need for efficient methods for prevention and treatment of heart disease.
Characterization of myocardial material properties has important clinical applications [4–6]. Mechanical stress has a key role in understanding the mechanics of myocardium . Myocardium stress can be evaluated if the mechanical properties of the tissue are quantified. The stress distribution within myocardium affects blood flow in coronary arteries  and also oxygen consumption within the myocardial tissue . Abnormal myocardial stress leads to hypertrophy of the left ventricle (LV) . However, experimental invasive determination of myocardial stress is not currently a reliable option because errors inherent to the experimental methods are so large that they make the measurements unreliable [11,12].
Computational modeling has been used to estimate myocardial mechanical properties noninvasively. We previously used tagged magnetic resonance imaging (MRI) to calibrate mechanical properties . Recently, we used ex vivo MRI data to reconstruct an anatomically accurate geometry of the heart with high resolution . However, the material properties were estimated based on global LV metrics, namely stroke volume and longitudinal strain. The computational circumferential and longitudinal strains were generally in agreement with experimental strains obtained from three-dimensional (3D) speckle tracking echocardiography, but the level of agreement, particularly for longitudinal strain, needs improvement .
The goals of this study were threefold. First, to validate our calibration method based on isight (Simulia, Providence, RI) and abaqus (Simulia, Providence, RI) by testing this method for LV models in which calibrated material parameters are known. Second, to explore the capability of isight to improve state-of-the-art material calibration for a porcine biventricular (BV) model based on segmental strain metrics . Third, to compare the BV ex vivo MRI-based models  to the in vivo 3D echo-based LV models. Using an anatomically accurate geometry of the LV, we estimated LV mechanical properties to minimize the difference between the computational and experimental strain metrics. Unlike the previous studies that used global LV metrics such as stroke volume and global longitudinal strain, we used the segmental strain metrics from 16 segments according to American Heart Association (AHA) guidelines. The experimental strain metrics in these segments were used to define and minimize the objective function. Moreover, we used strains provided by speckle tracking echocardiography of the heart. The abaqus living heart framework [14–16] was used to find computational strains using finite element (FE) modeling, whereas the whole optimization algorithm was conducted by isight. One of the possible applications of this project would be to quantify myocardial contractility after treatments with cardiac devices such as the left ventricle assistive device (LVAD) or MitraClip. For example, after patients with HF receive an LVAD, the methodology presented in this study can be used to determine when their LV function has recovered sufficiently for explant of the LVAD [17–19]. By providing a more realistic approach to the calibration of the mechanical properties of the LV, the findings of this paper should lead to more realistic cardiac models.
The animal experiment was 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, and the Animal Welfare Act, and an approved California Medical Innovations Institute IACUC protocol regarding the use of animals in research.
One Yorkshire domestic swine was used in this study. The animal was sedated with TKX (Telazol 10 mg/kg, Ketamine 5 mg/kg, and Xylazine 5 mg/kg, intramuscular (IM)), intubated, and maintained on surgical anesthesia with isoflurane (1–2%) and oxygen while on the ventilator. The neck and the right inguinal area were shaved and scrubbed with nolvasan, alcohol, and betadine. A 5Fr introducer sheath was placed percutaneously in the right jugular vein to administer fluids and drugs. Heparin 100 international units (IU)/kg body weight was administered intravenously before further instrumentation. A 6Fr introducer sheath was placed percutaneously into the right femoral artery. A hockey-stick guiding catheter was advanced over a 0.035 in. guide wire through the introducer sheath toward the LV to measure pressure. The geometry of the LV was measured using 3D echocardiography using a Philips (iE33) transesophageal probe (X7-2t). Strain was measured using speckle tracking (TomTec, Unterschleissheim, Germany) at 16 segments of the AHA segmentation. The LV pressure was recorded at the time of echo acquisition. After the procedures, the animal was euthanized with a bolus injection (60 mL) of potassium chloride under deep anesthesia (Isoflurane 5%).
The geometry of the LV in early diastole was reconstructed from echocardiography images using QLAB 10 (Philips, Andover, MA) after the echocardiography data were acquired . The reconstructed geometry was meshed using TrueGrid (XYZ Scientific Applications, Inc., Pleasant Hill, CA). An example echocardiography image and the mesh are shown in Fig. 1.
In addition to a single LV model, we analyzed our previously described BV model . The LV and BV models are from the same animal. We used the Simulia Living Heart Project framework for the single LV model [14,15]. In the BV model, atria were excluded but their effects were considered via lumped parameter models .
During diastole, the end-diastolic pressure (EDP) was applied. During systole, the myocardium contracted, and the interaction between the LV and the circulation system was considered using a Windkessel model [21,22]. Systolic blood flow between the LV and the circulation system was modeled by applying an initial pressure of 80 mmHg, and unidirectional fluid exchange was enforced.
It was assumed that the myofibers deviated from the circumferential direction smoothly, changing in direction from −60 deg in epicardium, aligned circumferentially at the midwall, and orientated to +60 deg in endocardium [25,26]. In line with experimental [27,28] and computational studies [24,29], we assumed sheet direction is normal to the epicardial and endocardial surfaces where the fiber directions are located.
To determine the passive material constants, we used the analytical end-diastolic pressure volume relation (EDPVR) curve reported by Klotz and colleagues . For this purpose, we used end-diastolic volume (EDV) (57.8 ml) and EDP (13.5 mmHg) to create the analytical curve, and then, determined the passive constants such that the difference between the numerical and analytical volumes was minimized at multiple points along the EDPVR curve (see Results section, Fig. 2).
The center of aortic endocardial annulus was fixed, and the average translational and rotational motions of the nodes at the endocardial annulus were coupled to the fixed point. The blood pressure was applied at the endocardial surface as the loading condition. The optimization was done in isight (Simulia, Providence, RI), which used abaqus Explicit (Simulia, Providence, RI) as the forward solver (Fig. 3). The element types in the LV and BV models were C3D8 (eight-node linear brick) and C3D10M (ten-node modified tetrahedron, with hourglass control), respectively (abaqus Analysis User's manual).
The calculation of computational strains for the BV model has been described in another paper from our group .
The reference configuration at which the model was initialized was assumed to be at the start of filling, when the ventricle cavity is both the smallest and the least loaded. The echo-based strain measurements were most accurate at the endocardium, especially the values in the longitudinal and in the circumferential directions [34,35]. These measurements were obtained as relative strains between two time points in the cardiac cycle and reported as nominal strains (relative changes in length). Since neither of the two time points corresponded to the reference configuration, to obtain the strain values equivalent to experimental strains, a tool had to be developed. The two measuring points in the FE model correspond to the end of diastole (ED) and the point right after the end of ejection, i.e., end systole (ES).
The single LV cavity model was coated with membrane elements that share their nodes with the solid model. These membranes are very thin, such that they do not interfere with the mechanics of the rest of the model. Rather, they are there as strain measuring devices, as they inherit the deformations of the parent surface. The material behavior of these membranes was modeled via an abaqus VUMAT user subroutine, which represents an isotropic hyperelastic material (we have chosen the simplest Neo-Hookean material). The most important aspect of the process is that the deformation gradient is available to the VUMAT user subroutine. Local orientations were defined in cylindrical coordinate systems (spherical for the apex segments) such that deformations would be represented roughly in the circumferential and longitudinal directions.
Consider the start of strain measurement configuration at time to which corresponds a deformation gradient . From that time on, the deformation gradient can be decomposed as , where is the deformation gradient of the relative deformation between the current configuration (ES) and the start configuration (ED). We then compute the relative left Cauchy–Green strain tensor from which we extract the principal values and directions. The principal values being the square of principal relative stretches allow the determination of the principal nominal strains followed by a projection back into the deformed configuration of the global coordinate system. We then have the relative nominal strains in the tangential and longitudinal directions. These strains are carried in the analysis database as state variables (abaqus User Subroutines Reference manual) and can be postprocessed. These computational radial and circumferential strains are defined on the endocardial surface. There were three components of strain but we used the two more accurate components .
We wrote scripts to monitor the minimum ventricle volume in order to select the end time-point that produces the frame from which the relative strains are obtained. These strains are integration point values. Therefore, the script averages the relative strain values over each one of the 16 AHA segments (Fig. 4) and feeds such averages to the optimization algorithm for comparison with measurements. The 16 AHA segments are consistent with those used by tomtec software (as described in TomTec-Arena™ 1.1 Operating Manual). It should be noted that the current recommendation from the AHA is to segment the LV into 17 segments, but 16 AHA segments in our paper refers to the TOMTEC convention.
The optimization was conducted for either LV or BV model, for either homogeneous or heterogeneous contractility. In all optimization scenarios, the downhill simplex optimization method was used to minimize the differences between the computational and experimental (target) strains [36,37]. The objective function was the sum of the difference between the experimental (target) and computational strains. Since the set of target experimental strains was different in the scenarios, different objective functions were used in each scenario (Table 1). To assess the reliability of the optimization process, in scenarios #3 and #4, the values of Tmax were altered and the optimization was repeated.
Calibration of Passive Properties for Left Ventricle Model.
Calibration of passive material properties was based on the end-diastolic pressure–volume relationship or EDPVR, and led to the passive properties (Table 2). The optimization scheme provided the EDPVR curve in agreement with the analytical curve (Fig. 2). The mean difference (the sum of the differences between the experimental and computational results at the points shown in Fig. 2 divided by the number of points) between the analytical and experimental LV volumes was approximately −0.04 ml. The computational and experimental LV EDVs were 57.1 and 57.8 ml, respectively (computational EDV 1.2% lower).
Calibration of Active Properties for the Left Ventricle Model.
The experimental circumferential and longitudinal strains for 16 segments changed from segment to segment (Figs. 5 and 6). The two strains had different patterns in basal (segments 1–6), midwall (segments 7–12), and apical (segments 13–16) segments.
The segmental values for Tmax in 16 segments were altered in the optimization process to calibrate mechanical properties based on segmental strain metrics (Table 2). The difference between the mean Tmax and segmental Tmax values changed for each segment (Fig. 7, Table 2). The largest difference occurred at midwall segment 9, where Tmax (680.4 kPa) was 352.1% larger than the mean Tmax (150.5 kPa) (the sum of Tmax values at all 16 regions divided by number of regions = 150.5 kPa). The computational and experimental circumferential strains were in correspondence; the mean difference between them in 16 segments was 0.046. The largest difference occurred at apical segment, segment 14 (0.182, 295.7% higher relative to the mean value). The computational and experimental longitudinal strains were generally in agreement. The mean difference (the sum of the difference between the computational and experimental strains at target regions divided by number of target regions) between them in 16 segments was 0.059. The maximum difference was at midmyocardium at segment 13 (0.169 and 186.4% higher relative to the mean value).
When only circumferential strains were used in the calibration process (Fig. 8), the agreement between the computational and experimental strains was better than when both circumferential and longitudinal strains were used. The mean difference between the computational and experimental strains was 0.040, and the maximum difference occurred at apical segment 14 (0.169 and 323.0% higher relative to the mean value).
Similarly, the agreement between the computational and experimental longitudinal strains improved when only longitudinal strains were used in the calibration process (Fig. 9). The mean difference at 16 segments between the two strains was 0.043, and the maximum difference occurred at midwall segment 8 (0.125 and 190.7% higher relative to the mean value).
Calibration of Active Properties for the Biventricular Model.
The calibration of material properties for the BV model led to converged objective function. The mean and maximum differences between the experimental and computational circumferential strains were 0.061 and 0.120 (segment 9, 100% higher relative to the mean value), respectively. Similarly, the mean and maximum differences between the computational and experimental longitudinal strains were 0.076 and 0.133 (segment 3, 75% higher relative to the mean value), respectively. The agreement between the experimental and computational strains improved after material properties were altered segmentally at segments 3–6 (Table 3, Figs. 10–13).
The differences between the experimental and computational strains were higher when uniform Tmax values were used. The mean differences for the LV model with uniform Tmax were 0.129 and 0.076 for the circumferential and longitudinal strains, respectively (Table 3; Figs. 5 and 6). For the BV model, the mean differences for the circumferential and longitudinal strains were 0.061 and 0.076, respectively, with uniform Tmax (Table 3, Figs. 10 and 11). After segmental Tmax calibration, these differences were 0.046 (LV, circumferential), 0.059 (LV, longitudinal), 0.061 (BV, circumferential), and 0.055 (BV, longitudinal).
After the segmental values for Tmax for models #2 and 3 were calculated, the starting values for the optimization were changed to ascertain the optimized values that are not those associated with a local minimum. The time to reach convergence was affected by the starting values for the parameters (segmental Tmax), but the calibration results did not change when the starting values were altered.
In this study, we used LV ED (EDV and EDP) and speckle tracking 3D echocardiography segmental strains from beating hearts to calibrate passive and active mechanical properties of the LV. Our results show that isight can be used as a platform to improve the state-of-the-art calibration of the LV mechanical properties . The calibration results were in line with the previous reports (Table 2, rows 1 and 7). The optimization algorithm was stable because the results did not change by altering starting Tmax values. The selection of segments where strains are measured affects the calibration accuracy. When circumferential and longitudinal strains were both used, the average difference between the predicted and experimental circumferential and longitudinal strains was higher (0.046 and 0.059, respectively) than when only circumferential (0.040) or longitudinal (0.043) strains were used (Table 3). The LV volume in our calibrated model was close to the experimental value. Although we did not use end systolic volume (ESV) as a calibration criterion, the predicted and experimental ESVs were close (29.9 versus 30.9 ml).
Rationale for Using abaqus and isight.
In a previous study, ls-dyna and ls-opt (Livermore Software Technology Corporation, Livermore, CA) were used to calibrate material properties of LV based on segmental strain metrics . Here, we used abaqus and isight for the same purpose. This study extends the work seen in our previous study , by utilizing more sophisticated methods, incorporating more modern strain tracking information, advanced cardiac modelling techniques (i.e., FE and coupled lumped circulatory models) and sophisticated models (i.e., realistic fiber orientations and biventricular modelling).
Global Versus Segmental Calibration Metrics.
The calibration of the mechanical properties of myocardium can be based on global or segmental metrics. From a mathematical perspective, using global metrics is more straightforward than using segmental metrics. Here, we used 16 circumferential and 16 longitudinal strains at 16 AHA segments to calibrate the myocardium. The convergence of the objective function took longer than when we used global metrics, or when we used 16 circumferential or 16 longitudinal strains. These findings indicate that the selection of calibration parameters plays a key role in an efficient calibration process.
The optimization platform presented here was used to improve the state-of-the-art calibration process . Our results showed that the agreement between the computational strains and the experimental strains improved after segmental variations in mechanical properties were implemented (Table 3, Figs. 12 and 13). The mean differences for the circumferential strains were reduced in the LV model (0.129 versus 0.046) and were nearly unchanged for the BV models (0.061 versus 0.061). Similarly, mean differences between experimental longitudinal strains decreased in the LV model (0.076 versus 0.059) and the BV model (0.076 versus 0.055).
Rationale for Segmental Variation of Tmax in the Left Ventricle Model.
The rationale behind segmental variations in values of Tmax could be limitations in geometry, fiber directions, and constitutive equations. Our results demonstrate the importance of muscle fiber orientation in the mechanics of LV. The single LV and BV models produced results that were in agreement with experimental data. However, the BV model used a simpler approach for calibration of mechanical properties than the single LV model, and had a more realistic geometry and fiber direction. Also, the reference configuration between the LV and the BV was different because of differences in fiber orientations, geometries, loads, and boundary conditions. This reference configuration difference also makes the BV model more realistic compared to the LV model. We suspect that if, for the single LV model, the LV geometry is reconstructed from MRI and it includes the diffusion tensor magnetic resonance imaging fiber data, the results will be closer to those of the BV model.
The values of Tmax at 16 segments deviated from mean value of Tmax (Fig. 7). Within the basal segments (1–6), the largest difference between the mean Tmax and the segmental Tmax happened at segment 2. Within the midwall segments (7–12), the largest difference was at segment 9, and within apical segments (13–16), the largest difference was at segment 14. It seems that at these segments (2, 9, and 14) larger errors happen in geometry approximation and fiber orientations in the LV model. It should be noted that the scenarios presented here for LV are complicated as we considered all 16 AHA segments. If a smaller number of LV regions are considered for calibration, the optimization will be faster, and the computational results will be closer to target strains (in calibrated regions).
Model Reconstruction From Echocardiography Versus Magnetic Resonance Imaging.
In this study, we did material calibration based on two types of imaging methods. The BV model was based on ex vivo MRI data, whereas the LV model was created from in vivo echocardiography data. The accuracy of the geometry and fiber orientation obtained from MRI was better than that from echocardiography data. However, processing MRI data needs more manual work and is more challenging than processing echocardiography images. Nevertheless, the LV model created from echocardiography images can provide rough information about the material properties of the LV. Although for an accurate material calibration, MRI is more suitable, echocardiography can be used to assess alterations in material properties (not to determine exact values of the material constants). In particular, for applications related to effectiveness of treatment strategies for HF, echocardiography may be a better choice because it is less challenging to analyze and can be used for patients who use cardiac implants such as devices (LVAD) and stents.
Our methodology could have applications in monitoring the recovery of myocardium contractility in patients with HF after treatment with cardiac devices such as LVADs. Although the calibration time was relatively long (on order of 24 h for LV and on order of 14 days for BV model), when linked to machine learning techniques, the material calibration could be achieved on the order of seconds. There are machine-learning finite element studies that successfully applied time-consuming computations into fast applied outcomes  that could be used in a clinical setup. In this project, we used echocardiography imaging data. As noted above, unlike tagged MRI, speckle tracking echocardiography can be used for patients with LVADs, stents, or other cardiac implants. Also, echocardiography data can be automated by software, whereas tagged MRI requires more manual operations. There are echocardiography-based handheld devices that can be used by physicians to monitor cardiac conditions . We plan to implement machine learning algorithms that use our calibration methodology in a hand-held device that could be used by physicians. Such a device can be used to find information about the health conditions of the heart in a more convenient, less expensive, and faster way than echo and MRI machines. Moreover, if required, the methodology presented in this paper works for MRI-based strain measurements.
Some of the differences between the model-predicted strains and the experimental strains were noticeably large (Table 3). The calibration of mechanical properties depends on several factors including Tmax, ns, and fiber orientation. The fiber orientation could be approximated from diffusion tensor magnetic resonance imaging data, and ns is mostly related to physiological data from cardiac muscle multiaxial stretching data . That is why in this work, we used Tmax in 16 segments to calibrate the model. Another parameter that has a profound role in the calibration process is variations of muscle contraction in the transmural direction, which was not considered in this study. In fact, a more realistic modeling approach would consider the electrical activity and link it to the transmural variation in the onset of muscle contraction and relaxation. Moreover, atrial contraction which was not considered in this study could noticeably affect the model-predicted strains.
The constitutive equations play a crucial role in calibration of the mechanical properties. In the LV model, we assumed that the parameters of active contraction in the longitudinal and circumferential directions are only determined by fiber directions. The segmental strains vary from segment to segment, but the trend in the strain from segment to segment is different for the longitudinal and circumferential directions (Figs. 5 and 6). This observation suggests that the muscle fiber activation is likely affected by other tissue characteristics, which cause different strain patterns in the circumferential and longitudinal directions. If the constitutive equation includes parameters that reflect these different trends in the circumferential and longitudinal strains, the calibration will be improved noticeably.
The discrepancies in experimental and computational strains in the LV model could be partly due to not taking into account the right ventricle (RV) pressure on the epicardial surface of the LV septum. A more realistic model would have a different transmural pressure for the LV septum that is less than that for the LV free wall (to take into account RV pressure).
Conclusions and Future Directions
A noninvasive solution to determine mechanical properties of myocardium may serve as a diagnostic or treatment tool for HF. The properties of myocardium change in the infarcted region. Moreover, stress, which may be the stimulus for remodeling in pathological conditions like HF, can be determined if mechanical properties are known. This study presented a methodology to determine mechanical properties of myocardium in beating hearts using isight. This methodology may lead to more accurate mechanical properties as it is based on segmental metrics (such as strains), rather than global metrics such as ESV and longitudinal strain. Our methodology could be implemented in assessing the effects of other parameters on the calibration outcome. For example, the distal boundary conditions representing circulation, the effects of RV on LV, the effects of errors introduced during reconstruction of geometry, and other modeling parameters could be studied using the methodology introduced in this paper.
We thank Pamela Derish in the Department of Surgery, University of California San Francisco, for assistance with proofreading the manuscript. We wish to especially thank Professor Yuan-Cheng Fung for his mentorship. Dr. Guccione and Dr. Kassab were his graduate students at UC San Diego in the late 1980s. The method used in this study is closely aligned with the method of approach described in the introductory chapter of Dr. Fung's classic monograph Biomechanics: Mechanical Properties of Living Tissues.
This work was supported by NIH (Grant Nos. R01-HL-077921, R01-HL-118627, and U01-HL-119578; Funder ID: 10.13039/100000002).
- and =
isotropic stiffness material parameters
- AHA =
American Heart Association
the relative left Cauchy–Green strain tensor
the right Cauchy–Green tensor
the peak intracellular calcium concentration
the maximum intracellular calcium concentration
the Lagrangian strain in the myofiber direction
- EDP =
end diastolic pressure
- EDPVR =
end diastolic pressure volume relation
- ESV =
end systolic volume
the deformation gradient
the deformation gradient of the relative deformation between the current configuration and the start configuration
the deformation gradient at start of strain measurement
- and =
vectors for the local myofiber and sheet directions, respectively
- HF =
- , and =
the invariants of
determinant of the deformation gradient
the sarcomere length with the stress-free condition
the sarcomere length with no active stress
- LV =
- LVAD =
left ventricle assistive device
used to partially transfer fiber stress to the sheet
- subscripts , , and =
material parameters associated with additional stiffness in the myofiber direction, sheet direction, and the connection between the myofiber and sheet directions
the isometric tension at the longest sarcomere length and highest calcium concentration
time to arrive at peak tension
the principal relative stretch