High-dimensional Bayesian Optimization of Personalized Cardiac Model Parameters via an Embedded Generative Model

by   Jwala Dhamala, et al.

The estimation of patient-specific tissue properties in the form of model parameters is important for personalized physiological models. However, these tissue properties are spatially varying across the underlying anatomical model, presenting a significance challenge of high-dimensional (HD) optimization at the presence of limited measurement data. A common solution to reduce the dimension of the parameter space is to explicitly partition the anatomical mesh, either into a fixed small number of segments or a multi-scale hierarchy. This anatomy-based reduction of parameter space presents a fundamental bottleneck to parameter estimation, resulting in solutions that are either too low in resolution to reflect tissue heterogeneity, or too high in dimension to be reliably estimated within feasible computation. In this paper, we present a novel concept that embeds a generative variational auto-encoder (VAE) into the objective function of Bayesian optimization, providing an implicit low-dimensional (LD) search space that represents the generative code of the HD spatially-varying tissue properties. In addition, the VAE-encoded knowledge about the generative code is further used to guide the exploration of the search space. The presented method is applied to estimating tissue excitability in a cardiac electrophysiological model. Synthetic and real-data experiments demonstrate its ability to improve the accuracy of parameter estimation with more than 10x gain in efficiency.



There are no comments yet.


page 9


Bayesian Optimization on Large Graphs via a Graph Convolutional Generative Model: Application in Cardiac Model Personalization

Personalization of cardiac models involves the optimization of organ tis...

Generative Modeling and Inverse Imaging of Cardiac Transmembrane Potential

Noninvasive reconstruction of cardiac transmembrane potential (TMP) from...

Fast Posterior Estimation of Cardiac Electrophysiological Model Parameters via Bayesian Active Learning

Probabilistic estimation of cardiac electrophysiological model parameter...

Reducing local minima in fitness landscapes of parameter estimation by using piecewise evaluation and state estimation

Ordinary differential equations (ODE) are widely used for modeling in Sy...

Bayesian multiscale deep generative model for the solution of high-dimensional inverse problems

Estimation of spatially-varying parameters for computationally expensive...

Accelerating Derivative-Free Optimization with Dimension Reduction and Hyperparameter Learning

We consider convex, black-box objective functions with additive or multi...

Adaptive Estimation of the Neural Activation Extent in Computational Volume Conductor Models of Deep Brain Stimulation

Objective: The aim of this study is to propose an adaptive scheme embedd...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Patient-specific simulation models of the heart have shown increasing potential in personalized treatment of cardiac diseases  [2, 15]. The estimation of patient-specific tissue properties, however, remains an unresolved problem. One significant challenge is that the unknown tissue properties (in the form of model parameters) are spatially varying at a resolution associated with the discrete cardiac mesh. To estimate these high-dimensional (HD) model parameters is not only algorithmically difficult given indirect and sparse measurements, but also computationally intractable in the presence of computing-intensive simulation models.

Numerous efforts have been made to circumvent the challenge of HD parameter estimation. Many works assume homogeneous tissue property throughout the myocardium, which can be represented by a single global model parameter [9, 11]. To preserve local information about the spatially distributed tissue properties, most existing works resort to dimensionality reduction through an explicit partitioning of the cardiac mesh. These efforts can be generally summarized in two categories. In the first approach, the cardiac mesh is pre-divided into 3-26 segments, each represented by a uniform parameter value [17, 18]. Naturally, this artificial low-resolution division of the cardiac mesh has a limited ability to represent the underlying tissue heterogeneity that is not known a priori. In addition, it has been shown that the initialization of model parameters becomes increasingly more critical as the number of segments grows [17]. In an alternative approach, the explicit partitioning of the cardiac mesh is done through a coarse-to-fine optimization along a pre-defined multi-scale hierarchy of the cardiac mesh, enabling spatially-adaptive resolution of tissue properties that is higher in certain regions than the others [4, 6, 7, 5]. However, the representation ability of the final partition is limited by the inflexibility of the pre-defined multi-scale hierarchy: homogeneous regions distributed across different scales cannot be grouped into the same partition, while the resolution of heterogeneous regions can be limited by the level of the scale the optimization can reach [6]. In addition, because these methods involve a cascade of optimizations along the coarse-to-fine hierarchy of the cardiac mesh, they are computationally expensive. In the context of parameter estimation for models that could require hours or days for a single simulation, these methods could quickly become computationally prohibitive.

