Pulmonary arterial hypertension was induced in rats as described in (Kwan et al., 2021) using the well-established sugen-hypoxia (SuHx) model, an animal model that recapitulates vascular remodeling found in PAH patients (Abe et al., 2010). Male Sprague-Dawley rats (7 weeks old and weighing 214 ± 23 g) were administered with 20 mg/kg of sugen, a vascular endothelial receptor blocker, and kept in 10% O₂ hypoxia for 3 weeks. The animals were then removed from the hypoxic chamber and returned to normoxia (21% oxygen) where the pulmonary arterial pressures continued to rise. Age-matched animals were kept in normoxia during the entire period to serve as the control group. The animals underwent invasive hemodynamic procedures at 4, 8, and 12 weeks post-SuHx induction. The three time points chosen to build three-dimensional biventricular models were based on a study of RV remodeling involving six time points along the disease progression (Kwan et al., 2021). Briefly, we found that after 4 weeks of sugen-hypoxia, rats had significant rise in end-systolic pressures but no changes in ejection fraction, attributed to significant RV hypertrophy. However, later in the disease or rats studied in later weeks of sugen-hypoxia, there were no more changes in hypertrophy. Instead, there was a sharp rise in end-diastolic pressure and end-diastolic elastance with preserved end-diastolic volume. By 12 weeks of sugen-hypoxia, animals show a small reduction in RV end-diastolic elastance, in ejection fraction, and a slight increase in volume. Here we sought to investigate these features. While we note a gradual decrease in the RV ejection fraction over 12 weeks (from 65% in the control to 50%, 47%, and 43% after 4, 8, and 12 weeks), the ejection fraction remained above heart failure thresholds of 35%–40% (Meyer et al., 2010).

Following previously described methods (Vélez-Rendón et al., 2018), all animals underwent invasive open-chest measurements of blood pressure and volume taken in the right and left ventricles while kept under 2.5% isoflurane. Pressure-volume (P-V) timeseries were aligned within the cardiac cycle and averaged. End of systole (ES) was determined by identifying the maximum pressure-to-volume ratio point in the P-V loop. The end-diastolic point was identified as the timepoint in the P-V loop where the pressure was at a minimum and the pressure rate of change (i.e., dp/dt) was at a maximum. After these hemodynamic measurements were taken, the heart was flushed and excised. RV free wall thickness measurements were taken ex vivo across the wall and averaged to obtain a representative RV wall thickness.

Rat-specific meshes were built using data from hemodynamic pressure-volume timeseries and morphological measurements from harvested hearts from four different rats. These rats were selected to represent normotensive and three distinct sugen-hypoxia groups, published in Kwan et al. (2021). The idealized three-dimensional biventricular meshes were built using Gmsh (Geuzaine and Remacle, 2009) and each rat-specific model was parameterized based on the animal hemodynamic and morphological data outlined in Table 1.

We adapted the RV wall thickness of each mesh to match the experimentally measured RV wall thickness and inflated both ventricles to match the cavity volumes presented in Table 1. The meshes were then partitioned into two regions representing the right ventricular free wall (RVFW) and the LV region, which consisted of the left ventricular free wall (LVFW) and the interventricular septum. The morphology, myocardium properties and hemodynamics of the LV were kept constant for all cases, as there was no indication of changes in the left ventricle in any of these animals (Kwan et al., 2021). The muscle-fibre architecture was registered on the meshes using the Laplace Dirichlet Rule-Based algorithm (Bayer et al., 2012). We prescribed the myofiber helix angle to vary transmurally from +60 deg at the endocardium to −60 deg at the epicardium. In Figure 1, we show the meshes built for each case, along with an example of the mesh partitioning and myocardial fiber registration on the control case mesh.

