Fast Parameter Inference in a Biomechanical Model of the Left Ventricle using Statistical Emulation

TL;DR

This paper introduces a fast, emulation-based approach to estimate myocardial properties from MRI data, significantly reducing computation time for clinical use in personalized heart modeling.

Contribution

It presents a novel application of statistical emulation to biomechanical LV models, enabling rapid parameter inference suitable for clinical settings.

Findings

Emulation methods drastically reduce computation time.

Comparison of two emulation strategies shows their effectiveness.

Different interpolation and loss functions impact inference accuracy.

Abstract

A central problem in biomechanical studies of personalised human left ventricular (LV) modelling is estimating the material properties and biophysical parameters from in-vivo clinical measurements in a time frame suitable for use within a clinic. Understanding these properties can provide insight into heart function or dysfunction and help inform personalised medicine. However, finding a solution to the differential equations which mathematically describe the kinematics and dynamics of the myocardium through numerical integration can be computationally expensive. To circumvent this issue, we use the concept of emulation to infer the myocardial properties of a healthy volunteer in a viable clinical time frame using in-vivo magnetic resonance image (MRI) data. Emulation methods avoid computationally expensive simulations from the LV model by replacing the biomechanical model, which is defined in terms of explicit partial differential equations, with a surrogate model inferred from simulations generated before the arrival of a patient, vastly improving computational efficiency at the clinic. We compare and contrast two emulation strategies: (i) emulation of the computational model outputs and (ii) emulation of the loss between the observed patient data and the computational model outputs. These strategies are tested with two different interpolation methods, as well as two different loss functions...

Tables

Table 1. Table 1 : Table giving the median (1st quartile, 3rd quartile) of the mean squared error (in parameter-space) in the prediction of all the model parameters. The considered interpolation methods (Interp. Meth.) are Low-Rank GPs and Local GPs, the target of the emulation (Emulation Target) is either the model output or the loss, and two loss functions are compared: Euclidean and Mahalanobis. The method with the best predictive performance, the output emulation method with local GP interpolation and the Euclidean loss function, is given in bold.

Interp. Meth.	Emulation Target	Euclidean	Mahalanobis
Low-Rank GP	Output	0.0048 (0.0012,0.0107)	0.0030 (0.0011,0.0062)
Low-Rank GP	Loss	0.6814 (0.2222,1.5234)	0.0113 (0.0041,0.0377)
Local GP	Output	0.0001 (0.0000,0.0003)	0.0009 (0.0003,0.0022)
Local GP	Loss	0.2201 (0.0588,0.6777)	0.0013 (0.0002,0.0063)

Equations

$$\begin{array}{rl}\Psi & =\frac{a}{2b}\{\text{exp}[b(I_1-3)]-1\}+\sum _{i\in \{f,s\}}\frac{a_i}{2b_i}\{\text{exp}[b_i(I_{4i}-1{)}^2]-1\}\\ & +\frac{a_{\text{fs}}}{2b_{\text{fs}}}[\text{exp}(b_{\text{fs}}{I}_{\text{8fs}}^2)-1]+\frac{1}{2}K(J-1{)}^2,\end{array}$$

$$I_1=\text{mbox}trace(\mathbb{C}),I_{\text{4f}}={\mathrm{f}}_0\cdot (\mathbb{C}{\mathrm{f}}_0),I_{\text{4s}}={\text{s}}_0\cdot (\mathbb{C}{\text{s}}_0),I_{\text{8fs}}={\mathrm{f}}_0\cdot (\mathbb{C}{\text{s}}_0)$$

$$a$$

$$a_{\text{f}}$$

$$b_{\text{f}}$$

$$a_{\text{fs}}$$

$$\ell (\mathbf{\theta }|\mathsf{m},{\mathbf{y}}_0)=\alpha d(\mathsf{m}(\mathbf{\theta }),{\mathbf{y}}_0)+Z,$$

$$\hat{\mathbf{\theta }}=\arg \underset{\mathbf{\theta }}{\min }\ell (\mathbf{\theta }|\mathsf{m},{\mathbf{y}}_0),$$

$$\ell (\mathbf{\theta }|\hat{\mathsf{m}},{\mathbf{y}}_0)=\frac{1}{2{\sigma }^2}\Vert \hat{\mathsf{m}}(\mathbf{\theta })-{\mathbf{y}}_0{\Vert }^2+Z(\sigma ).$$

$$\ell (\mathbf{\theta }|\hat{\mathsf{m}},{\mathbf{y}}_0)=\frac{1}{2}(\hat{\mathsf{m}}(\mathbf{\theta })-{\mathbf{y}}_0{)}^{\top }{\mathbf{\Sigma }}^{-1}(\hat{\mathsf{m}}(\mathbf{\theta })-{\mathbf{y}}_0)+Z(\mathbf{\Sigma }),$$

$${\hat{\mathbf{\theta }}}_0=\arg \underset{\mathbf{\theta }}{\min }\ell (\mathbf{\theta }\mid \hat{\mathsf{m}},{\mathbf{y}}_0).$$

$${\mathcal{D}}_{\ell }=\{({\mathbf{\theta }}_n,{\ell }_n):n=1,\dots ,N\},$$

$${\hat{\mathbf{\theta }}}_0=\arg \underset{\mathbf{\theta }}{\min }\mathbb{E}\{\hat{\ell }(\mathbf{\theta }|\mathsf{m},{\mathbf{y}}_0)\}.$$

$${\tilde{\mathbf{y}}}_j\sim {\beta }_j\text{1}+f(\mathbf{\theta })+f(\mathbf{\tau })+\mathbf{\epsilon }\text{mbox}forj=1,\dots ,J$$

$$\mathbf{Cov}[\mathbf{y}({\mathbf{\theta }}_i),\mathbf{y}({\mathbf{\theta }}_j)]=K({\mathbf{\theta }}_i,{\mathbf{\theta }}_j)\mathbf{A}$$

$$\mathrm{Cov}[y_k({\mathbf{\theta }}_i),y_l({\mathbf{\theta }}_j)]=K({\mathbf{\theta }}_i,{\mathbf{\theta }}_j)A({\mathbf{u}}_k,{\mathbf{u}}_l)$$

Full text

††journal: Journal of the Royal Statistical Society, Series C

1 Introduction

It is widely recognised that when integrated with in vivo data from cardiac magnetic resonance imaging (MRI), computational modelling of cardiac biomechanics can provide unique insights into cardiac function in both healthy and diseased states (Wang et al., 2015; Chabiniok et al., 2016; Gao et al., 2017a). For example, recent mathematical studies have demonstrated that passive myocardial stiffness is much higher in diastolic heart failure patients compared to healthy subjects (Xi et al., 2014). Similarly, myocardial contractility could be much higher in acute myocardial infarction (MI) patients than it is in healthy volunteers (Gao et al., 2017a). In particular, the myocardial passive properties not only affect left ventricular (LV) diastolic filling, but also influence the pumping function in heart chamber contractions (systole) through the ‘Frank-Starling’ law (Widmaier et al., 2016), the relationship between stroke volume and end diastolic volume.

