Assessing left ventricular pump function using an immersed boundary method combined with finite elements

Cardiovascular disease remains the leading cause of mortality and morbidity worldwide. Given the vital role of cardiac pump function, its precise assessment and prediction under altered conditions are essential for diagnosing and managing various heart diseases. Tremendous effort has been devoted to developing computational models of the heart over the last several decades,¹ ranging from lumped parameter models² to 3D four-chamber heart models.^3–5 Cardiac dynamics are inherently multiscale and multiphysics in nature, involving nonlinear soft tissue mechanics, ventricular blood flow, myocardial perfusion, electrophysiology, and their interactions. It is now widely recognized that computational cardiology, or in silico modeling, when integrated with clinical imaging, can provide profound insights into various heart conditions and ultimately improve patient management.^6,7 As a result, personalized computational modeling of cardiac function is becoming increasingly used in clinical practice, leading to a digital twinning era of precision cardiology.^1,8 Various cardiac models have been developed, including structure-only models with subject-specific geometry, anisotropic hyperelastic myocardium, and active contraction,^9,10 as well as electromechanical coupling,¹¹ valve dynamics,¹² and whole-heart models.^5,13 Moreover, computational models have been used to study a range of cardiac diseases, including heart failure,¹⁴ myocardial infarction,⁶ and mitral regurgitation,¹⁵ among others.

The primary function of the heart is to pump sufficient blood to the body through the contraction of the myocardium. One model for understanding cardiac function is the cardio-mechanics model. Mechanically, the myocardium is highly nonlinear and anisotropic, with complex microstructures,^16,17 making accurate modeling of myocardial behavior a significant challenge.^18,19 Two types of myocardial constitutive models have been described in the literature: phenomenological models describe the stress–strain relationship on the basis of observed macroscopic behavior without delving into the tissue’s microstructure,²⁰ whereas structural constitutive models are based on the tissue’s microstructural organization and incorporate components such as collagen and elastin fibers to provide a more detailed and accurate representation of myocardial mechanics at the microscopic level.²¹ Although structural constitutive models offer greater accuracy, they are more computationally expensive than phenomenological models, and consequently the latter remain widely used for personalized cardiac modeling. One very popular constitutive model is the H–O model,²⁰ which has been successfully applied to in vivo patient data in various studies.^22–24 

Predicting systolic stress inside the left ventricle (LV) further requires the incorporation of active contraction,^25,26 which is coordinated by the heart’s electrical activation through intracellular Ca²⁺ dynamics.^27,28 This process is known as excitation–contraction coupling.²⁹ Electrophysiological models, such as monodomain or bidomain models, have been coupled with ventricular mechanical models.³⁰ These coupled models have been used to study how different pacing sites affect cardiac pump function and to investigate the regulation of sarcomere length-dependent tension development on myocardial contractility.³¹ 

Consideration of fluid–structure interaction (FSI) is essential when assessing the heart’s pump function, as blood is ejected through the heart by myocardial contraction. However, it is still very challenging to simulate blood–myocardium interaction, owing to the large deformation of the heart wall. One widely used method is the arbitrary Lagrangian–Eulerian (ALE) approach,^32,33 which has the advantages of being able to handle large deformation, incorporate advanced constitutive laws, and make use of well-developed numerical solvers. The ALE approach has been widely used in modeling cardiovascular systems, such as heart valves,^34,35 cardiac dynamics,^36,37 and LV assist devices,³⁸ to name but a few. New computational methods have also been developed within the ALE framework through the introduction of stabilization methods. In particular, the variational multiscale method allows the use of computationally efficient elements, such as equal-order simplicial elements.³⁹ Nonetheless, the ALE approach, paired with the variational multiscale method, will continue to enjoy wide application, owing to its effectiveness and robustness in handling complex FSI problems.⁴⁰ With the ALE approach, however, the interface between the fluid domain and the solid domain needs to be tracked, and the computational meshes deform continuously to account for the moving interface, and thus large structure deformation might lead to severe mesh distortion. In such cases, frequent re-remeshing is needed, and this can become a problem.

