Ventricle stress/strain comparisons between Tetralogy of Fallot patients and healthy using models with different zero-load diastole and systole morphologies

Patient-specific in vivo ventricle mechanical wall stress and strain conditions are important for cardiovascular investigations and should be calculated from correct zero-load ventricle morphologies. Cardiac magnetic resonance (CMR) data were obtained from 6 healthy volunteers and 12 Tetralogy of Fallot (TOF) patients with consent obtained. 3D patient-specific CMR-based ventricle models with different zero-load diastole and systole geometries due to myocardium contraction and relaxation were constructed to qualify right ventricle (RV) diastole and systole stress and strain values at begin-filling, end-filling, begin-ejection, and end-ejection, respectively. Our new models (called 2G models) can provide end-diastole and end-systole stress/strain values which models with one zero-load geometries (called 1G models) could not provide. 2G mean end-ejection stress value from the 18 participants was 321.4% higher than that from 1G models (p = 0.0002). 2G mean strain values was 230% higher than that of 1G models (p = 0.0002). TOF group (TG) end-ejection mean stress value was 105.4% higher than that of healthy group (HG) (17.54±7.42kPa vs. 8.54±0.92kPa, p = 0.0245). Worse outcome group (WG, n = 6) post pulmonary valve replacement (PVR) begin-ejection mean stress was 57.4% higher than that of better outcome group (BG, 86.94±26.29 vs. 52.93±22.86 kPa; p = 0.041). Among 7 selected parameters, End-filling stress was the best predictor to differentiate BG patients from WG patients with prediction accuracy = 0.8208 and area under receiver operating characteristic curve (AUC) value at 0.8135 (EE stress). Large scale studies are needed to further validate our findings.