In order to comprehensively assess LV function, it is necessary to determine passive myocardial stiffness. Traditionally myocardial passive properties can be determined by a series of ex vivo or in vitro experiments (Dokos et al., 2002). The widely used Holzapfel-Ogden (HO) constitutive law (Holzapfel and Ogden, 2009) can give a detailed description of the myocardium response in passive state, including the effects of collagen fibre structure. However, determining the material parameters of this model is challenging for clinical applications, as one can not perform invasive experiments as in Dokos et al. (2002). One possibility of estimating these parameters non-invasively is by cardiac MRI. The biomechanical model used in this study describes the LV dynamics during the diastolic filling process, starting from early-diastole and finishing at end-diastole, which is the point of maximum LV expansion. Both early and end-diastolic states can be measured by MRI. We can then compare, for a given patient, these measurements to the predictions from the biomechanical model, which defines the likelihood. The biophysical parameters defining the myocardial properties (as described by the HO law) can then be inferred in an approximate maximum likelihood sense using an iterative optimisation procedure, as discussed in Gao et al. (2015). In the context of mathematical physiology, this procedure is referred to as solving the inverse problem.

The inverse problem itself can be solved using a variety of methods and many studies have demonstrated that it is possible to estimate constitutive material parameters using in vivo measurements even with very complex constitutive relations (Guccione et al., 1991; Remme et al., 2004; Sermesant et al., 2006; Sun et al., 2009). However, because of the strong correlation among the material parameters and sparse noisy data, the formulated inverse problem is highly non-linear (Xi et al., 2011; Gao et al., 2015). Furthermore, determining the unknown parameters in this way is very time consuming, with the process taking days or weeks to converge, even with a modern multi-core workstation (Gao et al., 2015; Nikou et al., 2016). The primary reason for this is the high computational expense of simulating from the biomechanical model, which requires a numerical integration of the underlying partial differential equations with finite element discretisation. This procedure has to be repeated hundreds or thousands of times during the iterative optimisation of the material parameters.

As a result of the high computational costs of simulating the biomechanical model, estimating myocardial properties using a process which uses this model as a simulator is not suitable for real-time clinical diagnosis. A potential approach to overcome this problem is emulation (e.g. Kennedy and O’Hagan (2001); Conti et al. (2009); Conti and O’Hagan (2010)), which has recently been explored in the closely related contexts of cardiovascular fluid dynamics (Melis et al., 2017), the pulmonary circulatory system (Noè et al., 2017) and ventricular mechanics (Achille et al., 2018).

Emulation methods are far more computationally efficient as most of the computation can be done in advance, making the in-clinic diagnosis faster. With emulation approaches, we simulate a large number of samples at different parameter specifications in advance and use these simulations combined with an interpolation method to replace the computationally expensive simulator in the optimisation procedure. The choice of parameter combinations from which simulations are taken can be determined effectively using a space filling design, in this case produced by a Sobol sequence (Sobol, 1967), to spread the parameter combinations chosen in a way that aims to maximise the information about the simulator for a given number of simulations via several uniformity conditions. Optimising this design is an active research area (see e.g. Overstall and Woods (2017)), which is beyond the remit of the present paper though.

The work presented here is designed as a proof of concept study to assess the accuracy of alternative emulation strategies for learning the material properties of a healthy volunteer’s LV myocardium based on only non-invasive, in vivo MRI data. To that end, we use a patient-specific model with a fixed, patient-specific LV geometry, and focus on the statistical methodology for biophysical parameter estimation. Additionally, we use a reduced parameterisation of the HO law with the biomechanical model based on the work of Gao et al. (2015) in MRI data. Based on this approach, we compare different emulation strategies, loss functions and interpolation methods.

The first of the emulation approaches we have tested is based on emulating the outputs of the simulator, Section 3.3.1, in this case the simulated clinical data based on the described biomechanical model. Here, individual interpolators are fitted to each of the simulator outputs, using our chosen interpolation technique. We can then calculate the loss function between the predicted output of the individual models and the observed new data points from which we wish to learn the underlying myocardial properties. Minimising this loss function via a standard optimisation routine then produces estimates of the material parameters of the new subject. A variety of loss functions can be used within our emulation methods and we have compared two different ones here. The first of these is the Euclidean loss function, which assumes independence between outputs and the second is the Mahalanobis loss function (Mahalanobis, 1936) which allows for correlations.

The second emulation approach involves emulating a loss function rather than the outputs directly, Section 3.3.2, where again we use both the Euclidean and the Mahalanobis loss functions. For new MRI data, we calculate the loss, which quantifies the discrepancy between the model predictions and the data. Statistical interpolation is then used to obtain a surrogate loss function over the biophysical parameter space, which can be minimised with standard iterative optimisation routines.

In addition to testing these two emulation paradigms, we test two interpolation techniques based on Gaussian Processes (Rasmussen and Williams, 2006). The first of these is a low rank Gaussian Process (GP) emulation method, which uses the complete dataset for interpolation, but uses a low rank approximation in order to scale to high dimensions (Wood, 2003). The second method uses a local GP, where the interpolation is based on the $K$-nearest neighbours closest to the current values of the material parameters. Using a reduced number of training points from the simulations at each stage of the optimisation procedure and thereby lowering the computational costs is important, as due to the cubic computational complexity in the number of training points a standard GP would not be suitable for clinical decision support in real time.

In this work, we firstly compare different combinations of emulation methods, interpolation methods and loss functions in order to determine which method provides the best estimate of the material LV properties. We do this via a simulation study, Sections 5.1, 5.2, 5.3 and 5.4, using additional independent simulations from the simulator as out-of-sample test data. Knowledge of the true parameter values allows us to assess the accuracy of the different combinations of methods. We then test the best combination of methods on real MRI data from the healthy volunteer from which we have taken the LV geometry, Section 5.5, to assess the accuracy of biomechanical constitutive parameter estimation in a time frame suitable for clinical applications.

2 Left-Ventricle Biomechanical Model

The LV biomechanical model describes the diastolic filling process from early-diastole to end-diastole. There are multiple different models that can be used to describe this process and these are reviewed in detail in Chabiniok et al. (2016). The model used here is similar to those used in Wang et al. (2013) and Gao et al. (2015). The biomechanical model initially described in Wang et al. (2013) can be thought of consisting of 5 parts: initial discretised LV geometry, the constitutive law (the HO law), the constitutive parameters, the finite element implementation, and corresponding initial and boundary conditions. Linking this biomechanical model to patient MRI data can allow the inference of unknown material parameters describing heart mechanics, potentially leading to improved disease diagnosis and personalised treatments (Gao et al., 2017b).

The mathematical model takes 3 inputs: the initial discretised LV geometry constructed from MRI images at early-diastole (Section 2.1), corresponding initial and boundary conditions (Section 2.2), and constitutive parameters (Section 2.3). Based on these inputs, the mathematical model, implemented in ABAQUS333ABAQUS (Simulia, Providence, RI USA), simulates the diastolic filling process using the HO law (Section 2.3) and a finite element implementation (Gao et al., 2015). The output of the mathematical model then gives a model of the LV state at end-diastole, which can be compared to the corresponding in-vivo MRIs. These MRIs at end-diastole are used to measure circumferential strains taken at 24 locations444These are based on the American Heart Association definition as in Gao et al. (2015). and the end diastolic volume. These measurements can be compared against those generated by the biomechanical model for various constitutive parameters in order to learn the parameters associated with the volunteer from which the MRI were taken.

Each simulation from the mathematical model without parallelisation takes about 18 minutes on our local Linux workstation555Intel(R) Xeon(R) CPU, 2.9GHz, 32G memory, or around 4.5 minutes with parallelisation on 6 CPUs. Note that the 18 or 4.5 minutes are required for just a single parameter adaption step of an iterative optimisation, or a single addition to the emulator.

