Robust and efficient fixed-point algorithm for the inverse elastostatic problem to identify myocardial passive material parameters and the unloaded reference configuration

Image-based computational models of the heart represent a powerful tool to shed new light on the mechanisms underlying physiological and pathological conditions in cardiac function and to improve diagnosis and therapy planning. However, in order to enable the clinical translation of such models, it is crucial to develop personalized models that are able to reproduce the physiological reality of a given patient. There have been numerous contributions in experimental and computational biomechanics to characterize the passive behavior of the myocardium. However, most of these studies suffer from severe limitations and are not applicable to high-resolution geometries. In this work, we present a novel methodology to perform an automated identification of in vivo properties of passive cardiac biomechanics. The highly-efficient algorithm fits material parameters against the shape of a patient-specific approximation of the end-diastolic pressure-volume relation (EDPVR). Simultaneously, an unloaded reference configuration is generated, where a novel line search strategy to improve convergence and robustness is implemented. Only clinical image data or previously generated meshes at one time point during diastole and one measured data point of the EDPVR are required as an input. The proposed method can be straightforwardly coupled to existing finite element (FE) software packages and is applicable to different constitutive laws and FE formulations. Sensitivity analysis demonstrates that the algorithm is robust with respect to initial input parameters.

Image-based computational models of the heart represent a powerful tool to shed new light on the mechanisms underlying physiological and pathological conditions in cardiac function and to improve diagnosis and therapy planning. However, in order to enable the clinical translation of such models, it is crucial to develop personalized models that are able to reproduce the physiological reality of a given patient. There have been numerous contributions in experimental and computational biomechanics to characterize the passive behavior of the myocardium. However, most of these studies suffer from severe limitations and are not applicable to high-resolution geometries. In this work, we present a novel methodology to perform an automated identification of in vivo properties of passive cardiac biomechanics. The highly-efficient algorithm fits material parameters against the shape of a patient-specific approximation of the end-diastolic pressure-volume relation (EDPVR). Simultaneously, an unloaded reference configuration is generated, where a novel line search strategy to improve convergence and robustness is implemented. Only clinical image data or previously generated meshes at one time point during diastole and one measured data point of the EDPVR are required as an input. The proposed method can be straightforwardly coupled to existing finite element (FE) software packages and is applicable to different constitutive laws and FE formulations. Sensitivity analysis demonstrates that the algorithm is robust with respect to initial input parameters.

