Model order reduction for left ventricular mechanics via congruency training

Summary

Computational models of the cardiovascular system and specifically heart function are currently being investigated as analytic tools to assist medical practice and clinical trials. To achieve clinical utility, models should be able to assimilate the diagnostic multi-modality data available for each patient and generate consistent representations of the underlying cardiovascular physiology. While finite element models of the heart can naturally account for patient-specific anatomies reconstructed from medical images, optimizing the many other parameters driving simulated cardiac functions is challenging due to computational complexity. With the goal of streamlining parameter adaptation, in this paper we present a novel, multifidelity strategy for model order reduction of 3-D finite element models of ventricular mechanics. Our approach is centered around well established findings on the similarity between contraction of an isolated muscle and the whole ventricle. Specifically, we demonstrate that simple linear transformations between sarcomere strain (tension) and ventricular volume (pressure) are sufficient to reproduce global pressure-volume outputs of 3-D finite element models even by a reduced model with just a single myocyte unit. We further develop a procedure for congruency training of a surrogate low-order model from multi-scale finite elements, and we construct an example of parameter optimization based on medical images. We discuss how the presented approach might be employed to process large datasets of medical images as well as databases of echocardiographic reports, paving the way towards application of heart mechanics models in the clinical practice.

Keywords:

Simulation and modeling – magnetic resonance imaging – Muscle cells – Hemodynamics – heart failure – Ejection fraction – Cardiac ventricles – Finite element analysis

Introduction

Multi-scale finite element (FE) models of cardiac mechanics are being investigated as novel tools for analyzing medical imaging data and assisting personalized diagnostics in cardiac resynchronization therapy and disease monitoring [1–4]. Despite recent advances, personalizing FE simulations of heart mechanics to a specific clinical case still constitutes an open research problem. While most state-of-the-art models can capture the anatomical details of atria and ventricles from 3-D medical images (e.g., from CT or from the more costly MRI [5, 6]), significant computational resources are necessary to inversely estimate, when at all possible, the many unknown model parameters. Applying FE models in large scale analyses of clinical trials is, therefore, currently intractable. At the same time, it is widely accepted that a more mechanistic understanding of trial results could be greatly beneficial. For example, direct applications of models could provide illustrations of certain biophysical correlates of trial trends, thus aiding the interpretation of complex drug effects, such as intropes on cardiac function in heart failure patients [7–9].

Clinical evaluations often rely on echocardiography as their main source of cardiac functional and anatomical assessments, with analysts typically not having access to the raw acquisitions from the study, but rather relying only on succint clinical reports summarizing the measurements as a list of discrete indices. To more completely leverage information from echocardiography, multi-scale models would additionally exploit simple descriptions of anatomy provided by these limited ultrasound assessments. In such a framework, anatomical parameters would be estimated, and uncertainty quantified statistically by sampling (e.g., in a generalized polynomial chaos framework [10]). Again, due to complexity, the large computational expense connected to statistical approaches cannot currently be absorbed into the clinical application when FE models are employed directly. Instead in this study, we lay groundwork for a more feasible approach employing model order reduction, and a novel strategy to reduce 3-D FE models is the subject of this paper.

When using state-of-the-art biophysical models of heart, simulating even just a few beats can require hours to days of supercomputer time. Simplified models are therefore desirable, but they typically lack accuracy or interpretability. A clear example is provided by image-based models of the heart (see [11, 12] for two recent reviews), for which low-order analogues have typically relied on restrictive assumptions that fail to accurately represent the anatomy of ventricles and atria [13–16]. The lack of a clear bridge between imaged anatomies and model representations makes it difficult to link unknown parameters to radiological or physiological measurements. Thanks to basis simplifications introduced directly at runtime, advanced numerical techniques such as proper orthogonal decomposition [17–20], reduced basis functions [21, 22], and hyper-reduction [23] are able to speed up high-resolution computations by simplifying both the temporal and spatial solution spaces while at the same time maintaining a strong correspondence with their computationally expensive counterparts. These techniques provide significant advantages as they can rely on strong theoretical proofs of congruency to high resolution simulations, and also because they do not require design changes such as modifications to the base model’s parameterization. Obstacles to their adoption have been, however, their implementation complexity and the relatively contained speed-ups that such techniques can achieve. Although less rigorously, further speedups can be achieved by selecting a priori the subset of main global deformation modes assumed to drive contraction, and then by restricting the solution space accordingly [24].

