Analysis of passive cardiac constitutive laws for parameter estimation using 3D tagged MRI

An unresolved issue in patient-specific models of cardiac mechanics is the choice of an appropriate constitutive law, able to accurately capture the passive behavior of the myocardium, while still having uniquely identifiable parameters tunable from available clinical data. In this paper, we aim to facilitate this choice by examining the practical identifiability and model fidelity of constitutive laws often used in cardiac mechanics. Our analysis focuses on the use of novel 3D tagged MRI, providing detailed displacement information in three dimensions. The practical identifiability of each law is examined by generating synthetic 3D tags from in silico simulations, allowing mapping of the objective function landscape over parameter space and comparison of minimizing parameter values with original ground truth values. Model fidelity was tested by comparing these laws with the more complex transversely isotropic Guccione law, by characterizing their passive end-diastolic pressure–volume relation behavior, as well as by considering the in vivo case of a healthy volunteer. These results show that a reduced form of the Holzapfel–Ogden law provides the best balance between identifiability and model fidelity across the tests considered.


Introduction
Cardiac imaging provides a powerful tool for assessing cardiac function and pathology, offering valuable insight into the kinematics and tissue characteristics of the heart. Its joint use with multiscale mathematical models of the heart (Costa et al. 2001;Guccione et al. 1995;Nash and Hunter 2000;Holzapfel and Ogden 2009; enables the quantification of model constitutive parameters, which can be used as clinical biomarkers of disease Xi et al. 2011b;Sermesant et al. 2006;Chabiniok et al. 2012;Imperiale et al. 2011). As a result, there is a strong need for reliable parameter estimates, an issue which depends on both the underlying cardiac constitutive model as well as the available clinical data.
A variety of noninvasive cardiac imaging techniques have been developed, which offer a powerful set of tools for personalization of mechanical parameters. Imaging methods such as echocardiography, computed tomography (CT), cardiac magnetic resonance imaging (MRI) have been successfully employed to accurately capture epicardial and endocardial motion, thus providing important bulk metrics such as ejection fraction and cavity volumes. In addition, the development of speckle tracking echocardiography (Meunier 1998;Craene et al. 2012) and cardiac MR tagging [SPAtial Modulation of Magnetization, SPAMM (Zerhouni et al. 1988;Axel and Dougherty 1989;Reichek 1999)] has enabled the quantification of regional cardiac motion in vivo (Young et al. 1995;Osman et al. 1999;Arts et al. 2010), by revealing local characteristics such as wall thickening, torsion and shear effects. The translation from 2D to 3D tagging techniques (Rutz et al. 2008) has enabled a direct extraction of the full 3D displacement field in the myocardium, leading to simultaneous quantification of radial, circumferential and longitudinal motion (Shi et al. 2012;Pan et al. 2005). This accurate "whole-ventricle" 3D deformation field creates an ideal setting for estimation of model-based parameters using tissue displacement observations.
The wealth and quality of information on myocardial motion has enabled patient-specific applications, where passive cardiac constitutive laws of varying complexity have been employed, ranging from simplified isotropic (Cheng et al. 2005) and transversely isotropic (Guccione et al. 1991) laws to orthotropic models accounting for the fiber anisotropy of the tissue (Holzapfel and Ogden 2009;Costa et al. 2001;Nash and Hunter 2000). However, as the model complexity and number of parameters increase to better approximate the tissue's complex behavior, estimating model parameters uniquely and accurately becomes an increasingly challenging task (Xi et al. 2011b). This raises the important question of structural identifiability for the various constitutive laws, i.e., whether it is possible to uniquely determine parameter values, given infinite well-defined noise-free data (Chis et al. 2011;Raue et al. 2009). Structural identifiability-a property of the model itself which does not depend on the available data-can be compromised by coupling between model parameters as in the case of the Guccione model Xi et al. 2011a, b;Augenstein et al. 2005) and nonlinear dependence of the model on the parameters. Lack of structural identifiability hinders the ability of any data assimilation method-mainly categorized into variational (Sun et al. 2009;Augenstein et al. 2005;Wang et al. 2009;Sermesant et al. 2006) and sequential (Moireau et al. 2008(Moireau et al. , 2009Chabiniok et al. 2012;Xi et al. 2011b;Wong et al. 2007;Liu and Shi 2009)-to accurately estimate parameter values.
In a clinical scenario, the estimation process is further compromised by limited data and measurement noise, leading to the issue of practical identifiability, i.e., whether we can determine unique parameter estimates given the limited amount and quality of data (Saccomani 2013). Absence of structural or practical identifiability in a cardiac law given a set of data leads to unreliable parameters, which can no longer provide clinically relevant information. The choice of an appropriate cardiac constitutive law should therefore balance the need for model fidelity, i.e., the ability of the model to accurately represent cardiac function, with the requirement for reliable identifiable parameters.
In this work, we aim to assist the choice of an appropriate constitutive law when the available data are 3D tagged MRI, by examining the practical identifiability and model fidelity of different cardiac mechanics models. Specifically, we look to compare progressively complex models to assess their capacity to both represent cardiac motion and be used reliably for parameterization. In order to gain insight into the parameter estimation process for these models, we investigate the behavior of a minimization criterion (J ) over the parameter space. This is first tested using synthetic 3D tags extracted from in silico simulations and performing parameter sweeps to obtain the parameter estimates. Following the work-flow described in Fig. 1, we characterize J and assess the error between estimated and actual parameters. The various models are further compared with respect to their ability to represent physiological cardiac deformation (model fidelity) and end-diastolic pressure-volume relation (EDPVR), in order to identify a constitutive law that would balance between practical identifiability and adequate representation of cardiac behavior. Our conclusions are then verified in an in vivo case, validating our in silico framework for characterizing practical identifiability using 3D tags.
Below, we expand on our approach to investigate practical identifiability and how it is influenced by the choice of constitutive law. The process for characterizing practical identifiability for each one of the considered models is reviewed in Sect. 2 and employed for in silico tests of diastolic filling using an idealized left ventricle (Sect. 3). The study is then extended to an in vivo case of a healthy volunteer, enabling the characterization of practical identifiability and model fidelity in a real-world scenario.