In this paper, we present a novel HD parameter optimization approach that replaces the explicitly defined low-dimensional (LD) or multi-scale representation of the parameter space with an implicit LD latent encoding, achieved by embedding within the optimization a stochastic LD-to-HD generative model that describes the generation of the HD spatially-varying tissue properties from a LD manifold. This generative model is obtained with a variational auto-encoder (VAE), trained from a large set of spatially-varying tissue properties reflecting regional tissue abnormality with various locations, sizes, and distributions. Once trained, the VAE is integrated with a Bayesian optimization (BO) [3] framework in two novel ways. First, the generative model (the VAE decoder) is embedded within the objective function to provide an implicit LD search space for the optimization of HD parameters. Second, the posterior distribution of the LD latent code as learned from the VAE encoder is used as prior knowledge within the BO for an efficient exploration of the LD manifold. To the best of our knowledge, this is the first work that utilizes a probabilistic generative process within an optimization framework for estimating HD patient-specific model parameters.

The presented method is applied to the estimation of local tissue excitability of a cardiac electrophysiological model using non-invasive electrocardiogram (ECG) data. On both synthetic and real data experiments, the presented method is compared against existing methods based on explicitly-defined LD [17] or multi-scale representation of the parameter space [6]. Experiments demonstrate that the presented method can achieve a drastic reduction in computational cost while improving the accuracy of the estimated parameters. Beyond the specific model considered in this paper, the presented method provides an efficient and reliable solution to a wider range of HD model parameter estimation problems.

2 Background: Cardiac Electrophysiological System

Cardiac Electrophysiology Model: Among the different types of cardiac electrophysiological models, phenomenological models such as the Aliev-Panfilov (AP) model [1] can explain the macroscopic process of cardiac excitation with a small number of model parameters and reasonable computation, making their use in parameter estimation widespread [15, 9, 6]. Therefore, in this study, the AP model given below is chosen to test the feasibility of the presented method:


Here, parameter controls the coupling between the normalized transmembrane action potential and the recovery current ,

is the diffusion tensor,

controls the repolarization, and controls the excitability of the cell. The transmural action potential is computed by solving the AP model (1) on a 3D myocardium discretized using the meshfree method [16]. Here, the output is most sensitive to the value of the parameter [6], which is associated to the ischemic severity of the myocardial tissue. Therefore, as an initial study, we focus on the estimation of the spatially distributed parameter .
Body-surface ECG Model: The propagation of the spatio-temporal transmural action potential U to the potentials measured on the body surface Y can be described by the quasi-static approximation of the electromagnetic theory [13]. Solving the governing equations on a discrete heart-torso mesh, a linear relationship between U and Y can be obtained as: , where

is the vector of local parameters

at the resolution of the cardiac mesh.

3 HD Parameter Estimation via an Embedded Generative Model

The parameter in the AP model (1) can estimated by maximizing the similarity between the measured ECG signal Y and those simulated by the combined cardiac electrophysiological and surface ECG model :


The direct estimation of at the resolution of the cardiac mesh is infeasible, both due to limited identifiability from limited ECG data and prohibitively high computational cost. The presented method enables this HD optimization by embedding within the Bayesian optimization framework a stochastic generative model that generates the HD parameter on the cardiac mesh from a LD manifold. As outlined in Fig. 1, the presented method includes two major components: 1) the construction of a generative model of HD parameters representing spatially-varying tissue properties at the resolution of the cardiac mesh, and 2) a novel Bayesian optimization method utilizing the embedded generative model.

Figure 1: The workflow diagram of the presented high dimensional Bayesian optimization via embedded variational auto-encoder for local parameter estimation.

