Efficient Ventricular Parameter Estimation Using AI-Surrogate Models

The onset and progression of pathological heart conditions, such as cardiomyopathy or heart failure, affect its mechanical behaviour due to the remodelling of the myocardial tissues to preserve its functional response. Identification of the constitutive properties of heart tissues could provide useful biomarkers to diagnose and assess the progression of disease. We have previously demonstrated the utility of efficient AI-surrogate models to simulate passive cardiac mechanics. Here, we propose the use of this surrogate model for the identification of myocardial mechanical properties and intra-ventricular pressure by solving an inverse problem with two novel AI-based approaches. Our analysis concluded that: (i) both approaches were robust toward Gaussian noise when the ventricle data for multiple loading conditions were combined; and (ii) estimates of one and two parameters could be obtained in less than 9 and 18 s, respectively. The proposed technique yields a viable option for the translation of cardiac mechanics simulations and biophysical parameter identification methods into the clinic to improve the diagnosis and treatment of heart pathologies. In addition, the proposed estimation techniques are general and can be straightforwardly translated to other applications involving different anatomical structures.

1. Introduction

Cardiovascular disease is the largest cause of death worldwide. Effective diagnosis and treatment are hampered by a lack of knowledge of the pathophysiological mechanisms underlying the development of the disease. Biomechanical factors, such as stiffness and stress, are known to have important influences on heart function, but are difficult to quantify. Patient-specific computer models of heart biomechanics allow intrinsic constitutive muscle properties, including stiffness, contraction, relaxation, stress and work, to be assessed using medical data from cardiac catheterisation and imaging (Wang et al., 2018). Such cardiac tissue indices provide a new dimension of diagnosis that can help to elucidate the mechanisms of heart disease, thus enabling more specific targeting of treatment and ultimately better outcomes for patients.

Simulating biomechanics via computational models of the heart is challenging due to the stress-strain non-linearities intrinsic of cardiac tissue undergoing large deformation. The finite element method (FEM) is the most commonly used approach to solve the finite elasticity governing equations to enable accurate predictions of these large deformations (Nash and Hunter, 2000). This requires specifying constitutive relations to describe the stress-strain behaviour of the myocardium. Calibration of patient-specific parameters of these constitutive relations typically involves minimising an objective function that quantifies kinematic discrepancies between the FE model's predictions and measurements from medical images. These kinematic measurements typically involve quantifying shape change or the displacement of the myocardium over the cardiac cycle. For example, segmentations of the endocardial and epicardial contours of the heart from Cine MRI (Chen et al., 2020) have been used as shape-based kinematic measurements for FE model calibration (Wang et al., 2018). Kinematic measurements that involve quantify the displacement of the heart wall have also been used for FE model calibration (Hadjicharalambous et al., 2017; Zhang et al., 2021). These approaches track the displacements of the tissue using techniques such as optical flow (Queiros et al., 2017) or cardiac magnetic resonance tagging (Aletras et al., 1999; Zhong et al., 2010; Ibrahim, 2011; Shi et al., 2012; Amzulescu et al., 2019).

A number of non-linear optimisation methods have been used to calibrate constitutive parameters of cardiac mechanics models. This includes the application of gradient-based methods (Gao et al., 2015; Wang et al., 2018) and gradient-free methods (Rumindo et al., 2020; Zhang et al., 2021). Each of these methods requires multiple evaluations of the objective function for each update of the parameters, and each evaluation involves a costly FE simulation of the mechanics model. The evaluation of the gradient is often obtained from finite difference approximations, which can involve a significant computational cost when increasing number of parameters need to be identified. Calibration procedures using gradient-based or gradient free methods are therefore computationally expensive, often taking many hours or days to complete [e.g., (Gao et al., 2015) reported calibration times of 63 h, and (Zhang et al., 2021) reported calibration times of 15 h]. This presents a significant barrier to clinical translation of cardiac mechanics models.

Surrogate models have been developed to improve simulation efficiency, replacing computationally expensive FE simulations with inexpensive surrogate simulations (Dabiri et al., 2019; Maso Talou et al., 2020; Cai et al., 2021). These surrogate models have also been used to accelerate calibration procedures (Di Achille et al., 2018; Davies et al., 2019; Noe et al., 2019; Longobardi et al., 2020; Cai et al., 2021). For example, Cai et al. (2021) presented an approach that used surrogate models with a trust-region-reflective gradient-based optimiser for estimating personalised constitutive parameter of the left ventricle. While these advancements are promising, existing calibration approaches do not exploit the full potential of neural networks-based surrogate models, which can directly provide the analytic gradients of objective functions of interest during the calibration procedure via automatic differentiation (Raissi et al., 2019). By using these analytic gradients, we only need to evaluate the model once for each iteration of gradient-based optimisation procedures. This can lead to an efficient personalisation of these cardiac models, and enable their future application to real-time continuous monitoring of cardiac function.

In this work, we propose a novel approach that uses AI-surrogate model and automatic differentiation to efficiently identify constitutive parameters or loading conditions of a biomechanical model given kinematic measurements from medical images. This builds upon our recent work in developing deep learning approaches that substantially reduce the computational cost of simulating cardiac biomechanics, by training an AI-surrogate model that accurately reproduces mechanics predictions with a fraction of the computational cost of numerical methods that solve the governing partial differential equations (Maso Talou et al., 2020). We demonstrate this approach for estimating passive stiffness or pressure of the left ventricle (LV). This approach involves performing an optimisation of the AI-surrogate inputs to best match its kinematic response against a set of given observations from medical images. We present two different strategies: (i) a full-field tracking approach, which requires the displacement field between two medical images; and (ii) a contour matching approach, which requires only the geometry of the ventricular surfaces in between two medical images. The latter approach relaxes the requirement for determining the displacement field of the tissue, and only requires contours describing the surface of the ventricle from medical images to determine the best matching kinematic response of the AI-surrogate.