Alternative approaches for modeling FSI have also been developed, including the immersed boundary (IB) method pioneered by Peskin in the 1970s,⁴¹ fictitious domain methods,⁴² and fully Eulerian methods.^43–45 By recasting the solid momentum equations into rate form, fully Eulerian methods solve the whole coupled system in Eulerian space, thereby avoiding the need for mesh deformation, while the fluid–solid interface needs to be tracked using the level-set method or phase-field methods. Fictitious domain methods treat the structure as part of the fluid domain, with its motion being constrained through penalty methods or Lagrangian multipliers, and thus the fluid–solid interface is usually diffused. In the IB method, the continuity, momentum, and viscosity of the coupled fluid–structure system are described using Eulerian quantities, while the solid deformations, stresses, and forces are described using Lagrangian quantities, and the fluid–solid interface is explicitly defined from the solid domain. Other FSI approaches also exist, such as the immersed particle method,⁴⁶ and smoothed particle hydrodynamics.⁴⁷ It is not our intention, however, to review the various FSI approaches in this study, and interested readers are referred to Refs. 48 and 49 for recent reviews of this aspect. The choice of a specific FSI approach will depend on its suitability for the research problem at hand. In this study, we will use the IB method because various studies have demonstrated its effectiveness and accuracy in modeling cardiovascular systems.^50,51 Conventionally, the immersed solid structures are modeled using pseudo-elastic “fibers”,^41,52 which makes it difficult to incorporates modern constitutive laws for soft tissue. To address this limitation, further developments of the IB method have been introduced,^53,54 allowing the use of experimentally informed hyperelastic constitutive models by incorporating finite-element-based structural models. The one that we use here is a hybrid finite difference/finite element IB method, referred to as IB/FE.⁵⁴ Using the IB/FE approach, we have successfully simulated cardiac dynamics with fiber-reinforced myocardium⁵⁵ and myocardial infarction,⁶ further coupled to valves,¹² the whole heart,⁵ and finally to blood perfusion within the myocardium.⁵⁶ 

In this study, we utilize an IB/FE human LV model^7,12 derived from in vivo cardiac magnetic resonance (CMR) imaging of a healthy human volunteer. The model includes geometrically idealized “inflow” and “outflow” tracts for the assessment of LV pump function. We inversely determine both the passive and active parameters of the myocardium by matching the predicted LV dynamics to the image-derived measurements. Using this IB/FE LV model, we investigate the effects of diastolic aortic pressure and intracellular Ca²⁺ transients on LV pump function.

An existing CMR study of a young healthy volunteer (28 years, male) was selected for constructing the LV geometry. The CMR study was approved by the local ethics committee, and written informed consent was obtained before the CMR scan. LV function was assessed using steady-state free precession cine images, including a stack of short-axis cine images of the LV from base to apex and a long-axis cine image obtained in the LV outflow tract (LVOT) view, as shown in Fig. 1. Typical cine imaging parameters were 10 mm for slice thickness, 1.3 × 1.3 mm² for in-plane pixel size, and 25 frames per cardiac cycle. Twelve planes along the LVOT view were further acquired to image the mitral and aortic valves. Typical parameters were 3 mm for slice thickness, 0.7 × 0.7 mm² for in-plane size, and 25 frames per cardiac cycle. The flow rate across the aortic valve was further measured using phase contrast imaging.

Given that the LV pressure is lowest at early-diastole, the early-diastolic phase is chosen to reconstruct the unloaded LV geometry, similar to Ref. 10. Seven short-axis slices from base to apex are used to reconstruct the LV geometry, and the inflow (left atrium) and outflow (aorta) tracts are reconstructed using long-axis slices along the LVOT view; see Fig. 1. The ventricular wall boundaries are manually extracted using an in-house MATLAB (The MathWorks Inc., Natick, MA, USA) code; see https://github.com/HaoGao/GlasgowHeart. Following the manual segmentation, we import LV wall boundaries into SolidWorks (SolidWorks Corp., Waltham, MA, USA) for generating the three-dimensional geometry using B-spline surface fitting, as shown in Fig. 2(a). The LV geometry is then meshed with 232 000 tetrahedral elements, as shown in Fig. 2(b). We further assume that the valvular region (between the two valve rings and the base) and the inflow and outflow tracts can only bear the load passively without any contractility. Thus, only the region below the base, denoted as the LV region, can actively contract. A rule-based method⁶ is used for generating myofiber architecture, including the fiber f, the sheet s and the sheet-normal n: the so-called f–s–n system. From the endocardium to the epicardium, the myofiber angle changes from −60° to 60°, and the sheet angle changes from −45° to 45°. Figure 2(c) shows the generated myofiber directions f in the LV region.