Introduction
Computational models of the heart are widely recognized as a powerful tool for the quantitative analysis of cardiac function. Their ability to explore mechanistic relationships of all variables of interest at high spatio-temporal resolution has turned computational models into an important, if not indispensable, adjunct in any basic research study. More recently though, driven by advances in medical imaging and simulation technologies, a translational trend has emerged that is geared towards turning cardiac modeling from a research tool into a clinical modality for diagnosis, stratification, and therapy planning [46,76,60].
However, unlike in basic research applications where a representative "one heart fits all" model is suitable to investigate generic mechanisms, in clinical scenarios the use of patient-specific models is of fundamental importance. Computational models must be built to represent the physiological reality of a given patient to provide a basis for therapeutic decisions for the specific patient, not for a representative average patient.
This process of minimizing the difference between the predictions of a computational model and clinical observations, referred to as "model personalization" or "digital twinning", is challenging to achieve [59]. Cardiac function emerges as a complex interplay between different physics -electrophysiology, mechanics, and hemodynamics. While accurate mechanistic computational models of each of these physical components exist and anatomically accurate multiphysics models of the heart, e.g., [6,50,80], are considered to be the state-of-the-art, they comprise a larger number of parameters, most of which cannot be observed in vivo or not at all. Thus, model parameters have to be inferred indirectly from clinical observable quantities which are, in general, afflicted with a significant degree of observational and residual uncertainties [52,51].
Integral model components in modeling cardiac function are constitutive relations which characterize the passive biomechanical behavior of myocardial tissue. The development and parameterization of these models is an active area of research, see [10] and references therein. As microstructural material components along local material axes determine the orthotropic behavior of the myocardium, recent experimental research has focused on multiaxial mechanical testing and microscopical investigations of the underlying structure [17,77,81]. These experiments are usually performed under artificial ex vivo conditions using excised tissue specimen from specific locations within the myocardium. While these experiments are an indispensable source for informing modeling, their inherent underlying limitations and uncertainties should be kept in mind when incorporating such data into models. The tissue is usually excised during autopsies within 12h of death and then stored in a suitable solution until testing. The testing should ideally take place immediately, but, in practice, may be postponed for weeks, although it has been shown for collageneous tissues that changes in microstructure and mechanics occur in as little as 48h after removal from the body [36]. Available material parameters are usually fitted to a normalized stress-strain curve obtained from multiple patients and used to represent the overall passive mechanical behavior of the ensemble. It is known though that tissue properties vary to a significant extent between individuals and throughout the myocardium. Additionally, even if orthotropic, microstructurally-based models are utilized, the underlying orthotropic fiber and sheet arrangement in the myocardium cannot be determined with certainty for individuals. To achieve translation it is therefore key to develop a methodology that facilitates the efficient, robust and, ideally, automated identification of model parameters from clinical observations in vivo [6].
A method towards in vivo passive parameter estimation was proposed by Augenstein et al. [5]. They utilized an ex vivo experimental method based on the combination of tagged magnetic resonance imaging (MRI) and simultaneous pressure recordings to estimate passive material properties in arrested pig hearts. The fiber architecture was incorporated from diffusion tensor MRI and the findings were validated against a silicon gel phantom with known material properties, proving that cardiac MRI can be used to extract meaningful material properties. Further work on isolated hearts combined with FE analysis was then conducted by e.g. Nair et al. [56], where an ex vivo rabbit left ventricular (LV) model was used for strain matching.
To take a step towards the translation of passive material identification into the clinics, in vivo MRI combined with ex vivo diffusion tensor MRI was carried out in heart failure patients [85] and, using a sequential quadratic programming optimization technique [12], passive material parameters were identified [86]. These studies used either transversely isotropic laws such as the Guccione law [32] or orthotropic laws such as the Holzapfel-Ogden (HO) model [39]. Further, initial parameter sets fitted to ex vivo experimental data (e.g. [33]) are taken and either only the isotropic parameter is varied [86,4,21] or two scaling parameters [58,46,72,64] are introduced to preserve the overall orthotropy of the material parameters. Additionally, material parameters of these laws are often correlated [70], and hence the unique identification poses a non-trivial problem.
Early studies attempted to identify material parameters based on non-invasive imaging of strain via MRI [57,92]. However, in vivo MRI strain measurements have major caveats as they vary among methods, modalities, and software version and mostly lack proper validation [2,69]. In addition, MRI strains are not usually acquired during routine clinical exams, and hence pressure-volume (PV) relations are more commonly used to estimate material parameters. There are two common fitting targets used in literature: fitting to displacement curves obtained from several frames during diastole or approximation of the empirical EDPVR by a power law as proposed by Klotz et al. [45].
Inverse estimation of patient-specific parameters is now widely used, where parameters of a forward problem are tuned to match PV-curves and motion fields available from clinical imaging [24], a method which was first used by Ghista et al. [28]. Xi et al. [89,90] utilized a reduced-order unscented Kalman filter [54] and later 3D tagged MRI to simulate unloading and estimate parameters in a reformulated Guccione law to overcome the problem of non-unique constitutive parameters. They used an early diastolic frame as the reference configuration and 29 frames per cycle for the fitting. Asner et al. [4] jointly estimated the reference state and passive parameters using deflation, considering a parameter sweep consisting of 25 simulations to assess the parameter space and find the minimum of the objective function. Gao et al. [24] used the first frame of early diastole as the unloaded reference configuration and assumed a population-based end-diastolic pressure, p ed , of 8 mmHg as no pressure recordings were available. They proposed a multi-step non-linear least-squared optimization procedure for the inverse estimation of all parameters of the HO law. The performance of multiple constitutive laws were compared by Hadjicharalambous et al. [34], who used 3D tagged MRI and a parameter sweep. The initial parame-ter set was chosen to match the Klotz relation. They found that the HO law provided the best balance between identifiability and model performance. Recently, Nasopoulou et al. [58] used 2D tagged MRI to improve the identifiability of parameters of the Guccione law by combining deformation and energy analysis to uniquely constrain the parameters. In more detail, the exponential scaling factor α was estimated through the minimization of an energy-based cost function. Then, a mechanical simulation was performed in a second step to optimize the linear isotropic parameter with a cost function based on differences in displacements. For this, the first frame and the two last frames of diastole were used to be compared to simulations of passive inflation. Finsberg et al. [20,21] used the backward displacement method to find the unloaded configuration and fitted one parameter of a reduced HO law with an initial guess of other parameters similar to Asner et al. [4].
The empirical EDPVR as proposed by Klotz et al. was utilized as a fitting target by, e.g. Nordsletten et al. [63], to estimate the unloaded geometry of human hearts. Simulations of passive inflation were conducted and compared to the Klotz relation, using least-square constraints on volume and pressure. In 2013, Krishnamurthy et al. [46] estimated the unloaded bi-ventricular geometry using a backward displacement method. They scaled parameters of a reduced HO law to fit the empirical Klotz curve by adjusting them manually and varying exponential parameters within a range of 15% of their initial values taken from experiments performed on canine tissue [17]. In 2018, Palit et al. [66] provided a summary of parameters of the HO law available in literature and examined the influence of p ed , fiber orientations and geometry on the estimated parameters. They used an objective function which minimized the differences between the simulated and real LV cavity volume and used the Klotz curve as a fitting target albeit they had no pressure recordings and assumed p ed to be 10 mmHg. The initial parameters were estimated by a Latin hypercube sampling (LHS) which generated 50 initial data sets. Sack et al. [72] used a similar approach, minimizing the error with respect to both pressure and volume. Initial parameters were fitted to data from [77] and scaled with one scaling parameter for the linear terms and one for the exponential terms to match the Klotz relation. For the passive filling calibration, they defined an objective function as the difference in pressure values along the PV-curve combined with a single measure of end diastolic volume, V ed . Genet et al. [26] calibrated two parameters of the Guccione law [32] using two nested loops: first, to optimize a parameter b, that defines the nonlinearity of the stress-strain relationship, such that the resulting loading curve is close to the Klotz curve and second, the scaling parameter of the Guccione law was optimized for a given b such that computed and prescribed end-diastolic pressures match.
In this paper, we present a novel model function-based fitting method (MFF) to find personalized material parameters for passive mechanical modeling of the heart. Only image data at one time point to create the anatomical model and one measured data point of the EDPVR are required as an input. The approach simultaneously personalizes passive material parameters and generates an unloaded reference configuration which is crucial for the image-based modeling of biomechanical LV diastolic function. In this regard, we also present an improved backward displacement algorithm where robustness is enhanced by a novel line search strategy. The MFF method was tested on a cohort of 19 LV meshes [49]. Excellent agreement with the empirical Klotz EDPVR was obtained for all cases under study. As the determination of an unloaded reference configuration and the identification of parameters is carried out simultaneously, the method is highly efficient, requiring only a small number of forward simulations (≤ 10 for all cases). Hence, computational costs were only about 2 to 3 times the cost of a standard passive inflation experiment. The versatile workflow is applicable to a large variety of constitutive models, model functions, and FE formulations and can be coupled to established FE software packages with relative ease. A thorough sensitivity analysis demonstrates the robustness of the method with regard to input uncertainty.

Patient data
Clinical data from the CARDIOPROOF study (NCT02591940), a recent clinical trial, was available, which includes data of N AS = 12 aortic valve stenosis (AS) patients and N CoA = 7 aortic coarctation (CoA) patients. As indicator for treatment, valve area and/or systolic pressure drop across the valve was taken into account for AS cases, whereas for CoA cases treatment indicators included an echocardiographically measured peak systolic pressure gradient across the stenotic region greater than 20 mmHg (2.66 kPa) and/or arterial hypertension. The institutional Research Ethics Committee approved the study following the ethical guidelines of the 1975 Declaration of Helsinki. From the participants' guardians written informed consent was attained. A detailed description of the data acquisition process and clinical protocols used in this study were reported in [19]. Pressure measurement in the LV, namely invasive catheterization, was routinely acquired in patients suffering from CoA only, making measured values of p dat ed solely available for this group. For AS cases p dat ed was determined empirically by statistically analyzing a reference data pool of N = 290 patient cases treated for AS. For more details on the reference data pool see the supplementary material of [49].