We represent the heart as a continuum body, where the coordinates in the reference configuration ( X ) are mapped to coordinates in the current configuration ( x ) via the deformation gradient: F = I + ∇ u , (1) with u = x − X denoting the displacement of a given point in the domain Ω (i.e., the myocardium). The displacement field is found by solving the quasi-static equilibrium equation given by: ∇ ⋅ P = 0 , (2) where P is the first Piola-Kirchhoff stress tensor, subject to imposed boundary conditions. The basal displacement of our biventricular domain was set to zero in the longitudinal (apex-to-outflow) direction. Basal movement in the other directions, as well as the movement of the epicardial surface, was restricted by a linear spring of stiffness k = 0.5 kPa/cm² as in a previous study (Finsberg et al., 2018b). Measured LV and RV pressures were applied as Neumann boundary conditions at the endocardial surfaces of the domain.

To model the passive behavior of the myocardium, we used the transversely isotropic form of the hyperelastic strain energy function proposed in Holzapfel and Ogden (2009): Ψ F = a 2 b e b I 1 − 3 − 1 + a f 2 b f e b f I 4 f 0 − 1 2 − 1 , (3) where a, a _f , b, b _f are material stiffness parameters, and the invariants are defined as: I 1 = t r C , I 4 f 0 = f 0 . C f 0 , (4) with C denoting the right Cauchy Green tensor and f ₀ the myocardial fiber direction. We assume the myocardium is incompressible and enforce this by adding an extra term p (J − 1) to the strain energy function with p being a Lagrange multiplier representing the hydrostatic pressure, and J = det( F ).

To model the active behavior of the myocardium, we applied the commonly used active stress formulation (Nash and Hunter, 2000): σ = σ p + σ a (5) where σ is the total Cauchy stress tensor which is decomposed into a passive stress contribution: σ p = J − 1 ∂ Ψ ∂ F F T , (6) and an active stress contribution due to the contraction of cardiomyocytes: σ a = T a f ⊗ f + η I − f ⊗ f . (7) The magnitude of the active stress is denoted by T _a and η controls the amount of active stress developed in the directions transverse to the fiber direction, as studies have shown that active stresses in the transverse direction (i.e., sheet and sheet-normal directions) are non-negligible (Lin and Yin, 1998). Similar to the study by Sundnes et al. (2014), we have assumed homogeneity in transverse active stress and consequently set η to a fixed value of 0.2 (or 20%). Notably, the total Cauchy stress tensor σ and the first Piola-Kirchhoff stress tensor P are related by the expression P = J σ ^T F ^−T .

The model was calibrated to match in vivo pressure-volume (P-V) time series data. Measured pressures were provided as input parameters to the model (endocardial boundary conditions), and the model parameters were adjusted until the calculated volumes agreed with the measured ones. The model-data volume mismatch was defined as V R V i − V ̃ R V i V R V i 2 , (8) where V ̃ R V and V _RV are the simulated and measured RV cavity volumes, respectively, and i denotes a specific time point.

Model calibration was carried out in two phases. In the first phase, the passive (isotropic) stiffness parameter a in (3) was estimated by fitting the model to P-V data in the passive filling phase of the cardiac cycle. Due to the sparsity of data used for the optimization, the remaining three material parameters in (3) were not estimated, but fixed to values (b = 5.0, a _f = 2.582kPa, and b _f = 5.0) from Finsberg et al. (2018b) for all simulations in this study. In the second phase, the optimized a parameter was held fixed at its fitted value from the first phase, and the active stress scaling parameter T _a in (7) was estimated by fitting the model to P-V data through the active phases of the cardiac cycle (i.e., isovolumic contraction, ejection and isovolumic relaxation). Since both the P-V values and muscle contraction varies throughout the cardiac cycle, T _a was allowed to vary in time, with a separate value estimated for each measured time point in the cardiac cycle. Figure 2 presents an overview of the parameter estimation pipeline.

End-diastolic elastance (E_ed) and end-systolic elastance (E_es), proposed as global indices of ventricular stiffness (Templeton et al., 1972) and ventricular contractility (Suga and Sagawa, 1974) respectively, were estimated in this study. These metrics were estimated by perturbing the loading conditions on the optimized model while keeping all other variables fixed. Specifically, at the ED point, the ED pressure P _ed was perturbed by incrementing it with a factor (P _ed+Δ = P _ed + ΔP), resulting in a change in ED volume (V _ed+Δ = V _ed + ΔV). The estimate of ED elastance was then obtained by dividing the change in pressure by the change in volume, given as E ̃ e d = Δ V Δ P , (9) where ΔP was set at 0.1kPa. The same approach was used for estimating E_es by perturbing the optimized model at the ES point.

