Improved identifiability of myocardial material parameters by an energy-based cost function

Myocardial stiffness is a valuable clinical biomarker for the monitoring and stratification of heart failure (HF). Cardiac finite element models provide a biomechanical framework for the assessment of stiffness through the determination of the myocardial constitutive model parameters. The reported parameter intercorrelations in popular constitutive relations, however, obstruct the unique estimation of material parameters and limit the reliable translation of this stiffness metric to clinical practice. Focusing on the role of the cost function (CF) in parameter identifiability, we investigate the performance of a set of geometric indices (based on displacements, strains, cavity volume, wall thickness and apicobasal dimension of the ventricle) and a novel CF derived from energy conservation. Our results, with a commonly used transversely isotropic material model (proposed by Guccione et al.), demonstrate that a single geometry-based CF is unable to uniquely constrain the parameter space. The energy-based CF, conversely, isolates one of the parameters and in conjunction with one of the geometric metrics provides a unique estimation of the parameter set. This gives rise to a new methodology for estimating myocardial material parameters based on the combination of deformation and energetics analysis. The accuracy of the pipeline is demonstrated in silico, and its robustness in vivo, in a total of 8 clinical data sets (7 HF and one control). The mean identified parameters of the Guccione material law were \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_1=3000\pm 1700\,\hbox {Pa}$$\end{document}C1=3000±1700Pa and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha =45\pm 25$$\end{document}α=45±25 (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b_f=25\pm 14$$\end{document}bf=25±14, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b_{ft}=11\pm 6$$\end{document}bft=11±6, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b_{t}=9\pm 5$$\end{document}bt=9±5) for the HF cases and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_1=1700\,\hbox {Pa}$$\end{document}C1=1700Pa and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha =15$$\end{document}α=15 (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b_f=8$$\end{document}bf=8, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b_{ft}=4$$\end{document}bft=4, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b_{t}=3$$\end{document}bt=3) for the healthy case.


Introduction
Left ventricular (LV) stiffness is proposed as a diagnostic indicator of cardiac function in heart failure (HF) patients (Westermann et al. 2008). Ventricular stiffness has been predominantly assessed in clinical practice through pressurevolume (p-V) analysis (Bermejo et al. 2013;Burkhoff et al. 2005;Zile et al. 2004). However, this approach is unable to distinguish between the anatomical and material contributions to LV stiffness. Specifically, an increment in ventricular size due to myocardial hypertrophy or an increase in collagen content with fibrosis may both lead to an equivalently stiffer LV behaviour using this methodology. Differentiating between these two components, anatomical and material, may improve the identification of HF aetiology in patients.
The development of biophysical models (Chen et al. 2016;Crozier et al. 2016;Krishnamurthy et al. 2013;Lee et al. 2015;Nordsletten et al. 2011;Plank et al. 2009) for the simulation of cardiac mechanics allows the distinct representation of the geometric and material components of stiffness. Using these models, the assessment of myocardial stiffness is posed as an inverse problem, where the material parameters are determined from known mechanical loads and deformations. Recent research in this field has focused on developing tractable pipelines for linking model parameters to data (Augenstein et al. 2005;Wang et al. 2009), evaluating the available material models (Criscione et al. 2002;Schmid et al. 2009) or improving the optimization strategies (Balaban et al. 2016;Moireau and Chapelle 2011;Moireau et al. 2008Moireau et al. , 2009Nair et al. 2007) and has led to the estimation of material parameters from clinical data sets Gao et al. 2015;Wang et al. 2013;Xi et al. 2013). An inherent limitation in these current methods is the intercorrelation of the material parameters in myocardial material laws (Augenstein et al. 2006;Gao et al. 2015;Remme et al. 2004), which results in multiple parameter combinations corresponding to equivalent solutions in the optimization process. The existence of multiple solutions for the inverse problem limits the interpretation of these parameters for characterizing patient pathology and understanding changes in material properties under conditions of HF.
In this paper we investigate the role of the cost function (CF) in parameter identifiability and develop a novel energy-based CF that allows us to uniquely constrain the myocardial material parameters. For our analysis we choose a popular material model in cardiac mechanics, the transversely isotropic constitutive equation proposed by Guccione et al. (1991), which is reported to suffer from parameter coupling (Augenstein et al. 2006;Xi et al. 2011). We examine how a group of CFs based on geometric attributes, and the energy-based CF, constrain the optimization in the search of the parameters that best explain the clinical data of pressure and deformation. After an evaluation on a synthetic data set, a novel parameter estimation pipeline emerges based on the combined use of the energy-based CF with one of the geometric CFs, and is tested in 8 clinical cases demonstrating its ability to identify unique material parameters from patient data.

Methods
In this section we summarize the synthetic and clinical data sets used (Sect. 2.1), the modelling framework (Sect. 2.2), the evaluated CFs (Sect. 2.3) and the proposed parameter estimation pipeline (Sect. 2.4). All data processing has been performed in MATLAB, and the meshes and simulation outputs have been visualized with cmGui 1 (Christie et al. 2002).