2.1 Initial discretised LV geometry

The initial discretised LV geometry can be obtained through constructing a 3D model based on the MRI scans (Wang et al., 2013). The scans consist of a series of 6 or 7 short-axis cine images which cover the ventricle666The MRI study was conducted on a Siemens MAGNETOM Avanto (Erlangen, Germany) 1.5-Tesla scanner with a 12-element phased array cardiac surface coil. Cine MRI images were acquired using the steady-state precession imaging protocol. Patient consent was obtained before the scan.. For each cardiac cycle there are usually around 35 frames from end-diastole to early-diastole. The images of the early-diastole are then used to create the initial discretised LV geometry, while the end-diastole images will provide the final measurements of the circumferential strains and the LV volume. To create the discretised LV model, the endocardial (inner) and epicardial (outer) boundaries of the LV are segmented from cine images at early-diastole as done in Gao et al. (2014b), e.g. Figure 1(a). A 3D model of the LV can then be constructed in Solidworks777Solidworks (Dassault Systems SolidWorks Corp., Waltham, MA USA), e.g. Figure 1(b). Finally, Figure 1(c) is constructed using a rule-based fibre-generation method, see Gao et al. (2014a), giving us the initial discretised LV geometry used in the biomechanical model. In the context of the present study, we consider this a fixed input and focus our work on developing parameter inference methods rather than a tool that can work for all possible subjects. Extensions to allow for different LV geometries is the subject of future work.

2.2 Initial and Boundary Conditions

The initial and boundary conditions, in particular LV pressure, play an important role in myocardial dynamics. Unfortunately, blood pressure within the cavity of the left ventricle can only be measured invasively, by direct catheter measurement within the LV cavity. Due to potential complications and side effects, these measurements are not available for healthy volunteers. We have therefore fixed the boundary conditions, including the pressure, at values considered sensible for healthy subjects, based on the work of Bouchard et al. (1971).

2.3 Constitutive Law

The final part of the biomechical model is the constitutive law for characterising the material properties of myocardium. In this study, we use the invariant-based constitutive law (Holzapfel and Ogden, 2009), based on the following strain energy function:

[TABLE]

in which $a$, $b$, $a_{\text{f}}$, $b_{\text{f}}$, $a_{\text{s}}$, $b_{\text{s}}$, $a_{\text{fs}}$ and $b_{\text{fs}}$ are unknown material parameters, $I_{1},I_{4i}$, and $I_{\text{8fs}}$ are the invariants corresponding to the matrix and fibre structure of the myocardium, which are calculated as

[TABLE]

in which $\bm{\mathrm{f}}_{0}$ and $\text{s}_{0}$ are the myofibre and sheet orientations, which are determined through a rule-based approach (Wang et al., 2013) and are known before the simulation (initial conditions). $\mathbb{C}$ is the right Cauchy-Green deformation tensor, defined as $\mathbb{C}=\mathbb{F}^{T}\mathbb{F}$, where $\mathbb{F}$ is the deformation gradient describing the motion of myocardium and hence how its shape changes in 3D with time. The term $\frac{1}{2}K(J-1)^{2}$ accounts for the incompressibility of the material, where $K$ is a constant ($10^{6}$) and $J$ is the determinant of $\mathbb{F}$. The HO law forms a major part of the biomechanical model, and the 8 constitutive parameters, $a$, $b$, $a_{\text{f}}$, $b_{\text{f}}$, $a_{\text{s}}$, $b_{\text{s}}$, $a_{\text{fs}}$ and $b_{\text{fs}}$, are unknown inputs into the model, which we wish to learn. The accuracy of parameter estimation for real data can be based on stretch-stress curves, as discussed in Section 1 of the online supplementary materials.

However, it has previously been found in Gao et al. (2015) that the 8 parameters are strongly correlated, which suggests that a model reduction is advisable to ensure identifiability. The authors further demonstrated that myofibre stiffness, the parameter most relevant for clinical applications, can be estimated from in vivo data with a reduced parameterisation; see Section 2 of the online supplementary materials. In fact, Hadjicharalambous et al. (2016) even estimated passive myocardial stiffness using a reduced form of the HO law with only a single unknown parameter. In the present study, similarly to Gao et al. (2015), we group the eight parameters of (1) into four, so that:

[TABLE]

where $\theta_{i}\in[0.1,5]:i=1,2,3,4$ are the parameters to be inferred from in vivo data, and $a_{0}$, $b_{0}$, $a_{\text{f0}}$, $a_{\text{s0}}$, $b_{\text{f0}}$, $b_{\text{s0}}$, $a_{\text{fs0}}$, and $b_{\text{fs0}}$ are reference values from the published literature (Gao et al., 2017a)888The reference values are, up to 2 decimal places: $a_{0}=0.22$, $b_{0}=1.62$, $a_{\text{f0}}=2.43$, $a_{\text{s0}}=0.56$, $b_{\text{f0}}=1.83$, $b_{\text{s0}}=0.77$, $a_{\text{fs0}}=0.39$, and $b_{\text{fs0}}=1.70$.. Our results obtained with this dimension reduction are consistent with the experimental results reported in Dokos et al. (2002).

3 Statistical Methodology

This section reviews the notion of a simulator and emulator, as well as establishing the notation that is used throughout the rest of the paper. It also provides details about the different emulation strategies that are going to be used in this paper, as well as the different interpolation methods considered.

3.1 Simulation

A simulator, $\bm{\mathsf{m}}$, is a mathematical model that relies on a computationally expensive numerical solution of the underlying systems equations. In the present study, the mathematical model is the soft-tissue mechanical description of the left ventricle based on the Holzapfel-Ogden strain energy function, as discussed in the previous section. The numerical procedure is the finite element discretization of the resulting partial differential equations. The inferential process, i.e. estimating the unknown inputs or parameters $\bm{\theta}_{0}$ underlying the observed clinical data $\bm{y}_{0}$, is computationally expensive and infeasible in settings where solutions are required within a short time frame, for instance in the context of clinical decision support. The prohibitive computational time that makes inference challenging is due to the time needed for a single (forward) simulation from the computational model, where by forward simulation we mean generating a (possibly multivariate) output $\bm{y}=(y_{1},\ldots,y_{J})=\bm{\mathsf{m}}(\bm{\theta})$ for a given parameter vector or input $\bm{\theta}$. In the context of the present study, $J=25$, and the outputs $y_{i}$ are the 24 circumferential strains and the LV volume at end of diastole, as predicted by the mathematical model.

Given our clinical data, $\bm{y}_{0}$, which are the measured circumferential strains and the end-of-diastole LV volume obtained from MRI, we can estimate the unknown parameter vector $\bm{\theta}_{0}$ by finding the corresponding input to the simulator which gives rise to an output which is as close as possible to the observed clinical data, $\bm{y}_{0}$. While our clinical data is assumed to come from the same data generating process $\bm{\mathsf{m}}$ for an unknown input, $\bm{\theta}_{0}$, in practise there will be a systematic deviation due to noisy measurement and model mismatch. The simplest approach to estimating the unknown input or parameter vector $\bm{\theta}$ is to choose the loss function as the negative log-likelihood:

[TABLE]

for a given metric function $d$ measuring the distance between a simulation $\bm{y}=\bm{\mathsf{m}}(\bm{\theta})$ and data $\bm{y}_{0}$, and some positive constants $\alpha$ and $Z$. We can then estimate the input to the model by minimising the true loss in (3):