Unloaded reference configuration
To find an unloaded reference configuration Sellier [75] proposed an iterative, fixed-point method for general elastostatic problems, see Fig. 1 for a schematic and Algorithm 3 in Appendix A. In each iteration k a forward problem φ(X k ) is solved, where a measured in vivo pressure is applied to an interim reference configuration with coordinates X k . This results in updated configurations with coordinates x k and, subsequently, the per node displacement vector between the , where a measured in vivo pressure p ed is applied to an interim reference configuration with coordinates X 0 = x dat . This generates an updated configuration with coordinates x 0 . Subsequently, the per node displacement vector between the updated deformed configuration and the target in vivo configuration ( R 0 ) is computed. Finally, the reference configuration is updated by subtracting this per node displacement vector (blue dashed line). This iterative procedure is repeated until a given error tolerance between the computed reference configuration and the given in vivo configuration is reached. Here, the algorithm converges to the unloaded reference configuration X * after four iterations. (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.) updated deformed configuration and the target in vivo configuration R k is computed. Finally, the reference configuration is updated by subtracting this per node displacement vector. This iterative procedure is repeated until a given error tolerance between the computed reference configuration and the given in vivo configuration is reached and hence the unloaded reference configuration X * is found. Rausch et al. [68] augmented this approach, see Algorithm 4 in Appendix A, based on Aitken's delta-squared process [1] that increases the convergence rate and also significantly accelerates the method. Both fixed-point approaches are simple and versatile and can be coupled to existing FE software packages with relative ease and have thus been applied to a number of biomechanical modeling problems, e.g., [41,43,87].
However, for some of our LV models the augmented approach still diverged and unfavorable updates of the reference vector even resulted in a failure of the algorithm. Since our parameter fitting approach described later in Section 2.3.2 is based on a robust unloading strategy for soft tissues, we further improved this iterative method by a line search strategy based on Armijo [3], see Algorithm 1. Here, unfavorable search directions are damped and with this improvement the unloading algorithm was robust for all cases and in most cases the procedure was sped up. See Table A.10 for a comparison of the robustness and computational times for the three specified unloading algorithms.
We successfully applied Algorithm 1 to all 19 LV models of the cohort from Marx et al. [49], as well as meshes of bi-ventricular, see Karabelas et al. [42], and four-chamber models from the cohort presented by Strocchi et al. [78].

Constitutive models
The myocardium is considered as a non-linear, hyperelastic, nearly incompressible, and orthotropic material with a layered organization of myocytes and collagen fibers that is characterized by a right-handed orthonormal set of basis vectors; see [10] for a great review on material properties and constitutive modeling of the passive mechanical behavior of myocardial tissue. The basis vectors consist of the fiber axis f 0 , which coincides with the orientation of the myocytes, the sheet axis s 0 and the sheet-normal axis n 0 . To enforce the condition of a nearly incompressible material, the strain energy function , is split into a volumetric U ( J ) and an isochoric part isc (C ): In this relation, C = F F is the right Cauchy-Green deformation tensor, with F the deformation gradient. U ( J ) is composed of the bulk modulus κ 0 kPa and a penalty term related to the Jacobian of the deformation gradient J = det(F ) and we choose U ( J ) = κ 2 ln( J ) 2 (2) for all considered material models. One of the simplest constitutive laws to model rubber-like materials is the isotropic Demiray model [16] Algorithm 1 Augmented Sellier's Inverse Method with Armijo strategy. 1: initialize X 0 = x dat ; R 0 = 0; k = 0; β = 1.0 2: initialize min which determines the maximal number of Armijo steps; we used min = 1 8 3: do 4: with the invariant I 1 = tr(C ) and parameters a > 0 kPa and b > 0. Here, C := J −2/3 C is the isochoric part of the right Cauchy-Green tensor resulting from the multiplicative split proposed by Flory [22] to model the nearly incompressible behavior of elastic materials.
Single Fung-type exponential models, e.g., Guccione et al. [31], Usyk et al. [82], can be generalized to where Here, are projections of the Green-Lagrange strain tensor E = 1 2 (C − I ) and parameter a > 0 kPa serves as a stress-like scaling. The dimensionless parameters b ff , b ss , b nn > 0 account for the mechanical behavior of the tissue along the fiber (f), sheet (s) and sheet-normal (n) direction, and b fs , b fn , b ns > 0 account for structural interactions. For parameter values of the different single Fung-type exponential models considered, see Table 1.
Holzapfel and Ogden [39] proposed a separated Fung-type exponential model which can be generalized to with parameters a (•) ≥ 0 kPa and b (•) > 0, invariant I 1 = tr(C ), unimodular fourth-invariants such that contributions of compressed fibers are excluded, and interaction-invariants Note that for the anisotropic contributions in Equation (6) the deformation gradient F and the right Cauchy-Green tensor C remain unsplit to avoid nonphysical results [35,74].   [80] ortho tr opic 0.88 5.0 Comparing to experimental data, Holzapfel and Ogden reduced the constitutive law with the full set of invariants (6) to Further, several papers [18,33,77] deal with in-plane and out-of plane dispersion of collagen fibers along the f 0 and s 0 direction. According to Eriksson et al. [18] this is modeled by modifying the unimodular fourth-invariants to Dispersion parameters have been identified previously by mechanical experiments on passive cardiac tissue by Sommer et al. [77] and are set to κ f = 0.08 and κ s = 0.09.
For parameter values of the considered separated Fung-type exponential models, also often referred to as HO-type models, see Table 2.