The whole computational domain Ω is set to be 16.5 × 16.5 × 16.5 cm³, as shown in Fig. 3, and the LV geometry is fully immersed inside this fluid box, with the inflow and outflow tracts attached to the top plane. We fix the inflow and outflow tracts during the simulations by using a spring-like tethering force. We also constrain the longitudinal and circumferential displacements of the LV base, but allow the radial expansion so that the LV region can contract freely in systole. Because the valvular region cannot contract, a spring-like tethering force is further applied to the valvular region in systole to avoid excessive expansion due to the high blood pressure, similar to what was done in Refs. 7 and 12.

We now introduce the fluid boundary conditions on the inflow and outflow tracts of this coupled LV model according to the four different phases of a cardiac cycle: (a) diastolic filling, (b) isovolumetric contraction, (c) systolic ejection, and (d) isovolumetric relaxation. In detail,

• Diastolic filling. A linearly increased pressure from 0 to a prescribed end-diastolic pressure (EDP) is applied to the endocardial surface of the LV region over 800 ms. We set EDP to 8 mmHg, a population-based value, because ventricular pressure measurement is not available for the healthy volunteer. We also set the pressure in the inflow tract P_in = 0, allowing free flow across the mitral valve, and the velocity in the outflow tract u_out = 0 (equivalently the flow rate Q_out = 0), because the aortic valve is closed in diastole.

• Isovolumetric contraction. We set both the inflow and outflow rates to be zero, Q_in = Q_out = 0, because both the mitral valve and the aortic valve are closed during isovolumetric contraction. The duration of the isovolumetric contraction will be determined by the LV contraction when the LV cavity pressure is higher than the aortic pressure.

• Systolic ejection. At the end of isovolumetric contraction, the LV pressure exceeds the aortic pressure, and then the aortic valve opens. The aortic pressure is initially set to be 85 mmHg, the measured cuff diastolic pressure of the healthy volunteer. When the aortic valve opens, a simplified systemic circulation model represented by a three-element Windkessel model (see Fig. 3) is connected to the outflow boundary to provide physiologically accurate pressure boundaries on the outflow tract. The stored pressure P^Wk(t) in the Windkessel model is initially set to be the same as the diastolic aortic pressure (85 mmHg), and thus the flow across R_p is zero when the aortic valve just opens. Subsequently, P_out(t) is determined by P^Wk(t) and the flow rate Q_out from the LV model. When Q_out reaches zero, closure of the aortic valve occurs; we then set Q_out = 0 and disconnect the systemic circulation model from the LV model.

• Isovolumetric relaxation. Because both the aortic valve and the mitral valve are closed in this phase, we set Q_in = Q_out = 0 again. The isovolumetric relaxation finishes when the total duration for one cardiac cycle reaches 1200 ms. During the whole cardiac cycle, zero pressure is applied to the remaining boundary of the fluid box ∂Ω\(∂Ω_in ∪ ∂Ω_out).

A time step size of Δt = 0.122 ms is chosen in the diastolic and relaxation phases. The large deformations in systole usually require much smaller time steps; thus, 0.25Δt is used in the first half of systole, and then 0.125Δt in the second half of systole. The delta function kernel in Eq. (4) is approximated by a standard four-point regularized version of the delta function. To achieve a density of approximately two quadrature points per Cartesian mesh width when approximating the integral transforms in Eq. (4), dynamically generated Gaussian quadrature rules are used. Further details of the spatial discretizations and time stepping scheme of this IB/FE framework can be found in Ref. 54.

The IB/FE LV model is implemented with the open-source IBAMR software https://ibamr.github.io/). IBAMR is a distributed-memory parallel implementation of the IB method with adaptive mesh refinement,^61,62 leveraging functionality provided by other open-source software, including SAMRAI⁶³ for Cartesian grid discretization management, PETSc⁶⁴ for linear solver infrastructure, and libMesh⁶⁵ for finite element discretization management. All simulations are run on a local Linux workstation using ten CPU cores (Intel Xeon CPU E5-2699, 2.30 GHz with 128 GB memory). The computed LV pressure, aortic flow rates, and active tension obtained in the second simulated cycle are within 1% compared with the third and fourth periods. Therefore, we use results from the second period in the following analysis.

In total, 12 simulations have been performed under different conditions. They are as follows:

• Baseline case. The passive and active parameters are determined by matching the diastolic and systolic LV dynamics measured from the CMR study of the healthy subject.