[TABLE]

effectively giving us the maximum likelihood estimate of $\bm{\theta}$. This method becomes prohibitive if a single simulation exceeds a certain amount of time, as it does with the biomechanical model considered in the present work. The numerical procedure based on finite element discretization requires approximately 18 minutes for a single simulation, or 4.5 minutes with parallelisation on 6 CPUs on our computing system999Dual Intel Xeon CPU E5-2699 v3, 2.30GHz, 36 cores and 128GB memory. Any optimisation of the true loss, (3), would require the evaluation of the simulator at every iteration of the optimisation routine, potentially hundreds or thousands of times, with each iteration taking between 4.5 and 18 minutes. This is computationally limiting if we wish to use the method for clinical decision support in real time.

3.2 Emulation

An emulator is a statistical model that is a cheap and fast approximation to the true computational model (simulator), $\bm{\mathsf{m}}$, in this case the biomechanical model. It is used to replace the simulator in order to speed up both computations and inference, and it is also referred to as a meta-model (or surrogate model) as it represents a model of a model. An emulator can be built using any interpolation technique such as regression splines, polynomial regression, Gaussian processes, etc; see Section 3.4 for more details. Once a method has been chosen and the emulator has been fitted to the training data, we will denote it as $\hat{\bm{\mathsf{m}}}$.

In order to fit a statistical model and replace the simulator, we need training data from the simulator itself in the form of simulations $\mathcal{D}=\{(\bm{\theta}_{1},\bm{y}_{1}),\dots,$ $(\bm{\theta}_{N},\bm{y}_{N})\}=\{\bm{\Theta},\bm{Y}\}$. In the context of the present application, the input vectors $\bm{\theta}_{i}$ are the biomechanical parameter vectors discussed in Section 2.3. These inputs into the simulator, $\bm{\Theta}$, are chosen based on a space filling design, using Sobol sequences. These so-called low-discrepancy sequences are known to lead to improved convergence in the context of quasi-Monte Carlo; see e.g. Gerber and Chopin (2015). A more efficient coverage of the input space is possible using more advanced statistical design methods, as e.g. discussed in Overstall and Woods (2017), but these explorations are beyond the remit of present work.

The outputs of the simulator, $\bm{Y}$, are the resulting clinical values based on the assumed data generating process, $\bm{\mathsf{m}}$. In the present application, the output vectors $\bm{y}_{i}$ are the vectors of 24 circumferential strains and LV volume at end of diastole. Whilst generating large numbers of simulations is computationally expensive, this can be massively parallelised in advance and before the arrival of the patient at the clinic.

Previously, given the clinical data, $\bm{y}_{0}$, and a simulator, $\bm{\mathsf{m}}$, we could not estimate the unknown input, $\bm{\theta}_{0}$, using the loss function (negative log-likelihood) given in (3) fast enough for effective use within a clinical environment. This was due to the high simulation time required for each single input. Now, however, we can replace the true loss function in (3) with a surrogate loss function, $\ell$, based on an emulation method; see Section 3.3 for details. Minimisation of the surrogate loss (surrogate negative log-likelihood) for any metric function $d$ will be fast and suitable for real-time precision medicine, as it does not involve any simulation from the computationally expensive model.

We can use a variety of different metric functions within our surrogate loss, $\ell$. The most obvious of these is the Euclidean norm, $\|\hat{\bm{\mathsf{m}}}(\bm{\theta})-\bm{y}_{0}\|^{2}$. Under the assumption of independent and identically (iid) normally distributed errors (i.e. deviations of the clinical data from the emulator outputs) with zero mean and variance $\sigma^{2}$, the Euclidean loss function is equivalent to the negative log-likelihood, up to a scaling factor and an additive constant $Z(\sigma)$:

[TABLE]

An extension of the Euclidean loss which allows for a correlations between the outputs is the Mahalanobis loss function:

[TABLE]

which is equivalent to the negative log-likelihood of a multivariate Gaussian distribution with covariance matrix $\bm{\Sigma}$ up to a constant, $Z(\bm{\Sigma})$. To minimise the computational costs at the clinic, the covariance matrix is pre-computed from the training data, $\bm{\Sigma}=\mathrm{cov}(\bm{Y})$, and then kept fixed. Its main purpose is to allow for the spatial correlations between the 24 circumferential strains at different locations on the LV.

3.3 Emulation Frameworks

3.3.1 Output Emulation

Emulating the outputs of the simulator, the LV model, involves fitting multiple individual models, one for each of the $J$ outputs of the simulator, $\bm{\mathsf{m}}$. These outputs, $y_{j}:j=1,\ldots,J$, are fitted using the inputs of the simulator, $\bm{\Theta}$, with an appropriate interpolation method; see Sections 3.4.1 and 3.4.2. Given the multiple independent models, $\hat{\bm{\mathsf{m}}}=(\hat{\mathsf{m}}_{1},\ldots,\hat{\mathsf{m}}_{J})$, estimates of the parameter vector, $\hat{\bm{\theta}}_{0}$, can be found for any new set of outputs $\bm{y}_{0}$ by minimising the difference between $\bm{y}_{0}$ and $\hat{\bm{\mathsf{m}}}(\bm{\theta})$ with a loss function:

[TABLE]

The loss function, $\ell$, in (7) can take a variety of forms, including the Euclidean and the Mahalanobis loss functions given in (5) and (6). An algorithmic description of the output emulation method is given in Algorithm 1.

The advantage of emulating the outputs is that the statistical models can be fitted in advance, before the data have been collected from the clinic, meaning that when a patient comes into the clinic, an estimation of the biomechanical parameter vector $\hat{\bm{\theta}}_{0}$ can be carried out relatively quickly. The disadvantage is that multiple potentially correlated model outputs must be fitted, leading to higher computational costs at training time than emulating the loss function directly.

3.3.2 Loss Emulation

An alternative strategy is loss emulation. This entails direct emulation of the losses $\ell_{n}=\ell(\bm{\theta}_{n}|\bm{\mathsf{m}},\bm{y}_{0})$ rather than the simulator outputs $\bm{y}_{n}=\bm{\mathsf{m}}(\bm{\theta}_{n})$, for $n=1,\dots,N$. Given simulations $\mathcal{D}=\{\bm{\Theta},\bm{Y}\}$ and clinical data $\bm{y}_{0}$, it is possible to evaluate the loss function at each of the training points and data $\bm{y}_{0}$ and record the obtained score. To follow this approach we fit a single real-valued emulator to training data:

[TABLE]

where $\ell_{n}=\ell(\bm{\theta}_{n}|\bm{\mathsf{m}},\bm{y}_{0})$ is the loss function, for a given metric $d$, evaluated at the $n$th design point from the corresponding simulation output, $\bm{y}_{n}=\bm{\mathsf{m}}(\bm{\theta}_{n})$. The metric $d$ should be chosen according to the problem, and it can capture the correlation between the model outputs. Now it is possible to fit a single real-value emulator $\hat{\ell}(\bm{\theta}|\bm{\mathsf{m}},\bm{y}_{0})$ of $\ell(\bm{\theta}|\bm{\mathsf{m}},\bm{y}_{0})$ based on the training data, $\mathcal{D}_{\ell}$, using a single statistical model instead of a vector of model outputs. Estimation of the parameters can now be done cheaply by minimizing the emulated loss function:

[TABLE]

where $\mathbb{E}$ denotes the conditional expectation predicted by the interpolation method, in our case the conditional mean of a Gaussian process. An algorithmic description of the loss emulation method is given in Algorithm 2. For further illustration, an additional example, on the Lotka-Volterra system, can be found in Noè (2019).

The advantage of loss emulation over output emulation is a reduction of the training complexity, as a multi-dimensional vector is replaced by a scalar as the target function. The disadvantage is that, as opposed to output emulation, the emulator can only be trained after the patient has come into the clinic and the training data have become available. This implies that on production of the training data, the emulator has to be trained and the resulting emulated loss function has to be optimized, leading to higher computational costs at the time a clinical decision has to be made. However, these computational costs are still low compared to running the simulator.

Loss emulation is closely related to Bayesian optimization, reviewed e.g. in Shahriari et al. (2016) and Noè (2019), which is a strategy to iteratively include further query points by trading off exploration versus exploitation via some heuristic or information-theoretic criterion. However, every additional query point requires a computationally expensive simulation from the mathematical model, which prevents fast clinical decision making in real time and renders Bayesian optimization infeasible for the purposes of our study.

3.4 Interpolation Methods

We have considered several interpolation methods, based on Gaussian processes (GPs). GPs have been widely used in the context of emulation; see e.g. Kennedy and O’Hagan (2001); Conti et al. (2009); Conti and O’Hagan (2010). For a comprehensive introduction to GPs, the reader is referred to Rasmussen and Williams (2006). Each of the interpolation methods can be used with both of the emulation paradigms described in the previous section, Section 3.3.

3.4.1 Local Gaussian Process

When the sample size $N$ is large, it is not feasible to use exact GP regression on the full dataset, due to the $O(N^{3})$ computational complexity of the $N\times N$ training covariance matrix $\bm{K}$ inversion. A possible approach is to use sparse GPs as in Titsias (2009), which considers a fixed number of $m$ inducing variables $\bm{u}=(u_{1},\dots,u_{m})$, with $m\ll N$, corresponding to inputs $\bm{Z}=[\bm{z}_{1},\dots,\bm{z}_{m}]^{\top}$. The locations of the inducing points and the kernel hyperparameters are chosen with variational inference, i.e. by maximizing a lower bound on the log marginal likelihood, which can be derived by applying Jensen’s inequality. The computational costs of this approach are $O(Nm^{2})$. Initially we tried sparse GPs with 100, 500 and 1000 inducing points but, using the code accompanying the paper by Titsias (2009), the prediction time was between 0.5 and 0.6 seconds for 100 inducing points, around one second for 500, and in the order of a few seconds for 1000 inducing points101010Dual Intel Xeon CPU E5-2699 v3, 2.30GHz, 36 cores and 128GB memory.. This means that minimization of the surrogate-based loss would still be slow as approximately 1 second is required for a single evaluation. The optimization time would exceed two and a half hours for 500 inducing points when using $10,000$ function evaluations. With the cost of variational sparse GP models with larger numbers of inducing points being so large, we can only use about 100 inducing points in order to keep to our goal of real-time in-clinic decision making. However, using such few inducing points was found to lead to around a quarter of the outputs of the biomechanical model being poorly predicted.

With the performance of the variational sparse GPs being poor when the number of inducing points are selected to give a clinically relevant decision time, we instead use a local GP approach based on the $K$-nearest-neighbours instead (Gramacy and Apley, 2015). This method uses the standard GP prediction formulas described in Rasmussen and Williams (2006), but subsetting the training data. Whenever we require a prediction at a given input, we find the training inputs representing the $K$-nearest-neighbours in input-domain, which will form the local set of training inputs, and the corresponding outputs will represent the local training outputs. Note that every time we ask for a prediction at a different input, the training sets need to be re-computed and the GP needs to be trained again. However, because of the small number of neighbours $K\ll 1000$ usually selected, this method is computationally fast and accurate; see Gramacy and Apley (2015) for a discussion.

Gramacy and Apley (2015) further discuss adding a fixed number of distant points in order to help in the estimation of the length scale parameters, but this comes with extra computational costs required by the iterative choice of which point to add to the set of neighbours. Given the time limitations required by our goal (real-time clinical decision support systems) we do not pursue this approach. Furthermore, this is mostly relevant when the interest lies in building predictive models able to make good predictions when the training data are distant from each other. Since we are working on a compact set which is densely covered by the Sobol sequence, this is less relvant. For generic training data $\mathcal{D}=\{(\bm{\theta}_{1},y_{1}),\dots,(\bm{\theta}_{N},y_{N})\}=\{\bm{\Theta},\bm{y}\}$, we give an algorithmic description in Algorithm 3.

In Algorithm 3, the $K\times K$ training covariance matrix is $\bm{K}=[k(\bm{\theta}_{i}^{\prime},\bm{\theta}_{j}^{\prime})]_{i,j=1}^{K}$, the $K\times 1$ vector of covariances between the training points and the test point is $\bm{k}(\bm{\theta}_{*})=(k(\bm{\theta}_{1}^{\prime},\bm{\theta}_{*}),\dots,$ $k(\bm{\theta}_{K}^{\prime},\bm{\theta}_{*}))$ and $\bm{m}=(m(\bm{\theta}_{1}^{\prime}),\dots,m(\bm{\theta}_{K}^{\prime}))$ is the $K\times 1$ prior mean vector. We consider a constant mean function $m(\bm{\theta})=c$. For the kernel $k(.,.)$ we choose the Automatic Relevance Determination Squared Exponential kernel (see e.g. Rasmussen and Williams (2006)), as widely used in the emulation of computer codes literature; see e.g. Fang et al. (2006); Santner et al. (2003). The kernel hyperparameters are the output scale (determining the function variance) and the input length scales, one length scale for each dimension. These hyperparameters are estimated by maximizing the log marginal likelihood using the Quasi-Newton method. The standard deviation of the additive Gaussian noise, $\sigma$, is initialized at a small value, $\sigma=10^{-2}$, to reflect the fact that the mathematical model of the LV is deterministic111111Even for deterministic models, a small non-zero value for $\sigma$ is usually assumed, to avoid numerical instabilities of the covariance matrix inversion.

The CPU time required to get a prediction from the local Gaussian process is approximately 0.18 seconds121212Dual Intel Xeon CPU E5-2699 v3, 2.30GHz, 36 cores and 128GB memory. using the $K=100$ nearest neighbours of a given point. The number of neighbours $K$ needs to be selected on the basis of the computational time allowed to reach a decision in a viable time frame, but keeping in mind that $K$ also controls the accuracy of the emulation. In our experiments we found that $K=100$ was sufficiently fast for the method to be applicable in the clinic while leading to accurate predictions at the test inputs, as discussed below in the Results section.

For this method, the surrogate-based loss and the emulated loss were optimized using the Global Search algorithm by Ugray et al. (2007), implemented in MATLAB’s Global Optimization toolbox.131313Available from https://uk.mathworks.com/products/global-optimization.html. We use the default choice of 2000 trial points and 400 stage one points. Consider running a local solver from a given starting point $\bm{\theta}_{0}$, ending up at the point of local minimum $\hat{\bm{\theta}}$. The basin of attraction corresponding to that minimum is defined as the sphere centred at $\hat{\bm{\theta}}$ and having radius equal to $\|\bm{\theta}_{0}-\hat{\bm{\theta}}\|$. All starting points falling inside the sphere are assumed to lead to the same local minimum $\hat{\bm{\theta}}$, hence no local solver is run and they are discarded. In simple words, stage one of the Global Search algorithm scatters initial points in the domain and scores them from best to worst by evaluating the function value and constraints. Then an interior-point local solver (Byrd et al., 2000) is run from each trial point, starting from the one that was scored best (lowest function value), and excluding points that fall into the basins of attraction of previously found minima. When all the stage one points have been analyzed, stage two generates more random points and the same procedure is run a second time.