Synthetic
To provide a known ground truth for the material parameters a synthetic data set was employed (see top panel in 1 www.cmiss.org/cmgui. Fig. 1). An in silico model of the LV diastolic mechanics was created from the passive inflation of a truncated confocal prolate spheroidal of typical human cardiac dimensions representing the myocardial domain (Evangelista et al. 2011;Ho 2009;Humphrey 2002) to an end-diastolic pressure of 1.5 kPa (Humphrey 2002). A mesh of 320 (4 transmural, 16 circumferential, 4 longitudinal and 16 in the apical cap) hexahedral elements and 9685 nodes was used for the passive inflation simulation (details on the interpolation schemes and solver used are provided in Sect. 2.2.4). The pressure was applied over 30 equal pressure increments of 0.05 kPa, keeping the nodes of the 'basal' plane fixed in all directions. The resulting 31 meshes (undeformed mesh and 30 deformed meshes from each pressure increment) and their corresponding cavity pressure values from the simulation compose the synthetic data set used for the in silico study.

Clinical
In this study 8 clinical data sets are utilized, obtained from 7 Cardiac Resynchronization Therapy (CRT) patients (denoted as PC1-PC7) and one healthy control (denoted as HC). The clinical profile of the 8 cases is shown in Table 1. PC1-PC7 were obtained according to the clinical protocols followed in St Thomas' Hospital, London, and consist of LV cavity pressure recordings and cardiac images covering the entire cardiac cycle. The data collection conforms to the principles of the Declaration of Helsinki and is guided by a local ethics committee approved protocol with patient informed consent. The healthy data set consists of pressure data and LV meshes covering diastole and were described previously (Xi et al. 2013).
The cardiac images of the CRT patient data sets (PC1-PC7) consist of 2-D short axis stacks of cine MRI with SENSE encoding (1.19 × 1.19 × 8 mm 3 to 1.45 × 1.45 × 10 mm 3 resolution), taken on a 1.5T-in six out of seven cases-or 3T-in one case-Achieva Philips Medical Systems MRI scanner. Each MRI sequence had 25 to 35 frames with a temporal resolution between 23 and 32 msec. The LV domain was manually segmented in itksnap 2 from the end-diastolic frame. Images were processed using a nonrigid registration ) which enables a spatially and temporally continuous description of the cardiac motion. Mesh personalization was performed on the segmented LV, as described previously . A set of deforming finite element (FE) meshes, consisting of 12 to 16 (4 circumferential, 3 to 4 longitudinal and 1 transmural) cubic hexahedral elements and 436 to 580 nodes, were created for each patient by warping the personalized end-diastolic anatomical mesh using the motion field corresponding to each frame of the cine MRI. As a result, correspondence Synthetic and clinical data sets used to evaluate CFs. a Synthetic data set was created by applying 30 equal pressure increments to an idealized finite element (FE) model to generate 30 deformed geometries 'frames', b the clinical data set combines an averaged pressure trace with FE meshes that capture the deformation calculated from registra-tion of the cine MRI frames. Combined, these provide a displacement and pressure measurement for each MRI frame recorded over the cardiac cycle. In our analysis only the diastolic frames, where a passive inflation approximation is relevant, are utilized (marked as diastolic window of interest) of material points between frames is obtained from the cine MRI images through the image registration and mesh personalization processes.
The LV cavity pressure transient was recorded during a catheterization procedure before the beginning of the CRT pacing protocols and separately from the MRI scans. For each patient an average pressure trace was calculated over 5-13 beats and then synchronized to the cavity volume trace estimated from the personalized FE meshes. The pressurevolume synchronization was based on the assumption that the inflection point in the pressure wave is approximately aligned with the R wave (acquisition time of the first frame of each MRI sequence) and finding the temporal offset that maximized the p-V loop area and was less than 5% of the R-R interval. The pressure transient was subsequently offset to ensure a zero pressure at the MRI phase that corresponded to the approximated reference configuration for the finite elasticity analysis, described below. The schematic of the steps followed for the processing of the clinical data sets is shown in the lower panel of Fig. 1.

Mechanical model
LV diastolic passive inflation is simulated using large deformation mechanics assuming that deformation is driven principally by the LV pressure, the myocardium has homogeneous material properties, is incompressible, inertia or viscoelastic effects are negligible, that the LV is in a stable relaxed state in late diastole and that the right ventricle

Cardiac microstructure
The myocardial microarchitecture requires the continuous description of the local fibre ( f ), sheet (s) and sheet normal (n) directions. In the model, local tissue microstructure was described by assuming a linearly varying preferential myocyte orientation from −60 • at the epicardium to 60 • at the endocardium based on the findings by Streeter et al. (1969).

Material description
The myocardium is modelled as a hyperelastic incompressible transversely isotropic material with the constitutive relation introduced by Guccione et al. (1991). The mathematical description of the Guccione law is given in Eqs. 1 and 2. The parameters b f , b t , b f t assign different mechanical responses to the tissue along the fibre ( f ) direction, across the transverse planes (t) and in the fibre-transverse shear planes ( f t), respectively. The f ,s,n indices in the strain components express projections of the Green-Lagrange strain tensor (E) along the fibre ( f ), sheet (s) and sheet normal (n) directions.
After the reformulation proposed by Xi et al. (2013), the strain energy density (Ψ ) is expressed as a function of the scaling (C 1 ) and bulk exponential (α) parameters, along which the primary coupling occurs (Eqs. 1, 3). (4) The r f , r f t , r t parameters are referred to as anisotropy ratios and range between 0 and 1, with the r f ratio obtaining usually the highest value in order to represent a stiffer behaviour along the fibre direction. Throughout our analysis the anisotropy ratios were kept constant at r f = 0.55, r f t = 0.25, r t = 0.2 while focusing on fitting the coupled C 1 and α parameters.

