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

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. Conclusion The new 2G model may be able to provide more accurate ventricular stress and strain calculations for potential clinical applications. Combining morphological and mechanical parameters may provide better predictions for post PVR outcome.


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][5][6][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.

Modeling active contraction and expansion by using different zeroload 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 beginsystole (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][5][6][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.

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 mm 2 ; 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][5][6][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 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.

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.

The anisotropic material models for RV tissues, fiber orientation and two-layer model construction
The governing equations for all material models were: 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][5][6][7]: where I 1 and I 2 are the first and second strain invariants given by, is the original position, c i and D i 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][6][7]): where I 4 = C ij (n f ) i (n f ) j , C ij is the Cauchy-Green deformation tensor, n f is the fiber direction, K 1 and K 2 are material constants. With parameters properly chosen, it was shown that stressstrain 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]: where E ff is fiber strain, E cc is cross-fiber in-plane strain, E rr is radial strain, and E cr , E fr and E fc are the shear components in their respective coordinate planes, C, b 1 , b 2 , and b 3 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

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

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.

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 3 rd period were extracted for analysis. Because stress and strain are tensors, for simplicity, maximum principal stress (Stress-P 1 ) and strain (Strain-P 1 ) were used as the representatives and referred to as stress and strain in this paper.

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.

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.

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).

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). 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 begindiastole 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.

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.

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.

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  Table 7. Comparison of RV wall thickness, circumferential and longitudinal curvature and stress/strain between Group 1 and Group 2 at begin-filling, end-filling, begin-ejection, and end-ejection showing stress is the only parameter with significant difference between the two groups. 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 beginejection 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 + Stres-sE-F with an area under the ROC curve of 0.855.

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 [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  [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][5][6][7].

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 firstorder 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.

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.

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.