^{1}

^{2}

^{3}

^{4}

^{4}

^{1}

^{2}

^{3}

^{*}

^{5}

^{1}

^{2}

^{3}

^{4}

^{5}

Edited by: Valeriya Naumova, Simula Research Laboratory, Norway

Reviewed by: Milan Žukovič, University of Pavol Jozef Šafárik, Slovakia; Joakim Sundnes, Simula Research Laboratory, Norway

This article was submitted to Computational Physics, a section of the journal Frontiers in Physics

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

A critical procedure in diagnosing atrial fibrillation is the creation of electro-anatomic activation maps. Current methods generate these mappings from interpolation using a few sparse data points recorded inside the atria; they neither include prior knowledge of the underlying physics nor uncertainty of these recordings. Here we propose a physics-informed neural network for cardiac activation mapping that accounts for the underlying wave propagation dynamics and we quantify the epistemic uncertainty associated with these predictions. These uncertainty estimates not only allow us to quantify the predictive error of the neural network, but also help to reduce it by judiciously selecting new informative measurement locations via active learning. We illustrate the potential of our approach using a synthetic benchmark problem and a personalized electrophysiology model of the left atrium. We show that our new method outperforms linear interpolation and Gaussian process regression for the benchmark problem and linear interpolation at clinical densities for the left atrium. In both cases, the active learning algorithm achieves lower error levels than random allocation. Our findings open the door toward physics-based electro-anatomic mapping with the ultimate goals to reduce procedural time and improve diagnostic predictability for patients affected by atrial fibrillation. Open source code is available at

