A Novel MRI-Based Finite Element Modeling Method for Calculation of Myocardial Ischemia Effect in Patients With Functional Mitral Regurgitation

Background: Functional Mitral Regurgitation (FMR) associated with coronary artery disease affects nearly 3 million patients in the United States. Both myocardial infarction (MI) and ischemia contribute to FMR development but uncertainty as to which patients will respond to revascularization (REVASC) of ischemia alone prevents rational decision making about FMR therapy. The aim of this study was to create patient-specific cardiac MRI (CMR) informed finite element (FE) models of the left ventricle (LV), calculate regional LV systolic contractility and then use optimized systolic material properties to simulate the effect of revascularization (virtual REVASC).

Methods: We describe a novel FE method able to predict the effect of myocardial ischemia on regional LV function. CMR was obtained in five patients with multi-vessel coronary disease and FMR before and 3 months after percutaneous REVASC and a single healthy volunteer. Patient-specific FE models were created and divided into 17 sectors where the systolic contractility parameter, Tmax of each sector was a function of regional stress perfusion (SP-CMR) and myocardial infarction (LGE-CMR) scores. Sector-specific circumferential and longitudinal end-systolic strain and LV volume from CSPAMM were used in a formal optimization to determine the sector based myocardial contractility, Tmax and ischemia effect, α. Virtual REVASC was simulated by setting α to zero.

Results: The FE optimization successfully converged with good agreement between calculated and experimental end-systolic strain and LV volumes. Specifically, the optimized T_max for the healthy myocardium for five patients and the volunteer was 495.1, 336.8, 173.5, 227.9, 401.4, and 218.9 kPa. The optimized α was found to be 1.0, 0.44, and 0.08 for Patients 1, 2, and 3, and 0 for Patients 4 and 5. The calculated average of radial strain for Patients 1, 2, and 3 at baseline and after virtual REVASC was 0.23 and 0.25, respectively.

Conclusion: We developed a novel computational method able to predict the effect of myocardial ischemia in patients with FMR. This method can be used to predict the effect of ischemia on the regional myocardium and promises to facilitate better understanding of FMR response to REVASC.

Introduction

Functional mitral regurgitation (FMR) associated with coronary artery disease (CAD) is a leading cause of valvular heart disease. FMR occurs in 1.2 to 2.1 million patients with CAD (Gorman et al., 2003) in whom it doubles the risk of heart failure and death (Aronson et al., 2006).

Both myocardial infarction (MI) and ischemia contribute to FMR development. Specifically, creation of inferior MI in large animals leads to left ventricular (LV) remodeling, papillary muscle displacement, restriction of mitral leaflets and subsequent development of FMR (Gorman et al., 1995; Liel-Cohen et al., 2000). Conversely, the evidence for ischemia comes from clinical studies where in approximately half of patients, FMR improves with coronary revascularization (REVASC), whereas in the remainder FMR persists or worsens (Aklog et al., 2001; Penicka et al., 2009; Kang et al., 2011). Uncertainty as to which patients will respond to REVASC impedes rational decision-making regarding FMR management.

Cardiac magnetic resonance imaging (CMR) allows LV remodeling, tissue characteristics including MI and non-ischemic fibrosis, ischemia and regional contractile function to be assessed with high precision. Cine-CMR has been employed as a reference standard for measuring LV size and function (Codella et al., 2008, 2012; Janik et al., 2008). Stress perfusion CMR (SP-CMR) has been shown to provide high diagnostic accuracy for ischemia as can occur with obstructive CAD (Klem et al., 2006). Late gadolinium enhancement CMR (LGE-CMR) provides near exact agreement with histopathology evidenced MI (Kim et al., 1999, 2000) and CMR with non-invasive tags CSPAMM-CMR has been widely used to measure systolic myocardial deformation and strain (Ryf et al., 2002).

The relationship between MI fibrosis, ischemia and function and recovery of function after REVASC in chronically ischemic hibernating myocardium in which myocardial metabolism and contractile function are down-regulated (Wijns et al., 1998) is complex (Canty and Suzuki, 2012). In that regard, we recently measured MI fibrosis, stress perfusion and regional strain with CMR in patients with FMR associated with CAD (Morgan et al., 2018). Briefly, CSPAMM-CMR measured circumferential and longitudinal strain were modestly reduced in sectors without MI or ischemia when compared to volunteers and then decreased in an approximately linearly fashion as the amount of MI fibrosis and ischemia increased (Morgan et al., 2018). In addition, chronically ischemic hibernating myocardium can partially recover function after surgical or percutaneous REVASC (Kim et al., 2000; Hocum Stone et al., 2017). However, fewer than 10% of dysfunctional segments with LGE hyperenhancement greater than 50% recover function (Kim et al., 2000).

Finite element (FE) modeling has the unique potential to (1) determine the effect of ischemia in FMR associated with CAD and to then (2) simulate the effect of REVASC. Prior FE modeling studies by our group (Walker et al., 2005; Wenk et al., 2011) and others (Bogen et al., 1980; Genet et al., 2015) have focused on the impact of MI but without incorporation of ischemia effect. Most previous FE modeling studies of the LV with ischemia have focused on arrhythmogenesis (Wang et al., 2013; Mendonca Costa et al., 2018) without incorporation of ischemia effect on LV function. There have been FE modeling studies of ischemia effect on regional LV function in which models were based on an idealized LV shape (Bovendeerd et al., 1996), cardiac imaging (XCAT) phantom (Veress et al., 2015) and physical measurements in an excised large animal heart (Mazhari et al., 2000). However, in those studies, regions of ischemia were assumed to have little (Mazhari et al., 2000) or no (Bovendeerd et al., 1996; Veress et al., 2015) regional contractile function. To date, FE models that incorporate CMR-based tissue characteristics to determine regional contractility have not been described.