Alternative strategies have proposed training purely statistical models (e.g., Gaussian Process regressions [25]) on a relatively limited number of high-resolution runs (e.g., see [26, 27] for applications in diastolic filling). Once fitted to a sampled training set, statistical models can be used as efficient surrogates for the computationally expensive simulations. These simulations can then be interrogated as needed by suitable optimization algorithms [28, 29]. Recent advances in co-kriging formulations [30] have also enabled to construct regressions that can leverage training on combined datasets comprising results from both low -⁠ and high-accuracy simulation results (e.g., see [31] for an application in cardiac electrophysiology). It is therefore apparent that while inferring model outputs from a statistical surrogate is computationally efficient, the accuracy of the approximations strongly depends on the quality of the training set, which might require a large computational cost to be built.

In this work, we demonstrate a method to circumvent some of the limitations of the above strategies by adopting a different hybrid approach. We introduce novel low-order (LO) models of ventricular mechanics that match global results from state-of-the-art FE simulations across multiple patients. The LO models comprise two classes of model components: first, “trained” black box transformations with parameters sized by machine learning to match model outputs to features of the detailed models, and therefore only indirectly linked to the underlying biophysics; and second, “physical” components, with parameters shared with the high-resolution biophysical model and therefore accessible to inference and prediction of underlying biophyiscal states, previously only available with the detailed models.

Fig 1 shows a schematic of the method. A typical FE model is composed of multiple elements, all of which in turn incorporate models of the biological units. Coupled interactions between FEs and the imposed boundary conditions ultimately result in the global outputs of the detailed model. The corresponding LO description accounts for only one or few biological units that are coupled to the global outputs through a series of transformations. In the simplest case, a transformation can be a linear operator (e.g., a scaling matrix or linear regression) defined by the set of “trained” parameters. Training is then performed by ensuring congruency between outputs of the LO and FE models upon a set of perturbed conditions.

As an example application, we construct LO models of 2 representative ventricles extracted from the publicly available Sunnybrook Cardiac MRI database [32]. After showing that the LO models are indeed congruent to their corresponding FE simulations upon varying training conditions, we test them under different simulated scenarios showing remarkable accuracy. Additional simulations performed on multi-unit LO models accounting for smoothly varying heterogeneities provide additional insights on how the LO models, despite their simplifications, can capture the global behaviors. We then explore the relationship between our “trained” phenomenological parameters and the anatomical features of the ventricles by training LO models for 106 “virtual” anatomies sampled from realistic distributions of geometric features. Finally, we employ the LO models as surrogates for the high-resolution FE simulations in an inverse optimization problem aimed at matching cardiac phase durations and ejection fraction measured from cardiac MRI.

Materials and methods

3-D FE models of left ventricular mechanics

By integrating the cellular contraction model into a 3-D finite element formulation of soft tissue biomechanics, we sought to explore how the cooperative behavior of myocardial myocytes would translate at the macroscopic organ scale. FE models predict global outputs such as time-varying ventricular volume and pressure while also naturally accounting for the detailed anatomies of the heart and for the complex boundary conditions acting on it (e.g., the behavior of the downstream vasculature). In this work, multi-scale coupling between cellular and tissue models was implemented according to a stabilized mixed u-P formulation with ad hoc preconditioning and adapted for P1/P1 elements [27, 33]. Briefly put, coupling between cellular and elemental stresses and strains occurred implicitly at the Gauss point level. Passive contributions from extracellular matrix and from myocyte structural proteins were described as a hyperelastic behavior following the constitutive model by Usyk et al. [34] (see S1 Supplemental Material for more details). To simulate also active contraction, calcium transients triggering myocyte response had to be synchronized to the pacing dictated by heart electrophysiology. Rather than simulating in detail the generation and propagation of action potentials, we followed a simplified decoupled approach where information on arrival timing of the electric impulse wavefront was computed by solving the Eikonal equation with anisotropic conductivities [35]. Specifically, myocardial fibers were assumed to have 2 times higher conductivity than myocardial sheets and 4 times higher conductivity than mutually normal cross-fibers. Studies on mapping electrical activation of the heart suggest that Purkinjie fibers might resurface at a small number of endocardial foci to initiate impulse propagation [36]. We selected, therefore, 4 small groups of endocardial elements as zero-time regions in the eikonal solution. In polar coordinates (z¯,θ), the set of 4 prescribed activation foci was (z¯,θ)={(0.6,70°),(0.3,140°),(0.6,−150°),(0.05,50°)}, where the first value of each tuple indicates whether the activation focus was closer to the base (i.e., z¯=0) or to the apex (i.e., z¯=1.0), while the second value represented the circumferential coordinate (in degrees) measured counterclockwise around the longitudinal axis of the ventricle starting from the approximate center of the left ventricular free wall (θ = 0°).

Computational domains and boundary conditions