Atrial fibrillation is the most common arrhythmia in the heart, affecting between 2.7 and 6.1 million people in the United States alone [

Deep learning has revolutionized many fields in engineering and medical sciences. However, in the context of activation mapping, we have to rely on a few sparse activation measurements. To address the limitations of deep learning associated with sparse data, recent techniques have emerged to incorporate the underlying physics, as governed by partial differential equations, into neural networks [

Here, we propose to use a physics-informed neural network to create activation maps of the atria [

This manuscript is organized as follows: In section 2 we introduce the physics-informed neural network, the method to estimate uncertainties and the active learning algorithm. In section 3, we show two numerical experiments to test the accuracy of the method. We end this manuscript with a discussion and future directions in section 4.

In this section we introduce a physics-informed neural network [

The electrical activation map of the heart can be related to a traveling wave, where the wavefront represents the location of cells that are depolarizing [

where

Further, we approximate both the activation time

where _{T} and _{V} are neural networks with parameters

where _{max} represents the maximum conduction velocity, specified by the user. Finally, we define a loss function to train our model:

_{T} activation time measurements available _{R} collocation points. The third and fourth terms serve as regularization for the inverse problem. The third term, which we evaluate at the _{R} collocation points, is a total variation regularization for the conduction velocity, which allows discrete jumps in the solution. We select this term to model slow regions of conduction, such as fibrotic patches. Finally, we use _{2} regularization on the weights of the activation time neural network, for reasons we explain in the following section. We solve the following minimization problem to train the neural networks and find the optimal parameters:

Physics-informed neural networks for activation mapping. We use two neural networks to approximate the activation time

We will be interested in quantifying the uncertainty in our predictions to both inform physicians about the quality of the estimates as well as to use active learning techniques to judiciously acquire new measurements. Given the large number of parameters in neural networks, using gold-standard methods for Bayesian inference, such as Markov Chain Monte Carlo, is prohibitively expensive. Instead, we borrow ideas from deep reinforcement learning and use randomized prior functions [

We draw the parameters

We take advantage of the uncertainty estimates described in the previous section to create an adaptive sampling strategy. We start with a small number of randomly located samples, fit our model, and then acquire the next measurement at the point that we estimate to reduce the predictive model error the most. We iteratively fit the model and acquire samples until we reach a user-defined convergence or until we exceed our budget or time to acquire new data. Since the exact predictive error is unknown, there are several heuristics to determine where to place the next sample. A common approach is to select the location where the uncertainty is the highest, which we can quantify by evaluating the entropy of the posterior distribution

Active learning algorithm to iteratively identify the most efficient sampling points

During electro-anatomic mapping, data can only be acquired on the cardiac surface, either of the ventricles or the atria. We thus represent the resulting map as a surface in three dimensions and neglect the thickness of the atrial wall. This is a reasonable assumption since the thickness-to-diameter ratio of the atria is in the order of 0.05. Our assumption implies that the electrical wave can only travel along the surface and not perpendicular to it. To account for this constraint, we include an additional loss term:

This form favors solutions where the activation time gradients are orthogonal to the surface normal _{R} collocation points as the centroids of each triangle in the mesh _{N}. If the gradient of the activation times were exactly orthogonal to the triangle normals, it would force a linear interpolation between mesh nodes in the neural network, which is not desirable and unlikely attainable.

We implement all models in Tensorflow [_{R} collocation points, we use a minibatch implementation, in which we use a subset of all available collocation points to compute the loss function and its gradient. For the two-dimensional benchmark problem, we randomly sample _{mb} points using a Latin hypercube design [_{mb}. For each iteration, we use the centroid locations of the triangles of one of these batches as collocation points and loop through them as the optimization progresses.

In this section, we explore our method using a two-dimensional benchmark problem and a three-dimensional personalized left atrium. We also quantify the effectiveness of the active learning algorithm.

To characterize the performance of the proposed model, we design a synthetic benchmark problem that analytically satisfies the Eikonal equation. We introduce a discontinuity in the conduction velocity and collision of two wavefronts in the following form:

with _{TV} and α_{L2} and then select them to _{mb} = 100 and then train with the L-BFGS method [

Benchmark problem activation times and conduction velocities. The top row shows the activation times, the bottom row the conduction velocity. The columns display the exact solution and the results of the physics-informed neural network (PINN), a neural network without physics, the Gaussian process regression, and the linear interpolation. The black circles indicate the sampling locations.

We compare our method against three other methodologies: a neural network with the same architecture and parameters except without including the physics, linear interpolation [

Benchmark problem x = y line. The dashed lines show the exact analytical solution, the solid lines the different interpolation methods. The physics-informed neural network (PINN) captures the gradient of the activation times produced by the collision of two wavefronts and closely predicts the conduction velocity. The neural network and Gaussian process regression produce a smooth interpolation that generates oscillating conduction velocities, especially near the collision of the two wavefronts. The linear interpolation produces an approximation that is highly sensitive to the position of the data points.

Normalized errors for benchmark problem.

0% | T | 0.49 | 0.48 | 1.56 | 2.05 |

V | 3.56 | 9.07 | 57.78 | 17.98 | |

1% | T | 1.91 (1.09–3.09) | 4.19 (2.74–7.93) | 1.92 (1.48–2.57) | 2.23 (1.97–2.56) |

V | 10.24 (6.09–17.86) | 63.93 (42.74–91.12) | 54.08 (47.82–67.72) | 27.37 (21.49–40.13) | |

5% | T | 3.42 (2.22–5.34) | 11.00 (6.90–18.42) | 3.34 (2.63–4.13) | 4.44 (3.58–5.76) |

V | 15.84 (10.40–23.20) | 87.33 (68.19–128.83) | 78.87 (53.77–102.81) | 66.30 (42.97–171.79) | |

10% | T | 6.70 (3.75–10.60) | 18.25 (12.84–33.39) | 6.73 (3.54–11.65) | 8.55 (6.02–11.55) |

V | 23.16 (10.81–40.78) | 90.09 (78.46–119.42) | 96.78 (50.49–177.69) | 81.06 (59.86–241.80) |

To conclude, we evaluate the performance of all four methods to noise. We introduce Gaussian noise with a standard deviation of 1, 5, and 10% of the maximum value of activation time and run all methods 30 times with the same datasets. In this case, we train the physics-informed neural network and the regular neural network for 100,000 ADAM iterations.

To study the accuracy of our uncertainty estimates, we vary the number of neural networks trained in our randomized priors ensemble. We set the noise parameter to σ_{N} = 0.01, based on reported uncertainties [_{mb} = 96. We define the entropy estimated by 100 neural networks as our ground truth and compute the root mean square error for the remaining cases. We also compute the time it takes to train the networks.

Uncertainty quantification and effect of number of neural networks. Entropy predicted by 10, 30, and 100 neural networks, and the tradeoff between entropy reduction and normalized training time for increasing number of neural networks. The black circles indicate the sampling locations.

We also test the hypothesis that the maximum entropy is co-localized with the predictive error of the model. If this is true, the active learning approach will be effective, since placing a new sample at the point of maximum entropy will reduce the error. We run the two-dimensional example with the same parameters specified in the previous sections with 30 samples drawn from a Latin hypercube design and train 30 neural networks in parallel.

Correlation of uncertainty and error. For the benchmark problem, trained with 30 samples, the computed entropy tends to be higher at regions where the error is higher and the points of maximum entropy and maximum error (⋆) are co-located. The black circles indicate the sampling locations.

To test the effectiveness of the proposed active learning method, we train models with 30 different initial conditions. We start with _{init} = 10 samples and follow algorithm 1, acquiring _{AL} = 40 additional samples. We also train 30 models where ^{−7}) and in conduction velocity (

Benchmark problem and active learning. We perform 30 simulations of active learning with different initial samples and compare them against a Latin hypercube design. In the middle and right box plots, we observe a significant reduction in activation time normalized root mean squared error (^{−7}) and in conduction velocity normalized mean absolute error (

To test our model in three dimensions, we studying the electrophysiology of the left atrium. We obtain the mesh from one of the examples of the MICCAI challenge dataset [^{2}/ms in the entire domain and one where it is heterogeneous such that half of the domain has a reduced conductivity of 0.05 mm^{2}/ms. In both cases, we initiate the activation at the center of the septum. We define the activation time as the time at which the transmembrane potential reaches 0 mV. We then compute the conduction velocity as _{max} = 1 m/s. We consider two experiments.

In the first experiment, we set the number of samples by the optimal density [^{2}. We randomly choose nodes in the mesh and use the activation times as data points. We repeat this process 30 times and compare the error to using a linear interpolation trained with the same data as detailed in section 3.1 with a Wilcoxon test [^{−5}) for both the homogeneous constant conductivity case (median 1.53, range 1.34–1.85 ms vs. median 3.92, range 3.42–5.34 ms) and the heterogeneous case with reduced conductivity in half of the domain (median 2.23, range 1.84–2.84 ms vs. median 4.77, range 3.90-6.02 ms).

Left atrium accuracy in activation times and conduction velocities. Performance of the physics-neural network (PINN) vs. linear interpolation, top left. The physics-neural network significantly reduces the root mean squared error (RMSE) in the homogeneous case with constant conductivity and in the heterogeneous case where the conductivity was reduced in half of the domain (^{−5}). Performance of the active learning algorithm for these two cases, top center and top right. Active learning decreases the root mean squared error of the activation times and the mean absolute error of the conduction velocity. The dashed line represents the median error achieved by linear interpolation at optimal sample density. Sample density required by the active learning approach to achieve the median error of the linear interpolation method at optimal sample density, bottom left. Comparison of active learning against random selection in the case of constant conductivity (bottom center) and reduced conductivity (bottom right). The active learning algorithm significantly reduces the error in all cases (

In the second experiment, we explore the performance of the active learning algorithm by running 30 cases that start with _{init} = 10 randomly selected samples, which corresponds to a density of 0.081 samples/cm^{2}. We acquire samples until we reach 90 measurements, which corresponds to a density of 0.734 samples/cm^{2}. ^{2} for the homogeneous case (1.01–1.88 ms) and the heterogeneous case (1.38–2.28 ms). The mean absolute errors in conduction velocity are relatively high for both cases (constant: median 0.086, range 0.083–0.101 m/s; half: median 0.085, range 0.081–0.095 ms). This can be explained, at least in part, by the difficulty of computing conduction velocities from the cardiac electrophysiology simulation. As ^{2}. This results in 3.92 ms in the constant conductivity case and 4.77 ms in the reduced conductivity case. The results show that a density of median 0.20, range 0.18–0.26 samples/cm^{2} is needed for the homogeneous case and a density of median 0.24, range 0.16–0.30 samples/cm^{2} for the heterogeneous case. Finally, we compare the performance of the active learning algorithm against randomly choosing samples by running 30 cases at 0.25 and 0.74 samples/cm^{2}.

Evolution of active learning for the homogeneous case with constant conductivity. The top two rows show the activation times for different sample densities and the ground truth for two different views. The bottom rows represent the conduction velocity.

Evolution of active learning for the heterogeneous case with partly reduced conductivity. The top two rows show the activation times for different sample densities and the ground truth for two different views. The bottom rows represent the conduction velocity.

In this work, we present a novel framework to create activation maps in cardiac electrophysiology by incorporating prior physical knowledge of the underlying wave propagation. To this end, we employ neural networks to represent the spatial distribution of activation times and conduction velocity, subject to physics-informed regularization prescribed by the Eikonal equation. In particular, we show that our method is able to capture the collision of wavefronts, a situation that is not captured by other state-of-the-art methods. This is a critical step toward reliably estimating conduction velocities and avoiding artificially high conduction velocity values. Our methodology directly predicts conduction velocities, without the need of creating

Even though our method shows promising results when compared to existing solutions, it displays some limitations. First, we have ignored the anisotropy in conduction of cardiac tissue [_{N} that we use in the randomized prior functions. Second, our uncertainty estimates are only approximations, since the true uncertainties depend on the geodesic distance between points on the manifold and not on the Euclidean distance in

In summary, we have presented a new and efficient methodology to acquire and create activation maps in cardiac electrophysiology. We believe that our approach will enable faster and more accurate acquisitions that will benefit the diagnosis of patients with cardiac arrhythmias.

The datasets generated for this study are available on request to the corresponding author. Source code is available at

FS, PP, DH, and EK designed the research and wrote the manuscript. FS, PP, and YY wrote the code and performed the simulations. FS and PP analyzed the data.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.