3.1 LD-to-HD Parameter Generation via VAE

Recently, generative models have shown promising potential in unsupervised learning of abstract LD representations from which complex images can be generated 

[10]. Inspired by this, here, we utilize a VAE to obtain a stochastic generative model of the HD spatially varying tissue properties from a LD manifold.

Generative VAE model: We assume that the spatially varying tissue properties at the resolution of a cardiac mesh

is generated by a small number of unobserved continuous random variables

z in a LD manifold. To obtain the generative process from z to , the VAE consists of two modules: a probabilistic deep encoder network with parameters that approximates the intractable true posterior density as ; and a probabilistic deep decoder network with parameters that can probabilistically reconstruct given z as . Both networks consist of three fully-connected layers as shown in Fig. 1.

To train the VAE, we generate consisting of configurations of heterogeneous tissue properties in a patient-specific cardiac mesh. The training involves optimizing the variational lower bound on the marginal likelihood of each training data with respect to network parameters and :


where we model

with a Bernoulli distribution. To optimize Eq. (


), stochastic gradient descent with standard backpropagation is utilized. Assuming that the approximate posterior

is a Gaussian density and the prior , their KL divergence can be derived analytically as:


where is the dimension of z, and and

are mean and variance from

. Because stochastic latent variables are utilized, the gradient of the expected negative reconstruction term during backpropagation cannot be directly obtained. The popular re-parameterization trick is utilized to express z as a deterministic variable as , where is noise [14, 10].

Probabilistic modeling of the latent code: After the training of VAE, in addition to the generative model provided by the decoder, the encoder provides an approximated conditional density of the LD latent code . This represents valuable knowledge about the probabilistic distribution of z learned from a large number of training data . To utilize this knowledge in the subsequent optimization, we integrate over the training data to obtain the density as a mixture of Gaussians , where and are mean and variance from . Because the number of mixture components in scales linearly with the number of training data, we obtain an approximation to by reducing the number of mixture components in two ways. First, we assume that can be represented with a single Gaussian density whose mean and variances is calculated as . Alternatively, we assume that can be represented by a mixture of Gaussians with components. To reduce the number of mixture components from to

, we use k-means clustering with the Bregman divergence 

[8] as a similarity metric on the -component Gaussian densities.

In this way, we obtain a generative model of HD tissue properties from an implicit LD manifold, and prior knowledge of the LD manifold from the probabilistic encoder. Both will be embedded into Bayesian optimization to enable efficient and accurate HD parameter estimation.

3.2 Bayesian Optimization with Embedded Generative Model

Representing the HD parameter with the expectation of the trained decoder , we embed the generative model into a revised objective function:


which allow us to optimize the HD parameter in an implicit LD manifold of z. For Bayesian optimization, we assume a Gaussian process (GP) – with a zero mean function and an anisotropic Mátern 5/2 co-variance function [3] – as a prior over the objective function ( 5). Then, the optimization consists of two iterative steps: first, by maximizing the acquisition function points in the LD manifold that allow the GP to both globally approximate Eq. (5) (exploration) and locally refine the regions of optimum (exploitation) is selected; second, the GP is re-trained with the recently selected point.

Acquisition function informed by VAE-encoded knowledge: To find optimal points in the LD manifold for updating the GP, we adopt the expected improvement (EI) function that selects point with maximum expected improvement over the current best objective function value [3]. For a GP posterior, it can be obtained analytically as:


where is the normal cumulative distribution, is the normal density function, and is the maximum of the objective function obtained so far. The first term here controls the exploitation (through high ) and the second term controls exploration (through high ). Because using only can lead to excessive local exploitation, a common practice is to augment with a constant trade-off parameter as:  [3]. Here, we utilize the VAE-encoded knowledge about the LD manifold

to enforce higher exploration in the regions of high probability density for

z, and lower elsewhere. In specific, we define , where , and are the weight, mean and variance of the Gaussian mixture components in .