While the FE formulation could be in principle applied to arbitrarily complex geometric domains such as detailed multi-chamber models, this study targeted only the behavior of left ventricles (LVs). Specifically, we employed a 6-parameter axisymmetric representation of ventricular anatomy that was recently employed to automatically process the Sunnybrook Cardiac MRI database [27, 32]. The parameterization scheme captured ventricular anatomy as a truncated prolate spheroid defined by the parameters: R_b, the outer LV radius at base; L, the LV thickness at base; Z, the main longitudinal axis dimension including thickness at the LV apex; H, the LV thickness at the apex; e, the sphericity/conicity parameter, with values between 0.5 (rotation cone) and 1 (perfect sphere); Ψ₀, the truncation angle, which typically assumes negative values. Fig 2 shows axial and longitudinal cross-section profiles of the two LV geometries that were selected as representative of the normal (N) and heart failure (HF) subject groups included in the database. To at least partially correct for the fact that all imaged ventricle configurations are subjected to non-negligible loads (e.g., due to intraventricular pressure and external boundary conditions), we first “unloaded” the geometries reconstructed from MRI using Gaussian Process regression under the assumption of a 10% mid-wall end-diastolic strain for the N geometry and a 15% mid-wall strain for the HF one (see leftmost and central columns for cross-section profiles at end-diastole and after unloading, respectively) [27]. The 2 idealized anatomies were then discretized into 3-D meshes of 34 590 (N) and 49 121 (HF) linear tetrahedral elements (see righmost column). Ventricle bases (marked in yellow) were prevented to move out of plane, while a 5 mm-wide stripe of epicardial elements (marked in red) closest to base was fully constrained to prevent rigid motions. Pressure and volume couplings between the left ventricle, the left atrium and the distal vasculature was weakly enforced at the endocardial surface as a Neumann boundary condition (marked in blue).

Low-order models of ventricular mechanics

Parallel schematics of the high-resolution FE models and the here proposed corresponding LO counterparts are reported in Fig 3. In terms of formulation, the FE models relied on a PDE-based description of mechanical equilibrium coupled via hemodynamic boundary conditions to a simplified representation of atrial pressure (P_at) and to a 3-element Windkessel model of the downstream vasculature with additional resistance to split aortic valve R₁, aorta R₂, and peripheral artery resistance R₃. In training the LO model, we used a simpler configuration with R₂ = 0, combining aortic and peripheral system. In the model, the capacitance of the peripheral system is denoted as C. The LO models maintained identical hemodynamic components, but also simplified the description of the computational domain (Ω) to a set of algebraic relationships linking cellular active tension (T_a) to left ventricular pressure (P), and cellular elongation (λ) to left ventricular volume (V).

Passive behavior

The hyperelastic passive mechanical behavior of ventricles was modeled by the 3-D strain energy function presented in [34] (see S1 Supplemental Material for more details). To reproduce an equivalent behavior, our LO models were endowed with empirical relationships mimicking the combined effects of nonlinear hyperelasticity and ventricular geometry. Specifically, we considered separate formulations for the responses observed at volumes higher and lower than the unloaded volume V₀,

where α, β, γ are coefficients modulating the pressure-volume exponential behavior upon compression; a, b, c are coefficients of the 3^rd-order polynomial describing the ventricle under passive tension. The differential treatment of the behavior under compression or tension was necessary to ensure good fits of the LO model behavior to the FE simulations.

Active behavior

Coupling between active myocyte response and time-varying ventricular hemodynamics is a pivotal component of the LO models introduced here. The relationships driving cell action can be summarized as

where S is a vector containing the time-varying state variables of an ODE model of myofilament contraction that depends on calcium concentration, Ca, and on strain of the myofilament, λ; state variables and strain also uniquely determine the active tension generated by the myofilament T_a (see S1 Supplemental Material for the specific f₁ and f₂ used for this work).

The relationships in (2) need to be opportunely transformed to incorporate the role played by the coordinated contraction of the left ventricle, with correct sizing of the transformations being of great importance for building LO models for a given ventricular anatomy of interest. The relationships linking cellular strain to ventricular volume and active tension to ventricular pressure are given by

where μ₁ and μ₂ are opportune scaling coefficients of the active component of the low order model, and P_pass follows what shown in (1).

Congruency training

We now describe the strategy used to find the overall best-fit coefficients for the linear transformations in (3) from high-resolution FE runs. To ensure that the LO models maintained correspondence to their FE counterparts over a wide range of myocyte lengths and downstream hemodynamic resistances, we considered a set of 3 active and 1 passive training simulations encompassing varying hemodynamic conditions. More specifically, to extract the passive coefficients {α, β, a, b, c} we fitted the LO model predictions to a FE pressure-inflation test performed using an intraventricular pressure range of 0-4 kPa. Note the estimation α, and β was done under the assumption of neutral role of γ (i.e., γ = 1), as we defined γ to capture differential behavior upon compression conditions (see (1)), which are not tested upon passive inflation. Best-fit values of the remaining parameters {μ₁, μ₂, γ} were instead found from active simulations accounting for different combinations of atrial pressure and Windkessel parameters (see Table 2). In other words, for any ventricular geometry defined by {R_b, L, Z, H, Ψ₀, e}, we solved two sequential problems to maximize correspondence between global outcomes of the FE and the LO models. A first minimization problem leveraged the results of the passive inflation simulation,

