Investigating the reference domain influence in personalised models of cardiac mechanics

A major concern in personalised models of heart mechanics is the unknown zero-pressure domain, a prerequisite for accurately predicting cardiac biomechanics. As the reference configuration cannot be captured by clinical data, studies often employ in-vivo frames which are unlikely to correspond to unloaded geometries. Alternatively, zero-pressure domain is approximated through inverse methodologies, which, however, entail assumptions pertaining to boundary conditions and material parameters. Both approaches are likely to introduce biases in estimated biomechanical properties; nevertheless, quantification of these effects is unattainable without ground-truth data. In this work, we assess the unloaded state influence on model-derived biomechanics, by employing an in-silico modelling framework relying on experimental data on porcine hearts. In-vivo images are used for model personalisation, while in-situ experiments provide a reliable approximation of the reference domain, creating a unique opportunity for a validation study. Personalised whole-cycle cardiac models are developed which employ different reference domains (image-derived, inversely estimated) and are compared against ground-truth model outcomes. Simulations are conducted with varying boundary conditions, to investigate the effect of data-derived constraints on model accuracy. Attention is given to modelling the influence of the ribcage on the epicardium, due to its close proximity to the heart in the porcine anatomy. Our results find merit in both approaches for dealing with the unknown reference domain, but also demonstrate differences in estimated biomechanical quantities such as material parameters, strains and stresses. Notably, they highlight the importance of a boundary condition accounting for the constraining influence of the ribcage, in forward and inverse biomechanical models.


Introduction
Cardiac biomechanics have proven particularly useful for improving our understanding of heart function in health and disease, providing a quantitative framework for studying cardiac pathophysiology. Myocardial wall stress, for instance, has long been hypothesised as a contributing factor to tissue growth and remodelling, reported to rise in heart failure patients or after myocardial infarction (Holmes et al. 2005). Hypokinetic hearts with heterogeneous deformation have been associated with cardiomyopathies (Chuang et al. 2010), rendering myocardial strains a valuable biomarker. Material properties such as stiffness and contractility are intrinsically linked to the tissue structural integrity and present abnormalities in ischemic conditions and diastolic heart failure (Zhang et al. 2010). Accordingly, evaluating biomechanical properties on a patientspecific basis holds significant potential to improve clinical diagnosis and identify individualised treatments.
Personalised cardiac mechanics models, combining mathematical modelling with clinical data, are becoming increasingly popular as a tool for noninvasively assessing biomechanical quantities. Latest finite element (FE) models incorporate patient-specific geometries , microstructure (Stoeck et al. 2018) and boundary conditions , and employ physiologically accurate material description (Holzapfel and Ogden 2009), providing enhanced accuracy in stress quantification. Model personalisation also enables the estimation of tissue stiffness and contractility, by fitting model outcomes to data-derived cavity volumes and pressures (Asner et al. 2015b;Costandi et al. 2006), threedimensional tissue displacements and strains (Finsberg et al. 2018). Substantial progress has also been achieved in the parameter estimation process, through studies focusing on acquiring physiologically relevant parameter estimates (Hadjicharalambous et al. 2014a;Nasopoulou et al. 2017) or proposing efficient schemes for parametrisation and heterogeneous parameter distributions (Xi et al. 2011).
Despite these advances considerable model uncertainties still remain, with one important confounder being the unknown zero-pressure domain. Due to the physiological blood pressure in cardiac cavities and the generation of active tension during systole, the heart is continuously under load throughout the cardiac cycle, making it impossible to capture the unloaded domain by routinely acquired clinical data. A typical approach is to instead generate the reference domain based on an in-vivo image frame. Common choices are the geometries at end systole , early diastole (Xi et al. 2013) and diastasis (Pfaller et al. 2019), as they correspond to cardiac phases with low endocardial pressures. These geometries are generally not stress-free but are instead associated with unknown in-vivo stress and strain fields, inevitably introducing inaccuracies in model-derived deformation, stress and material parameters. Other studies have employed inverse methodologies to acquire an approximation of the zero-pressure domain (Nordsletten et al. 2011;Nikou et al. 2016;Vavourakis et al. 2011) based on a deformed state and known loading conditions. Amongst them, iterative backward incremental schemes (Krishnamurthy et al. 2013) and inverse elastostatics approaches which reformulate the finite elasticity equations have drawn particular attention (Vavourakis et al. 2016;Rajagopal 2007;Peirlinck et al. 2018). Inverse approaches in general suffer from strong coupling between the estimated reference domain and the passive parameters, forcing studies to employ literature values or elaborate schemes for joint estimation (Nikou et al. 2016;Wang et al. 2020;Costandi et al. 2006). Furthermore, the uncertainties surrounding the unloaded domain extend to the boundary conditions that should be employed for its estimation, posing an additional challenge for inverse schemes.
Despite their shortcomings, these methods have shed light on the challenges associated with the selection of reference domain and on its influence on model estimates of cardiac biomechanics. For instance, usage of an image frame as reference domain or employment of an inverse approach has led to distinct differences in stresses (Nikou et al. 2016), strains (Behdadfar et al. 2017;Kallhovd et al. 2019) and material parameters (Finsberg et al. 2018;Nikou et al. 2016), while the choice of specific image frames also impacts parameter values Xi et al. 2013). While biomechanical model predictions are evidently impacted by the selection of zero-pressure domain, quantification of this influence is generally not possible due to the lack of suitable ground-truth data. In fact, very limited experimental data exist to describe the unloaded geometry and its relation to in-vivo known conditions. An exception is the study by Klotz et al. (2006), which proposes an empirical relation for estimating the unloaded state's cavity volume as a function of the end-diastolic pressure and volume, shown to correlate well with measurements in healthy and diseased hearts. In the absence of experimental data on the reference geometry and its position relative to in-vivo states, evaluating the accuracy of existing inverse schemes or determining appropriate boundary conditions would be unattainable. In parallel, without ground-truth data on the zero-pressure configuration it is infeasible to quantify the effect of selecting incorrect reference states on model-derived outcomes.
In this work, we aim to enhance our understanding of the importance of the unloaded configuration by employing a modelling framework relying on experimental data. Imaging data from in-vivo and in-situ experiments on a porcine heart are used to offer detailed information on the in-vivo behaviour of the individual heart, and more importantly, to provide a reliable approximation of the unloaded myocardial state. These novel datasets are used to personalise cardiac mechanics models using our previously developed data-model integration pipeline (Asner et al. 2015a;Hadjicharalambous et al. 2017), creating a unique opportunity for a validation study.
In particular, we develop three whole-cycle cardiac models, employing either the in-situ image, an in-vivo frame or an inversely estimated geometry as their zero-pressure domains. A thorough comparative analysis between the three models is performed, allowing a quantitative assessment of the effect of the unloaded domain on fundamental biomechanical quantities. The in-situ image of the reference geometry also enables the evaluation of existing inverse methodologies, identifying model aspects that should be carefully considered when aiming to develop high-fidelity patient-specific applications.