Reference configuration
The reference configuration represents an idealistic stressand strain-free geometry for the myocardium which is never reached within the cardiac cycle. For simplicity the LV geometry associated with the MRI frame corresponding to the minimum pressure was chosen as an approximation of the reference geometry.

Mechanical simulations and boundary conditions
The evaluation of the geometry-based CFs involves the performance of mechanical simulations where the LV was passively inflated to end-diastolic pressure applied on the endocardial surface of the LV mesh. The motion of the basal plane nodes was prescribed from the data, which for the case of the synthetic data set translates to maintaining a fully fixed basal plane. All boundary conditions (BCs) were applied in 30 equal increments. Figure 2 schematically shows where BCs are applied and how they are determined from the clinical data.
The finite elasticity problem was solved within a multifield variational principle approach, with incompressibility enforced through a Lagrange multiplier. Cubic and linear Lagrange interpolation were chosen for the displacement field and pressure, respectively (Hadjicharalambous et al. 2014). The mechanical simulations were performed in the CHeart 3 nonlinear FE solver following a Galerkin FE method (Lee et al. 2016).

Methodology to assess CF performance
To assess the parameter identifiability provided by the geometric and energy-based CFs, we computed the CF residual The evaluation of the energy-based CF is entirely data based, while the evaluation of the geometry-based CF requires the performance of mechanical simulations with sweeps over C 1 and α. The combination of the CFs ensures the unique estimation of these parameters across the α-C 1 parameter space. We then visualize the landscapes of the CF residuals, and locate the parameter subspaces that could potentially be identified as solutions to the inverse problem. This, for the case of the geometry-based CFs, requires the performance of mechanical simulations with parameter sweeps over C 1 and α. Conversely, the energy-based CF relies only on data analysis, as highlighted below.

Geometry-based CFs
We evaluate the ability of six geometrically based CFs to uniquely constrain the passive material parameters either independently or in combination. The geometry-based CFs are based on popular CFs from the literature, comparing displacement, strain and cavity volume Gao et al. 2015;Hadjicharalambous et al. 2015;Mojsejenko et al. 2015;Xi et al. 2013), and extended to include widely used clinical indices, such as wall thickening and apicobasal deformation. The required data for all these CFs can be readily provided from available imaging modalities. For the CF evaluation we consider the time at end diastole (ED), where the largest amount of deformation can be observed in the diastolic window. Specifically, the examined geometry-based CFs are: • L 2 displacement norm. The L 2 displacement norm CF is estimated by comparing the displacements between the simulated displacement (u sim ) and the clinically measured or synthetic data (u dat ): • L 2 strain norm. The L 2 norm of the difference in Green-Lagrange strains between simulated (E sim ) and synthetic or clinical data (E dat ): The L2 displacement and L2 strain norm CFs were estimated using 4 Gauss points per element direction. Increasing the Gauss points to 5 per element direction led to a maximum 5 10 −8 mm error for | u| and 2 10 −7 error for | E|, which is well within the expected magnitude of error due to data noise. • Cavity Volume. The Cavity Volume CF (| V |) describes the absolute difference between the LV cavity volumes in the clinical or synthetic data (V dat ) and model simulations (V sim ): • Wall Thickness. The wall thickness metric | d W T | compares the average wall thickness at the equatorial nodes between the simulation (d sim W T ) and the data (d dat W T ): where n 1 , . . . , n n are the node pairs at the equator used for the wall thickness measurements.
apicobasal distance metrics estimate the average difference between the distance of the endocardial (or epicardial) basal nodes to the endocardial (or epicardial) apical node at the data d dat ABendo (or d dat ABepi ) and simulation d sim ABendo (or d sim ABepi ) meshes: where m 1 , . . . , m n are the basal nodes whose distance from the apex is calculated.