To solve the set of partial differential equations described in (2), we implemented a Galerkin finite element method, which involved discretizing the variational form of (2) using Taylor-Hood tetrahedral finite elements (Hood and Taylor, 1974). Specifically, we used piecewise quadratic basis functions for the displacement field and piecewise linear basis functions for the hydrostatic pressure field. We used a previously developed cardiac mechanics software (Finsberg, 2019) implemented in the FEniCS finite element framework (Logg et al., 2012) to solve the numerical problem. For a detailed derivation of the variational form of (2), interested readers can refer to the work by Finsberg et al. (2018b).

For the minimization of the objective function (8) we used the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm (Broyden, 1970; Fletcher, 1970; Goldfarb, 1970; Shanno, 1970) implemented in the SciPy library (v 1.11.3) (Virtanen et al., 2020) in Python.

The simulated and measured P-V loops for the RV are shown in Figure 3 for the different simulated cases. The results show a very good fit between our model and the experimentally measured data. Time traces of the optimized active stress parameter T _a are also presented in Figure 3. In addition, we present the optimized passive material parameter a for the control, week 4, week 8 and week 12 SuHx cases in Table 2. As previous studies have shown no significant change in LV hemodynamics in early-stage PAH (Kwan et al., 2021), the same P-V data and material parameters were used for the LV in all simulated cases (Figure 4; Table 2). The LV active and passive material parameters were calibrated using P-V data from the control animal.

We conducted a mesh convergence analysis based on the control geometry and hemodynamics to find the optimal mesh resolution needed for accurate model predictions. Specifically, we calculated the average Cauchy stress and Green strain employing four distinct mesh resolutions, ranging from a low-resolution mesh comprising approximately 5,500 elements to a high-resolution mesh comprising approximately 61,000 elements. As depicted in Figure 5, the derived metrics exhibited a very low sensitivity to mesh resolution beyond the medium-high resolution (36,473 elements) threshold. This observation implies that a mesh resolution of 36,473 elements is adequate to ensure the accuracy of our model predictions. The chosen mesh sizes for the control, 4, 8, and 12 weeks post-PAH analyses are detailed in Table 3, including the respective average evaluation time for the cost functional (8), the number of cost-functional evaluations to fit one P-V point, and the total run times for each optimization process. All computational analyses were conducted on a computing cluster utilizing a single node with 32 cores.

End-diastolic elastance (E_ed) and end-systolic elastance (E_es) serve as valuable indices for evaluating ventricular stiffness and contractility, respectively (Templeton et al., 1972; Suga and Sagawa, 1974). The gold standard for computing these metrics is by transiently varying RV preload, generating multi-beat pressure-volume loops subsequently used for E_ed and E_es measurements (Suga et al., 1973; Maughan et al., 1979). This method, while also applied in clinical studies (Dell’Italia and Walsh, 1988; Hsu et al., 2020), is less established in the clinical setting compared to experimental studies. Here, we implemented a different approach that is comparable to the gold standard. We altered RV loading conditions by perturbing end-diastolic and end-systolic pressure, keeping all other parameters fixed. We then computed E_ed and E_es as the slope of the resulting pressure-volume relationship, i.e., ΔP/ΔV. This approach has been used previously for calculating E_es (Finsberg et al., 2018a; Finsberg et al., 2018b). We observed that the model predicts a slight increase in diastolic stiffness from week 8 to week 12, although the material stiffness parameter displayed in Table 2 is reduced. This apparent inconsistency is most likely the result of the increased wall thickness at week 12, which directly impacts the overall chamber elastance.