Materials and methods
This section describes the specifics of the experimental procedure and model development followed to examine the reference domain influence. Exploiting data from innovative invivo and in-situ experiments, three modelling scenarios were considered to elucidate the differences in model-derived indices caused by the deployment of different geometries as the model unloaded domain: (1) the in-situ model, which employs the unloaded geometry imaged in the in-situ experiment as the reference domain, (2) the end-systolic model which uses the end-systolic frame of the in-vivo image as the reference domain, and finally (3) the inverse model, which employs an inverse finite deformation approach to acquire an estimate for the unloaded geometry.

Experimental procedure
The data used in this study were acquired over prospective in-vivo and in-situ experiments on a porcine heart. During the initial in-vivo phase of the experiments, imaging data were acquired, while the animal was under anaesthesia. Subsequently, the heart of the pig was arrested through administrating an overdose of Pentobarbital. High doses of Pentobarbital are known to inhibit the heart's contractile function, causing substantial decreases in the systolic pressures and the maximal rate of change of cavity pressures (Jiang et al. 2011). Therefore, the arrested heart was at a purely passive state during the in-situ phase of the experiment. Importantly, as there was no blood flowing through the left-ventricular cavity, the blood pressure load on the endocardial wall was near zero. With the blood remaining in the cavity (i.e. not drained), there was no risk of the heart wall collapsing. All in-situ images were acquired within 10 minutes of arresting the heart, to ensure that the myocardial tissue was still alive. The precautions taken, along with ensuring the myocardial tissue was at a passive state with a near-zero pressure in the left-ventricular cavity, enabled a reasonable approximation to the animal's myocardial unloaded domain. Accordingly, all images recorded at this state provided views of the subject-specific zero-pressure domain (in-situ reference images).
The in-vivo and in-situ experiments were performed on a healthy porcine heart (Swiss Edelschwein breed, 50 kg, female), during a study approved by the local Committee for Experimental Animal Research (Cantonal Veterinary Office Zurich, Switzerland) under the License numbers ZH219/2016. The animal was anaesthetised by intramuscular injection of Ketamine (15 mg/kg bodyweight (BW)), Azaperone (2 mg/kg BW) and Atropine sulfate (0.05 mg/kg BW). Analgesia was ensured by application of Buprenorphine (0.01 mg/kg BW). Throughout the experiment anaesthesia was maintained with isoflurane (2-3%) by positive pressure ventilation with an inspired oxygen fraction (FiO2) of 100%, a tidal volume of 10 mL/kg and a positive end expiratory pressure (PEEP) of 5 cmH2O.
Imaging was performed on a clinical 1.5T Philips Achieva System (Philips Healthcare, Best, The Netherlands) with 5-channel dedicated cardiac receiver coil. The protocol consisted of CINE Magnetic Resonance Imaging (MRI) in four-(4CH), three-(3CH) and two (2CH)-chamber long-axis (LA) views and multi-slice two-chamber short-axis (SAX) view covering the entire heart. Imaging parameters were as follows: in-plane resolution: 1.5 × 1.5 mm 2 , slice thickness 8 mm, number of heart phases 25 (SAX)/50 (LA), number of slices 13 (SAX)/1 (LA). ECG electrodes (Quatrode patch, In-Vivo Corporation, Gainesville, USA) were placed onto the chest between the front limbs for retrospective cardiac gating.
During the in-vivo experiment, the pressure of the left ventricle (LV) was recorded at 2 ms intervals, using a pigtail catheter inserted into the femoral artery and advanced to the LV cavity. For the in-situ phase, the heart was arrested by administration of an overdose of Pentobarbital (75 mg/Kg) and the imaging protocol was repeated.

Image registration
A critical step for model personalisation was the spatiotemporal registration of available images, to minimise 1 3 misalignment due to changes in pig's position or breathing cycle. Rigid registration between the SAX and LA images, for both the in-vivo and in-situ images, was performed using MIRTK 1 . Registration of the static in-situ image to the dynamic in-vivo images required a more elaborate procedure, described in "Appendix Registration of reference domain to in-vivo images".

Construction of personalised FE meshes and fibre architecture
LV meshes were constructed based on manual segmentations of the in-situ image and the end-systolic (ES) and end-diastolic (ED) frames of the in-vivo images. For each case, segmentations of all LA and SAX images created with ITK-SNAP 2 (Yushkevich et al. 2006), were combined into a single mask image, which was further refined and smoothed.
The smoothed masks for the three models were subsequently truncated at corresponding planes below the mitral valve, to avoid the need for modelling the basal tissue which is often poorly resolved in standard MRI. Additionally, wall volumes were almost identical for the three truncated masks (less than 5% difference), in accordance with the myocardial incompressibility assumption. Volumetric tetrahedral meshes ( Fig. 1, Table 1) were then created based on the mask images, using SimModeler. 3 To represent the anisotropic behaviour of myocardial tissue, a fibre architecture was constructed for each model. In the absence of personalised data providing insight into the microstructure on a patient-specific basis (Stoeck et al. 2014), a rule-based fibre distribution was employed, consisting of circumferentially symmetric fibres, with the angle between fibres and circumferential direction varying linearly between 60 • at the endocardium and −60 • at the epicardium (Streeter et al. 1969).

