Figures
Abstract
Background
Accurate calculation of ventricular stress and strain is critical for cardiovascular investigations. Sarcomere shortening in active contraction leads to change of ventricular zero-stress configurations during the cardiac cycle. A new model using different zero-load diastole and systole geometries was introduced to provide more accurate cardiac stress/strain calculations with potential to predict post pulmonary valve replacement (PVR) surgical outcome.
Methods
Cardiac magnetic resonance (CMR) data were obtained from 16 patients with repaired tetralogy of Fallot prior to and 6 months after pulmonary valve replacement (8 male, 8 female, mean age 34.5 years). Patients were divided into Group 1 (n = 8) with better post PVR outcome and Group 2 (n = 8) with worse post PVR outcome based on their change in RV ejection fraction (EF). CMR-based patient-specific computational RV/LV models using one zero-load geometry (1G model) and two zero-load geometries (diastole and systole, 2G model) were constructed and RV wall thickness, volume, circumferential and longitudinal curvatures, mechanical stress and strain were obtained for analysis. Pairwise T-test and Linear Mixed Effect (LME) model were used to determine if the differences from the 1G and 2G models were statistically significant, with the dependence of the pair-wise observations and the patient-slice clustering effects being taken into consideration. For group comparisons, continuous variables (RV volumes, WT, C- and L- curvatures, and stress and strain values) were summarized as mean ± SD and compared between the outcome groups by using an unpaired Student t-test. Logistic regression analysis was used to identify potential morphological and mechanical predictors for post PVR surgical outcome.
Results
Based on results from the 16 patients, mean begin-ejection stress and strain from the 2G model were 28% and 40% higher than that from the 1G model, respectively. Using the 2G model results, RV EF changes correlated negatively with stress (r = -0.609, P = 0.012) and with pre-PVR RV end-diastole volume (r = -0.60, P = 0.015), but did not correlate with WT, C-curvature, L-curvature, or strain. At begin-ejection, mean RV stress of Group 2 was 57.4% higher than that of Group 1 (130.1±60.7 vs. 82.7±38.8 kPa, P = 0.0042). Stress was the only parameter that showed significant differences between the two groups. The combination of circumferential curvature, RV volume and the difference between begin-ejection stress and end-ejection stress was the best predictor for post PVR outcome with an area under the ROC curve of 0.855. The begin-ejection stress was the best single predictor among the 8 individual parameters with an area under the ROC curve of 0.782.
Citation: Tang D, del Nido PJ, Yang C, Zuo H, Huang X, Rathod RH, et al. (2016) Patient-Specific MRI-Based Right Ventricle Models Using Different Zero-Load Diastole and Systole Geometries for Better Cardiac Stress and Strain Calculations and Pulmonary Valve Replacement Surgical Outcome Predictions. PLoS ONE 11(9): e0162986. https://doi.org/10.1371/journal.pone.0162986
Editor: Miklos S. Kellermayer, Semmelweis Egyetem, HUNGARY
Received: February 10, 2016; Accepted: August 31, 2016; Published: September 14, 2016
Copyright: © 2016 Tang et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: The study involves human subjects. The Boston Children’s Hospital Committee on Clinical Investigation approved the study. The Boston Children’s Hospital IRB approval number is: IRB-CRM09-04-0237. Written informed consent was obtained from participants. Data available on request. Data cannot be made publicly available for ethical or legal reasons (public availability would compromise patient privacy). Interested researchers and individuals should contact the Data Access Committee made of DT (dtang@wpi.edu), Michael Malone (mmalone@wpi.edu), Rhonda Podell (rpodell@wpi.edu), and Tal Geva (Tal.Geva@CARDIO.CHBOSTON.ORG).
Funding: This work was supported by: National Heart, Lung and Blood Institute, 1R01-HL 089269, PIs: del Nido, Tang, Geva; National Heart, Lung and Blood Institute, R01-HL63095, del Nido; National Heart, Lung and Blood Institute, 5P50HL074734, Geva, del Nido; National Sciences Foundation of China 11171030 Yang; National Sciences Foundation of China, 1571691, Tang; and National Sciences Foundation of China, 11672001, Tang. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: Dr. Geva is a consultant to Medtronic.
1. Introduction
Accurate ventricle stress and strain calculations are of fundamental importance for cardiovascular research and investigations. From mechanical point of view, zero-stress ventricle geometry information is required for its stress/strain calculations. Ventricle modeling, especially ventricle active contraction modeling based on in vivo data, is extremely challenging because of complex ventricle geometry, dynamic heart motion and active contraction and relaxation where the reference geometry (zero-stress geometry) changes constantly in a cardiac cycle. As a first-order approximation, an approach using two zero-load geometries (2G) is proposed to model ventricle cardiac motion: one zero-load ventricle geometry is used to model the diastole phase where sarcomere has its relaxed zero-stress length, another zero-load ventricle geometry is used to model the systole phase where sarcomere has its contracted zero-stress length (therefore the zero-load systole geometry is smaller than the zero-load diastole geometry). Essentially, we are using two models to model the cardiac cycle to handle the active contraction and relaxation which are caused by zero-stress sarcomere length changes. It should be noted that zero-stress and zero-load are two different concepts. Zero-load geometries are used as an approximation since zero-stress state is really hard to get, and zero-load geometries are what we need for model construction purposes. More details are given later.
A brief review of ventricle modeling is given in the Discussion section. Active contraction is caused by sarcomere shortening which leads to increased strain and stress. A brief description is given with the use of some related terminologies: a) at the beginning of active contraction, the zero-stress sarcomere length is shortened in a very short time duration while ventricle volume has no change (isovolumic); b) since the ventricle volume does not change, while there are local SL variations and ventricle shape changes [1], the average in vivo sarcomere length under pressure (referred to as in vivo SL) may not change much when the zero-stress sarcomere length shortens (not visible in vivo) which leads to ventricle strain increase; c) the increased strain leads to the “added” stress, equivalent to the active tension in models in [2,3]. The same could be said about “active relaxation” which happens at the end of systole when the zero-stress systolic sarcomere length changes back to its zero-stress diastolic length.
The active-tension models in [2,3] is well-accepted which adds an active stress term to ventricle total stress. Since it does not change its reference geometry, strain calculation was still based on one geometry and does not reflect the effect caused by zero-stress SL changes. A new modeling approach using two different zero-load geometries (diastole and systole) was introduced in this paper to properly model active contraction and relaxation and provide ventricle diastole and systole stress and strain calculations based on their respective zero-load geometries. 2G models were constructed for 16 patients with repaired Tetralogy of Fallot (TOF) (data from a clinical trial and available for our research) and results were compared with our previously published one-geometry (1G) models [4–7]. The new morphological and stress/strain results from the new 2G models were used to identify potential predictors for post pulmonary valve replacement (PVR) outcome for ToF patients.
2. Methods
2.1. Modeling active contraction and expansion by using different zero-load diastole and systole geometries
As a first order approximation, the right ventricle cardiac cycle can be divided into 4 phases involving two different RV zero-stress geometries (diastole and systole). A short description of the 4 phases is given below since this is the base for our 2-geometry models:
Phase 1.
Filling (diastole phase). The right ventricle starts with its minimum volume under minimum pressure with minimum stress and strain. One zero-load geometry (diastole geometry) is used for this phase, corresponding to diastole zero-stress sarcomere length (SL). It should be noted that zero-stress status is a concept for stress/strain calculations. It is not observable in a living heart under in vivo conditions. At beginning-of-filling, tricuspid valve opens; RV volume increases, pressure increases, in vivo SL expands; strain and stress increases. Phase 1 ends when RV reaches its maximum volume under end-diastole pressure (denoted by P-dia) which is lower than the maximum pressure condition.
Phase 2.
Isovolumic contraction: Both tricuspid and pulmonary valves are closed; RV volume has no change; zero-stress SL shortens (changing from diastole zero-stress length to systole zero-stress length); however, this sarcomere shortening is not physically observable. Roughly, average in vivo SL does not change much (small local SL changes are possible) since RV volume does not change. So zero-stress SL shortening leads strain and stress increase (This is similar to the active tension in other models, but our model have both strain and stress increase); increased stress pushes pressure to maximum. This phase is short. This phase involves dynamic change of zero-stress sarcomere length which is very difficult to implement. It was skipped in our model.
Phase 3.
Ejection (systole phase): This phase starts from max volume, pressure, stress and strain. One zero-load geometry (for systole phase) is used for this phase, corresponding to systole zero-stress SL. At begin-systole (referred to as begin-ejection as well), pulmonary valve opens up and ejection starts; RV volume drops; in vivo SL shortens and strain decreases; pressure drops; stress drops. At end-systole (end-ejection), RV volume reaches its minimum, pressure drops to the end-systole pressure denoted as P-sys, which is greater than minimum pressure. Pressure will continue to drop in Phase 4 when systole zero-stress SL changes to diastole zero-stress SL.
Phase 4.
Isovolumic relaxation: Pulmonary valve closes (both valves closed); zero-stress SL relaxes from systole zero-stress length to diastole zero-stress length (non-contracted length); similar to the comments made in Phase 2, roughly, average in vivo SL does not change much since volume does not change; zero-stress SL relaxation leads to strain and stress decreases; pressure drops to minimum. This phase is short. It was also skipped in our model.
Our 2G model includes the diastole and systole phases described above with their respective zero-load ventricle geometries reconstructed from patient-specific MRI data. In fact, the 2G model includes two sub-models, systole and diastole models, each with its own zero-load geometry and pressure conditions. They are based on the same CMR data with continuous volume variation in a cardiac cycle. The two isovolumic phases were omitted from our model. So stress and strain values have discontinuities going between the two phases. RV volume, pressure, stress and strain achieve their minima and maxima at begin-diastole (BD) and begin-systole (BS), respectively. End-diastole and end-systole pressure, stress and strain are also available from our 2G model which were not available from our 1G models [4–7] since end-diastole and begin-systole were made identical in 1G models, as well as end-systole and begin-diastole, respectively. Therefore, 2G model represents an improvement over 1G model also in that regard.
2.2. Data acquisition and 3D geometry reconstruction
The Boston Children’s Hospital Committee on Clinical Investigation approved the study. The Boston Children’s Hospital IRB approval number is: IRB-CRM09-04-0237. Written informed consent was obtained from participants. Cardiac magnetic resonance (CMR) data before and 6 months after PVR surgery were obtained from 16 patients (8 male, 8 female, mean age 34.5) who were previously enrolled in the RV surgical remodeling trial [8]. For this analysis, we selected the 8 best (Group 1) and 8 worst (group 2) responders based on their change in RV EF from pre- to post-PVR. RV EF was chosen due to its strong association with adverse clinical outcomes in patients with repaired TOF. Table 1 describes the basic patient data and RVEF from study cohort. Briefly, ventricular assessment was performed by an electrocardiographically-gated, breath-held steady state free precession cine CMR in ventricular short-axis planes (12–14 equidistant slices covering the atrioventricular junction through the apex). The pressure data were obtained from pre-PVR cardiac catheterization procedures. The use of CMR in evaluating RV function and size in clinical practice has been well established [9–10]. The resolution of the CMR technique used in our patients was 1.6 × 1.6 × 6–8 mm reconstructed to 1.0 × 1.0 × 6–8 mm (field of view 260 × 260 mm2; matrix 160 × 160 reconstructed to 256 × 256; 30 reconstructed frames per cardiac cycle). The valve and patch positions were determined with cine MR imaging, flow data, and delayed enhancement CMR to delineate location and extent of scar/patch, and were further verified by the surgeon (PJdN, 30 years of experience) who performed PVR for those patients. Ventricular volumes and function were measured by manual tracing of endocardial and epicardial borders on each short-axis steady-state free precession cine slices throughout the entire cardiac cycle. Analyses were performed using commercially available software (QMass, Medis Medical Imaging Systems, Leiden, the Netherlands). Simpson’s method was applied to calculate end-diastolic volume (EDV), end-systolic volume (ESV), EF, stroke volume, ventricular mass, and mass-to-volume ratio. Three-dimensional RV/LV geometry and computational meshes were constructed as previously described [4–7]. Fig 1 shows sample pre-operative CMR images from a TOF patient with segmented contour plots, reconstructed 3D zero-load diastole and systole geometries (explained in Section 2.3), a front view showing the patch, scar, and the RVOT. Fiber orientation and the two-layer model construction were also shown (Fig 1d–1g) [3,11]. Since patient-specific fiber orientation data were not available, data from the current literature was used in our models.
(a) Selected CMR slices from a patient, end of systole; (b) segmented contours; (c) zero-load diastole geometry; (d) zero-load systole geometry; (e) fiber orientation from a human heart; (f) model with marked fiber orientations; (g) model with patch and scar; (h) two-layer structure.
2.3. A pre-shrink process to obtain no-load diastole and systole geometries
Under the in vivo condition, the ventricles were pressurized and the zero-load ventricular geometries were not known. In our model construction process, an iterative pre-shrink process was applied to the in vivo minimum volume ventricular geometry to obtain the two zero-load geometries so that when in vivo pressure was applied, the ventricle would regain its in vivo geometry. Shrinking is achieved by shrinking each slice (short-axis direction) by a shrinking rate and by reducing the slice distances (long-axis direction). However, if the slice was shrunk uniformly, the ventricle wall volume (the muscle) would become smaller, which should not happen. So the inner contour (inner wall of the ventricle) was shrunk more, the outer contour (ventricle outer wall) was shrunk a little less (rate was determined by volume conservation). To get the zero-load diastole geometry, we started with a 2% shrinkage, construct the model, and apply the minimum pressure to see if the pressurized RV volume matches the CMR data. If not, we adjust the shrinkage, re-made the model, pressurize it and check again. The process is repeated until RV volume matches CMR volume with error < 0.5%. For the zero-load systole geometry, assuming a 10–15% sarcomere shortening, we started with a 15% shrinkage. The same process was repeated until the pressurized RV volume under end-systole pressure (higher than the minimum pressure) matched the CMR volume data. LV geometries were handled in similar way, with proper shrinkages determined corresponding to LV pressure conditions.
2.4. Direct biaxial testing of human myocardium tissue material properties
Tissue mechanical properties are essential to computational ventricle models. However, human heart tissue material properties are not readily available from the literature. Based on the methods of Sacks and Choung [12–13] for canine hearts and informed by our previous biaxial testing [14] and the methods of Humphrey et al. and Novak et al. [15–16]. We generated the first complete multiaxial mechanical data set for ventricular tissues using a cadaveric human heart sample [7]. A detailed description of the custom planar biaxial testing device and method has been previously described [12,14]. It should be noted that the direct measurement of biaxial material data was used to verify that the modified Mooney-Rivlin model is able to fit the recorded stress-strain data. The parameter values in the Mooney-Rivlin model were adjusted for each patient to match CMR-measured volume data.
2.5. The anisotropic material models for RV tissues, fiber orientation and two-layer model construction
The governing equations for all material models were:
(1)
(2)
(3)
where σ is the stress tensor, ε is the strain tensor, v is displacement, and ρ is material density. The normal stress was assumed to be zero on the outer RV/LV surface and equal to the pressure conditions imposed on the inner RV/LV surfaces. Structure-only RV/LV models were used to optimize model computing time. These models provided RV volume, ejection fractions, and RV stress/strain values for analysis.
The RV and LV materials were assumed to be hyperelastic, anisotropic, nearly-incompressible and homogeneous. The patch and scar materials were assumed to be hyperelastic, isotropic, nearly-incompressible and homogeneous. The nonlinear Mooney-Rivlin model was used to describe the nonlinear anisotropic and isotropic material properties. The strain energy function for the isotropic modified Mooney-Rivlin model is given by Tang et al. [4–7]:
(4)
where I1 and I2 are the first and second strain invariants given by,
(5)
C = [Cij] = XTX is the right Cauchy-Green deformation tensor, X = [Xij] = [∂xi/∂aj], (xi) is the current position, (ai) is the original position, ci and Di are material parameters chosen to match experimental or patient-specific CMR measurements. The strain energy function for the anisotropic modified Mooney-Rivlin model was obtained by adding an additional anisotropic term in Eq 4 (Tang et al. [5–7]):
(6)
where I4 = Cij (nf)i (nf)j, Cij is the Cauchy-Green deformation tensor, nf is the fiber direction, K1 and K2 are material constants. With parameters properly chosen, it was shown that stress-strain curves derived from Eq 6 agreed very well with the stress-strain curves from the anisotropic (transversely isotropic) strain-energy function with respect to the local fiber direction given in McCulloch et al. [2]:
(7)
(8)
where Eff is fiber strain, Ecc is cross-fiber in-plane strain, Err is radial strain, and Ecr, Efr and Efc are the shear components in their respective coordinate planes, C, b1, b2, and b3 are parameters to be chosen to fit experimental data. Even though two different zero-load geometries were used for diastole and systole phases, parameter values in 6–8 still needed to be adjusted to fit CMR-measured RV volume data. It should be noted that Eqs 7 and 8 were used because it is desirable to use local coordinate system to identify material parameters which are independent of fiber directions. Stress-stretch curves from our old one-geometry model and the new two geometry model for the patient given in Fig 1 are given in Fig 2. Imposed RV pressure conditions are given by Fig 3. Fig 4 shows good agreement between computational and CMR-measured volume data (error < 2%).
(a) Recorded RV pressure profile; (b) RV pressure condition used in the model with Pmin (begin-filling), Pdia (end-filling), Pmax, (begin-ejection) and Psys (end-ejection) marked; (c) recorded aorta pressure profile; (d) recorded LV pressure profile.
(a) Stress-stretch curves for patch, scar and ventricle tissue used in 1G model; (b) stress-stretch curves used for the diastole phase in the 2G model; (c) stress-stretch curves used for the systole phase in the 2G model. Tff: stress in fiber direction; Tcc: stress in cross-fiber direction.
(a) 1G model agreement with CMR data; (b) 2G model agreement with CMR data.
2.6. Fiber orientation
As patient-specific fiber orientation data was not available, we chose to construct a 2-layer RV/LV model and set fiber orientation angles using fiber angles given in Hunter et al. (see Fig 1) [3]. Fig 1 shows fiber orientations from a human heart and how the 2-layer RV/LV model was constructed (Fig 1(h)) [4–7].
2.7. Mesh generation and geometry-fitting technique for patient-specific CMR-based models
Because ventricles have complex irregular geometries with patch and scar tissue components, a geometry-fitting mesh generation technique was developed to generate mesh for our models. Fig 1(h) illustrates RV/LV geometry between 2 slices. Each slice was first divided into geometry-fitting areas called “surfaces”. The neighboring slices were stacked to form volumes. Using this technique, the 3D RV/LV domain was divided into multiple “subvolumes” to curve-fit the irregular ventricle geometry with patch as an inclusion. 3D surfaces, volumes and computational mesh were made manually, slice by slice. Mesh analysis was performed by decreasing mesh size by 10% (in each dimension) until solution differences in stress/strain predictions were less than 2%. The mesh was then chosen for our simulations.
2.8. Solution methods and simulation procedures
The unsteady periodic RV/LV computational models were constructed for the 16 patients and the models were solved by ADINA (ADINA R&D, Watertown, MA, USA) using unstructured finite elements and the Newton-Raphson iteration method. Simulations were carried out for 3 periods until the solutions became period and stress/strain distributions from the 3rd period were extracted for analysis. Because stress and strain are tensors, for simplicity, maximum principal stress (Stress-P1) and strain (Strain-P1) were used as the representatives and referred to as stress and strain in this paper.
2.9. Data collection for statistical analysis
For each CMR data set, we divided each slice into 4 quarters, each with equal inner wall circumferential length. Ventricular wall thickness (WT), circumferential curvature (C-cur), longitudinal curvature (L-cur), and stress/strain were calculated at all nodal points (100 points/slice). The “slice” values of those parameters were obtained by taking averages of those quantities over the 100 points for each slice and saved for analysis. We checked the data for the proper assumptions of statistical methods. In particular, for statistical studies using the 16 patients (sample size n = 16), the assumption of normality was checked and satisfied. Pairwise T-test and Linear Mixed Effect (LME) model were used to determine if the differences from the new and old models were statistically significant, with the dependence of the pair-wise observations and the patient-slice clustering effects being taken into consideration [17]. For group comparisons, continuous variables (RV volumes, WT, C- and L- curvatures, and stress and strain values) were summarized as mean ± SD and compared between the outcome groups by using an unpaired Student t-test. Associations between pre-PVR RV parameters and the outcome (change in RV EF) were explored using Pearson correlation analysis. At the patient level (assuming independence), logistic regression analysis was used to identify pre-PVR parameters that best predicted the primary outcome—RV EF response to PVR. Prediction performance was evaluated based on 20 repeats of 2-fold cross-validation in order to stabilize the accuracy measurements for all combinations of these parameters [7,18]. Repeated cross-validation is a standard technique to reduce errors and improve prediction stability, especially when sample size is relatively small [19,20]. For each predictor (it can be a combination of the risk factors) and for each test run, the data set was divided randomly into two parts, one part (model fitting data) was used to fit the statistical predictive model, and the other (model testing data) was used to test the model and determine the sensitivity and specificity of the predictor for this test run. The test run was repeated 20 times for improved stability [7]. The whole process was performed automatically used our developed codes and the R Package. The optimal sensitivity and specificity (the largest summation of the two values) of these parameters and their area under the receiver operating characteristic curve (AUC) were determined.
3. Results
The differences of 1G and 2G models were presented first, then the differences between the two patient groups using 2G models. Post PVR outcome prediction results were presented in 3.10.
3.1. Terminologies
The purpose of this paper is to introduce the new model with 2 zero-load geometries (2G model), compare the results with our previous model which used 1 zero-load geometry (1G model), and then use the 2G models of the selected 16 ToF patients to investigate possible associations between morphological and mechanical parameters and post PVR surgical outcome. For the 1G model, results at begin-filling (BF) and begin-ejection (BE) corresponding to minimum and maximum pressure and RV volume were obtained for comparison. For the 2G model, results at begin-filling (BF), end-filling (EF), begin-ejection (BE), and end-ejection (EE) were obtained for comparison. The traditional end-systole, end-diastole, and begin-diastole conditions correspond to our begin-ejection, end-ejection, end-filling, and begin-filling, respectively. They may be used interchangeably as needed.
3.2. Begin-ejection stress from the 2G model was 28% higher than that from the 1G model
Table 2 summarizes the average stress values of the 16 patients from the 2 models. Fig 5 shows the stress plots from the 2 models using the patient given by Fig 1 as an example. According to the total average values in Table 2, begin-ejection stress from the 2G model was 28% higher than that from the 1G model (108.4 kPa vs. 84.7 kPa). The begin-filling stress values from the two models were about the same (7.17 kPa vs. 7.32 kPa, 2% difference).
BF: Begin-Filling; BE: Begin-Ejection; EF: End-Filling; EE: End-Ejection.
(a) 1G model, begin-ejection; (b) 2G model, begin-ejection; (c) 2G model, end-filling; (d) 1G model, begin-filling; (e) 2G model, end-ejection; (f) 2G model, begin-filling. Note: (a)-(c) all with maximum RV volume; (d)-(f) all with minimum RV volume.
3.3. 2G model provides end-systole and end-diastole stress conditions
It should be noted that the 2G model now provides end-systole and end-diastole stress conditions which were not available from 1G model. The right ventricles had the same volume at end-diastole and begin-systole, but RV begin-systole (peak-systole) stress average value was 151.5% higher than the end-diastole average stress (108.41 kPa vs. 43.10 kPa). Similarly, the right ventricles had their minimum volumes at both end-systole and begin-diastole, but RV end-systole stress average value was 300% higher than the begin-diastole average stress (28.89 kPa vs. 7.17 kPa). These stress conditions are of fundamental importance for many cardiovascular investigations.
3.4. Begin-systole strain from the 2G model was 39.6% higher than that from the 1G model
Table 3 summarizes the average strain values of the 16 patients from the 2 models. Fig 6 shows the strain plots from the 2 models using the same patient as Fig 5. Fig 7 gives plots of average stress/strain variations in a cardiac cycle from both models. According to the total average values in Table 3, begin-systole strain from the 2G model was 39.6% higher than that from the 1G model (0.606 vs. 0.434). The begin-diastole strain value from the 2G model was 23% lower than that from 1G model (0.048 vs. 0.062).
(a) 1G model, begin-ejection; (b) 2G model, begin-ejection; (c) 2G model, end-filling; (d) 1G model, begin-filling; (e) 2G model, end-ejection; (f) 2G model, begin-filling. Note: (a)-(c) all with maximum RV volume; (d)-(f) all with minimum RV volume.
Sudden increase at the end of diastole and sudden decrease at the end of systole reflected our omission of the two isovolumic phases. (a) Average stress from the old 1G model; (b) average stress from the new 2G model; (c) Average strain from the 1G model; (d) average strain from the 2G model.
3.5. 2G model provides end-systole and end-diastole strain conditions not available from 1G model
Table 3 shows that RV peak-systole strain average value was 42.6% higher than the end-diastole average strain (0.606 vs. 0.425), while the ventricles had their maximum volumes under both conditions. Similarly, the right ventricles had their minimum volumes at both end-systole and begin-diastole, but RV end-systole strain average value was 242% higher than the begin-diastole average strain (0.164 vs. 0.048). We are able to get those values because we used different zero-load geometries for diastole and systole phases, respectively.
3.6. Comparison of RV wall thickness and curvatures between the new and old models
Fig 5 demonstrated that both new and old models matched CMR measured RV volumes very well. Table 4 gave RV wall thickness and curvature results from the two models. The begin-diastole RV wall thickness from the new and old models were 0.575 cm and 0.564 cm, respectively (2% difference). The begin-systole RV wall thickness from the new and old models were 0.491 cm and 0.505 cm, a 3% difference. The begin-diastole RV circumferential curvatures from the new and old models were 0.638 1/cm and 0.630 1/cm, respectively (1% difference). The begin-systole RV circumferential curvatures from the new and old models were 0.524 1/cm and 0.514 1/cm, a 2% difference. For the longitudinal curvature, the begin-diastole and begin-systole values (1.216, 1.186) from the old model were very close to the begin-diastole and end-diastole values (1.193, 1.155) from the new model. The peak-systole and end-systole longitudinal curvatures from the new models were 1.263 1/cm and 1.389 1/cm, respectively. The higher longitudinal curvatures in the systole phase were due to the larger shrinkage in the longitudinal direction linked to our selected fiber orientations.
3.7. Statistical significance of the reported model differences
Pairwise T-test and Linear Mixed Effect (LME) model were used to determine if the differences from the new and old models were statistically significant with the dependence of the pair-wise observations the clustering effects taking into consideration. Patients were assumed independent. The p-values of our comparisons are given in Table 5, which indicates that the differences > 5% reported were statistically significant.
p = 0.00000 indicates that p < 0.00001. Begin-ejection = begin-systole; end-ejection = end-systole in 2G model.
3.8. Correlations between RV EF change (ΔEF) and the morphological and mechanical stress/strain parameters
Correlation analyses were performed to determine whether changes in RV EF from pre- to post-PVR were associated with RV volume, WT, C-curvature, L-Curvature, or stress/strain data and results are given in Table 6. Using the 2G model results, RV EF change correlated negatively with stress (r = -0.609, P = 0.012) and with pre-PVR RV end-diastole volume (r = -0.60, P = 0.015), but did not correlate with WT, C-curvature, L-curvature, or strain.
3.9. Group comparison: RV stress of Group 2 was much higher than that of Group 1
Table 7 summarizes the comparison of RV WT, C-curvature, L-curvature, and stress and strain values between the outcome groups at begin-filling, end-filling, begin-ejection, and end-ejection showing stress is the only parameter with significant difference between the two groups. At begin-ejection, mean RV stress of Group 2 was 57.4% higher than that of Group 1 (130.1 ±60.7 vs. 82.7±38.8 kPa; P = 0.0042). Differences at other three time points were similar. It should be noted that stress was the only parameter that showed significant differences between the two groups.
3.10. Combination of morphological and mechanical stress parameters may provide better prediction for post PVR outcome
The logistic regression method was applied to all 255 possible combinations of the 8 candidate predictors to calculate their prediction accuracy for patient’s group category. The 8 predictors are WT, C-cur, L-cur, RV volume, and stress at begin-ejection, plus three stress variations from one time point to another: StressE-D is stress difference between begin-ejection and end-filling; StressE-F is stress difference between begin-ejection and end-ejection; and StressE-C is stress difference between begin-ejection and begin-filling. Table 8 shows the 6 best combinations (out of 255) of RV parameters that correctly assigned patients to their ultimate outcome group and the prediction accuracy and ranking of the single predictions. Pre-PVR RV stress at begin-ejection was the best single predictor among the 8 individual parameters with an area under the ROC curve of 0.782. The best combination of parameters was C-cur + RV volume + StressE-F with an area under the ROC curve of 0.855.
4. Discussions
It should be noted that our 1G model starts from RV minimum volume, with minimum pressure, stress and strain conditions.
For ventricle modeling, Peskin pioneered heart modeling effort and simulated blood flow in a pumping heart with his immersed boundary method [21]. McCulloch et al. and Hunter et al. conducted comprehensive investigations for cardiac mechanics, which included passive and active ventricle modeling, the Physiome Project and the Continuity package [2–3]. Kerckhoffs et al. introduced a multi-scale approach starting from the cellular scale and building to the tissue, organ and system scales [22–23]. As heart contraction is triggered by electrical activation, Pfeiffer et al. investigated cardiac mechanics with electromechanical coupling and mechanoelectric feedback [24]. Guccione et al. presented a detailed left ventricle active contraction model with parameter values determined from canines [25]. Costa et al. studied ventricle laminar fiber architecture and 3D systolic mechanics using canine model [26]. Humphrey et al. and Novak et al. investigated ventricle tissue material properties using animal tissues [15–16]. Gan et al. used MRI-based left ventricle models and Cine-MRI to perform strain rate analysis and indicated that strain rate may be used to differentiate diabetic patients from normal controls [27]. Early magnetic resonance imaging (MRI)-based ventricle models were introduced by Saber et al. for mechanical analysis and investigations [28]. Choi et al. used MRI-based models which coupled electromechanics with hemodynamics to compare normal and diseased canine left ventricle cardiac functions [29]. Das et al., used patient-specific phase-contrast MRI (PC-MRI) to quantify flow velocity boundary conditions for ventricle models to get better flow predictions [30]. Asner et al. used tagged MRI to estimate passive and active properties in the human heart [31]. Nguyen et al. patient-specific CFD models to investigate left ventricle (LV) diastolic dysfunctions [32]. Krishnamurthy et al. described a unique approach mapping a 3D bi-ventricular model obtained from a fixed heart to patient-specific geometric models using large deformation diffeomorphic mapping [33]. These methods were tested in five heart failure patients and their results showed good agreement with measured echocardiographic and global functional parameters (ejection fraction and peak cavity pressures). Fomovsky et al. used models in design of mechanical therapies for myocardial infarction [34]. Holmes et al. indicated that image-based cardiac mechanical models could provide useful information for clinical and surgical applications [35]. In our previous papers, patient-specific MRI-based computational right ventricle/left ventricle (RV/LV) models with fluid-structure interactions were introduced to assess outcomes of various RV reconstruction techniques with different scar tissue trimming and patch sizes [4–7].
4.1. Significance of the new 2G models
Correct stress/strain calculation is of fundamental importance for many cardiovascular research where mechanical forces play a role. Ventricle remodeling, disease development, tissue regeneration, patient recovery after surgery and many other cell activities are closely associated with ventricle mechanical conditions. The 2G modeling approach is setting up the right stage for diastole and systole stress/strain calculations using proper zero-load geometries. 1G models do not use different reference geometries for systole and diastole phases, therefore have difficulties in giving right strain calculations. It should be noted that direct measurements of stress, strain, and zero-load sarcomere length are either extremely difficult or even impossible. Even by using tagging, the strain determined uses in vivo references and could not account for zero-stress SL changes. Actual ventricle contraction and relaxation are very complex. Our model is only a first-order approximation, an improvement over the 1G models. Lack of in vivo data and model construction cost are also considerations. Data from the literature or from ex vivo experiments have to be used to complete the computational models. We are in need of patient-specific data such as fiber orientation, sarcomere length contraction rate, regional material properties, etc.
4.2. Predictors for post PVR outcome
Survival of patients with tetralogy of Fallot (TOF) has steadily increased since the introduction of open-heart surgery, with operative mortality currently <2% [36]. Survival past the first two decades of life has also improved with recent reports showing a 30-year survival rate nearing 90% [37]. Since this operation was first performed in the mid-1950s, a conservative estimate projects that the number of survivors of TOF repair in the United States exceeds 100,000 and increases by 3,000–4,000 patients every year [38]. As a result of the surgical reconstruction of the right ventricular (RV) outflow tract and other operative sequelae, patients are exposed to chronic pulmonary regurgitation that leads to progressive RV dilatation and dysfunction. The current surgical approach to address chronic pulmonary regurgitation includes pulmonary valve replacement/insertion (PVR) with or without RV remodeling. However, while most patients demonstrate a variable degree of decrease in RV size, many do not experience an improvement in RV function and some show a decline after PVR [38–39]. In our previously reported randomized clinical trial of surgical remodeling of the RVOT in 64 patients with repaired TOF undergoing PVR, using available clinical data and CMR, we found no significant difference between groups with and without RV remodeling in post PVR RV EF changes [8]. It has remained unclear why some patients had experienced an improvement in RV EF whereas in others RV function had deteriorated. Our modeling added mechanical stress analysis to this study and the identified predictors provided a potential approach to identify patients and factors for possible surgical outcome improvements.
4.3. Model limitations and future directions
Several improvements can be added to our models in the future for better accuracy and applicability: a) fluid-structure interaction to obtain both flow and structural stress/strain information for complete mechanical analysis; b) direct measurements of tissue mechanical properties which will be very desirable for improved accuracy of our models; c) patient-specific fiber orientations and better estimate of zero-stress sarcomere fiber contraction rate; d) valve mechanics would enable us to investigate regurgitation and other valve-related diseases.
Small sample size (n = 16) is a limitation for the prediction method. Repeated cross-validation was used which is a standard technique to reduce errors and improve prediction stability, especially when sample size is relatively small [19, 20]. Our prediction method was also used in our previous one-geometry paper [7] where consistent results were reported. It should be noted that the only real solution for the small sample size issue is to increase the patient numbers. Since getting new patients to obtain data with follow-up data after surgery is extremely difficult and requires time and resources, and model construction takes long time (each model takes 2–4 weeks to construct and adjust; each patient needs practically 3 models: 1G, 2G systole and 2G diastole), adding more patients will be our future effort.
Author Contributions
- Conceptualization: DT PJdN CY TG.
- Data curation: DT PJdN CY HZ TG KLB.
- Formal analysis: DT PJdN CY ZW TG.
- Funding acquisition: PJdN DT TG.
- Investigation: DT PJdN CY HZ XH RHR VG AT ZW TG.
- Methodology: DT PJdN CY TG HZ XH RHR VG AT ZW KLB.
- Project administration: DT PJdN TG.
- Resources: DT PJdN KLB TG.
- Software: DT CY HZ XH AT.
- Supervision: DT PJdN TG CY.
- Visualization: CY HZ XH TG.
- Writing – original draft: DT TG PJdN ZW HZ.
- Writing – review & editing: DT PJdN TG.
References
- 1. Goetz WA, Lansac E, Lim HS, Weber PA, Duran CM. Left ventricular endocardial longitudinal and transverse changes during isovolumic contraction and relaxation: a challenge. Am J Physiol Heart Circ Physiol. 2005 Jul;289(1):H196–201. pmid:15708963
- 2. McCulloch AD, Waldman L, Rogers J, Guccione JM. Large-scale finite element analysis of the beating heart, Critical Rev. in Biomed Eng, 1992; 20(5,6): 427–449.
- 3. Hunter PJ, Pullan AJ, Smaill BH. Modeling total heart function, Annu Rev Biomed Eng., 2003; 5:147–77. pmid:14527312
- 4. Tang D, Yang C, Geva T, del Nido PJ. Image-Based Patient-Specific Ventricle Models with Fluid-Structure Interaction for Cardiac Function Assessment and Surgical Design Optimization, Progress in Pediatric Cardiology, 30:51–62, 2010. pmid:21344066
- 5. Tang D, Yang C, Geva T, Gaudette G, del Nido PJ., Multi-Physics MRI-Based Two-Layer Fluid-Structure Interaction Anisotropic Models of Human Right and Left Ventricles with Different Patch Materials: Cardiac Function Assessment and Mechanical Stress Analysis, Computers & Structures, 89, 1059–1068, 2011.
- 6. Tang D, Yang C, Geva T, del Nido PJ. Patient-specific MRI-based 3D FSI RV/LV/Patch models for pulmonary valve replacement surgery and patch optimization, J. of Biomech. Engineering, 2008;130(4) 041010.
- 7. Tang D, Yang C, Del Nido PJ, Zuo H, Rathod RH, Huang X, Gooty V, Tang A, Billiar KL, Wu , Geva T. Mechanical stress is associated with right ventricular response to pulmonary valve replacement in patients with repaired tetralogy of Fallot. J Thorac Cardiovasc Surg. 2015 Oct 3. pii: S0022-5223(15)01815-2.
- 8. Geva T, Gauvreau K, Powell AJ, Cecchin F, Rhodes J, Geva J, del Nido P. Randomized trial of pulmonary valve replacement with and without right ventricular remodeling surgery, Circulation. 2010 Sep 14;122(11 Suppl):S201–8. pmid:20837914
- 9. Geva T, Is MRI the preferred method for evaluating right ventricular size and function in patients with congenital heart disease?: MRI is the preferred method for evaluating right ventricular size and function in patients with congenital heart disease, Circ Cardiovasc Imaging. 2014;7:190–197. pmid:24449548
- 10. Geva T, Repaired tetralogy of Fallot: the roles of cardiovascular magnetic resonance in evaluating pathophysiology and for pulmonary valve replacement decision support, Journal of Cardiovascular Magnetic Resonance 2011, 13:9.
- 11. Sanchez-Quintana D, Anderson R, Ho SY, Ventricular myoarchitecture in tetralogy of Fallot, Heart, 1996;76:280–286. pmid:8868990
- 12. Sacks MS, Chuong CJ, Biaxial mechanical properties of passive right ventricular free wall myocardium, J Biomech Eng, 1993; 115:202–205. pmid:8326727
- 13. Valdez-Jasso D, Simon MA, Champion HC, Sacks MS. A murine experimental model for the mechanical behaviour of viable right-ventricular myocardium. J Physiol. 2012 Sep 15;590(Pt 18):4571–84. Epub 2012 Jul 30.
- 14. Billiar KL, Sacks MS. Biaxial mechanical properties of the natural and glutaraldehyde treated aortic valve cusp—Part I: Experimental results. J Biomech Eng. 2000;122(1):23–30. pmid:10790826
- 15. Humphrey JD, Strumpf RK, Yin FC. Biaxial mechanical behavior of excised ventricular epicardium. Am J Physiol. 1990 Jul;259(1 Pt 2):H101–8. pmid:2375396
- 16. Novak VP, Yin FC, Humphrey JD. Regional mechanical properties of passive myocardium. J Biomech. 1994 Apr;27(4):403–12. pmid:8188721
- 17. Tang D, Yang C, Geva T, del Nido PJ. Right ventricular local longitudinal curvature as a marker and predictor for pulmonary valve replacement surgery outcome: an initial study based on preoperative and postoperative cardiac magnetic resonance data from patients with repaired tetralogy of Fallot. J Thorac Cardiovasc Surg. 2014 Jan;147(1):537–8. pmid:24100105
- 18.
Hastie T, Tibshirani R, Friedman J, The Elements of Statistical Learning: Data Mining, Inference, and Prediction, Second Edition, New York, Springer, 2009.
- 19. Kim JH, Estimating classification error rate: Repeated cross-validation, repeated hold-out and bootstrap. Computational Statistics & Data Analysis, 2009; 53(11):3735–3745.
- 20. Braga-Neto UM, Dougherty ER, Is cross-validation valid for small-sample microarray classification? Bioinformatics, 2004; 20(3):374–380. pmid:14960464
- 21. Peskin CS. Numerical analysis of blood flow in the heart. J Com Phys. 25:220–252, 1977.
- 22. Kerckhoffs RC, Neal ML, Gu Q, Bassingthwaighte JB, Omens JH, McCulloch AD. Coupling of a 3D finite element model of cardiac ventricular mechanics to lumped systems models of the systemic and pulmonic circulation. Ann Biomed Eng. 2007 Jan;35(1):1–18. pmid:17111210
- 23. Kerckhoffs RC, Campbell SG, Flaim SN, Howard EJ, Sierra-Aguado J, Mulligan LJ, McCulloch AD. Multi-scale modeling of excitation-contraction coupling in the normal and failing heart. Conf Proc IEEE Eng Med Biol Soc. 2009;2009:4281–2. pmid:19963818
- 24. Pfeiffer E, Tangney J, Omens J, McCulloch AD, Biomechanics of Cardiac Electromechanical Coupling and Mechanoelectric Feedback. J Biomech Eng. 2014, 136(2):021007. pmid:24337452
- 25. Guccione JM, Waldman LK, McCulloch AD. Mechanics of active contraction in cardiac muscle: Part II—Cylindrical models of the systolic left ventricle, J Biomech Eng. 1993;115(1):82–90. pmid:8445902
- 26. Costa KD, Takayama Y, McCulloch AD, and Covell JW. Laminar fiber architecture and three-dimensional systolic mechanics in canine ventricular myocardium, Am J Physiol. 1999; 276(2 Pt 2):H595–607. pmid:9950861
- 27. Gan Y, Chen Q, Zhang S, Ju S, Li ZY. MRI-based strain and strain rate analysis of left ventricle: a modified hierarchical transformation Model, BioMedical Engineering OnLine 2015, 14, S1:S9. Epub 2015 Jan 9, 2015.
- 28. Saber NR, Gosman AD, Wood NB, Kilner PJ, Charrier CL, Firman DN. Computational flow modeling of the left ventricle based on in vivo MRI data: initial experience, Annals of Biomed. Engng., 2001; 29:275–283.
- 29. Choi YJ, Constantino J, Vedula V, Trayanova N, Mittal R. A new MRI-based model of heart function with coupled hemodynamics and application to normal and diseased canine left ventricles. Front Bioeng Biotechnol. 2015 Sep 23;3:140. pmid:26442254
- 30. Das A, Wansapura JP, Gottliebson WM, Banerjee RK. Methodology for implementing patient-specific spatial boundary condition during a cardiac cycle from phase-contrast MRI for hemodynamic assessment. Med Image Anal. 2015 Jan;19(1):121–36. pmid:25461332
- 31. Asner L, Hadjicharalambous M, Chabiniok R, Peresutti D, Sammut E, Wong J, Carr-White G, Chowienczyk P, Lee J, King A, Smith N, Razavi R, Nordsletten D. Estimation of passive and active properties in the human heart using 3D tagged MRI. Biomech Model Mechanobiol. 2015 Nov 26. [Epub ahead of print]
- 32. Nguyen VT, Wibowo SN, Leow YA, Nguyen HH, Liang Z, Leo HL. A patient-specific computational fluid dynamic model for hemodynamic analysis of left ventricle diastolic dysfunctions. Cardiovasc Eng Technol. 2015;6(4):412–29. pmid:26577476
- 33. Krishnamurthy A, Villongco CT, Chuang J, Frank LR, Nigam V, Belezzuoli E, Stark P, Krummen DE, Narayan S, Omens JH, McCulloch AD, Kerckhoffs RC. Patient-specific models of cardiac biomechanics. J. Comput. Phys. 2013; 244: 4–21. pmid:23729839
- 34. Fomovsky GM, Macadangdang JR, Ailawadi G, Holmes JW. Model-based design of mechanical therapies for myocardial infarction. J Cardiovasc Transl Res. 2011 Feb;4(1):82–91. Epub 2010 Nov 19. pmid:21088945
- 35. Holmes JW, Costa KD. Imaging cardiac mechanics: what information do we need to extract from cardiac images? Conf Proc IEEE Eng Med Biol Soc. 2006; 1: 1545–1547. pmid:17946901
- 36. Ooi A, Moorjani N, Baliulis G, Keeton BR, Salmon AP, Monro JL, Haw MP. Medium term outcome for infant repair in tetralogy of Fallot: indicators for timing of surgery. European Journal of Cardio-thoracic Surgery 2006;30:917–922. pmid:17052914
- 37. Chiu SN, Wang JK, Chen HC, Lin MT, Wu ET, Chen CA, Huang SC, Chang CI, Chen YS, Chiu IS, Chen CL, Wu MH. Long-Term Survival and Unnatural Deaths of Patients With Repaired Tetralogy of Fallot in an Asian Cohort. Circulation Cardiovasc Qual Outcomes. 2012;5:120–125.
- 38. McKenzie ED, Khan MS, Dietzman TW, Guzmán-Pruneda FA, Samayoa AX, Liou A, Heinle JS, Fraser CD Jr. Surgical pulmonary valve replacement: a benchmark for outcomes comparisons. J Thorac Cardiovasc Surg. 2014;148(4):1450–3. pmid:24703628
- 39. Waien SA, Liu PP, Ross BL, Williams WG, Webb GD, McLaughlin PR. Serial follow-up of adults with repaired tetralogy of Fallot. J Am Coll Cardiol. 1992;20: 295–300. pmid:1634663