Nevertheless, the general trend of elevated stiffness in disease compared to control aligns with previous measurements in sugen-hypoxia (Kwan et al., 2021), monocrotaline (Vélez-Rendón et al., 2018), and pulmonary artery banding (Rain et al., 2016) animal models of PAH. Table 4 displays results from the study by Kwan et al. (2021), supporting this general trend. Likewise, Rain et al. (2016) demonstrated increased RV stiffness in rats with mild and severe RV dysfunction, attributing the rise in mild dysfunction to myofibril-mediated stiffness and in severe dysfunction to both increased myofibril stiffness and fibrosis. In a monocrotaline rat model, Vélez-Rendón et al. (2018) reported stiffening of the passive myocardium after 4 weeks post-PAH induction, observing an initial decrease in stiffness at week 1, possibly due to a temporary increase in myocardial compliance required to preserve RV stroke volume. However, by 2 weeks post-PAH induction, passive stiffness began trending upward. Our results also align with human clinical (Rain et al., 2013; Trip et al., 2015; Hsu et al., 2018) and computational (Finsberg et al., 2019) studies of PAH. In all four studies, they noted a progressive increase in RV passive stiffness in PAH patients compared to controls. This increase was found to be influenced by the degree of remodeling or the severity of the disease.

In our study, we observed an initial twofold increase in contractility from control to week 8, followed by a decrease at week 12, indicating a downregulation of RV contractility by 12 weeks post-PAH induction. This initial increase in contractility, reported in both animal (Blaudszun and Morel, 2012; Vélez-Rendón et al., 2018; Kwan et al., 2021) and human (Rain et al., 2013; Hsu et al., 2018; Finsberg et al., 2019) PAH studies, is linked to adaptive hypertrophy (Bogaard et al., 2009; Kwan et al., 2021), believed to preserve systolic function in early PAH stages (Naeije et al., 2014; Naeije and Manes, 2014). However, increased ventricular contractility is unsustainable, and a downregulation of contractile force is commonly observed as the disease progresses (Fan et al., 1997; Bogaard et al., 2009). Our model results indicate that this downregulation starts after 8 weeks post-PAH induction in a SuHx-rat model. Few PAH studies have explored the time course of contractility changes; most only distinguish between control and PAH subjects. In a relevant study, Finsberg et al. (2019) noted a 20% initial increase in RV contractility in mild RV remodeling, later down-regulating in severely remodeled cases to values below the control. Their study, focusing on human subjects, did not specify the time-course of the significant decrease in contractility but distinguished between mild and severe remodeling based on RV end-diastolic volume to LV end-diastolic volume ratios. Similarly, Kwan et al. (2021) reported downregulation at 10 weeks post-PAH induction in a sugen-hypoxia animal model, with values still higher than control at the 10-week timepoint. This suggests a peak in RV contractility at 8 weeks post-PAH induction, followed by a progressive downregulation, aligning with our observations.

Ventricular wall stress, particularly at end-systole (ES) and end-diastole (ED), plays a pivotal role in systolic and diastolic cardiac function. Its significant correlation with myocardial oxygen consumption and adverse cardiac remodeling has been well-established by previous studies (Strauer, 1979; Bogaard et al., 2009; Alter et al., 2016; Haque and Wang, 2017). Pressure overload directly impacts wall stress, adhering to the law of Laplace, and an elevation in wall stress can hinder myocardial oxygen supply by compressing coronary circulation (Chin et al., 2005). Consequently, deviations from normal or baseline wall stress can detrimentally affect oxygen availability to cardiomyocytes, potentially leading to cardiac ischemia, adverse remodeling, and ultimately, heart failure.

The inverse linear correlation identified by Quaife et al. (2006) between ES wall stress and RV ejection fraction, supported by Alter et al. (2016) for the LV, emphasizes the critical relationship between wall stress and cardiac function. This underscores the potential of wall stress as a diagnostic index for evaluating heart performance in disease. However, the lack of a direct method to measure RV wall stress necessitates reliance on mathematical approximations, with several proposed methods in the literature and no universally recognized gold standard. Consequently, the accuracy and agreement of these approximations carry significant experimental and clinical implications.