• Effect of diastolic aortic pressure. High blood pressure is considered to be a major risk factor for myocardial infarction and stroke. Favorable clinical outcome has been demonstrated in hypertensive patients by lowering systolic and diastolic arterial pressures.⁶⁶ Here, we run the IB/FE LV model with four additional diastolic aortic pressures: 70, 80, 90 and 100 mmHg.

• Effect of Ca²⁺ profiles. In a failing myocyte, mishandling of Ca²⁺ is a central cause of both contractile dysfunction and arrhythmias.²⁹ To maintain sufficient cardiac output, the failing myocytes undergo a remodeling process that will result a prolonged action potential and altered Ca²⁺ handling. Changes in Ca²⁺ handling include among others, reduced Ca²⁺ level, increased Ca²⁺ transient duration, and prolonged Ca²⁺ decay time.⁶⁷ To mimic such remodeling in intracellular Ca²⁺ transients, we consider seven additional intracellular Ca²⁺ transients in two groups as shown in Fig. 4: for Ca²⁺ profiles in group 1 [Fig. 4(a)], the peak Ca²⁺ increases from 0.7 to 1.1 μM; for Ca²⁺ profiles in group 2, they differ mainly in the time-to-peak compared with the baseline. An extreme condition, case 7, has both a lower peak and a much longer decay time. The decay time is defined as the duration in which the Ca²⁺ level decreases to 0.1 μM after reaching the peak. The peak level, time-to-peak, and decay time of different Ca²⁺ profiles are summarized in Table I.

The optimized baseline parameters for the passive and active responses of the healthy subject are a = 0.19 kPa, b = 5.08, a_f = 1.2 kPa, b_f = 4.15, a_s = 0.70 kPa, b_s = 1.6, a_fs = 0.24 kPa, b_fs = 1.3, and T_scale = 4. Using these optimized passive and active parameters, the simulated end-diastolic LV cavity volume is found to be 142.6 ml, which is 0.4 ml less than the measured value (<1%), and the ejection fraction is 52%, which agrees well with the measured value of 57%. Figure 5 further shows the deformed LV geometries at end-diastole and end-systole superimposed with CMR images in short-axis and long-axis views. Overall, the deformed LV geometry agrees well with the CMR measurements of the subject, although it can be noticed that the predicted apex moves slightly beyond the CMR-imaged apex as shown in Fig. 5(c) (indicated by the blue arrow), which could be explained by the fixed basal plane along the long-axis direction.

Figure 6 shows vector plots of the ventricular blood flow in diastole (top panels) and systole (bottom panels) superimposed on the corresponding LVOT CMR images. In diastole, vortices form inside the LV cavity and moving toward the apex to mix with the blood from previous pump cycles and push it farther up toward the outflow tract in late-diastole. In systole, it can be seen that the blood is ejected out through the outflow tract with a strong jet, and the velocity in the outflow tract is much higher than the velocity inside the LV cavity.

Figure 7 compares the simulated aortic flow rate during systole with the CMR measurements. The maximum flow rate from the IB/FE LV model is 491 mL/s, slightly less than the measured peak value 498 mL/s. The simulated peak value from the IB/FE LV model also arrives earlier compared with the measured data, which is possibly caused by the lack of a realistic aortic valve in our model. On the basis of the simulated aortic flow rate, we find that the systolic ejection starts 51 ms after the myocardium begins to contract, and ends at 307 ms with an ejection duration of 256 ms. The ejection duration is slightly shorter by about 30 ms compared with the in vivo measurement, which could be explained by the mismatch between the individual intracellular Ca²⁺ transient and the prescribed profile in Eq. (12). The peak LV cavity pressure in systole is 161 mmHg, which is comparable to the measured systolic cuff pressure of the healthy volunteer (150 mmHg). In Fig. 8, we further compare the circumferential strains in mid-ventricle with the values derived from the CMR cine image at six segments in line with the standard classification of the American Heart Association:⁶⁸ the mean difference is −0.008 ± 0.03. These close agreements with the CMR measurements suggest that the baseline model is functioning physiologically with the applied boundary conditions and estimated parameters.