Methods
In this section, we describe the process followed in this work in order to assess the practical identifiability of various laws, focusing on the creation of synthetic tags, the motion extraction algorithm used, and the parameter sweeps performed (Sect. 2.1). We then present the cardiac model of LV diastolic filling used, as well as the various cardiac constitutive laws considered (Sect. 2.2). Finally, we review a general theoretical framework for the inverse problem of parameter estimation using 3D tags (Sect. 2.3), focusing on the concepts of structural and practical identifiability, and the factors that influence them (observations, constitutive laws, objective function).

In silico tagging and assessment protocol
A primary goal of this study was to assess the potential of using 3D tagged MRI in parameter estimation applications. Even though 3D tagged MRI offers a rich dataset for parametrization, the process may be compromised by low-resolution or noisy data and error introduced during the motion-tracking Workflow followed for the study of practical identifiability using 3D tags. The in silico testing protocol is presented in blue, while in red is the pipeline followed in the in vivo case procedure. In order to investigate this issue, we have created synthetic 3D tagged images directly from simulation results. Within this controlled environment, the actual parameters of the heart model are known, allowing for an assessment of the error between actual and estimated parameters. Further, as the synthetic tags approximate real 3D tagged images (see Fig. 2), within this framework, we can quantify the error associated with various aspects of 3D tags such as resolution, number of tag lines, noise in the data, and error introduced by the tracking algorithm.
Initially, we ran a simulation of LV diastolic filling, choosing parameters that produce a physiological enddiastolic volume. Synthetic 3D tags were then generated from the resulting deformation, and the myocardial motion was extracted and propagated on a mesh. These deformed meshes were then treated as data and were used directly for comparisons with simulation results. By performing parameter sweeps, computing and minimizing J over a bounded parameter space, we obtained parameter estimates and quantified the behavior of J .

LV diastolic filling
Several constitutive laws (see Sect. 2.2.2) were employed to model the passive behavior of the tissue (step "SIMULA-TION" in Fig. 1), and the simulated diastolic deformation in each case was compared with diastolic data to provide parameter estimates. The LV was modeled as a thick-walled truncated ellipsoid of typical cardiac dimensions. The domain was discretized using a mesh composed of 56 hexahedral ele-ments, with two elements transmurally, and a quadratic linear interpolation scheme was employed for the displacement and pressure variables, respectively. A generic heterogeneous fiber field was applied, with the fiber angle varying linearly between 60 • and −60 • from endocardium to epicardium. A zero traction condition was enforced on the epicardial surface, while the base plane was fully fixed. Finally, the endocardial surface of the LV model was passively inflated to a typical end-diastolic pressure of 1.5 kPa, using 50 uniform loading steps.

In silico assessment protocol
As we are interested in the passive cardiac parameters, we have generated synthetic tags (step "SYNTHETIC TMRI" in Fig. 1) from a passive inflation simulation of a model left ventricle (LV) (see Sect. 2.1.1) using the various constitutive laws, which will be described in Sect. 2.2.2. Using rasterization, a binary mask of the mesh was created, and tag planes were generated within the image (Sermesant et al. 2003;Duan et al. 2007), resulting in a final 3D tagged image of resolution 1×1×1 mm. Simulated deformations were then mapped and interpolated within the image, producing a 3Dtagged representation of our passive inflation simulation. For the remaining steps of the parameter estimation study, these synthetic tags served as our data and were treated as real 3D tags.
In order to assess the effect of data noise, Gaussian noise was added to the simulation results, prior to the in silico tagging (step "UNBIASED NOISE" in Fig. 1). The mean value of the added Gaussian noise was zero, and the variance was a percentage (usually 5-20 %) of the maximum deformation of the diastolic simulation.

3D tagged MRI/in silico tagged motion extraction
Myocardial motion was extracted from 3D tagged and in-silico tagged images (step "MOTION TRACKING" in Fig. 1) using the Image Registration Toolkit. 1 This software uses a nonrigid registration technique proposed by Rueckert et al. (1999); Schnabel et al. (2001) and subsequently extended to tracking of cardiac motion (Chandrashekara et al. 2004;Shi et al. 2012). The registration algorithm which is based on freeform deformations and optimizing the similarity between two subsequent images allows tracking any point within the myocardium throughout the cardiac cycle and provides the deformation field with respect to the beginning of systole. The obtained deformation fields were then applied on an initial mesh and propagated through time, resulting in a set of deformed meshes, which were used as the "observations" within our parameter estimation process. 1 http://www.doc.ic.ac.uk/~dr/software.

Mechanical simulations and J characterization
The parameter estimates for the constitutive laws considered were obtained by parameter sweeps (step "PARAMETER SWEEP" in Fig. 1). Specifically, for each law, we performed passive inflation simulations ( see Sect. 2.1.1) for parameter combinations within a neighborhood of the actual parameters. The parameter estimates were then obtained as the set of values that minimized the objective function J over the parameter space. Within this process, 10 synthetic tagged frames were used as observations. The objective function employed-defined and discussed in Sect. 2.3.2-is discerning and thus able to provide a unique minimum, assuming that the constitutive law is practically identifiable. Note that as the same constitutive law is used for the generated data and simulations, the estimation process does not suffer from model fidelity issues, leading to safe conclusions about practical identifiability.