The manuscript is structured as follows. In section 2, we introduce the FE mechanical model for the left ventricle of the heart, and its AI-surrogate. Then, we present the proposed parameter identification strategies followed by the setup of the optimisation scheme. In section 3, we study the performance of both strategies for the identification of constitutive parameters and haemodynamic loading conditions. Finally, we discuss the contributions and limitations of this work in section 4 and outline our final remarks in section 5.

2. Methods

2.1. Mechanical Model

Kinematics of the LV are simulated using a patient-specific FE model (Wang et al., 2018). This involves solving the finite elasticity equilibrium equations during the diastolic phase of the cardiac cycle under an endocardial pressure boundary condition to simulate passive filling. Patient-specific geometrical models of the LV are constructed at the diastasis frame of the cardiac cycle for a range of individuals, which are assumed to be in a load-free configuration. Cubic Lagrange basis functions are used for constructing the FE mesh of the geometry. A typical mammalian description of the myocyte orientation through the LV wall (Nielsen et al., 1991) is incorporated into the geometry through a material fibre field. The LV myocardium is modelled as an ideally-incompressible transversely isotropic material by means of a Fung-type exponential constitutive model (Guccione et al., 1991) with the following strain energy density function

$\begin{array}{ll}\begin{array}{l}\text{Ψ}=\frac{c_1}{2}(e^Q-1)\\ Q=c_2{E}_{ff}^2+c_3(E_{cc}^2+E_{rr}^2+E_{cr}^2)+2c_4(E_{fc}^2+E_{fr}^2)\end{array}& (1)\end{array}$

where the c₁ parameter scales the overall stiffness of the myocardium, and c₂, c₃, and c₄ control the material non-linearity and anisotropy in the fibre (f), cross-fibre (c), and radial (r) directions, respectively. Incompressibility of the myocardium was enforced through a mixed formulation that uses linear Lagrange basis functions for describing the hydrostatic pressure (Nash and Hunter, 2000). Homogeneous Dirichlet boundary conditions were applied on nodes of the FE mesh at the epicardial perimeter of the basal surface of the model. All FE simulations were performed using the OpenCMISS-Iron open-source computational modelling software package (Bradley et al., 2011).

2.2. Surrogate Network Model

From the FE mechanical model defined in section 2.1, we derive a surrogate model as described in Maso Talou et al. (2020). The AI-surrogate model predicts the displacement of a material point x = (x_d, y_d, z_d) for a given intra-ventricular pressure p, domain description g = (g₁, g₂) (PCA weights as explained in Maso Talou et al., 2020) and the constitutive parameters (c₁, c₂). Particularly, we fixed parameters c₃ = 3.67 and c₄ = 25.77 as their identifiability from macro-scale observations is low.

To train the AI-surrogate, we minimised the squared displacement error with the FE predictions given by training the following loss function

$\begin{array}{ll}l({B})={\displaystyle \sum _{{\text{u}}_d\in {{B}}_d}}\parallel {\text{u}}_d-{\tilde{\text{u}}}_d{\parallel }_{L_2}^2+\alpha {\displaystyle \sum _{{\text{u}}_b\in {{B}}_b}}\parallel {\text{u}}_b-{\tilde{\text{u}}}_b{\parallel }_{L_2}^2& (2)\end{array}$

where ${{B}}_b$ and ${{B}}_d$ define sets of points (training batches) on the basal/endocardial boundary (where Dirichlet or Neumann boundary conditions were applied) and inside the domain, respectively, u_i and ${\tilde{\text{u}}}_i$ are the displacements predicted with the neural network and the FE model, respectively, α = 4.5 is the penalty factor to impose the boundary conditions, and ${B}={{B}}_b+{{B}}_d$ is a given training batch. For further details about the training of the AI-surrogate, refer to (Maso Talou et al., 2020).

2.3. Parameter Identification Strategies

In this work, we propose two strategies to identify inputs of the surrogate model from a given set of observations of ventricular kinematics from medical images. As presented in our previous work (Maso Talou et al., 2020), we can encode boundary conditions, applied tractions, domain geometry and constitutive parameters as inputs of these networks, thus the proposed techniques can be used to characterise any of these inputs.

The estimation of parameters is studied using two approaches, which differ based on the image data available to quantify the kinematics of the heart. One approach considers a displacement field description of the cardiac wall throughout the cardiac cycle, which can be derived from medical images using motion tracking methods (Wang and Amini, 2011; Shi et al., 2012), or post-processing functions for MR sequences that provide effective tracking of material points (e.g., CMR tagging). For the second approach, the input data describe the kinematics of the cardiac wall surfaces during the heart cycle. For each approach, two different time-points, corresponding to the initial and final positions of the cardiac walls, were combined with the measured intra-ventricular pressures that obtained from cardiac catheterisation for the analyses. Quantification of the cardiac wall surfaces can be obtained through segmentation of clinical images, such as 3D echocardiography or cardiovascular magnetic resonance imaging (CMR) (Chen et al., 2020).

2.3.1. Full-Field Tracking Approach

The full-field tracking approach involves finding the parameters Θ (inputs of the surrogate model) that minimise the error between the surrogate model's predictions of tissue motion and displacement data derived from medical images, i.e.,

$\begin{array}{l}\Theta ={\displaystyle \underset{\widehat{\Theta }}{arg min}}{L}(\widehat{\Theta },\widehat{\text{u}})\\ ={\displaystyle \underset{\widehat{\Theta }}{arg min}}\parallel \text{u}(\widehat{\Theta })-\widehat{\text{u}}{\parallel }_{L_2}& (3)\end{array}$

where $u(\widehat{\Theta })$ is the displacement field predicted by the surrogate model for parameters $\widehat{\Theta }$, and $\widehat{\text{u}}$ is the displacement field from medical images.

2.3.2. Contour Matching Approach

Let us define a contour s_t at the time-point t composed of P points as ${\text{s}}_t=({\text{x}}_t^0,\dots ,{\text{x}}_t^P)$. The endocardial and epicardial contours are extracted from medical images at the reference (diastasis) and pressure loaded (end-diastolic) time-points. We denote these contours as s^endo and s^epi for endocardial and epicardial contours, respectively, and in the following we use subscripts i or f to indicate the absence or presence, respectively, of pressure load.