3.4.2 Low-Rank Gaussian Processes

Along with local GPs based on the $K$-nearest-neighbours, described in Section 3.4.1, we report results for another type of statistical approximation: low-rank GPs, as described in Section 5.8.2 of Wood (2017), whose main ideas are summarized here for generic training data $\mathcal{D}=\{(\bm{\theta}_{1},y_{1}),\dots,(\bm{\theta}_{n},y_{n})\}=\{\bm{\Theta},\bm{y}\}$.

Let $\bm{C}=\bm{K}+\sigma^{2}\bm{I}$ be the $n\times n$ covariance matrix of $\bm{y}$ and consider its eigen-decomposition $\bm{C}=\bm{U}\bm{D}\bm{U}^{\top}$ with eigenvalues $|D_{i,i}|\geq|D_{i+1,i+1}|$. Denote by $\bm{U}_{k}$ the submatrix consisting of the first $k$ eigenvectors of $\bm{U}$, corresponding to the top $k$ eigenvalues in $\bm{D}$. Similarly, $\bm{D}_{k}$ is the diagonal matrix containing all eigenvalues greater than or equal to $D_{k,k}$. Wood (2017) considers replacing $\bm{C}$ with the rank $k$ approximation $\bm{U}_{k}\bm{D}_{k}\bm{U}_{k}^{\top}$ obtained from the eigen-decomposition. Now, the main issue is how to find $\bm{U}_{k}$ and $\bm{D}_{k}$ efficiently enough. A full eigen-decomposition of $\bm{C}$ requires $O(N^{3})$ operations, which somewhat limits the applicability of the rank-reduction approach. A solution is to use the Lanczos iteration method to find $\bm{U}_{k}$ and $\bm{D}_{k}$ at the substantially lower cost of $O(N^{2}k)$ operations, see Section B.11 in Wood (2017). Briefly, the algorithm is an adaptation of power methods to obtain the truncated rank $k$ eigen-decomposition of an $N\times N$ symmetric matrix in $O(N^{2}k)$ operations. However, for large $N$, even $O(N^{2}k)$ becomes prohibitive. In this scenario the training data are randomly subsampled by keeping $n_{r}$ inputs and an eigen-decomposition is obtained for this random selection with $O(n_{r}^{2}k)$ computational cost.

We used the implementation found in the R package mgcv by Wood (2017), with the following settings: $n_{r}=2000$ (the package default), $k=2000$ for output emulation, while $k=1000$ for loss emulation. The kernel used was an isotropic Matérn 3/2 kernel, with lengthscale set to the default of Kammann and Wand (2003): $\lambda=\max_{ij}\|\bm{\theta}_{i}-\bm{\theta}_{j}\|$. The remaining model hyperparameters are estimated by maximizing the log marginal likelihood. The final model used an interaction term between each of the 4 model parameters, as well as a second interactive term between the inverses of the model parameters:

[TABLE]

where $\bm{\tau}=1/\bm{\theta}$, $f(\bm{\theta})\sim\mathsf{GP}_{\mathsf{LR}}(m(\bm{\theta}),K(\bm{\theta},\bm{\theta}^{\prime}))$, $f(\bm{\tau})\sim\mathsf{GP}_{\mathsf{LR}}(m(\bm{\tau}),K(\bm{\tau},\bm{\tau}^{\prime}))$ and $\mathsf{GP}_{\mathsf{LR}}(\cdot)$ denotes a low rank GP. The model specification with the two interaction terms was found to reduce the variation in the predictive accuracy as the volume increases and the strains decrease. This can be seen in the predictions of the test and training data in Figure 2 and 3 of the online supplementary materials.

Minimization of the surrogate-based loss $\ell(\cdot\mid\hat{\bm{\mathsf{m}}},\bm{y}_{0})$ and the emulated loss $\hat{\ell}(\cdot\mid\bm{\mathsf{m}},\bm{y}_{0})$ is performed by the Conjugate Gradient method implemented in the R function optim (Nash, 1990), with maximum number of iterations set to 100. To avoid being trapped in local minima, 50 different starting points from a Sobol sequence were used. The best minimum found was kept as the estimate, discarding the remaining 49 optima.

3.4.3 Multivariate-output Gaussian Processes

The previous two subsections have focussed on single-output GPs, while potentially correcting for the correlation structure of the outputs via a modified objective function, using the Mahalanobis distance defined in (6). One can model the correlation structure between the outputs directly via

[TABLE]

where $K(\bm{\theta}_{i},\bm{\theta}_{j})$ is the covariance between $y_{k}(\bm{\theta}_{i})$ and $y_{k}(\bm{\theta}_{j})$ for any output $k$, and ${\bf A}$ is a matrix of the covariances between the outputs, i.e. the circumferential strains and the LV volume. Various approaches have been proposed in the literature. The approach taken in Conti and O’Hagan (2010) and Conti et al. (2009) is to place a non-informative prior on ${\bf A}$ and integrate ${\bf A}$ out in the likelihood. This leads to a closed-form solution in terms of a matrix-normal distribution; see Conti and O’Hagan (2010) and Conti et al. (2009) for explicit expressions. However, we found that in combination with Algorithm 3 – to deal with the $O(N^{3})$ computational complexity – the computational costs of running the emulator were in the order of hours, rather than minutes, which renders this approach not viable for clinical decision support in real time.

An alternative approach is to explicitly model the correlation structure of the outputs via

[TABLE]

taking into account covariates ${\bf u}_{k}$ and ${\bf u}_{l}$ associated with the $k$th and $l$th outputs, $y_{k}$ and $y_{l}$, respectively. Roberts et al. (2013) pursue this approach in the context of time series analysis, where ${\bf u}_{k}$ and ${\bf u}_{l}$ are scalar variables indicating different time points. In our application, ${\bf u}_{k}$ and ${\bf u}_{l}$ are vectors indicating the locations on the surface of the left ventricle associated with the circumferential strains. Due to the highly non-Euclidean geometry of this space, the choice of kernel is not obvious. A naive approach that we tried is to project the locations onto a linear space defined by the first principal component (Huang, 2016). The results were not encouraging, due to the information loss incurred by the map. Future work could try projections onto nonlinear maps, like Hilbert curves (Hilbert, 1891; Hamilton and Rau-Chaplin, 2007), generative topographic maps (Bishop et al., 1998), or self-organising maps (Kohonen, 1982).