Figure 9(a) shows the active tension at peak systole when the LV cavity pressure is highest. It is found that the active tension is almost uniformly distributed within the LV wall, except for the apical region. Figure 9(b) plots the mean active tension and the myofiber stress throughout the LV region during systole. In early-systole, the myofiber stress is slightly higher than the active tension, suggesting that the myofibers are being stretched, whereas in late-systole, the active tension is higher than the total myofiber stress, suggesting that myofibers are being compressed.

As the diastolic aortic pressure rises, the average myocardial active tension and LV cavity pressure increase, as shown in Fig. 10. It can be seen that the aortic flow rate is slightly lower, whereas the active tension is increased. With decreased diastolic aortic pressure, both the mean active tension and the LV cavity pressure decrease [Figs. 10(a) and 10(c)]. On the contrary, the aortic flow rate increases slightly [Fig. 10(b)]. Consequently, decreasing the diastolic aortic pressure will increase the LV cardiac output with a smaller end-systolic LV volume: see the pressure–volume loops in Fig. 10(c). From Table II, it can be seen that decreasing the diastolic aortic pressure will reduce the peak active tension, the peak LV cavity pressure, and the stroke work (the area enclosed by the pressure–volume loop), but with increased stroke volume. This implies that the LV works more efficiently when the diastolic aortic pressure is lower.

Figure 11(a) shows the LV active tension in systole triggered by the Ca²⁺ transients from group 1 [Fig. 4(a)]. The LV active tension drops significantly in systole for cases 1 and 2, but increases in case 3 with a higher peak Ca²⁺. Similar results can be found for the aortic flow rates shown in Fig. 11(b). Case 1 has a much lower aortic flow rate with the shortest ejection duration (216 ms), which is 40 ms less than the baseline case. The pressure–volume loops in Fig. 11(c) suggests that for case 1, the LV pumping capability may be inadequate, as reflected by its stroke volume of 44 mL and a much smaller stroke work (65 cJ) that is half of the baseline case (139 cJ). Both the stroke volume and stroke work increase in case 3 with a higher peak Ca²⁺. Figures 11(d)–11(f) show the effects of the different Ca²⁺ transients from group 2 on the LV pump function. For cases 4, 5, and 6, the active tension and the aortic flow rates are similar in early-systole because the three Ca²⁺ transients have same peak values, and are slightly different only in the time-to-peak. Noticeable differences can be found in late-systole, which may be explained by the different decay time. The longer decay time will maintain a relatively higher level of Ca²⁺, which will result in a relatively higher active contraction and pressure, and a longer systolic ejection duration. In case 6, the ejection duration is increased by 20 ms compared with the baseline case. In case 7, because of the lower Ca²⁺ peak, both the peak active tension and the peak aortic flow rates are much smaller than the baseline case. However, since the Ca²⁺ decay time is much longer, the LV pump function seems to be preserved; see Fig. 11(f). As a result, the stroke volume is only 4 ml less than the baseline case, suggesting that the long decay time could compensate for the lower Ca²⁺ peak to maintain the cardiac output in case 7. A summary of LV pump function with different Ca²⁺ transients from Fig. 4 can be found in Table IV.

Table III further summarizes the changes in the LV pump functions resulting from variations in diastolic aortic pressure and the intracellular Ca²⁺ transients. Decreasing the diastolic aortic pressure not only increases the stroke volume, but also lowers the LV active tension and LV pressure. Changes in the Ca²⁺ transients have more complex effects on the LV pump function. High peak Ca²⁺ will generally decrease the duration of isovolumetric contraction and increase the LV stroke volume. The time-to-peak Ca²⁺ correlates positively with isovolumetric contraction duration. A longer decay time of the Ca²⁺ profile will help to maintain the stroke volume, even with low levels of intracellular Ca²⁺.

We have developed an IB/FE human LV model and used it to study LV pump function under physiological and pathological conditions, including variations in diastolic aortic pressure and intracellular Ca²⁺ levels. The extended inflow and outflow tracts enable simulations of more realistic blood–myocardium interactions. The LV geometry is derived from a clinical CMR study of a healthy subject, with passive and active parameters determined by matching the measured cardiac function, including circumferential strains and LV cavity volume. The simulated LV dynamics agree well with in vivo LV dynamics (Figs. 5 and 7), even with a simplified electromechanical coupling for active contraction.