GP Update: Once a new point is selected by maximizing EI, the objective function (2) is evaluated at the HD parameter given by the mean of obtained from the embedded generative model. The GP is then re-fitted to the updated data obtained by adding the new pair of and the corresponding value of the objective function by maximizing its log marginal likelihood with respect to kernel parameters: length scales and covariance amplitude.

4 Experiments

Figure 2: Comparison of BO-VAE EI post-1 (blue bar) with: 1) FH and FS (green bars); and 2) BO-VAE using standard EI, EI prior, and EI post-K (yellow bars) in terms of DC, RMSE, and number of model evaluations (from left to right).

Synthetic Experiments:

We include 27 synthetic experiments on three CT-derived human heart-torso models. In each case, an infarct sized of the heart was placed at differing locations using various combination of the AHA segments. The value of the parameter in the infarcted region and the healthy region is set to and , respectively. 120-lead ECG is simulated and corrupted with 20dB Gaussian noise as measurement data. To evaluate the accuracy in estimated parameters with two metrics: 1) root mean square error (RMSE) between the true and estimated parameters; and 2) dice coefficient (DC) = , where and are the sets of nodes in the true and estimated regions of infarct which is determined by using a threshold that minimizes the intra-region variance on the estimated parameter values [12].

VAE architecture and training: For each heart, we generate a training dataset corresponding to tissue properties with various heterogeneous infarcts. Each infarct is generated by random region growing in which, starting with one infarct node, one out of the five closest neighbors of the present set of infarct nodes is randomly added as an infarct node. This is repeated until an infarct of desired size is attained which is then added to the training dataset. Because infarcts generated in this fashion tend to be irregular in shape, we also add to the training data infarcts generated by growing the infarct using the closest neighbor to the infarct center. In total, we extracted 123,896, 155,099, and 116,459 training data for each heart. On these dataset, VAE with an architecture as shown in Fig. 1

that consists of an encoder and a decoder network with two hidden layers of 512 hidden units each, a pair of two-dimensional units for the mean and log-variance of the latent code, and sigmoid activation function is trained using

Adam optimizer. The training time for each dataset using Titan X Maxwell GPU was 9.77, 13.96, and 9.00 minutes, respectively.

Comparison with existing methods: The presented method (termed as BO-VAE) is compared against two common approaches based on explicitly defined LD representation of the cardiac mesh: 1) pre-fixed 18 segments (termed as fixed-segment (FS) method); and 2) partitions along a fixed multi-scale representation of the cardiac mesh obtained during coarse-to-fine optimization (termed as fixed-hierarchy (FH) method). Fig. 2

(a)-(b) summarizes the accuracy, demonstrating that BO-VAE (blue bar) is more accurate than the other two methods (green bars) in both DC and RMSE (paired t-tests, p

). Fig. 2(d) shows that this improvement in accuracy is achieved at a reduction of computation cost by for the FS method and for the FH method.

Figure 3: Left: Examples of estimated parameters with BO-VAE, FH, and FS. Right: Progression of FH on the multi-scale hierarchy for parameter estimation of case (a) (green leaf: homogeneous tissues and red leaf: heterogeneous tissues).