The contour matching approach involves finding the parameters Θ (inputs of the surrogate model) that minimise the error between the initial contours ($s_i^{·}$) displaced by surrogate model predictions, and the end-diastolic contours identified from the medical images, namely

$\begin{array}{l}\Theta ={\displaystyle \underset{\widehat{\Theta }}{arg min}}{L}(\widehat{\Theta },{\text{s}}_i^{\text{endo}},{\text{s}}_f^{\text{endo}},{\text{s}}_i^{\text{epi}},{\text{s}}_f^{\text{epi}})\\ ={\displaystyle \underset{\widehat{\Theta }}{arg min}}\left(d({\text{s}}_i^{\text{endo}}+\text{u}(\widehat{\Theta }),{\text{s}}_f^{\text{endo}})+d({\text{s}}_i^{\text{epi}}+\text{u}(\widehat{\Theta }),{\text{s}}_f^{\text{epi}})\right)& (4)\end{array}$

where the surface-field addition yields a predicted surface ${\text{s}}_t+\text{u}=({\text{x}}_t^0+\text{u}{|}_{{\text{x}}_t^0},\dots ,{\text{x}}_t^p+\text{u}{|}_{{\text{x}}_t^p})$ with u|_x being the displacement of the cardiac wall at the spatial position x, and

$\begin{array}{ll}d({\text{s}}_1,{\text{s}}_2)=\frac{1}{|{\text{s}}_1|}{\displaystyle \sum _{\text{x}\in {\text{s}}_1}}{\displaystyle \underset{\text{y}\in {\text{s}}_2}{min}}\parallel \text{x}-\text{y}\parallel & (5)\end{array}$

where ∥·∥ is the Euclidean norm and |s| is the cardinality of s.

2.4. Optimisation Scheme

A surrogate model of the ventricle is obtained as described in Maso Talou et al. (2020) by minimising a displacement error metric (see Equation 2) between the model predictions and a finite element model. During the training of the AI-surrogate, both the training and testing errors decreased monotonically across epochs, reaching a plateau at 4.43 × 10⁻⁴ mm² and 1.30 × 10⁻³ mm², respectively. The resulting surrogate model presents an absolute displacement error of 0.0499 ± 0.0374 mm in ranges of intraventricular pressure (p ∈ [0.15, 1.5] kPa) and myocardial elasticity (c₁ ∈ [2, 5] kPa and c₂ ∈ [4, 40]) within the physiological ranges.

As the parameters to be estimated are inputs of surrogate model's network, we freeze the weights of the network layers and optimise the inputs by solving Equations (3) or (4). By using Tensorflow v2.1, we implemented the network and the objective function in a Tensorflow graph. Via this execution graph and automatic differentiation, we obtain the analytic derivatives of Equations (3) or (4) with respect to the parameters to be estimated. Finally, we employ the derivatives in an ADAM method with exponential decay of the learning rate given by

$\begin{array}{ll}{\tau }_i={\tau }_{0}0.985^{i/10}& (6)\end{array}$

where τ_i is the learning at the ith epoch, and τ₀ is the initial learning rate. We empirically choose τ₀ = 1 for both full-field tracking and contour matching approaches. By using learning decay, the optimisation proceeds more rapidly in the neighbourhood of the target value, and the overshooting oscillations are damped, as the learning rate diminishes, leading to convergence. The optimisation process stops when estimates of all of the parameters show a relative change lower than 10⁻⁵ in the last epoch, i.e.,

$\begin{array}{ll}{\Vert \frac{{\Theta }_{n+1}-{\Theta }_n}{{\Theta }_{ref}}\Vert }_{\infty }<10^{-5}& (7)\end{array}$

where Θ_i is the estimate of the parameters at the ith epoch, Θ_ref is a vector with reference values for each parameter, and vector division is the element-wise division of the vector components.

3. Results

In this section, we analyse the ability of both formulations (see Equations 3 and 4) to identify myocardial constitutive parameters, and the intra-ventricular pressure. In our assessment, we disregard the representation error of the domain associated with the PCA and learning space presented in Maso Talou et al. (2020), as it can be reduced by including additional PCA modes as inputs of our surrogate model. Without loss of generality, we fixed the geometry to g = 0 (i.e., mean ventricular shape across our population, see Figure 1). The analyses reported here will be analogous for different geometrical variations of the left ventricle as well as for different Dirichlet boundary conditions. All observations $\widehat{\text{u}}$, ${\text{s}}_f^{\mathrm{\text{endo}}}$ and ${\text{s}}_f^{\mathrm{\text{epi}}}$ in this work have been generated using the same FE model (see section 2.1) to train the surrogate model, thus there is no model error in this study.

FIGURE 1

Figure 1. (left) AI-surrogate depicting its inputs and outputs (Maso Talou et al., 2020); (right) Example of the geometry of the ventricle at diastasis for g = 0.

Additionally, we studied the effect of tracking and segmentation errors (observation errors) on the estimates of interest. The observation error is modelled as

$\begin{array}{ll}{\text{u}}_{ϵ}=\frac{\text{r}}{\parallel \text{r}\parallel }\eta & (8)\end{array}$

where r is a randomly oriented vector with components $r_i\in {U}(0,1)$ and $\eta \in {N}(0,\sigma )$ is the magnitude of u_ϵ. Note that the magnitude of our error is normally distributed with zero-mean, and we choose σ = 1 mm to represent errors of similar magnitude to the image resolution of MRI images for myocardial structures. This error is added to $\widehat{\text{u}}$, ${\text{s}}_f^{\mathrm{\text{endo}}}$, and ${\text{s}}_f^{\mathrm{\text{epi}}}$ for the cases with noise.

