Improving cardio-mechanic inference by combining in vivo strain data with ex vivo volume–pressure data

Cardio-mechanic models show substantial promise for improving personalised diagnosis and disease risk prediction. However, estimating the constitutive parameters from strains extracted from in vivo cardiac magnetic resonance scans can be challenging. The reason is that circumferential strains, which are comparatively easy to extract, are not sufficiently informative to uniquely estimate all parameters, while longitudinal and radial strains are difficult to extract at high precision. In the present study, we show how cardio-mechanic parameter inference can be improved by incorporating prior knowledge from population-wide ex vivo volume–pressure data. Our work is based on an empirical law known as the Klotz curve. We propose and assess two alternative methodological frameworks for integrating ex vivo data via the Klotz curve into the inference framework, using both a non-empirical and empirical prior distribution.


INTRODUCTION
The last decade has seen impressive progress in the mathematical modelling of soft-tissue mechanics (Chabiniok et al., 2016;Mangion et al., 2018), and the specific application to cardio-mechanics promises to improve our understanding of cardiac physiology and pathophysiology (Gao et al., 2017b;Hadjicharalambous et al., 2015;Krishnamurthy et al., 2013;Maso Talou et al., 2020;Peirlinck et al., 2019;Wang et al., 2013). In a recent case-control study involving patients diagnosed with myocardial infarction (MI) and healthy volunteers, we have demonstrated that the biophysical parameters of a cardio-mechanic model of the left ventricle (LV) of the heart can be non-invasively estimated from cardiac magnetic resonance (CMR) images, and that the parameters thus inferred allow an accurate classification of the patients' disease status (Gao et al., 2017a). However, the cardio-mechanic equations have no-closed-form solution and require computationally expensive numerical procedures based on finite element discretization (Bathe, 2006). During optimization or sampling, these procedures have to be repeated iteratively, leading to excessive computational run times that prevent uptake of the methodology in the clinical practice. To address this challenge, we have demonstrated in two recent proof-of-concept studies that the original cardio-mechanic model can be emulated with a statistical surrogate model (Conti et al., 2009;Sacks et al., 1989), allowing a reduction of the computational complexity by three orders of magnitude at negligible loss in prediction and estimation accuracy Noè et al., 2019). While these studies are a promising first step on the pathway to clinical translation, they are limited by the fact that the LV geometry was assumed to be known. To apply the method to new patients, variations in the LV geometry need to be taken into consideration. Since there is no time to train a new emulator with a different LV geometry once a patient has been admitted to the clinic, the LV geometry needs to be included as an additional input to the emulator, along with the cardio-mechanic parameters. The LV geometry is typically represented by a mesh of in the order of 10 4 finite element nodes, so a substantial dimension reduction is needed. In previous work (Romaszko et al., 2019) we have shown that a projection of this high-dimensional ventricular mesh space onto a space spanned by the first five to ten leading principal components gives an accurate reconstruction of the LV geometry, and that dimension reduction with more advanced methods based on convolutional neural networks does not give any significant improvement over dimension reduction with principal component analysis (PCA). However, the prediction accuracy of the cardiac emulator is adversely affected by the increase of the dimension that the emulator now has to cover. In addition, we have the problem of intrinsic poor identifiability of certain biomechanical parameters when trying to infer them from the circumferential strains typically extracted from CMR scans (Gao et al., 2015), while longitudinal and radial strains are difficult to extract at high precision. The purpose of the present article is to investigate to what extent the loss in parameter estimation accuracy can be compensated by the inclusion of prior knowledge derived from complementary ex vivo data.