Data-derived displacement and pressure-volume loop
Whole-cycle simulations for the three models were personalised using data-derived displacements and pressure-volume loop, which were incorporated into the models through suitable boundary conditions (Sect. 2.3.3). Displacements and cavity volumes were derived from the three reconstructed geometries (end-systolic and and in-situ reference (right) geometries, with the boundaries of the personalised tetrahedral meshes overlaid, in short-axis (top) and 4CH (middle) views. Bottom row demonstrates the spatial field R for each mesh, which describes the relative influence of the ribcage on the epicardial surface based on the distance between them ( R = 1 where the epicardial surface is almost in contact with the ribcage (in red) and R = 0 in regions not affected by the ribcage (in blue)) LV cavity pressure measurements were available throughout the duration of the experimental procedure. However, these measurements had consistently high ED and low ES values compared to typical literature values. This suggested the presence of a systematic measurement error, potentially caused by an inaccurate positioning of the pressure catheter relative to the heart or by the anaesthetic administered. To avoid discarding the entire dataset, the pressure curve was scaled to achieve physiological porcine end-diastolic and end-systolic LV pressures (Kolipaka et al. 2010), thus retaining personalised temporal information and relative pressures through the cycle. Combined with the volume trace, a subject-specific pressure-volume loop was generated, enabling full-cycle simulations.

Kinematics and governing equations
An indispensable part of cardiac mechanics models is the reference domain, the state of the myocardium under zero load, which is of particular interest in this work. Let 0 ⊂ ℝ 3 denote the unloaded domain with coordinates X ∈ 0 , and let the domain's boundary 0 be subdivided into the basal b 0 , epicardial e 0 and endocardial 0 surfaces. The heart deforms continuously during the cardiac cycle, with its deformed domain (t) ⊂ ℝ 3 characterised by the physical coordinates x ∈ (t) , for each time t > 0 , expressed as x = X + u , u being the tissue displacement. The deformation gradient, defined as F = ∇ X u + I , is used to relate points in the deformed state to their reference configuration, while its determinant J is used to enforce the incompressibility constraint ( J = 1 ). Other kinematic quantities of interest are the right Cauchy-Green tensor C = F T F and its material invariants I C = C ∶ I and I C f = f 0 ⋅ (Cf 0 ) (Holzapfel and Ogden 2009), with f 0 denoting the fibre direction throughout the domain.
For the models considered the total potential energy of the heart is given as a sum of its internal and external energies (Asner et al. 2015a), (u, p, ) = int (u, p) + ext (u, ), with int determined by material properties and ext depending on external forces introduced through boundary conditions. Here p and denote the hydrostatic pressure and boundary Lagrange multipliers, respectively (see Sect. 2.3.3). Based on the principle of stationary potential energy (Bonet and Wood 2008), the primary variables ( u, p, ) are determined as the saddle-point solution of the potential energy: minimising the internal energy of the myocardium, while ensuring adherence to the incompressibility and boundary constraints. Here, U × W × Λ denote suitable Sobolev function spaces ).

Constitutive laws
In the models considered, the myocardium is treated as an incompressible, hyperelastic material, with its internal energy described as a sum of passive ( W p ) and active ( W a ) strain energy functions: Material laws were selected to combine physiologically accurate myocardial function with good parameter identifiability characteristics (Hadjicharalambous et al. 2014a;Asner et al. 2015a). In particular, a reduced form of the Holzapfel-Ogden law (Holzapfel and Ogden 2009), shown to provide unique parameter estimates, physiological cardiac deformation and end-diastolic pressure-volume relation (EDPVR) (Hadjicharalambous et al. 2014a), was employed: Here a, b, a f and b f are material constants: a characterises the bulk tissue stiffness, a f the stiffness along the fibre direction, while b and b f modulate the form of the nonlinear tissue response. The exponents are kept constant throughout simulations ( b = 5 , b f = 5 ) as they cannot be uniquely identified based on the physiological deformations observed during a cardiac cycle (Hadjicharalambous et al. 2014a). The specific values selected for the exponents allow the model to reproduce physiologically relevant cardiac function, while ensuring unique parametrisation for a and a f reproduce physiologically relevant cardiac function (Hadjicharalambous et al. 2014a).
The active behaviour was described with a law entailing a single, time-dependent, contractility parameter (t) applied across the domain (Kerckhoffs et al. 2003;Asner et al. 2017): The specific form of iso was selected to enable isotropic activation near the apex, to avoid known singularities in the apical region , while the activation occurs along the fibre direction in most of the domain: Here, D(X) is a distance metric which is one near the apex X a and rapidly decays to zero away from it: where L la is the long-axis length. The activation function g is selected to account for length dependence, with l 0 = 0.8 the lower bound for the compressive strain of myocytes (Kerckhoffs et al. 2003):