A further alternative is the method of Alvarez and Lawrence (2009, 2011), who have proposed sparse convolved GPs for multi-output regression. Their method assumes that there is an underlying process which governs all of the outcomes of the model and treats it as a latent process. Modelling this latent process as a GP leads to a GP prior over the outputs, inducing cross covariance between the outputs and effectively introducing correlations between them. We can use the interpolation method of Alvarez and Lawrence (2009, 2011) within either of the emulation frameworks introduced in Section 3.3. There are however problems with doing this: training a convolved GP with $N$ training points requires the inversion of a $DN\times DN$ matrix (where $D=25$ is the number of outputs) which is currently infeasible with all of the training data ($N=10,000$), even when choosing the number of inducing points using the method proposed in Alvarez and Lawrence (2009, 2011). Instead we can choose a strategy similar to that proposed in Section 3.4.1. This again, however, proves to be computationally expensive as fitting a single local emulator requires more than 15 minutes141414Intel Xeon CPU E5-606,2.13GHz, without consideration of the computational costs of the subsequent optmization of the LV model parameters. When this is included within either of the emulation methods (Algorithms 1 and 2), the time becomes too large for a clinical decision support system, as it is infeasible to make a prediction within a clinically relevant time frame.

Since the focus of our study is to develop an emulation framework for a clinical decision support system that can work in real time, we have restricted our analysis to the univariate methods described in Sections 3.4.1 and 3.4.2.

4 Data and Simulations

For training the emulator, we used 10,000 parameter vectors generated from a Sobol sequence (Sobol, 1967) in a compact 4-dimensional parameter space, with $\theta_{1},\ldots,\theta_{4}\in[0.1,5]^{4}$, where the parameter bounds reflect prior knowledge available from Gao et al. (2015). The 4-dimensional parameter vectors are then transformed to the original 8-dimensional parameter space using the transformation (2). The 8-dimensional parameter vectors are then inserted into the HO strain energy function (1). Following the finite element discretisation method described in Wang et al. (2013), the soft-tissue mechanical equations are numerically solved to produce a 25-dimensional output vector associated with each parameter vector; these are 24 circumferential strains and the LV volume at end of diastole. The Sobol sequence is extended to generated an independent test set of additional 100 parameter vectors, for which the same procedure is followed to associate them with output vectors of circumferential strains and LV volume. As a real data set, we used 24 circumferential strains and the LV volume at end of diastole obtained from the cardiac MRI images of a healthy volunteer, following the procedure described in Gao et al. (2015).

5 Results

To summarise, we have introduced two emulation frameworks which can be used to infer the parameters of the LV biomechanical model; see Sections 3.3.1 and 3.3.2. We have applied these methods with two different loss functions, the Mahalanobis loss function and the Euclidean loss function, and two different interpolation methods, low rank GPs and local GPs; see Sections 3.4.1 and 3.4.2. Testing each combination of these methods means that there is a total of 8 different alternative procedures.

We have applied and assessed the proposed methods in a two-pronged approach. Firstly, in Sections 5.1, 5.2, 5.3 and 5.4, we have tested the 8 different combinations of methods on synthetic data, where the true parameter values of the underlying biomechanical model are known; see the previous section for details on how the training and test data were generated. We compare the methods using the Mean Square Error (MSE).151515Note that the likelihood is computationally expensive and intractable. Hence, we do not compare the methods using the log posterior of the parameters, as this would involve approximations, using e.g. variational methods or expectation propagation, and this would have to be repeated 100 times (the number of test data) at high computational costs. The distribution of 100 MSEs is given in Figure 2 and summarised with the median and the (1st, 3rd) quartiles in Table 1, representing 3 of Tukey’s five number summary.161616We do not present plus or minus the interquartile range as this can lead to the wrong impression of a negative MSE.

Finally, we have applied the method with the best performance in Section 5.5 to clinical data generated from a healthy volunteer’s cardiac MRI scan, where we can compare our performance against the gold standard results of Gao et al. (2015).

5.1 Comparison of Interpolation Methods

Looking at the two interpolation methods, the local GP method (boxplots 5-8 in Figure 2) outperforms the low rank GP method (boxplots 1-4 in Figure 2). The reason for the difference in performance between the two methods is the size of the noise variance that is estimated. With the low rank GP method, a larger noise variance is estimated as the interpolation must fit to the entire dataset. The larger variance of the errors is in mismatch with the deterministic nature of process that we are modelling. Instead of estimating the variance with maximum likelihood, one could consider a Bayesian approach that discourages larger values of the variance with a restrictive prior. However, besides the confounding effect of switching statistical inference paradigms (from maximum likelihood to Bayesian inference), the available code for this is not available in the mgcv package in R.

Conversely, with the local GP method, a much smaller error variance is estimated, which more closely matches the deterministic data generation method. This is a result of there only being a small number of points that the interpolant must fit. These points are local, giving more detail of the local surface than the low rank GP method, which uses a selected number of points from across the whole dataset.

5.2 Comparison of Emulation Frameworks

Out of the two emulation frameworks, the output emulation method (boxplots 1, 2, 5 and 6 in Figure 2) gives the most accurate parameter estimates, outperforming the loss emulation method (boxplots 3, 4, 7 and 8 in Figure 2) for all interpolation methods and loss functions. The output emulation method provides accurate estimates for all the different combinations of interpolation methods and loss functions, while the loss emulation method provides poor estimates in some cases. The improved parameter estimation of the output emulation method is a result of using multiple separate emulators. These multiple emulators better model the complex non-linear relationships between the parameters and the outputs than it is possible with the single emulator used with the loss emulation method. In the loss emulation method, the differences between the patient data and the simulations are summarised in one loss function, reducing the performance of the method.

5.3 Comparison of Loss Functions

In terms of the accuracy of the parameter inference, the Euclidean loss and Mahalanobis loss perform differently in different emulation methods. Firstly, for the loss emulation method the Mahalanobis loss function (boxplots 3 and 7 in Figure 2) clearly outperforms the Euclidean loss function (boxplots 4 and 8 in Figure 2) in all cases. The reason for the difference is that the loss function summarises how similar the patient data is to the simulations and this is done more realistically by the Mahalanobis loss function in this case. This is because there are spacial correlations between the outputs due to measuring the circumferential strains at different neighbouring locations on the left ventricle. The Mahalanobis loss function accounts for this through including a correlation estimate, whereas the Euclidean loss function does not.

In comparison to the loss emulation method, for the output emulation method it is less clear which loss function gives the best results. The Mahalanobis loss function is marginally better for the low rank GP method (boxplot 1 is better than boxplot 2 in Figure 2), while the Euclidean loss function gives the best performance for the local GP method (boxplot 6 is better than boxplot 5 in Figure 2). The reason for the Euclidean loss function performing best for the local GP method is because of potential inaccuracies in the correlation matrix used for the Mahalanobis loss function. Firstly, the correlation matrix is a global measure based on the whole dataset and may not accurately represent the true correlations between the local points due to limited numerical precision171717Using a local correlation matrix was also tested, but limited accuracy and numerical stability of the correlation matrix due to using only a small number of local points meant that the performance did not improve over the global correlation matrix.. Secondly, this is aggravated by lack of numerical stability when inverting the covariance matrix. Thirdly, the loss function minimised in the output emulation method is based on the errors between the emulators and patient data, whereas the correlation matrix has been calculated based on only the patient data.

5.4 Overall Best Method in Simulation Study

In conclusion, the results of our simulation study show the following. (1) The local GP method outperforms the low rank GP method and is the better of the two interpolation methods. (2) The best emulation method is the output emulation method and this outperforms the loss emulation method in all the different combinations of interpolation method and loss function tested. (3) The Mahalanobis loss function gives the best performance for the loss emulation method. (4) For the output emulation method, the Mahalanobis method is marginally better for the low rank GP method, but for the local GP method the Euclidean loss function gives the best parameter estimates. (5) Overall, the simulation study results show that the best performing combination of methods is the output emulation method, using the local GP as the interpolation method and the Euclidean loss function (boxplot 6 in Figure 2). This combination of methods will be used on the cardiac MRI data of the healthy volunteer in Section 5.5.