Energy-based CF
Based on the modelling assumptions described in Sect. 2.2, the energy conservation dictates the equality of the external work W ext (the work performed by the external loads acting on the tissue) to the internal energy W int (the work of the internal stresses and strains), giving: The internal energy for hyperelastic materials can be expressed via the strain energy function Ψ , which with the chosen constitutive law (Eqs. 1 and 3) yields the internal energy expression in Eq. 14 as a function of the Green-Lagrange strain tensor E.
The external work is estimated as the increase in LV cavity volume (V ) from the reference configuration (V 0 ) to a given volume V D caused by the LV pressure p.
Defining two points in diastole with corresponding volume, Green-Lagrange strain, external work and internal energy of V 1 , E 1 , W ext1 and W int1 and V 2 , E 2 , W ext2 and W int2 , respectively, we can write : The ratios of the external work and internal energy at these two points must be equal and the difference of these two ratios should tend to zero. This provides the energy-based CF, f : Substituting in the definition of W int (Eq. 14) and W ext (Eq. 15) then gives : We can see in Eq. 18 that the constant C 1 of the constitutive Eq. 14 is cancelled out from the numerator and denominator of the right part of Eq. 17. Note that here the strain field is directly derived from the deformation field extracted from the medical images, without any forward simulation involving a choice of C 1 , giving a CF dependent only on the material parameters in Q.
In implementing the energy CF we select two time points. In the clinical study these correspond to diastolic frames (D F) of the MRI sequence. To avoid potential artefacts from slow decaying active tension we choose to use frames from the MRI that corresponded to the last two frames of end diastole (this choice is reviewed in Appendix 7). We define D F 2 as the end-diastolic frame, and D F 1 as the frame prior to D F 2 (see also Fig. 2). For consistency, we also chose to use the last two 'frames' in the analysis of the synthetic dataset. These correspond to the last two increments of the simulation used to generate it (D F 1 corresponds to the solution after inflation to 1.45 kPa and D F 2 to 1.5 kPa ).
The external work for each time point is calculated by integrating the product of the pressure and change in volume between sequential MRI frames, giving: In Eq. 19, R F corresponds to the index of the cine sequence that corresponds to the reference frame, M R I F corresponds to the index of the MRI frame of interest (for example the index for D F 1 or D F 2 ), p n is the pressure at MRI frame n and V n is the LV cavity volume at MRI frame n.
The internal energy W int is calculated solely from the Green-Lagrange strain field which is derived from the displacement field between the geometries of the D F under consideration and the MRI frame employed as the reference frame. This tensor field can be calculated directly from the image registration algorithm without any further requirement for mechanical simulations of the LV model.
The energy-based CF is only dependent on the parameters in Q in the Guccione law. Assuming constant anisotropy ratios then allows the α parameter (Eq. 3) to be uniquely inferred form the energy-based CF.

Proposed parameter estimation workflow
The proposed workflow relies on the combination of the energy-based and the L 2 displacement norm (| u|) CFs, inferring the unique C 1 − α parameter set from the point of intersection between the minimum residual contours of the two CFs (Fig. 2). Thus, the steps are: -Step 1. Estimate α through minimization of the energybased CF from analysing the data. -Step 2. Perform mechanical simulations in order to optimize C 1 from the | u| CF.
The L 2 displacement norm is chosen as the geometric CF (Sect. 2.3.2) for the C 1 parameter estimation due to its robustness (see Sect. 3.2.1) and the comprehensive data model deformation comparison it provides compared to simpler metrics. Note that the optimization of C 1 can be achieved by setting α to the value obtained in step 1 and sweeping over C 1 , resulting in 1D searches of the C 1 value that minimizes | u|. However, C 1 estimation through 1D optimization may not always be possible, as for certain C 1 , α combinations the nonlinear mechanical solver may not converge (Land et al. 2015b). To overcome this, simulations with 2D sweeps over both C 1 and α can be performed in order to allow an exponential fit of the parameter combinations yielding the minimum | u| residual (for a justification on the choice of the exponential fit, see Xi et al. (2013)). C 1 can then be uniquely estimated from the intersection of this curve with the flat line corresponding to the α solution from step 1 (see Fig. 5c for an example). The latter approach was followed in our study for parameter estimation from both synthetic and clinical datasets.

Parameter estimation in the synthetic data set
To determine if geometric CFs or the energy-based CF can constrain passive stiffness parameters, in the absence of data noise and under conditions of absolute model fidelity to the data, we evaluate the CF performance on a synthetic data set with baseline Guccione constitutive law parameters set to α = 30, C 1 = 1000 Pa, r f = 0.55, r f t = 0.25, r t = 0.2 (Eqs. 1, 3).

Identifiability of the geometry-based CFs
The reported coupling between the C 1 -α parameters (Xi et al. 2013) is confirmed for the L 2 displacement norm (Fig. 4a) and extended for the remaining geometry-based CFs (Fig. 4b-f). Fitting an inverse exponential function to the parameters with the minimum residual for each CF reveals that the CF minimization contours are highly coincident (Fig. 6). This shows that for the in silico case the geometrybased CFs, independently or in combination, are unable to uniquely constrain the parameters of the Guccione law.

Identifiability of the energy-based CF
The landscape of the energy-based CF residual in the C 1 -α parameter space is shown in Fig. 3. Due to its formulation the energy-based CF is independent of the C 1 parameter, as is evident by the fact that its minimization contour is parallel to the C 1 axis and the minimum occurs for a unique value of α. Combining the energy-based CF with the L 2 displacement norm the ground truth C 1 , α parameters of the synthetic dataset were recovered (Fig. 3).

Parameter estimation in the clinical data sets
Following the in silico analysis we investigated the CF performance in 8 clinical cases.

Evaluating Geometric CFs on Clinical Data
The energy-based CF must be paired with a geometric CF to constrain both the C 1 and α parameters. To determine the geometric CF to pair with the energy-based CF we evaluated the six proposed geometric CFs on the 8 clinical data sets. The identifiable parameter combinations for each CF for each clinical data set are presented in Fig. 5 as summary plots of the exponential fits to the CF residual minimization parameter contours. This figure confirms that the C 1 -α parameter coupling exists in vivo for all the geometric CFs. However, the minimization contours are not always coincident in the clinical setting, with some of the CF producing discordant parameter solutions as in cases PC2 and PC7 in Fig. 5.
The L 2 norm of displacements was selected as the geometric CF to pair with the energy CF, as it is based on a thorough global comparison of the agreement of the deformation field between model and data and consistently accorded well with the majority of the other geometric CFs across cases.
. The parameter grid used is shown as the empty black circles and is in the range of 200-5000 Pa for C 1 and 5-300 for α. The parameter combinations yielding the minimum CF residual are plotted in blue. The white patches in the plots a-f signify parameter combinations that resulted in simulations that could not solve with the defined loading paradigm