In vivo J characterization
The study was then extended to an in vivo case of a healthy volunteer, to allow for the assessment of practical identifiability and model fidelity in a real-world setting.
The data used in this study were acquired from a healthy 28-year old male volunteer with a normal heart function. The geometry and cardiac motion were characterized using cine MRI images in retrospective ECG gating (short-axis view, acquired spatial resolution 2 × 2 × 8 mm, temporal resolution ∼20 ms), and 3D tagged MRI images in prospective ECG triggering (acquired in 3 breath-holds and reconstructed to spatial resolution of 1 × 1 × 1 mm, temporal resolution ∼30 ms). The myocardial motion was extracted from the 3D tagged MRI images using the algorithm described in Sect. 2.1.3.
The cubic Lagrange hexahedral mesh used for the simulations was based on the end-systolic frame of the data ) and was composed of 72 elements, with two elements transmurally. This mesh was created by first constructing a cubic hexahedral end-diastolic mesh ) from a short-axis CINE MRI stack registered to the 3D tagged MRI images using the IRTK imaging toolkit (see Sect. 2.1.3). The mesh was then uniformly refined into a linear end-diastolic (ED) mesh consisting of 96,768 hexahedral elements. Nodes of the linear ED mesh were then tracked through the cardiac cycle, providing the deformed LV geometry at end systole. This deformed ED mesh was used as a template for the coarser 72 element end-systolic (ES) cubic mesh, which was generated by least-square fitting. Observations-in this case displacements from end systole-were then required to compute the objective function (see Eq. 15). These were extracted from the linear ED mesh by subtracting the displacements observed at the endsystolic frame from the diastolic displacements. These computed displacements at the linear ED mesh nodes were then mapped onto a linear version of the ES mesh by doing a nearest point search (mean point search error 0.8 mm). A large number of elements were used in the linear ED mesh to minimize potential errors due to the mapping of data. The final projected displacements on the linear ES mesh served as observations within our parameter estimation process.
The cubic end-systolic mesh was inflated by prescribing the data-derived cavity volume at each time step, instead of inflating by pressure as used in the in silico tests. This was due to the lack of cavity pressure measurements for this in vivo test. The volume constraint was enforced weakly through a Lagrange multiplier. The base motion was prescribed directly from the observations-13 diastolic frames based on increasing cavity volume-and zero traction was applied on the epicardial surface. The myocardial model was assumed to be incompressible (see Sect. 2.2.1), even though potential compressibility of the extracted myocardial motion was not examined. This is due to the fact that a possible degree of compressibility [2-12 % (Iwasaki et al. 1984)] would be within the level of noise of the tagged MRI data. Running the described simulation with parameter sweeps and comparing with observations, we could then characterize the behavior of J in a real-world scenario.

Finite elasticity
The passive diastolic filling of the LV considered in this work was simulated within the finite elasticity framework due to the large deformation of the myocardium during the cardiac cycle (Holzapfel and Ogden 2009).
We consider here a body defined over a reference domain Ω 0 , which deforms under the action of a traction t (such as the endocardial pressure) on a subset Γ N of the boundary ∂Ω, where Ω is the current configuration. Given a set of parameters θ related to the employed constitutive law, the mechanics problem can be written as: Find the deformation and hydrostatic pressure pair where In the finite elasticity framework considered here, σ denotes the deviatoric Cauchy stress tensor. In this setting, x = (u, p) and y = (v, q) represent the state solutions and test functions, respectively.
The Cauchy stress tensor introduced in Eq. 1 depends on the passive behavior of the material and the constitutive law chosen to describe it. As the myocardium is most typically modeled as a hyperelastic material, the mechanical behavior is expressed using a strain energy function, whose deviatoric component is denoted here by Ψ . The deviatoric component of the Cauchy stress tensor can then be obtained through the expression Here F denotes the deformation gradient defined as which relates a point in the deformed state x ∈ Ω to its reference configuration X ∈ Ω 0 (Bonet and Wood 2008). For volume-preserving materials, incompressibility is enforced using the determinant of F through the constraint J = det(F) = 1. Additionally, C represents the right Cauchy-Green deformation tensor, defined as C = F T F. In cardiac mechanics, the finite elasticity problem (1) is commonly solved using the finite element method (FEM) due to the nonlinearities introduced by the constitutive laws and cardiac geometry. The FEM framework is based on discretization of the continuous domain and function spaces (Hadjicharalambous et al. 2014).

Constitutive laws
We begin with the Neo-Hookean law, a well-known isotropic hyperelastic law, which has also been employed in cardiac models (Cheng et al. 2005). The strain energy function Ψ for the Neo-Hookean law is defined as where μ is the stiffness of the material,Ĉ = J − 2 3 C is the isochoric component of C and IĈ is the first invariant ofĈ (IĈ =Ĉ: I). The deviatoric stress can then be expressed as where b = F F T denotes the left Cauchy-Green deformation tensor.
A more structurally accurate model is then examined, by augmenting the Neo-Hookean law with a fiber-dependent component (Humphrey 2002). This enhanced version which we will refer to as the Neo-fiber law is defined with respect to a fiber coordinate system, where the axes are aligned with the fiber f 0 , sheet s 0 and sheet-normal n 0 unit vectors. The strain energy function for the Neo-fiber law is defined as where a = 1 or 2. C 1 and C 2 are material parameters of the Neo-fiber law, corresponding to the fiber-dependent and isotropic terms, respectively. In this definition, IĈ f represents the first invariant ofĈ in the fiber direction, defined as The deviatoric Cauchy stress for the Neo-fiber law is then expressed as To better capture the exponential response of the myocardial tissue, our study is then extended to the structurally based orthotropic law by Holzapfel and Ogden (2009). The deviatoric form of the strain energy function for a 3-dimensional body is defined as In this definition IĈ s =Ĉ: s 0 ⊗ s 0 and IĈ f s =Ĉ: f 0 ⊗ s 0 denote invariants associated with the sheet and cross-fiber directions.
In what follows, we use a reduced version of the Holzapfel-Ogden law, where a s and a f s are set to zero and the exponents b and b f are kept constant, as we are restricting this study to constitutive laws with a small number of parameters to allow for better identifiability. A similar formulation has also been applied in (Caruel et al. 2014) in 0D and 1D models, demonstrating its ability to fit experimental data and thus reproduce physiological cardiac behavior. The values of the exponents (b = 5 and b f = 5) are chosen so that the model is able to provide a physiological EDPVR (see Sect. 3.3 and Fig. 12), as we note that for several combinations of b and b f a physiological EDPVR could not be produced for any values of the scaling parameters α and α f . Nevertheless, there is a range of values that would be appropriate for this choice, as we can assume interdependence between the exponents and scaling constants similar to that of the Guccione law (Xi et al. 2011a). The added value of this formulation over the Guccione law is that due to its structure as a sum of individual exponential terms, it can be reduced into a form with more than one uncoupled parameters. For this reduced version, the deviatoric Cauchy stress is derived as follows: Finally, we examine the well-known transversely isotropic exponential law by Guccione et al. The strain energy function Ψ is defined with respect to a fiber-oriented Green-Lagrange strain tensor E F where the rotation tensor Q is defined as Q = [ f 0 , s 0 , n 0 ]. The strain energy function is then defined as where a is a matrix of constants describing the degree of anisotropy in each component: The Cauchy stress tensor can then be expressed as

