The workflow of the model characterization is outlined in Figure 1. In summary, the cardiac model is customized to an individual rat heart by first fitting a generic template mesh to a segmented three-dimensional DT-MRI image. The primary eigenvectors of the diffusion tensor at each image voxel specify the fiber directions, which are fitted and incorporated into the mesh. Contraction is then stimulated by applying an intracellular calcium transient simultaneously and ubiquitously throughout the muscle tissue. The numerical solution of the ordinary differential equations describing the system mechanics yields the pressure-volume characteristics over a full heart beat, that can be compared directly to hemodynamic and echocardiographic measurements performed on the individual rats. The tissue stiffness, aortic pressure, and calcium-binding sensitivity are then tuned to achieve maximal consistency with the rat-specific measurements.

In some instances, it was not possible to reproduce the measured diastolic and contractile properties using any choice of model parameters because the diameter of the segmented heart was too similar to the measured end-diastolic diameter. An analysis of these hearts indicated the likely presence of residual LV pressure during the ex vivo fixation process, performed prior to DT-MRI, resulting in an overestimate of the LV dimensions in diastasis, and hence precluding reliable model parameterization. Those hearts were therefore discarded. In total, one SHAM heart and three AB hearts were successfully simulated and considered for further analysis.

Having thus parameterized the heart models, we investigated the impact of aortic banding on cardiac function. Conceptually, our methodological approach associates the remodeling process with a hypothetical trajectory through a multidimensional parameter space U with coordinates that describe the cardiac system:

where c1˜ measures the passive tissue stiffness, [Ca²⁺]₅₀ is the calcium-binding sensitivity of troponin, p_a is the aortic pressure, z is the aortic-valve impedance, PCa and DCa are respectively the peak and diastolic intracellular calcium concentrations, and u_geom and u_fiber are generalized coordinates representing the geometrical and fiber configurations. From this perspective, the end points of the trajectory represent the real heart models (i.e., the individual SHAM and AB hearts), while the trajectory represents the continuous evolution from a “healthy” heart (SHAM) to the remodeled heart, as a result of aortic banding. At each intermediate state considered along the trajectory, a new hypothetical heart model was fitted again using the process outlined in Figure 1 and assuming linearly interpolated values for the phenotypes.

The geometrical and fiber parameters u_geom and u_fiber, which represent the geometrical and fiber structures, are defined as continuous variables bounded by ugeom(SHAM)=ufiber(SHAM)=0 and ugeom(AB)=ufiber(AB)=1. Other points along the trajectory are expressed by intermediate values of u_geom and u_fiber, representing hypothetical heart structures generated by linear interpolation of, respectively, the node coordinates or the fiber angles of the end-point models.

Within this formalism, we quantified the “compensatory ability” of a given model, with respect to each variable, in terms of the sensitivity of EF, considered as a function over U:

This formulation ensures the non-dimensionality of S and considers changes in a given parameter u in proportion to the phase-space distance spanned by the overall trajectory. The derivatives ∂EF(u)/∂u were calculated by applying changes to u that were sufficiently small to remain in the local linear regime around the model considered. As discussed below, we interpret S[EF, u] as a measure of the ability of each parameter u to control EF and maintain the heart in a compensated state.

An initial cohort of 20 Sprague-Dawley rats was initially considered for echocardiography (University of California, San Diego). Eight weeks (±4 days) after birth, 10 of the rats were applied an AB on the transverse aorta with a 0.5 mm clip, while the other 10 underwent a sham surgery. All animal use followed NIH guidelines and was approved by the Institutional Animal Care and Use Committee (IACUC) at the University of California, San Diego.

The heart rate of each rat was monitored every 7 days post intervention. Echocardiograms were measured immediately prior to the operation, and subsequently after 2 and 4 weeks. The rats were rested on an IACUC-approved heating pad and anesthetized with 1.5% isoflurane in 100% O₂ by nose cone. M-mode echocardiographic images were acquired using a GE vivd/I imaging system. The entire procedure took approximately 15–20 min per animal. The animals were transferred to their cages and quickly regained full consciousness.

Four weeks post operation, the rats underwent a hemodynamic study to determine the maximum aortic pressure during LV ejection using a manometer-tipped catheter (Millar, Houston, TX, US). A subset of six rats from the initial cohort were then imaged by diffusion-tensor MRI (DT-MRI). Two of the SHAM and four of the AB hearts were extracted and fixed in diastasis by arrest using a high-potassium solution, followed immediately by fixation to approximate the relaxed state. Following a method described elsewhere (Benson et al., 2011), the muscle-fiber orientations were determined from DT-MRI measurements performed on these fixed hearts.