Regarding spatial discretisation of the observations (i.e., the number of samples for $\widehat{\text{u}}$, ${\text{s}}_f^{\mathrm{\text{endo}}}$, and ${\text{s}}_f^{\mathrm{\text{epi}}}$), we sampled 1, 109 points across the cardiac wall and 7, 579 points over each of the ventricular surfaces to describe $\widehat{\text{u}}$ and ${\text{s}}_f^{·}$, respectively. In turn, for the spatial discretisation of the AI-surrogate inputs, we sampled 1, 109 points across the cardiac wall for the full-field tracking approach and 1, 072 points to describe each of the surfaces ${\text{s}}_i^{·}$. By oversampling ${\text{s}}_f^{·}$ with respect to ${\text{s}}_i^{·}$, we can compute the distance between these surfaces by substituting ${\text{s}}_1={\text{s}}_f^{·}$ into Equation (5). This approximation yields a reasonable trade-off between computation time and discretisation error (note that the discretisation error decreases as |s₁| increases).

3.1. Identification of Constitutive Parameters

We first studied the precision of simultaneously estimating the constitutive parameters c₁ and c₂ for a known intra-ventricular pressure p. The estimates were generated for both approaches under three experimental scenarios: (i) a two time-point experiment without observation error, i.e., only two images are available during the cardiac cycle; (ii) a two time-point experiment with observation error; and (iii) a multiple time-point (N = 10) experiment with observation error. The first two experiments illustrate the degradation of c₁ and c₂ estimates for both approaches as measurements become less reliable due to observation errors. This indicates the robustness of the approaches with respect to observation error. A comparison of the last two experiments shows the benefits of including additional observations (time-points), which diminish the impact of the stochastic component of the observation error.

3.1.1. Experiment 1

Both approaches showed a similar accuracy in the absence of observation errors. Since variations in the predicted displacements are more sensitive with respect to c₁ than c₂, we observed better estimates of c₁ (see Figure 3). The distribution of the error in the parametric space (c₁, c₂) does not show an association between the displacement magnitudes and the estimation error (compare Figures 2, 3).

FIGURE 2

Figure 2. Mean displacement of the ventricular wall (denoted ${\mu }_{\widehat{\text{u}}}$) for different combinations of the constitutive parameters c₁ and c₂, and an intra-ventricular pressure load of p = 0.9 kPa. For details of mean displacements for other loads, refer to the Supplementary Material.

FIGURE 3

Figure 3. Relative error of the two-parameter estimation. Estimates of the constitutive parameters c₁ (left) and c₂ (right) are presented using the full-field tracking and contour matching approaches with a load of p = 0.9 kPa in the absence of noise in the observations. For details of estimations for other loads, refer to the Supplementary Material.

3.1.2. Experiment 2

To analyse the effect of observation errors, we repeated the previous study 10 times while adding different independent instances of u_ϵ (see Equation 8) to the observations $\widehat{\text{u}}$, ${\text{s}}_f^{\mathrm{\text{endo}}}$ and ${\text{s}}_f^{\mathrm{\text{epi}}}$. To quantify the effect of the noise on the estimates across these different samples, we computed the mean and standard deviation of the relative error when estimating the parameters c₁ and c₂ (see Figure 4).

FIGURE 4

Figure 4. Mean and standard deviation of the relative error for the two-parameter estimation of c₁ (left) and c₂ (right) constitutive parameters using the full-field tracking and contour matching approaches with a load of p = 0.9 kPa and the presence of Gaussian noise (σ = 1mm ≈ 2 pixels) in the observations. For details of estimations for other loads, refer to the Supplementary Material.

The contour matching approach provided more accurate estimates than the full-field tracking approach in the presence of observation errors (see Figure 4). This is shown by the lower overall error of the contour matching approach, which is insensitive to displacements within the muscular heart wall, and within the surfaces of the ventricles. As only a small tangential component of displacement is expected in the endocardial and epicardial surfaces from the loading process in our FE model, the insensitivity toward such component is beneficial for the estimation process. Note that this may not be true in practical scenarios, where the quantification of such components is yet to be thoroughly explored.

The full-field tracking approach showed larger estimation errors for c₁ in stiffer materials. In such cases, the displacement-to-u_ϵ ratio (analogous to the signal-to-noise ratio in signal processing) was smaller, hindering the estimation due to the increasing effects of noise in the observations. c₂ presented lower identifiability due to the plateau in the objective function, ${L}(\Theta )$, associated with the observation error (see Equation 3 and Figure 9). The sensitivity of c₂ with respect to the displacements was higher than that of c₁ due to the form of the constitutive equation (Equation 1). The identifiability of c₂ is susceptible to surrogate approximation errors, i.e., discrepancies between the AI-surrogate and the FE model. An example of this, is the increase in error of c₂ as shown in Figure 4. The AI-surrogate underestimates ∂u/∂c₂ in its predictions contributing to the plateau in ${L}(\Theta )$ described previously. Note that, as the pressure is increased the AI-surrogate error has a diminishing effect on the predicted displacement (see Supplementary Figures 7, 8).

3.1.3. Experiment 3

We analysed the effect of adding multiple time-points for the estimation of the constitutive parameters when observation errors were present. We assumed that the observation errors, u_ϵ, for all time-points, were independent and identically distributed (uncorrelated).

We observed that both approaches showed an improvement in reducing the mean relative error in the estimation of c₁ and c₂ (see Figure 5), closer to values from error-free estimates reported in Figure 3. As the observation error is normal, zero-mean and independent between the time-points, both formulations efficiently cope with the uncertainty, because the Euclidean norm involved in both formulations (Equations 3 and 5) optimises toward the mean of the error distribution.

FIGURE 5

Figure 5. Mean and standard deviation of the relative error for the two-parameter estimation of c₁ (left) and c₂ (right) constitutive parameters using full-field tracking and contour matching approaches, a pressure trace with 10 time-points, and independent Gaussian noise (σ = 1mm ≈ 2 pixels) in the observations.

The contour matching approach is slightly more accurate, mainly in the estimation of c₂ (maximum error of 3%, in comparison to the full-field tracking approach, with a maximum error of 5 ± 4 %).