Boundary conditions
Selecting boundary conditions is a critical step in model development, as they provide a means of representing the individual loads affecting the heart's surface and can incorporate data-derived motion into the model. Boundary conditions were introduced through boundary energies, using Lagrange multipliers: Whole-cycle simulations were driven using cavity volumes, whereby the mesh lumen volume V(t) was set equal to the data-derived volume trace V d (t) at each timestep: In this setting, is equivalent to the cavity pressure (Asner et al. 2015b (the weak form of this formulation can be shown to be equivalent to the weak form of a problem employing the commonly used pressure boundary condition).
For the basal boundary, a relaxed boundary condition was used to avoid strict adherence to noisy data and limit nonphysiological stress peaks associated with Dirichlet conditions ): .
is a penalty matrix, which forces displacements along the base normal vector n b to be equivalent to data, while allowing a relaxed adherence through the base plane. The penalty parameter b ( 0 < b ≪ 1 ) controls the degree of enforcement of the constraint, depending on data quality or numerical considerations.
The data-derived displacement on the basal boundary was computed by first constructing an affine mapping to transform the reference basal domain to its deformed configuration at ED (or ES), as extracted from the ED mesh (or ES mesh). In practice, the affine mapping A (REF to ED) transforms the reference base plane X b 0 to its best approximation of the ED base plane is the base boundary of the mesh employed as reference domain in each model.
In terms of epicardial boundary conditions, the zero external energy constraint ( e ext = 0 ) was considered following the majority of studies. However, the heart is not in isolation, but instead is surrounded by the pericardial sac and is in close proximity to the lungs, the ribcage and the diaphragm. In fact, due to the porcine anatomy, the ribcage is found particularly close to the heart, potentially affecting myocardial motion and stress fields. Therefore, additional constraints were also considered which account for the ribcage's influence on the LV, by restricting the displacement of the epicardial wall in a region adjacent to the ribcage: Here e is a Lagrange multiplier defined over the epicardial surface and R ∶ e 0 → [0, 1] is a spatial field signifying the influence of the ribcage on the epicardial surface based on their distance, as derived from the segmentations of SAX images (Fig. 1). R varies smoothly between 0 and 1, with 1 marking the epicardial regions most adjacent to the ribcage. Through R, the motion of the epicardial areas closer to the ribcage is restricted, i.e. u d = 0 . The penalty matrix follows the basal boundary condition ( K e = I − n r ⊗ n r ), enforcing a stronger adherence of the constraint along the normal vector ( n r ) of the rib-adjacent region. This constraint was employed for the end-systolic and inverse models.
In the in-situ model, however, constraints were also considered for the apical region, due to substantial variations in the positions of the reference and in-vivo domains (see "Appendix Registration of reference domain to in-vivo images"). The epicardial surface was further subdivided into an apical region ( e 0 = a 0 ∪ e�a 0 ), to enable individual constraints on the apical region: where K a = I − n a ⊗ n a , with n a the normal vector to the apical region. In this case, u d was computed as the displacement of the apex from the in-situ reference image to the ED (or ES) image. The scalars w e and w a were either 0 or 1, allowing for selective enforcement of the two epicardial constraints, leading to the following combinations of boundary conditions for the in-situ model: B: Only basal boundary condition while e ext = 0 ( w e = 0 , w a = 0), BA: Data-derived deformation prescribed on basal and apical regions ( w e = 0 , w a = 1), BR: Data-derived deformation prescribed on basal plane, restricted motion in ribcage-adjacent region ( w e = 1 , w a = 0), BAR: Data-derived deformation on basal and apical regions and restricted motion near the ribcage ( w e = 1 , w a = 1).

Parameter estimation
For the in-situ and end-systolic models, passive material parameters were estimated by minimising the error between the simulated ED state for each model and the in-vivo ED geometry, through volume-driven simulations. Three error metrics were considered, namely the root-mean-square error (RMSE), and the mean and maximum distance errors ( d mean and d max ). Error metrics were calculated with respect to nodal distances between the simulated ED mesh ( x sim ED ) for each model and the data-derived ED mesh ( x ED ): Here d i represents the distance of node i of the simulated ED mesh to the closest node in the data-derived ED mesh, and N n is the total number of nodes for the former mesh.
As previously described (Asner et al. 2015a), for the estimation of passive parameters a and a f in volume-driven simulations, we can exploit the linear dependence of the law on the parameters, to reduce the parameter space to only the parameter ratio a∕a f . In particular, scaling both passive parameters by an amount m would lead to the same scaling of the endocardial multiplier or cavity pressure, i.e. m × , without altering the simulated deformation ) (i.e. volume-driven simulations with the same parameter ratio a∕a f produce identical deformations).
Taking advantage of the proportionality of parameters to cavity pressure, we can uniquely estimate the parameter ratio-which is independent of the cavity pressure-and retrieve the absolute values of a and a f by scaling by the ratio between the data-derived ED pressure P ED and its simulated value sim : The ratio was estimated through parameter sweeps and considering the aforementioned error metrics, by keeping a sim f constant and varying a sim . The described estimation process was, however, not appropriate for the inverse model. Reference domain estimation is directly dependent on the material parameters, which, in turn, require an assumed unloaded domain for quantification. To circumvent this coupling, the reference domain and material parameters were jointly estimated using sweeps over passive parameters and inverse simulations. As the application of the cavity volume constraint in an inverse setting is not straightforward, a cavity pressure boundary condition was applied on the endocardial wall. Assuming no knowledge of the geometry of the reference domain-to properly evaluate inverse methodologies-joint estimation was performed using a cavity volume error metric V ref d . This error metric relied on the empirical estimate of the cavity volume of the unloaded domain (Klotz et al. 2006). Based on experiments on excised animals and human hearts, Klotz et al. proposed an expression for calculating the cavity volume of the unloaded domain, V Klotz , using a single pair of end-diastolic pressure P ED and volume V ED measurements: Here, V Klotz = 44.27 mL was calculated based on P ED = 10 mmHg and V ED = 81.98mL. The cavity volume metric used for the joint estimation of parameters and reference domain was then defined as the difference between the cavity volume of the simulated reference domain V ref (a,a f ) and its Klotz estimate: However, it was not possible to uniquely identify both passive parameters using only estimation by cavity volumes ("Appendix Investigation of passive parameters estimation for inverse model"). To circumvent this, a combination of existing schemes (Nikou et al. 2016;Xi et al. 2013) was used, whereby the end-systolic frame provided an initial guess for the zero-pressure geometry, and an initial estimate of passive parameters was acquired through forward simulations (i.e. these coincide with parameters estimated for the end-systolic model). Subsequently, inverse simulations with parameter sweeps were run, whereby a f was set to its initial estimate and a was allowed to vary. The new estimate of the reference domain INV and a was then obtained by minimising V ref d . By construction, the estimated combination of material parameters and unloaded domain should provide an almost identical match to the ED geometry in forward diastolic simulations, while ensuring a good agreement with two points in the EDPVR, (V Klotz , 0) and (V ED , P ED ) . The process for jointly estimating the zero-pressure domain and material parameters is illustrated in Fig. 2.
Finally, once passive characteristics were configured, whole cycle simulations for all models were run by prescribing data-derived cavity volumes. Throughout the cycle, active characteristics were determined at each timestep, by matching model cavity pressures to the pressure measurements described in Sect. 2.2.3. Specifically, for each timestep t, the simulated cavity volume V(t) was enforced to match the measured V d (t) , while the contractility parameter (t) in Eq. 4, which is constant throughout the LV domain, was estimated by minimising a cavity pressure error functional: where P(t) is the cavity pressure measurement at time t. The model-derived cavity pressure (t) coincides with the Lagrange multiplier used to enforce the volume constraint in Eq. 9, as described in our earlier works Asner et al. 2017).