Our model results demonstrate a progressive increase in ED and ES fiber stress in disease compared to control. Circumferential and apex-to-outflow stress also exhibit an increase at ED and ES. Comparison with existing literature that reported wall stress in PAH, consistently revealed increased RV wall stress at ED and ES in disease compared to control (Quaife et al., 2006; Vélez-Rendón et al., 2018; Gold et al., 2020; Kwan et al., 2021), despite variations in geometric model assumptions and stress computation methods across these studies. Vélez-Rendón et al. (2018) and Kwan et al. (2021) assumed a spherical RV geometry and used the thin-walled Laplace law to estimate fiber stress in the RV free wall. Quaife et al. (2006) also assumed a spherical RV but employed a modified Laplace law suitable for non-circular cross-sections (Janz et al., 1989). Gold et al. (2020) utilized a two-dimensional model based on short-axis echo images and estimated wall stress using the von Mises formula, providing an equivalent stress over the entire RV rather than individual stress components. Finsberg et al. (2019), utilizing a modeling technique and stress approximation method similar to our study, did not report ED and ES wall stress, but instead reported peak fiber stress, which does not necessarily coincide with end-diastole or end-systole. Interestingly, they reported an increase in peak fiber stress only in severely remodeled cases which they defined as a ratio of RV ED volume (RVEDV) to LV ED volume (LVEDV) greater than 1.5. In the mildly remodeled case (RVEDV/LVEDV ≤1.5), peak fiber stress was at control level. This finding aligns with our observation that peak fiber stress, which occurred at a time point between ED and ES in our model, only surpassed control values at 12 weeks post-PAH induction, although it exhibited an increasing trend in the disease cases (Figure 6). This suggests that peak fiber stress increases as PAH progresses, and RV geometric remodeling transitions from wall thickening to dilation.

At the onset of PAH, both geometric and material remodeling occur simultaneously, and manifest as compensatory mechanisms. Our study aimed to isolate these remodeling mechanisms and explore their individual impacts on RV mechanics, as indexed by wall stress and strain. The results revealed that material remodeling has a minimal effect on wall stress and strain, while geometric remodeling plays a predominant role in reducing these parameters as the disease progresses. This agrees with existing knowledge that attributes the primary role of geometric remodeling to the reduction of wall stress (to normal physiological values) in the presence of pressure overload (Grossman et al., 1975).

Our study provides quantifiable insights into these established findings. Specifically, we found that if the RV wall had not thickened, ED fiber stress would have exhibited a substantial increase, reaching six-fold by week 4, eleven-fold by week 8, and a remarkable twenty-fold increase by week 12 compared to baseline (control) levels. However, due to observed geometric remodeling, the maximum increase in ED wall stress was limited to seven-fold, occurring at week 12. Similarly, ES wall stress would have experienced a three-fold increase by week 4 and a five-fold increase at weeks 8 and 12 without RV wall thickening. Again, the presence of RV wall thickening significantly mitigated these increases, with ES wall stress nearly fully normalized by week 12.

From these findings, two key takeaways emerge. Firstly, while geometric remodeling alone fell short of fully normalizing wall stress in the 12 weeks time course of this study, it effectively mitigated the majority of stress increases associated with rising RV afterload. Secondly, the presented magnitudes of wall stress increase, in the absence of RV thickening, shed light on what can be expected when RV geometric remodeling transitions from wall thickening to wall thinning in the later stages of PAH (in an attempt to maintain stroke volume and cardiac output). These excessive wall stress levels may lead to a detrimental sequence of cardiac ischemia, further compromised RV contractility, causing additional RV dilation, and ultimately culminating in RV dysfunction.

There are a number of important limitations to the present study, which should motivate further experimental and computational research. We applied the computational framework to a limited cohort of four animals, which included one normotensive and three hypertensive animals at distinct time points. Although the animals chosen for this study were representative of the control and disease groups, they do not represent a full description of the longitudinal and progressive remodeling of the RV in PAH. While not comprehensive, the study aimed to analyze mechanical changes in the RV during the initial 12 weeks post-PAH induction. Future research could apply the framework to a larger cohort over an extended period for a better understanding of RV behavior in PAH progression.

