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 [Formula: see text] and [Formula: see text] ([Formula: see text], [Formula: see text], [Formula: see text]) for the HF cases and [Formula: see text] and [Formula: see text] ([Formula: see text], [Formula: see text], [Formula: see text]) 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 pressure–volume (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. 2008, 2009; Nair et al. 2007) and has led to the estimation of material parameters from clinical data sets (Asner et al. 2015; 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).

Data sets

Synthetic

To provide a known ground truth for the material parameters a synthetic data set was employed (see top panel in 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.

Fig. 1

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 registration 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)

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).

Table 1

Summary of the patient cases (PC1-PC7) and healthy data set (HC) used

Case	Age	Sex	EF (%)	ESV (ml)	EDV (ml)	EDP (kPa)
PC1	61	M	13.5	266	307	2.59
PC2	61	M	6.2	348	371	1.21
PC3	70	M	19.5	174	216	4.44
PC4	76	F	32.3	86	127	1.84
PC5	57	F	19.3	214	265	2.99
PC6	65	M	29.7	122	173	1.17
PC7	39	M	19.7	176	219	2.98
HC	36	M	-^a	-^a	134	1.89

Data not available

The abbreviations used are as follows: EF ejection fraction, ESV end-systolic volume, EDV end-diastolic volume all corresponding to the LV

Volumes were estimated from the personalized meshes following the p–V synchronization

The cardiac images of the CRT patient data sets (PC1–PC7) consist of 2-D short axis stacks of cine MRI with SENSE encoding ( to 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 non-rigid registration (Shi et al. 2013) which enables a spatially and temporally continuous description of the cardiac motion. Mesh personalization was performed on the segmented LV, as described previously (Lamata et al. 2014). 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 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 pressure–volume 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 (RV), atria, pericardium and other neighbouring structures have secondary roles.

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 at the epicardium to 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 , , assign different mechanical responses to the tissue along the fibre (f) direction, across the transverse planes (t) and in the fibre-transverse shear planes (ft), respectively. The f,s,n indices in the strain components express projections of the Green–Lagrange strain tensor () 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 () and bulk exponential () parameters, along which the primary coupling occurs (Eqs. 1, 3). The parameters are referred to as anisotropy ratios and range between 0 and 1, with the 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 while focusing on fitting the coupled and parameters.

Reference configuration

The reference configuration represents an idealistic stress- and 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.

Fig. 2

Overview of the proposed parameter estimation pipeline. The required input consists of LV meshes and corresponding pressure values at the defined diastolic window of interest (covering the frames from minimum pressure until before the beginning of contraction, see also Fig. 1). 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 and . The combination of the CFs ensures the unique estimation of these parameters

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).

Examined CFs and their evaluation

Methodology to assess CF performance