Identifiability of the energy-based CF
The energy-based CF was estimated for the 8 clinical data sets. Its independence on C 1 is verified in clinical data, as the CF minimizing parameter combinations form a horizontal line parallel to the C 1 axis (see Fig. 5c, where the fitted line to the minimum residual contour is overlain on top of the exponential fits to the geometry-based CF minimums).

Estimated Parameters from the proposed pipeline
Following the proposed pipeline (Sect. 2.4), the passive material parameters for the 8 clinical data sets were determined from the intersection of the fits to the minimum residuals of the energy and L 2 displacement norm CFs (Fig. 5). The identified parameters are shown in Table 2 along with the L 2 displacement norm residual for each case. A certain level of variability is evident in the estimated parameters. In HF patient cases, C 1 ranges from 820 Pa for case PC2 to 5300 Pa for PC1, and α from 5 for PC4 to 66 for PC3 and PC7. The healthy volunteer data set yielded 1700 Pa for C 1 and 15 for α. To provide context for the fitted parameters in the parameter estimation results of this study, previous estimates of the Guccione parameters fitted to human data are presented in Table 3. Fig. 4 Lines of minimal residual for all the geometry-based CFs in the synthetic dataset, after an exponential fitting (each line is an exponential fitting to the points of minimum residual, see Fig. 5c for an example). The result shows that the lines of minimal residual are nearly identical for all geometry-based CFs in silico

Comparison with Previous Methods
We test if the proposed energy CF method for unique parameter estimation predicts significantly different local stresses as compared to one of the previous approaches, specifically where C 1 was fixed to 2000Pa (Xi et al. 2013). The difference in stresses developed at the ED frame using the two methods is presented in Table 4. The large discrepancies observed in  In case PC5 there is no available residual as the forward simulation with the identified parameters did not converge. The α parameter is estimated by the energy-based CF and the C 1 parameter from the | u| CF  Stress values are the deviatoric second Piola-Kirchhoff stress in the fibre direction (Ŝ f f ) computed at end diastole. Estimated α from fixing C 1 (α B ) are also reported for completeness some cases, notably PC2, can be explained by the existence of large local strains that amplify differences due to the exponential term in the strain energy function. In addition, Fig. 7 illustrates the difference in stress-strain curves corresponding to a 1-D fibre stretch for the parameter pairs estimated with both methods.

Discussion
We have shown that unique identification of myocardial material parameters is possible with a suitable choice of the CF. To our knowledge, this is the first time that the two coupled parameters in the Guccione model have been uniquely constrained by clinical data; previously this issue was addressed by fixing part of the parameter set Hadjicharalambous et al. 2015;Wang et al. 2013;Xi et al. 2013).

Identifiability by an energy-based CF
The core methodological contribution of this work is the proposal of a CF that removes the parameter coupling. The energy-based CF identifies α due to its independence to the C 1 parameter. Its accuracy was tested in silico, where it estimates the correct α parameter and combined with the L 2 displacement norm provides accurate estimates of the ground truth parameter values. Results in 8 real clinical data sets with the complete pipeline demonstrate that the CF is robust to the inherent noise in clinical data and finite model fidelity. The novel energy-based CF is also a data driven metric. Only the data of deformation (strain and cavity volume) and pressure are required to compute it, without the need of computational simulations or data assimilation pipelines. This has three main benefits. Firstly, the data derived deformation field employed in the CF is unaffected by the C 1 -α coupling that arises from the simulation. Secondly, the computationally expensive search over the full parameter space involved in current data assimilation schemes has one dimension of the parameter space reduced since the α parameter is fixed. Thirdly, the reduction of methodological complexity to obtain the α parameter opens the possibility for a quicker and easier clinical adoption.
It is important to note that the energy-based CF raises the demands on data quality and quantity, since it requires strain data of the entire myocardium at two time points during diastole and the pressure-volume information covering the filling phase of the cycle. The importance of data quality on parameter estimation is demonstrated in a sensitivity study, provided in Appendix 6. In the absence of high fidelity strain data it is possible to recast the energy CF in terms of a pressure CF. This allows unique parameter estimates from pressure and volume transient data alone. The efficacy of this approach is presented in Appendix 8.
In our calculations the external work is estimated using a pressure-volume approach (see Eq. 19) which is fully accurate for the case of a deformation field consistent to the passive inflation assumption we have adopted. However, the image driven Dirichlet boundary conditions applied on the basal plane in the clinical data sets contribute to external work. This contribution is quantified as a mean 5% of the elastic energy in the clinical cases based on forward simulations with the identified parameters.
The efficiency of the energy-based CF was demonstrated for the Guccione material law, but can be extended to other exponential constitutive relations for reducing the mate-rial parameter redundancy by one, such as the Holzapfel-Ogden law (Holzapfel and Ogden 2009) as demonstrated in Appendix 9 or the pole-zero (Nash and Hunter 2000) among others.