We only fit a single passive material parameter, while the others are set to values from the literature. This modeling decision was necessitated by the limited data available for the optimization. Attempting to fit all four parameters in (3) to in vivo P-V data would have substantially increased the computational cost. More importantly, it would have resulted in a scenario where multiple sets of parameters would minimized the cost function (8), rendering the chosen optimal parameter set overly sensitive to the initial guess provided to the optimization algorithm. While our approach aligns with previous research suggesting that a unique solution can be achieved by optimizing a single parameter (or at most two of the parameters) of (3) when fitting to in vivo P-V data (Hadjicharalambous et al., 2015; Balaban et al., 2017), it is crucial to recognize its inherent limitations, specifically the dependence on parameter values sourced from the literature. More detailed experimental data, for example, local ventricular strains estimated from magnetic resonance images, or biaxial stress-strain data extracted from myocardial tissue mechanical testing will be important for more accurate characterization of material properties, as well as for validation of the model results. As shown in Figure 3, the passive filling phase - which is the part used for optimizing a _RV - is nearly flat for several of the P-V loops used in this study. As such, small deviations in the pressure can lead to large variations in the fitted parameters. In addition, our decision to fit the isotropic material parameter in 3, while keeping the parameter that describes stiffness in the fiber direction, a _f , constant in both control and disease cases, may have contributed to the absence of a discernible trend in active stress development, as depicted in Figure 3B. To address these issues, future model development should be based on more comprehensive datasets which should enable improved characterization of the passive tissue stiffness.

Due to the unavailability of imaging data of sufficient resolution to build meshes for the rats used in this study, we employed idealized biventricular geometries in the mechanical model. Despite this limitation, the idealized geometry gives insights on the effect of wall thickening on RV mechanics with PAH progression. It is important to note that the altered septal wall positions depicted in Figure 1 as the disease progresses stems from the methodology used in generating the biventricular meshes. Specifically, it is a result of the inflation of the mesh geometries to match the initial cavity volumes outlined in Table 1. While the magnitude of septal wall flattening was not validated against imaging data, qualitatively, it is consistent with cardiac magnetic resonance imaging data obtained from a different rat at time points similar to those used in this study.

We modeled the myocardium as transversely isotropic. The passive myocardium is inherently complex, displaying strong non-linearity and anisotropy, and would require an orthotropic model to fully characterize its mechanical behavior. However, it has been documented that a transversely isotropic formulation of the Holzapfel-Ogden law strikes a good balance between parameter identifiability–the ability to determine a unique parameter set given limited amount and quality of the data–and model fidelity–the ability of the model to adequately represent cardiac deformation and function (Gjerald et al., 2015; Hadjicharalambous et al., 2015). Given that the data available for model fitting in this study was quite limited, and the extensive use of a transversely isotropic model to approximate orthotropic cardiac tissue properties in the literature (Guccione et al., 1991; Xi et al., 2011; Balaban et al., 2017; Finsberg et al., 2018b; Avazmohammadi et al., 2019), we deemed the use of this model to be appropriate.

We did not include rat-specific measurements of fiber orientation or regional strain. Instead, we implemented the same fiber orientation across all simulations using a rule-based method and assumed a transmural variation in the fiber direction. This modeling choice is a candidate for further refinement as some studies have reported changes in myofiber architecture in PH (Avazmohammadi et al., 2017; Mendiola et al., 2023). However, it remains unclear to what degree myofiber angle remodeling occurs and whether previous studies are confounded by the animal model, or the stage of disease. In addition, our assumption of transmural (through-thickness) variation in the fiber direction is a common assumption in computational studies of cardiac mechanics that employ a rule-based method to assign myocardial fibers to the mesh geometry and has been validated by histology (Streeter Jr et al., 1969; Hill et al., 2014) and DTI-based studies (Hsu et al., 1998; Holmes et al., 2000; Scollan et al., 2000) of cardiac fiber architecture. Nevertheless, future research could benefit from incorporating diffusion tensor magnetic resonance imaging (DT-MRI) data from Sprague-Dawley rats to create a more accurate representation of the fiber architecture.

Finally, while isoflurane and other anesthetics can alter cardiovascular function (Loushin, 2005), isoflurane was carefully regulated at 2% (in 100% O₂), below reported doses that drastically alter heart rate and hemodynamics (Yang et al., 2014; Pang et al., 2018). The effects on hemodynamics due to isoflurane were likely small compared with the differences due to the PAH treatment, as seen previously (Hill et al., 2014; Vélez-Rendón et al., 2019; Kwan et al., 2021).