Fitting with model function
Due to limited availability of clinical data representing the EDPVR, the passive biomechanical material model is fitted to the empirical Klotz EDPVR estimated from a single measured PV-pair [45]. A more detailed description about the method proposed by Klotz et al. [45] can be found in Appendix B.
Fitting is done by adjusting material parameters of the constitutive relation during the unloading procedure (Algorithm 1) such that the cavity volume of the unloaded reference configuration X * matches V klotz 0 , see Equation (B.1), and such that the simulated PV-curve matches the EDPVR as predicted by Equation (B.2). In particular, the parameters are varied in each unloading step as follows: consider In a second step the procedure is repeated for the simulated PV-curve in the current step k of the unloading algorithm with .., n ls } being the volumes at the different loading points and x 0 the cavitary volume of the current reference configuration X k , to obtain the parameter set [a sim fit , b sim fit ]. The parameter fitting is an extension of the unloading algorithm described above, see also Algorithm 2. Thus, we use the same index k to specify the current iteration step while index i specifies the current loading step of the forward simulation. Note that a sim fit > 0 holds, due to the fact that an increase in pressure leads to an increase in volume, and b sim fit = 0, due to the design of the model function (8). We compute the scalings Here, the absolute value of b sim fit is chosen, as this fitting parameter might be negative in rare circumstances when the simulated PV-curve is concave. Also, for the unlikely case of a linear PV-curve, b sim fit is close to 0, leading to a very large scaling parameter b scale . Hence, to avoid large jumps in the material parameters, we define the bounding interval I = [ 1 5 , 5] with a scale and b scale taking the value of the closest bound if a scale / ∈ I or b scale / ∈ I, respectively. Finally, all parameters of the material model in Equations (3) to (6) The aim of this scaling step is to get the current simulated loading curve to converge against the empiric EDPVR. When the material is too stiff (V klotz 0 < V sim 0 ), then stress-like a (•) parameters are scaled with a value less than one (a scale < 1) to soften the material; vice versa when the material is too soft (V klotz 0 > V sim 0 ), then the stress-like a (•) parameters are scaled with a value greater than one (a scale > 1) to increase stiffness. Further, when the simulated loading curve is too straight, parameters b (•) are increased (b scale > 1); and vice-versa when the simulated loading curve is too bent, then b (•) are decreased (b scale < 1). These changes to the current parameter setting are ensured by the independent fitting of the model function first to the Klotz-curve and second to the current loading curve and the scaling as described in Equations (9) and (10). Interactions between simultaneous changes of a (•) and b (•) parameters are settled with the iterative approach extending the unloading algorithm as described above. Upon convergence of the algorithm in step k, see Section 2.3.3 for convergence criteria, we get the unloaded configuration X * = X k and final parameter set a * (•) = a k , see Algorithm 2 and Fig. 2.
Note that all discussed material models require that a k (•) > 0 and b k (•) > 0 for all k and hence a scale and b scale have to be positive, which is guaranteed by design. Additionally, we point out that for all simulations where the initial loading function was concave the presented algorithm still converged in a satisfactory manner.
To verify the proposed MFF approach, a more costly cost function based fitting (CFF) approach was established and applied to the available patient data. For details see Section 3.4.

Error estimates and goodness of fit
To check for convergence, detect potential stagnation, and measure goodness of fit of the proposed method, several error estimates were introduced, which are described in detail below. The difference of V 0 and V ed between the simulated and the empirical Klotz EDPVR was quantified by the difference of the initial volumes where V sim 0 is the cavity volume of the reference configuration X k at the current unloading step k, and the difference of end-diastolic volumes where V dat ed is the measured end-diastolic volume and V sim ed is the cavitary volume of the inflated configuration x k . For convergence we chose an error tolerance of vol = 0.5% of V dat ed such that r ed < vol and r V 0 < vol .
The difference between the simulated and measured end-diastolic geometry was calculated by consideration of the maximal nodal error where • 2 stands for the 2 -norm and N nodes represents the total number of nodes in the geometry. To obtain convergence, r x ∞ was required to be smaller than 0.1 mm. Finally, we define the error in parameter update between unloading iterations as which is used to detect stagnation of the algorithm and required to be smaller than 0.001 for convergence.
To define the goodness of fit for the simulation outcome in terms of the fitted curve and the Klotz curve, the relative difference of initial volumes and the relative area difference of the normalized curve to the area under the Klotz curve A klotz were calculated. The absolute area difference r A n between two curves in the 2D space is computed by the method introduced by Jekel et al. [40], which positions quadrilaterals q between two curves and subsequently calculates the sum of their areas r A n = q A q , Algorithm 2 Fixed-point algorithm to identify myocardial passive material parameters and the unloaded reference configuration. 1: initialize geometry, X 0 = x dat ; k = 0; β = 1.0 2: initialize PV-pair from measured data p dat ed , V dat ed 3: initialize default parameters with values from the literature, a 0 4: compute empirical EDPVR as proposed by Klotz et al. [45], see Appendix B, using PV-pair p dat ed , V dat ed 5: fit model function Equation (8)   Starting from an initial geometry X 0 obtained from image data and initial parameters a init (•) , b init (•) from the literature, we perform an iterative procedure to find estimates for the unloaded reference configuration X * and personalized passive material parameters a * (•) , b * (•) . While iterating, the simulated loading curve (in green) gets closer to the empiric EDPVR [45] (in red) which is determined by a measured pressure volume pair (p dat ed , V dat ed ). Convergence of the algorithm is based on error estimates given in Section 2.3.3.
using Gauss' area formula where (x i , y i ) represent the vertices of a quadrilateral.

