1 Introduction
Patientspecific simulation models of the heart have shown increasing potential in personalized treatment of cardiac diseases [2, 15]. The estimation of patientspecific 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 highdimensional (HD) model parameters is not only algorithmically difficult given indirect and sparse measurements, but also computationally intractable in the presence of computingintensive 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 predivided into 326 segments, each represented by a uniform parameter value [17, 18]. Naturally, this artificial lowresolution 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 coarsetofine optimization along a predefined multiscale hierarchy of the cardiac mesh, enabling spatiallyadaptive 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 predefined multiscale 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 coarsetofine 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 lowdimensional (LD) or multiscale representation of the parameter space with an implicit LD latent encoding, achieved by embedding within the optimization a stochastic LDtoHD generative model that describes the generation of the HD spatiallyvarying tissue properties from a LD manifold. This generative model is obtained with a variational autoencoder (VAE), trained from a large set of spatiallyvarying 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 patientspecific model parameters.
The presented method is applied to the estimation of local tissue excitability of a cardiac electrophysiological model using noninvasive electrocardiogram (ECG) data. On both synthetic and real data experiments, the presented method is compared against existing methods based on explicitlydefined LD [17] or multiscale 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 AlievPanfilov (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:
(1) 
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 .Bodysurface ECG Model: The propagation of the spatiotemporal transmural action potential U to the potentials measured on the body surface Y can be described by the quasistatic approximation of the electromagnetic theory [13]. Solving the governing equations on a discrete hearttorso 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 :
(2) 
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 spatiallyvarying tissue properties at the resolution of the cardiac mesh, and 2) a novel Bayesian optimization method utilizing the embedded generative model.
3.1 LDtoHD 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 fullyconnected layers as shown in Fig. 1.To train the VAE, we generate consisting of configurations of heterogeneous tissue properties in a patientspecific cardiac mesh. The training involves optimizing the variational lower bound on the marginal likelihood of each training data with respect to network parameters and :
(3) 
where we model
with a Bernoulli distribution. To optimize Eq. (
3), 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:(4) 
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 reparameterization 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 kmeans 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:
(5) 
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 covariance 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 retrained with the recently selected point.
Acquisition function informed by VAEencoded 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:
(6) 
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 tradeoff parameter as: [3]. Here, we utilize the VAEencoded 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 refitted 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
Synthetic Experiments:
We include 27 synthetic experiments on three CTderived human hearttorso 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. 120lead 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 intraregion 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 twodimensional units for the mean and logvariance 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 BOVAE) is compared against two common approaches based on explicitly defined LD representation of the cardiac mesh: 1) prefixed 18 segments (termed as fixedsegment (FS) method); and 2) partitions along a fixed multiscale representation of the cardiac mesh obtained during coarsetofine optimization (termed as fixedhierarchy (FH) method). Fig. 2
(a)(b) summarizes the accuracy, demonstrating that BOVAE (blue bar) is more accurate than the other two methods (green bars) in both DC and RMSE (paired ttests, 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.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 multiscale 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, BOVAE 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 VAEencoded knowledge about the LD manifold:
To study the effect of incorporating the VAEencoded 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 Post1); and 3) approximated with a mixture of Gaussian with 10 components (EI PostK). 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 Post1 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 higherdimensional latent code z. Fig. 5(c)(d) give examples of the estimated parameters using a fivedimensional (5d) vs. twodimensional (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).
Real Data Experiments:
Realdata studies are conducted on two patients who underwent catheter ablation of ventricular tachycardia due to previous myocardial infraction. The patientspecific geometrical models of heart and torso are constructed from axial CT images. Using 120lead 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 invivo 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 CTderived cardiac meshes.
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 BOVAE 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 BOVAE 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.
References
 [1] Aliev, R.R., Panfilov, A.V.: A simple twovariable 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 electroanatomical 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.: Spatiallyadaptive multiscale optimization for local parameter estimation in cardiac electrophysiology. IEEE TMI (2017)
 [6] Dhamala, J., Sapp, J., Horacek, M., Wang, L.: Spatiallyadaptive multiscale 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 processbased 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] GiffardRoisin, 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.: Autoencoding variational bayes. arXiv preprint arXiv:1312.6114 (2013)
 [11] Lê, M., Delingette, H., KalpathyCramer, 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 graylevel histograms. Automatica 11(285296), 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.: Patientspecific 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.: Physiologicalmodelconstrained 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.: Velocitybased cardiac contractility personalization from images using derivativefree optimization. Journal of the mechanical behavior of biomedical materials 43, 35–52 (2015)
 [18] Zettinig, O., Mansi, et al.: Fast datadriven calibration of a cardiac electrophysiology model from images and ecg. In: MICCAI. pp. 1–8. Springer (2013)
Comments
There are no comments yet.