Estimation of passive and active properties in the human heart using 3D tagged MRI

Advances in medical imaging and image processing are paving the way for personalised cardiac biomechanical modelling. Models provide the capacity to relate kinematics to dynamics and—through patient-specific modelling—derived material parameters to underlying cardiac muscle pathologies. However, for clinical utility to be achieved, model-based analyses mandate robust model selection and parameterisation. In this paper, we introduce a patient-specific biomechanical model for the left ventricle aiming to balance model fidelity with parameter identifiability. Using non-invasive data and common clinical surrogates, we illustrate unique identifiability of passive and active parameters over the full cardiac cycle. Identifiability and accuracy of the estimates in the presence of controlled noise are verified with a number of in silico datasets. Unique parametrisation is then obtained for three datasets acquired in vivo. The model predictions show good agreement with the data extracted from the images providing a pipeline for personalised biomechanical analysis.


Introduction
Integrated imaging and mathematical modelling hold significant potential for augmenting and improving current clinical practice. The ongoing advancement in computational modelling enables fast and efficient simulation of complex physiological systems (Casoni et al. 2014;Gurev et al. 2015;Lafortune et al. 2012). Biomechanical models provide the distinct advantage in that they enable clinicians and researchers to estimate characteristics of the heart that are otherwise difficult or impossible to measure in vivo, such as tissue stiffness and contractility, local strains and stresses. These model-derived metrics yield additional, potentially more sensitive, indicators of health and disease. Consequently, personalised cardiac modelling has generated significant clinical interest (Ge and Ratcliffe 2009;Sermesant et al. 2012;Yacoub and Terracciano 2011), actively stimulating research into model developments Krishnamurthy et al. 2013;Marchesseau et al. 2013;Neal and Kerckhoffs 2010) that may improve diagnosis, facilitate prognosis and assist in treatment planning.
Advances in medical imaging and image processing are paving the way for the use of cardiac biomechanical modelling in clinical applications. Magnetic resonance imaging (MRI), computed tomography and echocardiography have advanced to the state where anatomy, motion (kinematics) and flow can all be routinely imaged in patients. Various types of cardiac MR sequences, such as high-resolution threedimensional (3D) anatomy at multiple cardiac states (Uribe et al. 2008), cine images of cardiac motion (Uribe et al. 2007;Usman et al. 2013Usman et al. , 2015, 4D phase contrast MRI (Markl et al. 2011) and 3D tagged MRI (Rutz et al. 2008), provide an extensive set of data for modelling. Motion tracking in 3D tagged images (Chandrashekara et al. 2004) in particular enables the extraction of three-dimensional displacements of individual points in the myocardium and can therefore serve as basis for quantification of the mechanical properties of cardiac tissue (Mojsejenko et al. 2014;Xi et al. 2013). Moreover, non-invasive approaches to estimating pressure from imaging data (Buyens et al. 2005;Donati et al. 2015;Ebbers et al. 2001;Tyszka et al. 2000) or peripheral cuff measurements (Brett et al. 2012;Chen et al. 1997) are becoming more advanced and pervasive in clinical assessments.
Leveraging this quantitative data for construction of cardiac mechanics models, it is essential to provide reliable and robust models and parameterisation pipelines for processing patient data to producing results that enable medical interpretation. In this context, model selection is of critical import, requiring a balance between the objectives of the analysis and input extractable from the data. This mandates selection of a model that is inherently a compromise, providing sufficient richness to represent the physiological processes of interest (i.e. sufficient model fidelity), while limiting the number of patient-specific parameters to ensure identifiability. Many modelling approaches have been reported in the literature (Camara et al. 2015), using a variety of laws to describe the passive (Costa et al. 2001;Guccione et al. 1995;Holzapfel and Ogden 2009) and active Kerckhoffs et al. 2003;Niederer et al. 2011) behaviour of the myocardium. Quantification of mechanical properties from imaging data has also been pursued (Augenstein et al. 2005;Chabiniok et al. 2012;Chandrashekara et al. 2004;Gao et al. 2015;Göktepe et al. 2011;Imperiale et al. 2011;Mojsejenko et al. 2014). However, model complexity often limits the ability to identify patient-specific parameter values uniquely (Hadjicharalambous et al. 2014b;Xi et al. 2013), thus adding uncertainty to the quantitative characterisation the model can provide.
The aim of this work is to develop a relatively simple model characterising full-cycle left-ventricular (LV) mechanics with sufficient accuracy, and an associated pipeline for reliable estimation of a limited number of patient-specific parameters, based on a non-invasively acquired dataset including 3D tagged MRI. Successful implementation of this pipeline would provide a basis for linking the modelderived values, such as passive and active parameters, to the underlying pathophysiology. The computational model and the associated parametrisation protocol introduced here are building on our previous work for passive parameter estimation from 3D tagged data (Hadjicharalambous et al. 2014b) and full-cycle parameter estimation from pressure-volume data (Asner et al. 2015). A reduction in the Holzapfel-Ogden material model (Holzapfel and Ogden 2009) is used that was shown to balance identifiability and model fidelity across a range of passive behaviours (Hadjicharalambous et al. 2014b). Novel external boundary conditions are applied, so that deformations are driven by LV volumes and basal longaxis motion extracted from processed image data. Finally, a simplified active constitutive law is integrated into the model, providing an estimate of the active tension in the myocardium over the cardiac cycle.
In passive parameter estimation, optimal values are selected based on the L 2 (Ω 0 )-norm difference between 3D motion fields extracted from the data and model predictions. Active tension through the cardiac cycle is estimated based on a combined functional incorporating relative L 2 -norm displacement errors and cavity pressure errors. Parameter identifiability and model sensitivity are verified in a series of tests based on in silico data emulating acquired data and processing errors. These tests provide a useful tool for evaluating the estimation pipeline under perfect or nearperfect model fidelity and allow us to assess the errors in the estimated quantities. With the view of using the process in patient-specific modelling, we then apply the model and the parameterisation protocol to three in vivo datasets including MR images and pressure data as an important verification step. Model fidelity is reduced for in vivo data, and the feasibility of the estimation process, as well as the reliability of computed parameter values, needs to be re-evaluated for potential use in clinical cases. Parameter identifiability is verified in the in vivo cases, and good agreement between the data and the model can be observed.
Section 2 gives the details of the full-cycle LV mechanics model used in the pipeline and describes the datasets and the process for model personalisation. The results for all of the in silico and the in vivo tests are presented in Sect. 3. The outcomes of the testing, along with limitations and future work, are discussed in Sect. 4.