M-mode echocardiograms were measured in vivo along the vertical diameter of a basal slice, displaying the time course of the LV inner and outer walls over five heart beats (Figure 2A). Mean diastolic and systolic LV diameters (LVD) and wall thicknesses (LVWT), were determined for each rat at time points t = 0, 2, and 4 weeks from the AB or SHAM intervention. To assess the potential impact of aortic banding, the average rate of change of these parameters y, over the 4-week period, was quantified using the metric

The EF for each heart was estimated by treating the M-mode echocardiography images as cross sections through a cylindrical LV slice (Figure 2A). The change in the internal volume of these cylinders reflects the ventricular EF in so far as the blood ejected from the basal region of the LV constitutes the majority of the total ejection. The incompressibility of the wall tissue implies that any change in the cross-sectional area of the cavity wall is accompanied by an extension of the cylinder height h in the base-apex direction to conserve the tissue volume π(r + LVWT)²h − πr²h, where r = LVD/2. Hence,

To validate this approach, we compared the result of Equation (4) with a more explicit measurement of the rat LV volume using MRI data available in the literature (Wise et al., 1998). The good agreement in the two approaches confirmed the suitability of Equation (4) for our present purpose. The details of this analysis are presented in Appendix 1.

Diffusion-weighted images were acquired ex vivo in 15 non-collinear directions on a 7 T scanner equipped with Avance II hardware (Bruker) using a DTI-EPI sequence. A gradient insert allowed for a 1 T/m maximum gradient strength and maximum slew rate of 11,250 T/m/s. Imaging parameters were as follows: FOV 20 × 20 × 20 mm, isotropic 200 μm resolution, gradient duration 2 ms, gradient separation 11 ms, nominal b-value 500 s/mm², number of non-diffusion weighted images 1. The total scan time was approximately 1 h. This scan protocol was repeated five times and the resultant images were manually averaged in AFNI (NIH, Bethesda, MD).

The resulting image stacks were segmented (ITK-SNAP) and used to reproduce the LV geometry and muscle-fiber configuration in the computational model, as explained below.

We complemented the above imaging results with measurements of intracellular calcium stimuli, performed in separate experiments at physiological temperature (37°C) in single myocytes derived from sham and aortic-banded Wistar rat hearts 6 weeks following the the placing of a constriction on the ascending aorta (Oslo University, Norway) (Røe et al., 2017). Samples from the aortic-banded rats required that the hearts display wall thickening in excess of 1.9 mm with no signs of systolic heart failure (left atrial diameter less than 5.0 mm and no left ventricular dilation). The overall phenotypes of the animals were similar to those described in section 2.2. The EF was preserved in the AB rats, compared to SHAM-operated rats, consistent with the rats in section 2.2.

LV myocytes were isolated using a standard enzymatic dispersion technique (Louch et al., 2011). Calcium measurements and calibrations were done as described by Røe et al. (2017).

The cardiac cycle was simulated using the protocol described in detail elsewhere (Land and Niederer, 2011; Land et al., 2012) to reproduce the heart cycle (pressure-volume loop in Figure 1). Briefly, a phenomenological model of diastolic filling governed LV inflation toward a set end-diastolic pressure of 1 kPa. Upon stimulating the muscle with a calcium transient, the LV volume was fixed to simulate isovolumetric contraction. The calcium transient, combined with the myocyte contraction model, generated active tension within the tissue along the fiber directions, resulting in a sharp LV pressure rise. Upon reaching the pre-set aortic pressure p_a, the Windkessel model governed blood ejection from the LV until volume flow was reversed. Fixing the LV volume again in this state, isovolumetric relaxation was then initiated, terminating when the LV pressure attained the chosen diastolic pressure (1 kPa). Diastolic filling was then reactivated and the process was repeated until the pressure-volume loop converged to a steady-state limit cycle.

To minimize the reliance of the simulation on unmeasured parameters, we limited the analysis to consider a simplified phenotype space defined by three variables: the mean LV diameter at end diastole (LVD, measured from the echocardiography data), the maximum pressure p_max during LV ejection (obtained from the hemodynamic measurements), and EF (determined from the echo data using Equation 4). Each of these quantities can be compared directly with the simulation output (section 2.7).

Preliminary simulations were conducted to explore the parameter space. Variations in model parameters yielded smooth monotonic changes in the phenotypes, thereby justifying the use of a gradient-base method of solution. A Newton method was therefore used to fit each rat LV model, to reproduce the three corresponding phenotypes, in terms of three model parameters: (1) the effective tissue stiffness c1˜, which predominantly constrains the end-diastolic volume; (2) the aortic pressure p_a; and (3) the dissociation constant [Ca²⁺]₅₀ between Ca²⁺ ions and troponin C, which controls the extent of the contraction (Land et al., 2012). Convergence of the solution was generally achieved after one or two iterations and the results showed no significant dependence on initial values.