3.2. Identification of Intra-Ventricular Pressure

The previous section studied the accuracy of estimating the constitutive parameters, c₁ and c₂, given a known intra-ventricular pressure, p. In this section, we analogously analyse the accuracy in the estimation of p, assuming that the constitutive parameters are known. Two experiments were conducted: (i) a two time-point experiment without observation error where only two images are available; and (ii) a two time-point experiment with observation error. Analogous to the previous section, we studied the degradation in the ability to recover the intra-ventricular pressure for both approaches as the observations become less reliable.

3.2.1. Experiment 1

Both approaches estimated intra-ventricular pressure p with a similar accuracy (relative error ≤ 1%, see Figure 6). The results did not present evidence of a correlation between the magnitude of the displacements and the error in the estimation of p.

FIGURE 6

Figure 6. Relative error in estimating intra-ventricular pressure p using full-field tracking and contour matching approaches with a load of p = 0.9 kPa in the absence of noise in the observations. For details of the estimation of other loads, refer to the Supplementary Material.

3.2.2. Experiment 2

Introducing observation errors slightly reduced the accuracy of the full-field tracking approach (see Figure 7), especially for stiffer materials. The interpretation of these results is analogous to the analysis presented for the estimation of c₁ in Experiment 2 of section 3.1. As reported in the previous section, with increasing material stiffness, the displacements of the ventricle for the same intra-ventricular pressure are smaller and the same intensity of noise will have more detrimental effects on the estimation process (due to less displacement-to-u_ϵ ratio in the observations).

FIGURE 7

Figure 7. Mean and standard deviation of the relative errors for recovering ground truth intra-ventricular pressure using the full-field tracking and contour matching approaches and an intra-ventricular pressure of p = 0.9 kPa with Gaussian noise (σ = 1 mm ≈ 2 pixels) in the observations. Refer to the Supplementary Material for details of the relative errors for other intra-ventricular pressures.

3.3. Convexity of the Objective Function

Both approaches are based on measuring the displacement error by means of the Euclidean distance between the ground truth and the AI-surrogate prediction. As the mechanical model monotonically increases the ventricular displacements with the increase of p and the decrease of c₁ and c₂, the minimisation problem defined in Equations (3) or (4) with respect to p, c₁ or c₂ is convex if the AI-surrogate sufficiently approximates the FE model response.