Numerical framework
The unloading and parameter estimation scheme was implemented in the FE framework Cardiac Arrhythmia Research Package (CARP) [83]. The simulations were constrained by applying spring-like boundary conditions at the rim of the clipped aorta and over the epicardial surface at the apex of the LV, respectively, see [49,79]. For all the passive inflation simulations we rely on solver methods established previously [9] which have been verified in a N-version benchmark study [47]. In brief, the linear systems were solved by a generalized minimal residual method (GMRES) with a relative error reduction of = 10 −8 . Efficient preconditioning was based on the library PETSc (https://www.mcs .anl .gov /petsc/). For all cases we used 100 loading steps to ramp up the pressure from p = 0 to p = p ed . To speed up computations, we limited the number of Newton steps to two for the passive inflation simulations during the unloading scheme. This did not alter simulation results and is justified as we were only estimating the unloaded reference configuration and material parameter sets. For the MFF with the Levenberg-Marquardt least-squares algorithm we used Python V. 3.8 and packages NumPy V. 1.16.6 and SciPy V. 1.2.3; see [55] for details on the implementation of this algorithm.
To validate our results, we performed a passive inflation experiment that starts with the found unloaded reference configuration and the determined parameter set. For this passive inflation we used a fully converging Newton algorithm Table 3 Fitting results for all N = 19 cases are listed in terms of final scaling parameters, a scale and b scale after convergence of the algorithm, and measures of goodness of fit, r V 0 ,rel and r An ,rel , along with input parameters, V dat ed and p dat ed , used to compute the Klotz EDPVR. Mean values and standard deviation (SD) of goodness of fit measures were computed separately for each etiology. The whole set of fitted parameters for these cases is given in Table C with a relative 2 -norm error reduction of the residual of = 10 −6 and 100 loading steps. All goodness of fit measurements in the sections below were calculated using the PV-curve of this validation experiment.

LV patient cohort results
The MFF method as described in Section 2.3.2 was executed for all N = 19 patient cases. Parameters of the reduced HO constitutive law were fitted using initial values from the literature, see Table 2. P1-P0-elements with the bulk modulus κ = 650 kPa were utilized for all cases and simulations were run on the Vienna Scientific Cluster 4 (VSC-4) using 96 cores. Given the convergence criteria as presented in Section 2.3.3, the MFF algorithm converged for all 19 cases and fitting results are summarized in Table 3. Cases with best (06-CoA) and worst (03-AS) fit were distinguished by calculating the relative area difference r A n ,rel , see Equation (16), and used for further analysis in the sections below. The Klotz curve and the simulated PV-curve for these two cases are visualized in Fig. 3. Even for the worst case the fitted model approximated the Klotz curve remarkably well, with an almost exact match in terms of V klotz 0 and V dat ed . Indeed, the MFF method achieved excellent goodness of fit for all cases.
Run-times t ul for the MFF method which include unloading and parameter fitting were between 10.85 min and 74.18 min. For the subsequent passive inflation experiment for validation the run-times t val were between 6.31 min and 24.50 min. As expected, computational cost increased with mesh size, see Fig. 4.
Further, the influence of the input parameters V dat ed and p dat ed on goodness of fit, more precisely r A n ,rel , is shown in Fig. 5. No trend can be recognized for different end-diastolic volumes V dat ed while it appears that lower end-diastolic pressures p dat ed might lead to better fits; this is investigated further in Section 3.3.

Results for different constitutive laws
The default constitutive law in Section 3.1, i.e., the reduced HO law with initial material parameters from Guan et al. [30], was chosen simply because this study is one of the most recent papers with fits to data from human ventricular myocardium. In this section, we show that the method described above works equally well for a large variety of constitutive laws, in fact most of the laws that are currently used to model passive myocardium.
To this end, we applied the MFF method using the different constitutive laws listed in Tables 1 and 2 for the case with best (06-CoA) and the worst (03-AS) fit in Section 3.1. P1-P0-elements were used and all simulations were run on VSC4    5. Influence of input parameters V dat ed (left) and p dat ed (right) on relative area difference r An,rel is shown. Data from CoA cases is marked as dots, whereas data from AS cases is marked as crosses.
using 96 cores. Normalized fitting results for cases 06-CoA and 03-AS are given in Fig. 6 and Table 4 showing excellent fits and almost similar PV-curves for all material laws. Due to its simplicity and isotropic nature, the Demiray material model resulted in the lowest number of iterations and shortest run-times for both example cases. More complex orthotropic models needed more iterations of the MFF algorithm and thus also run-times increased. Overall, MFF simulations were tractable with a computational run-time below 2.5 hours for all experiments.

Results for locking-free finite elements
It is well known that simple P1-P0-elements may suffer from locking effects and hence other FE formulations are required for certain applications where accurate stresses are essential. To show the capabilities of the MFF method in this   scenario, we applied the algorithm to cases 06-CoA and 03-AS using stabilized, locking-free P1-P1-elements [44] and an incompressible material, i.e., 1/κ = 0. Normalized fitting results for the two example cases are presented in relation to results obtained above using P1-P0-elements in Fig. 7 and Table 5. The PV-curves and the measures of goodness of fit r V 0 ,rel and r A n ,rel show that the method works equally well for stabilized P1-P1-elements. Looking at the parameter values in Table 5 we can see that a scale is smaller and b scale is larger for P1-P0-elements compared to the locking-free elements. This is anticipated as the fitting compensates for a certain degree of locking in the P1-P0-element formulation and thus it gives parameters that correspond to a softer material. Unsurprisingly, run-times were significantly longer when a stabilized P1-P1 FE formulation was used: 32 min vs 246 min for case 06-CoA and 61 min vs 364 min for case 03-AS, whereas iteration numbers of the MFF algorithm were similar. We could achieve similar results with locking-free MINI elements as introduced by Karabelas et al. [44]; only run-times varied a bit compared to the stabilized P1-P1 formulation.

Sensitivity analysis
To show robustness of the proposed MFF methodology, a sensitivity analysis was performed to investigate the influence of clinical data uncertainty in terms of the end-diastolic PV-point used for computation of the Klotz EDPVR, changes in myocardial fiber orientation, altered model parameters used as initial guess and variations of the bulk modulus κ, respectively.
The sensitivity analysis was executed for example case 03-AS using the reduced HO material law and P1-P0-elements if not mentioned otherwise.

Influence of variations in the end-diastolic PV-point on fitting outcome
The behavior of the Klotz curve was investigated changing the inputs V dat ed and p dat ed independently and subsequently comparing the normalized curves. For that purpose, n = 20 evenly spaced values in the range of 93-282 mL for V dat ed and 0.61-2.80 kPa for p dat ed were chosen; range boundaries correspond to the respective minimum and maximum values in the patient cohort data. For the variable held constant the mean values (V dat ed = 175 mL and p dat ed = 2.26 kPa) obtained from the cohort data were taken. First, V dat ed was varied while keeping p dat ed constant and vice versa. Normalized results of the generated curves are visualized in Fig. 8. It shows that curves generated with alternating V dat ed match when normalized, whereas normalized curves resulting from alternating p dat ed differ. For that reason, varying V dat ed will have little effect on the goodness of fit of the MFF method. We could also observe this behavior previously in Fig. 5. Hence, the sensitivity of the proposed MFF method was studied for varying p dat ed only, considering deviations of ±10%. Results in Fig. 9 and Table 6 suggest that fitting results are better for Klotz EDPVR that have higher normalized pressure values, especially in the lower normalized volume region. Material parameters change for different values of p dat ed , but no particular trend is recognizable. Sheet and sheet-normal directions were adapted accordingly to preserve the orthogonal system of local fiber coordinates. Results of fitting performance are visualized in Fig. 9 and values of fitted parameters and measures of goodness of fit are summarized in Table 7. As expected, parameter values are changing with the fiber orientation, as myocyte directions have a great impact on local material stiffness. Further, we can observe that the goodness of fit of the MFF algorithm is almost indifferent to changes in fiber orientation, leading to almost the same PV-curves in Fig. 9.

Influence of initial model parameters on fitting outcome
To show the robustness of the method, the influence of variations in initial model parameters of the reduced HO material model was investigated. LHS from pyDOE V. 0.3.8 was used to create n = 10 different sets of initial scaling parameters {a scale , b scale } in the interval (0, 1). For all executed simulations the exact same parameter set and goodness of fit were obtained, also matching results with the default initial parameter set {a scale = 1, b scale = 1}, see case 03-AS in Table 3. While this shows the high robustness of the MFF algorithm, we cannot prove uniqueness and the run-times of the simulations varied significantly: t ul = 39.48−302.25 min; depending on how close the initial parameters are to the final parameter set.

Influence of variations in bulk modulus on fitting outcome
Finally, sensitivity of the method to changes of the bulk modulus κ was studied using both P1-P0-elements and lockingfree, stabilized P1-P1-elements. Simulations were run for example case 03-AS for κ = {1000, 3000, 5000} kPa. Additionally, default values were chosen for P1-P0-elements (κ = 650 kPa) and for P1-P1-elements (1/κ = 0). Results are visualized in Fig. 10 and summarized in Table 8. The forward simulations were not robust for P1-P0-elements and κ = 5000 kPa resulting in an irregular curve in Fig. 10(left). This is a shortcoming of the penalty formulation and P1-P0-elements, not the MFF method itself. We can see that for P1-P0-elements the material parameters vary greatly with the bulk modulus κ, compensating locking effects with material parameters that correspond to a softer material. However, when using P1-P1elements, the method was not sensitive to changes in κ. Material parameters and goodness of fit are similar, which is a further indication of the robustness of the MFF algorithm.

Comparison to fitting with cost functional and parameter space sampling
To compare and verify the results of the MFF approach we implemented an ad hoc optimization procedure, based on cost functional fitting (CFF), to find sets of material parameters. To penalize deflections at lower and higher pressures and Here, the first term is the sum of squared normalized volume differences, the second term is the sum of squared normalized pressure differences, and the third term serves to further ensure that V klotz 0 and V sim 0 are close, with γ > 0 a weighting parameter; in all presented simulations we chose γ = 1. A downhill-Simplex based fitting (Nelder-Mead method) from SciPy was used to carry out the optimization, where we chose an absolute error of 0.001 in the parameter update between iterations to be acceptable for convergence. We tried different varieties of CFF, the above described method with Equation (17) proved to be the most efficient.

Results
The CFF approach was executed for the two example cases (06-CoA and 03-AS) with the reduced HO constitutive law, P1-P0-elements, and the bulk modulus κ = 650 kPa. We used default values for the material parameters as in Table 2 scaled by a factor of 0.5 as we already expect from results in Table 3 that material parameters from the literature are in general too stiff; this scaling sped up computational times considerably.
Using the Nelder-Mead algorithm we optimized (i) n opt = 1 scaling parameter (a scale ) which was applied to all material parameters a (•) , whereas all b (•) parameters of the model were scaled with 0.5; (ii) n opt = 2 scaling parameters (a scale , b scale ), with a scale applied to all a (•) parameters and b scale applied to all b (•) parameters of the model; and (iii) n opt = 3 scaling parameters (a scale , b iso scale , b aniso scale ), where a scale is the scaling parameter that was applied to all a (•) parameters; b iso scale is the scaling parameter that was applied to b in the isotropic contribution of the constitutive law; and b aniso scale is the scaling parameter that was applied to all b (•) parameters of the anisotropic contribution.
We used bounds such that all scaling parameters are > 0.1 to prevent that parts of the constitutive law are eliminated. The case n opt = 4, where additionally a from the isotropic contribution and all a (•) from the anisotropic contribution were scaled by two different parameters, was not converging for the patient-specific cases within a reasonable amount of time (72 hours). We performed the optimization using values from the literature, see Table 2, as starting values (no init) as well as using the MFF approach to generate initial guesses for the parameter scalings (init).  11. Normalized fitting results of case 06-CoA (left) and case 03-AS (right) executed using the reduced HO constitutive law with initial guesses from literature (no init) and with a MFF initialization step (init) as initial guess. Outcomes of the CFF method, optimizing one, two, or three parameters, n opt = {1, 2, 3}, are compared.

Table 9
Fitting results in terms of fitted scaling parameters and measures of goodness of fit for the two example cases executed using the reduced HO constitutive law with parameters from literature [39] (no init) as initial guess and with an MFF initialization step (init). Outcomes of the CFF method, optimizing one, two, or three parameters, n opt = {1, 2, 3}, are listed with number of iterations (it.) of the Nelder-Mead algorithm.

Fitted Parameters Goodness of Fit Fitted Parameters Goodness of Fit
Mode Results are summarized in Table 9 and in Fig. 11 and we see that the worst fit was obtained for case 06-CoA optimizing two parameters and using initial values from the literature (no init), where the Nelder-Mead method got stuck in a local minumum. For n opt = 1 similar results were reached when using values from the literature (no init) and initial scalings generated by a MFF (init). Best fits were acquired for n opt = 2 with MFF initialization and n opt = 3 in both variants. For case 03-AS obtained fitted parameters and goodness of fit were very similar for both initialization variants. Overall, results did not deviate much when one, two, or three parameters were optimized. In general, simulations run with default parameters from the literature needed more Nelder-Mead iterations compared to simulations with MFF iterations and thus also runtimes were longer (>1.5 times as long). (Table 3) and CFF (Table 9) agree well in terms of fitted parameters and measures of goodness of fit. Also the PV-curves of the MFF and CFF approach match almost exactly in Fig. 12 where the MFF result is compared to the best fit of the CFF method. However, since a full unloading step has to be solved in each Nelder-Mead iteration, the number of passive inflation simulations for the CFF approach is significantly larger compared to the MFF method: 79/156/516 vs 8 for 06-CoA and 110/460/835 vs 10 for 03-AS for CFF no init and n opt = 1/2/3 fitted parameters vs. MFF; thus also the run-time of the CFF approach was considerably longer compared to the MFF approach: 330/546/2095 min vs 32 min for 06-CoA and 541/2197/3974 min vs 61 min for 03-AS for CFF no-init and n opt = 1/2/3 fitted parameters vs. MFF.

Comparison to other methods and the state-of-the-art
This study describes a novel methodology to fit passive mechanical parameters of the myocardium. Additionally, the algorithm computes an unloaded reference configuration where a new line search strategy was implemented that improves the robustness of a previously introduced backward-displacement method. We have demonstrated that the presented MFF method is an efficient, robust, and versatile approach for the automated identification of model parameters from clinical image data using the empirical EDPVR as proposed by Klotz et al. [45] as a fitting target. With this, the MFF algorithm only requires image or mesh data from one time instance during diastolic filling and a single measured data point of the EDPVR as inputs. Hence, the MFF method does not depend on a known approximation of the reference configuration which is a limitation of previous works [34,58]. This is of particular interest as often only a single anatomical snapshot taken within the diastatic window is available from clinical studies. Alternatively, previously generated geometries can be used for the MFF algorithm and thus no additional segmentation work is required to estimate material parameters for already existing meshes. Note that in many cases pressure and volume at end-diastole and image acquisition are not synchronized as both involve different modalities, e.g., catheter-based pressure measurement vs. echocardiography/CT/MRI. However, the presented method will still work if the image is scanned during early diastole or diastasis given an estimate of the pressure, p image , in the image-based geometry. Then, the first step of the unloading algorithm will be performed using p image and the maximal nodal error comparing simulated and measured geometries Equation (13) is computed at this pressure. See also Finsberg et al. [21] for a similar workflow.
Contrary to other studies that jointly estimate material parameters and the unloaded reference configuration [4,21,20,62], we did not require to fix any parameters of the material law. To the best of our knowledge this is thus the first study that allows to fit not only for the predicted volume of the reference configuration V 0 but also for the shape of the EDPVR. In opposition to studies above, we could apply our algorithm to much finer meshes (up to 1.5 million elements) and still keep simulations tractable with regard to computational cost. This is essential for many translational application of computational models in industry or in the clinic where time constraints apply.
Parameter estimation results are in line with previous studies where experimental parameters from the literature were also often considered as too stiff [41,72,84]. Fittings are in a similar range compared to [26, Table 1], [34, Figure 10], and [62, Table 2] for the Guccione law; and [4, Table 6] and [21, Table A1] for the one-fiber HO law; in these studies the authors were also able to reproduce the Klotz curve. Comparison of parameter values presented in other studies is difficult as these did not take the Klotz relation into account [58,90] or used different ratios between parameters [72].
To verify our novel MFF approach we implemented an ad hoc CFF based optimization procedure. Here, MFF and CFF approaches yield similar fitting results distinguishing the MFF approach as the superior method due to significantly shorter run-times and a higher robustness. While the MFF approach converged for all examined patient-cases, the more expensive CFF approach did not. This was due to time constraints on clusters and numerical problems stemming from instabilities of the forward solver owing to very soft material parameters chosen by the Nelder-Mead algorithm. For the latter case, we returned a very high cost, however, for some patient-specific geometries this still led to convergence problems. We further showed that even if a CFF approach is pursued, the MFF method can be used as a fast approach to generate an initial guess to improve convergence and efficiency. Note that for the CFF approach with n opt = 3 scaling parameters, we fitted the material parameters of the isotropic and the anisotropic contribution independently. Nevertheless, we are aware that it is hardly possible to estimate the degree of anisotropy based solely on pressure-volume data. Here, experimental setups -as, e.g., described in papers mentioned above for default material parameters (Tables 1 and 2) -are by far better suited. For the MFF approach we thus keep the degree of anisotropy as it was given by default parameter sets and we do not distinguish between isotropic and anisotropic parts. For the sake of comparison, we tried the CFF approach also with 3 and 4 fitting parameters. As we can see in Table 9 the CFF approach with 3 fitting parameters did not improve results. Further, the CFF approach with 4 fitting parameters did not converge within a reasonable amount of time. In regard of the methods discussed in this paper we thus conclude that it is not worth pursuing an independent fitting of the isotropic and anisotropic contribution and, consequently, it is better to keep the degree of anisotropy as it was determined in experiments.
There are several other works that discuss the determination of the reference configuration from a loaded state, e.g., [29,48,67]. Each of these algorithms, as well as Algorithms 3 and 4, can be used instead of Algorithm 1 as a basis for the presented MFF approach to fit passive mechanical parameters of the myocardium. We only need the loading curve generated by inflating the unloaded reference configuration to a prescribed pressure as an input.
Finally, compared to some above-mentioned studies, the proposed method can be coupled with relative ease to established FE software packages. It only requires a limited number of passive inflation simulations, a standard experiment in all continuum mechanics simulators, and basic least-squares fitting tools which are available through open-source packages such as SciPy.

Computational costs
For all patient-specific cases, run-times for the MFF method (between 10.85 min and 74.18 min) were only 2 to 3 times higher compared to a standard passive-inflation experiment. Note that a further speed-up is possible by using only one Newton step during unloading, making the MFF method sometimes even faster than a single passive inflation experiment. However, in this case we observed stair-casing effects in the simulated PV-curve that might affect results to a minor degree. Also, a considerably lower number of loading steps will accelerate simulation times. In this scenario, the least-squares fitting to the piece-wise linear PV-curve could be less accurate.
The high efficiency of the parameter fitting is of utmost importance as the computational burden imposed by highresolution patient-specific models demands fast numerical methods to keep simulations tractable.

Versatility of the workflow with respect to material laws and FE formulations
The MFF approach works for various material models including the most widely used passive cardiac tissue models, e.g., HO-type and Fung-type materials as introduced in Section 3.2. Run-times and iteration numbers of the MFF algorithm varied to a small degree with complexity and anisotropy of the model. Convergence with criteria discussed in Section 2.3.3 was reached for all constitutive laws within 10 iterations. It is worth mentioning, that the MFF using the Klotz curve as a target did not work for the isotropic neo-Hookean and Mooney-Rivlin materials. Though sometimes used for cardiac mechanics [15,73], also Hadjicharalambous et al. [34] reported previously that with these constitutive laws the Klotz curve could not be reproduced. This is suggestive of such materials being inappropriate to model cardiac tissue. The problem for the fitting in our case is that the PV response with these materials results in a concave curve -as also observed in [34, Figure 12] -while the Klotz curve is convex. Hence, the presented MFF failed, however, our approach might still work with these materials for other applications than cardiac mechanics with target functions that are not necessarily convex.
We could further show that the MFF approach can be used with different FE formulations such as P1-P0 elements and locking-free stabilized P1-P1 and MINI elements. Here, the locking-free elements resulted in material parameters that correspond to a stiffer material as the fitting compensates for locking effects in the P1-P0-element formulation. As expected, run-times are significantly higher for locking-free elements; still the efficiency of the MFF approach allowed for unloading and parameter estimation within reasonable time frames.

Goodness of fit and sensitivity analysis
We were able to show that our method allows to match the predicted volume of the unloaded geometry V klotz 0 and the given end-diastolic volume V dat ed almost perfectly while reproducing the shape of the PV-curve within uncertainty of the empiric Klotz relation (Table 3). It is important to note that the EDPVR cannot be accurately determined in vivo for a number of fundamental limitations. First, the LV volume during diastolic filling depends on the pressure difference between LV cavity and intra-pericardial pressure and is also influenced by the presence of time-varying non-zero active stresses. Time-varying intra-pericardial pressure is non-negligible and influenced by intra-thoracic pressure modulated by breathing and, as such, cannot be recorded easily beyond specifically designed experiments [38,37]. Further, tails of the active stress transients during diastole depend on cytosolic calcium, the state of Troponin-C buffering as well as the spatially varying sarcomere geometries. None of these factors can be monitored. As such, the time course of active stresses during diastole must be considered unknown, albeit efforts have been invested to estimate these [91]. The Klotz relation has been shown to be robust across species and a range of pathological conditions under ex vivo conditions where both volume and intra-cavitary pressure can be controlled with higher accuracy, but, owing to its empirical nature, must be considered an approximation of the EDPVR in an individual patient under in vivo conditions.
Independently of the uncertainties a Klotz-based surrogate EDPVR is afflicted with, the application to the patient cohort resulted in excellent fittings to the Klotz curve for all cases with a goodness of fit measure r A n ,rel ranging between 7.38% and 14.77%. Also r V 0 ,rel was very low for all cases, although a small trend towards better fits in the CoA cases was observed. Note that only for these cases the measured pressure was available, see also Section 3.3.1.
The influence of the input pressure p dat ed was noticeable, resulting in different shapes of the Klotz EDPVR (Figs. 8 and 9) and therefore in different material parameters. This influenced the goodness of fit to some extent, hence, the sensitivity to this measure should be considered, especially in the case when invasively measured pressures are not available.
Fiber orientations are known to have a great impact on cardiac mechanics simulations [25,71], thus, the influence of variations of the fiber angle α on fitting performance is an integral point of the sensitivity analysis. Variations of fiber angles resulted in very similar outcomes in terms of goodness of fit as well as computational cost. The altered stiffness due to altered fiber directions was compensated by the variation of material parameters, hence resulting in similar PV-curves.
Changes in initial guesses for model parameters were not reflected in the fitting outcome at all and resulted in the same fitted material parameters and curves. Only the run-times for unloading varied when using initial guesses lying further from the final result. This robustness of the approach is essential for the interpretation of material parameters for characterizing patient pathology and understanding changes in material properties under HF conditions. The fitting outcome was not sensitive to variations of the bulk modulus κ when using locking-free elements. However, for P1-P0-elements the shape error increased for higher values of κ, compensating for the increased locking of the P1-P0elements. As in all simulations modeling nearly-incompressible behavior of soft tissue, special care has to be taken when setting this parameter for simple linear elements.
Finally, the choice of the model function (8) was motivated by the mathematical structure of the constitutive laws and the Klotz relation. In Appendix D, we could demonstrate that several other model functions can be used for the fitting, all leading to excellent results. As the MFF approach is applicable to many kinds of target loading curves other than the Klotz relation, the definite choice of the model function is specific to the problem.

Limitations
First, the workflow presented in this paper relies on the empiric Klotz relation. However, the approach also works in the exact same fashion for other experiment-based PV-curves, discrete data points obtained from loading experiments, or experimental/clinical in vivo measurements of the EDPVR.
Second, as residual strain is not considered, in general, this method does not generate a uniform fiber-stretch field which is generally assumed at end-diastolic state [14,23,65]. Heterogeneity in fiber stretch in an end-diastolic state impairs the Frank-Starling mechanism, as shown in a recent modeling study [8]. In future studies, the presented workflow will be extended by using ideas based on growth and remodeling [27,88] to address this issue.
Third, we presented the methodology only for single-chamber LV models, while solid-and electromechanical wholeorgan simulations of the heart are becoming feasible [9,11,72,78]. Nevertheless, the workflow was successfully tested for bi-ventricular models [42] as well as 4-chamber models. Here, we used one homogeneous material for the ventricles and default values from the literature for the atria, e.g., [7], and applied patient-specific pressures in the different chambers for the forward solve in Algorithms 1 and 2. We did not present these simulations as an empirical estimation for the EDPVR was only available for the LV given as the Klotz law. To the best of our knowledge, no comparable approximation of EDPVR exists for the right ventricle or the atria.
Finally, we cannot provide a rigorous proof that the resulting unloaded configuration and material parameter sets are unique. However, we could show that the method was robust with respect to initial values leading to the same result for all simulations, see Section 4.4. Further, as we are fitting a strictly convex function (8), we can assume that the Levenberg-Marquardt algorithm results in a unique parameter set [13]. Since the MFF method requires that the estimated parameters are close to this solution, see Equation (14), we can expect that the MFF method is robust with regard to material parameters.

Conclusion
We report on a novel MFF approach combining the fitting passive mechanical parameters of the myocardium to the shape of the EDPVR and the simultaneous estimation of the unloaded reference configuration. The algorithm only requires image or mesh data from one time instance during diastolic filling and a single measured PV data point of the EDPVR as inputs. The MFF is efficient, robust, and versatile and can be applied to reproduce clinically-relevant PV relationships for patient-specific LV anatomical models within clinically easily feasible time frames in a fully automated fashion. Thus, the method constitutes a further step forward towards a realistic representation of LV passive mechanical function. As such, the MFF is a pivotal building block in workflows for building computational digital twin models of human cardiac EM function at scale, to facilitate the generation of virtual cohorts in translational applications.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix B. Klotz relation
Considering the limited availability of clinical data of the EDPVR, the computational method proposed by Klotz et al. [45] is utilized, which enables prediction of the EDPVR by a single measured PV-pair. According to this seminal work the volume of the unloaded geometry V klotz 0 of the LV can be empirically determined by where V dat ed and p dat ed is a measured PV-pair at end-diastole. Further, the EDPVR is described by the power law p = αV β , (B.2) where p is the cavity pressure in mmHg, V is the cavity volume in mL, and the constants α and β are defined by the relations α = , which applies to all patient cases in the CARDIOPROOF cohort (see Table 3). Fitted parameters for all N = 19 cases for the reduced Holzapfel-Ogden law; for this default material parameters given in Table 2 are multiplied according to Equation (10) with final scaling parameters given in Table 3.

Appendix D. Choice and comparison of model functions
As default we chose the model function in Equation (8) for its similarity to the constitutive laws and to keep b = 0.
Nevertheless, the fitting method also worked for other model functions, e.g., those listed below: (energy-Demiray) V klotz (p) − V sim (p) . (D.1) Hence, for certain cases it can be a good alternative to the (energy-Demiray) function. Model functions i given above are only a small subset of functions that we tried for our fitting to the Klotz curve. All of them showed satisfying results and convergence rates but other functions might work as well for the procedure. Since the Klotz curve resembles an exponential function, it is not surprising that the (energy-Demiray) and (exp-function) worked best. However, as our method would work for all kinds of target EDPVR other than the Klotz law, the choice of the model function is specific to the problem.