A mathematical description of the diastolic filling process
The heart consists of four chambers that work together with the circulatory system during the cardiac cycle to provide blood rich in oxygen and nutrients to the rest of the body. Each chamber operates in unison, with the left ventricle responsible for pumping the blood out of the heart on its way to the body's organs. This process is permitted by the pump function of the LV, which occurs in two stages. The first is called diastole, where the chamber fills with blood and the second is systole, where the blood is ejected out of the cavity. The effectiveness of these processes is largely determined by the mechanical properties of the myocardium, which is a muscular tissue surrounding the cavity. During diastole, this muscle is passive and acted upon by the blood within the LV. The adherence of the muscle to this diastolic filling process contributes to the ability of the LV to fill with blood. Corresponding to the stiffness of the material, this adherence can be described by the mechanical properties of the tissue and provides knowledge of the level of LV function or dysfunction. Traditionally, these properties were learnt invasively through ex vivo testing by harvesting myocardial tissue from the body and subjecting it to a series of stretch tests. Of interest in this paper is their estimation from a set of in vivo measurements, using a non-invasive methodology that would be relevant for clinical decision making. Non-invasive estimation of the mechanical properties involves estimating the parameters of a mathematical expression-called a constitutive law-that describes the build-up of energy in the myocardium under deformation. In our mathematical model, this is provided by the Holzapfel-Ogden (HO) law: where a, b, a i , b i , (i = f, s, fs) are the material parameters. When properly estimated, these provide a summary of the mechanical properties of the tissue. A full description of the HO law is outwith the scope of this work, and interested readers should refer to the literature (Holzapfel & Ogden, 2009). Motivated by a local sensitivity study carried out in the literature (Gao et al., 2015), we reduce the dimension of the parameter space from 8 to 4 by adopting the following transformation: This reparameterization was implemented in previous studies Noè et al., 2019). We will focus on determining the parameters of this low dimensional mapping, 1 , 2 , 3 and 4 with values of a 0 , b 0 , a f0 , b f0 , a s0 , b s0 , a fs0 and b fs0 taken from the literature (Gao et al., 2017a). Given these constants, we assume each i is bounded by [0.1, 5] as used by Davies et al. (2019) and Noè et al. (2019).
The Holzapfel-Ogden law can be used in combination with a set of momentum balance equations to simulate the dynamics of the LV in diastole. This simulator is treated as a black box for the work in this study: a computer code that takes a set of inputs and produces some output. In our context, the inputs are the material parameters, an LV geometry and the end diastolic pressure, as visualized in Figure 1. After a simulation time of around 8 min (this varies depending on the material parameters and the hardware used), we obtain a predicted left ventricle at end of diastole, from which we can calculate the end diastolic volume of the LV cavity and a set of circumferential strains at different regions on the LV wall. The values that are predicted by this simulator can also be measured from CMR images in the clinic, allowing us to optimize the material parameters in order to match the measured data. To do so, we require a fixed patient-specific LV geometry and end diastolic pressure. For the in vivo behaviour, we assume this end diastolic pressure is 8 mmHg, in line with other studies . The patient-specific geometry, an example of which is provided on the right of Figure 1, is obtained by segmenting the left ventricle in the CMR scans. The final LV geometry is in the form of a series of nodes, of which there are 2896 on both the inner and outer surfaces (endocardium and epicardium respectively) of the left ventricle. Note that this constitutes a 2 × 3 × 2896 = 17,376 dimensional representation in three dimensional Cartesian coordinates. The high computational costs for solving the system of equations using the finite element method conflicts with the objective of in clinic implementation of cardiac models, leading to the recent introduction of machine learning and computational statistics methods to the field of cardiac modelling (Dabiri et al., 2019;Maso Talou et al., 2020;Peirlinck et al., 2019). These approaches seek to replicate the behaviour of the expensive mathematical models by replacing the expensive simulations with cheap statistical approximations. Methods considered vary from boosted tree based approaches (Dabiri et al., 2019) to Gaussian processes (Peirlinck et al., 2019) and deep neural networks (Maso Talou et al., 2020). The reduced computational costs are particularly relevant in the context of parameter estimation where optimization or statistical sampling schemes require repeated evaluation of the cardiac models, and approximate techniques have already been shown to achieve results similar to brute force optimization of an expensive simulator loss function Noè et al., 2019). The importance of successful parameter inference is emphasized in the work done by Gao et al. (2017a) to successfully classify MI based on mechanical properties of the LV.
Recently, the subject of LV geometry representation in cardio mechanical modelling has seen increased relevance. Specifically, PCA (Shlens, 2014) has drawn particular interest. Wang et al. (2020) implemented this dimension reduction technique for the purpose of load-free LV geometry estimation in the context of inferring the passive parameters. Perhaps most similar to the method of dimension reduction considered in the current paper, Maso Talou et al. (2020) considered a reduced geometry representation to represent the LV in a statistical model approximating the diastolic filling process, which is replicated by deep neural networks. We build on this latter piece of work by quantifying the performance of the reduced LV geometry cardiac model when attempting to learn the material properties while also properly quantifying any uncertainty in our parameter estimation through the use of computational Bayesian statistics based on Markov Chain Monte Carlo (MCMC).

Emulation of the in vivo likelihood function
Of interest is the efficient estimation of the parameters = ( 1 , 2 , 3 , 4 ) from Equation (2). The data available from in vivo CMR images are K = 24 circumferential strains at well-defined positions on the endocardial wall and the LV cavity volume, all measured at end diastole. In what follows, we denote the measured quantities by where y * 0 is the measured LV cavity volume at end of diastole (LVV), and y * 1∶ = (y * 1 , … , y * 24 ) are the measured circumferential strains. For the inverse problem, we require the end diastolic pressure of the blood within the LV cavity. In principle this quantity can be measured, but the procedures to do so are invasive, meaning this value might only be available in the case of diseased patients. Instead, we take the approach of Davies et al. (2019) and assume a fixed end diastolic pressure (EDP) of 8 mmHg. In particular, this means that y * 0 is the measured LVV at an EDP of 8 mmHg. The true quantities, representing the underlying values of the data from Equation (3), are modelled by the emulator: where the indices have the same meaning as above. Note that these outputs depend on the cardio-mechanic parameters and the LV geometry representation 1 , , which when made explicit in the notation gives us: In Gao et al. (2015) the cardio-mechanic parameters for a given LV geometry  were estimated by minimizing the L2 norm of the difference between y and z: Under the assumption that the measurement noise is iid additive Gaussian with variance 2 m : we obtain the log-likelihood z( , ): Maximizing the likelihood with respect to , for a given LV geometry , is equivalent to minimizing the original loss function (6), where the model outputs z( , ) have been set to the outputs of the emulator. This approach was adopted by Noè et al. (2019) and shown to outperform emulation of the loss function by Davies et al. (2019). Unlike the current work, they considered statistical emulation and parameter estimation for a fixed LV geometry.

Introduction of ex vivo information and the aims of this study
The behaviour of the myocardium is governed by components operating at a microscopic scale within the tissue whose influence changes with the level of stretch in the muscle (Voorhees & Han, 2015). This behaviour is replicated by the HO law from Equation (1), where different mathematical terms take effect at different stretches. In particular, some of the parameters in this expression relate to the behaviour at small stretch and some to the behaviour at larger stretch. This property is useful in the forward problem, allowing us to simulate the behaviour of the tissue at stretches beyond those typically observed in vivo. Unfortunately, this introduces identifiability problems to the inverse problem since the in vivo data do not inform us about the behaviour of the material at higher stretches. To overcome this problem, we aim to introduce ex vivo information to the inference routine through a relationship known as the Klotz curve (Klotz et al., 2006). Klotz et al. (2006) investigated the dependence between pressure and end-diastolic cavity volume (LVV) over a large population of various diseased patients and healthy volunteers. In their study, the end-diastolic volume-pressure curves varied substantially and spanned a wide range of volumes. However, they found that on standardizing the LVV, the curves become nearly identical, and that despite different disease status, the shape of the curve is nearly the same (figure 2 in Klotz et al. (2006)), with a functional form provided by: V P is the volume at pressure P, V 30 is the volume at pressure 30 mmHg and V 0 is the volume of the unloaded left ventricle, at pressure 0 mmHg. We will seek to introduce this information to the inference problem by penalizing deviations from this mathematical relationship. The Klotz curve has already been explored in the context of cardio mechanical modelling (Hadjicharalambous et al., 2015;Krishnamurthy et al., 2013;Palit et al., 2018;Peirlinck et al., 2019). We will expand on the current literature in the following ways: firstly, we will introduce a prior distribution based on the Klotz curve, allowing us to correctly introduce this information in a Bayesian inference scheme for the material parameters. These priors will be implemented using statistical emulators, enabling the introduction of the Klotz information to the parameter inference procedure while avoiding any large computational costs. Secondly, we will seek to use these priors to improve parameter inference in the presence of LV geometry errors and information loss caused by the dimension reduction, reducing the uncertainty and improving the accuracy of the parameter estimation. Thirdly, we carry out a comparative performance assessment of two alternative prior representations, using both synthetic data and real CMR images, and we quantify the effect that the inclusion of the empirical Klotz law has on the inference. This effect will be assessed using the marginal posterior distribution, which we take as representing the practical identifiability of the parameters.