To assess the parameter identifiability provided by the geometric and energy-based CFs, we computed the CF residual across the - 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 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 (Asner et al. 2015; 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:

displacement norm. The displacement norm CF is estimated by comparing the displacements between the simulated displacement () and the clinically measured or synthetic data ():

strain norm. The norm of the difference in Green–Lagrange strains between simulated () and synthetic or clinical data (): 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 mm error for and 2 10 error for , which is well within the expected magnitude of error due to data noise.

Cavity Volume. The Cavity Volume CF () describes the absolute difference between the LV cavity volumes in the clinical or synthetic data () and model simulations ():

Wall Thickness. The wall thickness metric compares the average wall thickness at the equatorial nodes between the simulation () and the data (): where are the node pairs at the equator used for the wall thickness measurements.

Apicobasal distance. The endocardial (or epicardial ) 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 (or ) and simulation (or ) meshes: where 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 (the work performed by the external loads acting on the tissue) to the internal energy (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.  as a function of the Green–Lagrange strain tensor .The external work is estimated as the increase in LV cavity volume () from the reference configuration () to a given volume caused by the LV pressure p.Defining two points in diastole with corresponding volume, Green–Lagrange strain, external work and internal energy of , , and and , , and , 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 (Eq. 14) and (Eq. 15) then gives :We can see in Eq. 18 that the constant 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 , 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 (DF) 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 as the end-diastolic frame, and as the frame prior to (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 ( corresponds to the solution after inflation to 1.45 kPa and 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, RF corresponds to the index of the cine sequence that corresponds to the reference frame, corresponds to the index of the MRI frame of interest (for example the index for or ), is the pressure at MRI frame n and is the LV cavity volume at MRI frame n.

The internal energy is calculated solely from the Green–Lagrange strain field which is derived from the displacement field between the geometries of the DF 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 displacement norm () CFs, inferring the unique parameter set from the point of intersection between the minimum residual contours of the two CFs (Fig. 2). Thus, the steps are:The displacement norm is chosen as the geometric CF (Sect. 2.3.2) for the 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 can be achieved by setting to the value obtained in step 1 and sweeping over , resulting in 1D searches of the value that minimizes . However, estimation through 1D optimization may not always be possible, as for certain , combinations the nonlinear mechanical solver may not converge (Land et al. 2015b). To overcome this, simulations with 2D sweeps over both and can be performed in order to allow an exponential fit of the parameter combinations yielding the minimum residual (for a justification on the choice of the exponential fit, see Xi et al. (2013)). 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.

Fig. 5

Proposed CF combination for unique parameter estimation emerging from the in silico analysis. a The landscape of the CF (residual in mm), which is the chosen geometric CF (see Sect. 2.4 for more details). b Landscape of the energy-based CF. A straight line pattern emerges that demonstrates the independence of the CF on and its ability to uniquely identify . c Unique parameter estimation with the combined use of the energy-based CF for parameter identification and the displacement norm CF. A superimposed black curve is fitted to the parameter combinations corresponding to the minimum CF residual values. The horizontal black line corresponds to the minimum energy-based CF residual. The identified parameters (indicated by the black circle at the intersection of the energy-based and geometry-based CF minimization lines) coincide with the ground truth values (, ). (Minimum and maximum CF residuals are shown in blue and red, respectively.)

Step 1. Estimate through minimization of the energy-based CF from analysing the data.

Step 2. Perform mechanical simulations in order to optimize from the CF.

Results

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 , , , , (Eqs. 1, 3).

Identifiability of the geometry-based CFs

The reported coupling between the - parameters (Xi et al. 2013) is confirmed for the 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 geometry-based CFs, independently or in combination, are unable to uniquely constrain the parameters of the Guccione law.

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

Fig. 6

Lines of minimal residual from the geometry-based CF superimposed to the solution for the parameter resulting from the energy-based CF (flat line) in each of the 8 cases studied: a PC1, b PC2, c PC3, d PC4, e PC5, f PC6, g PC7, h HC. i Legend

Identifiability of the energy-based CF

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

Fig. 3

Plots of landscapes of the examined geometry-based CF residuals over the and parameter space for the synthetic data set: a CF (residual in mm). b CF. c CF (residual in ml). d CF (residual in mm). e CF (residual in mm). f CF (residual in mm). The parameter grid used is shown as the empty black circles and is in the range of 200–5000 Pa for 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

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 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 - 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 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 energy-based CF was estimated for the 8 clinical data sets. Its independence on is verified in clinical data, as the CF minimizing parameter combinations form a horizontal line parallel to the 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 displacement norm CFs (Fig. 5). The identified parameters are shown in Table 2 along with the displacement norm residual for each case.

Table 2

Parameter estimation results from the application of the proposed pipeline to the clinical data sets

Case	α_sol	C_1 (Pa)	|Δu| (mm)
PC1	61	5300	0.95
PC2	61	820	1.96
PC3	66	1960	2.56
PC4	5	4780	5.61
PC5	24	3140	-^b
PC6	29	1460	1.91
PC7	66	3300	1.18
HC	15	1700	3.35

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 parameter from the CF

A certain level of variability is evident in the estimated parameters. In HF patient cases, 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 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.

Table 3

Estimated , parameters for human myocardium from previous studies

Case	C_1 (Pa)	α	|Δu| (mm)
Human
Healthy^a	2000^e	43	1.78
Patient 1^a	2000^e	105	1.58
Patient 2^a	2000^e	95	1.39
Healthy^b	3600±1200^e	38	–
HT^b,c	12000±2600^e	38	–
NI-HF^b,d	11800±3400^e	38	–

Xi et al. (2013). Wang et al. (2013). Patients with hypertrophic LV. Patients with non-ischaemic HF with reduced EF. In these studies the scaling constant is defined as half the parameter in Eq. 1, and therefore, values reported here have been doubled for consistency in the results

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 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 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.

Table 4

Differences in stress calculated with parameters estimated by proposed method (A) and a previous one (Xi et al. 2013), where was fixed at 2000 Pa (B)

Case	Mean S^_ff^A-S^_ff^B (Pa)	Standard deviation (Pa)	α^B
Synth	31.7	102.8	17
PC1	-461.8	47103.7	142
PC2	-577767.5	34820092.3	25
PC3	3127.3	52912.1	65
PC4	-7.3	284.9	10
PC5	-0.5	82.2	35
PC6	32.6	114.1	22
PC7	18.7	598	92
HC	45.7	972.7	13

Stress values are the deviatoric second Piola–Kirchhoff stress in the fibre direction () computed at end diastole. Estimated from fixing () are also reported for completeness

Fig. 7

Cauchy stress vs stretch curves for an idealized 1-D extension along the fibre direction (up to 120% stretch) of an incompressible cube with the Guccione law using the material parameters estimated with the proposed methodology (A: energy-based & CFs) and a previous one where one of the parameters was assigned a certain value (B: fixed at 2000 Pa & CF as in Xi et al. (2013))

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 (Asner et al. 2015; 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 parameter. Its accuracy was tested in silico, where it estimates the correct parameter and combined with the 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 - 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 material 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 - 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 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 and parameters, where the coupling occurs (Xi et al. 2013), are found while fixing two of the less intercorrelated ratios , , (the third is bound by Eq. 6). Once and are found, the ratios can be uniquely found (as reported in Xi et al. (2013)). The impact of a wrong initial choice of , , on the estimation of 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 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 - 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 (Asner et al. 2015), 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 (Asner et al. 2015; 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 —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, 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 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) 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 (Wang et al. 2013), 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 (Lamata et al. 2016). 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 , 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.

Table 5

Sensitivity analysis results with respect to the identified parameter

Data/model modification	α_sol
Pressure
+10% p¯ offset	22
−10% p¯ offset	40
E noise
STD 5% E¯_ij	28
STD 10% E¯_ij	26
STD 20% E¯_ij	37
Ref. Fr.
1	26
2	22
5	3
Fibres
−/+ 90 ^o	29
r_f-r_ft-r_t
0.85-0.1-0.05	29
0.34-0.33-0.33	28

The ground truth value for is 30

Table 6

Results of the estimation of with an energy cost function extended to consider the last 3 frames of diastolic filling (), instead of the 2 of the original formulation (referred to here as )

Case	Synth	PC1	PC2	PC3	PC4	PC5	PC6	PC7	HC
α_3fr.	30	78	79	52	1	-^a	41	77	17
α_3fr.-α_2fr.	0	17	18	12	-4	-^a	12	11	2

PC5 had data inconsistent with passive inflation, see text for further details

Table 7

Parameter estimation results from the pressure-based CF (pCF), compared to the energy-based CF (eCF)

Case	α^pCF	C_1^pCF	α^pCF-α^eCF	C_1^pCF-C_1^eCF
Synth	30	1000	0	0
PC1	18	21630	-43	16330
PC2	6	8335	-55	7515
PC3	>120^a	<885^a	>54^a	<-1075^a
PC4	50	240	45	-4540
PC5	60	1070	36	-2070
PC6	130	275	101	-1185
PC7	>400^a	<210^a	>334^a	<-3090^a
HC	50	310	45	-1390

Note that here is obtained in both cases by the L displacement norm

Results in PC3 and PC7 are bounded since the passive inflation simulation for larger values did not converge