This LV model has allowed us to quantitatively investigate how different factors, ranging from the cellular level to the tissue level, affect LV pump function. We have found that increasing diastolic arterial pressure leads to increases in peak active tension and stroke work, in agreement with findings by Veress et al.⁶⁹ Clinical outcomes have demonstrated that lowering systolic and diastolic arterial pressure for patients with heart conditions is very beneficial.⁷⁰ In this study, we have shown that when the diastolic aortic pressure is decreased (from 85 to 70 mmHg), both the peak active tension and LV cavity pressure are slightly decreased, but the stroke volume is increased; see Table II and Fig. 10. On the other hand, when the after-load increases, the LV needs to contract further to generate a higher pressure. The aortic flow rate is then reduced, and systolic ejection ends up with a higher end-systolic pressure. These changes result in a reduced stroke volume. For example, the case with the highest aortic diastolic pressure experienced a 5 ms longer duration of isovolumetric contraction, 7 ms shorter ejection duration, and 5 kPa higher peak active tension than the values in the baseline case. This seems consistent with the “Anrep effect” that occurs when the aortic blood pressure increases to maintain appropriate stroke volume.⁷¹ 

In the failing heart, mishandling of Ca²⁺ is a central cause of both contractile dysfunction and arrhythmias.²⁹ To maintain sufficient cardiac output, the failing myocytes undergo a remodeling process that results in a prolonged action potential and altered Ca²⁺ handling.⁷² By mimicking changes in Ca²⁺ transients under different pathological conditions, such as the effect of inotropic drugs and mishandling of Ca²⁺ in failing myocytes, we are able to study the effects of different Ca²⁺ profiles on LV pump function. In case 3 from Fig. 4(a), the peak Ca²⁺ is 10% higher than the baseline case, and the peak active tension is increased by 5 kPa (5%), the peak LV cavity pressure is increased by 7 mmHg (4%), and the stroke work is increased by 15 cJ (11%), with slightly increased stoke volume (2 ml). If the Ca²⁺ level is substantially reduced [see case 1 from Fig. 4(a)], the stroke volume decreases to 44 ml with an ejection fraction of 31% and a much lower LV cavity pressure than the baseline case. Owing to the dramatic drop in stroke volume when a low-level Ca²⁺ transient is present, negative inotropic drugs may precipitate or exacerbate heart failure.⁷³ Considering case 7 in Fig. 4(b), the low Ca²⁺ level leads to a low level of active tension, and consequently a low LV cavity pressure, but because of the increased Ca²⁺ transient duration and prolonged decay time, the ejection duration is 30 ms longer than the baseline case. Therefore, the final stroke volume is reduced by only 4 mL. However, if the failing myocyte is unable to remodel the Ca²⁺ transient with a longer decay time, the stroke volume will drop significantly.

Although our IB/FE LV model provides some interesting results, its limitations should be mentioned. First of all, estimating material parameters of the myocardium from in vivo clinical data remains a major challenge in subject-specific cardiac modeling. In this study, we have used only the end-diastolic volume and the circumferential strains to infer the passive myocardial parameters without uncertainty quantification, which is similar to Refs. 6 and 22. We find that this approach could determine the mechanical response that matches the in vivo LV dynamics, but the uniqueness of inferred material parameters cannot be ensured. This is partially due to the correlations among these parameters, and the limited in vivo measurements that we can obtain. To overcome the challenge posed by these limited measurements, studies have used 3D tagged CMR for estimating myocardial motion and strain,^10,74 although 3D tagged CMR requires additional scan time and much more complex post-processing. Furthermore, crucial ventricular pressure data can only be measured invasively. Noninvasive approaches for estimating ventricular filling pressure using medical imaging techniques are under investigation⁷⁵ and could be useful in reducing the uncertainty imposed by the assumed end-diastolic pressure. Xi et al.⁷⁶ reported that an average error of 1.3 mmHg in the pressure measurements resulted in a maximum error of 11% in the estimated myocardial stiffness. They therefore suggested that noninvasive assessment of pressure using medical imaging data could be used for estimating myocardial stiffness, although further validation work is needed. Recently, Lazarus et al.²⁴ performed a comprehensive sensitive study on the H–O model making use of Gaussian surrogate models, including forward uncertainty quantification using global sensitivity and inverse uncertainty quantification in a Bayesian manner. They found that a and a_f could be reliably inferred from end-diastolic LV cavity volume and strains, and that b and b_f could be inferred if high-pressure data were available, such as the so-called Klotz curve.⁷⁷ By taking advantage of machine-learning-based approaches, an appropriate statistical inference pipeline could be integrated into personalized cardiac models to account for the uncertainties arising from both the data and the model.