METHODOLOGY
Two alternative prior distributions will be formulated for the introduction of information from the Klotz curve, each derived from the parametric form suggested by Klotz et al. (2006) and provided in Equation (9). The first distribution is called the nonempirical Klotz prior, and its evaluation depends on the data obtained from the study in Klotz et al. (2006) y ex . The second prior depends on y ex and the measured LV cavity volume y * 0 . Due to its dependence on the measured in vivo data, this prior is called the empirical Klotz prior. To summarize: the important difference between the empirical Klotz curve prior and the nonempirical Klotz prior is the dependence of the former on the in vivo measured LVV. The two priors are described in Table 1.
Both of the Klotz priors will require samples of the parameters from Equation (9), for which a rearrangement of the original function is considered, expressing the normalized volume as a function of pressure. For now, consider this relationship to be of the formṼ = f (p, ) where f depends on some parameters . The form of this function is deliberately not made explicit at the moment so that the discussion generalizes to both of the Klotz prior distributions. Assuming a heteroscedastic noise model:Ṽ where log ( the parameters , c 1 , c 2 and c 3 are sampled from the posterior distribution conditional on the data from Figure 2. The posterior mean values of c 1 , c 2 and c 3 are used to provide a heteroscedastic noise model estimate. Prediction of the Klotz curve is shown in Figure 2, along with uncertainty bands obtained as ±2 (P).

Method A: The nonempirical Klotz prior
The emulated form of the cardio-mechanic model from Section 2.2 is used to model the in vivo data, represented by the right branch in Figure 3 with an arrow from the cardio-mechanic model to y * . To include extra information, obtained from the ex vivo study performed by Klotz et al. (2006), information from the relationship in Equation (9) will be introduced to the inference problem. To do so, we require a volume prediction at an EDP other than 30 mmHg. Using a low EDP value makes little sense because this information already comes into the model through the

Prior Description
Nonempirical Klotz prior This uses the normalized form of the Klotz function from Equation (9), not depending on any in vivo end of diastole measurements Empirical Klotz prior This uses a rearrangement of (9) where we evaluate the Klotz predicted V 30 using the measured LVV y * Observed quantities are in grey circles, functions are in boxes and variables are in white circles. In addition to using the cardio-mechanic model for modelling in vivo data (y * ) at the in vivo reference pressure (right branch), it is also used to model the LVV at higher ex vivo reference pressures, providing a distribution on the normalized volume, which is labelled asṼ e 20 . In addition, an empirical law based on the function from Equation (9) fitted to ex vivo data (y ex ) provides a distribution on the normalized volume at the same reference pressure (left branch), with noise model assumed to be heteroscedastic (see details below (10)). Despite the fact that the normalized volume resulting from the emulator and the Klotz relationship are actually the same variable, we choose to distinguish between these based on their distributions (Ṽ e 20 for emulator induced distribution andṼ Kl 20 for Klotz induced.) The two distributions are combined via a product-of-experts approach (dashed line) likelihood function. Since the use of the Klotz prior was intended to introduce information about the high stretch behaviour of the material, a pressure of 20 mmHg seems appropriate.
Eliminating from the Klotz curve function in Equation (9) (see section 5.1 of the SM), we can evaluate the following function for some parameter value , providing a prediction of the normalized volume at pressure 20 mmHg based on the Klotz relationship: Recalling the general form of the Klotz noise model in Equation (10), we have the following function: so that = . Information from the ex vivo data is introduced by inferring the parameters of this model , c 1 , c 2 and c 3 (the latter three parameters provide the heteroscedastic noise), from the Klotz data plotted in Figure 2. Assuming priors of the form ∼ u(0, 10), 2 c 1 , c 2 , c 3 ∼ u(−100,100), where the bounds are deliberately overly conservative, these parameters were sampled using MCMC. The posterior distribution on is then approximated with a Normal distribution: the validity of which was confirmed using the Shapiro-Wilk test.
Estimating the mean and standard deviation using the samples from the posterior, we can use the delta method (see section 5.1 of the SM) to propagate the uncertainty in through to the normalized volume via (12): giving a distribution onṼ 20 induced by the Klotz data and function. This part of the model is represented by the left branch of Figure 3, where we relate the normalized volume with the ex vivo data, y ex , through the Klotz curve model. Instead of predicting the normalized volume using the Klotz curve, we can predict this quantity for a given set of material parameters using the simulator from Figure 1, linking the normalized volume from Equation (9) and the material parameters. Training Gaussian process emulators, we learn a probabilistic representation of the LVV at pressures 20 mmHg and 30 mmHg. Using the delta method (see section 10 of the SM), we can propagate the uncertainty in V 20 and V 30 -induced by the probabilistic nature of the emulators-through to the normalized volume: where Here,V i ( , ) and̂2 i ( , ) are the mean and variance of the predictive distribution of the pressure i volume emulator. This distribution is represented by the branch leading toṼ e 20 in Figure 3, where a dependence of the normalized volume on the material parameters is introduced through the underlying cardio-mechanic model.
We now have two different distributions forṼ 20 : Kl (Ṽ 20 ), from fitting the Klotz function to the ex vivo data, and e (Ṽ 20 ) from the cardio-mechanic emulators. Intuitively, we want to reward cardio-mechanic parameters for which the agreement between these two distributions is high, and penalize parameters for which the agreement is poor. There are different ways to quantify the agreement between the distributions. In the present study, the so-called 'product of experts' approach is used, which is widely used in the Machine Learning literature; see e.g. Calderhead et al. (2008) and Dondelinger et al. (2013). In the present case, it leads to a Gaussian convolution integral with a closed-form solution. Define: inserting Equations (14) and (15) into (17) and using standard results for Gaussian integrals, we obtain: where 2 K is the variance of the distribution in Equation (14). This quantity, for which a larger value represents better compliance with the Klotz curve, is represeted by the dashed line connecting the two ex vivo reference volume terms in Figure 3. Note that while  interpreted as a function of is a normal distribution,  interpreted as a function of the cardio-mechanic parameters is a complex non-linear function where: does not have a closed-form solution. We obtain, for a given LV geometry representation , the prior distribution of via normalization: where MCMC sampling techniques have to be applied due to the intractability of Z. We have now obtained a prior distribution of the cardio-mechanic parameters by fitting an empirical law for the pressure-volume dependence to ex vivo data from the literature, combining it with the cardio-mechanic model via the emulation of the LV cavity volumes at two ex vivo reference pressures -at 20 mmHg and 30 mmHg; see (16)-and taking explicitly the uncertainty predicted by the emulator into consideration.

Method B: The empirical Klotz prior
The first Klotz model was developed under the framework of the normalized Klotz curve, incorporating directly the relationship between normalized volume and pressure in Equation (9) through a prior that measures the adherence of the predicted LV dynamics to the Klotz model. An alternative model, referred to as the empirical Klotz prior, will now be considered where the level of fit with the Klotz curve is measured as a function of the measured LVV. Recall that the measured LVV is the cavity volume at end of diastole, measured from CMR scans and denoted as y * 0 in Equation (3). The empirical Klotz prior formulation is more in line with the use of the Klotz curve in the literature, which usually depends on an in vivo volume measurement.
Consider a rearrangement of the original Klotz function from Equation (9), where P ∼  (0, (P)) (see (10)) arises from an assumption of heteroscedastic Gaussian noise on the Klotz data:Ṽ This equation expresses the volume at pressure 30 mmHg as a function of the volume at some pressure P. Making the volume at pressure 30 mmHg the dependent variable is a natural consequence of the Klotz normalization, where V 30 appears in the denominator. Choosing P = 20 mmHg and plugging in the volume prediction from an emulator will result in a prior the same as the one from Section 3.1. Instead, P = 8 mmHg will be used, which means plugging in the LVV from the measured dataset, denoted as y * 0 . From the above relationship, and with reference to Equation (10), we have: Similar to the model from Section 3.1, ex vivo information can be introduced to the empirical Klotz model by inferring the parameters , and (P) conditional on the data in Figure 2 (recall from Equation (11) that estimating (P) involves inferring the parameters c 1 , c 2 and c 3 ). For estimation, the priors on the parameters were as follows: , ∼ u(0, 10) and c 1 , c 2 , c 3 ∼ u(−100, 100).
Using the posterior samples of and , the uncertainty in and can be propagated through to V 30 via the delta method (see the final section of the SM) to give a distribution on V 30 : with parameters: where the variable V 8 is replaced by the measured LVV y * 0 . The suitability of the approximation with the delta method was assessed by comparing with a distribution obtained using the full MCMC sample of and , which was found by evaluating (21) at the MCMC samples. The comparison is found in Figure 4, where we see good agreement between the two distributions.
In the clinical setting we must also account for noise in the LVV measurement y * 0 . Assuming that y * 0 is the mean of a distribution on measured LVV values with standard deviation 8 = 5 ml (this value was obtained from an empirical inter-operator study by Lazarus (2022)), the uncertainty on the volume measurement is propagated through the model via the delta method and 2 Kl includes the extra term F I G U R E 4 A visual comparison of V 30 distributions. The first distribution is obtained by evaluating (21) over MCMC samples of and , then using a Gaussian kernel density estimator with the bandwidth selected according to Silverman's rule to produce the density given by the orange dashed line. The delta method, which was used to evaluate the prior densities in this study, was used to produce the density shown by the solid blue line [Colour figure can be viewed at wileyonlinelibrary.com] F I G U R E 5 Visualizing the model with the empirical Klotz prior. The node containing y * 0 represents the measured left ventricle (LV) cavity volume at end diastole and y * 1∶ is the corresponding circumferential strain measurements. The part of the graph leading to these nodes represents the in vivo noise model, where we match the simulator predictions with the in vivo measured volume (y * 0 ) and strains (y * 1∶ ), which are assumed to be obtained at end diastolic pressure 8 mmHg. Based on the cardiac model, an emulator has been trained to predict the volume at pressure 30 mmHg and the predictive distribution provides a distribution on the volume, indicated by V e 30 . Inferring the Klotz relationship from the ex vivo data, y ex , and estimating V 30 from the in vivo volume using (21) gives another distribution on V 30 , denoted by V Kl 30 in the diagram. Using a product of experts approach, we assess the similarity of the two different distributions, indicated by the dashed line. Although the formulation of this model is different from the nonempirical Klotz prior, including a dependence on both and in a rearranged expression for predicting V 30 from the in vivo volume (compare (21) with (12)), the important difference with the model in Figure 3 is the introduction of a connection between the measured volume and the Klotz curve model, represented by the horizontal arrow exiting the node y * 0 which, unlike Kl (⋅), depends explicitly on the material parameters. This distribution is represented by the right branch of Figure 5 (leading to V e 30 ), which shows its dependence on the material parameters and the cardio mechanical model. We now have two distributions on V 30 : Kl (V 30 ) from the Klotz function in Equation (23) and the other, e (V 30 ), from the V 30 emulator with mean 30 ( , ) and variance 2 30 ( , ). Considering a statistic similar to that from Equation (17): we obtain: ) .
Letting Z = ∫  EK ( , )d , we retrieve a prior distribution: where the intractability of the denominator, Z, does not affect the inference if we use MCMC.

Klotz priors with model mismatch
Reconsidering Figure 2, we see that there is fairly substantial variance in the data. Due to the restrictive nature of the model in Equation (12) and the large quantity of data, much of the variation is inferred as noise during inference instead of an increased uncertainty in the parameters. This could lead to a slight underestimation of the true variation in the Klotz curve so an alternative form of the Klotz prior will be considered where this noise variance is included in calculations. The derivation of this prior in the nonempirical and empirical cases is outlined in sections 5.1 and 5.2 of the SM respectively. From this point forward, the priors implemented without the inclusion of this noise variance will be referred to as the Klotz prior without model mismatch and the prior with the noise variance included will be referred to as the Klotz prior with model mismatch.

IMPLEMENTATION
Codes to reproduce all results in this paper can be found at: https://github.com/lazarusal/klotzcodes

Construction of the in clinic emulators
We describe two in-clinic emulators, both described in Table 2. Our multi-geometry emulator emulates the cardio-mechanical model, taking as input an LV geometry at early diastole and the four material property parameters, and predicting the same outputs as the cardio-mechanical model (circumferential strains and LV cavity volume at end diastole). The LV geometry is internally mapped to a lower dimensional space defined by the leading principal components obtained from our reference population (see SM, section 1 for more details on the reference population and section 3 of the SM for more details on using PCA in the surrogate model). This is required to reduce the high-dimensional LV configuration space (about 17 k) to a dimension that can be more densely covered for training the emulator. We chose 5 PC components, based on our previous study (Romaszko et al., 2019), leading to a configuration space of dimension 9 overall (5 PCs components and 4 cardio-mechanic parameters). This space was covered with a space-filling design generated from a Sobol sequence, with a training set size of 16,000 (to train the emulator) and an independent validation set of slightly over 500 points to monitor potential overfitting, triggering early stopping of training (machine learning parlance for parameter adaptation

Generalized Klotz emulator
Emulator constructed over an input space of four material properties, early diastolic cavity volume and early diastolic wall volume In clinic prediction of volumes at ex vivo pressure ranges for implementation of Klotz priors

Multi geometry emulator
The emulator of the cardio mechanic model, generalized over five dimensional LV geometry representations obtained with PCA In clinic efficient approximation of the simulator model during an iterative optimisation process). Details on the selection of our emulation method can be found in section 2 of the SM, while the training procedure is discussed in section 4 of the SM. In addition, we generated two in-clinic emulators for predicting the volumes V 20 and V 30 , required for evaluation of the Klotz priors and referred to as generalized Klotz emulators. The input space here is six dimensional: four dimensions for material properties, one for cavity volume and one for wall volume at early diastole.

Exploring the effects of LV geometry approximation
Using 6 different LV geometries that were excluded from the original PCA analysis during emulator training, we will assess the ability of the empirical and nonempirical Klotz priors to reduce the error incurred by the LV geometry approximation in the emulator, as well as the intrinsic identifiability problems in the model. To do so, we require the gold standard parameters so we will generate synthetic data that are consistent with the Klotz priors to assess the improvement in inference accuracy that can be obtained under the assumption that the Klotz law is physiologically realistic (which we base on the evidence provided in the literature cited in Section 2.3). In the generation of these test data, we fix 1 and 4 at population based averages (this corresponds to fixing them equal to 1). There will be five different parameter configurations for each LV geometry, giving 30 test cases in total. Further details on the sampling procedure for test data generation can be found in section 6 of the SM.

Measured data study
As a second study, we consider the task of inferring the material properties from measured CMR scans. In this case, we do not possess the ground truth parameters but we can instead quantify the effect the Klotz curve priors have on cardio-mechanic parameter estimation and assess the physiological plausibility of this effect. For five different healthy volunteers, we choose four positions of short-axis cine images from the base to the middle ventricle for strain measurements. A B-spline deformable registration approach is used for estimating segmental averaged circumferential strains at each location (Mangion et al., 2016). Following the clinical convention, 6 peak diastolic strains are calculated at each short-axis position with respect to the end-diastolic phase. The end-diastolic LV cavity volume is calculated from short-axis cine images at the end of diastole after manually delineating the endocardial boundaries. The LV cavity volume at early diastole is calculated from the reconstructed LV geometry using the MATLAB function convhull. In total, we obtain 24 peak diastolic strains and 2 LV volumes (one at early diastole, and the other at end of diastole).

Simulations
Sampling of the material properties is carried out using Hamiltonian Monte Carlo with NUTS with dual averaging step size adaptation (Betancourt, 2017;Hoffman & Gelman, 2014). The likelihood function is evaluated using the statistical emulator so no forward simulations are required F I G U R E 6 Comparing the shapes of the nonempirical and empirical Klotz priors, from Equations (18) and (28) respectively, (top row) and the RSS functions (bottom row) for the in vivo data: nonempirical (evaluation of (18)) on the left, and empirical (evaluation of (28)) on the right. The orientation of the mode in the empirical Klotz prior (top right), along with identifiability of 2 based on the in vivo data, should allow for a reduction in the uncertainty of the parameters in the model. In the case of the nonempirical Klotz prior (top left), the high density ridge is aligned with the 3 direction. In each plot, the horizontal and vertical axes show the entire physiological range for the parameters, which is [0.1, 5] [Colour figure can be viewed at wileyonlinelibrary.com] during inference. We use a burn-in phase of length at least 5000, continuing until a potential scale reduction factor (Gelman & Rubin, 1992) of 1.1 is reached and then sample until an effective sample size (see section 11.5 of Gelman et al., 2013) of 100 is achieved. The burn-in is then discarded and the sample is used for analysis of the different methods.

Initial prior exploration with high fidelity emulators
First, we wish to investigate the ability of the Klotz priors to constrain the space of relevant constitutive parameters. The top row of Figure 6 provides visualizations of the prior surfaces, obtained by evaluation of (18) for the nonempirical Klotz prior and (28) for the empirical Klotz prior. Along with these, we plot the corresponding residual sum of squares (RSS) functions obtained based on the in vivo data in the bottom row of the figure. As we have discussed before, for the estimation of the cardio-mechanic parameters the data extracted from the CMR scans are informative to different degrees, generally allowing us to accurately estimate 2 but not 3 . 3 In light of this, the contour plot in the top left panel of Figure 5 suggests that with the nonempirical Klotz prior we will achieve little improvement in our inference, in that on estimating 2 from the CMR data, the conditional probability of 3 conditional on 2 remains diffuse. The problems witnessed in the nonempirical Klotz prior are largely alleviated by consideration of the empirical Klotz prior on the right of Figure 6. Here we introduce the measured EDV, meaning that only V 30 is predicted using an emulator (see left branch of Figure 5). We see that the ridge in the prior distribution is re-oriented such that conditioning on 2 (which can usually be easily estimated from the data as shown in the bottom right panel of this figure) substantially reduces the uncertainty for 3 . The most likely reason for this improvement is the fact that one of the two LV volume emulators required for the nonempirical Klotz prior has become obsolete as a consequence of explicitly including a measured LV volume. This in particular avoids the correlation problem mentioned in the previous paragraph.
Combining the priors from Figure 6 with a likelihood function obtained using the multi geometry emulator, we sample from the two dimensional posterior of 2 and 3 while fixing 1 and 4 to their true values. As a benchmark, we compare with posterior samples when adopting a non-informative  (0.1, 5) prior in each dimension. The posterior density, when using the nonempirical Klotz prior, is shown in Figure 7a. For data generated in the way described in Section 4.2, the nonempirical Klotz prior appears to reduce the uncertainty in the parameters compared with the non-informative prior in Figure 7c but some bias is present in the posterior samples. Using the empirical Klotz prior, we observe reduced uncertainty in 3 -shown in Figure 7b-compared with the non-informative prior case that is presented in Figure 7d. In this case, bias is not introduced to the posterior samples. This initial study demonstrates the potential the empirical Klotz prior has for improving cardio-mechanic parameter inference.

Sampling in four dimensional space
A key question of interest is the following: can the introduction of ex vivo information help to alleviate any error that is introduced to the inference problem by the LV geometry approximation while reducing uncertainty in the parameters? With test data generated as described in Section 4.2.1, we will answer this question with ex vivo information introduced by the generalized Klotz emulators and the simulator model approximated by a multi geometry emulator, both described in Table 2. The results will be assessed both in parameter space and in stress-stretch space. Beginning with parameter space, we will evaluate the inference results using squared error, which is defined as follows: wherêi is an estimate of material parameter i and̃i is the ground truth value of this parameter. For each of the six LV geometries, we have generated synthetic data from five different material parameter configurations. Using MCMC, we obtain a sample from the posterior distribution of the material parameters for each of these 30 test cases, providing a distribution of squared errors for each of the four material parameters. Given that the performance is measured in absolute terms via the squared error, we choose to combine the results, separately for each LV geometry, into one summary in the boxes of Figure 8. That is to say that, in total, each box provides the distribution of SEs obtained from 5 (the number of parameter configurations per LV geometry used for synthetic data generation) × 100 samples of the material parameters, where each subset F I G U R E 7 Kernel density estimates obtained using samples from the posterior of the material properties when ex vivo information is introduced in (a) the nonempirical Klotz prior (b) the empirical Klotz prior. The red dashed lines show the true parameter values. Panels (c) and (d) show posterior density plots obtained from samples when assuming a uniform prior on each of the parameters. In two dimensional plots, 2 is on the horizontal axis and 3 is on the vertical axis, with the entire physiological range plotted. For the 1D density plots, the horizontal axis shows the entire physiological range [Colour figure can be viewed at wileyonlinelibrary.com] of size 100 is the posterior distribution of squared errors corresponding to different test parameter configuration. For LV geometry i, these distributions are labelled as Gi, with separate plots for each of the four material parameters. We consider the empirical Klotz prior in Figure 8. The top row shows the performance when not accounting for the mismatch in the Klotz model (discussed in Section 3.3), compared against results with a non-informative prior on the four parameters. We observe improved inference results for 2 and 3 while sampling of 4 is largely unchanged. Information from the in vivo data is already enough to identify 1 , so we observe little difference in this parameter. At the bottom of Figure 8, we observe less improvement upon inclusion of the Klotz model mismatch term. This is exactly as expected given the weaker form of the prior distribution. Figure 9 displays the squared error in parameter samples when adopting the nonempirical Klotz prior. Similar to the empirical Klotz prior case, the inference for 2 and 3 is improved by introducing the Klotz information but this time the change is less pronounced. Again, we experience marginal improvement in 1 while 4 is unchanged. These results are mirrored in the F I G U R E 8 Squared error in samples from the posterior distribution of the material properties when using the empirical Klotz prior without the Klotz model mismatch term as discussed in Section 3.3. in the top row and with this term in the bottom row. A larger version is available in section 7.1 of the SM. The boxplots are arranged into the six different left ventricle (LV) geometries, denoted by Gi for the ith geometry, with samples obtained while adopting the empirical Klotz prior and a non-informative prior shown side by side in each pair. For each LV we have five parameter configurations and the squared errors from the 100 posterior samples for each configuration are combined into one distribution of 500 squared errors. The physiological range for each parameter is [0.1, 5] bottom row, where the Klotz model mismatch term was included in the prior. We note that the lack of improvement observed with the nonempirical Klotz prior could be partially attributed to errors observed in the approximation of the EDV. In the empirical Klotz prior this approximation is only required for the volume at pressure 30 mmHg, but in the nonempirical Klotz prior we make an approximation twice since V 20 is also required (see section 7.2 of the SM for more details).
In light of the superior results with the empirical Klotz prior, we now assess this method in more detail. For each test case, we have obtained a sample from the four-dimensional posterior distribution of the model parameters from Equation (2). These parameters physiologically characterise the tissue stiffness, which is more easily viewed in the stress-stretch space, that is, by virtually stretching a myocardial strip along a specific direction uni-axially (Gao et al., 2015). This provides a metric with which to assess the health of the cardiac tissue that could be used in the clinical setting. Note that stress here means the Cauchy stress (Holzapfel, 2000, chapter 3), which is the actual stress that the myocardium experiences. From our parameter samples, we obtain a distribution of stress-stretch curves with one curve for each sample from the posterior distribution of the material parameters. Similarly, knowledge of the ground truth parameters provides a ground-truth stress-stretch curve for each test case. We can measure the distance between each sampled curve and the ground truth curve at each stretch level, providing a distribution of distances that measures the accuracy with which the sampled material parameters can replicate the tissue behaviour of the ground truth parameters at a given stretch. In Figure 10, we determine the proportion of the 30 test cases (recall we have 6 LV geometries each with 5 test parameter configurations) where the distribution of curves at a given stretch, obtained from material parameters sampled under an assumption of the empirical Klotz prior, is significantly closer than the distribution of curves obtained from material parameters sampled under an assumption of a uniform prior, where closeness is measured by distance to the ground truth stress-stretch curve. A significant difference between the distributions of distances, which is assessed via the one-sided Mann-Whitney U Test at a 5% level (Mann & Whitney, 1947), signifies that the stress-stretch curves obtained from the samples under an assumption of the empirical Klotz prior are significantly closer to the ground truth stress-stretch curve than the stress-stretch curves obtained under an assumption of a uniform prior on the material parameters (note that a lack of significance does not imply a deterioration). We repeat this same procedure-checking for significant improvements in the stress stretch curve estimation-to check for improvement when using the uniform prior with the proportion of significant results also plotted in Figure 10. At the majority of stretch locations we observe an improvement upon introduction of the empirical Klotz prior but we do notice a dip in the proportion of significant improvements around stretch 1.12 for the myocyte direction and 1.2 for the sheet direction. Looking at the number of significant improvements when adopting the uniform prior we see a reverse of the pattern but the proportion of significant improvements with the uniform prior is low, reaching a peak of around 0.3 in the no model mismatch case. The corresponding stretch, where proportion of improvements with the empirical Klotz prior is low, is at a level where behaviour of the myocardium can be inferred from F I G U R E 10 Using stress-stretch curves to compare the performance of the empirical Klotz prior, with and without model mismatch as discussed in Section 3.3, and a uniform prior. First we tested for improvement in two different cases: 1. if the empirical Klotz prior without model mismatch improves upon a unif(0.1, 5) prior (K no mismatch) and 2. if the empirical Klotz prior with model mismatch improves upon a unif(0.1, 5) prior (K mismatch). The plots show the proportion of test cases where we find a significant improvement in our inference when using the empirical Klotz prior, assessed by comparing the distribution of distances of the sample of stress-stretch curves obtained under an assumption of the different priors. To complete the study, we also tested improvement in the other directions: 1. if a unif(0.1, 5) prior improves upon the empirical Klotz prior without model mismatch (U no mismatch) and 2. if a unif(0.1, 5) prior improves upon the empirical Klotz prior with model mismatch (U mismatch), assessing improvement as before. The online version is in colour [Colour figure can be viewed at wileyonlinelibrary.com] the in vivo data alone. Moreover, the fact that we do see a few cases where the inference is better with a uniform prior is not surprising since the effect of the Klotz prior distribution is to act as a regularizer, potentially introducing a small bias to the parameter estimates at the in vivo location while allowing for improved estimation in the ex vivo domain. This will be further discussed in the discussion section.

Inferring parameters from real data
Finally, we consider the task of inferring material properties from the real measured data of five healthy volunteers, outlined in Section 4.2.2. This time, the strains and LV volume are subject to measurement error, inevitably affecting the parameter estimation both with and without the Klotz prior (see section 7.4 of the SM for an additional exploration of parameter estimation in the presence of different levels of noise). Results for two of these are shown in Figure 11 (one subject in (a) and one in (b)). In the final section of the SM we provide plots with the full sample of curves for all five real data cases. Again, we only consider the empirical Klotz prior and, given that we use real data, we do not possess ground truth material properties so inference will only be

F I G U R E 11
Stress-stretch samples for two real test cases (left and right). From a sample of material parameters, obtained from the posterior distribution of the material parameters for each real test case, we obtain a sample of stress-stretch curves. We then plot these stresses as points at the different stretch locations along the myocyte (top) and sheet direction (bottom). We compare our inference with prior knowledge expressed through the uniform prior on material parameters (grey dots) and the empirical Klotz prior on the material parameters (white dots). The physiological range for the stress is less than 100 kPa assessed in stress-stretch space. In the distributions shown in Figure 11 we see that with the introduction of the Klotz prior without Klotz model mismatch, the stress in the upper stretch region is reduced, signifying a softer material than is described when we do not introduce the Klotz curve. This is true both along the myocyte (top row) and along the sheet direction (bottom row). This can be partially explained by the large uncertainties in the strain measurement, and its underestimation will lead to a much stiffer property prediction (see final paragraph of section 1 of the SM for more detail). Without the Klotz prior, mismatch between the measured strains and volumes can lead us to sample material properties that give smaller V 30 than the observed V 8 , violating the 'Frank-Starling' law in the heart (Moss & Fitzsimons, 2002). Introducing the Klotz information helps to alleviate this problem by balancing the weights between the measured strains and volumes. Thus the material property prediction with Klotz information will be physiologically more realistic within a wide range of LV pressures covering from in vivo to ex vivo.

DISCUSSION AND CONCLUSION
The work in this paper has adopted an approximation to the in vivo likelihood function in the form of a statistical emulator. Currently, statistical emulation provides the most hopeful route to implementation of biomechanical models in clinical decision making. The emulator implemented in this work takes a further step towards this goal by generalizing a fixed LV geometry emulator with the introduction of a low dimensional LV geometry component to the input space, allowing the statistical model to interpolate between behaviours of different LVs. This provides a computationally efficient approximation to the true model, but with the introduction of some error due to the required LV geometry approximation. These errors exacerbate the weak-identifiability of F I G U R E 12 Change in the shape of the simulated LVV function over parameter space as the end diastolic pressure P varies. The horizontal axis is 2 and the vertical axis is 3 . Increasing P from 8 to 30 mmHg aligns the contours more with the 2 axis. The final plot shows that by normalizing the volume at P = 20 mmHg with the volume at P = 30 mmHg we end up with a completely new alignment of the surface that is often aligned with the 3 direction the parameters in the Holzapfel-Ogden model. More fundamentally, the weak identifiability of the parameters is a result of the absence of information from the inference problem about the high stretch behaviour of the tissue. In the parameterization considered in this paper, this is of particular relevance to the parameter 3 .
A key motivation for the work performed in this paper is found in Figure 12, which shows the simulated volume as a function of 2 and 3 . This plot is also of more general relevance, giving a more precise explanation for the improvements seen with high pressure data. At low pressures (8 mmHg in this case) the contours of the volume function are aligned with the 3 axis. If we were to form a loss function based on this simulated function then the value of the loss function would be similar along the contours, and identifiability of 3 would be poor. As we increase the pressure, first to 20 mmHg and then 30 mmHg, the contours begin to align with the 2 axis, which allows for better estimation of 3 . This is the effect that is included in the empirical Klotz prior and justifies its success in the inference scheme. The final plot shows why the nonempirical prior was less successful. By normalizing the simulated V 20 with the simulated V 30 we have an alignment of the loss function that varies over the input space, but is largely in the direction of the 3 axis.
The first Klotz paradigm presented in this paper considered the normalized form of the Klotz curve in a prior named the nonempirical Klotz prior. The presence of the ratio of V 20 and V 30 in the Klotz function in Equations (12) via (9) leads to identifiability issues, as a result of a potential correlation between the two volumes. This observation, visualized in Figure 12, is the main motivation for our development of the empirical Klotz prior, where we express V 30 as a function of the empirical EDV measurement at pressure 8 mmHg. We have shown that in the context of inference with a statistical emulator, the inclusion of an in vivo EDV measurement in the Klotz prior proves to be important for achieving a performance improvement in the inference scheme.
Our evaluation has been based on synthetic data that are consistent with the empirical Klotz law, and real data extracted from the CMR scans of five healthy volunteers.
In reality, any mismatch between the behaviour of the in vivo and ex vivo left ventricle-or the consideration of a left ventricle that violates the Klotz curve relation-will lead to a bias in the parameter estimates, particularly those parameters that relate to the behaviour of the tissue in the high stretch regime such as 3 in our parameterization. It is important to note, however, that the effect of the Klotz prior is also to reduce the variance in the posterior samples. Considering the bias-variance decomposition of the generalization error of the model, the prior will act to reduce the variance term while incurring a small increase in the bias (see section 3.2 of Bishop (2006) for more details). In this sense, the overall effect of the Klotz prior will be to improve the generalization error in much the same way as any regularization method from statistics.
There are some limitations in our methodology. Currently, the statistical emulator does not account for differing patient pressures in the input space. This limits the applicability of the model to healthy volunteers, since this pressure assumption is only valid for this population. It is likely that the generalization of the input to different end diastolic pressures would increase the uncertainty in the material properties since the uncertainty in the pressure measurement must also be accounted for in the parameter inference scheme. A second limitation of the work in this study is the use of the LV geometry at early diastole as the unloaded state (with zero pressure in the LV cavity). This approximation, which is required for evaluation of the Klotz priors in Equations (18) and (28) where both i ( , ) terms depend on V 0 , also introduces error to the in vivo likelihood for measured data and future work will aim to estimate the unloaded configuration, as has been done in the recent literature (Wang et al., 2020). Another issue we encounter is the informativeness of the measured data, particularly with reference to parameters 3 and 4 . The absence of radial and longitudinal strain measurements from the likelihood function is partially responsible for the identifiability issues that we observe in the model. An interesting future line of research would be to introduce these extra strain data to the likelihood function and infer the material properties when adopting a Klotz prior distribution. The combination of the two should lead to better parameter identifiability than is observed in the case of a non-informative prior with a likelihood that only depends on circumferential strain data. However, obtaining high precision measurements of these data is currently not possible.
With the implementation of the Klotz curve, we are able to partially alleviate errors introduced by the approximated in vivo likelihood function while reducing uncertainty in the weakly identifiable 3 parameter. We have shown for synthetic data that are consistent with the empirical Klotz law that this ex vivo information leads to a significant bias and uncertainty reduction for one of the two Klotz prior models proposed (the empirical Klotz prior). For measured datasets, we have shown the ability of the Klotz curve to alter the inferred behaviour of the myocardium in the ex vivo domain, reducing our uncertainty about this behaviour. Expansion on the current work is currently underway, with attempts to improve identifiability of the model parameters by improved LV geometry representation or improved HO law parameterization. Another interesting course of study is the inclusion of end-diastolic pressure in the emulator input space, allowing further generalization in the clinical setting.