This study describes a novel FE method able to predict the effect of myocardial ischemia on regional LV function. The method leverages CMR imaging data obtained from patients with multivessel CAD and FMR before and after percutaneous REVASC. It includes patient-specific FE models in which the regional systolic contractility parameter is a function of regional CMR-based tissue characteristics including stress perfusion and MI. The aim of the study is to calculate systolic material properties of ischemic myocardium and then use optimized systolic material properties to simulate the effect of REVASC (virtual REVASC).

Materials and Methods

Five patients were selected from a larger group of patients prospectively enrolled in a protocol entitled “CMR Myocardial Tissue Based Prediction of Ischemic MR Revascularization Response” at Weill Cornell Medical College (New York, NY, United States) examining REVASC effect in patients with FMR associated with CAD. The protocol included CMR and transthoracic echocardiography before and 3 months after percutaneous revascularization (PCI). One healthy volunteer without cardiovascular risk factors was also studied. The Cornell Institutional Review Board approved this study, and written informed consent was obtained at time of enrollment.

Computational Modeling Pipeline

The proposed method is summarized in Figure 1. The CMR images were used to create LV and right ventricular (RV) surfaces (see section “Image Analysis”), to calculate 3D regional strains (see section “Constitutive Laws”), and to configure LV mechanical parameters (see section “Sector-Based Myocardial Material Parameters”). LV mechanics models were constructed using the patient-specific LV and RV surfaces, LV and RV pressures, rule-based fiber angles, diastolic and systolic myocardial constitutive relationships and LV myocardial material parameters. A formal parameter optimization framework was then constructed based on the LV mechanics models and the in vivo sector-specific circumferential and longitudinal end-systolic strain and LV volume.

FIGURE 1

Figure 1. Flowchart of the proposed method. Note that the raw data are in yellow, the derived pressure in blue, the CMR-derived geometry and strains in green, and the FE modeling process in gray.

Image Acquisition

CMR was performed with a 3.0 Tesla MRI scanner (General Electric, Waukesha, WI, United States) (Morgan et al., 2018) both before and 3 months after PCI. CMR included 4 pulse sequences. Specifically, Cine-CMR was performed with a steady-state free precession pulse sequence. Images were acquired in contiguous LV short axis and standard (2-, 3- and 4-chamber) long axis orientations (Figure 2A). SP-CMR was assessed in accordance with established methods previously validated by members of our group (Klem et al., 2006; Heitner et al., 2019). In brief, pharmacologic stress was induced with regadenoson (0.4 mg), during which gadolinium was infused (0.1 mmol/kg) and LV short axis images (4–5, evenly distributed from base-apex) were acquired using a gradient echo pulse sequence (Figure 2B). Perfusion CMR was repeated 5-min thereafter under baseline (non-stress) conditions. Myocardial infarction (LGE-CMR) was assessed using an inversion recovery pulse sequence, which was acquired 10–30 min post-gadolinium (0.2 mmol/kg) infusion and acquired in spatial orientations matched to Cine-CMR (Figure 2C). CSPAMM in contiguous LV short and long axis slices (8 mm tag spacing, 10 mm slice thickness, no gap) was performed to measure myocardial deformation and strain (Figure 2D).

FIGURE 2

Figure 2. Examples of the 4 CMR image sequences used in this study: (A) Short axis Cine-CMR with LV epicardial and endocardial and RV endocardial contours; (B) A stress perfusion image with arrows indicating the ischemic region; (C) A late gadolinium enhancement (LGE-CMR) image with red arrows indicating the infarcted region, and (D) A short axis tagged (CSPAMM) image.

Image Analysis

LV and RV Segmentation

Left ventricular short and long axis and RV short axis Cine-CMR images were contoured at LV early diastolic filling (EDF; MeVisLab, version 2.7.1, Bremen, DE) where the EDF phase was defined as the time step immediately prior to opening of the mitral valve. The LV and RV contours included the papillary muscle bases but excluded the trabeculae (Figure 2A).

To eliminate short axis image alignment breathing artifacts, an LV surface was created based on the three epicardial long axis contours. LV short axis contours were moved in plane so that the contour centroid was over the centroid of the corresponding plane from the long axis surface. After alignment, LV epicardial and endocardial and RV endocardial surfaces were created using the aligned short axis contours.

Scoring of Cine-, SP- and LGE-CMR Images

Left ventricular function, MI scar and stress perfusion were localized using the AHA/ACC 17 segment model (Cerqueira et al., 2002). LV wall thickening was measured on Cine-CMR using a 5-point scale (0 = normal (at least 50% wall thickening), 1 = mild hypokinesis, 2 = moderate hypokinesis, 3 = akinesis (no wall thickening), 4 = dyskinesis).

Myocardial infarction was identified on LGE-CMR, for which transmural extent was graded using a 5-point segmental scale (0 = no hyperenhancement; 1 = 1–25%; 2 = 26–50%; 3 = 51–75%; 4 = 76–100%) (Sievers et al., 2007). LGE in non-coronary arterial patterns (mid-wall or epicardial) was not included in MI analyses.

Ischemia was assessed semi-quantitatively on SP-CMR based on peak myocardial signal intensity during first-pass gadolinium infusion. Perfusion deficits were defined as persistent hypo-enhancement on ≥4 consecutive heart beats (Klem et al., 2006) and the severity of stress perfusion defects was graded semi-quantitatively using an established 4-point scale (0 = absent, 1 = mild, 2 = moderate, 3 = severe) based on the magnitude of hypo-enhancement (Kim et al., 2019).

CMR-Measured Strains

Regional 3D circumferential (E_cc) and longitudinal (E_ll) strain for each sector were calculated from CSPAMM images. The HARmonic Phase (HARP) method developed by Osman et al. (2000) was implemented. Implementation in MeVisLab has been previously reported (Morgan et al., 2018). Please note that all the strain presented in this study is LV end-systolic strain.

LV FE Model Generation