Numerical solution
Biomechanical models were solved using a standard Galerkin FE scheme. Displacement and hydrostatic pressure variables were interpolated using quadratic and linear tetrahedral elements, respectively, to ensure the mixed formulation's stability (Hadjicharalambous et al. 2014b), while quadratic triangular elements were used for epicardial and basal multipliers. All numerical problems were solved using C Heart, a multiphysics FE solver (Lee et al. 2016).

In-situ model ( ˝0 =˝I N−SITU )
In the in-situ modelling scenario, the model reference domain was provided by the in-situ MR image of the unloaded myocardial configuration (Sect. 2.1). To simulate the passive mechanics, the reference domain was inflated to the ED cavity volume, through Eq. 9. In terms of basal and epicardial boundary conditions, numerous perturbations of Eq. 12 were considered, to identify a combination able to sufficiently match the ED geometry. Passive and active parameters were evaluated as described in Sect. 2.3.4 allowing for full-cycle simulations.

End-systolic model ( ˝0 =˝E S )
For this model scenario, the segmentation of the ES image provided the unloaded domain. Base displacements were prescribed through Eq. 10, while deformation in the region adjacent to the ribcage was constrained through Eq. 11.

Inverse model ( ˝0 =˝I NV )
The final model scenario employed an inverse methodology to estimate the unloaded domain. The aim was to evaluate the influence of inverse approaches on model outcomes; therefore, any information derived from the in-situ reference image was ignored. Model personalisation was thus performed based only on cavity pressures P(t), ED geometry and V Klotz . The basal surface was kept constrained, as is common in most inverse cardiac models (Eq. 10, u d = 0 ). Two constraints were considered for the epicardial surface: a zerotraction condition, used almost exclusively in inverse cardiac mechanics applications, and a ribcage condition restricting the deformation of the rib-adjacent epicardium (Eq. 11).
INV was estimated using the approach of Rajagopal et al., whereby the Jacobian matrix is computed with respect to reference state parameters (Rajagopal 2007). Inverse simulations were run with pressure inflation of the LV, summed up as estimating the reference domain INV which would provide the ED domain when inflated to P ED . Joint estimation of INV and passive parameters was performed as described in Sect. 2.3.4, followed by active tension estimation.

Boundary conditions for the in-situ model
Several combinations of epicardial and basal boundary conditions (Sect. 2.3.3), along with perturbations of relaxation parameters ( b , e and a ), were considered for the in-situ model. The suitability of each combination was assessed using the distance error metrics in Sect. 2.3.4. A representative illustration of perturbations of relaxation parameters is  Fig. 3, when the BAR boundary condition was employed. Good identifiability characteristics were observed for each error metric, with a unique minimum in each case. Notably, the combinations of ( b , e ) providing minimum errors were similar in all three metrics, suggesting that ( b = 5 × 10 −6 , e = 5 × 10 −7 ) would provide a reasonable choice. Figure 4 illustrates the overall effect of the employed combination of boundary conditions on model accuracy. For all error metrics, the maximum value was observed when only the basal boundary condition was employed, while the introduction of the apical constraint caused a substantial reduction. The decrease in error was even more striking when the ribcage condition was employed, leading to a more than 50% reduction in the maximum error.

Boundary conditions for the end-systolic model
The choice of boundary conditions was also examined for the case of a data-derived geometry employed as the reference domain. Several combinations of relaxation parameters b and e were run, with the combination ( b = 1 × 10 −6 , e = 5 × 10 −6 ) found to produce consistently good results across three error metrics (Fig. 3).

Boundary conditions for the inverse model
Numerous perturbations were also considered for the inverse model to investigate the suitability of boundary conditions. However, simulations employing the relaxed basal boundary condition and zero epicardial energy exhibited considerable numerical issues, with the vast majority failing to converge to a solution. In particular, a wide range of passive parameters were employed ( a sim ∈ [50, 10000] Pa and a sim f ∈ [50, 10000] Pa), with all cases with a, a f < 7 kPa presenting simulation failure. Certain combinations of larger parameters would converge to a solution; however, such stiff parameters would result in estimated reference domains with very high cavity volumes (e.g. V ref = 57 mL for a = a f = 10 kPa) which differ substantially from V Klotz = 44.27mL. Additionally, several values were considered for the penalty parameter of the basal constraint ( b ∈ [10 −9 , 10 −2 ] ). Although a weaker enforcement of the constraint would allow the simulations to run to a larger loadstep, simulation failure would still occur. Finally, to ensure that these numerical difficulties were not simply an issue of the solution being difficult for the solver to find, different loading steps were considered (in parallel with perturbations of passive parameters and b ). Accordingly, the loading steps for the cavity pressure would increase from 10 to 100 and 1000. Although this helped the simulations run to a lower cavity volume, simulations would still fail to run through all loading steps. Visualisations of model deformation prior to simulation failure revealed nonphysiological deformation including surface bending and pinching (Fig. 5), with the irregular indentation observed mid-length of the heart's long axis, being present in all failing simulations.
On the contrary, similar computational issues were not present when inverse simulations accounted for the ribcage influence on LV motion. The majority of these simulations ran without any issues (see "Appendix Investigation of passive parameters estimation for inverse model", Fig. 11), for most combinations of passive parameters (with the exception of very compliant material, e.g. a, a f < 100 Pa). When the constraining influence of the ribcage was accounted for, inverse estimation provided reasonable geometries as reference domain estimates (Fig. 5), avoiding any spurious deformation modes.