5.5 Application to Cardiac MRI Data

Figure 2 and Table 1 show that the method which gives the most accurate parameter prediction is the emulation of the outputs method combined with the local GP interpolation and the Euclidean loss function. We have applied this strategy to estimate the material parameters for the heart model of a healthy volunteer described in Section 2, using the set of 24 circumferential strains and the LV cavity volume extracted from cardiac MRI images, as described in Section 4. The true model parameters are not known in this case, so as opposed to the simulation study we do not have a proper ‘gold standard’ for evaluation. We therefore use the following alternative procedure. We first estimate the constitutive parameters with the method of Gao et al. (2015, 2017a), that is, with the method using the computationally expensive simulator. From these parameters, we calculate the stretch-stress relationships along the directions of the sheets and the myocytes, following the procedure described in Holzapfel and Ogden (2009). We use these graphs as a surrogate ‘gold standard’, which we compare with the corresponding graphs obtained from the parameters obtained with our emulation approach.

Figure 3 shows, as dashed lines, the estimate of the stretch-stress relationship for the healthy volunteer using the ‘gold standard’ method of Gao et al. (2015, 2017a). For comparison, the solid blue lines show the estimates of the stress-stretch relationship obtained from the best emulation method identified in the previous sections, Sections 5.1–5.4, the emulation of the outputs method combined with the local GP interpolation method and the Euclidean loss function.

For uncertainty quantification, we numerically estimated the Hessian at the minimum surrogate loss (5). Its inverse represents an approximate lower bound on the variance-covariance matrix in parameter space181818The Hessian is the empirical Fisher information matrix. The lower bound would be exact (Cramer-Rao lower bound) if we could take an expectation with respect to the data distribution. Recall that saying that matrix A is a lower bound on matrix B means that B - A is positive semi-definite.. The uncertainty in the estimate can then be obtained by sampling from a multivariate normal distribution, with the covariance set to the inverse of the Hessian, $\mathrm{MVN}(\hat{\bm{\theta}},\textbf{H}(\hat{\bm{\theta}})^{-1})$, and calculating the corresponding confidence intervals.

The results in Figure 3 show that the emulation method accurately estimates the stretch-stress relationship in the myocyte direction. The agreement between the ‘gold standard’ and the prediction with our emulation method is nearly perfect, with a deviation that is less than the predicted single-standard deviation width. For the stretch-stress relationship in the sheet direction, the agreement is also very good, although the deviation exceeds the predicted standard deviation in this case. A possible explanation is that parameter sensitivity in the sheet directions is very low when only using regional circumferential strains and the LV cavity volume to formulate the objective function, as reported in Gao et al. (2015), thus the uncertainty of estimating the stiffness in the sheet direction will be higher than that in the myocyte direction. It is expected that higher accuracy will be achieved when radial (transmural) strains are included when inferring the parameters. While the differences between the stretch-stress curves obtained with the simulator and our emulator are minor, there is a substantial difference in the computational costs. For the simulator, that is, the original procedure described in Gao et al. (2015, 2017a), the computational costs are in the order of over a week. The estimation procedure with the proposed emulator, on the other hand, could be carried out in less than 15 minutes191919Dual Intel Xeon CPU E5-2699 v3, 2.30GHz, 36 cores and 128GB memory., giving us a reduction of the computational complexity by about three orders of magnitude.

Hence, while the former procedure is only of interest in a pure research context, the latter procedure gives us estimation times that are acceptable in a clinical decision context. This is an important first step towards bringing mathematical modelling into the clinic and making a real impact in health care.

6 Discussion

We have developed an emulation framework that can be used to infer the material properties of the LV of a healthy patient in a clinically viable time frame. We have focused on developing an emulation framework that can be used in future more generalised work and have therefore tested 2 emulation methods, 2 interpolation method and 2 loss functions; see Section 3. Each combination of these methods has then been evaluated in a simulation study in order to determine the best method. The best method was found to be the output emulation method, using the local GP as the interpolation method and the Euclidean loss function; see Table 1.

We have then applied the proposed emulation method to cardiac MRI data and demonstrated that it is able to accurately estimate the stretch-stress relationship along the myocyte and sheet directions of the LV from a healthy volunteer. Our method provides a notable improvement in computational time with a speed-up of approximately 3 orders of magnitude. In particular, while conventional parameter estimation based on numerical simulations from the mathematical LV model, following e.g. the approach of Gao et al. (2015), leads to computational costs in the order of weeks, the proposed emulation method reduces the computational complexity to the order of a quarter hour, while effectively maintaining the same level of accuracy. This is an important step towards a clinical decision support system that can assist a clinical practitioner in real time.

A limitation of the current approach is the fact that the LV geometry is fixed. This LV geometry varies from patient to patient, and these variations need to be taken into consideration for applications to wider patient populations. We discuss how to potentially address this challenge in the next section.

7 Future Work

The next step for this work is to design a method that is capable of fast parameter inference for multiple patients on whom we have not directly trained the emulator. For each new patient we would need to replace the single geometry used here as an input, with the new patient’s data on arrival at the clinic. With no time limits on the inference, we could simply replicate this study with a different input geometry. However, in order to treat patients in a clinically viable time frame we must be able to train the emulator for the unobserved patient before they enter the clinic. We can do this by using simulations from multiple different LV geometries as our training data. Low-dimensional representations of each geometry can then be included as variables in the interpolation method of the emulator and we can learn how these changes affect the output of the biomechanical model. When new patient data then arrives, these low dimensional representations can be calculated and included in the loss function, which must be minimised in the emulation method.

A straightforward approach for achieving this low-dimensional representation is principle component analysis (PCA), illustrated in Figure 4(a), where the high-dimensional LV geometries are mapped onto a low-dimensional space that captures the maximum variation in the population. A variation along the PCA directions can be mapped back into the high-dimensional LV geometry space to illustrate typical organ deformations, as illustrated in Figure 4(b). However, while fast and easy to implement, the limitation of PCA is its restriction to linear subspaces. If the LV geometries extracted from the patient population are grouped along a non-linear submanifold in the high-dimensional LV geometry space, as illustrated in Figure 4(c), PCA is suboptimal. A variety of non-linear extensions of and alternatives to PCA have been proposed in the machine learning and computational statistics literature. The most straightforward extension is kernel PCA (Scholkopf et al., 1998), which conceptually maps the data non-linearly into a high-dimensional vector space and makes use of Mercer’s theorem, whereby the scalar product in this high-dimensional space is equivalent to a kernel in the original data space and therefore never has to be computed explicitly. Alternative non-linear dimension reduction methods to be explored are generative topographic maps (Bishop et al., 1998), self-organising maps (Kohonen, 1982), and variational auto-encoding neural networks (Kingma and Welling, 2014).

Acknowledgement

This work was funded by the British Heart Foundation, grant numbers PG/14/64/31043 and RE/18/6134217, and by the UK Engineering and Physical Sciences Research Council (EPSRC), grant number EP/N014642/1, as part of the SofTMech project. Benn Macdonald is supported by The Biometrika Trust, Fellowship number B0003. Umberto Noè was supported by a Biometrika Scholarship. Alan Lazarus is partially funded by a grant from GlaxoSmithKline plc.