The AI-surrogate used in this manuscript, satisfies such a condition (see Figures 8, 9). We can observe convexity in the loss function with respect to p and the constitutive parameters c₁ and c₂. Because we use a stochastic gradient descent optimiser, the identifiability of the parameters is related to the rate of change of the loss function with respect to the parameter of interest. We observed a decrease in $\frac{\partial {L}}{\partial c_1}$ and $\frac{\partial {L}}{\partial c_2}$ for stiffer materials. Nonetheless, both approaches presented good accuracy in recovering physiologically realistic ranges of parameters (2.0 kPa ≤ c₁ ≤ 5.0 kPa and 4.0 ≤ c₂ ≤ 40.0, even in the presence of noise (see Figures 8, 9). The identifiability of p does not exhibit degradation in the physiological range, even in the presence of noise (see Figure 8).

FIGURE 8

Figure 8. Loss function for the full-field tracking and contour matching approaches with respect to p in the presence (right) and absence (left) of noise in the observations. The surface plot in the p-$\widehat{p}$ plane and the opaque manifold represent the loss function value ${L}(p)$ when $\widehat{p}$ is the ground truth (optimal) value. Level curves in ${L}$-$\widehat{p}$ and ${L}$-p planes correspond to parallel cuts of the opaque manifold for fixed values of p and $\widehat{p}$, respectively (blue to red shades indicate values of p and $\widehat{p}$ that correspond to the ticks on their respective axes). Note that for a given $\widehat{p}$, ${L}(p)$ is convex with the minimum $p=\widehat{p}$ resulting in precise and well-behaved formulations for gradient-based optimisers.

FIGURE 9

Figure 9. Loss function for the full-field tracking and contour matching approaches with respect to c₁ and c₂ in the presence (right) and absence (left) of noise in the observations. The first and second rows present the loss function values when c₁ = 3.5 kPa, c₂ = 22.0 and c₁ = 2.5 kPa, c₂ = 10.0 are the ground truth parameters, respectively. Surface plots in c₁-c₂ plane and the opaque manifolds represent the loss function value ${L}(c_1,c_2)$. Level curves in ${L}$-c₁ and ${L}$-c₂ planes correspond to parallel cuts of the opaque manifold for fixed values of c₂ and c₁, respectively (blue to red shades indicate values of c₁ and c₂ that correspond to the ticks on their respective axes). Note that the objective functions are more sensitive to perturbations of c₁ and c₂ for softer materials (as highlighted by the blue level curves, which represent lower values of c₁ and c₂).

3.4. Computational Cost

We quantify the performance of the formulations by measuring the wall-clock time taken to estimate four different set of parameters: (i) intra-ventricular pressure; (ii) constitutive parameter c₁; (iii) constitutive parameter c₂; and (iv) constitutive parameters c₁ and c₂. For each test, we executed the estimation 50 times with different observation errors (using the same noise distribution presented in section 3). During these tests, the observation $\widehat{\text{u}}$ corresponded to a ground truth FE model prediction obtained using p = 0.9 kPa, c₁ = 3.5 kPa and c₂ = 22.0. We then estimated the mean and standard deviation of the wall-clock times over the 50 executions (see Table 1). The networks and optimisation schemes were implemented in TensorFlow 2.1 with GPU-support, using an NVIDIA Quadro P6000 and CUDA v11.0 for their executions.

TABLE 1

Table 1. Wall-clock times for parameter estimation using the proposed strategies.

Both approaches were able to estimate the different set of parameters in less than 18 s. In particular, the full-field tracking approach executed 6.2, 7.3, 6.4, and 3.6 times faster for experiments (i)-(iv), respectively, than the contour matching counterpart. The slower response of the contour matching is due to the oversampling of the ventricular surface ${\text{s}}_f^{·}$ that is necessary for obtaining a good approximation of Equation (5). Such oversampling increases the number of points in s₂ requiring more model evaluations and thus computational expense.

In terms of computational complexity, the full-field tracking approach has a cost of ${O}(|\widehat{\text{u}}|)$ where $|\widehat{\text{u}}|$ is the number of points across the ventricular wall, and the contour matching approach has a cost of ${O}(|{\text{s}}_1||{\text{s}}_2|)$ where |s_i| is the number of points at the ventricular surfaces s_i. This analysis shows a higher computational complexity for the contour matching approach, which detrimentally impacts the scaling of the approach with respect to the discretisation of the ventricle. Specifically, the computational cost grows linearly and quadratically for the full-field tracking and contour matching approaches, respectively.

Another contribution to the lower performance of the contour matching approach is given by the number of network evaluations. The total number of evaluations of $\text{u}(\widehat{\Theta })$ scales linearly with respect to the discretisation of the ventricle for both formulations. However, contour matching requires almost twice the number of evaluations when compared to the full-field tracking approach (1109 evaluations for each contour, vs. 1072 evaluations for tracking). Thus, less refined representations of the ventricular surfaces may help reduce the computational expense if needed for clinical translation.

It is worth noting that we optimised Equations (3) and (4) using the ADAM algorithm, because of its direct support toward neural network optimisation. It is possible that the use of efficient optimisers for convex problems, such as L-BFGS (Liu and Nocedal, 1989) or even Newton's method, may further enhance convergence and reduce the computational effort in the proposed formulations.

4. Discussion

The two techniques proposed for parameter estimation feature appealing properties such as low computational cost, simple implementation, and no need for analytical derivatives of the objective function. To achieve this, we used automatic differentiation, already implemented in neural networks frameworks (such as Tensorflow and PyTorch), for solving Equations (3) and (4). This allow us to use exact gradient information from complex AI-surrogate models without any additional effort. Thus, we were able to assess the sensitivity of the objective function residuals with respect to its inputs, allowing for an efficient convex optimisation of the inputs. Additionally, the same neural network frameworks are endowed with GPU-efficient implementations, accelerating the evaluation of models by orders of magnitude with respect to a cost equivalent CPU infrastructure.

Regarding the estimation of the intra-ventricular pressure and constitutive parameters of the left ventricle, we conclude that both of the proposed approaches can provide accurate predictions of parameters, even in the presence of measurement noise. The noise was modelled as a normal zero-mean distribution with a standard deviation of two pixels. This noise represented reconstruction errors of the ventricle wall displacement (full-field tracking approach), or the ventricular surface geometry (contour matching approach). Note that in our analysis, the errors present no biases. If a bias were present (e.g., consistent segmentation errors due to mis-identification of the structures, or assimilating experimental data with a model with a significant modelling error), this may lead to larger estimation errors than those reported here.

Both approaches demonstrated good identification properties for the physiological range considered in our experiments. Particularly, we focused on the analysis of relatively stiffer materials (c₁ > 2.0 kPa) and lower intraventricular pressures (p = 0.9 kPa) because it present a more challenging scenario to assess the parameter estimation task. As displacements are smaller for lower pressures and stiffer materials, the noise and the displacement field are within the same displacement magnitudes, and a lower performance for the estimation is expected. However, the results only showed this degradation under those conditions for the full-field tracking approach, increasing its error estimate from 0.6 to 1.8% in c₁, from 1.9 to 4.6% in c₂ and from 0.3 to 2.7% in p (see Figures 5, 7). The contour matching approach only presented a slight degradation in parameter estimates of c₂ when c₁ reached the upper bound of its physiological range. Note that, after estimating the constitutive parameters, the displacement errors between the AI-surrogate predictions and the target observations were visually negligible when the corresponding ventricle contours were overlaid. The case with the largest disagreement across all simulations reported here is shown in Figure 4 for the full-field approach with target parameters c₁ = 3.5 kPa and c₂ = 10.0 (the specific sample estimated c₁ = 3.76 kPa and c₂ = 7.05). In this case, the displacements of the AI-surrogate had an error of 0.06 ± 0.02 mm with respect to the FE ground truth. In particular, when not using information from multiple frames, the full-field tracking approach failed to precisely identify c₂ with errors reaching 16 ± 9.6%. Nonetheless, in practical applications, clinical MR datasets often contain enough temporal resolution to perform a multi-frame kinematic assessment.

For normal human physiological cases (diastasis pressures of approximately 1.2 kPa), both methods present a slight improvement in performance due to the larger displacements in the observations (see Supplementary Figures 3, 7, 11, 15). As the observed displacements increase, the approximation error of the AI-surrogate to the FE model has a smaller contribution, leading to more accurate estimates. Specifically, we can observe this for the estimates of stiffer materials, where the larger displacements improved the parameter identification in comparison with the 0.9 kPa pressure case.

Regarding clinical translation, this approach offers rapid and efficient estimation of the mechanical properties using commodity computational resources, e.g., a standard computer with state-of-the-art GPU. The generated training datasets used in section 3 assume realistic data constraints (i.e., observation errors and resolution) expected in medical data. In this study, we assumed segmentation errors of two pixels, and a temporal resolution of 10 frames during diastole (see section 3), which are both attainable using a 3T clinical MRI scanner. FE models (which are used to train our AI-surrogate) have demonstrated clinical utility (Wang et al., 2018; Hasaballa et al., 2021), evidencing the suitability of this approach for mechanical characterisation of the heart. Nonetheless, robust uncertainty quantification analysis should be performed to analyse all sources of error in the specific clinical environment. For the assessment of mechanical properties of the ventricular wall, it is important to quantify the uncertainty in the pressure measurements, segmentation error, and geometric representation error. The characterisation of such uncertainties is out of the scope of the present work.

In terms of the number of parameters to be estimated, our study demonstrated reasonable efficiency for the simultaneous estimation of two parameters. Compared to the single parameter problem, the computational time for the two parameter problem was more than two-fold, but still presenting time ranges compatible with clinical practice. This is due to the coupled effect of the parameters in the model response which is a problem shared by all parameter estimation techniques. The extension of our technique for the simultaneous identification of additional parameters is possible, as long as the objective function remains strictly convex with respect to the parameters of interest.

While the estimation approaches considered in this work involved solving convex optimisation problems, the use of AI-surrogates can offer advantages for non-convex problems. In such cases, non-convex optimisation solvers, such as genetic algorithms, may benefit from the use of AI-surrogates due to the computationally inexpensive model evaluations. This feature enables a more efficient exploration of the parameter space with reduced computational intensity.

Finally, both of the proposed estimation techniques are general, and can be translated straightforwardly to other applications (e.g., to estimate constitutive properties of other tissues, such as the breast, lung or liver) as long as an AI-surrogate can be generated from the appropriate models. The full-field tracking approach is limited to applications where material point tracking measurements are available (e.g., using CMR tissue tagging, or image registration techniques). On the other hand, the contour matching approach can be applied to any applications where the surfaces of the tissues or organs of interest can be quantified experimentally.

5. Conclusions

This study proposed two approaches for parameter estimation using AI-surrogates, depending on whether (i) tracking kinematic measurements, or (ii) only surface measurements are available. We focused our application on the estimation of left ventricular constitutive properties and its intra-ventricular pressure during the passive filling phase of the cardiac cycle.

We conclude that: (i) both approaches are robust with respect to Gaussian noise when the measurement data for multiple loading conditions were combined; and (ii) estimates of one or two constitutive parameters could be obtained in less than 9 or 18 s, respectively. We found that the contour matching approach was more robust toward Gaussian noise, recovering the ground truth parameters with high accuracy even when only one loaded configuration was available. Conversely, the full-field tracking approach was more efficient than its counterpart by a factor of ≈4, while providing the possibility of further improving scalability as medical imaging resolution improves.

Data Availability Statement

The data analyzed in this study is subject to the following licenses/restrictions: The de-identified data is available from the authors upon reasonable request, with permission of the institutional review board and under ethical approval. Requests to access these datasets should be directed to Martyn P. Nash, bWFydHluLm5hc2gmI3gwMDA0MDthdWNrbGFuZC5hYy5ueg==.

Author Contributions

GM co-designed the study, implemented, validated and executed the optimisation techniques, and prepared the manuscript. TB co-designed the study, developed the biomechanical modelling framework, and contributed to the manuscript text. MN contributed to the design of the study, the development of the biomechanical modelling framework, and edited the manuscript. All authors contributed to the discussions presented in this manuscript.

Funding

This work was supported by research funding grants 13/317 and 17/608 from the Health Research Council of New Zealand, and funding grant 9077/31/8402 from the Li Ka Shing Foundation.

Conflict of Interest

GM and TB are scientific advisors at HeartLab (NZ) Ltd. MN is Chief Scientific Officer at HeartLab (NZ) Ltd.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2021.732351/full#supplementary-material

References

Aletras, A. H., Ding, S., Balaban, R. S., and Wen, H. (1999). DENSE: displacement encoding with stimulated echoes in cardiac functional MRI. J. Magnet. Reson. 137, 247–252. doi: 10.1006/jmre.1998.1676

PubMed Abstract | CrossRef Full Text | Google Scholar

Amzulescu, M. S., De Craene, M., Langet, H., Pasquet, A., Vancraeynest, D., Pouleur, A. C., et al. (2019). Myocardial strain imaging: review of general principles, validation, and sources of discrepancies. Eur. Heart J. Cardiovasc. Imaging 20, 605–619. doi: 10.1093/ehjci/jez041

PubMed Abstract | CrossRef Full Text | Google Scholar

Bradley, C., Bowery, A., Britten, R., Budelmann, V., Camara, O., Christie, R., et al. (2011). OpenCMISS: a multi-physics & multi-scale computational infrastructure for the VPH/Physiome project. Prog. Biophys. Mol. Biol. 107, 32–47. doi: 10.1016/j.pbiomolbio.2011.06.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Cai, L., Ren, L., Wang, Y., Xie, W., Zhu, G., and Gao, H. (2021). Surrogate models based on machine learning methods for parameter estimation of left ventricular myocardium. R. Soc. Open Sci. 8:201121. doi: 10.1098/rsos.201121

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, C., Qin, C., Qiu, H., Tarroni, G., Duan, J., Bai, W., et al. (2020). Deep learning for cardiac image segmentation: a review. Front. Cardiovasc. Med. 7:25. doi: 10.3389/fcvm.2020.00025

PubMed Abstract | CrossRef Full Text | Google Scholar

Dabiri, Y., Van der Velden, A., Sack, K. L., Choy, J. S., Kassab, G. S., and Guccione, J. M. (2019). Prediction of left ventricular mechanics using machine learning. Front. Phys. 7:117. doi: 10.3389/fphy.2019.00117

PubMed Abstract | CrossRef Full Text | Google Scholar

Davies, V., Noe, U., Lazarus, A., Gao, H., Macdonald, B., Berry, C., et al. (2019). Fast parameter inference in a biomechanical model of the left ventricle by using statistical emulation. J. R. Stat. Soc. Ser. C 68, 1555–1576. doi: 10.1111/rssc.12374

PubMed Abstract | CrossRef Full Text | Google Scholar

Di Achille, P., Harouni, A., Khamzin, S., Solovyova, O., Rice, J. J., and Gurev, V. (2018). Gaussian process regressions for inverse problems and parameter searches in models of ventricular mechanics. Front. Physiol. 9:1002. doi: 10.3389/fphys.2018.01002

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, H., Li, W. G., Cai, L., Berry, C., and Luo, X. Y. (2015). Parameter estimation in a Holzapfel-Ogden law for healthy myocardium. J. Eng. Math. 95, 231–248. doi: 10.1007/s10665-014-9740-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Guccione, J. M., McCulloch, A. D., and Waldman, L. K. (1991). Passive material properties of intact ventricular myocardium determined from a cylindrical model. J. Biomech. Eng. 113, 42–55. doi: 10.1115/1.2894084

PubMed Abstract | CrossRef Full Text | Google Scholar

Hadjicharalambous, M., Asner, L., Chabiniok, R., Sammut, E., Wong, J., Peressutti, D., et al. (2017). Non-invasive model-based assessment of passive left-ventricular myocardial stiffness in healthy subjects and in patients with non-ischemic dilated cardiomyopathy. Ann. Biomed. Eng. 45, 605–618. doi: 10.1007/s10439-016-1721-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Hasaballa, A. I., Babarenda Gamage, T. P., Wang, V. Y., Zhao, D., Mauger, C. A., Gilbert, K., et al. (2021). “Sensitivity of myocardial stiffness estimates to inter-observer variability in LV geometric modelling,” in Functional Imaging and Modeling of the Heart, eds D. B. Ennis, L. E. Perotti, and V. Y. Wang (Cham: Springer International Publishing), 287–295. doi: 10.1007/978-3-030-78710-3_28

CrossRef Full Text | Google Scholar

Ibrahim, E.-S. H. (2011). Myocardial tagging by Cardiovascular Magnetic Resonance: evolution of techniques-pulse sequences, analysis algorithms, and applications. J. Cardiovasc. Magnet. Reson. 13:36. doi: 10.1186/1532-429X-13-36

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, D. C., and Nocedal, J. (1989). On the limited memory BFGS method for large scale optimization. Math. Program. 45, 503–528. doi: 10.1007/BF01589116

CrossRef Full Text | Google Scholar

Longobardi, S., Lewalle, A., Coveney, S., Sjaastad, I., Espe, E. K., Louch, W. E., et al. (2020). Predicting left ventricular contractile function via gaussian process emulation in aortic-banded rats. Philos. Trans. R. Soc. A 378:20190334. doi: 10.1098/rsta.2019.0334

PubMed Abstract | CrossRef Full Text | Google Scholar

Maso Talou, G. D., Babarenda Gamage, T. P., Sagar, M., and Nash, M. P. (2020). Deep learning over reduced intrinsic domains for efficient mechanics of the left ventricle. Front. Phys. 8:30. doi: 10.3389/fphy.2020.00030

CrossRef Full Text | Google Scholar

Nash, M. P., and Hunter, P. J. (2000). Computational mechanics of the heart. J. Elasticity Phys. Sci. Solids 61, 113–141. doi: 10.1023/A:1011084330767

CrossRef Full Text | Google Scholar

Nielsen, P. M., Le Grice, I. J., Smaill, B. H., and Hunter, P. J. (1991). Mathematical model of geometry and fibrous structure of the heart. Am. J. Physiol. Heart Circ. Physiol. 260, H1365–H1378. doi: 10.1152/ajpheart.1991.260.4.H1365

PubMed Abstract | CrossRef Full Text | Google Scholar

Noe, U., Lazarus, A., Gao, H., Davies, V., Macdonald, B., Mangion, K., et al. (2019). Gaussian process emulation to accelerate parameter estimation in a mechanical model of the left ventricle: a critical step towards clinical end-user relevance. J. R. Soc. Interface 16:20190114. doi: 10.1098/rsif.2019.0114

PubMed Abstract | CrossRef Full Text | Google Scholar

Queiros, S., Vilaca, J. L., Morais, P., Fonseca, J. C., D'hooge, J., and Barbosa, D. (2017). Fast left ventricle tracking using localized anatomical affine optical flow. Int. J. Num. Methods Biomed. Eng. 33:e2871. doi: 10.1002/cnm.2871

PubMed Abstract | CrossRef Full Text | Google Scholar

Raissi, M., Perdikaris, P., and Karniadakis, G. E. (2019). Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. J. Comput. Phys. 378, 686–707. doi: 10.1016/j.jcp.2018.10.045

CrossRef Full Text | Google Scholar

Rumindo, G. K., Ohayon, J., Croisille, P., and Clarysse, P. (2020). In vivo estimation of normal left ventricular stiffness and contractility based on routine cine MR acquisition. Med. Eng. Phys. 85, 16–26. doi: 10.1016/j.medengphy.2020.09.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Shi, W., Zhuang, X., Wang, H., Duckett, S., Luong, D. V. N., Tobon-Gomez, C., et al. (2012). A comprehensive cardiac motion estimation framework using both untagged and 3-D tagged MR images based on nonrigid registration. IEEE Trans. Med. Imaging 31, 1263–1275. doi: 10.1109/TMI.2012.2188104

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, H., and Amini, A. A. (2011). Cardiac motion and deformation recovery from MRI: a review. IEEE Trans. Med. Imaging 31, 487–503. doi: 10.1109/TMI.2011.2171706

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Z. J., Wang, V. Y., Bradley, C. P., Nash, M. P., Young, A. A., and Cao, J. J. (2018). Left ventricular diastolic myocardial stiffness and end-diastolic myofibre stress in human heart failure using personalised biomechanical analysis. J. Cardiovasc. Transl. Res. 11, 346–356. doi: 10.1007/s12265-018-9816-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Wang, V. Y., Morgan, A. E., Kim, J., Tafreshi, R., Wallace, A. W., et al. (2021). Finite-element based optimization of left ventricular passive stiffness in normal volunteers and patients after myocardial infarction: Utility of an inverse deformation gradient calculation of regional diastolic strain. J. Mech. Behav. Biomed. Mater. 119:104431. doi: 10.1016/j.jmbbm.2021.104431

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhong, X., Spottiswoode, B. S., Meyer, C. H., Kramer, C. M., and Epstein, F. H. (2010). Imaging three-dimensional myocardial mechanics using navigator-gated volumetric spiral cine DENSE MRI. Magnet. Reson. Med. 64, 1089–1097. doi: 10.1002/mrm.22503

PubMed Abstract | CrossRef Full Text | Google Scholar