Comparison of passive mechanics
For the in-situ and end-systolic models, parameter estimation was performed through parameter sweeps (Sect. 2.3.4). Simulations with varying values of the isotropic parameter were run ( a sim ∈ [50, 2000] Pa), while the fibre parameter was set ( a sim f = 1 kPa). The range of parameter ratios a∕a f ∈ [0.05, 2] considered was based on in-vivo estimates of human heart models, employing the reduced Holzapfel-Ogden law . Absolute values of a and a f are provided in Table 2, while the errors achieved for the estimated parameters for each model are presented in Table 3. In addition to the nodal distances, model errors were also assessed with respect to the Dice similarity coefficient D (Zou et al. 2004). Specifically, D was used to compare images generated from the simulated ED geometry of each model versus the segmented ED image. D values close to one denote a good agreement between model and data, whereas values closer to zero denote a poor match.
For the inverse model, joint estimation of passive parameters and reference domain was performed. The fibre parameter was not estimated but instead set to an initial value ( a f = 1326 Pa), and only the isotropic parameter was allowed to vary ( a sim ∈ [50, 5000] Pa). The range of passive parameters considered was based on in-vivo estimates for healthy and diseased human hearts . A unique minimum was observed for a in this Fig. 5 Views of simulated inverse geometries are provided and are overlaid onto the insitu SAX and 4CH views for comparison. Top row: inverse simulation employing a zero epicardial boundary condition, with a = 0.2 kPa, a f = 2 kPa and e b = 10 −3 , which failed at loadstep 52 out of 100, at a cavity volume of 59.6 mL. Middle row: inverse simulation employing a zero epicardial boundary condition, with a = 5 kPa, a f = 10 kPa and e b = 10 −6 , which failed at loadstep 68 out of 100, at a cavity volume of 52.4 mL. Bottom row: inverse simulation employing the ribcage epicardial boundary condition, with a = 0.8 kPa, a f = 1.326 kPa and e b = 5 × 10 −6 , which ran without issues and provided a reference domain with cavity volume 44.29 mL  Table 2 and to the estimated reference domain depicted in Fig. 5. It should be, however, noted that this estimate of a is strongly coupled to the value selected for a f , as discussed in "Appendix Investigation of passive parameters estimation for inverse model".
The three models were also compared with respect to their ability to model the nonlinear stress-strain properties of the passive myocardium. Figure 6 presents the stress-strain response in the fibre direction over the passive inflation of the LV, as it is derived by the three models. The three curves present good agreement, suggesting that the three models produce consistent stress-strain properties and are thus able to consistently capture the nonlinear tissue response, despite differences in the parameters employed.

Comparison of estimated active tension
Active mechanics were characterised following the procedure in Sect. 2.3.4. Parameter sweeps were performed for each timestep, with (t) varying between [0, 150] kPa, following in-vivo estimates of active tension (Asner et al. 2015a) . For all modelling scenarios, a unique minimum of the pressure error was achieved for each timestep. Estimated active tension curves are depicted in Fig. 7, maximum values are presented in Table 2, while the agreement with the ES geometry is assessed in Table 3.

Comparison of stresses and strains
For each modelling scenario, full-cycle simulations employing the estimated passive and active parameters were run. with overlaid the end-systolic mesh (in white) and the simulated endsystolic geometries when employing the in-situ (in green), end-systolic (in red) and inversely estimated (in blue) unloaded domains The Cauchy stress and Green-Lagrange strain along the fibre direction were computed, both at ED and ES. For each model, strains were computed with respect to the reference domain employed. Average values are presented in Table 4, while their spatial distributions are illustrated in Figs. 8 and 9.

Assessment of boundary conditions
Exploiting the invaluable data on the myocardial reference domain and its relation to in-vivo geometries, we were able to develop a computational framework for reliable forward and inverse simulations, transitioning between the unloaded and ED domains. Error quantification revealed the importance of accounting for the surrounding tissues, with a remarkable improvement in model errors reported when the ribcage boundary condition was employed for all models. Importantly, inverse simulations employing only a basal constraint exhibited numerical issues and predicted geometrical abnormalities, suggesting that the ribcage constraint was not only important for model accuracy, but also necessary for simulation convergence. The ribcage boundary condition led to more than 50% reduction in errors, while errors were lowest when all constraints were enforced. While this is not surprising as larger parts of the geometry were penalised, the relative impact of constraints is interesting. Inclusion of the apical constraint resulted in improved agreement with the data-derived geometries, yet this effect was not as critical as neglecting the ribcage influence which led to nonphysiological cardiac deformations, penetrating the ribcage surface.