Geometry-based CFs
Geometry-based CFs, and their combination, were not able to identify unique myocardial material parameters, agreeing with previous reports (Augenstein et al. 2006;Xi et al. 2011).
We investigated if a combination of geometrical CFs could improve parameter identifiability. The C 1 -α parameters would then be identified by an intersection of the minimization contour of two or more CFs. Nevertheless, in the in silico data set the minimization contours are almost identical for the different CFs, suggesting the low complementary value of the CFs (Fig. 4). On the contrary, results with real data report a large variability in the agreement between the different CFs in half of the cases (see the offsets between lines that identify the coupling in Fig. 6), suggesting that this strategy does not lead to a unique set of parameters in practice. The disagreement between CFs with real data, and not with simulated data, is interpreted as a reflection of the mismatch between model and real data, caused by a combination of lack of model fidelity and data quality.
Note that while the geometric cost functions are based on a single frame, in contrast to the two frames used in the energy CF, the addition of an adjacent frame is not expected to improve the identifiability of parameters from geometric cost functions when working with clinical data due to the presence of noise that is sufficiently large to obscure the global minimum, as reported in Xi et al. (2013).

Parameter estimation workflow
The proposed parameter estimation pipeline lead to a unique estimation of the Guccione material parameters in the 8 clinical data sets analysed in this study, and with | u| residuals (Table 2) comparable to previously reported errors (Table 3). It is important to note that the set of kinematic BCs in this study was lighter (only constraining the base, and not also the apex as in Xi et al. (2013)), thus making the task of reproducing the clinical observation more challenging.
In the proposed methodology the unique C 1 and α parameters, where the coupling occurs (Xi et al. 2013), are found while fixing two of the less intercorrelated ratios r f , r t , r f t (the third is bound by Eq. 6). Once C 1 and α are found, the ratios can be uniquely found (as reported in Xi et al. (2013)). The impact of a wrong initial choice of r f , r t , r f t on the estimation of C 1 and α parameters was evaluated in the sensitivity study (Appendix 6) and was found to be relatively low. An important remark in the methodology is the existence of challenges associated with the convergence of nonlinear mechanics solvers in incompressible applications (Land et al. 2015b). The lack of convergence can often obstruct the calculation of the geometry CF residual for the whole C 1 range of interest under the known α. Knowing that the set of coupled parameters that lead to very similar minimum costs draw a line in the α-C 1 log-scale space (Xi et al. 2013) allows the problem to be reformulated as finding the parameters of this exponential line.

Model assumptions
There are a series of model assumptions which, although not affecting the contribution of the proposed pipeline to parameter identifiability, may have an impact on the parameters found.
One important factor determining diastolic filling is the residual active tension (AT), which is known to be present in diastole (Bermejo et al. 2013;Xi et al. 2013). In this study the parameters are estimated from late diastolic instants and the inclusion of an earlier frame suggested the presence of remaining AT as detailed in Appendix 7. End-diastolic events, where the AT can be assumed to be limited and its contribution to the work estimation negligible, are the most suitable observations. Following this approach, any remaining AT at end diastole leads to an apparent increased myocardial stiffness , specifically in the fibre direction (Xi et al. 2013).
The most important element in the proposed methodology was revealed to be the choice of the reference frame (see Appendix 6), in concordance with previous studies (Xi et al. 2013). The reference geometry directly affects the observed myocardial stiffness as it dictates the measured strain under a given cavity pressure. In this study the LV geometry at the lowest pressure frame was chosen to describe this geometry following a popular approach Gao et al. 2015;Land et al. 2012;Nikou et al. 2016;Wang et al. 2009) in order to simplify the workflow. Inclusion of a reference frame estimation (Krishnamurthy et al. 2013) can possibly enhance the parameter estimation in future applications.
Assumptions are also made regarding the definition of the myocardial microstructure and material. Conforming to the majority of FE studies in the field of cardiac mechanics, the myocardium is assumed to be incompressible although capillary and coronary flow are known to locally violate this assumption (Ashikaga et al. 2008;Yin et al. 1996). This assumption affects both the simulated deformation fields from the mechanics solver, as well as the novel energy-based CF, where deformation due to cardiac perfusion (increment of volume during diastole) is assimilated to contribute to the tissue strain energy. Also in the definition of myocardial microstructure, the inclusion of a more realistic fibre field might improve the accuracy in the estimation of the projected strain components (Eq. 3) and thus that of the energy-based CF (Eq. 18). Nevertheless, results in the sensitivity study (Appendix 6) suggest that the impact of this assumption is very small, in accordance to previous studies (Land et al. 2015a).
A central assumption in the modelling approach followed here is myocardial material homogeneity, which by reducing model complexity facilitates the parameter estimation procedure. However, this currently restricts the application of the proposed method to disease where this assumption is valid. Although rendering our method suitable for cardiac disease with localized stiffness alterations -such as myocardial infarction-is within our future plans, our workflow is readily suitable for applications on disease, such as dilated cardiomyopathy, diffuse fibrosis or heart failure with normal ejection fraction (HFnEF), where tissue properties are expected to be more homogeneous.
One last set of assumptions are needed to define the BCs of the model. First, a homogeneous pressure load is assumed to act on the endocardial boundary during ventricular filling. This is a reasonable simplification based on the reported cavity pressure variations in the literature (de Vecchi et al. 2014), and is commonly taken for computational efficiency (Gao et al. 2015;Hadjicharalambous et al. 2015;Mojsejenko et al. 2015;Nikou et al. 2016;Wang et al. 2013). However the impact of the RV, atria and pericardium on achieving more physiological deformations is known (Belenkie et al. 2001;Tyberg and Smith 1990;Williams and Frenneaux 2006), and thus the use of more advanced mechanical models (Fritz et al. 2014) is anticipated to improve model fidelity and therefore parameter estimation accuracy. We would also expect that more realistic BCs will enable the model to better reproduce the recorded myocardial deformation, thus leading to smaller residual | u|-and this should be especially beneficial in 3 of our cases (PC3, PC4 & HC, see Table 2).