In the absence of experimental data relating to diastolic filling, we assumed a constant diastolic pressure p_ed = 1 kPa throughout the filling phase. Consequently, the parameter c1˜ encapsulates both the physical stiffness c₁ (Equation 5) and the assumed p_ed, i.e., c1˜~c1/ped.

Figure 3B shows a representative mesh constructed from an LV segmentation (Figure 3A). The fiber unit vectors, computed at each voxel contained within the mesh domain, were projected onto the local basis vectors of the mesh to yield elevation (α) and imbrication angles (β). The cross-sectional maps of α and β shown in Figures 3C,D are consistent with histological observations of fibers being predominantly parallel to the cavity walls (|β| ≲ 10°), and with α varying linearly from ~ +50° at the endocardial wall to ~ −40° at the epicardial (Figure 3E). For the purpose of the simulations, we therefore imposed β = 0 throughout the tissue, leaving α as a unique scalar field characterizing the fiber configuration. The fiber vectors were determined at the mesh nodes by cubic-Hermite fitting (Figure 3B).

Our main results, summarized in Figure 4, identify an evolution in the subset of cardiac properties with the greatest ability to regulate EF. Whereas the EF sensitivities associated with the LV geometry and tension generation (via Ca²⁺ binding to troponin C) are of similar magnitude in the hearts with no aortic banding, the latter sensitivities increase significantly at the expense of the former following aortic banding. In other words, a healthy heart can be expected to compensate better for the onset of hypertension by a thickening of its cavity wall, than a heart in an advanced state of concentric hypertrophy; at that stage, remodeling of calcium signaling becomes a more effective compensatory mechanism. The complementary simulations summarized in Figure 5 and Appendix 2 further identify the wall thickness as the variable that most significantly governs the evolution in compensatory ability observed in Figure 4. The consistency of the observed behavior, over a range of phenotype values, supports the generality of the effect, beyond the few experimental cases available in the present study.

Another natural concern for this interpretation is that the hypertrophic process necessarily occurs much more slowly than the sudden impact of aortic banding. However, our simulations probed only the experimentally accessible end points of the overall remodeling process. Therefore, insofar as aortic banding serves as an experimental model for the more gradual onset of hypertension, our interpretation remains independent of the evolution of the rat hearts between these end points.

