Personalized computer models of the heart have shown promise in various clinical tasks such as risk stratification [2, 13] and predicting treatment response . With advances in medical imaging technologies, personalization of high-resolution anatomical models is now feasible. By contrast, organ tissue properties vary spatially over the three-dimensional anatomical domain but have to be estimated from indirect, sparse and noisy data, resulting in a difficult ill-posed optimization problem with a high-dimensional (HD) unknown.
Numerous works have been presented in estimating these spatially-varying tissue properties in the form of high-dimensional model parameters. Most existing methods choose to represent spatially-varying tissue properties via a low-dimensional (LD) partitioning of the underlying geometrical model, either as pre-defined segments , or iteratively optimized in a coarse-to-fine fashion [5, 14]. This LD-to-HD definition directly exploits the spatial proximity and hierarchical composition of the underlying geometry. However, it is of such limited expressiveness that the number of partitioning is either too low to faithfully represent high-resolution tissue properties, or too high to allow effective optimization. In contrast, recent work presented the use of a data-driven generative model of HD tissue properties, via a variational auto-encoder (VAE), to embed the optimization into a LD latent space . Being more expressive, this VAE-based generative model is able to represent high-resolution tissue properties with a latent code sufficiently small for effective optimization. However, as the regular VAE is defined over Euclidean data, it does not take into account the valuable geometry information in the data, nor does it allow transferring among different geometry without first establishing point-by-point correspondence.
If we view organ tissue properties over a 3D geometrical model as an image, convolutional neural networks (CNNs) are a natural choice to incorporate knowledge of the spatial proximity and hierarchical composition of the image. However, standard CNNs have been most successful on data with an underlying Euclidean structure (i.e., image grids). Generalizing CNNs to non-Euclidean domains is an emerging area of research , where significant efforts have been presented on addressing the challenges of defining convolution , pooling, and up-sampling operations . However, most developments to non-Euclidean CNNs are focused on supervised discriminative networks. To date, very limited work have been presented to enable generative modeling of non-Euclidean data.
In this paper, we present a novel VAE architecture that allows generative modeling of data over non-Euclidean domains, and utilize this generative model to embed Bayesian optimization of large graphs into a LD manifold. The presented approach bridges the gap of previous works by introducing an expressive generative model that is able to represent high-resolution tissue properties with a small latent code, while incorporating the geometrical knowledge in the data and being transferable across geometries. We evaluate the presented method in synthetic and real-data experiments of estimating tissue excitability in a cardiac electrophysiological model, where we compare the expressiveness of the generative model and the accuracy of the subsequent parameter optimization with those obtained by using a linear reconstruction model based on principal component analysis (PCA) and a regular fully-connected VAE. We further demonstrate the feasibility of transferring the presented non-Euclidean VAE across patients. To our knowledge, this is the first introduction of a graph convolutional VAE and its use to enable Bayesian optimization over large graphs.
2 Background: Models of Cardiac Electrophysiology
Cardiac Electrophysiological Model: Among various computational models of cardiac electrophysiology, phenomenological models are widely used in parameter personalization as they are able to express key macroscopic properties of cardiac excitation with a small number of parameters [14, 5]. We thus adopt the phenomenological two-variable Aliev-Panfilov (AP) model :
where is the normalized transmembrane potential, and is the recovery current. Parameter controls the coupling between and ,
is the diffusion tensor,controls the repolarization, and controls tissue excitability. As is most sensitive to parameter in the AP model (1) , we focus on its estimation and use standard literature values for the remaining parameters . Solving the AP model on a 3D discrete cardiac anatomy of meshfree nodes , we obtain a 3D electrophysiological model of the heart that describes the temporal evoluation of 3D transmembrane potential .
Measurement Model: is measured on the body surface following the quasi-static electromagnetic theory, solving which on a discrete subject-specific heart-thorax mesh gives a linear relationship between and its surface potential measurement as: .
3 Personalizing HD Parameters on Unstructured Meshes
We seek parameter that minimizes the sum of squared errors between model output and patient’s measurements as:
Directly solving (2) for the spatially-distributed is difficult . Below we describe how we learn a LD-to-HD generation of that accounts for the underlying geometry, and embed the HD optimization into the expressive LD manifold.
3.1 Graph Convolutional VAE
model the generation of spatially-distributed with a VAE . A VAE consists a probabilistic encoder network with parameters that approximates the intractable true posterior density as ; and a probabilistic decoder network with parameters that represents the likelihood as . For data defined on a Euclidean grid, structural information is incorporated in VAE through CNNs. Here, data resides over the 3D heart geometry, on which standard convolution and pooling operations are not applicable. Below we present a novel VAE architecture that enables convolution, pooling, and unpooling over non-Euclidean geometry of the heart.
Local Connectivity & Graph Convolution: We model the cardiac mesh as a graph: , where vertices consist of all meshfree nodes and edges exist between each meshfree node and its nearest neighbors. consists of edge attributes , calculated as the normalized if an edge exists between vertices at and and 0 otherwise. On this graph, we use a convolution operator based on spatial continuous convolution kernels because it was shown to allow better generalization to similar graphs . In specific, given the graph and -dimensional input features , the -th convolution kernel is:
where denotes open B-spline bases of degree
based on equidistant knot vectors with-dimensional kernel size of , is the Cartesian product of the B-spline bases, and are the trainable parameters. Given kernel functions and input features , the spatial convolution operator for each vertex with a neighborhood based on its edge connectivity is then defined as:
Hierarchical Composition & Pooling: To define pooling and unpooling operations necessary for the encoding-decoding architecture, a hierarchical representation of the graph is needed. We obtain this by an efficient multilevel graph clustering method based on minimizing the normalized cuts (Graclus) , which reduces the graph size by half the number of vertices at each coarsening.
We store hierarchical graph representation in matrices which reduces pooling/unpooling operations to efficient matrix multiplications. In specific, if is a graph with vertices and is its coarsened graph with vertices, we populate a binary matrix , where if the vertex in was grouped to the vertex in and otherwise. Given feature maps over the vertices of graph and over graph , the average pooling in the encoder can be obtained by and unpooling in the decoder by , where is column normalized from .
Graph Convolutional VAE: Using these building blocks along with global pooling and standard convolution layers, we construct a VAE architecture as shown in Fig. 1. It is trained by optimizing the variational lower bound on the marginal likelihood of the training data : . We set and to be Gaussian parameterized by the graph convolutional networks. The prior is set to be an isotropic Gaussian, producing an analytical form for the KL divergence. Using the re-parameterization trick , to express the random as a deterministic variable, standard stochastic gradient methods can be used to optimize .
3.2 Bayesian Optimization on Large Graphs
Bayesian optimization is a popular choice in optimizing complex objective functions such as (2) . It begins by defining a surrogate over the objective function. The optimization then consists of two iterative steps: 1) actively find a point that optimizes a utility function based on the surrogate, and 2) update the surrogate with the newly-selected point. Direct Bayesian optimization over HD space is difficult and its use over large graphs has not been reported [3, 9]. To enable this, we reformulate the original objective function in (2) by replacing the unknown with the expectation of its non-Euclidean generative model:
This allows us to embed surrogate construction and active selection of training points in a LD manifold. We initialize the Gaussian process (GP) surrogate of (5) with a zero mean function and an anisotropic Mátern 5/2 kernel.
Active Selection of Training Points: To select a training point, we maximize the expected improvement (EI) utility function that favours a point with the highest expected improvement over the current optimum [3, 15]:
are the predictive mean and standard deviation of the GP,
is the cumulative normal distribution, andis the normal density function. The first term promotes exploitation by favoring regions of high , while the second term promotes exploration by favoring the regions with high .
4 Synthetic experiments
We evaluate the presented graph convolutional VAE (termed as gVAE) by: 1) its reconstruction accuracy, and 2) optimization accuracy of the gVAE-based Bayesian optimization, both in comparison to existing methods. Accuracy is evaluated by the sum of squared error (SSE) in , and the dice coefficient (DC) of the abnormal region obtained by thresholding with Otsu’s method .
We synthetically generate data of heterogeneous tissue excitability via random region growing in each cardiac model. Beginning with a single meshfree node as abnormal, we randomly grow the abnormal region by adding one of the nearest neighbors of the abnormal nodes, until the abnormal region reaches a desired size ( to of the heart volume). On average, we generated data for training and data each for validation and testing. The various layers and sizes of feature maps in each layer of the presented gVAE architecture are detailed in Fig. 1. We use B-spline basis degree of with kernel size of in all graph convolution layers. All models are trained with a learning rate of with Adam optimizer .
gVAE as a Generative Model: On five different human heart models constructed from CT images, we first evaluate the ability of the presented gVAE model to reconstruct tissue excitability in comparison to: 1) PCA-based linear reconstruction, and 2) fully-connected VAE (termed as fVAE)  with three to five hidden layers (termed as fVAE-3h, -4h, and -5h, respectively). As summarized in Table 1, PCA being a linear model has the lowest accuracy. Compared to fVAE, the reconstruction accuracy of gVAE is consistently higher in both DC and SSE, even when its number of trainable parameters is similar to or lower than fVAE. However, gVAE is more expensive to train: 37.21 hrs vs. 7.51 mins for fVAE-3h in TITAN Xp GPU. Nevertheless, note that the number of trainable parameters for gVAE does not increase for larger meshes.
We further compare gVAE with two-dimensional (2d) latent codes vs. fVAE and PCA with various latent dimensions. As shown in Fig. 2(a), to achieve the reconstruction accuracy of gVAE with 2d latent codes, at least 13 principles components are required with PCA: this increase in dimension will make the subsequent optimization difficult. With fVAE, a similar SSE is attained with four latent dimensions, while a similar DC could not be attained with even 50 latent dimensions. This may be because fVAE does not consider the geometry underlying the spatial distribution of tissue excitability. Note that this increase in the number of latent dimensions will increase both the difficulty and the computational cost of the subsequent optimization. Fig. 2(b) and (c) show that the latent code learned with gVAE are clustered by the location of the abnormal tissue and its radial direction encodes the size of the abnormal tissue.
gVAE-based Parameter Optimization: On 40 synthetic cases on two different heart geometries, we conduct experiments on estimating unknown tissue excitability. In each case, we set an abnormal region by using various combinations of AHA segments which are very different from the training set. Measurement data was simulated, with models described in Section 2, and then sub-sampled and corrupted with 20dB Gaussian noise. We compare the accuracy of gVAE-based Bayesian optimization with three previous approaches: 1) optimization on 17 fixed segments (termed as FS) [17, 14], 2) coarse-to-fine optimization along a fixed multi-scale mesh hierarchy (termed as FH) , and 3) optimization on a LD manifold obtained with fVAE-3h .
Result in Fig. 3
(a) shows that gVAE-based Bayesian optimization is more accurate than all other approaches in both DC and SSE (paired t-test,). This indicates that gVAE was better at capturing the LD generative factors necessary for accurate optimization. Fig. 3(c)-(f) shows visual comparison on a few examples. The computational cost of gVAE-based optimization is much lower than that of FH (x) and FS (x), and similar to that of fVAE-based optimization.
Feature Sharing across Geometries: To demonstrate the feasibility of transferring the presented gVAE across different geometries, we take a pre-trained gVAE, fix the learned features in the encoder’s graph convolution layers, and fine-tune the remaining layers for a different anatomy. We compare this training strategy to training a gVAE from scratch. Results in Fig. 4(a) show that a pre-trained model can be fined-tuned with as small as 1088 new examples. In comparison, gVAE could not be trained from scratch with samples, shown both by the low reconstruction accuracy (Fig. 4(a): DC = 0.24; SSE = 17.68) and a flat test loss plot (Fig. 4(b)). Test loss plots in Fig. 4(b) also show that a pre-trained model starts with a lower loss and a larger size of training data leads to faster convergence. Parameter optimization via a gVAE fine-tuned with 7360 data on 20 cases achieved an average DC and SSE of 53.10 and 11.01 respectively. Fig. 4(c) shows some examples of the estimated parameters. More detailed analysis and experimentation is left to future work.
5 Real Data Experiments
We conduct real-data studies on two patients who underwent catheter ablation of ventricular tachycardia due to chronic myocardial infarction. Patient-specific heart-thorax models are obtained from axial CT images. Using 120-lead ECG as measurements, we evaluate the performance of the presented gVAE in estimating tissue excitability in comparison to the fVAE , FH , and FS  methods. Training dataset and network architectures are as described in Section 4. We qualitatively evaluate the results with in-vivo catheter mapping data which, as shown in Fig. 5, provides a reference for the location of the abnormal (red, voltage ) and healthy (purple, voltage ) regions.
Case 1 has a large heterogeneous abnormal region in the lateral LV region (Fig. 5(a)). All methods are able to localize this region, but with varying degree of heterogeneity. Optimizations based on gVAE and fVAE are much faster requiring only 100 model evaluations, in comparison to FH and FS that required 4056 and 1058 model evaluations, respectively. By contrast, case 2 has a smaller but dense abnormal region in the lateral LV (Fig. 5(b)). While all methods identify the general location of this abnormality, gVAE more accurately differentiates the region of dense core and border. In comparison, fVAE and FH estimate a larger border region; and the abnormal region revealed by FS is less accurate. Again, in this case, gVAE and fVAE required only 100 model evaluations, whereas FH and FS required 5798 and 1501 model evaluations, respectively.
We presented a novel graph convolutional VAE model that allows generative modeling of data defined over non-Euclidean structures and integrated it with the Bayesian optimization to enable optimization on large graphs. Experiments on estimating tissue properties distributed on a 3D cardiac mesh shows higher accuracy in terms of both reconstructing the tissue excitability and estimating them from indirect measurements, both in comparison to existing baselines. In future, we will incorporate realistic data from high resolution 3D images and investigate efficient transfer learning methods for models trained on such data.
-  Aliev, R.R., Panfilov, A.V.: A simple two-variable model of cardiac excitation. Chaos, Solitons & Fractals 7(3), 293–301 (1996)
-  Arevalo, H.J., Vadakkumpadan, F., Guallar, E., Jebb, A., Malamas, P., Wu, K.C., Trayanova, N.A.: Arrhythmia risk stratification of patients after myocardial infarction using personalized heart models. Nature communications 7 (2016)
-  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)
Bronstein, M.M., Bruna, J., LeCun, Y., et al.: Geometric deep learning: going beyond euclidean data. IEEE Signal Processing Magazine34(4), 18–42 (2017)
-  Dhamala, J., Arevalo, H., Sapp, J., et al.: Spatially-adaptive multi-scale optimization for local parameter estimation in cardiac electrophysiology. IEEE TMI (2017)
-  Dhamala, J., Ghimire, S., Sapp, J.L., Horáček, B.M., Wang, L.: High-dimensional bayesian optimization of personalized cardiac model parameters via an embedded generative model. In: MICCAI. pp. 499–507. Springer (2018)
Dhillon, I.S., Guan, Y., Kulis, B.: Weighted graph cuts without eigenvectors a multilevel approach. IEEE TPAMI29(11), 1944–1957 (2007)
-  Fey, M., Eric Lenssen, J., Weichert, F., Müller, H.: Splinecnn: Fast geometric deep learning with continuous b-spline kernels. In: CVPR. pp. 869–877 (2018)
Kandasamy, K., Schneider, J., Póczos, B.: High dimensional bayesian optimisation and bandits via additive models. In: International Conference on Machine Learning. pp. 295–304 (2015)
-  Kingma, D.P., Ba, J.: Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980 (2014)
-  Kingma, D.P., Welling, M.: Auto-encoding variational bayes. arXiv preprint arXiv:1312.6114 (2013)
-  Otsu, N.: A threshold selection method from gray-level histograms. Automatica 11(285-296), 23–27 (1975)
-  Prakosa, A., Arevalo, H.J., Deng, D., Boyle, P.M., Nikolov, P.P., Ashikaga, H., et al.: Personalized virtual-heart technology for guiding the ablation of infarct-related ventricular tachycardia. Nature Biomedical Engineering (2018)
-  Sermesant, M., Chabiniok, R., Chinchapatnam, P., 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)
-  Shahriari, B., Swersky, K., Wang, Z., Adams, R.P., De Freitas, N.: Taking the human out of the loop: A review of bayesian optimization. Proceedings of the IEEE 104(1), 148–175 (2016)
-  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)
-  Wong, K.C., Relan, J., Wang, L., Sermesant, M., Delingette, H., Ayache, N., Shi, P.: Strain-based regional nonlinear cardiac material properties estimation from medical images. In: MICCAI. pp. 617–624. Springer (2012)
-  Ying, Z., You, J., Morris, C., Ren, X., Leskovec, J.: Hierarchical graph representation learning with differentiable pooling. In: NeurIPS. pp. 4805–4815 (2018)