The schematic of our physics-based model in the reference state is shown in the left panel of Figure 2. The right panel shows an example of deformation with end-diastolic ventricular wall thinning. Conversely, the wall thickens during systolic contraction. Matching the simplicity of the single curve for either EDPVR or ESPVR, the cardiac shape is approximated by a thick-walled sphere. The use of such simplified geometries dates back to the early days of cardiac modeling and can be found still in modern research (Anani and Rahimi, 2016; Kalhöfer-Köchling, 2020). The reference geometry Ω₀ of the LV is shown on the left-hand side, whereas the right-hand side depicts the geometry in a deformed state. The inner and outer radius of the sphere for the reference geometry are R _endo and R _epi , respectively. The wall thickness is thus R _epi − R _endo . While in the deformed geometry, the inner and outer radius become r _endo and r _epi , correspondingly.

A spherical coordinate system is adopted so that its origin is located at the center of the sphere. The three orthogonal basis vectors of coordinate systems are (e _r ,e _θ ,e _φ ). A deformation maps a point X in the reference geometry Ω₀ to point x in the deformed geometry Ω. For a given point X in the reference geometry, the radial coordinate is R, while the corresponding radial coordinate of the point x in the deformed geometry is r. The second and third coordinates of the point, θ and φ, keep unchanged under deformation due to the assumption of centrosymmetry, which will be clarified later. The deformation gradient tensor is defined as F = ∂ x ∂ X . (3) Under the spherical coordinate system, it is straightforward to get F = diag λ ρ , λ θ , λ φ , (4) where λ _ρ is the radial strain; λ _θ and λ _φ are two tangential strains.

We assume that the myocardium is incompressible, resulting in the volumetric strain J = det(F) = 1, yielding the relation λ ρ λ θ λ φ = 1 . (5) The sphere only has expansion and contraction deformations, which means that the points in the domain only have radial displacement. We have the symmetry constraint λ θ = λ φ . (6) The strain λ _θ can be calculated by the ratio of the perimeter l of the cross-section through the spherical center (which is represented by the red dotted circle in Figure 2) on the deformed geometry and L on the reference geometry λ _θ = l/L, yielding λ θ = r R . (7) Substituting Eqs 6, 7 into Eq. 5, we can get the radial strain λ ρ = R 2 r 2 . (8)

Numerically, the right Cauchy-Green strain tensor C is a better choice for solving the balance equation than the deformation gradient tensor F, since the former is symmetric and positive definite for all points X ∈ Ω₀, which reduces computational costs. Said right Cauchy-Green strain tensor is defined as C = F T F . (9) Substituting Eqs 4–6 into Eq. 9, the right Cauchy-Green deformation tensor reads C = diag λ ρ 2 , 1 λ ρ , 1 λ ρ . (10) The first invariant of the right Cauchy-Green deformation tensor I ₁ can be expressed as I 1 = λ ρ 2 + 2 λ ρ . (11)

Due to the incompressibility of the myocardium, the volume between the inner surface and the red dotted spherical surface (see Figure 2 for reference), stays constant. It follows straightforward that R 3 − R endo 3 = r 3 − r endo 3 . (12) The geometrical parameters are further normalized by the inner radius R _endo and the ventricular volume V ₀ of the reference geometry as follows R ̂ = R R endo , r ̂ = r R endo , δ = R − R endo R endo , Δ = R epi − R endo R endo , V ̂ = V V 0 , (13) where Δ is the normalized thickness; V is the ventricular volume at deformed geometry. With Eqs 12, 13, the normalized radius can be expressed as r ̂ = R ̂ 3 + V ̂ − 1 1 / 3 . (14)