Parameter estimation
In patient-specific mechanics simulations, we often wish to tune or parameterize our models based on measurement data (observations). Supposing we have N parameters which govern the material behavior, it is a common approach to try and parameterize based on objective function minimization. For example, we aim to find a set of N parameters θ min which satisfies, for an objective function J θ , where P ⊆ R N is a subset of vectors of real numbers which constitute the admissible parameters for the problem. The behavior of the model, the observations (or data) over which the behavior is considered, and the objective function itself all play an important role in the behavior of the minimization problem and uniqueness of the minimizer. This is particularly important in clinical contexts where the obtained set of parameters θ min is used to, in some sense, provide an indicator of the health/state of the myocardium. In this section, we examine how these factors-the model, observations and objective-can influence the identifiability of θ min .

Model identifiability
To understand the behavior of the minimization problem, we first aim to better understand the behavior of the model and its parameters. Suppose we consider N s loading conditions imposed on our model (shown in Eq. 1). In this case, we can write the total model problem using the operator F s where we sum each quasi-static equilibrium state defined in Eq. 1, In this notation {t 1 , . . . t N s } denotes the set of N s loading conditions (boundary tractions) and denote the set of state solutions and test functions for each load state k. We can also compose the solution spaces set- where spaces W d 0 and W p correspond to displacement and pressure variables, respectively.
Using this notation, our quasi-static mechanics problem is (given a set of parameters θ ): find an X θ ∈ W s 0 such that, Here X θ constitutes the state solution composed of displacements and pressures at each of the N s loading conditions. We can further collect all solutions to (P1) and construct a space of solutions V ⊂ W s 0 , i.e., From the definition of V, we observe that we can identify pairings between a subset of P V ⊆ P and the space of state solutions V. These pairings, in general, have no well-defined properties. Indeed, V and P V may be empty. Here we suppose that (P1) induces a surjective mapping on the parameter space P to the space of state solutions V (see Fig. 3), i.e., This assumption is equivalent to assuming that there exists a unique X θ ∈ V for every θ ∈ P. If ϕ: P → V, it implies that for any θ ∈ P there is an X θ ∈ V. In other words, it implies there exists a (X θ , θ ) satisfying (P1). Moreover, nonuniqueness of solutions would imply that for some θ ∈ P, the best we could do to write the mapping ϕ is to write This possibility is precluded, however, by the surjective assumption (implying uniqueness). Though it is not proven for general cardiac mechanics boundary value problems, it is often assumed in these applications that, given the admissible Fig. 3 Schematic representation of the objective function J over the solution space W 0 for a given parameter set θ . The bijectivity of mapping ϕ between parameter space (P) and space of solutions (V) ensures practical identifiability set of loading states and parameters, the solution X θ to (P1) exists and is unique.
However, this condition is insufficient to guarantee unique parameter identifiability in Eq. 11. A stronger condition, which we show can ensure unique parameter identifiability, occurs when the mapping ϕ is in fact bijective, i.e., there exists a ϕ −1 where In the case where ϕ is bijective, we can ensure that the transition from state to parameters is well defined.
Two important determinants of the properties of ϕ stem from the behavior of the constitutive law itself and the set of loading states. The parameter dependence of the constitutive law-whether it be linear or nonlinear-can significantly influence a model's ability to uniquely identify parameters. It may also influence the range of loading states (and, as a result, deformations) which must occur to elucidate parameter dependencies. The final set of N s loading states {t 1 , . . . t N s } then defines this scope of deformations. For example, in the case of pure tension of the orthotropic Costa law, while fiber, sheet and sheet-normal parameters influence the end solution, the non-diagonal strain terms have no influence.
These considerations lead to a basic property required of a model referred to as structural identifiability (Raue et al. 2009(Raue et al. , 2011. A model is said to be structurally identifiable if there exists an arbitrarily large set of N T loading states {t 1 , . . . t N T } such that the mapping ϕ: P → V is bijective. As we show below (see Theorem 2 and Appendix 1), this property is easily proven for constitutive laws with linear parameter dependence, but becomes more complex when this dependence is nonlinear.
However, in most in vivo scenarios, model parameterization is limited to a given set of loading states which cannot be arbitrarily selected, leading to the concept of practical identifiability. A model is said to be practically identifiable if, for a given set of N s loading states {t 1 , . . . t N s }, the mapping ϕ: P → V is bijective. The key difference here is that limited observations comprise a set of loading states that ensure identifiability of all parameters, θ. In our case, these represent the in vivo states observed through the cardiac cycle which must be sufficient to yield practical identifiability of all parameters of the law.
For general nonlinear cardiac models, practical identifiability of passive parameters is difficult to prove a priori as it fundamentally depends on the structural identifiability of the model and the set of loading states and observations provided by the data. However, these considerations dictate the choice of model best suited for a given set of material deformations.
In general, bijectivity can be ensured by a coercivity assumption, i.e., Theorem 1 Suppose ϕ: P → V is a surjection (i.e., the solution to (P1) exists and is unique). Then if for any X ∈ V a pair θ 1 , θ 2 ∈ P, Proof The proof follows from contradiction. Suppose that θ 1 , θ 2 ∈ P both happen to satisfy (P1) when paired with the solution states X . Then by (P1), Dividing both sides by Y W u 0 and taking the absolute value and supremum, the coercivity assumption ensures, 0 ≥ θ 1 − θ 2 P or that the parameters θ 1 and θ 2 are, in fact, the same. Hence, any solution X has a single pair θ ∈ P.
A much simpler case occurs when the model depends linearly on N p parameters, the properties of ϕ are easier to decipher. In this case, the model may be written as: where σ n is the stress tensor (which may nonlinearly depend on u) scaled by the nth parameter. In this case, the bijectivity of ϕ can be ensured by guaranteeing that it is possible to construct N -constraints by using different Y 's in (P1). Due to the linear parameter dependence, the constraints may then be written as a matrix system, where the invertibility of the matrix ensures ϕ is bijective (see Theorem 2).
Proof The proof may be shown, again, by contradiction. Suppose there are two sets of nonidentical parameters θ , ψ ∈ P which result in the same state X . Then, by (P1), for any Y ∈ W u 0,Div × (W p ∩ 0). Hence, choosing {Y 1 , . . . Y N p }, Eq. 14 may be rewritten as Theorem 2 depends on a sufficient number of deformation states such that A gains linearly independent rows. We then rely on any test functions in W u 0,Div to further accentuate differences in material response, providing a flexible source from which to select constraints. Using this theorem, we can prove the structural identifiability for the Neo-Hookean, Neofiber and reduced Holzapfel-Ogden law which have a linear dependence on their parameters (see Appendix 1).

Objective function-based minimization
In data-based parameter estimation procedures, we often rely on some objective function to guide the choice of parameters. Since the parameters are not observed in most cases, we rely instead on comparing states with observations. In these cases, it is necessary that the objective function J : V → R obtains a unique minimum (a discerning objective).
Using 3D tagged data where the states are usually displacements, the natural choice of objective function is to use the L 2 norm over all states, i.e., where X d = {u 1 , . . . u N s } are observations on the displacements in the myocardium and we divide through by the overall scale of displacements so that J gives a percentage error.
In this case, ||| · ||| is a norm on W d 0 defined as, where (·, ·) is the L 2 −inner product on the reference domain Ω 0 . We then look to minimize the objective, finding X min ∈ V where As |||·||| acts as a norm on displacements in V (and a seminorm on the entire space), if the observed displacements in the state X d constitute a set of displacements X ∈ V, then J is automatically a discerning objective as the norm has a unique zero by definition (and is strictly nonnegative). This is, however, unlikely to occur in a real context due to two dominant factors: (1) data noise and resolution, (2) model fidelity. The introduction of noise, or degradation in data due to image resolution, introduces offsets which make the likelihood of X d being a member of V minimal. In addition, the fidelity of the model can strongly influence whether or not the model can capture the behavior observed in the data, making it possible that one or more than one minima exist.
Supposing that the model is a good representation of the tissue in vivo, we can then write X d = X + ε. In this case, if we assume that ε is, in fact, some random unbiased noise which satisfies then we observe that the noise does not bias our minima, but instead introduces a constant offset in J , i.e., The assumed relation in Eq. 17 is a reasonable assumption when the noise fluctuations occur over a small spatial scale compared with the change of the state solutions X near the minima.
The offset in Eq. 18 does not influence the minimization on V and, as a result, J remains a discerning objective. Obtaining a unique minimum in V is essential as even if a model is practically identifiable based on loading constraints, multiple minima for J guarantee multiple minima in parameter space. However, with a discerning objective, we then rely on the bijectivity (or practical identifiability) of ϕ so that the objective formed through composition, also obtains a unique minimum In practice, we can see that a discerning objective and a set of load states giving practical identifiability are suffi-cient conditions to ensure that the set of parameters θ min are uniquely identifiable.

Parameter coupling
Characterizing the behavior of the objective function J over the parameter space is critical for the performance of the parameter estimation process. The behavior of J around the minimum value (a distinct localized minimum instead of a wide valley) indicates a discerning objective function, which would allow data assimilation methods to retrieve the parameter estimate. Further, the landscape of J over the parameter space provides important information regarding the practical identifiability of the constitutive law, revealing the presence of a unique or multiple minima or possible interparameter coupling.
Coupling can also be deduced by the Hessian matrix of the objective function at the obtained minimum. Using the Taylor expansion of J around the obtained minimum θ min , Due to the gradient being zero at the minimum θ min , where H denotes the Hessian matrix. As can be deduced by this expression, the Hessian matrix can provide important information as it allows one to relate growth in J locally to local perturbations in the parameters. Further, to allow comparison between laws with varying parameters' scale, we use a scaled HessianH defined as where θ i and θ j correspond to the i-th and j-th components of θ min . Then using =˜ • θ min , Eq. 21 can be expressed as where now we are dealing with parameter percentages, which enables comparison between the different laws. The scaled HessianH can then characterize the sensitivity of J to the parameters and demonstrate possible inter-parameter coupling. Specifically, the minimum diagonal value ofH indicates that J is least sensitive to the corresponding parameter as a large error in the parameter can result in an insignificant change in J . Similarly, the minimum eigenvalue ofH indicates the parameter combination that J is least sensitive to. Accordingly, the ratio of the minimum diagonal value ofH over the minimum eigenvalue λ(H) demonstrates the degree of coupling between parameters. Specifically, large values indicate that there is a parameter combination whose possible error will cause a smaller change in J than error in each parameter separately, suggesting inter-parameter coupling. Similarly, coupling ratios close to 1 suggest that there is no significant coupling between parameters.

Comparison of practical identifiability using 3D tags
For each of the considered constitutive laws, the behavior of J over the parameter space is examined and the error between actual and estimated parameters is quantified. For each law we select a ground truth set of parameters which gives physiologically reasonable end-diastolic volume and pressure and generate synthetic 3D tags from an LV inflation simulation. The extracted myocardial motion applied on meshes is then used as our data and compared with simulations with varying parameter combinations to provide the landscape of the objective function and assess the error in the parameter estimates. Through this process we are able to characterize the practical identifiability of each law and assess its potential use in patient-specific applications.
All tests under consideration were implemented in CHeart-a multi-physics software tool based on (Nordsletten 2009; Nordsletten et al. 2010) and expanded by the CHeart team at KCL. The problems were solved on a Dell OPTIPLEX 990, quad-core (Intel Core TM i7-2600 CPU @ 3.40 GHz), on a quad-core (Intel 4th Generation Core TM i7-4790 CPU @ 3.60 GHz) and on an 2.1 GHz AMD Opteron TM Interlagos 32 processor.

J characterization of the Neo-Hookean law
We begin by investigating the practical identifiability of the Neo-Hookean law, by examining the behavior of J over a range of stiffness values. Due to the simple structure of the law and its linear parameter dependence, we expect that given some deformation, we should be able to get good identifiability characteristics. Specifically, we expect J to have a unique and distinct minimum, and we anticipate that the incorporation of unbiased noise should not affect the behavior of J and the estimated parameter, but only cause a shift toward larger J values.
Indeed, as illustrated in Fig. 4a, the objective function J has a unique and distinct minimum at the initial stiffness value (μ = 10 kPa). Further, the actual stiffness value is 10 × 10 and 11 × 11 simulations were performed for the initial and refined parameter sweeps, respectively, with an average computational time of 51.120 s retrieved even in the case of noisy data (5 and 20 % Gaussian noise), and the overall behavior of J remains the same, with just a small shift toward bigger values. These results suggest that using 3D tags we can uniquely and accurately estimate the stiffness value. The Neo-Hookean in silico test was extended to study the influence of the actual parameter value on the estimation process. Using the same inflating pressure, increased stiffness would result in smaller deformation that might be insufficient to allow for parameter estimation. We therefore performed four passive inflation simulations, where in each case the inflating pressure was adjusted to provide the same end-diastolic volume. We note that as can be observed in Eq. 1 (where the inflating pressure is introduced through traction as the product of pressure and deformed surface normal vector), due to the linear dependence of the law on the parameter the inflating pressure required to produce the same deformation was simply scaled by the ratio between the stiffness values. Figure 4b presents the behavior of the objective function J over scaled stiffness (μ over the initial stiffness for each case), showing consistent behavior for any initial stiffness value. This fact confirms practical identifiability of the Neo-Hookean law using 3D tags for any initial stiffness. Neo-fiber J over a the isotropic parameter C 2 and b fiber parameter C 1 , for a = 1 and a = 2. The actual parameters used are marked by the blue and red circles, for α = 1 and α = 2, respectively

J characterization of the Neo-fiber law
Figures 5 and 6 illustrate the behavior of the objective function J over the parameter space of the Neo-fiber law, for a = 1 and a = 2, respectively. As can be deduced from Fig. 5, the Neo-fiber law (a = 1) maintains the practical identifiability of the Neo-Hookean law (distinct minimum) and provides accurate parameter estimates. Even in the case of noisy data, the landscape of J remains similar and maintains a clear distinct minimum, with a small error (3.3 %) in the fiber parameter.
When the exponent, a, is increased to 2 (see Fig. 6) however, the practical identifiability of the Neo-fiber law is compromised (valley) and the error between actual and estimated parameters increases significantly (2 and 16 % for the isotropic and fiber parameters, respectively). Note that as no noise is added in this case, this error is created during the tagging process due to the combination of limited resolution of the tags and higher nonlinearity of the law.
Further, Fig. 7 examines the behavior of J over a range of values for each parameter separately. The steeper slope of J around the minimum value in the case of parameter C 2 suggests that the isotropic parameter has better identifiability than the fiber parameter. This issue is even more prominent in the a = 2 case due to the increased nonlinearity in the fiber dependence. Based on Theorem 2, the Neo-fiber law is structurally identifiable due to its linear parameter dependence, suggesting that the deterioration of the identifiability of the fiber parameter is due to insufficient deformation in the data. In fact, when the end-diastolic pressure is increased to 3 kPa instead of 1.5 kPa, the error in the fiber parameter decreases from 16 to 12 %.

J characterization of the reduced Holzapfel-Ogden law
The practical identifiability of the reduced Holzapfel-Ogden model (see definition in Eq. 8) was tested setting the values of the exponents (b = 5 and b f = 5) to provide a physiological EDPVR (see Sect. 3.3 and Fig. 12). Figure 8 presents the landscape of the objective function J over the parameter space, indicating a clear and unique minimum for the objective function. Even though the fiber parameter α f presents deteriorated identifiability characteristics compared with the isotropic parameter α (as also observed in the Neo-fiber case), we are still able to retrieve the parameter values with small relative errors (2.5 and 1 % for the isotropic and fiber parameters, respectively). Further, Fig. 9 which illustrates J over the parameter ratio α f /α for varying α, indicates the presence of a unique distinct minimum at the actual parameter. In this case, the behavior of J is also examined when only the enddiastolic frame is taken into account, illustrating that identifiability is preserved even when only one diastolic frame is used.

J characterization of the Guccione law
The practical identifiability of the transversely isotropic Guccione law was tested by choosing the parameters to fit an empirical end-diastolic pressure-volume relation (EDPVR), proposed by Klotz et al. (Klotz et al. 2006) (see Fig. 12). In order to assess the effect of noise in the data, 5 % Gaussian error was added in the simulation results, prior to in silico tagging. Table 1 presents the 5 parameter combinations with smallest J values, with and without 5 % Gaussian noise in the data. These combinations vary significantly, suggesting the presence of multiple minima. Indeed, the presence of 5 % noise in the data results in a different estimate of parameters compared with the non-noisy data. Note that this estimate which has a larger difference from the actual simulation parameters, compensates for the increase in C with a decrease in b f , suggesting an inter-parameter dependence. In order to further examine this issue, we reformulate the Guccione law as suggested by Xi et al. (2011a): where parameter α denotes the sum of b f , b t , b f s . We performed parameter sweeps over parameters C and α while keeping the ratios between b f , b t , b f s and α constant (r 1 = 0.5, r 2 = 0.3, r 3 = 0.2). Figure 10 illustrates the behavior of the objective function J over a range of the parameters C and α. The exponential shape of the blue valley representing model parameters resulting in small J values, verifies coupling between C and α as previously reported by Xi et al. Xi et al. (2011a, 2013. The presence of inter-parameter dependence is also demonstrated in Table  2 which shows a significantly larger coupling ratio R (see Sect. 2.3.3) for the Guccione law. Coupling may significantly deteriorate the parameter estimation process, as any noise in the data is likely to result in a large error in the estimated parameters (in this case 5.6 % for both α and C). The coupling in the Guccione law therefore suggests that we cannot guarantee unique and reliable parameter estimates using 3D tags.

Effect of noise in the 3D tagged data
The practical identifiability of the constitutive laws considered may be significantly compromised by the presence of noise in the data. In order to assess this effect, we have considered noisy data, where unbiased Gaussian noise was added in the simulation results prior to in silico tagging. Due to the limited resolution of the data, the addition of noise is expected to deteriorate the behavior of J , especially for parameters with very low sensitivity. The presence of noise resulted in increased J values (see Figs. 4a,5) and larger errors in the parameter estimates as indicated in Tables 1 and 3. However, the landscape of J was not significantly altered due to the uniform noise used, as indicated by the representative case of Neo-fiber a = 1 in Fig.  5. Nonetheless, unbiased noise caused a minor change to the topology of the objective function in parameter space as can be deduced by the increase in the coupling ratio in Table  2.

Comparison of models' fidelity
Keeping in mind that the choice of an appropriate constitutive law should balance between parameter identifiability and model fidelity, the constitutive laws described in Sect. 2.2.2 are tested with respect to their ability to represent physiological cardiac deformation. Further, the behavior of the objective function for any constitutive law is also influenced by model fidelity (see Sect. 2.3.2), as a model which cannot provide a good approximation to data can lead to unreliable parameters. We note that even though model fidelity does not affect the in silico tests in Sect. 3.1, where the same constitutive law is used for generated data and simulations, it is an issue for in vivo cases where the model should represent cardiac deformation.
For the purposes of this comparison, Guccione deformations were considered as the ground truth for physiological cardiac deformation. In order to cover a range of possible deformation modes during the cardiac cycle, 18 widely varying parameter combinations were used. For each of these combinations, a parameter sweep was performed for each law to provide the parameter that minimizes the difference between simulations and the ground truth cardiac deformation. The comparison was performed on simulations with the same end-diastolic volume through volume-prescribed loading to reduce the parameter space by 1. Figure 11 compares the minimum, maximum and average errors between the various laws considered and ground truth deformation. The Neo-Hookean law exhibits a significant error compared with ground truth deformation. This is mainly due to the inability of the Neo-Hookean law to produce adequate elongation and twist, which are important characteristics in cardiac deformation. On the other hand, the added fiber dependence in the Neo-fiber law allows a more accurate approximation to physiological cardiac motion, as the lower average and maximum errors indicate that the Neo-fiber can on average reproduce most of the deforma-  First row corresponds to the gradient with respect to the first parameter for each law. Coupling ratios R are also presented, for data with and without 20 % unbiased noise tion modes considered. The approximation to cardiac motion is further improved as the exponent α increases. Finally, the reduced form of the Holzapfel-Ogden law presents the smallest average and minimum error between the deformations produced using a given constitutive law and the Guccione law, confirming that the values used for the exponents b and b f are appropriate and allow for physiological cardiac deformation. This is also illustrated in the meshes in Fig.  11 which present the maximum difference from Guccione combinations. For the Neo-fiber and reduced Holzapfel-Ogden models, even the maximum difference from the ground truth is still quite small as indicated by the close match between model and ground truth cardiac deformation.  Fig. 12 Typical end-diastolic pressure-volume curves for the constitutive laws considered and ground truth Klotz curve. While the reduced Holzapfel-Ogden and Guccione laws are able to reproduce the Klotz curve, the Neo-Hookean and Neo-fiber laws cannot produce a physiological pressure-volume response

Comparison of models' EDPVRs
End-diastolic pressure-volume relation (EDPVR) is an important determinant of cardiac function; therefore, the considered constitutive laws were compared with respect to their ability to reproduce a physiological EDPVR. For this comparison, the empirical EDPVR proposed by Klotz et al. (2006) was chosen as the ground truth physiological EDPVR, as it is considered capable to represent healthy and diseased cases. The EDPVR is derived from a single set of end-diastolic pressure (EDP) and volume (EDV) measurements, which for the in silico tests were chosen as EDP = 11 mmHg, EDV = 140 ml. Figure 12 illustrates typical EDPVRs for the various constitutive laws considered. As indicated by the curve, the Neo-Hookean and Neo-fiber laws were not able to reproduce a physiological EDPVR, even though case a = 2 gives a better approximation for the Neo-fiber law. On the contrary, the exponential Guccione and the reduced Holzapfel-Ogden laws were able to provide a physiological EDPVR as indicated by the close match to the Klotz curve.

In vivo behavior of J
The in silico tests performed in Sects. 3.1, 3.2 and 3.3 suggest that the reduced Holzapfel-Ogden law combines practical identifiability with good representation of cardiac deformation and EDPVR. Therefore, the reduced Holzapfel-Ogden is suitable for patient-specific applications as it is an accurate cardiac model with reliable-thus potentially clinically important parameters.
In this section the reduced Holzapfel-Ogden law is employed in an in vivo case of a healthy volunteer. The behavior of the objective function J is examined in this setting as well, to allow for the characterization of practical identifiability in a real-world scenario. Within this setting, we can also assess the effect of model fidelity on the parameter estimation process, as the simulations are now compared with actual cardiac deformation data.
The extracted 3D displacement field was compared with passive inflation simulation results to provide parameter estimates for the reduced Holzapfel-Ogden law. As the LV pressure trace was not part of the available data, we were not able to obtain unique values for each parameter separately. However, taking advantage of the linear dependence of the law on the parameters and accordingly their proportionality to inflating pressure, we were able to uniquely estimate the ratio α f /α, which is independent of the inflating pressure. Note that if the end-diastolic pressure is available, we can retrieve the actual values of α f and α by multiplying by the ratio between the actual pressure values and the values used in the estimation process. Figure 13 illustrates the behavior of J over a range of ratios α f /α, where the fiber parameter α f was kept constant. Even though in this in vivo case we cannot assess the error in the parameter estimates, we can still infer that the practical identifiability of the reduced Holzapfel-Ogden law observed in in silico tests 3.1.3 is maintained when actual 3D tags are used, based on the distinct and unique minimum. The identifiability of the law using all or only the final diastolic frames is examined as well, indicating that J presents distinct and unique minima in both cases, with a 12 % difference between the two estimated values. The similarity in behavior of J over the parameter ratio between in silico and in vivo tests (see Figs. 9, 13) confirms model validity for the reduced Holzapfel-Ogden law. Finally, it offers a validation of our in silico testing protocol, thus allowing for the conclusions for the considered constitutive laws to be extrapolated to parameter estimation applications with real 3D tags.

Study limitations
A relatively coarse, lower order mesh was used for the in silico example. This choice was based on the small computational time per simulation, which allowed for the large number of simulations performed to provide the landscape of the objective function over the parameter space. Nevertheless, as illustrated in Appendix 2, the mesh resolution used was sufficient to examine parameter identifiability characteristics in our models, for both the in silico and in vivo tests.
Further, only one objective function has been considered, even though other objective functions might be able to elucidate other characteristics for the behavior of the various laws. Other objective functions are not considered here, as such a study would require proving that a potential objective function is discerning. If non-discerning, proving that unique parameter identification is achievable becomes challenging as it then must rely on the observations not exposing non-uniqueness of the objective function. Nevertheless, the chosen J uses an L 2 norm on the displacements, which is generally considered a robust criterion and should be able to provide adequate information and accurately describe the identifiability characteristics of each law.
Ten diastolic frames were used as observations in the in silico test. This choice was based on the number of diastolic frames in the available 3D tagged MRI data (e.g., where the cavity volume is increasing). However, the number of diastolic image frames used in other studies is variable with authors considering all or part of diastole. This variability is due to assumptions on residual active tension, the presence of which is confirmed by decreasing cavity pressures even after the opening of the mitral valve (Pasipoularides et al. 1986). As a result, early diastolic frames do not contain purely passive tissue behavior, but also contain residual active stress. We observed that parameter identifiability for the reduced Holzapfel-Ogden model was preserved in both in silico and in vivo scenarios, using only the last diastolic frame (see Figs. 9 and 13). J presented unique and distinct minima for all cases, with a 12 % difference between the two estimated values in the in vivo case. However, this robustness is likely due to the single parameter (α f /α ratio) estimated over the entire left ventricle. Incorporation of spatially varying parameters would increase sensitivity to noise and would likely require additional passive diastolic frames to ensure identifiability.
Even though the effect of unbiased noise is examined, several aspects of the process or the data that may have a biased influence are not considered. For instance, the resolution and number of tagged lines in the data along with the tracking algorithm used are likely to incorporate consistent error in the parameter estimation process. The boundary conditions used in the simulations are also likely to influence the identifiability and model fidelity results. Understanding these attributes is important for patient-specific personalization; therefore, further work is required to clarify the influence of these factors on the landscape of the objective function and the estimation process in general.
The reference configuration is assumed to be known and correspond to a specific frame of diastole in the in vivo example. Even though this is a common approach, we have verified that the identifiability of model parameters is not sensitive to the frame used as a reference configuration for the in silico case (see Appendix 3). Inter-estingly, the parameter ratio is relatively consistent (20 % maximum error between actual and estimated ratio) irrespective of the reference domain. We note, however, that this might be an artifact of the idealized geometry used or other simplifications inherent in the model produced data.
Finally, only one in vivo case is considered; thus, further work is needed to examine effects in vivo, such as noisy or low-resolution data, or diseased cases where the cardiac deformation is likely to differ significantly. However, the in silico tests performed provide a standard, giving the "best case scenario," which can be anticipated when using real data. Further, the in vivo and in silico model fidelity of the reduced Holzapfel-Ogden law encourage the application of the law to diseased cases as well.

Conclusions
In this paper, we examine the practical identifiability and model fidelity of a range of cardiac constitutive laws using 3D tagged MRI as the available data. In order to investigate the practical identifiability of the laws considered and examine the potential of using 3D tags in parameter estimation applications, we generate synthetic 3D tags directly from simulation results and assess the behavior of the objective function over the parameter space through parameter sweeps. The laws considered are also compared with respect to their ability to represent physiological cardiac motion and EDPVR, elucidating the primary components that should guide the choice of an appropriate cardiac constitutive law-namely reliable parameters and adequate representation of cardiac deformation and function.
Our results verify the reported coupling of the transversely isotropic Guccione law, suggesting the need for a law with better identifiability characteristics that would allow for reliable parameter estimates. The Neo-Hookean law is shown to have good identifiability characteristics, due to linear parameter dependence. The stiffness parameter is identifiable, provided adequate deformation is present in the available data. However, due to its isotropy, Neo-Hookean deformation misses key characteristic deformation modes, mainly long-axis elongation and twist. Further, it cannot reproduce physiological pressure-volume response.
Building on the Neo-Hookean model, the Neo-fiber law maintains the good identifiability characteristics, while reproducing physiological cardiac deformation. Both parameters are identifiable, even though sufficient deformation is required to allow identifiability of the fiber parameter due to the structure of the constitutive law. However, using the Neo-fiber law leads to inaccurate pressure-volume response, which cannot match the empirical Klotz curve.
The reduced Holzapfel-Ogden law, however, combines all important attributes considered, offering a balance between practical identifiability and adequate representation of cardiac deformation and EDPVR. Its use in an in vivo case where good identifiability characteristics are maintained supports the conclusion that the reduced Holzapfel-Ogden law offers a sensible choice in patient-specific applications with 3D tagged data.
For the case of the Neohookean law, due to its isotropy, pure tension in any of the three directions is sufficient to ensure structural identifiability. Specifically, for elongation in the X direction the deformation gradient F and left Cauchy-Green deformation tensor can be expressed as vivo test could be justified by the fact the finer mesh was also used to propagate the extracted myocardial motion, and its cavity volume was used to drive the simulation.

Appendix 3: Effect of reference configuration on J behavior
This section investigates the effect of choosing a different diastolic frame as a reference configuration on the behavior of the objective function. Specifically, we have used an in silico test to examine the effect of the reference configuration on the identifiability of the parameter ratio of the reduced Holzapfel-Ogden law, in a volume-prescribed diastolic filling simulation. Five different reference configurations were considered for the simulations performed over parameter sweeps: the initial reference mesh and meshes corresponding to loading steps 10, 20, 30 and 40 of the simulation used for the creation of the synthetic tags. In each case, the number of loading steps  Fig. 17 J over the parameter ratio α f /α for the reduced Holzapfel-Ogden law. Five different reference configurations are used, which correspond to different diastolic phases and thus observations was adjusted, so that each observation would correspond to 5 loading steps. Figure 17 compares the behavior of the objective function J when different diastolic phases are used as the reference configuration. Based on the similar behavior of the curves, we can deduce that the identifiability of the parameter ratio is not affected by the assumed reference configuration.