where the subscripts _FE and _LO are used to label the left ventricular volumes predicted by the FE and LO models, and p is a set of left ventricular pressures considered in the passive simulation Ⓟ. The optimal best-fit coefficients for the passive material behavior were then used in a second minimization problem seeking the active “trained” parameters,

where the discrepancies between global outputs P and V are computed for all of the active training simulations (i.e., {①, ②, ③}, see Table 2). To solve the problem in (4) we used a standard curve fitting procedure based on the Levenberg-Marquardt algorithm available in SciPy. Solving (5) was less straightforward, as it required iterative runs of the LO model to compare global outputs during active contraction. To exploit parallel execution on multicore machines, we used a custom implementation of the OrthoMADS algorithm [37], with search steps interrogating an iteratively updated kriging surrogate to speed up convergence [28, 29].

Multi-element low-order models

In a simulated heartbeat, the many cell units of a high-resolution model undergo deformations that may vary based on their anatomical location (e.g., across the cardiac wall or along the apex-base direction). As a result, even for the same computational domain, cellular models of active contraction may operate at heterogeneous sarcomere lengths that cannot be explicitly accounted for a single-cell LO description. To incorporate heterogeneity and evaluate approximation errors, we then extended the LO models to account for multiple interconnected elements that share the same left atrium and Windkessel models. Fig 4 shows schematics of the two considered multi-element configurations, with LO submodules connected either in parallel (see panel A), or in series (see panel B). Outputs from each configuration were then compared to a single-element LO model that was adapted via congruency training to maximize correspondence between global outputs of the single-element and multi-element models. Congruency training was carried out solving (5) and the error analysis was repeated for 3 different atrial pressures, P_at = 0.3, 0.6, and 1 kPa.

Parallel configuration

The parallel model can be viewed as the combination of LO subcomponents that share the same ventricular volume V but operate at different sarcomere lengths because of non-uniform μ₁ values. In this case, Eq (3) takes the form

where λ and T_a are the stretch ratio and the active tension of the LO submodels, respectively, and ξ ∈ [0, 1] indicates the discretization of the system into subcomponents. To include heterogeneity, the μ₁ parameter was imposed to vary linearly within the [μ₁(ξ = 0), μ₁(ξ = 1)] = [0.01, 0.179] range. We overall considered 256 LO subcomponents.

Series configuration

The series model combines LO subcomponents each contributing a volume portion, V(ξ), to the global left ventricular volume, V, while operating under the same ventricular pressure, P. The weights w(ξ) introduce heterogeneity in the series multi-element model by modulating the effects of a subset of LO parameters, {α, a, b, c, μ₂}. Accordingly, Eq (3) takes the form

where ξ ∈ [0, 1] indicates the discretization of the system into subcomponents. The w(ξ) weights were imposed to vary linearly within the range [w(0), w(1)] = [0.6, 1.0]. Overall, the model was discretized into 256 LO subcomponents connected in series.

Regression for congruency coefficients

While the minimization problems (4) and (5) were solved for any given geometry on the basis of 4 high resolution runs, the computational burden would scale up quickly in situations where many anatomies might be of interest. For this reason, we built a Gaussian Process regression model mapping the geometric parameters of the left ventricle (i.e., {R_b, L, Z, H, Ψ₀, e}) to the 8 parameters of the LO module (i.e., {α, β, γ, a, b, c, μ₁, μ₂}). To train the regression model, we randomly sampled 150 virtual anatomies uniformly distributed over the range of geometric features observed for the patients included in the Sunnybrook Cardiac MRI datasets [27]. 3-D high-resolution simulations were run for all geometries and under all training conditions (see Table 2). To ensure meaningful fits also for the γ parameter, which acts only when LVV is smaller than V₀, we excluded from subsequent training the cases for which the left ventricle did not contract enough to reach 85% of the unloaded volume V₀ under hemodynamic conditions ②, which imposed a minimal afterload constraint. Efficacy of the Gaussian Process regression model in learning the behavior of the LO indices was evaluated via 5-fold cross-validation.

Inverse problem solution: Fit to MRI data