Estimated parameters in vivo
Few studies have reported passive material properties for multiple patients. The identified parameters for the 8 cases in this study (Table 2) are within the reported range in the literature (Table 3). In our results, C 1 falls within a range of 800-5300 Pa for human HF patients, as opposed to higher values obtained when α was held fixed. The α values fitted here span from 5 to 66 and lie within the literature range, while the higher values corresponding to HF patients are significantly lower to the ones previously reported when fixing C 1 at 2000 Pa (Table 3).
In our results we found that in four out of seven patient data sets (PC1, PC3, PC5, PC7) both identified parameters were higher than those of the healthy data set, in two cases (PC2, PC6) only the α parameter was increased while in one case (PC4) C 1 was increased instead. The sample of 7 HF subjects already reveals functional differences by the linear and exponential components of the material law, and we can speculate that this may have links with the aetiology of the disease. While the general trend is for the HF patient data sets processed to have higher stiffness than the healthy case in agreement with previous clinical studies , further cases need to be investigated to confirm this observation.

Cardiac Mechanics Application
Computational models are increasingly used in clinical applications. Inferring material properties from these models offers three potential applications with each benefiting from uniquely constrained parameters. Firstly, material parameters may provide a more sensitive descriptor of patient pathology . To test the utility of this application first requires a method for inferring constitutive parameters from clinical data. Assuming that these parameters reflect some underlying material property, and hence reflect the patients pathology, then the inferred parameter values should be unique and independent of the fitting method.
Secondly, the use of biophysical models, that are constrained by physical laws, increases the capacity of models to predict outside of the data used to constrain them. This is particularly important for patient-specific models that could be used to predict response to treatments from pre procedure data. The use of non-unique parameters will make these predictions dependent on the parameter inference method, increasing the uncertainty in the model predictions.
Finally, models can be used for predicting values of interest that can not easily be measured, including material stress, regional work and local mechanical efficiency. We have shown that an alternative approach for a unique estimation of parameters, by fixing C 1 , leads to significant discrepancy in model predictions (see Fig. 7).

Conclusions
A novel and clinically tractable pipeline for passive myocardial material estimation is proposed, which manages to decouple the material constitutive law parameters and guarantee reliable material estimation. This is an important step towards the use of myocardial stiffness as a reliable tool for the understanding of cardiac pathophysiology and the development of biomechanically relevant biomarkers. This work highlights the central role of CFs in the identifiability of material parameters.

Sensitivity analysis
This section tests the performance of the novel energybased CF in the presence of data quality issues or arbitrary modelling assumptions that are necessarily made in clinical models. For this purpose we use the in silico workbench with access to ground truth parameter values, and then perform a sensitivity analysis. We focus on the impact in the estimation of the α parameter (ground truth is 30), since it is provided by the novel energy-based CF.

Sensitivity to pressure offset errors
Artefacts in the pressure data were studied introducing an offset, which can be caused by calibration errors and are hypothesized to be the main source of material parameter errors (see the large impact reported in Xi et al. (2014)). The effect of offset pressure values is investigated by shifting all pressure data points by +/−10% of the mean pressurep (corresponding to a 75 Pa increase/decrease in the pressure values).

Sensitivity to strain data quality
The quality of the strain field in the imaging data can affect the accuracy of the estimated α, as the energy-based CF relies on the integration of strains in two diastolic timepoints (Eq. 18).
To assess the anticipated error in parameter estimation from strain data noise, we inserted white noise to the strain tensor at the gauss points used for the integration of strain energy density function in Eq. 14. The noise was independently applied to each strain component in fibre coordinates, and corresponded to a standard deviation of 5, 10 and 20% of the mean strainĒ i j , which varied from −0.15 to 0.10 in the employed diastolic frames (corresponding to increments 29 and 30 of the original simulation).

Sensitivity to modelling assumptions
Two modelling assumptions are challenged, the choice of the reference configuration and the mathematical description of the myocardial microarchitecture (fibre field). The reference frame, which represents the geometrical configuration at which the myocardium is unloaded and the developing stresses and strains are zero, is difficult to calculate. In the clinical application of the parameter estimation pipeline, the reference frame is estimated by taking an observed geometry at the beginning of diastole, similar to previous approaches Mojsejenko et al. 2015;Wang et al. 2013). To estimate the effect of an inaccurate reference configuration on the parameter solu- tion, later frames (frames 1 to 5, where frame 0 is the original reference frame) were used as reference configurations in the mechanical model. The cardiac microarchitecture is approximated by rulebased mathematical models (Bayer et al. 2012;Nielsen et al. 1991) with −60 • to +60 • angle variation. To assess the effect of this choice of fibre microarchitecture, a set of extreme vertical fibres was selected (fibre angles varying from −90 • to +90 • ).