Influence of reference domain on model outcomes
The employed reference domain was shown to affect passive parameter quantification. Estimates obtained with the endsystolic model suggested more anisotropic material properties compared to in-situ and inverse models, in agreement with earlier studies (Nikou et al. 2016). The parameters estimated within the inverse framework were in better agreement with the parameters of the in-situ model; nevertheless, their values are not unique and depend on the joint estimation scheme employed ("Appendix Investigation of passive parameters estimation for inverse model"). Despite variations between models, passive parameters were of the same order of magnitude. For comparison, a and a f exhibited fivefold increase in patients with dilated cardiomyopathy compared to healthy volunteers . Therefore, while the reference domain influences the absolute parameter values, the correct order of magnitude was captured by all three models. Notably, all models were able to capture the data-derived diastolic part of the cycle, with satisfactory agreement to the ED geometry and pressure-volume loop, while producing consistent stress-strain response.
The unloaded domain had little qualitative effect on active characteristics, with active tension curves following the same temporal evolution. However, quantitative differences of 15.2% and 8.8% were observed in max for the end-systolic and inverse models, respectively, compared to in-situ values. Although active tension estimation neglected any displacement observations, there was a fair agreement between the simulated and data-derived ES geometries, for all cases. Notably, despite differences in the estimated active tension, all modelling scenarios were able to reproduce the pressure-volume loop, suggesting that such gross metrics are insensitive to the selection of reference domain.
Substantial variations were observed in ED stresses. The end-systolic model gave higher ED stresses, in agreement with previous studies reporting lower ED stresses when an unloading algorithm was employed (Nikou et al. 2016). Apart from geometrical discrepancies in the employed reference configurations, the different passive parameters of the three models might contribute to ED stresses' variations, as passive parameters have a significant correlation with ED stresses (Behdadfar et al. 2017). On the contrary, ES stresses did not vary markedly in the three modelling scenarios, with both the global average values and their spatial distributions showing considerable similarities.
ED strains presented similar spatial patterns, particularly between the in-situ and inverse models. The end-systolic model exhibited larger average and maximum values, potentially due to its larger cavity volume difference to the ES mesh. End-systolic strains presented very similar spatial distributions, suggesting their insensitivity to the employed reference domain.
Notably, the variations in model outcomes resulted from reference domains with very similar cavity and wall volumes (Table 1), highlighting the importance of incorporating patient-specific geometries and boundary conditions.

Evaluation of common choices for the unknown reference domain
Usage of in-vivo frames as the unloaded configuration has been very popular in cardiac mechanics, due to its ease in implementation using routinely acquired data and low computational cost compared to inverse schemes. This approach is particularly attractive for clinical applications, as it does not require any invasive clinical data for modelling the unloaded domain. However, based on the findings of this study, certain biomechanical properties derived from the ES model differ from in-situ estimates in terms of absolute values (e.g. material parameters). Models employing image frames as the reference domain can nevertheless provide reliable clinical tools for investigating the comparative effect of a treatment or the impact of a pathology, provided the same data frame is employed for all cases (i.e. reliable intracomparative analysis). Such models can also be used reliably for inter-comparisons, for quantities which are insensitive to the choice of reference domain such as ES stresses and strains or cavity pressures. For instance, they could be valuable for personalised therapy planning in the case of Cardiac Resynchronization Therapy (CRT), whereby the peak LV pressure time derivative is used as a clinical index (Sermesant et al. 2012). Furthermore, models employing an imagederived reference domain could be useful for disease stratification, focusing on quantities whereby any deviations from in-situ values are much smaller than the difference observed between health and disease. In the case of dilated cardiomyopathy for example, the stiffness difference between patients and normals (fivefold increase)  is much larger than the variation between the ES and in-situ parameters observed in this work, suggesting that stiffness estimates obtained with the ES model can be reliably used to distinguish between normals and heart failure patients. However, this is not the case for the maximum active tension max , whereby the difference observed between patients and healthy volunteers (~ 5 to 30%) ) is similar to the error in max estimates between ES and in-situ models (15.2%). Accordingly, models employing image frames as the unloaded domain have numerous applications as reliable clinical tools, yet, care should be taken to prevent biased inter-comparative analyses.
Alternatively, inverse methodologies have been employed in many cardiac mechanics models to provide a geometry which-unlike in-vivo frames-is under no load. The approach followed within led to a reference estimate with a cavity volume very close to the ES frame (Table 1), yet very different geometries (see Figs. 1 and 5), which in turn resulted in substantial variations in results. The passive and active parameters along with ED stresses obtained with the inverse model were in closer agreement with in-situ values, compared to the end-systolic model. Studies interested in such quantities-which appear to be sensitive to the choice of reference domain-might therefore consider an inverse approach. It should be noted, however, that inverse methodologies incur an additional computational cost and are less straightforward to implement in existing FE solvers compared to image-driven forward FE methodologies, factors which might limit their clinical applicability. In that case, researchers should weigh the benefits of enhanced accuracy over the additional computational cost.
Finally, both approaches could be clinically useful in cases whereby the comparative effect of a treatment or intervention is under consideration. Additionally, both cardiac modelling techniques could be reliably used to provide a mechanistic understanding of heart function and pathology, by linking mechanisms of action to disease progression or by examining the interactions of a medical device with the mechanics of the heart.