To test the efficacy of our LO models as convenient surrogates of high-resolution 3-D FE counterparts, we employed the overall model reduction framework to solve an inverse parameter identification problem. More specifically, we sought values for the model parameters that would ensure an optimal match with the ejection fraction and cardiac phase durations measured from the Sunnybrook Cardiac MRI scans [27, 32]. Knowing that heart failure might likely affect both calcium handling and overall contractility, we chose to optimize the τ₂ and S_a parameters, which modulate the duration of the calcium transient in the cellular unit and the overall efficiency of contraction, respectively. At the same time, we also penalized deviations of central aortic blood pressure from the normal range 85-110 mmHg by adapting the hemodynamic parameters of the Windkessel submodel, i.e., the aortic resistance R₂, the distal resistance R₃, and the capacitance C. We used the same value for the aortic valve resistance as during the congruency training simulation (i.e., R₁ = 0.0015 kPa ml⁻¹ s). Finally, we assumed that the left atrium would contribute to diastolic relaxation by providing a progressively decreasing blood pressure. For this, we employed an electrical circuit analogy and represented the atrium as a capacitor. The capacitance value, C_at, was considered unknown and optimized together with {τ₂, S_a, R₂, R₃, C}.

As in 5, the inverse optimization problem was solved via OrthoMADS, leveraging parallel execution. Comparison analyses were finally carried out with respect to the pressure and volume traces computed for the optimal set of parameters.

Results

In our approach, the simulated global behavior is reproduced at a reduced computational cost by transforming the outputs of a modeled cellular subcomponent. To determine a suitable functional form for the cell-to-global transformations, our first analyses investigated the relationships between local and global biomechanical measures in active FE simulations. Fig 5 shows LVV plotted against mid-wall stretch in the fiber direction (λ_mw) in simulations run under the same boundary conditions used for congruency training (i.e., ①, ② and ③ in Table 2). Probing locations were taken across the ventricular wall at half of the apex-to-base distance and averaged over 40% and 60% of thickness. As expected, simulated perturbed conditions led to pronounced changes in local stretch for both the N (see black lines) and the HF (see gray lines) cases, with stretch ranges within 0.88-1.21. A least-square error fit revealed that the relationship between LVV and λ_mw could be reasonably assumed to be linear for both considered geometries (R² = 0.99).

This relationship between simulation outputs at the global (LVV) and local (λ_mw) levels caused us to adopt simple linear forms for the transformations in (3). Fig 6 shows the outputs from high-resolution FE (solid lines) overlayed with outputs from LO models (dashed lines) after parameter optimization via congruency training. Best-fit values for the trained parameters {μ₁, μ₂, γ, α, β, a, b, c} are reported in Table 3 for the N and HF ventricular geometries.

The perturbed hemodynamic parameters considered during training affected significantly global performance both in terms of ejection fraction (observed ranges within 44%-60% and 40%-60% for the N and HF geometries, respectively), and peak systolic LVP (observed ranges within 77-158 mmHg and within 91-174 mmHg for the N and HF geometries, respectively). Note that LO and FE models employed the same values for all the “physical” parameters modulating myocyte contraction and hemodynamic coupling. Outputs from the two classes of models matched closely, with small discrepancies appreciable between pressure traces (see top row of Fig 6) and almost complete overlap occurring at the volume level (see central row in Fig 6). The largest, but still contained, mismatch was appreciable for simulation ③ in the HF case (see gray lines in the top right panel Fig 6). For both classes of models, simulated pressure-volume loops are shown in the bottom row overlayed to results from the passive inflation tests run to optimize the LO coefficients {α, β, a, b, c}. Both relations in (1) resulted in good matches during passive filling, but the polynomial form provided better accuracy at lower pressures.

Despite supporting our choice for the functional forms (3), the results of Fig 6 did not provide any indications on whether the LO models would show congruent behavior also in scenarios not directly considered during training. To assess performance under more general conditions, we then compared global outputs of the two classes of models when independently varying sets of parameters modulating cellular action as well as a hemodynamic scenario not accounted for during training. Meanwhile, the trained parameters of the LO models were maintained fixed to the optimal values obtained via congruency training without further tuning (see Table 3). Fig 7 compares LVP (first and third rows) and LVV traces (second and fourth rows) obtained in the additional simulated scenarios. The first column shows effects of varying the maximum velocity of muscle shortening, V_max, by simultaneously modulating the attachment and detachment rates of crossbridges, i.e., f_XB and g_XB, respectively. The 3 shown simulations account for 50% rate variations above and below the baseline values used for training. Altering crossbridge dynamics resulted either in faster (for increased f_XB and g_XB) or slower (for decreased f_XB and g_XB) contractions during systole, as appreciable from the different LVV trace slopes observed during ejection. The LO results (see dashed lines) followed closely their FE counterparts (see solid lines) for both the N (top two rows) and HF (bottom two rows) cases, and with mismatches limited to the isovolumic relaxation phase, when ventricles operated under compression behavior.