As expected, the FS method shows the lowest accuracy in which, as illustrated in Fig. 2, the estimated parameters either completely miss the infarct or are associated with a large region of false positives. This could partly be because direct optimization of 18 unknown model parameters without good initialization is difficult, while many of these dimensions are wasted at representing a region of homogeneous tissue. The FH method overcomes this issue to some extent, although it suffers from limited accuracy in many cases primarily due to an inflexible multi-scale hierarchy. An example is shown in the right panel of Fig. 3: several dimensions are wasted at representing homogeneous healthy regions (green nodes) distributed across different scales, which increasingly limits the ability of the optimization to go deeper along the tree to represent heterogeneous tissues at the necessary scale. In contrast, BO-VAE is not limited by such explicitly imposed structure of the parameter space, which allows it to attain higher accuracy with only 2 latent dimensions and a and fraction of the computation time of the FS and FH methods, respectively.
The effect of VAE-encoded knowledge about the LD manifold: To study the effect of incorporating the VAE-encoded knowledge of the LD manifold in the EI acquisition function, we compare the standard EI with EI augmented with three types of distributions on z: 1) (EI isotropic), 2) approximated with a single Gaussian density (EI Post-1); and 3) approximated with a mixture of Gaussian with 10 components (EI Post-K). As shown in Fig. 2, the estimation accuracy using all three distributions is higher than that without using any knowledge about z. In particular, the estimation accuracy with EI Post-1 is the highest. Fig. 4 illustrates how the knowledge from enables a more efficient search of the latent manifold. As shown, when is utilized, the exploration gradually proceeds from the region of high probability density to the region of low probability density (Fig. 4(b)). In comparison, without anything to guide the placement of the points, they are spread in an attempt to reduce overall variance (Fig. 4(a)). As a result, it could result in incorrect or suboptimal solution as shown in Fig. 4 (c) and (d), respectively.

We also experimented with utilizing a higher-dimensional latent code z. Fig. 5(c)(d) give examples of the estimated parameters using a five-dimensional (5d) vs. two-dimensional (2d) latent code, where only a marginal improvement of accuracy is observed with the increase in dimensions of z. This could be because, given the focus of the training data on tissue properties resulting from local infarcts, a 2d latent code is sufficient to account for the necessary latent generative factors. Fig. 5(a)-(b) shows the plot of these 2d latent codes corresponding to infarcts of different sizes and locations. Interestingly, it appears that the radial direction in the LD manifold accounts for the infarct size (a), while latent codes cluster by infarct location in the LD manifold (b).

Figure 4: Comparison of points selected by EI and EI post-1 shows that with EI post-1 the regions of higher is explored before the regions of lower .
Figure 5: LD manifold based on: (a) infarct size, and (b) infarct location shows that these information are encoded in the latent code. (c)(d) Examples in which estimated parameters with 5d latent code is more accurate than that with 2d latent code.

Real Data Experiments:

Real-data studies are conducted on two patients who underwent catheter ablation of ventricular tachycardia due to previous myocardial infraction. The patient-specific geometrical models of heart and torso are constructed from axial CT images. Using 120-lead ECG as measurement data, we evaluate the performance of the presented method in estimating local parameters in comparison with the FH and FS methods. The accuracy of the estimated parameters is evaluated using the bipolar voltage data from in-vivo catheter mapping. Note that since the voltage maps are not a direct measure of tissue properties, it is used as a reference rather than ground truth. The first two columns of Fig. 6 show the reference catheter mapping data (red: dense scar 0.5mV, purple: healthy tissue mV, and green: scar border mV) and the same data registered to CT-derived cardiac meshes.

Figure 6: Model parameter estimated with BO-VAE, FH, and FS on real-data study.

The catheter map in case 1 (Fig. 6(a)) shows a highly heterogeneous infarct spread over a large region in the lateral LV region. The estimated parameters by all three methods capture this region of infarct. To attain this accuracy, the FH and FS methods required 4056 and 1058 model evaluations, whereas BO-VAE required only 105 model evaluations. By contrast, as shown in Fig. 6(b), case 2 has a smaller region of dense scar in the lateral LV. The estimated model parameters by BO-VAE and FH correctly reveal this region of abnormal tissue, whereas although FS reveals a scar it has less accurate overlap with the low voltage region shown by the voltage map. To attain this level of accuracy, the presented method required only 105 model evaluations in comparison to the FH and FS methods that required 5798 and 1501 model evaluations, respectively.

5 Conclusion