Study limitations
The experimental procedure followed enabled us to reliably capture the heart's unloaded domain, creating a unique opportunity for quantitatively assessing its influence on model outcomes. However, the pressure in the postmortem setting was not likely to be exactly zero, which would exert some load in the myocardium-albeit much lower than invivo. Due to the animal's positioning, additional loading stems from gravity which is also likely to induce deformation in our imaged reference domain. However, this effect should be similar in the in-vivo images and is thus unlikely to affect the comparative analysis.
Furthermore, although studies have employed different image frames as the reference domain (e.g. end systole, early diastole or diastasis), only one data frame (the end systolic) was examined. However, considering different data frames as the unloaded domain might lead to differences in the absolute values of the data-derived outcomes. In fact, the use of a different diastolic frame as reference was found to influence the estimate of the parameter ratio a∕a f ) as well as affect the active tension estimation (Xi et al. 2013). Accordingly, examining different image frames would provide an invaluable future contribution, as it might reveal a specific image frame providing better agreement to in-situ values and thus more suitable to be used as the unloaded domain.
Similarly, despite many iterative schemes for the joint estimation of the reference domain and material parameters, only one approach was considered herein. The evaluation of the current approach was valuable in highlighting important aspects to consider when running inverse simulations, the associated numerical difficulties, and more importantly the necessity for choosing suitable boundary conditions. However, the specific joint estimation scheme employed might had introduced additional uncertainty. While the in-situ and end-systolic models should provide unique passive parameters (Hadjicharalambous et al. 2014a), parameter uniqueness cannot be guaranteed for the inverse model, whereby an almost arbitrary value was selected for a f . Nevertheless, different passive parameters ("Appendix Investigation of passive parameters estimation for inverse model") led to similar active tension, stresses and strains, suggesting that our conclusions would hold for similar schemes.
Parameter estimates are also affected by the boundary conditions employed ). This dependence was observed in the passive estimates of the in-situ model with different basal constraints ("Appendix Investigation of passive parameters estimation for inverse model"); nonetheless, the effect on active tension was less pronounced. Boundary conditions are also likely to affect stress and strain fields, both directly and indirectly through their influence on passive parameters. Nevertheless, similar boundary conditions were employed in all models, enabling reliable comparisons. Additional boundary conditions could also be considered in the future, to account, for instance, for the effect of the right ventricle on the LV, which is expected to be substantial Kolipaka et al. 2010).
In the present study, all models assumed negligible residual stresses (owing to growth and remodelling) within the heart tissue. Nevertheless, residual stresses are likely present in the myocardium throughout the cardiac cycle (Genet et al. 2015), as evidenced by opening angle experiments on excised hearts which revealed compressive circumferential endocardial stresses and tensile epicardial stresses (Omens and Fung 1990). However, in order to properly quantify residual stresses in a subject-specific manner, opening angle experiments would have to be conducted for the heart under consideration, following the in-situ image acquisition. Such experiments are very lengthy and particularly challenging, as they require the excised muscle to remain viable for the duration of the experimental procedure. An additional challenge is posed by the need to model residual stresses in a personalised manner. Modelling studies have sought to characterise residual stresses through appropriate constitutive laws (Wang et al. 2014) or as a result of heterogeneous tissue growth (Genet et al. 2015), yet they often involve additional parameters that would need to be estimated for the individual heart. However, such modelling studies have shown that accounting for residual stresses has a smallyet non-negligible-effect on diastolic (Wang et al. 2014) and systolic characteristics (e.g. decrease in the peak fibre systolic stress (Taber 1991)). Particularly, residual stresses were seen to lead to higher epicardial ED stresses (Wang et al. 2014), while the compliance of the ventricle was found to increase when a higher level of pre-stress was assumed (Genet et al. 2015). While it is likely that such biases were present in our analysis, understanding the impact residual stresses would incur on our findings warrants further investigation and would provide a valuable future development.
Finally, the conclusions of this study stem from a porcine subject. While direct extrapolation to human hearts cannot be guaranteed, important similarities in cardiac function between the two species suggest that similar conclusions might hold for humans. Importantly, these experimentswhich would not had been possible on humans-enabled a thorough analysis of common choices for unloaded domain, and a systematic assessment of their influence on model biomechanics.

Conclusions
This work deals with a major limitation in models of cardiac mechanics, namely the unknown reference domain. Exploiting validation data from in-vivo and in-situ experiments on a porcine heart, this study has employed an in-silico modelling framework to quantitatively assess the effect of the unloaded domain on model-derived biomechanics. Personalised whole-cycle cardiac models were developed, employing common approaches for selecting the reference domain (image-derived, inversely estimated), and were systematically compared against in-situ data. Our results give credit to both approaches, but also elucidate important discrepancies in their estimates of material parameters, stresses and strains, during both the passive and active phase of the cardiac cycle. Notably, they draw attention to the importance of patient-specific image-driven boundary conditions, in both forward and inverse biomechanical models.
However, an important advantage of this study is that the reference domain is actually available. Accordingly, estimation does not need to rely only on the gross estimate V Klotz but could also take into account the actual geometry of the unloaded state. In this case, parameter estimates of the inverse model were acquired by matching the simulated reference domain to the in-situ unloaded domain. Estimation was performed by comparing the reference geometries resulting from parameter sweeps with the data-derived reference domain. However, these simulations were run with zero motion applied on the basal boundary, as it was assumed that the relative position between the reference and in-vivo images was not known. Accordingly, the distance errors were calculated after first registering each estimated unloaded geometry to the reference mesh, through point set registration with ICP. Figure 11 demonstrates the behaviour of the RMSE over the parameter space, which presented a unique minimum. Interestingly, the parameters estimated through this inverse process ( a = 500 Pa, a f = 1700 Pa) were not very close to the parameters estimated with the forward process ( a = 792 Pa, a f = 1200 Pa) in the in-situ model. This discrepancy in parameters could be attributed to the difference in the boundary conditions employed. Specifically, the insitu parameters were obtained with a prescribed motion on the basal plane, while inverse parameters were obtained with a fixed base plane.
To further investigate this issue, a separate set of forward simulations employing IN−SITU as the unloaded domain was run. Specifically, following the boundary conditions of the inverse process, LV inflation simulations of the in-situ reference domain were run with zero displacement enforced weakly on the basal plane. The passive parameter estimates ( a = 187 Pa, a f = 1876 Pa) were then obtained using parameter sweeps and ICP registration to the end-diastolic domain, comparing more favourably to the parameters estimated with the inverse process ( a = 500 Pa, a f = 1700 Pa), as illustrated with a yellow circle in Fig. 11. Interestingly, when employing similar boundary conditions, the forward and inverse simulations led to similar model outcomes, strengthening the credibility of the modelling framework employed. It is also worth noting that the active tension was not substantially affected by the basal boundary condition in the in-situ model, giving a maximum value of max = 133 kPa for zero basal motion versus max = 125 kPa when displacement was prescribed.
Funding Financial support was received from the University of Cyprus, the Engineering and Physical Sciences Research Council (EP/ N011554/1 and EP/R003866/1) and from Swiss national science foundation (PZ00P2_174144).
Availability of data and material Imaging data are available upon request.
Code availability All numerical problems were solved using C Heart (Lee et al. 2016).

Conflict of interest
The authors declare no conflict of interest.

Fig. 11
Behaviour of objective functions relying on cavity volumes (left) and distance (right) metrics, over the parameter space, in inverse simulations. The absolute value of the difference between the simulated cavity volume and V Klotz is presented in mL, while the RMSE between the estimated unloaded domain and the in-situ refer-ence mesh is given in mm. White circles depict the combination of passive parameters providing the lowest error values, while the yellow circle corresponds to the parameter estimates obtained in the forward simulations' framework