For a second scenario, we compared LO and FE simulation results for varying length-tension relationships. When progressively increasing n₀, i.e., the minimum stretch required for myocytes to generate active force, ventricles generated progressively less force and contracted less. Maximum relative mismatch (i.e., absolute pointwise difference normalized by pointwise mean value) in LVP and LVV traces between FE and LO was ∼ 3% on average. The third scenario is again depicted in Fig 7 (see third column) and targeted the n_A parameter, which defines coperativity in the interaction among fraction of troponin with bound calcium, tropomyosin state, and crossbridge dynamics (see S1 Supplemental Material for more details). Increasing n_A resulted in smaller forces of contraction (i.e., larger end-systolic volumes) combined to faster ejection and isovolumic relaxation phases. These findings can be attributed to crossbridges being recruited at lower calcium levels when n_A is low, while increasing n_A resulted in crossbridges more quickly reaching their maximum generated force. Even for this set of simulations, relative errors between LO and FE models were small, with a maximum relative error of ∼5% observed for the HF geometry at the lowest value of n_A. For the fourth scenario (see second column from right), we varied τ₂, which impacts the relaxation time of the intracellular calcium transients and, hence, contraction duration. Higher τ₂ values resulted in prolonged ejection and isovolumic relaxation times, and a maximum mismatch of 5% between LO and FE results was observed for the LVV trace of the HF ventricle at the longest considered τ₂ (i.e., for τ₂ = 285 ms). For the final scenario (see rightmost column), we progressively increased R₁, the aortic resistance, to simulate the large afterloads that ventricles might need to overcome, for example, due to stenosis of the aortic valve. We considered 3 sets of simulations where R₁ was progressively increased up to 0.1 kPa ml⁻¹ s, which is half of the distal arterial resistance. As expected for a 3-element Windkessel model, increasing the ratio between proximal and distal resistances resulted in smoother transition in ventricular pressure at the opening and closing of the valves. Again, as expected, maximum contraction was achieved for the lowest R₁ value (i.e., R₁ = 0.025 kPa ml⁻¹ s), while maximum ventricular pressure was reached for the largest R₁ (i.e., R₁ = 0.1 kPa ml⁻¹ s).

Overall, the good matches between FE and LO global outputs observed both under training conditions (see Fig 6) and for the simulated scenarios of Fig 7 suggested that, despite accounting for just one active cell unit, the LO models were able to effectively capture the global behavior of 3-D FE models, which instead results from the cooperation of many units operating under heterogeneous strain conditions. For further investigation, we then decided to extend the LO modeling framework to account for heterogeneity by allowing for the global ventricular dynamics to be captured by multiple LO sub-components connected either in parallel or in series and varying for their material properties (see section Multi-element low-order models). More specifically, for units connected in parallel, we imposed variations in μ₁ within the range [μ_1,min, μ_1,max] = [0.01, 0.179]. Similarly, heterogeneity was achieved for the series configuration cases by modulating 6 LO parameters, {a, b, c, α, β, μ₂}, with weights varying within the range [w_min, w_max] = [0.6, 1.0]. To analyze the capability of single-element models to approximate multi-element configurations, we then trained LO models with a single cell unit to maximize the match in global ventricular volumes and pressures. Average mismatches, quantified as C (see (5)), between single and multi-element models are reported in Table 4 upon variations of the same testing parameters considered for Fig 7. In all cases, the single element models were able to capture well the multi-unit behaviors, with a maximum mismatch of 3.5% observed for the parallel configuration at the lowest V_max considered. Maximum mismatch with the series configuration was observed instead for the largest τ₂ value (C = 2.59%).

As an example, Fig 8 shows pressure (top row) and volume (bottom row) traces from multi-element models (connected in series) overlayed on the corresponding traces from best-fit single-element models. The results shown derived from two different levels of myofilament cooperativity (i.e., n_A = 10 and and n_A = 5, in the left and right panel, respectively) demonstrating almost complete overlap between outputs from the two modeling schemes (i.e., compare dashed and dotted black lines). Heterogeneity within multi-element models was evident when comparing volume traces from the two most different subunits of the model (i.e., those modulated by weights at the range extremes). As weights modulated also passive stiffness, the w_max (w_min) subunit was least (most) compliant and therefore contributed the smallest (largest) volume subportion to the resulting total trace. As expected, no differences could be appreciated among pressure traces, as subunits connected in series were designed to sustain the same left ventricular pressure.