The total elastic energy W stored in the myocardium can be calculated by integrating the energy density function Ψ over the domain Ω₀ W = ∫ Ω 0 Ψ d Ω . (15) The sphere experiences an inner pressure P, representing the blood pressure inside the LV. While the outer pressure is set to zero, indicating a free boundary condition. Based on classical mechanics, it is known that any mechanical work W performed on the sphere due to a given internal pressure P follows the relation P = d W d V . (16) Substituting Eqs 9–15 into Eq. 16, we get an important relationship P = − 2 ∫ 0 Δ λ ρ 2 r ̂ d Ψ d λ ρ d δ . (17) Therein, the ventricular pressure is expressed as the integral of the energy density function Ψ over the domain defined by the normalized thickness. The EDPVR and ESPVR can be further deducted based on this relationship. Besides such mechanical work approach, another approach based on stress analysis can also be found in Kalhöfer-Köchling (2020).

The myocardium is considered to be a homogeneous, incompressible, anisotropic, and fiber-reinforced soft material that generates active forces. Inspired by experimental data, several constitutive laws (energy density function Ψ) for the myocardium have been developed in the last decades, including the orthotropic Holzapfel-Ogden model (Holzapfel and Ogden, 2020) and our recently developed squared generalized structure-tensor (SGST) models (Kalhöfer-Köchling et al., 2020), in which the fiber dispersion of the myocardium is taken into account. Considering the microstructure of the myocardium, these constitutive laws contain different contributions of isotropic, fibrous, and laminar structures as well as their coupling, and are used for simulations at tissue or organ level (Baillargeon et al., 2014; Chabiniok et al., 2016; Kalhöfer-Köchling, 2020).

As a first approximation, by neglecting the anisotropy of the cardiac tissue, the isotropic energy function is considered in this work Ψ = a 2 b e b I 1 − 3 − 1 , (18) where a and b are mechanical parameters representing the material property, i.e., the stiffness, which can be obtained from experiments.

Incorporating the constitutive law Eq. 18 into Eq. 17, one gets the passive contribution of the myocardium on the ventricular pressure, implying the EDPVR P E D = 2 a ∫ 0 Δ 1 − λ ρ 3 r ̂ e b I 1 − 3 d δ . (19) Therein, the input parameters are a, b and Δ, while the output is a function indicating the relation between the ventricular volume and the ventricular pressure at the ED state.

There are two contributions to the stresses and strains in cardiac muscle tissue: the passive and the active components. As for the passive contribution, the tissue generates resistive stress when it is deformed. This tension contributes to the ventricular pressure of the LV. On the other hand, the tissue actively contracts and the active force generated inside the tissue also contributes to the ventricular pressure. Active contraction of the cardiac muscle can be modelled in two different ways: using the additive stress approach or the multiplicative strain approach, as described in Pezzuto (2013). We have adopted the former due to its simplicity and intuitive nature. Thus, we have the total pressure P = P _p + P _a and the energy function Ψ = Ψ _p + Ψ _a . The indices p and a represent the passive and active parts, respectively.

To model the ESPVR, an active contribution, which indicates the active force generated by the myocardium during the ES state, is added onto the passive contribution Eq. 18. Such an active contribution reads Ψ a = T a λ 2 2 − λ 0 λ , (20) where T _a is the maximum active stress; λ ₀ = l ₀/l _r with the sarcomere smallest length l ₀ = 1.58 μm and rest length l _r = 1.85 μm (Guccione et al., 2001). The strain of the sarcomere λ is defined as λ = C : H a . (21) The active force structure tensor H _a is defined such that contractile forces act in the tangential plane of the myocardium H a = I − e r ⊗ e r , (22) Substituting Eqs 21, 22 into Eq. 20, the energy function for active force reads Ψ a = T a 1 λ ρ − λ 0 2 λ ρ . (23) Incorporating Eq. 23 into Eq. 17, we obtain the active contribution of pressure P a = 2 T a ∫ 0 Δ 1 − λ 0 λ ρ 2 r ̂ d δ . (24) The pressure at the ES state contains both the passive and active parts. Adding Eq. 24 onto Eq. 19, we get the ESPVR P E S = 2 ∫ 0 Δ a 1 − λ ρ 3 e b I 1 − 3 + T a 1 − λ 0 λ ρ / 2 r ̂ d δ . (25)