We presented a novel approach for estimating HD cardiac model parameters by embedding within the Bayesian optimization framework a generative model of the HD tissue properties from an implicit LD manifold. Experiments show a gain in accuracy with drastically reduced computational cost. Future works include two direction: 1) to incorporate more realistic training data from high resolution 3D imaging for a more expressive generative model and potentially an improved accuracy in estimating highly heterogeneous tissues; and 2) to improve the efficiency by investigating novel ways to incorporate the knowledge of latent manifold to guide the active selection of points during BO.


  • [1] Aliev, R.R., Panfilov, A.V.: A simple two-variable model of cardiac excitation. Chaos, Solitons & Fractals 7(3), 293–301 (1996)
  • [2] Arevalo, H.J., Vadakkumpadan, F., Guallar, E., et al.: Arrhythmia risk stratification of patients after myocardial infarction using personalized heart models. Nature communications 7 (2016)
  • [3] Brochu, E., Cora, V.M., De Freitas, N.: A tutorial on bayesian optimization of expensive cost functions, with application to active user modeling and hierarchical reinforcement learning. arXiv preprint arXiv:1012.2599 (2010)
  • [4] Chinchapatnam, P., Rhode, K.S., Ginks, M., et al.: Estimation of volumetric myocardial apparent conductivity from endocardial electro-anatomical mapping. In: Engineering in Medicine and Biology Society, 2009. EMBC 2009. Annual International Conference of the IEEE. pp. 2907–2910. IEEE (2009)
  • [5] Dhamala, J., Arevalo, H., Sapp, J., Horacek, M., Wu, K., Trayanova, N., Wang, L.: Spatially-adaptive multi-scale optimization for local parameter estimation in cardiac electrophysiology. IEEE TMI (2017)
  • [6] Dhamala, J., Sapp, J., Horacek, M., Wang, L.: Spatially-adaptive multi-scale optimization for local parameter estimation: application in cardiac electrophysiological models. In: MICCAI. pp. 282–290. Springer (2016)
  • [7]

    Dhamala, J., Sapp, J., Horacek, M., Wang, L.: Quantifying the uncertainty in model parameters using gaussian process-based markov chain monte carlo: An application to cardiac electrophysiological models. In: IPMI. pp. 223–235. Springer (2017)

  • [8]

    Garcia, V., Nielsen, F., Nock, R.: Levels of details for gaussian mixture models. In: Asian Conference on Computer Vision. pp. 514–525. Springer (2009)

  • [9] Giffard-Roisin, S., Jackson, T., Fovargue, L., et al.: Noninvasive personalization of a cardiac electrophysiology model from body surface potential mapping. IEEE Trans. on Biomed. Eng. 64(9), 2206–2218 (2017)
  • [10] Kingma, D.P., Welling, M.: Auto-encoding variational bayes. arXiv preprint arXiv:1312.6114 (2013)
  • [11] Lê, M., Delingette, H., Kalpathy-Cramer, J., et al.: Mri based bayesian personalization of a tumor growth model. IEEE TMI 35(10), 2329–2339 (2016)
  • [12] Otsu, N.: A threshold selection method from gray-level histograms. Automatica 11(285-296), 23–27 (1975)
  • [13] Plonsey, R.: Bioelectric phenomena. Wiley Online Library (1969)
  • [14] Rezende, D.J., Mohamed, S., et al.: Stochastic backpropagation and approximate inference in deep generative models. arXiv preprint arXiv:1401.4082 (2014)
  • [15] Sermesant, M., Chabiniok, R., Chinchapatnam, P., Mansi, et al.: Patient-specific electromechanical models of the heart for the prediction of pacing acute effects in crt: a preliminary clinical validation. Medical image analysis 16(1), 201–215 (2012)
  • [16] Wang, L., Zhang, H., Wong, K.C., Liu, H., Shi, P.: Physiological-model-constrained noninvasive reconstruction of volumetric myocardial transmembrane potentials. IEEE Trans. on Biomed. Eng. 57(2), 296–315 (2010)
  • [17] Wong, K.C., Sermesant, M., Rhode, K., et al.: Velocity-based cardiac contractility personalization from images using derivative-free optimization. Journal of the mechanical behavior of biomedical materials 43, 35–52 (2015)
  • [18] Zettinig, O., Mansi, et al.: Fast data-driven calibration of a cardiac electrophysiology model from images and ecg. In: MICCAI. pp. 1–8. Springer (2013)