Concentric hypertrophy, clearly exhibited in the wall thickening of the AB hearts, has been interpreted as a passive compensatory response to pressure overload, insofar as it reduces the wall stress in the presence of enhanced LV pressure (Laplace's law) (Grossman et al., 1975; Olivetti et al., 1988). Its supercession by remodeling of the calcium dynamics, through PCa and [Ca²⁺]₅₀, 6 weeks after the intervention, marks a qualitative transition in the heart's capacity to achieve compensation. The peak calcium level PCa encapsulates the underlying electrophysiological machinery of the muscle cells, which expends energy at every heart beat. The increase in S[EF, PCa] in the AB hearts can be explained by the lowering of PCa (Figure 3F). In other words, a given absolute change in calcium becomes proportionally more significant in an AB than in a SHAM heart. Our dynamic model fits indicate a concomitant decrease in [Ca²⁺]₅₀, the intracellular Ca²⁺ concentration required to produce 50% occupancy of the troponin C regulators of actomyosin cross-bridge formation. Biochemically, [Ca²⁺]₅₀ is controlled by the phosphorylation of troponin I through various kinases (Layland et al., 2005). There was no a priori reason for the [Ca²⁺]₅₀ values listed in Table 2 to vary in the same manner as the measured PCa values (1.53 μM for SHAM-heart cells, 0.93 μM for AB). The model fitting, however, implies a strong correlation between those parameters that arguably reflects a functional benefit by maintaining contraction dynamics near the middle of its dynamic range.

The observed depression of the Ca²⁺ transient in cases of heart failure has been associated with impaired contractile behavior (Siri et al., 1991). However, when accompanied by a proportional decrease in [Ca²⁺]₅₀, this change may benefit compensation by magnifying the impact of small changes in the calcium concentration on EF. This evolution may become increasingly beneficial with the decline of compensation through hypertrophy.

By limiting the scope of our simulations to the principal features of diastole and systole, we reduced the parameter space to maximize the objectivity of the model parameterization. We lumped the passive mechanical tissue properties into the single phenomenological parameter c1˜, amalgamating the physical stiffness c₁ and the unknown end-diastolic pressure p_ed. The measured end-diastolic diameter was hence uniquely fitted. In systole, tension generation is controlled by both the intracellular calcium dynamics and the strength of crossbridge contraction. Aortic banding can in principle affect either—or both—of these contributions. As outlined in section 2.6.3, contractile strength is itself influenced by several mechanisms, whose detailed characterization lies outside the scope of the available measurements. We therefore chose to lump the sensitivity to tensile behavior in the calcium-binding parameter [Ca²⁺]₅₀, and assumed a constant T_ref in all the hearts. This simplification assumes that an impact on T_ref by aortic banding would likely involve a substantial remodeling of the myofibril structure. On the other hand, a change in [Ca²⁺]₅₀ through the modified expression of kinases arguably involves a simpler regulation, at the protein level, that can more reasonably be expected to occur in the early stages of hypertrophy. Variations in calcium sensitivity in cases of heart failure have been reported, as a result of the phosphorylation of myofibrillar regulatory proteins (Wolff et al., 1995; Velden et al., 2009).

DT-MRI is a costly and delicate technique that affords a considerably lower throughput than echocardiography. Thus, whereas the latter experiments provided statistically significant samples of geometric measurements, only a small proportion of those hearts underwent DT-MRI processing. As explained above, the main appeal of performing both experiments successively on the same hearts was to increase the consistency of the models with the biological systems. One performance criterion for the computational simulations was their ability to reproduce the morphological dynamics inferred from the echocardiograms (see section 2.8). Of the initial cohort of hearts (ten SHAM and ten AB), two SHAM and four AB hearts proceeded to DT-MRI. Within this subgroup, only one SHAM heart and three AB hearts successfully simulated the expected diastolic size and EF. In both of the unsuccessuful cases, blood ejection into the aorta stalled at a fraction of the expected EF, even when the strength of contraction (T_ref) was increased in the extreme. We deduced that this failure to achieve the expected dynamics resulted from the nominal LV size in diastasis (derived from the DT-MRI measurements) being almost identical to the diastolic diameter (to within 10% or less). In contrast, the successfully simulated hearts had an equilibrium size that typically amounted to 60–70% of the end-diastolic volume, which provided scope for sufficient contraction. One potential cause of error in the erroneous DT-MRI measurement is the possibility of even a small intramyocardial pressure being exerted during perfusion with the high-K⁺ arrest solution and fixative. The near coincidence of the LV diameter in this state with the measured end-diastolic size is consistent with an unintended inflation of the LV. We therefore considered these anomalies as artifacts and discarded those hearts from the analysis, retaining one SHAM and three AB hearts.

The model characterizations would no doubt have benefitted from additional measurements, for example of tissue deformation by MRI or echography. However, these measurements were not part of the initial experimental protocol. Although future exploration of these parameters would be desirable, we evaluated LV torsion in the simulated contractions (angle of twist per unit length of ventricular wall) and found no systematic trend distinguishing the SHAM and AB hearts. The available measurements thus provide no evidence that this is a significant parameter for the present purpose.

To increase confidence in the validity of our conclusions despite the paucity of available hearts, the complementary simulations of Figure 5 and Appendix 2, done with the canonical semi-ellipsoidal heart meshes, explored a range of the phenotype space around the SHAM heart (see Table 1: SHAM). All the cases considered display the same qualitative behavior, with the sensitivity S[EF, LVWT] decreasing with increasing LVWT while S[EF, PCa] and S[EF, [Ca²⁺]₅₀] both increasing. Taken together, these results demonstrate the general pattern that, as LVWT increases, EF becomes less sensitive to variations in LVWT and more sensitive to those of [Ca²⁺]₅₀ and PCa.

To our knowledge, the present work constitutes a first attempt to provide a global view of cardiac compensation, within the context of ventricular pressure overload, by considering the anatomical shape and the passive and active mechanical behavior on an equal footing. We conclude that intracellular behavior becomes a more efficient compensatory mechanism than organ-level remodeling as the rat heart responds to the onset of pressure overload. This phenomenon was observed with regard to EF, a traditional indicator of cardiac dysfunction routinely measured in clinical protocols. Although other metrics of cardiac function (e.g., cardiac output and the stroke volume per body weight) directly reflect cardiac efficacy, those metrics are notoriously dependent on temporal and biological factors, and may vary significantly among healthy subjects of the same species, or even, over time, in a given individual (Haxhe and Lammerant, 1964; Ovsyshcher et al., 1993; Collis et al., 2001). Such incidental influences are arguably normalized in EF, allowing a more equitable comparison of different hearts. Furthermore, the EF is preserved at its healthy level in nearly half of heart-failure patients (Sharma and Kass, 2014). One major obstacle in fully understanding this increasingly prevalent condition is the limitation of existing experimental models for highlighting the distinct facets of the disease. Toward meeting this challenge, computational modeling provides an invaluable basis for unifying experimental results and investigating subsystems within an integrative perspective.

AM and JO are co-founders, equity holders and scientific advisors to Insilicomed Inc., a licensee of UCSD software used in this research. This relationship has been disclosed to the University of California San Diego and is overseen by an independent conflict of interest management subcommittee appointed by the university. The other authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.