To validate the physics-based EDPVR model, we fit different models to the dataset from Klotz et al. (2006) and Klotz et al. (2007) and compare them in Figure 3. The dataset contains ex vivo EDPVR data for 80 human hearts. The fitting was implemented by minimizing the mean squared error (MSE) for the models. The full dataset with pressure up to 30 mmHg, and a subset with a physiologically reasonable pressure range (up to 20 mmHg), were considered respectively.

As we can see in Figures 3A, B, the green solid lines representing the physics-based model show a good fit to the experimental data. We compared our model with two other widely used EDPVR models, i.e., the exponential model (Artrip et al., 2001; Rodriguez et al., 2015; Shimizu et al., 2018) and the Klotz curve (Klotz et al., 2007). The Klotz curve (Eq. 2) belongs to the family of polynomial power functions, while the former (Eq. 1) is classified in the form of exponential functions. Our physics-based model entails the combination of an exponential energy function (Eq. 18) with a volume integral, hence resulting in exponential behaviour. It should be noted that during the curve fitting the original exponential model was adapted to the same normalized form as the Klotz curve. This was done by replacing the term V − V ₀ with V _n , which yields P E D = A ( e B V n − 1 ) .

The optimized parameters for the three models are given in Table 1. Contrary to previous models, the parameters in our model have a physical meaning. For example, a and b together reflect the stiffness of the material. The fitted values of a and b of our model are close to most values from experiments and simulations in the literature (Baillargeon et al., 2014; Nikou et al., 2016), although parameter estimates themselves often vary considerably across different datasets and experimental protocols. These two parameters are exactly the same as those in the isotropic constitutive law (Eq. 18). Δ is the normalized thickness of the LV wall, which is an important measure of the ventricular geometry.

The resulting curves of the three models in the full dataset are shown in Figure 3A. The MSEs of the original exponential model, the Klotz curve, and the physics-based model are 2.82, 3.19, and 3.08, respectively. Since our physics-based model uses an exponential constitutive law, it is fundamentally similar to the original exponential model. This is why the curves of the physics-based model and the original exponential model are almost identical, and both perform better than the Klotz curve. Figure 3C presents the residuals of the exponential model and the Klotz curve, using the physics-based model as a reference. It can be seen that the exponential model is much closer to the physics-based model, especially in the region with small volumes. It is also worth mentioning that the Klotz curve is not consistent with its definition at the maximum volume or pressure. When the normalized volume V _n is equal to 1, the resulting pressure in the Klotz curve is the same as the value of A _n , which, by design, may be different from the expected 30 mmHg.

The physics-based model shows its strong utility for small volumes. The normalized volume V _n of this region ranges from 0 to 0.8, corresponding to the pressure 0–20 mmHg, which covers the EDP of the human heart. To better evaluate these three models within this physiologically reasonable range, we generated a sub-dataset with pressures no more than 20 mmHg. The results are shown in Figures 3B, D. Here, the physics-based model shows the best fit with a MSE of 2.46. The MSE for the exponential model is 2.48. The Klotz curve has the worst fit with a MSE of 2.69, indicating its weakness at small ventricular volumes.

The Klotz curve is often used as a reference when estimating material parameters of myocardium, like stiffness, in numerical simulations (Krishnamurthy et al., 2013b; Hadjicharalambous et al., 2017; Sack et al., 2018). These simulations mostly use exponential constitutive laws to describe the mechanical properties of the myocardium. Based on the above-mentioned comparison, the new physics-based model shows the capacity to replace the Klotz curve in similar simulations in the future.