Using LV epicardial and endocardial surfaces as input, TruGrid (XYZ Scientific Applications, Inc., Pleasant Hill, CA, United States) was used to create patient-specific meshes of the LV composed of hexahedral elements as shown in Figure 3. Meshes with different density were created and optimal mesh density was determined as a function of parameter calculation accuracy and calculation time.

FIGURE 3

Figure 3. (A) RV endocardial surface (triangle mesh), LV epicardial surface (hexahedral mesh), and the septum (in green) determined using a ray-casting method; (B) the final 17 sector FE model with the septal sectors in green.

Each element was assigned a material direction with a rule-based fiber angle method (Bayer et al., 2012) where the myofiber helix angle varying transmurally from 60° at the endocardium to −60° at the epicardium.

Custom software (C#, Visual Studio 2017, Microsoft, Redmond, WA, United States) that included a ray casting method (Möller and Trumbore, 1997) to locate the interventricular septum (Figure 3A) was used to break the mesh into 17 AHA sectors (Cerqueira et al., 2002; Figure 3B)¹.

Loading and Boundary Condition

LV pressure (LVP) at end systole was obtained from a blood pressure cuff and LVP at end-diastole (ED) and RV end-diastolic and systolic pressures were estimated from concomitant transthoracic echocardiography data (Ommen et al., 2000; Table 1). Patient-specific LV pressure was applied to the LV endocardial surface and RV pressure was applied to the LV septal epicardial surface.

TABLE 1

Table 1. Patient-specific LV and RV pressures and LV ejection fraction (EF).

Note that these models contract from the apex to the base. Specifically, nodes at the LV base (top surface) were fixed in the valve plane (Z direction movement = 0) but able to slide in the X and Y directions during passive filling and systolic contraction.

Constitutive Laws

Passive and active myocardial constitutive laws described by Guccione et al. (1991, 1993) were used in this study. Specifically, the passive myocardium is modeled by a strain energy function, W, that is anisotropic relative to the local fiber direction:

$$W=0.5C\left(e^Q-1\right)(1)$$

where

$$Q=[b_f{E}_{11}^2+b_t(E_{22}^2+E_{33}^2+E_{23}^2+E_{32}^2)+(2)$$

$$b_{fs}(E_{12}^2+E_{21}^2+E_{13}^2+E_{31}^2)]$$

E₁₁ is fiber strain, E₂₂ is cross-fiber in-plane strain, E₃₃ is radial strain, E₂₃ is shear strain in the transverse plane, and E₁₂ and E₁₃ are shear strain in the fiber-cross fiber and fiber-radial planes, respectively, and where b_f = 49.25, b_t = 19.25 and b_fs = 17.44 (Walker et al., 2005).

A time-varying elastance model is used to simulate the active contraction of cardiac muscles (Guccione et al., 2001). The time-varying elastance model has the following form:

$$T_0=T_{max}\frac{C{a}_0^2}{C{a}_0^2+EC{a}_{50}^2}{C}_t(3)$$

where T_max is the maximum isometric tension achieved at the longest sarcomere length and maximum peak intracellular calcium concentration(Ca₀)_max. Ca₀ is the intracellular calcium concentration and ECa₅₀ is length-dependent calcium sensitivity. C_t is a time-varying variable defined as follows:

$$C_t=\frac{1}{2}\left[1-cos\left(\mathrm{\varpi }\right)\right](4)$$

where ϖ is a time and sarcomere length dependent variable that increases from 0 at the start of contraction to π during peak contraction and then decreases to 0 during diastole.

Material laws were implemented with a user-defined material subroutine in the explicit FE solver, LS-DYNA (Livermore Software Technology Corporation, Livermore, CA, United States).

Sector-Based Myocardial Material Parameters

The LGE score was an integer from 0 to 4 and the SP score was an integer ranging from 0 to 3. In each case, a score of 0 represented the absence of disease. Further LGE and SP scoring details are described above (see section “Scoring of Cine-, SP- and LGE-CMR Images”). Passive stiffness, C, was determined using the following:

$$C_n=C_H+9\cdot C_H\cdot \frac{LG{E}_n}{4}(5)$$

where C_n represents the passive stiffness for each sector and C_H the stiffness of a “healthy” sector with LGE = 0. Note that when LGE = 4 (transmural MI), C_n=10C_Hwhichis consistent with our prior work (Wenk et al., 2011).

On the other hand, both LGE and SP were assumed to effect regional contractility, Tmax according to the following linear relationship:

$$Tma{x}_n=Tma{x}_H\cdot \left(1-\frac{LG{E}_n}{4}\right)\cdot \left(1-\mathrm{\alpha }\cdot \frac{S{P}_n}{3}\right)(6)$$

where Tmax_n represents the contractility for each sector, Tmax_H the contractility for healthy sector with LGE and SP = 0, α the ischemia effect and LGE_n and SP_n the LGE and SP scores for each sector. It’s known that Tmax_n≥0 and LGE_n is ∋[0, 4] and SP_n is ∋[0, 3] then it can be determined α∋[0, 1]. Since we assumed that the LGE and SP scores were 0 for the healthy volunteer, C_n and Tmax_n for the volunteer were equal to C_H and Tmax_H, respectively.

Model Optimization

A formal optimization of C_H,Tmax_H and α were performed where the objective function for the optimization was taken to be the mean-squared errors (MSE) (Guccione et al., 2001). C_H was determined such that the FE model predicted LV ED volume matched the patient-specific in vivo measured volume. Tmax_H was estimated by minimizing the MSE between FE model-predicted and in vivo MRI-measured end-systolic longitudinal and circumferential strains and the LV end-systolic volume. The goal of the optimization is to minimize the MSE as follows:

$$Objectivefunction$$

$$=\frac{1}{2N}\left(\sum _{n=1}^N\left({\left(E_{cc,n}-{\overline{E}}_{cc,n}\right)}^2+{\left(E_{ll,n}-{\overline{E}}_{ll,n}\right)}^2\right)\right)$$

$$\mathrm{\hspace{1em}}+W{\left(\frac{V_{ES}-{\overline{V}}_{ES}}{{\overline{V}}_{ES}}\right)}^2,(7)$$

where n is the in vivo average strain at each sector (note the apex sector was excluded so n = 16), E_cc,n the calculated FE circumferential strain, E_ll,n the calculated longitudinal strain, V_ES the LV end-systolic volume and W is the weight applied to the volume term. The overbar represents the experimental in vivo measurements.

W = 10 was applied to the volume term to make the strain and volume effects more balanced. For each case, the Tmax_H was initially = 350 kPa. α was initially = 0.5 in the FMR patient models and = 0 in the healthy volunteer model. No constraints were applied to either Tmax_Hor V_ES.

Mesh Convergence Study

A mesh convergence study was performed on models of patients 2 and 5 to find the minimum number of elements needed to obtain stable calculations of C_H, Tmax_H, and α within the fastest computation time. Four sets of meshes with different density were created with the transmural elements from 3 to 4, circumferential elements from 40 to 64, and longitudinal elements from 22 to 38.

Testing With Synthetic Data

Method accuracy was determined using idealized input data. Briefly, multiple simulations of diastolic inflation and systolic contraction for Patients 2 and 5 were firstly conducted by setting α to be 0 and 1, respectively (note that C_H and Tmax_H were set as the optimized values as presented in section “Testing With Synthetic Data” under “Results”). For each case, the simulated strain and LV end-systolic volume (ESV) were used as input in our optimization framework as “experimental data.” In addition, for Patient 2, the optimization was initiated with α set at 0, 0.5, and 1 while the Tmax_H was set at 100 kPa, 350 kPa, and 600 kPa to test whether the current optimization scheme is sensitive to the initial guess of the parameters.

Virtual REVASC

Virtual REVASC was performed on patient-specific models that had an ischemia effect (α > 0) by setting α=0. As seen in Eq. 6, by setting α=0, Tmax_n becomes a function of LGE only and SP has no effect.

Virtual REVASC was validated by comparing sector specific radial strain, E_rr calculated from baseline models and after virtual REVASC with a Cine-CMR wall motion based estimation of revascularization effect. Briefly, the following relationship was developed:

$$E_{rr}\mathrm{\_}WM=E_{rr}\mathrm{\_}BL\cdot \left(1+\mathrm{\Delta }WM\cdot 50\%\right),(8)$$

where E_rr_BL is E_rr calculated from baseline models, ΔWM the difference between Cine-CMR wall motion scores at baseline and after revascularization and E_rr_WM is an estimated E_rr after revascularization based on the change in wall motion score. Eq. 8 is consistent with a wall motion score of 0 having a wall thickening of at least 50%. For instance, if E_rr_BL = 0.2 and wall motion improved by 2 points, E_rr_WM would be 0.4.

Results

Patient-specific LV and RV pressures and LV ejection fraction (EF) are shown in Table 1. Patient-specific LGE and SP scores at baseline and wall motion scores at baseline and follow-up studies are shown in Table 2. Four of the 5 patients had perfusion defects in sectors without MI/fibrosis (have over 3 sectors with SP score – LGE score >2). One of the patients had MI but no significant ischemia. Note that the wall motion scores were normal and LGE and SP scores were assumed = 0 for the healthy volunteer.

TABLE 2

Table 2. Patient-specific LGE and SP scores at baseline (BL), and wall motion (WM) scores at BL and follow up (FU) studies.

Mesh Convergence Study

The mesh convergence study determined that 7952 elements (Table 3) are required and further mesh refinement only results in an average of 2% change in C, < 2% change in Tmax_H and no change in α. Note that only the optimization results for using 7952 elements will be presented in the following sections.

TABLE 3

Table 3. Mesh convergence study.

Testing With Synthetic Data

Figure 4 shows excellent agreement between the synthetically generated “experimental strain” and the strain after optimization. The target and optimized parameters are summarized in Table 4 with a maximal difference of 1.4%, which shows the proposed method has great capability of predicting Tmax_H and α accurately. The synthetic test also demonstrated that the current optimization scheme is insensitive to the selection of the initial guess of the parameters.

FIGURE 4

Figure 4. Comparison of synthetic E_cc with FE model predicted strain in patients 2 and 5.

TABLE 4

Table 4. Results of testing with synthetic data.

Prediction of C_H, Tmax_H and α

C_H was determined such that the FE model predicted LV EDV matched the patient-specific in vivo measured EDV (Table 5).

TABLE 5

Table 5. CMR-based in vivo volumes and the FE-predicted volumes after optimization.

Figure 5A shows excellent convergence of the objective function (OF) during a representative model optimization. Figure 5B shows a surface plot of MSE of the strain with respect to Tmax_H and α, where the valley of the surface indicates lowest MSE value relating to the best match between the experimental and FE model-predicted strain. The optimized Tmax_H and α are summarized in Table 3. The optimized Tmax_H and α displayed good convergence as shown in Figure 6. The final OF value of 0.297 (an average calculated using all six subjects) was obtained indicating generally good agreement between the FE model-predicted systolic strain and the patient-specific in vivo measured strain. The LV EDV and ESV were summarized in Table 5, where the results showed that the difference between the predicted LV ESV and volume measured with CMR was on average less than 2%.

FIGURE 5

Figure 5. (A) Optimization history of the overall objective function, and (B) surface plot of mean-squared-error of strain with respect to Tmax_H and α.

FIGURE 6

Figure 6. The convergence of the parameter (lines in blue represent the upper and lower bounds) for (A) Tmax_H for healthy myocardium, and (B) ischemia effect, α, on the LV myocardial contractility, with each parameter resulted in a precise final converged optimum.

The average of E_cc along with its standard error of the mean for all sectors with LGE score <3 from overall subjects derived from the experimental data and calculated with the FE model are −0.135 ± 0.044 and −0.127 ± 0.042. Figure 7 shows the regional comparisons of E_cc calculated from CSPAMM and FEM simulations. It can be seen that the model-predicted E_cc generally has good agreement with the E_cc derived from CSPAMM.

FIGURE 7

Figure 7. Comparison of sector-specific CSPAMM and FE model predicted E_cc across all 6 subjects. E_cc derived from CSPAMM and calculated from FE simulations are in gray and red bars, respectively. The average sector specific MI/fibrosis (LGE) and stress perfusion (SP) scores are in yellow circle and blue diamond, respectively.

Virtual REVASC

Figure 8 shows the E_rr_BL, E_rr_VR and E_rr_WM for all 3 patients that were determined with positive α. In all studies, virtual REVASC underestimated the actual E_rr effect. In patients 1 and 2 the difference was mild. However, patient 3 had a small predicted alpha (0.08) and minimal improvement in E_rr with virtual REVASC in spite of the fact that at baseline there was significant ischemia and the improvement in actual wall motion after revascularization was significant.

FIGURE 8

Figure 8. Comparison of the sector specific E_rr before and after virtual REVASC for the 3 patients with positive ischemia effect, α. Sector specific MI/fibrosis (LGE) and stress perfusion (SP) scores are plotted for reference. E_rr BL represents E_rr calculated using the FE simulations at baseline, E_rr VR is the E_rr obtained from the FE simulations after virtual REVASC, and E_rr WM is an “actual” E_rr effect estimated using the qualitative wall thickening scores.

Discussion

In this study, we demonstrated a novel method to estimate the effect of ischemia on regional LV myocardial contractility. The proposed method takes the advantage of multiple CMR sequences and the FE modeling to formally optimize the LV mechanics parameters. To the best of our knowledge, this study is the first study that incorporated 17-segment-based patient-specific perfusion and LGE scores into constitutive model, which allowed regional mechanical properties to be estimated with better physiological constraints. A linear relationship was proposed to describe the effect of LV infarct and ischemia on LV active contraction. Our approach was tested with CMR data from five patients with multi-vessel CAD and moderate FMR and a healthy volunteer. Good agreement between the FE model-predicted systolic strains and the patient-specific in vivo measured strains were observed, which suggests that the optimized model was faithful to the experimental data.

Ischemia Effect

The effect of ischemia, α, on LV contractility was found to be >0 for Patients 1, 2, and 3, but = 0 for Patients 4 and 5. This is consistent with the patient-specific wall motion scores seen in Table 2 where the cumulative wall motion score in patients 1, 2 and 3 decreased (improved) after PCI but remained unchanged or increased (worsened) in patients 4 and 5.

It should be noted that stress induced perfusion defects are associated with either normal or depressed regional myocardial function and in the latter case the myocardium would be described as hibernating (Wijns et al., 1998). First, by definition, our method identifies only hibernating myocardium and simple demand ischemia without dysfunction is not considered. As a corollary, if a patient with multi-vessel CAD had a large amount of demand ischemia without LV dysfunction, our method would not identify an ischemia effect. Specifically, if the SP defect is showing simple demand ischemia without underlying contractility deficit, then α would likely be 0. If the SP defect is associated with hibernating myocardium where contractility is reduced, then α will likely be >0.

MI Stiffness

Our assumption that infarct stiffness is ten times that of normal myocardium is consistent with our prior studies (Walker et al., 2005; Wenk et al., 2011) and others (Genet et al., 2015). However, there is reason to believe that better measurement of infarct stiffness is necessary in patients with multi-vessel coronary disease and FMR. First, model accuracy was much better when sectors with transmural infarction were excluded. Second, the optimized Tmax_H was higher in Patients 1, 2, and 5 (523.3, 336.8, and 401.4 kPa respectively) than the healthy volunteer (259.4 kPa). The reason for this discrepancy is unclear but might be explained if Eq. 6 is overestimating the effect of infarct and ischemia on LV contractility. It is anticipated that future studies that measure infarct stiffness in patients with MI but without ischemia will help to better determine the relationship between infarct and myocardial contractility.

Diastolic Stiffness in Patients With Multi-Vessel CAD

As seen in Table 3, there is a large variance in optimized C_H between the healthy volunteer and patients with multi-vessel CAD and FMR. We suspect this is likely a function of myocardial fibrosis (LGE scores >0) and, consistent with this, prior clinical studies have shown the extent of LGE to correspond to impaired LV relaxation (Moreo et al., 2009). However, the cumulative LGE scores of patients 1 and 2 are similar and further studies are therefore necessary to establish our hypothesis. On the other hand, the lack of diastolic strain in this study limits our ability to better optimize passive LV material parameters.

It is interesting to note that in Figure 7 in the sectors with large LGE scores (i.e., sectors 5 and 11), the model-predicted E_cc is much smaller than that derived from CSPAMM, while this is not seen in the sectors with large SP scores (i.e., sectors 4 and 7). This may be due to the effect of MI fibrosis/LGE on either regional diastolic stiffness or systolic contractility or interaction between the two. In the current study we assumed that systolic contractility was affected by both LGE and SP scores but that diastolic stiffness was only affected by the LGE scores. Those assumptions about the form and amount of LGE effect on diastolic stiffness and systolic contractility may be overestimated. Furthermore, it should be noted that fibrosis (LGE) -based effect on diastolic stiffness will affect systolic E_cc through its effect on strain at ED and sarcomere stretch. The mismatch between the model-predicted E_cc and E_cc derived from CSPAMM in the regions with large LGE scores shows that a better understanding of the effect of LGE scores on the diastolic stiffness is important.

Model Accuracy

Results of a mesh sensitivity study suggest that the optimized LV material parameters are relatively insensitive to the mesh density, although slight differences were observed between FE models with different mesh density. Table 2 indicates that a mesh with 7952 elements may be the optimal choice for the current study because of the balanced calculation time and parameter optimization accuracy. Furthermore, synthetic test results suggest that the proposed method is capable of precisely determining the LV passive and active material properties and the ischemia effect on LV active contractility.

Virtual REVASC

Validation of virtual REVASC using Cine-CMR wall thickening shows that the magnitude of ischemia effect is underestimated. This could be due to an incorrect T_max and LGE relationship. For example, we assumed a sector with LGE = 2 had the T_max as 50% of the T_max for a healthy sector and this may be over-estimating the effect of infarct on the myocardial contractility. As a result, according to Eq. 6, the model-predicted ischemia effect might be smaller than reality to compensate the over-estimated LGE effect on myocardial contractility. A non-linear relationship between T_max and LGE and a scaling factor in front of LGE score is under development.

Figure 8 shows that patient 3 had a small predicted alpha (0.08) and minimal improvement in E_rr with virtual REVASC in spite of the fact that at baseline there was significant ischemia and the improvement in actual wall motion after revascularization was significant. The reason for this is unclear but possibilities include that measured CSPAMM strain was inaccurate. In addition, patient 3 had a dilated LV and the subsequent decrease in systolic strain magnitude might make it hard to discriminate the ischemia (SP) effect. Regardless, this points out the need for accurate strain measurement with good spatial resolution as simulation input.

Limitations

First, only systolic strains were calculated since the CSPAMM data were acquired during the systolic phase, which resulted in a fact that only systolic myocardial material parameters can be optimized. Because of the lack of diastolic strains, the main diastolic material parameter, C_H was calibrated to match the measured LV end-diastolic volume. Second, the border zone effect was not taken into consideration in the current study. Third, only a linear relationship between the infarct and ischemia and the LV contractility were employed in this study. Fourth, boundary conditions near the valve plane of the LV may affect longitudinal strain calculations in that local area. However, longitudinal strain in the body of the LV is realistic. Last, the path from end-diastole to end-systole is not specified, and our models cannot therefore provide accurate simulation of isovolumic contraction and ejection. However, simulation of diastolic filling and end-systole are accurate.

Conclusion

This study proposed a novel FE method to predict the effect of ischemia on the LV contractility in patients with multi-vessel coronary disease and FMR. The proposed method has good agreement with patient-specific in vivo measured strains and is able to predict REVASC effect.

Future studies will explore the stress perfusion and contractility relationship to determine if non-linear relationships allow greater model accuracy. Flow structure models will be implemented so that mitral regurgitation effect can be directly simulated. Last, a study with larger numbers of patient-specific models will better establish the utility of our FE method.

Data Availability Statement

All datasets generated for this study are included in the article/supplementary material.

Ethics Statement

The Cornell Institutional Review Board approved this study, and written informed consent was obtained at the time of enrollment.

Author Contributions

YZ, JW, and MR designed the study. YZ conducted simulations, analyzed results, and wrote the initial draft of the manuscript. MR created the FE models. VW, AM, and JK did data analysis. LG contributed to development of methodology. JW collected the clinical data. LG, JW, and MR supervised the study. YZ, VW, JG, JW, and MR contributed to manuscript writing/review.

Conflict of Interest

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.

Footnotes

1. ^ https://github.com/mratcliffe/ElementSearch/, available on request

References

Aklog, L., Filsoufi, F., Flores, K. Q., Chen, R. H., Cohn, L. H., Nathan, N. S., et al. (2001). Does coronary artery bypass grafting alone correct moderate ischemic mitral regurgitation? Circulation 104(12 Suppl. 1), I68–I75.

PubMed Abstract | Google Scholar

Aronson, D., Goldsher, N., Zukermann., Kapeliovich, M., Lessick, J., Mutlak, D., et al. (2006). Ischemic mitral regurgitation and risk of heart failure after myocardial infarction. Arch. Intern. Med. 166, 2362–2368. doi: 10.1001/archinte.166.21.2362

PubMed Abstract | CrossRef Full Text | Google Scholar

Bayer, J. D., Blake, R. C., Plank, G., and Trayanova, N. A. (2012). A novel rule-based algorithm for assigning myocardial fiber orientation to computational heart models. Ann. Biomed. Eng. 40, 2243–2254. doi: 10.1007/s10439-012-0593-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Bogen, D. K., Rabinowitz, S. A., Needleman, A., McMahon, T. A., and Abelmann, W. H. (1980). An analysis of the mechanical disadvantage of myocardial infarction in the canine left ventricle. Circ. Res. 47, 728–741. doi: 10.1161/01.res.47.5.728

PubMed Abstract | CrossRef Full Text | Google Scholar

Bovendeerd, P. H., Arts, T., Delhaas, T., Huyghe, J. M., van Campen, D. H., and Reneman, R. S. (1996). Regional wall mechanics in the ischemic left ventricle: numerical modeling and dog experiments. Am. J. Physiol. 270(1 Pt 2), H398–H410. doi: 10.1152/ajpheart.1996.270.1.H398

PubMed Abstract | CrossRef Full Text | Google Scholar

Canty, J. M. Jr., and Suzuki, G. (2012). Myocardial perfusion and contraction in acute ischemia and chronic ischemic heart disease. J. Mol. Cell Cardiol. 52, 822–831. doi: 10.1016/j.yjmcc.2011.08.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Cerqueira, M. D., Weissman, N. J., Dilsizian, V., Jacobs, A. K., Kaul, S., Laskey, W. K., et al. (2002). Standardized myocardial segmentation and nomenclature for tomographic imaging of the heart. A statement for healthcare professionals from the cardiac imaging committee of the council on clinical cardiology of the American Heart Association. Circulation 105, 539–542. doi: 10.1161/hc0402.102975

PubMed Abstract | CrossRef Full Text | Google Scholar

Codella, N. C., Lee, H. Y., Fieno, D. S., Chen, D. W., Hurtado-Rua, S., Kochar, M., et al. (2012). Improved left ventricular mass quantification with partial voxel interpolation: in vivo and necropsy validation of a novel cardiac MRI segmentation algorithm. Circ. Cardiovasc. Imaging 5, 137–146. doi: 10.1161/CIRCIMAGING.111.966754

PubMed Abstract | CrossRef Full Text | Google Scholar

Codella, N. C., Weinsaft, J. W., Cham, M. D., Janik, M., Prince, M. R., and Wang, Y. (2008). Left ventricle: automated segmentation by using myocardial effusion threshold reduction and intravoxel computation at MR imaging. Radiology 248, 1004–1012. doi: 10.1148/radiol.2482072016

PubMed Abstract | CrossRef Full Text | Google Scholar

Genet, M., Chuan Lee, L., Ge, L., Acevedo-Bolton, G., Jeung, N., Martin, A., et al. (2015). A novel method for quantifying smooth regional variations in myocardial contractility within an infarcted human left ventricle based on delay-enhanced magnetic resonance imaging. J. Biomech. Eng. 137:081009. doi: 10.1115/1.4030667

PubMed Abstract | CrossRef Full Text | Google Scholar

Gorman, R. C., Gorman, J. H. III, and Edmunds, L. H. Jr. (2003). “Ischemic mitral regurgitation,” in Cardiac Surgery in the Adult, eds L. H. Cohn, and L. H. Edmunds, Jr. (New York: McGraw-Hill), 1751–1770.

Google Scholar

Gorman, R. C., McCaughan, J. S., Ratcliffe, M. B., Gupta, K. B., Streicher, J. T., Ferrari, V. A., et al. (1995). Pathogenesis of acute ischemic mitral regurgitation in three dimensions. J. Thorac. Cardiovasc. Surg. 109, 684–693. doi: 10.1016/s0022-5223(95)70349-7

CrossRef Full Text | Google Scholar

Guccione, J. M., McCulloch, A. D., and Waldman, L. K. (1991). Passive material properties of intact ventricular myocardium determined from a cylindrical model. J. Biomech. Eng. 113, 42–55. doi: 10.1115/1.2894084

PubMed Abstract | CrossRef Full Text | Google Scholar

Guccione, J. M., Moonly, S. M., Moustakidis, P., Costa, K. D., Moulton, M. J., Ratcliffe, M. B., et al. (2001). Mechanism underlying mechanical dysfunction in the border zone of left ventricular aneurysm: a finite element model study. Ann. Thorac. Surg. 71, 654–662. doi: 10.1016/s0003-4975(00)02338-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Guccione, J. M., Waldman, L. K., and McCulloch, A. D. (1993). Mechanics of active contraction in cardiac muscle: part II–Cylindrical models of the systolic left ventricle. J. Biomech. Eng. 115, 82–90. doi: 10.1115/1.2895474

PubMed Abstract | CrossRef Full Text | Google Scholar

Heitner, J. F., Kim, R. J., Kim, H. W., Klem, I., Shah, D. J., Debs, D., et al. (2019). Prognostic value of vasodilator stress cardiac magnetic resonance imaging: a multicenter study with 48000 patient-years of follow-up. JAMA Cardiol. 4, 256–264. doi: 10.1001/jamacardio.2019.0035

PubMed Abstract | CrossRef Full Text | Google Scholar

Hocum Stone, L. L., Swingen, C., Holley, C., Wright, C., Chappuis, E., Ward, H. B., et al. (2017). Magnetic resonance imaging assessment of cardiac function in a swine model of hibernating myocardium 3 months following bypass surgery. J. Thorac. Cardiovasc. Surg. 153, 582–590. doi: 10.1016/j.jtcvs.2016.10.089

PubMed Abstract | CrossRef Full Text | Google Scholar

Janik, M., Cham, M. D., Ross, M. I., Wang, Y., Codella, N., Min, J. K., et al. (2008). Effects of papillary muscles and trabeculae on left ventricular quantification: increased impact of methodological variability in patients with left ventricular hypertrophy. J. Hypertens. 26, 1677–1685. doi: 10.1097/HJH.0b013e328302ca14

PubMed Abstract | CrossRef Full Text | Google Scholar

Kang, D. H., Sun, B. J., Kim, D. H., Yun, S. C., Song, J. M., Choo, S. J., et al. (2011). Percutaneous versus surgical revascularization in patients with ischemic mitral regurgitation. Circulation 124(11 Suppl.), S156–S162. doi: 10.1161/CIRCULATIONAHA.110.011254

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, J., Alakbarli, J., Yum, B., Tehrani, N. H., Pollie, M. P., Abouzeid, C., et al. (2019). Tissue-based markers of right ventricular dysfunction in ischemic mitral regurgitation assessed via stress cardiac magnetic resonance and three-dimensional echocardiography. Int. J. Cardiovasc. Imaging 35, 683–693. doi: 10.1007/s10554-018-1500-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, R. J., Fieno, D. S., Parrish, T. B., Harris, K., Chen, E. L., Simonetti, O., et al. (1999). Relationship of MRI delayed contrast enhancement to irreversible injury, infarct age, and contractile function. Circulation 100, 1992–2002. doi: 10.1161/01.cir.100.19.1992

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, R. J., Wu, E., Rafael, A., Chen, E. L., Parker, M. A., Simonetti, O., et al. (2000). The use of contrast-enhanced magnetic resonance imaging to identify reversible myocardial dysfunction. N. Engl. J. Med. 343, 1445–1453. doi: 10.1056/NEJM200011163432003

PubMed Abstract | CrossRef Full Text | Google Scholar

Klem, I., Heitner, J. F., Shah, D. J., Sketch, MH Jr, Behar, V., Weinsaft, J., et al. (2006). Improved detection of coronary artery disease by stress perfusion cardiovascular magnetic resonance with the use of delayed enhancement infarction imaging. J. Am. Coll. Cardiol. 47, 1630–1638. doi: 10.1016/j.jacc.2005.10.074

PubMed Abstract | CrossRef Full Text | Google Scholar

Liel-Cohen, N., Guerrero, J. L., Otsuji, Y., Handschumacher, M. D., Rudski, L. G., Hunziker, P. R., et al. (2000). Design of a new surgical approach for ventricular remodeling to relieve ischemic mitral regurgitation: insights from 3-dimensional echocardiography. Circulation 101, 2756–2763. doi: 10.1161/01.cir.101.23.2756

PubMed Abstract | CrossRef Full Text | Google Scholar

Mazhari, R., Omens, J. H., Covell, J. W., and McCulloch, A. D. (2000). Structural basis of regional dysfunction in acutely ischemic myocardium. Cardiovasc. Res. 47, 284–293. doi: 10.1016/s0008-6363(00)00089-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Mendonca Costa, C., Plank, G., Rinaldi, C. A., Niederer, S. A., and Bishop, M. J. (2018). Modeling the electrophysiological properties of the infarct border zone. Front. Physiol. 9:356. doi: 10.3389/fphys.2018.00356

PubMed Abstract | CrossRef Full Text | Google Scholar

Möller, T., and Trumbore, B. (1997). Fast, minimum storage ray-triangle intersection. J. Graphics Tools 2, 21–28. doi: 10.1080/10867651.1997.10487468

CrossRef Full Text | Google Scholar

Moreo, A., Ambrosio, G., De Chiara, B., Pu, M., Tran, T., Mauri, F., et al. (2009). Influence of myocardial fibrosis on left ventricular diastolic function: noninvasive assessment by cardiac magnetic resonance and echo. Circ. Cardiovasc. Imaging 2, 437–443. doi: 10.1161/CIRCIMAGING.108.838367

PubMed Abstract | CrossRef Full Text | Google Scholar

Morgan, A. E., Zhang, Y., Tartibi, M., Goldburg, S., Kim, J. J., Nguyen, T. D., et al. (2018). Ischemic mitral regurgitation: abnormal strain overestimates nonviable myocardium. Ann. Thorac. Surg. 105, 1754–1761. doi: 10.1016/j.athoracsur.2018.01.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Ommen, S. R., Nishimura, R. A., Appleton, C. P., Miller, F. A., Oh, J. K., Redfield, M. M., et al. (2000). Clinical utility of doppler echocardiography and tissue doppler imaging in the estimation of left ventricular filling pressures: a comparative simultaneous doppler-catheterization study. Circulation 102, 1788–1794. doi: 10.1161/01.cir.102.15.1788

PubMed Abstract | CrossRef Full Text | Google Scholar

Osman, N. F., McVeigh, E. R., and Prince, J. L. (2000). Imaging heart motion using harmonic phase MRI. IEEE Trans. Med. Imaging 19, 186–202. doi: 10.1109/42.845177

PubMed Abstract | CrossRef Full Text | Google Scholar

Penicka, M., Linkova, H., Lang, O., Fojt, R., Kocka, V., Vanderheyden, M., et al. (2009). Predictors of improvement of unrepaired moderate ischemic mitral regurgitation in patients undergoing elective isolated coronary artery bypass graft surgery. Circulation 120, 1474–1481. doi: 10.1161/CIRCULATIONAHA.108.842104

PubMed Abstract | CrossRef Full Text | Google Scholar

Ryf, S., Spiegel, M. A., Gerber, M., and Boesiger, P. (2002). Myocardial tagging with 3D-CSPAMM. J. Magn. Reson. Imaging 16, 320–325. doi: 10.1002/jmri.10145

PubMed Abstract | CrossRef Full Text | Google Scholar

Sievers, B., Elliott, M. D., Hurwitz, L. M., Albert, T. S., Klem, I., Rehwald, W. G., et al. (2007). Rapid detection of myocardial infarction by subsecond, free-breathing delayed contrast-enhancement cardiovascular magnetic resonance. Circulation 115, 236–244. doi: 10.1161/CIRCULATIONAHA.106.635409

PubMed Abstract | CrossRef Full Text | Google Scholar

Veress, A. I., Fung, G. S., Lee, T. S., Tsui, B. M., Kicska, G. A., Paul Segars, W., et al. (2015). The direct incorporation of perfusion defect information to define ischemia and infarction in a finite element model of the left ventricle. J. Biomech. Eng. 137:051004. doi: 10.1115/1.4028989

PubMed Abstract | CrossRef Full Text | Google Scholar

Walker, J. C., Ratcliffe, M. B., Zhang, P., Wallace, A. W., Fata, B., Hsu, E. W., et al. (2005). MRI-based finite-element analysis of left ventricular aneurysm. Am. J. Physiol. Heart Circ. Physiol. 289, H692–H700. doi: 10.1152/ajpheart.01226.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, D., Kirby, R. M., Macleod, R. S., and Johnson, C. R. (2013). Inverse electrocardiographic source localization of ischemia: an optimization framework and finite element solution. J. Comput. Phys. 250, 403–424. doi: 10.1016/j.jcp.2013.05.027

PubMed Abstract | CrossRef Full Text | Google Scholar

Wenk, J. F., Sun, K., Zhang, Z., Soleimani, M., Ge, L., Saloner, D., et al. (2011). Regional left ventricular myocardial contractility and stress in a finite element model of posterobasal myocardial infarction. J. Biomech. Eng. 133:044501. doi: 10.1115/1.4003438

PubMed Abstract | CrossRef Full Text | Google Scholar

Wijns, W., Vatner, S. F., and Camici, P. G. (1998). Hibernating myocardium. N. Engl. J. Med. 339, 173–181. doi: 10.1056/NEJM199807163390307

PubMed Abstract | CrossRef Full Text | Google Scholar