Second, the active tension model⁵⁸ used in this study allows us to compute the active force from Ca²⁺ transient profiles, rather than relying on constant Ca²⁺ levels or action potentials.²⁶ One advantage of using the Ca²⁺ transient in the active tension calculation is that the model developed here can be readily linked with electrophysiological models. However, the parameters associated with the active tension model, including the Ca²⁺ transient, cannot be easily measured or inferred from in vivo patient data. Therefore, we prescribe the Ca²⁺ transient, rescale the active tension generated by myocytes through T_scale, and take other parameters from Ref. 58. This is because the original active tension model was developed using rat data. Our results show that by changing the peak and decay time of Ca²⁺, the LV pump function changes significantly, as shown in Table IV, and thus the rescaling procedure seems to be necessary. In the simulations, T_scale is determined by matching the model predictions with the measured ejection fraction and the aortic flow rate. However, it is challenging to match all the measurements by tuning only one parameter in the active contraction model, as can be seen in Fig. 7, where the simulated aortic flow rate peaks earlier than the measured value.

A further limitation of our work is that the structure above the two valvular rings, as shown in Fig. 2, is not physiologically derived. Although this is needed to implement the flow boundary conditions, the flow patterns near the two trunks will be different from in vivo conditions, where heart valves are present. Despite this, the flow patterns within the LV away from the boundary largely agree with physiological observations. We also assume that the valvular region can only passively bear the load, which also affects modeling accuracy, since there could be myocytes between the fibrous valvular rings and the basal region, and a transition from contractile myocardium to stiff fibrous tissue needs to be considered in the future studies. However, CMR cine images are unable to differentiate these. Therefore, simplified displacement boundary conditions are applied to the LV base, and the valvular region is described with a soft isotropic material model without contraction and kept in place through a spring-like tethering force during systole.

Finally, we have not considered myocardial growth and remodeling. Cardiac tissue can adapt its structure and function to maintain its normal pump function in response to external and internal environmental changes,⁷⁸ such as increased afterload or impaired contractility as we have studied here. There are two main frameworks in the literature that have been developed for modeling soft tissue growth and remodeling: kinematic growth theory⁷⁹ and constrained mixture theory.⁸⁰ Both approaches have been successfully applied to cardiac growth and remodeling under physiological and pathological conditions,^81,82 and other applications can also be found in the literature, such as to membranes⁸³ and skin.⁸⁴ 

We have presented an IB/FE cardiac model that simulates blood–myocardium interaction in the human left ventricle. Model parameters are inversely inferred to match measured LV dynamics in vivo. Using this LV model, we have studied the LV pump function under various physiological and pathological conditions. Our findings indicate that increasing diastolic aortic pressure will reduce LV pump function, while the intracellular Ca²⁺ transient plays an essential role in regulating LV pump performance. Specifically, a reduced Ca²⁺ transient may result in insufficient blood supply to the body, whereas a prolonged decay time could compensate for a low level of Ca²⁺, thereby maintaining stroke volume. This personalized modeling approach will deepen our understanding of cardiac pump functions across multiple scales and hold promise for improving risk stratification and treatment planning for various heart diseases.

H.G. acknowledges funding provided by the British Heart Foundation (No. PG/22/10930) and by the UK EPSRC (No. EP/S030875/1). H.G. and P.V. are also grateful for funding provided by the British Council the Going Global Global Partnerships - Springboard Programme (No. 1270441886). D.G. is grateful for funding provided by the National Natural Science Foundation of China (No. 12402353) and by Shandong Provincial Natural Science Foundation (Nos. 2024HWYQ-008 and ZR2024QA054).

The authors have no conflicts to disclose.

Hao Gao: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Funding acquisition (lead); Investigation (lead); Methodology (lead); Project administration (lead); Software (lead); Supervision (lead); Validation (equal); Visualization (equal); Writing – original draft (lead); Writing – review & editing (equal). Debao Guan: Conceptualization (equal); Data curation (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Validation (equal); Writing – review & editing (equal). Pierre-Frédéric Villard: Conceptualization (equal); Investigation (equal); Methodology (equal); Resources (equal); Validation (equal); Writing – review & editing (equal).

The simulation was implemented using the IBAMR software, freely available on GitHub (http://ibamr.github.io), Ref. 62. The left ventricular model will be made available on request.