Methods
The specifics of the proposed pipeline for personalisation of full-cycle LV mechanics models are described below. First, we present the cardiac mechanics model, including the governing equations, the passive and active constitutive laws, and the boundary conditions. A minimal clinical dataset required is outlined next, and the procedure for personalisation of the model based on this dataset is described. We also present the steps involved in generating the in silico datasets.

Governing equations
The equations used to model myocardial deformation over the cardiac cycle are derived from the basic conservation principles (Bonet and Wood 2008). We employ a standard Galerkin finite element formulation of the solid mechanics problem, while choosing specific constitutive laws and boundary conditions to ensure good model fidelity and parameter identifiability based on the available data. The details of the model are described below.
Let Ω 0 ∈ R 3 be the reference configuration of the domain of interest (the unloaded myocardium) with X denoting the corresponding reference coordinate. Let also Ω t ∈ R 3 be the deformed domain with x the corresponding physical coordinate at a given time t > 0. The mapping between the deformed and the reference coordinates is assumed to be diffeomorphic at all times. At any time t > 0, the displacements of material points u = x − X and the hydrostatic pressures p together with the Lagrange multiplier vector λ (see Sect. 2.1.3) are found as the saddle point of the energy potential functional Π (Bonet and Wood 2008): where stable combination of Sobolev spaces for any t ≥ 0 (Brenner and Scott 2008). The solution must coincide with a critical point of the functional: for any (v, q, μ)(t) ∈ X × Y × Λ at a given t > 0. Eq.
(2) provides the weak form equations for solving the solid mechanics problem.

Constitutive relations
The precise form of the energy functional depends on the choice of models that represent the myocardium. The total energy can be split into internal and external energy terms: The former is determined solely by the properties of the material, and the latter comes from the boundaries where external forces are applied. We assume that myocardial tissue behaves like an incompressible hyperelastic material, where the internal strain energy can be represented as the sum of the total passive (W p ) and active (W a ) strain energies of the body and an incompressibility penalty: is the determinant of the deformation gradient tensor for a deformation v: The constitutive laws defining the strain energy functions are usually derived from experimental data on material response to different loading conditions. It is established that the passive stress-strain relationship for the myocardium is highly nonlinear (Nash and Hunter 2001), and several laws have been used to model such behaviour (Costa et al. 2001;Guccione et al. 1995;Holzapfel and Ogden 2009). A good combination of model fidelity and parameter identifiability in estimation based on 3D tagged data is observed in the reduced form (Hadjicharalambous et al. 2014b) of the Holzapfel-Ogden law (Holzapfel and Ogden 2009): where a, a f , b and b f are constant parameters, I C v = C v : I is the first invariant of the right Cauchy-Green deformation with f denoting the myocardial fibre direction at a specific point in the domain. The values of the exponents are not estimated, but instead fixed at b = b f = 5, ensuring good agreement of the diastolic pressure-volume curve with the experimentally derived Klotz curve (Hadjicharalambous et al. 2014b;Klotz et al. 2006).
Similarly, a number of constitutive laws can be used to model active behaviour of the myocardium Kerckhoffs et al. 2003;Niederer et al. 2011). We sought to employ a simple version of these laws that would limit parameter uncertainties. Prescribing a single timedependent active tension α(t) over the whole myocardium led to non-physiological contractile behaviour of the model, but incorporating the effect of cell length dependence φ iso , as described in Kerckhoffs et al. (2003), allowed us to resolve the abnormality while keeping the model sufficiently simple:

Boundary conditions
Any external loading affecting the system can be incorporated into the model using the Lagrange multiplier method (Babuska 1973), with individual multipliers representing pressures on each of the boundaries. We can prescribe boundary conditions on the endocardium and the base of the left ventricle, so that Π ext = Π ext endo + Π ext base . The full Lagrange multiplier vector combines the endocardial and the base components: λ = (λ endo , λ base ).
The endocardial energy term penalises deviations of the model from given LV cavity volumes: where V (t) is the volume of the cavity computed from model displacements v (Asner et al. 2015), and V data (t) the volume computed from 3D tagged data, both at a given time t in the cycle. While it is common to drive cardiac mechanics models with cavity pressures (Nash and Hunter 2001;Wang et al. 2009;Xi et al. 2013), measuring these pressures non-invasively over the cardiac cycle is currently impossible. At the same time, cavity volumes can be easily derived from the time-resolved images. A volume-driven model can therefore be more easily personalised using routinely acquired patient data. Alternatively, a lumped parameter model can be coupled into the system in order to drive the cardiac cycle (Sainte-Marie et al. 2006;Kerckhoffs et al. 2007;Krishnamurthy et al. 2013). However, any such model requires additional personalisation based on the data that are already available, which is complicated by the 3D-0D coupling. The base energy term penalises deviations of the model from given displacements at the base of the ventricle: where λ base ∈ R 3 , u base are the deformations extracted from the images, and K b the penalty matrix. By choosing is the normal vector to the base plane at time t, we enforce strict adherence to u base in the n b (t) direction and relaxed adherence in the base plane controlled by ε. A weak Dirichlet condition on the base can be applied by setting ε = 0.

Finite element discretisation
A standard Galerkin finite element method is employed to find a numerical solution to Eq. (2), as described in Asner et al. (2015), Hadjicharalambous et al. (2014a) and Nash and Hunter (2001). The domain is discretised to produce a computational mesh Ω h , and the spaces in the weak form are substituted for their discrete counterparts: The discrete spaces are quadratic for displacements and linear for pressures in all tests to preserve stability in the mixed formulation.
The combination of the governing equations, constitutive relations and boundary conditions fully defines the cardiac mechanics model used in this work, with the finite element method providing a practical tool to solve the model equations. Personalisation of this model relies on a non-invasively acquired dataset, which is described below.