As the the above results show, after training, even single cell-to-global transformations (2) were able to capture the structural organization and heterogeneity accounted for in FE models and that we imposed in our multi-element LO models. To characterize better the links existing between ventricular anatomy and corresponding trained LO parameters, we then repeated congruency training for a set of virtual ventricles representative of the anatomies collected by the Sunnybrook Cardiac MRI database. Similar to [27], we employed an axisymmetric 6-parameter geometric description and used latin hypercube sampling to extract 150 ventricles from uniform distributions constructed over the entire range of the dataset anatomies. From this initial sampling, 44 simulations were discarded as either not compatible (e.g., due to unrealistically large ventricular thickness as compared to radius), or because they did not contract enough (i.e., LVV did not reach 85% of V₀ at peak contraction, see section Regression for congruency coefficients for more details). For the remaining 106 simulations, congruency training was achieved with mean square errors always lower than 1%. Fig 9 shows scatter plots and overall trends of the LO parameters after training from the active simulations (μ₁, μ₂ and γ) with respect to the variations of the two most influential anatomical parameters (R_b and the ventricular thickness L). When considering the trends separately (see panel A), all trained parameters were most affected by L. Some of this dependence on L can be expected, as thicker ventricles are able to accomodate more force-generating myocytes, thus explaining the direct proportionality between L and μ₂ (see central row). The opposite was true for μ₁, which transforms LVV into myocardial strain, indicating that, in thicker ventricles, changes in volumes tend to correspond to smaller changes in effective myocardial strain (see top row). Effects on γ (see bottom row) were more complex with a minimum γ = 0.92 observed for L = 13 mm and higher values observed both at smaller and larger thicknesses. The role played by R_b emerged more clearly in contour plots showing expected coupled effects when varying also L in a 2-D Gaussian Process regression (see panel B). In general, R_b had opposite effects with resepect to L, which resulted in the maximum expected parameter values being located in the vicinity of either the top (for μ₁) or bottom (for μ₂ and γ) left corners. In addition, the combined variations assumed locally complex behaviors (e.g., see isocontours in the lower half the μ₁ plot), which would difficult to capture without a statistical regression model. Cross-validation results of the regression are reported in the S1 Supplemental Material.

Fig 10 shows simulated global outputs after the model parameters were optimized to match cardiac phase durations and ejection fraction from MR imaging for the N (left panels) and HF (right panels) cases. In addition, we show volume measurements extracted from MRI that provided the target for the optimization (see solid circles in the two bottom panels). As expected, the MRI analysis revealed a reduced ejection fraction (-19%) and a prolonged ejection time (+21% after normalizing to a 60 bpm heart rate) for the HF case than to the N one (see Table 3 in the S1 Supplemental Material). Accordingly, the optimization algorithm achieved optimal matches by assigning smaller S_a (-46%) and larger τ₂ (+31%) values to the HF case compared to the N one. Given our simplified description of the atrium, we did not correctly reproduce the dynamics of diastoling filling. Accordingly, our target variables (i.e., ejection fraction, ejection time, isovolumic relaxation time, and aortic pulse pressure) did not include diastolic filling metrics, as they would not be helpful without proper accounting of atrial and mitral valve function. Average mismatches to the target variables were smaller than 3% in both cases, suggesting that the combined tuning of hemodynamic parameters, calcium handling via τ₂ and contraction efficiency via S_a might be sufficient to constrain the model to global contraction measurements from images.

Discussion

Models of cardiac mechanics promise to provide significant insights into the state of health of patients with heart disease by providing quantitative descriptions of the biophysics driving cardiac contraction at different space and time scales. Image-based mechanistic reconstructions such as state-of-the-art FE models allow for capturing high levels of detail in anatomic complexity and spatial heterogeneity (e.g., to account for the varying preferential directions of myofibers), but are also computationally expensive. Not surprisingly, significant obstacles are added when large parameter searches and inverse problem solutions are required to adapt these models to experimental or clinical observations. To at least partially address such challenges, we have demonstrated a procedure to train efficient LO models of ventricular mechanics and reproduce the global behavior of FE models by scaling the outputs of just a single modeled cell. The computational costs of cardiac models were therefore radically reduced. By implementing the LO model to exploit GPU computations via CUDA, we were able to run 48 simulations in parallel (20 heartbeats each) in 5 seconds on a workstation mounting an nVidia Quadro K5200. In comparison, 20 heartbeats in a FE model would require, depending on resolution, between 1 and 40 CPU hours on high performance servers.