Effect of fixed parameter values
Throughout the parameter estimation procedure fixed values for the exponential stiffness parameter ratios were assumed. To evaluate the effect of this arbitrary choice, two different stiffness ratio combinations were used: one pertaining to a highly anisotropic material (r f = 0.85, r f t = 0.10, r t = 0.05) and one to an isotropic material (r f = 0.34, r f t = 0.33, r t = 0.33).

Results and discussion
Results (see Table 5) demonstrate a satisfactory robustness of the approach to possible sources of disagreement between model and data. The largest effect on parameter estimation accuracy derives from an unsuitable choice of reference configuration (choice of frame 5, α sol = 3). The offset in pressure recordings has the second largest impact on the parameter estimation, with smaller parameters found with larger pressure values and vice-versa. The anticipated effect of the quality of the strain field on the estimated parameter values is also confirmed. The microstructural representation in the model seemed to have a negligible effect on the identified values. The stiffness anisotropy ratios did not bias the estimation of α (despite the pronounced anisotropy vs isotropy scenarios tested) but had a larger effect on the estimation of C 1 , which was 2100 Pa for the highly anisotropic assumption and 620 Pa for the isotropic one. In both cases the resulting | u| residual of the forward simulations with the identified parameters was increased from 3 10 −9 mm for the ground truth parameters to 1.15 and 1.07 mm, respectively.

The effect of using additional frames in the energy-based CF on α
In our analysis the two last frames of diastole were used as they provide larger strain signals with lowest probability for presence of residual AT. Alternatively, additional late diastolic frames can be included in the analysis, leading to a modified functional using the ratios of the paired combinations of the last 3 diastolic frames f 3DF (Eq. 20): Results in Table 6 show that the change in α is generally small (up to a 25% except for a 40% in PC6). In most of the cases the additional earlier frame led to a larger value in α, suggesting the presence of residual AT. The data of the additional earlier diastolic frame in PC5 indicated a decrease in strain energy, despite pressure increasing, which is incompatible with underlying model assumptions and results in a negative α value.

An alternative approach for unique parameter estimation
In this work a novel CF was proposed that can isolate α, allowing for C 1 to be uniquely determined from previous geometrical CF. The proposed energy CF relies on strain fields derived directly from the images, and thus on the availability of high quality data and robust image registration tools (Shi et al. 2012). A similar approach can be followed by exploiting Note that C 1 here is obtained in both cases by the L 2 displacement norm a Results in PC3 and PC7 are bounded since the passive inflation simulation for larger α values did not converge the fact that the cavity pressure that inflates the myocardium to a given volume scales linearly with the parameter C 1 . This leads to the definition of an alternative pressure based CF, which relies on the good quality of the pressure data and their synchronization to the volume trace: where p dat1 , p dat2 are the measured cavity pressures at two timepoints in diastole and p sim1 , p sim2 are the esti at the two last diastolic frames (DF2, DF1 are again the end-diastolic and pre-end-diastolic frames from the MRI sequence, respectively) and p sim DF 1 , p sim DF 2 are the estimated cavity pressures at the respective timepoints, which are required to inflate the reference mesh to the measured cavity volumes. For consistency with the energy CF analysis, these diastolic timepoints are defined as frames DF 1 and DF 2 corresponding to the pre-end-diastolic and end-diastolic frames of the MRI sequence (Sect. 2.3.3). As p sim scales with C 1 , the ratio p sim1 p sim2 is independent of C 1 and therefore f p depends only on α (for fixed anisotropy ratios). Note that the pressure CF can only be evaluated through forward simulations, while in contrast the energy CF can be evaluated directly from the clinical data.
To evaluate the pressure form of the cost function we find the α parameter for the 8 cases of this study. For this purpose, parameter sweeps over α for an arbitrary fixed C 1 at 2000 Pa were performed. For each α the myocardium was inflated to an elevated cavity pressure (3 times the end-diastolic pressure) in 90 increments. The values p sim1 and p sim2 were found by linear interpolation between the two simulation increments with volumes that bound the cavity volume from the data. Following estimation of α from Eq. 21, C 1 was estimated from the already fitted exponential curves to the minimum L 2 displacements norm CF residual (as explained in Sect. 2.4) and the results are presented in Table 7.
Despite the fact that both the energy CF and the pressure CF do find the correct unique set of parameters in an in-silico test, parameters show large differences depending on the CF when working with real data (see Table 7). This is a similar finding as when working with the collection of geometrical CFs, and reflects the incongruity between data and model accentuated by the different origins and types of data errors.
The assumption that these patients, who share a common pathology, have similar material properties would suggest that for the available data the energy-based CF provides a more robust estimate of model parameters. Nevertheless, there is no ground truth methodology available for the evaluation of the performance of the different CFs. Further research in data acquisition and analysis, on model fidelity, and on alternative in-vivo strategies may improve the agreement between the different CFs and lead to the true material parameters of the beating myocardium.