Non-invasive datasets
For a given set of patient-specific parameter values a, a f and α(t) the data required for simulating cardiac mechanics using the above model includes the following: 1. reference geometry Ω 0 , 2. fibre orientations f , 3. LV cavity volumes over the cycle V data (t), 4. displacements of the LV base plane u base .
The estimation of a, a f and α(t) relies on additional information: 5. local displacements throughout the entire LV and 6. LV cavity pressure estimates: a single value or a full-cycle trace, as discussed below.
Obtaining the reference geometry from image data is a non-trivial task, since the myocardium is not observed in its unloaded state during the cardiac cycle. In practice, endsystolic (Wang et al. 2009) and early (Mojsejenko et al. 2014;Xi et al. 2013), or even end-diastolic (Dokos et al. 2002) geometries have been used as reference for simulation purposes. The LV cavity volume at end systole and early diastole is normally close to the refence volume estimates. Moreover, passive parameter estimates were shown to be minimally affected by changing the reference state from end-systolic to early-diastolic geometries in an idealised LV in Hadjicharalambous et al. (2014b). At the same time, the effect of selecting the reference as an early-diastolic motion state, where the myocardium is known to experience residual active tension, needs to be quantified.
An alternative approach to obtaining the reference state is the solution of an inverse mechanics problem (Krishnamurthy et al. 2013;Xi et al. 2014). However, any such model is directly dependent on the values of material parameters estimated in the forward problem, as well as the boundary tractions, and it is not clear that the coupled problem is well posed. We investigated the potential for joint estimation of passive parameters and the reference state in "Appendix 1" and were unable to obtain satisfactory results for the available dataset. The observed lack of identifiability was indicative of strong coupling between the estimated values.
In the present pipeline, the end-systolic geometry was used as the reference state, ensuring consistent selection of the reference state between cases and balancing the accuracy of approximation with resulting parameter identifiability.
A commonly used linear fibre angle distribution between 60 • on the endocardial and −60 • on the epicardial surface (Spotnitz 2000;Streeter et al. 1969) was used in the model. In practice, it is possible to generate personalised fibres based on diffusion tensor MRI data (Rohmer et al. 2006;Nagler et al. 2013;Stoeck et al. 2014;Toussaint et al. 2013). However, at the moment scan duration as well as spatial resolution and accuracy of the acquired data is still prohibitive for clinical use.
The reference displacements of the computational mesh through the cardiac cycle were obtained by tracking the motion of the myocardium in 3D tagged images. The process is based on a non-rigid registration algorithm (Chandrashekara et al. 2004;Rueckert et al. 1999;Shi et al. 2012Shi et al. , 2013 and carried out with the Image Registration Toolkit 1 (IRTK). The tracking procedure provides a transformation for each voxel of the image. This can be applied to the vertices of a computational mesh positioned in the physical coordinates, producing the deformed states at the time points in the cardiac cycle where 3D tagged frames were available.
As discussed above, the model is driven by LV cavity volumes, since the corresponding pressure curve cannot be measured non-invasively. The volumes can be computed with reasonable accuracy for each deformed state of the computational mesh, and the pressures λ endo are then computed in the simulations. Due to the linear nature of the chosen parametrisation, scaling all of a, a f and α(t) by a constant factor does not affect the displacements, and scales the pressures λ endo and p by the same factor. In order to recover the correct absolute values of pressures and parameters, we can use a single known cavity pressure at a specific time point in the cardiac cycle. An example would be the peak systolic LV pressure (SP) estimated from cuff measurements using a Centron cBP301 device, which transforms acquired peripheral pressures to compute central pressures (Brett et al. 2012). In systole, when the aortic valve is open, the central pressure estimate can be used as an accurate LV pressure estimate. Alternatively, an end-diastolic pressure (EDP) estimate can be obtained from the E/E a ratio (Nagueh et al. 1997(Nagueh et al. , 2009). E is the peak early-diastolic flow velocity through the mitral plane and can be measured in phase contrast MR sequences. E a (sometimes denoted by E , or e ) is the early-diastolic velocity of the mitral annulus on the lateral side of the base and can be computed using the tracked motion of the mesh.
Moreover, both the EDP and the SP estimates, along with the timings of mitral and aortic valve opening and closing, can be used to generate a full-cycle patient-specific LV pressure curve from a normalised trace obtained in Russell et al. (2012), shown to be consistent across a range of cardiac disorders. Even though the accuracy of such pressure data, and hence its use as a driving force in the model, is limited (compared to the accuracy of image-derived cavity volumes), matching pressures along with displacements in the active estimation process proves beneficial (see Sect. 4).
The types of clinical data used in the simulation and estimation process in the in vivo tests in Sect. 3.2 were: 1 http://www.doc.ic.ac.uk/~dr/software/.
1. short-axis cine to produce the computational mesh, 2. 3D tags to extract mesh deformations, 3. pressure cuff measurement to estimate the SP, 4. mitral flow velocity to estimate the EDP.
The datasets were acquired as part of the BHF New Horizons project Integrated Mathematical Modelling and Imaging for Dilated Cardiomyopathy (DCM). Specifically, we processed one dataset from a healthy volunteer and two datasets from patients being treated for moderate DCMrelated heart failure. The cine, 3D tagged and flow images were acquired on a 1.5T Philips Achieva system with the following specifications: -cine bSSFP in retrospective ECG gating, spatial resolution 2 × 2 × 8 mm, temporal resolution ∼20 ms, FOV 350 × 350 mm. -3D tagged MRI in prospective ECG triggering, spatial resolution 3.4 × 7.7 × 7.7 mm, temporal resolution ∼30 ms, FOV 100 × 100 × 100 mm, reconstructed interpolated image with spatial resolution 1 × 1 × 1 mm. -4D flow 2 in prospective ECG triggering using an MRI breathing navigator, spatial resolution 2.3×2.3×2.3 mm, temporal resolution ∼35 ms, velocity encoding range 100-150 cm/s.
An end-diastolic LV mesh was created using the CGAL library 3 (The CGAL Project 2015) from a CardioViz3D 4 (Toussaint et al. 2008) segmentation of the end-diastolic frame of the short-axis cine image. The mesh was deformed using the results of IRTK motion tracking  from the 3D tagged image set covering the majority of the cardiac cycle. The mesh with the lowest LV cavity volume was used as the reference state for the simulations. The distances between the deformed meshes and this reference mesh at all time steps served as reference displacement data.
The end-diastolic pressure was estimated using the E/E a ratio: λ 0 = 1.9 + 1.24E/E a (Nagueh et al. 1997), with the peak early-diastolic flow velocity through the mitral plane E estimated using GyroTools GTFlow software 5 from 4D flow images, and the early-diastolic velocity E a computed at the lateral side of the base from reference displacements (assumed to be sufficiently close to the mitral annulus). The peak systolic pressure was estimated from cuff measurements using a Centron cBP301 device 6 (Brett et al. 2012). The fullcycle pressure curve was obtained by scaling the normalised pressure curve between λ max at maximum and λ 0 at the time corresponding to the first frame of the 3D tagged sequence. The time scaling was based on valve events extracted from the mesh volume curve, so as to obtain a physiological shape of the reference P-V loop.
Building the personalised models for the described in vivo datasets allows us to establish the potential applicability of the pipeline to clinical data and the feasibility of interpreting the estimated material parameters as characteristic of tissue properties, and consequently cardiac health and disease.

In silico datasets
The evaluation of the proposed pipeline using typical clinical datasets is the ultimate test of its applicability. However, such evaluation is often obscured by variable data quality. The lack of ground truth for most computed variables as well as model parameters means that validation using medical images is somewhat limited in scope. In silico testing, whereby the datasets are generated using the model, provides a useful tool for bridging the gap between model and clinical measurement. In this case, when all of the reference parameters are known and the model has perfect fidelity, we can investigate the accuracy of the estimation procedure as well as its sensitivity, specifically the effects of noise in the data and motion tracking errors on the quality of parameter estimates.
An idealised left ventricle was modelled as an ellipsoid cropped at an angle at the base, as shown in Fig. 1a. This geometry can be easily generated for testing purposes and provides a straightforward point for comparison. The idealised fibre field was produced as described for the in vivo datasets.
The displacement data were produced in three stages, as illustrated in Fig. 2. First, we obtained clean displacements ( Fig. 2a) by running a full-cycle simulation using a given volume curve with given reference parameters a ref , a ref f and a given active scaling function α ref (t). The input data were chosen so that the model produced physiological cavity pressures, P-V loop and LV deformations. The mesh underwent characteristic twist and shortening (Fig. 1d) in systole before relaxing (Fig. 1e) and inflating (Fig. 1c) in diastole.
Second, random zero-mean uniform noise with a given standard deviation was added to clean displacements to produce noisy displacements (Fig. 2b, c). The random nature of the added noise means that its spatial distribution is not in any way related to any local geometric features of the mesh. This ensures that the symmetry of the idealised geometry does not artificially improve parameter identifiability or the accuracy of estimation. Third, we generated artificial 3D tagged images (Fig. 2d-f) using clean and noisy displacements and performed motion tracking as described above to extract processed displacements from the images (Fig. 2g-i).
The full pressure curve was produced as simulation output together with clean displacements.

Model personalisation process
The full-cycle estimation process is carried out via the following steps: 1. estimation of the passive parameter ratio from diastolic displacements, 2. scaling the passive parameters and estimated pressures using the reference EDP, 3. estimation of the active tension from a combination of the displacements and LV cavity pressures through the cycle. Comparison of in silico datasets in mid-systole (long-and short-axis views). The surfaces show the domain deformed using clean displacement data, while the lines crossing the surfaces show the mesh deformed using the six different datasets produced for testing: clean unprocessed, unprocessed with noise at 10 and 20 % standard deviation, clean processed and processed noisy at 10 and 20 % standard deviation. Surface shading represents the distance between the clean displacements and the displacement for each of the cases. a Clean, b noisy, 10 % std, c noisy, 20 % std, d clean tags, e noisy tags, 10 % std, f noisy tags, 20 % std, g processed clean, h processed noisy, 10 % std, i processed noisy, 20 % std The procedure is aimed at finding parameters that produce simulation results closest to the data, as measured by objective functions J .
Let u ref n and λ ref n , n = 0, . . . , N , be the reference displacements and cavity pressures, respectively, n = 0 corresponding to end diastole and n m to end systole.
Since the passive parameters a and a f are assumed to be constant in time, the objective function in the passive stage is defined as the relative total displacement error: with · the standard L 2 (Ω 0 )-norm. The effect of residual active tension at the start of diastole is neglected due to our inability to reliably deduce its presence in vivo. This assumption can be corrected for based on active estimation: if the active scaling is estimated at a significant positive value in the initial diastolic frames, then the passive parameters can be adjusted based on estimation in the remaining diastolic frames (see Sect. 4).
As discussed above, in volume-driven simulations the pressure solutions scale together with a and a f , so that displacements alone do not identify both passive parameters, only the ratio between them γ = a/a f . In practice, a fixed value of a sim f can be used to find a sim which minimises the objective function. 7 The correct absolute values of parameters a and a f , as well as estimated cavity and myocardial pressures λ endo and p, are recovered by scaling the respective values used or computed in the simulations by the ratio of the reference to the estimated EDP λ ref 0 /λ sim 0 . In the current procedure, this scaling will also be consistent with the peak systolic pressure, since the personalised full-cycle pressure curve is scaled to both the EDP and the SP. The correct absolute values of the base pressure λ base can be obtained by running inflation with the correct absolute parameter values.
In the active stage, the parameter α is estimated through time, or, more precisely, at times when displacement data for n = 0, . . . , N . Equal weighting is given to relative displacement and pressure errors in the combined functional in order to balance the uncertainty in the two variables, but this can be adjusted if deemed necessary. The displacement errors are taken in relation to the maximum over all frames, not current displacement, in order to neutralise the effect of zero and low displacements at and around end systole. In addition to computing objective functions which measure relative errors in problem variables, for the in silico tests, we can compute the relative errors in the parameter estimates: These metrics are useful in assessing the quality of the estimation and the sensitivity of the model.

Results
First, we describe test specifications (including geometry, input data, parameters and generation of clean reference data) and present results for the in silico cases. Model outputs were compared to reference data, and the effects of noise and tags processing were assessed. Next, we describe the in vivo datasets processed and present pipeline outputs for these cases. All simulations were performed using CHeart (McCormick et al. 2011(McCormick et al. , 2013, 8 a multiphysics solver package developed at the Biomedical Engineering Department, King's College London.

In silico cases
The computational mesh was composed of 56 quadratic hexahedral elements. 9 The cycle was taken to last 1000 ms, driven by the cavity volume curve shown in Fig. 3a Fig. 3a, and the reference P-V loop in Fig. 3c.
Reference simulation results were obtained at 5-ms intervals to achieve convergence of the nonlinear solver. However, only solutions at 30, 70, . . . , 830 ms were taken as data points, making a total of 21 "frames" to emulate a 3D tagged image sequence typically obtained using the in vivo imaging protocol by means of temporal sampling and initial offset given by the prospective ECG triggering.
Two sets of unbiased uniform noise were generated for each set of clean data: with standard deviations of 10 and 20 % of maximum displacement length, where the maximum was taken for each displacement to end-diastolic state (ranging between 6 and 12 mm).

Estimation of the passive parameters
All diastolic frames (t = 510, 550, . . . , 830 ms) were used for the estimation of the passive parameter ratio. In addi-8 http://cheart.co.uk/. 9 The structured hexahedral mesh consisted of two conforming layers: the majority of the myocardium was subdivided into 8 elements circumferentially, 3 in the long-axis direction and 2 transmurally, and the apical patch subdivided into 2 × 2 × 2 elements (see Fig. 1a). This resolution was sufficient to characterise objective function behaviour, as observed in Hadjicharalambous et al. (2014b).  Unprocessed (0 ± 1) % (7 ± 1)% (11 ± 1)% Processed (13 ± 1) % (12 ± 1)% (14 ± 1)% tion, the first frame (t 0 = 30 ms) was assumed to be sufficiently close to the end-diastolic state to provide a data point at t = 1000 ms. Since the reference volume changes between frames were too large to apply in single simulation steps, additional linearly interpolated volumes were used. Even though the exact reference volumes were known at intermediate steps between frames, they were not used in estimation, so that the dataset resembled the in vivo case.
For each of the six datasets, a parameter sweep consisting of ∼25 simulations was carried out in order to assess the behaviour of the objective function J p over the γ = a/a f space (see Fig. 4) and find its minimum (expected at or near γ ref = 0.1). The fibre parameter was fixed at a sim f = 1 kPa in all simulations.
The estimated values of e p [γ ] are given in Table 1 with uncertainty ranges determined by sweep intervals.
Scaling the estimated a sim , a sim f , λ sim endo and p sim by λ ref 0 /λ sim 0 = 1.6 allowed us to recover the exact EDP.

Estimation of the active tension
The estimated passive parameters for each dataset were used in active scaling estimation. We performed parameter sweeps for each time frame in order to map the behaviour of the active objective function. In the systolic frames with nonzero reference active tension α ref  Fig. 5b. In this case, the minimum was always expected at or near α = 0.
The mean and maximum errors in the active scaling estimated using the combined objective function J a are given in Table 2 for the systole and Table 3 for the diastole.
Sample comparison plots (worst case, processed displacements with 20 % noise) for the reference and estimated active scalings and P-V loops are shown in Fig. 7a, b.

In vivo cases
We processed three datasets acquired in a healthy volunteer (V) and two patients with dilated cardiomyopathy (P1 and P2). Case details are presented in Table 4.
The meshes were composed of 7795 (V), 10740 (P1) and 17047 (P2) linear tetrahedral elements. The reference volumes were computed from meshes deformed according to 3D tagged image tracking, and reference pressure curves were scaled between estimated end-diastolic and peak systolic pressures, and in time based on volume-derived valve events in order to obtain physiological shapes of the reference P-V loops (see Fig. 8).
Weak adherence to reference base displacements was required in all simulations (ε = 10 −7 −10 −8 ), and the reference mesh volumes were used in the endocardial boundary condition.

Estimation of the passive parameters
As before, estimation of the passive parameter ratio was carried out as a parameter sweep, refining the sweep intervals   Estimated active tension does not exceed 1 kPa, which is less than 1 % of the peak active tension, allowing a differentiation between the passive and the active stages of the cycle  HR heart rate, EDV/ESV end-diastolic/end-systolic volume, EF ejection fraction, SP peak systolic pressure, EDP end-diastolic pressure in the area of lower J p . The fibre parameter was fixed at a f = 1 kPa. All diastolic frames were used in computing the objective function, as well as the end-diastolic state obtained from the first frame of the 3D tagged images. Figure 9 shows the variation of J p over a range of passive parameter ratios γ .
The minimising ratios γ along with the passive parameter values a and a f scaled by the ratios of the reference and the simulated EDPs are shown in Table 5.

Estimation of the active tension scaling
A parameter sweep was carried out at each time step corresponding to a frame of the 3D tagged image. At a given time step, a simulation was carried out whereby the active tension scaling increased from 0 to 200 kPa using 1 kPa increments (i.e. each simulation took 200 parameter increment steps). The variation of J a over the parameter space for each time step is shown in Fig. 10. The lowest and highest absolute values of the combined displacement and pressure objective function in systole together with peak estimated active tension and the relative time to peak activation are given in Table  5 for all cases. (a) (b) (c) Fig. 9 Variation of the objective function J p over the passive parameter ratio space in the in vivo cases, with minima shown as circles Table 5 Estimation results for the in vivo cases. α max denotes the peak active tension,t max the scaled time to peak activation:t max = t max · (HR)/60, and min t and max t J a refer to the lowest and highest absolute values of the objective function across all systolic frames at the minimum over α  Figure 11 shows the reference mesh configurations for case V1 at different stages of the cardiac cycle based on 3D tagged motion tracking alongside simulation results obtained using estimated parameters.

Discussion
Patient-specific cardiac modelling holds significant promise for providing details on the performance, function and properties of the heart. However, the process of personalising mechanical models needs to be reliable and robust in order to become a viable tool for model-based analysis of patient cases. The modelling pipeline proposed in this work is based entirely on a non-invasive dataset, which includes a shortaxis cine MRI stack (used to construct the LV mesh) and a 3D tagged MRI sequence (providing mesh deformations over time) as well as end-diastolic and end-systolic pressure estimates derived from flow MRI and a cuff pressure measurement. The main goal of the study was to ensure the existence of unique sets of parameters that provide a good fit to the three-dimensional cardiac motion over the cardiac cycle. The presented model aims to balance model fidelity with parameter identifiability, producing mechanical equations that depend linearly on the parameters. Identifiability was verified and validated by sweeping parameter spaces in a series of tests where in silico and in vivo data were matched to model predictions. Unique parametrisation was possible in all cases, providing two constant stiffnesses (defining passive behaviour) and a time-dependent active tension (driving the active contraction) estimated based on a rich set of displacements across the myocardium and an LV cavity pressure curve estimate ensuring physiological behaviour of the model. While the constitutive laws used were highly nonlinear with respect to problem variables, their dependence on the estimated parameters was linear, resulting in each of the objective functions J having a single global minimum. In all cases, model fidelity was sufficient to retain this property of J ensuring our ability to estimate reliable parameters.

In silico cases
Testing with in silico datasets provided a useful insight into how the proposed model and parameterisation performed. The in silico tests (see Fig. 2) enabled a decoupling between the issues of model fidelity and parameter identifiability, since the model was tested under the circumstances where it was able to replicate reference deformations. The knowledge of the reference parameters used for generating the data allowed the assessment of the sensitivity of the model parameters to different types of errors (in addition to the sensitivity of the problem variables).
In active estimation, the identification of active tension scaling throughout the cardiac cycle was uniformly achieved (see Fig. 5). The inclusion of random displacement noise of 10-20 % resulted in a minimal increase in error, while errors in pressure values translated directly to active parameter errors. Activation through diastole was also minimal, showing a bias towards activation using processed data (though this was below 0.01 % peak tension). It is interesting to note that identification of myocardial activation was achievable from displacement data alone (see Fig. 6), suggesting that, given sufficient model fidelity, a full-cycle pressure curve would not be required for estimation. However, the objective function based on both displacement and pressure errors provided improved accuracy (max 5 %, average 0.5 % error for 20 % noise in data) compared to that based on the displacement term alone (max 25 %, average 20 % error for 20 % noise in data). In the latter case (which would allow fullcycle estimation based on a single pressure measurement), peak systolic pressure was underestimated by 25 %. While scaling using the measured to estimated pressure ratio would adjust the peak active tension to the correct value, it would also overestimate both passive parameters by an additional ∼33 %, as well as overestimate the active tension curve away from the peak up to ∼33 %.

In vivo cases
Following the pipeline introduced in Sect. 2.2, a unique parametrisation of the model was obtained for each of the volunteer and patient datasets acquired non-invasively. The myocardial geometry was obtained from the short-axis cine stacks registered to 3D tagged images and tracked through the cycle using IRTK, providing observations for estimation as well as boundary condition data, namely LV cavity volumes (see Fig. 8) and basal deformations. Following the parameterisation process, Figs. 9 and 10 show that unique identification of the passive ratio γ and active tension scaling α(t) was achieved. Moreover, Fig 11 illustrates that qualitative model motion corresponded well with extracted motion from the patient data.
The active estimation results indicated that residual active tension was present in the myocardium in 4-5 frames past end systole. This additional information could be used to correct passive parameter estimates. For example, the objective function J p was recomputed for case V based on the frames with no residual tension (15-24), and the minimum was achieved at γ = 0.08. After scaling by end-diastolic pressure ratio, the corrected parameter values were a = 0.304 kPa and a f = 3.8 kPa, including relative corrections of 10 % for a and 3 % for a f .
Unlike the in silico tests, where good model fidelity allows estimation of the active tension based solely on the tracked myocardial data, the in vivo cases relied on the pressure curve in producing a viable active tension estimate. The minimum of the displacement component of the active functional was achieved at near-zero active tension throughout the cardiac cycle for all datasets due to model errors dominating the term. The use of non-invasive pressure estimates ensured that a reasonable pressure-volume loop was produced by the model (see Fig. 8).
While we cannot make reliable interpretation of parameters based on three datasets only, some differences can already be observed in the results in Sect. 3.2. The estimated passive parameter ratio is significantly lower for the healthy volunteer, and consequently we estimate that the tissue is much stiffer in the fibre direction than isotropically, which does not appear to be the case for the DCM patients. Case P1 activates significantly slower, and case P2 significantly higher than V. The higher total relative error J p and the shallower curves near the minimum in the passive estimation for P1 and P2 as compared to V could be caused by reduced motion of the diseased myocardium and therefore a stronger effect of model errors. Identifiability is not affected in the active estimation due to reliance on the same normalised pressure curve. A detailed analysis of estimation results for a wider range of healthy and diseased cases would be fundamental for building an understanding of the link between the estimated quantities and physiology.

Limitations and future work
In the current study, a range of in silico datasets and three in vivo cases are considered in order to validate unique identifiability of parameters in the personalisation process. A comparative analysis of parameter estimates for a wider range of cases including healthy volunteers and patients is necessary for determining the diagnostic and prognostic potential of the model parameters as biomarkers. While personalisation of multiple in vivo datasets using parameter sweeps would incur a high computational cost, the results of this work suggest that commonly used optimisation methods, such as the reduced-order unscented Kalman filter Moireau et al. 2008;Xi et al. 2011), would perform well in the same parameter spaces using the outlined objective functions.
A direct comparison of the in vivo parameter estimates with results available in the literature for the Holzapfel-Ogden law would be problematic due to differences in some of the basic assumptions on fixed parameter values taken in different studies. Further effort in this direction can be made based on recent experimental results in Sommer et al. (2015).
The objective functions used to determine the optimal parameters for a given case took advantage of the local information on the myocardial position over the cardiac cycle extracted from 3D tagged images. To our knowledge, this is the first instance when the relative L 2 (Ω 0 )-norm of the error between myocardial deformations obtained using image processing and model predictions has been assessed and reported. In the passive parameter estimation, the error is accumulated over time as well as space, leading to the values of the objective function min J p = 45-63 %. Active estimation relies on a combined displacement and pressure error functional, with the pressure component ensuring parameter identifiability throughout the cycle.
The strictness of the objective function accumulating all local errors can partly explain the relatively high values of J . Nevertheless, further adjustment of the model should be sought to improve fidelity and specifically to reduce displacement errors. In particular, local errors were largest in the septum and at the right-ventricular attachment points, suggesting that the influence of the right ventricle should be accounted for when considering left-ventricular mechanics.
Additional tests have been carried out to assess the potential for reducing model errors (thus improving model fidelity) by using different fibre distributions in the myocardium, or introducing a transverse activation law (see "Appendix 2"). As may be expected, the absolute values of the parameters obtained using these modified models change to some extent. More importantly, the identifiability of parameters is not dramatically affected by either modification. The viability of using regional passive parameters and non-homogeneous activation should be investigated.
The choice of the end-systolic geometry as the reference state is likely to have an effect on displacement errors, as well as passive estimation results. An investigation of the bias introduced by this choice would be of interest. While the present set-up did not permit simultaneous estimation of the passive parameters and the reference state, it is possible that refinement of the model and/or enrichment of the dataset would make it a viable option.
The need for an estimate of the LV pressure curve for active estimation means that further validation of the employed scaled normalised data would be beneficial. While the curve in Russell et al. (2012) was shown to be similar for patients with a number of different heart conditions (and consistent between diseased and healthy canines), its application to personalised models of healthy volunteers is yet to be fully verified. Any errors in the pressure will translate linearly into errors in all parameter estimates due to their linear dependence prescribed by the model.
Experimental validation of these results would be an important step in ensuring that model-based estimates and predictions can be relied upon in clinical practice.

Conclusions
The proposed pipeline for personalisation of the full-cycle LV mechanics model produced unique passive and active parameters for all of the in silico and in vivo datasets tested. Moreover, the behaviour of the objective functions over the respective parameter spaces implied that an efficient optimisation or filtering method can be successfully used to find the patient-specific parameters of the model. The choice of the L 2 (Ω 0 )-norm of displacement error as the objective function in parameter estimation allowed us to take advantage of the rich dataset extracted from 3D tagged MR images. In the in silico tests, the presence of random noise of up to 20 % and typical processing errors caused a moderate deterioration of model fidelity, without a significant effect on the objective function behaviour. The associated errors in parameter estimates were within the data noise/error level. In vivo estimation results showed good agreement with global measurements, such as the P-V loop and the overall motion of the ventricle.   Three sets of fibres were generated for the in vivo case V, all linearly varying through the wall, with angles ranging between 55 • and −55 • , 60 • and −60 • , and 65 • and −65 • going from the endo-to the epicardial surface.
The estimation procedure was carried out as described in Sect. 2.4. The variation of the passive objective function for the three tests is shown in Fig. 14. The passive parameter values scaled by end-diastolic pressures are shown in Table 6.
We observed that the total displacement error J p at the minimum reduced as the fibre angles were increased, as did the optimal value of the passive ratio. Increasing the fibre angles in diastole increased dilation of the LV in short axis and reduces elongation. The reduced ratio γ balanced this by reducing dilation and increasing elongation. Identifiability of Fig. 15 Estimated active tension using fibre-only activation, and fibre and transverse activation for case V parameters was not dramatically affected, and the change in error at the minimum of J p was <2 %.
The estimated active tension profile and identifiability in J a were not significantly affected by the chosen fibre distribution. The characteristics of the optimal α(t) and the associated J a are shown in Table 6, and we also note that the highest difference between the estimated tensions for the different fibre orientations was 3 kPa at peak activation. A personalised fibre distribution, such as that based on diffusion tensor MRI, might have a stronger effect on model fidelity; however, this type of data was not available in the current datasets.

Transverse activation
The simple active model used in this work and a lot of the more complex models assume that active systolic stress is developed mostly in the direction of the fibres. However, it has been suggested through experiment and computational modelling (Usyk et al. 2001(Usyk et al. , 2002 that activation at the level of ∼30 % of fibre stress occurs in the directions transverse to the fibres. We have carried out a test for case V in order to assess whether the inclusion of the transverse components in the active constitutive law has the potential to improve model fidelity. As a first approximation, the tension in each of the transverse directions was taken to equal 30 % of the tension in the fibre direction. A comparison of the estimated active tension scaling α(t) obtained using the two methods is shown in Fig. 15.
The range of total errors J a in systole is reduced at the minimum over α from 24-35 to 18-26 %. The identifiability using the objective J a remains unaffected, and identification of a physiological active tension based on the displacement component of the active functional is still impossible.