Introduction
Patients with repaired tetralogy of Fallot (TOF) account for the majority of cases with late onset right ventricular (RV) failure. Surgical repair of TOF usually involves incision into the right ventricular outflow tract (RVOT), pulmonary valve replacement/insertion (PVR), removal of obstructive cardiac muscle, removal of obstructive pulmonary valve tissue, and placement of a patch made of non-contracting tissue or synthetic material to widen the RVOT and pulmonary valve annulus. The current surgical approach, which includes PVR and often RV remodeling, has yielded mixed results with many patients not recovering RV function after pulmonary valve insertion with or without concomitant RV remodeling surgery [1,2,3]. Therefore, there is a need to identify potential predictors which could identify patients with better post-PVR outcome from patient with worse post-PVR outcome [4,5].
It is well accepted that mechanical stress and strain conditions play an important role in cardiovascular disease initiation, development, treatment, organ remodeling, and surgical recovery process. Accurate ventricle stress and strain calculations are of fundamental importance for ventricle disease assessment, surgical design, and post-surgical recovery. McCulloch et al., Hunter et al., Holmes et al., Costa et al., and many other authors have made great contributions to passive and active ventricle modeling, including the Physiome Project and the Continuity package [6][7][8][9][10][11]. Early cardiac magnetic resonance imaging (CMR)-based ventricle models were introduced by Axel and Saber et al. for mechanical analysis and investigations [11,12]. We have introduced patient-specific CMR image-based computational right and left ventricle (RV/LV) models with fluid-structure interactions (FSI) to assess outcomes of various RV reconstruction techniques with different scar tissue trimming and patch sizes [4,5,13,14,15]. A recent review can be found in [15].
From mechanical point of view, zero-stress ventricle geometry information is required for its stress/strain calculations. Computational models often need to start from geometries with zero stress/strain conditions. Due to active contraction and myocardium sarcomere zero-stress length shortening, ventricle zero-load geometries for diastole and systole phases are different, with zero-load systole geometry being smaller than that of diastole phase. In this paper, different zero-load systole and diastole ventricle geometries were used to obtain more accurate stress/strain calculations: 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. Essentially, we are using two models (called 2G 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.
In this paper, patient-specific 2G models were constructed for 6 healthy volunteers and 12 TOG patients based on cardiac magnetic resonance (CMR) data to obtain ventricle stress/ strain conditions. Results from 2G models were compared with that from previous models with one zero-load geometry to observe the improvements by the 2G models. Differences between systole and diastole stress/strain conditions and between healthy and TOG patients were identified. Ventricle stress/strain and other morphological parameters were used to differentiate patients with better post-PVR outcome from patients with worse post-PVR outcome to identify potential possible predictors. Ventricle stress and strain of healthy and TOG patients calculated using different zero-load diastole and systole morphologies will be good contributions to the literature.

Data acquisition and 3D geometry reconstruction
The study was approved by Boston Children's Hospital Committee on Clinical Investigation. The IRB approval number is: IRB-CRM09-04-0237. Written informed consent was obtained from participants. CMR data from 6 healthy volunteers and 12 TOF patients before and 6 months after PVR were obtained. Demographic and CMR data for the 18 participants (11 male, average age 32.1) are given by Table 1. The 18 participants were divided into healthy group (HG) and TOF group (TOG) for comparison purposes. Depending on the post-PVR outcome, 12 TOF patients were further divided into Better-outcome group (BG) and Worseoutcome group (WG) corresponding to the post-PVR outcome for prediction analysis. Ventricle ejection fraction (EF) is commonly used as a measure for ventricle cardiac functions. In this paper, EF variance after PVR (ΔEF, difference between post-PVR EF and pre-PVR EF) was used as the criterion for BG and WG. Since we only have a small group, for the 12 TOF patients, the "best" way for us is to divide them into two even groups, 6 in BG and 6 in WG. This is purely for better comparisons with no clinical reasons. Ideally, we should use positive and negative ΔEF as a criterion for BG and WG. That would give us 4 in BG and 8 in WG, leading to poorer results and uneven comparisons. For our even division (BG = 6, WG = 6), the median ΔEF (-5%) was used as our criterion, i.e., TOF patients with ΔEF>-5% were classified into better group, and the rest were classified into worse group. P-values from Wilcoxon rank sum test in Table 1 indicated that there was no significant age difference between BG and WG, and between HG and TG (TOF patient group). Data acquisition techniques were previously published [14]. Image segmentation were performed with commercial software (QMass, Medis Medical Imaging Systems, Leiden, the Netherlands). Simpson's method was used to calculate end-diastolic volume (EDV) and end-systolic volume (ESV). Stroke volume and ejection fraction were determined from EDV and ESV. Location of patch and valve of TOF patients were determined with flow data, cine MR imaging, and delayed enhancement CMR, and were further verified by the surgeon (PJdN, over 30 years of experience) who performed PVR for those patients. Ventricle pressure data from TOF patients were obtained from cardiac catheterization procedures. RV pressure in healthy controls was derived from Doppler measurement of the tricuspid regurgitation jet velocity by echocardiography. Fig 1 shows 3 selected CMR slices, corresponding segmented contour plots, zero-load systole and diastole ventricle geometries, reconstructed 3D RV/LV geometry with scar, patch and myocardium fiber orientation, and measured ventricle diastole and systole pressure conditions from a TOF patient.

Stress/strain definition and calculations and the pre-shrink process to obtain zero-load diastole and systole geometries
Mechanical stress and strain are defined and should be calculated using zero-stress reference frames. Fig 2 gave an illustration of the concept using a one-dimensional bar. Fig 2 (A) gives a bar with its zero-stress (un-stretched) length L 0 and stretched length L. The 1D strain and stress (for simplification) are defined by: where YM (Young's modulus) is a commonly-used material stiffness parameter. Knowing zero-stress length of the bar is essential for its stress/strain calculations. Fig 2(B) and 2(C) gave bar illustrations of myocardium fiber (sarcomere) zero-stress diastole and systole length and its length when stressed in vivo. Sarcomere shortening is reflected by the shorter length (L 0-sys ) of the bar in Fig 2(C). This bar plot also shows how sarcomere shortening leads to increased ventricle stress and strain values (note the bar length L in (b) and (c) remains the same). Ventricle material models and 3D stress/strain calculations are given in 2.3 and 2.4. Since ventricle is under pressure, its accurate zero-stress geometry is unknown under invivo conditions. Our model construction needs to start from zero-stress geometries (the reference geometries) with zero pressure and zero stress and strain assigned. It should be noted that zero-load and zero-stress are two different concepts, since the residual stress exists in myocardium tissues [16]. However, it is near impossible to obtain ventricle zero-stress geometry from in vivo CMR data. Thus, zero-load geometries were used in our model construction process as approximations to zero-stress geometries. Due to myocardium active contraction and relaxation, zero-stress sarcomere lengths are different in systole (shortened) and diastole (relaxed) phase. Two different zero-load (diastole and systole) geometries are needed to simulate ventricle diastole and systole phases, respectively. These models are called 2G models since two zero-load geometries are used [13]. The 2G model is our way to change ventricle zeroload geometry between systole and diastole phases. A 2G model is practically using two models (each with a different zero-load geometry) to complete the cardiac cycle. In our 2G model construction process, a pre-shrink process was applied to in vivo ventricle CMR data with minimum volume to obtain the approximate two zero-load geometries so that when in vivo pressure was applied, the ventricle would regain its in vivo geometry [13]. Shrinking is achieved by shrinking each slice (short-axis direction) with an initial short-axis shrinking rate and by reducing the slice distances (long-axis direction) with an initial long-axis shrinking rate. 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 so that the ventricle wall volume was conserved. To get the zero-load diastole geometry, we started with a 2% shrinkage, construct the model, and apply the minimum pressure (begindiastole) 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. The zeroload systole geometry is smaller reflecting shortened zero-stress sarcomere length. Fig 1 provided the zero-load diastole and systole geometries from one TOF patient.
In the 1G model construction process, only one zero-load geometry was used to stimulate the cardiac cycle, with maximum pressure corresponding to maximum ventricle volume and minimum pressure corresponding to minimum volume, respectively. In 1G models, end-filling state and begin-ejection state were the same with identical pressure, volume, and ventricle stress and strain conditions; the same were true for begin-filling and end-ejection. That means 1G models used wrong pressure conditions for end-filling and end-ejection. Results for endfilling and end-ejection from 1G models could not be right. Those were 1G model limitations.

Governing equations and material models
Ventricle stress and strain were solved from the governing equations given below: where σ is the stress tensor (which has 6 components), ε is the strain tensor (6 components), v is displacement vector (3 components), p is ventricle pressure, 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. Eqs (3) and (4) contain 15 unknown functions (stress: 6; strain: 6; displacement: 3), but only 9 equations (Eq (3) has 3 scalar equations; Eq (4) has 6 scalar equations). The other 6 equations come from material properties linking stress and strain together. To that end, 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,[13][14][15][17][18][19]: 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 (6) (Tang et al. [4,5,[13][14][15][17][18][19]): 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 the strain energy function determined, the stress-strain relationship is given by: Since stress tensor is symmetric, Eq (9) only has 6 scalar equations. Now Eqs (3), (4) and (9) include 15 scalar equations to determine stress, strain and displacement, altogether 15 scalar functions.
Because tissue mechanical properties are essential for computational ventricular modeling, we generated the first complete biaxial mechanical data set for ventricular tissues using a cadaveric normal human heart sample obtained from NDRI-National Disease Research Interchange following required procedures and protocol approved by NDRI (see Fig 2) [14]. Detailed description of the custom biaxial testing device and method has been previously described [20,21]. Myocardium tissue came from a 30-year-old female who was deceased because of head trauma. Approximately 2 cm by 2 cm full-thickness sections of the right free ventricular wall and thin left ventricular wall specimens were dissected from the cadaver heart obtained within 24 hours of harvest from the donor from National Disease Research Interchange (NDRI, Philadelphia, PA). Upon arrival the heart was perfused and stored in chilled cardioplegic solution to eliminate contraction of the muscle. For the right ventricular free wall, one epicardial sample was obtained due to the thinness of the wall. The preferred fiber direction was determined by placing four fiduciary (graphite chip) markers in the center of each sample, mounting it on the biaxial test machine via 16 stainless steel hooks with the base-apex direction aligned with the machine axes, then loading it equibiaxially and recording the displacement of the markers. The stiffest material axis will be determined by the deformation gradient, its orientation relative to the base-apex direction recorded, and the sample was trimmed to align with the material axes (~15mmx15mm), and the sample re-mounted on the biaxial test device. Five tests with different fiber-direction vs. cross-fiber-direction stress ratios were run and data were recorded for model fitting. We were able to choose parameter values in our modified Mooney-Rivlin model to fit our direct measurement of biaxial stress-strain data (see Fig 3).
It should be noted that the ex vivo biaxial testing data was only used to support the choice of our material model, i.e., the modified anisotropic Mooney-Rivlin model is able to represent ventricle tissue anisotropic material properties. The ex vivo biaxial testing data was not directly used in this paper. Patient-specific ventricle parameter values in the material model for each patient were determined using CMR-measured RV volume data and an iterative procedure [4,5,[13][14][15][17][18][19]. The parameter determination procedure and parameter values for all participants are given in the appendix.
The orientations of myofibrils cause the myocardial tissue to exhibit mechanical anisotropy [22]. Since patient-specific fiber orientation data was not available, we chose to construct a 2-layer RV/LV model and set fiber orientation angles using the fiber angles published by Hunter et al. (Fig 1) and available human data [7,23]. In our 2-layer models, left ventricle (LV) fiber orientation was approximately -60 o (relative to circumferential direction) at the epicardium and +80 o at the endocardium. RV fiber orientation was set at -45 o at the epicardium and +40 o at the endocardium [7,23]. Fiber orientations could be adjusted when patient-specific data becomes available.

Solution methods
The RV/LV computational models were constructed for the 18 participants and solved by ADINA (ADINA R&D, Watertown, MA, USA). Model construction and solution procedures were reported previously [13][14][15]. Parameters in the Mooney-Rivlin models were adjusted iteratively until good agreement between the computational and CMR-measured volume data was found (error<0.2%). Mesh analysis was performed by reducing mesh density by 10-15% incrementally in each dimension until solutions became mesh independent, i.e., differences of solutions (stress/strain) between two consecutive meshes became smaller than the set tolerance (2%). Mesh density is the step size for an edge (line) of a volume ADINA uses to generate mesh for the given volume. Since each RV/LV model consisted of hundreds of manually-constructed small volumes and mesh was generated by ADINA with specified mesh density (or number of divisions for each edge of the volume), mesh density was refined by increasing number of divisions in each edge for a given small volume. That was done for all volumes of the given RV/LV model.  strain error plots showing convergence. L1 norm was to calculate stress and strain errors between two consecutive meshes. The average of mesh density for all volumes for a given mesh was used to make Fig 4D. When mesh was settled for a model, simulation procedures were continued for 3 periods when differences in the solutions between the last two periods became less than 0.1%. The solutions for the last period were accepted for analysis. Three periods were used since a) our simulation started from zero pressure with zero stress/strain so Period 2 solutions differed from Period 1 solutions considerably; b) the solutions became periodic from the second period so Period 2 and Period 3 solutions were almost identical.

Data extraction and preparation for comparative and statistical analysis
Since stress and strain are tensors, maximum principal stress and strain (both are scaler functions) were used to as the representative stress/strain variables for comparative analyses. For each ventricle model, 100 evenly-spaced points were assigned for each slice. Maximum principal stress and strain values on those points from all the slices were extracted from the model and their mean values were recorded for analysis [14]. Average right ventricle stress and strain values at begin-filling (BF), end-filling (EF), begin-ejection (BE) and end-ejection (EE) were recorded for analyses. The Wilcoxon rank sum test for unpaired observations was used to compare stress and strain differences between different groups. The Wilcoxon signed rank test for paired treatments was selected to compare stress and strain values at the 4 selected key time points (BF, EF, BE, EE). Spearman rank correlation was used to calculate correlation between computational stress/strain and CMR-measured volume data.
For a given patient, it is of clinical importance to predict if the patient would have better post-PVR outcome based on pre-surgery data. For our study, this comes to determine if a patient belongs to BG (better outcome group) or WG (worse-outcome group) based on patient's baseline data, without knowing post-PVR outcome. To predict if a patient belongs to BG or WG using pre-PVR data, we used the logistic regression model: where y i denotes the ith patient and W i denotes predictor W's value from that patient. y i = 1 (positive outcome) indicates that the model predicted that the patient belongs to BG, while y i = 0 (negative outcome) indicates that the model predicted that the patient belongs to WG. RV stress, strain, RVEDVi, EVESVi, gender, EF, and age was used as the predictors to check their prediction accuracy. For "sex" variable, male was represented with 1 and female was represented with 0 in the analysis. To stabilize the result, a 5-fold cross-validation procedure was adopted and repeated 200 times automatically with a random partition of training and testing groups each. The gradient descent method was used to obtain the coefficients by the training data, and the logistic regression model obtained from the training set was applied to the testing set to make predictions. Prediction results were compared to actual post-PVR outcome to calculate prediction accuracy, sensitivity, and specificity defined as: where TP is the number of true positive outcomes, TN is the number of true negative outcomes, FP is the number of false positive outcomes, and FN is the number of false negative outcomes. Receiver operating characteristic (ROC) curve together with false positive rate and true positive rate were determined to qualify classification efficiency.

2G model stress/strain results showed considerable improvements over 1G results
Being able to provide end-diastole and end-systole stress/strain values is one of the major improvements of 2G models over 1G models since 1G models did not have those results. Results from 2G models included stress/strain values at begin-filling (BF, corresponding to minimum volume with minimum pressure), end-filling (EF, corresponding to maximum volume with end-diastole pressure), begin-ejection (BE, corresponding to maximum volume with maximum pressure) and end-ejection (EE, corresponding to minimum volume with end-systole pressure). Results from 1G models included stress/strain values at begin-filling (BF, corresponding to minimum volume with minimum pressure) and begin-ejection (BE, corresponding to maximum volume with maximum pressure). 1G models treated EF and BE (and EE and BF as well) as identical since 1G models did not use end-diastole and end-systole pressure conditions at end-diastole and end-systole, respectively. Tables 2 and 3 summarized the stress and strain values from the 18 participants at BF, EF, BE, and EE, respectively. At end-ejection, 2G mean stress value was 321.4% higher than that of 1G models (14.54±7.41vs. 3.45±1.62, p = 0.0002). 2G mean strain values was 230% higher than that of 1G models (0.099±0.015 vs. 0.030±0.009, p = 0.0002). At begin-ejection, mean stress of 2G models was 4.9% higher than that of 1G models (61.49±26.83kPa vs. 58.60±25.41, p = 0.0010). Mean strain of 2G models was 23.1% than that of 1G models (0.548±0.126 vs. 0.445 ±0.125, p = 0.0002). The differences at end-ejection was much greater because when pressure and deformation were smaller, the stress/strain differences caused by the zero-load geometries became far more noticeable. Differences in different groups can be observed from Tables 2 & 3.

Mean RV stress value at begin-ejection is higher than that at end-filling
1G models treated end-filling and begin-ejection as identical (they were the same time point connection systole and diastole phases) so 1G models could not provide the stress/strain differences between end-filling and begin-ejection. In real cardiac cycle, active contraction happens between end-filling and begin-ejection (it is called the isovolumic phase), sarcomere zerostress length shortens, ventricle stress/strain and pressure increase while its volume is kept unchanged (both valves are closed). 2G models actually have stress/strain values at end-filling and begin-ejection from their diastole and systole phases, respectively. Mean stress of all 18 cases at begin-ejection was 150% higher than end-filling mean stress value (61.49±26.83kPa vs. 24.64±10.12kPa, p = 0.0002). Begin-ejection mean strain value was 23.26% higher than that of end-filling (0.548±0.126 vs. 0.444±0.126, p = 0.0002). The higher stress and strain at begin-systole is caused by sarcomere shortening in iso-volume contraction (active contraction). The average wall stress and strain values reach their maximum values in begin-ejection due to the high pressure and smaller zero-load configurations reflecting sarcomere shortening.

Mean RV stress value in begin-filling is lower than that of end-ejection
Similar to isovolumic contraction, isovolume relaxation happens between end-ejection and begin-filling from, 2G models are able to calculate RV stress and strain at end-ejection and begin-filling with their respective systolic and diastolic zero-load geometries and corresponding pressure conditions. Mean stress value at end-ejection was 321.45% higher than that of begin-filling (14.54±7.41kPa vs. 3.45±1.62kPa, p = 0.0002). Mean strain value in begin-filling was 233.3% higher than that of begin-filling mean strain value (0.030±0.009 vs. 0.099±0.015, p = 0.0002). The average wall stress and strain values reach their minimum values at begin-filling due to the lower pressure and sarcomere relaxation. Table 4 and Fig 5 provided stress/strain comparison between different groups at BF, EF, BE and EE with their respective p-values. In general, mean stresses from TOG group were higher than that from the healthy group. Statistically significant stress differences between HG and TG were found at begin-filling and end-ejection. At begin-filling, TG mean stress value was 94.34% higher than that of HG (4.12±1.59kPa vs. 2.12±0.39kPa, p = 0.0414). At end-ejection, TG mean stress value was 105.3% higher than that of HG (17.54±7.42kPa vs. 8.54±0.92kPa, p = 0.0245). Mean strain values of healthy individuals is 40.7% higher than that of TOF patients at end-filling (0.550±0.055 vs 0.391±0.119, p = 0.0008) and 28% higher than that of    strain at BF, EF, BE and EE all included for comparisons. RV stress had the best AUC (0.8135) and prediction accuracy (0.8208). Among the begin-filling, end-filling begin-ejection and endejection time points, end-filling stress had the best prediction accuracy at 0.8208, noticeably better than its accuracies at BE (0.7638), BF (0.8054), and EE (0.7879). Other the other hand, end-ejection stress had best AUC value at 0.8135. Prediction accuracy of age was 0.6542, just below the stress predictors. RVEDVi and RVESVi had prediction accuracies at 0.6296 and 0.5675, respectively. Prediction accuracy of ejection fraction was only 0.5000. Ejection fraction, sex and strain all had accuracies at about 0.5, indicating poor predicting power.

Correlations between RV stress/strain and MRI-derived parameters
Correlation analyses were performed to determine whether RV stress/strain were associated with MRI-derived parameters including ejection fraction change (4EF), RV end-diastole volume and end-diastole volume index. Results are given in Table 6. Using the 2G model results, RV EF change correlated negatively with end-filling stress (r = -0.650, p = 0.026) and beginejection stress(r = -0.608, p = 0.040). Other correlations were not found. It is worth noting that while end-filling and begin-ejection stress/strain mean values differ considerably, their correlation results were similar.

2G models can provide more accurate stress/strain calculations with proper zero-load ventricle geometries, especially end-diastole and endsystole stress/strain conditions
It is self-evident from the results presented that 2G models could give end-diastole and endsystole stress/strain conditions that 1G models could not produce due to model limitations.
Using correct reference frame is a basic requirement for correct stress/strain calculations [13][14][15]. 2G model could be considered an improvement over 1G model based on stress/strain theoretical definitions in biomechanics. The ultimate gold standard to prove that 2G models are Ventricle stress/strain of TOG patients & healthy more accurate than 1G models requires the knowledge of zero-stress reference ventricle diastole and systole geometries which is extremely difficult to have if after all possible. One limitation of 1G models is that not only 1G models do not change zero-load geometry between systole and diastole phases, they also use wrong pressure conditions at end-filling and end-ejection. As a result, end-filling and end-ejection stress/strain calculations using 1G models would have large errors and could be fairly misleading.  Various circumstances and applications may require accurate ventricle stress/strain calculations, such as measuring how hard the ventricle is working, energy consumed, and designing ex vivo experiments studying myocardium material properties. 2G models could provide more realistic stress/strain calculations in those cases. Table 7 gave post-PVR outcome prediction results using 1G models. It should be noted that 1G models had only two extreme time points in a cardiac cycle (begin-filling and begin-ejection) since (end-filling, begin-ejection) and (end-ejection, begin-filling) were modeled as identical points. Our 1G and 2G begin-filling conditions were the same. So Table 7 had prediction results for only one time point. 1G prediction results were in line with 2G results given in Table 5. In particular, prediction accuracy (0.8141) using 1G Begin-Ejection stress was better than that (0.7638) using 2G Begin-Ejection stress, and slightly lower than the best 2G accuracy (0.8206) using 2G End-Filling stress. It should be noted that whether stress from 2G model was a better predictor or not does not prove or disprove if 2G models gave better stress calculators. Those two were separate issues.

Our mechanical stress/strain study compared to other TOF patient studies in the literature
Our previous publications indicated that mechanical stress could be used to differentiate TOF patients with better post-PVR outcome from TOF patients with worse post-PVR outcome [4,5]. While our purpose of this modest paper was mainly to provide some benchmark stress/ strain data for TOF patients and healthy and their comparisons, many investigators have made interesting studies, often based on morphological features with aims for clinical applications. Leonardi et al. used an unbiased shape analysis framework and linear regression to analyze 38 rToF patients [24]. They found that regurgitation severity was significantly associated with RV dilatation (p = 0.01) and associated with bulging of the outflow tract (p = 0.07) and a dilatation of the apex (p = 0.08). Toussaint et al. pushed out an integrated platform which aims at helping researchers and clinicians to visualize and process dynamic image data, as well as evaluate simulation results. They used personalized simulation of the Tetralogy of Fallot to illustrate their platform [25]. Geva et al. performed a study using clinical and laboratory data of 100 consecutive patients with repaired TOF (median age 21) to identify independent factors associated . Some other studies also states the correlation between ventricle shape and pulmonary regurgitation indexed with body surface area [28]. Compared to healthy ventricles, the motion and shape features obtained from TOF ventricles also showed high variability and abnormality, indicating the potential use as TOF progression indicators such as RV dilation [29]. Interestingly, similar phenomenon was also found within our results: compared to healthy volunteers, TOF patients RV stress showed higher variability. In our analysis, the morphology features such as RV EDVi and ESVi of TOF patients showed weak capability in post-PVR outcome prediction, which, however, didn't deny the whole usage of ventricle morphology features. In our TOF population, RV EDVI and ESVI in BG and WG had no statistical significant differences (p = 0.937, p = 0.818). More importantly, compared to TOF population discussed in the references [24][25][26][27][28][29], 12 tested subjects in our analysis had much higher RV EDVi and ESVi and lower RV EF, indicating TOF patients in our analysis suffering from more severe RV failure and dilation. The severity in pathology should be considered.

Potential clinical applications
It has been a challenge for cardiovascular researchers to know before-hand which patients would have better post-PVR outcome. Our finding that ventricle stress had the best prediction accuracy (0.8208) in differentiating patients with better post-PVR from patients with worse post-PVR outcome leads to the potential in using mechanical analysis in patient screening and treatment strategies. Knowing that better-outcome patients had lower RV stress level, novel surgical plans could be designed in reducing RV stress level for possible better outcome.
Comparison between healthy and TOF patients revealed that TOF patients have in general higher stress and lower strain levels. That suggested in surgical treatment considerations, the surgeon could introduce novel strategies to reduce ventricle stress and increase ventricle strain (linked to ventricle contractivity) levels. Dr. del Nido proposed a novel PVR surgical strategy to use smaller patch with aggressive scar trimming to improve post-PVR outcome [30,31]. In our clinical trial (NIH 5P50HL074734, Geva and del Nido), 64 patients with repaired tetralogy of Fallot, and who fulfilled defined criteria for PVI were randomly assigned to undergo either PVI alone (n = 34) or PVI with surgical RV remodeling (n = 30). However, no significant difference was observed in the primary outcome (change in RV ejection fraction, -2±7% in PVI alone (Fig 1A), vs. -1±7% in PVI with RV remodeling; P = 0.38) or in any secondary outcomes [31]. Our modeling study indicated that PVR with a smaller patch and more aggressive scar removal led to reduced stress conditions and may lead to improved recovery of RV functions (2-3% EF improvement) [15,32]. However, results from modeling study were with simplified and idealized conditions and need to be validated by large-scale patient studies.

Model limitations and future directions
Limitations in our work include the following: a). Iso-volumic contraction and relaxation were skipped since they involve continuous reference frame change caused by sarcomere zero-stress length change that our current method is not able to handle; b). valve mechanics were not included; c) Patient-specific fiber orientation were unavailable to us. Patient-specific TOF RV/ LV fiber orientations should be acquired and used in our future models; d). a pre-shrink method was applied to in vivo ventricle geometry obtained from CMR data to obtain approximate zero-load geometries and used in our models, which could alter the shape of ventricle when in-vivo pressure was applied. Ohayon et al. and Speelman et al. have developed good method to obtain zero-load geometries which will be considered in our future effort [33,34].
Direct measurements of ventricle in-vivo stress and zero-load sarcomere length are nearly impossible. Ventricle strain measured by some medical devices are often calculated based on some in vivo ventricle reference frames (with either maximum or minimum geometry). Those strains are different from the strain calculated in this paper based on zero-load geometries. Our model is an improvement over the old 1G models [13].