Because of its bottom-up, physical nature, the model can be used to predict the EDPVR of a ventricle with given information, such as mechanical properties and thickness of the myocardium. In Figure 4A, Dokos et al. (2002) represents the parameters for pig hearts, Demiray (1972) and Marx et al. (2022) are for the human hearts. We further performed finite element simulations using COMSOL Multiphysics (version 6.0.0.405) with the same parameter sets. The geometry of the finite element model is the same as that of the physical model, i.e., a thick-walled sphere. The normalized inner radius and thickness are 1.0 and 0.27, respectively. The inner surface was subjected to a constant pressure mimicking the ventricular pressure from blood, while the outer surface was free. Due to rotational symmetry, only 1/8 of the sphere with 26,392 tetrahedral elements (119,614 degrees of freedom) was simulated using the symmetry boundary condition. Equation 18 was used for the passive myocardium. The volumetric energy function Ψ _v = κ((J ² − 1)/2 − ln J) was applied to ensure tissue incompressability. The bulk modulus κ was 1 GPa; J = det(F) was the volumetric strain. Finally, in Figure 4B the normalized volume V/V ₀ and pressure obtained from the simulation are compared with the algebraic solutions of our physics-based EDPVR model. The perfect match between the finite element simulation and the theoretical model is not surprising as they are conceptually identical and only the calculation scheme differs. This is a demonstration of the relation between constitutive models commonly employed in cardiac mechanics simulations and our algebraic equation.

By changing the thickness Δ in our model, we studied how myocardium thickness affects the EDPVR of the LV, as shown in Figure 5A. An increased myocardium thickness leads to an upward lift of the EDPVR curve. To keep the same volume, higher pressure is needed when increasing the thickness of the LV, as shown in Figure 5B. This reflects the strong utility of a physics-based model over that of simply fit, i.e., the physics-based model is predictive over a wide range of parameters while a fit cannot.

The ESPVR describes the relationship between ventricular pressure and volume at the ES state of the LV as loading conditions change. It is composed of two parts, as shown in Eq. 25. The first part is identical to the EDPVR and the second part comes from the active force generated by the myocardium. The two contributions are shown in Figure 6A. Therein, the solid blue line is the overall ESPVR, while the dashed green one is from the active contraction. A large proportion of the pressure in the LV during the contraction is due to the active force generated by the myocardium. The shape of the ESPVR, especially when the normalized ESV is less than 1.0, is mainly determined by the passive resistance due to the deformation of the myocardium. In addition, the ESPVR with small volume is roughly linear, while the overall curve is almost bilinear. This implies that both the linear or bilinear forms of the ESPVR used in the literature have some validity. Please also notice that for the passive part of the ESPVR, the pressure is shown as a negative value. Negative pressure means that the tissue is resisting the deformation caused by its own active contraction. The greater the deformation at ES state, the more negative pressure (resistance stress) is generated, so that the required ventricular pressure is lower. It should be noted that the resistance stress increases in a nonlinear fashion with the decrease of the ESV, due to the nonlinear stress/strain behavior of the tissue. This results in the nonlinear shape of the ESPVR. Our model also shows that the positive slope of the ESPVR is mainly due to the decrease of the resistance stress as ESV increases.

To validate the proposed ESPVR model, we performed an additional finite element simulation, with the same geometry, mesh, boundary condition and volumetric energy function as previously considered for Figure 4B. The active force generated by the myocardium was determined by the energy function Eq. 20. For the passive response of the myocardium, the energy function was chosen according to Eq. 18. Parameters used in both the physics-based ESPVR model and numerical simulation are a = 1.15 kPa, b = 3.82, and T _a = 76.90 kPa. In Figure 6B, the orange dots present simulation results while the blue solid line is from the physics-based ESPVR model. In a large region of the normalized volume, the ESPVR predicted by the physics-based model agrees with the simulation results very well.

To further check the validity of our physics-based ESPVR model, we compared it with experimental ESPVR data. The experimental data for a pig heart was extracted from Moriz et al. (2021), by changing afterload pressure. The stress-free volume of the LV, V ₀, was assumed to be 50 mL based on Figure 3B of Moriz et al. (2021). Since the mechanical property and contractility of that pig heart were unknown, a = 2.52 kPa and b = 6.79 (Dokos 2002 in Figure 4) were adopted from Dokos et al. (2002), which performed shear test for cardiac tissue of the same specie; T _a = 81.89 kPa was obtained by fitting our ESPVR model to the experimental points and this value is physiologically reasonable. Other parameters used in our physics-based ESPVR model are Δ = 0.27 and λ ₀ = 0.85 as usual. As shown in Figure 7, our model shows good agreement with the experimental data with a maximum relative error of 6.77%.