Our strategy was motivated by our finding that the relationship between local fiber strain and ventricular volume could be well linearized over a large range of stretches (see Fig 5), and surprisingly also when accounting for material and geometric nonlinearities [33, 34], and spatially-varying distribution of fiber orientation and heterogeneous activation times [35]. By lumping the nonlinearities introduced by the 3-D architecture of the myocardial tissue at the passive level (e.g., during filling) into an additive pressure term (see (1)), we were able to reproduce the global behaviors observed at the organ level by just linearly scaling the active contributions of the cellular model (see (3)). In effect, our approach exploited the links between length dependence active force and Frank-Starling mechanism [38] to generate a concise model of global heart function, which can be computed with minimal cost. Coupling between local cell behavior and global outputs was achieved by transformation modules incorporating the linear relationships (3), which were trained to ensure congruency with FE results upon varying hemodynamic conditions. Given the empirical nature of the transformation relationships, the “trained” parameters {μ₁, μ₂, γ, α, β, a, b, c} did provide only an indirect physiological meaning. Although we did not conduct a systematic investigation to minimize the number of high-resolution simulations necessary for training, our results suggested that even just 3 active and one passive high-resolution simulations were sufficient to achieve good fits of the LO model to the FE model (see Fig 6).

We submit that the presented model order reduction strategy presents significant advantages over standard machine learning approaches where training is conducted directly on high-resolution simulations (e.g., [27, 31]). While efficient at the inference stage, statistical models tend to require extensive training to achieve sufficient accuracy. By maintaining a biophyisical component (i.e., the cell unit), our LO models can be more easily trained, especially given that probing simulations can be chosen to explore large variations in operating sarcomere lengths (see Fig 5 and Table 2). After training, LO models then achieve remarkable computational performance, with the possibility of simulating thousands of heartbeats every second on commodity hardware. This level of computational efficiency is comparable to alternative order reduction methods [13–16], but, unlike these, our approach does not require restrictive assumptions on ventricular anatomy. The popular CircAdapt model [13], for example, achieves computational efficiency by assuming sphericity of ventricles, and then by analytically exploiting symmetry. The congruency training proposed here can handle any anatomy that the FE model can handle. The LO model also outperform the speed-ups typically reported for more sophisticated numerical approaches, such as reduced-basis functions and proper orthogonal decomposition (e.g., [17, 18, 21, 22]), although, in the current form, our reduction technique does not maintain details on spatially-varying fields.

To investigate further the properties of the LO reduction approach, we also considered extensions of the single-element LO architectures, and built configurations where multiple LO subcomponents are connected either in series or in parallel (see Fig 4). Even when introducing heterogeneity among material properties of different subcomponents, the global outputs could still be recapitulated by a trained single-cell analogue. This finding was also backed analytically under the assumptions of a simplified model of active contraction and of linear passive material properties (see section “Approximation of multi-element model with single-element LO model” in the S1 Supplemental Material). Our empirical and analytical analyses demonstrate that the global aspects of the biomechanics simulated by FE models could indeed be captured by these simplified descriptions.

Given their computational efficiency and demonstrated versatility, the presented LO models are thus ideal for solving difficult inverse parameter estimations problems. To showcase this strategy, we adapted model parameters to capture contraction features measurable from cardiac imaging. More specifically, we relied on our previous analysis of the Sunnybrook Cardiac MRI database [27] to select 2 representative patients from the normal and heart failure groups, and found model parameterizations that would match ejection fraction and cardiac phase durations from MRI, while also maintaining central aortic pressure within a reasonable 85-110 mmHg range. As expected for this type of complex inverse estimation problem, parameter optimization required running the LO models for several thousands of parameter combinations and for tens of beats at each evaluation to reach steady-state. Tackling such an inverse problem would be impractical with a standard FE technique, even when considering efficient surrogate-based optimization algorithms [28, 29, 37], as we did in this study. By exploiting multi-thread execution on GPU, we were instead able to optimize LO model parameters in just a few minutes on commodity hardware, and achieved good matches (see Fig 10). Interestingly, the optimization algorithm did not reveal large differences in hemodynamic lumped parameters, but it attributed instead the discrepancies in ejection fraction and cardiac phase duration to the S_a and τ₂ parameters, which modulate cardiac contraction efficiency and decay time of calcium transients, respectively. More specifically, our analysis suggested that the HF MRI measurements could be explained by smaller S_a values and longer τ₂ times, which is in agreement with experimental findings on heart failure myocyte phenotypes, which are less effective in contracting partially due to dysfunctional calcium handling [39, 40].

Given the encouraging preliminary results in extracting mechanistic insights from clinical assessments of cardiac function, we submit that, after validation, our approach could prove helpful to assist statistical analyses of large databases of clinical measurements and provide interpretations for variable associations. By providing an efficient model-order reduction strategy that does not compromise anatomical details, our work might also facilitate future studies involving high resolution simulations, where FE and LO models could be employed in combination for effective inverse estimation of parameter values.

Supporting information

S1 Supplemental Material [pdf]
Myofilament contraction model, optimized model parameters, cross-validation of GP regression and approximation of multi-element LO model with single-element LO model.