In this work, we use previously developed highly detailed human atrial model of atrial fibrillation (AF) (Potse et al., 2018; Gharaviri et al., 2020). We briefly summarize here the relevant aspects of the model. The anatomy, including heart and torso, is based on MRI data. Several key features (bundles, fibers) are based on histological studies and added manually. The atrial wall is 3-dimensional with variable thickness.

In the numerical experiments for this study, we consider different combinations of fibrosis patterns and ablation lines, for a total of 9 scenarios. Firstly, we consider three fibrosis patterns (Figure 2), one case with moderate fibrosis, corresponding to 50% of fibrotic tissue, and two cases with severe fibrosis, corresponding to 70% of fibrotic tissue. We consider endomysial fibrosis, which is modeled by formally imposing zero cross-fiber intracellular conductivity in fibrotic regions. Secondly, we implement two standard-of-care ablation strategies, pulmonary veins isolation (PVI) and PVI with roof lines (PVI+BOX), see Figure 2D. Ablation lines are non-conductive tissue.

The electrical activity is modeled with the monodomain system (Colli Franzone et al., 2014), which reads as follows (1){χ(Cm∂tVm+Iion(Vm,w)     +Istim(x,t))=∇·(Gm∇Vm),(x,t)∈Ω×(0,T],∂tw=g(Vm,w),(x,t)∈Ω×(0,T],Gm∇Vm·n=0,(x,t)∈∂Ω×(0,T],Vm(x,0)=V0, w(x,0)=w0,x∈Ω, where V_m(x, t) is the transmembrane potential, w(x, t) is a vector of ion gating and concentration variables, Ω is a domain describing the active myocardium, C_m = 1 μF cm⁻² is the membrane capacitance, χ = 800 cm⁻¹ is the membrane surface-to-volume ratio, I_stim is the current stimulus, G_m(x) is the monodomain conductivity tensor, and I_ion and g describe the ionic model. In particular, we consider the Courtermanche-Ramirez-Nattel model (Courtemanche et al., 1998) adapted to an AF phenotype, with minor adaptations to guarantee numerical stability when evaluating the gating parameters for certain values of V_m (Potse, 2019). The initial condition (V₀, w₀) corresponds to the resting state.

The conductivity tensor G_m is defined as Gi(Gi+Ge)-1Ge, where G_i and G_e are, respectively, the intra- and extra-cellular conductivity tensors, both assumed transversely isotropic with respect to the local fiber direction. The intracellular longitudinal and cross conductivity are, respectively, set 3 and 0.3 mS cm⁻¹, while the extracellular conductivities are 3 and 1.2 mS cm⁻¹, respectively. The resulting conduction velocity in the fiber direction is 55.6 cm s⁻¹. In the Bachmann's bundle, faster conduction is obtained with a longitudinal intracellular conductivity of 9 mS cm⁻¹. Finally, the region between the superior and inferior vena cava is assumed isotropic, with all conductivities set to 1.5 mS cm⁻¹.

The numerical solution of Equation (1) is based on a second-order finite difference scheme for the spatial discretization, and a fully explicit first-order Euler scheme for time stepping (Potse et al., 2006). The Rush-Larsen scheme is adopted to update the gating variables. The computational domain is discretized using a uniform mesh with hexahedral elements of side h.

For the high fidelity simulations, we consider a fine mesh with h = 0.2 mm and a time step of Δt = 0.01 ms. For the low fidelity simulations, we double the discretization parameters, with h = 0.4 mm and Δt = 0.02 ms. The coarsening of the grid is performed by employing a majority rule to determine the tissue type and fiber orientation of the coarse hexahedral elements from the eight sub-elements of the fine mesh. Moreover, the coarse model assumes a reduced surface-to-volume ratio χ = 450 cm⁻¹ to balance out the expected reduction in conduction velocity due to a coarser space discretization (Pezzuto et al., 2016).

All simulations are performed with the Propag-5 software (Potse et al., 2006; Krause et al., 2012) on the Swiss National Supercomputing Centre (CSCS). For one simulation with T = 4 s, the compute time of the high fidelity model is 1 h40 min with 8 nodes, whereas the compute time of the low fidelity model is 14 min with 4 nodes. This means that the low fidelity model is approximately 16 times faster than the high fidelity model.

The stimulation protocol, encoded in the function I_stim(x, t), is defined by a point x_stim ∈ Ω and a vector of distinct times τstim={τj}j=1Nstim through the expression (2)Istim(x,t;xstim,τstim)={Imax,(x,t)∈Br(xstim)×⋃j=1Nstim[τj,τj+Δτ],0,otherwise, where B_r(x_stim) = {x ∈ Ω:x_stim ≤ x ≤ x_stim + r} is a r-neighborhood (the ≤ is meant component-wise) of the stimulation site, and Δτ > 0 is the stimulus duration. In this study, the vector τ_stim is fixed as in Figure 3 (middle panel), which consists in a series of N_stim = 14 stimuli with decreasing temporal distance, whereas x_stim varies for each simulation. Each stimulus lasts Δτ = 4 ms, and has a strength I_max = 800 μA cm⁻² with a fixed radius r = 0.8 cm, which is enough to maximize the chance that the tissue correctly captures it (Potse et al., 2018). The induction of AF is not successful when (V_m, w) asymptotically approaches the resting state after the delivery of the last stimulus. Otherwise, if a self-sustained activity is still present at the end of the simulation, the induction of AF is successful. The idea is summarized in Figure 3. For sake of simplicity, in this work there is no distinction between a true AF episode and an atrial flutter, which could be understood as a periodic solution of the monodomain system.

The objective of this work is to learn the set A⊂Ω, such that if xstim∈A a sustained episode of AF is observed. In particular, we are interested in approximating the indicator function of A, denoted by F:Ω → {0, 1} such that F-1(1)=A. The overall inducibility, which reflects the fraction of the tissue where AF can be initiated, follows immediately from the definition of A as I=|A||Ω|=1|Ω|∫ΩF(x) dx. Interestingly, the formula generalizes to the case of non-uniformly distributed ectopic foci. Let ρ(x) be the probability density function of the distribution of foci, then the inducibility can be obtained as Iρ=∫ΩF(x)ρ(x) dx. In this way, for instance, it is possible to account for a higher density of ectopic activity around the pulmonary veins and fibrotic regions. In this work, we will only consider a uniform distribution of foci, equivalent to select ρ(x) = |Ω|⁻¹.

Next, we present the proposed methodology for learning the inducibility function F from a limited set of simulations. We start by assuming that we have a data-set of N input/output pairs D={(xi,yi)}i=1N, where x_i ∈ Ω and y_i ∈ {0, 1}. Since the atrial wall is thin, we constrain the points to belong to a mid-wall smooth atrial surface S⊂Ω. We remark however that there is no loss of generality in the following presentation, as the methodology applies to the volumetric domain Ω in the same manner. Moreover, since y_i takes only binary values, we also restrict the scope of this work to binary classification. We also note that it is straightforward to extend this framework to the multi-class classification setting.

The classical formulation of Gaussian process classification defines an inter mediate variable which is computed from a latent function f(x) (Rasmussen and Williams, 2006). Throughout this article, we will assume standardized data-sets and work with zero-mean Gaussian process priors of the form f~GP(0,k(x,x′;θ)). Here, k(·, ·; θ) is a covariance kernel function, which depends on a set of parameters θ. We adopt a fully Bayesian treatment and prescribe prior distributions over these parameters, which we will specify later (Neal, 1999). To obtain class probability predictions we pass the Gaussian process output f through a non-linear warping function σ:ℝ → [0, 1], such that the output is constrained to [0, 1], rendering meaningful class probabilities. We define the conditional class probability as π(x) = ℙ[y = 1|x] = σ(f(x)). A common choice for σ(f) is the logistic sigmoid function σ(f) = (1 + exp(−f))⁻¹, which we will use throughout this work. We assume that the class labels are distributed according to a Bernoulli likelihood with probability σ(y) (Nickisch and Rasmussen, 2008).

A crucial step in building a Gaussian process classifier is the choice of the kernel function. A popular choice is the Matérn kernel, which explicitly allows one to encode smoothness assumptions for the latent functions f(x) (Rasmussen and Williams, 2006). In a Euclidean space setting, the kernel function has the form (Rasmussen and Williams, 2006): (3)k(x,x′,θ)=η221-νΓ(ν)(2ν||x-x′||ℓ)νKν(2ν||x-x′||ℓ), where Γ is the gamma function, and K_ν is the modified Bessel function of the second kind. The parameter η controls the overall variance of the Gaussian process, the parameter ℓ controls the spatial correlation length-scale, and ν controls the regularity of the latent functions f(x) (Rasmussen and Williams, 2006). When ν → ∞, we recover the popular squared exponential kernel, also known as radial basis function, that yields a prior over smooth functions with infinitely many continuous derivatives.

The form presented in Equation (3) is not suitable to be used on manifolds, as the atrial surface. A naive approach is to replace the Euclidean distance between points with the geodesic distance on the manifold surface. Even though this approach may work for some cases, there is no guarantee that the resulting covariance matrix between input points will be positive semi-definite (Pezzuto et al., 2019; Borovitskiy et al., 2020), a key requirement for a kernel function. As a matter of fact, the choice of the kernel is problematic in this case. For instance, the Matérn family does not yield positive definite kernels even on the sphere, except for a few exceptional choices of the parameters (Gneiting, 2013). Here, we follow an alternative approach, implicitly based on the solution of the following stochastic partial differential equation (SPDE) (Whittle, 1963; Lindgren et al., 2011): (4){(κ2I-Δ)ν/2+d/4u=W,x∈Ω,n·∇(κ2I-Δ)ju=0,x∈∂Ω,j=0,…,⌊ν-12+d4⌋, where −Δ is the Laplace-Beltrami operator on the d-dimensional manifold, and W is the spatial Gaussian white noise on Ω. When Ω = ℝ^d, the solution of the fractional SPDE is a Matérn random field with κ=2νℓ (Lindgren et al., 2011). However, compared to Equation (3), the SPDE in Equation (4) trivially generalizes to manifolds with no loss of positive definiteness of the correlation kernel, thanks to the properties of the pseudo-differential operator (Borovitskiy et al., 2020). The correlation function can be explicitly written as follows. Let {(λi,ψi)}i=0∞ be the eigenvalue/eigenfunction pairs of the Laplace-Beltrami operator with pure Neumann boundary conditions, that is (5){-Δψi=λiψix∈Ω,-n·∇ψi=0,x∈∂Ω, for all i ∈ ℕ. Then, we can represent Matérn-like kernels on manifolds as (Coveney et al., 2019; Borovitskiy et al., 2020). (6)k(x,x′;θ)=η2C∑i=0∞(1ℓ2+λi)-ν-d2ψi(x)ψi(x′) where C is a normalizing constant. This eigen-decomposition also enables a direct solution of the SPDE, providing the following representation of the Gaussian process prior: (7)f(x)≈η2C∑i=0∞wi(1ℓ2+λi)-ν2-d4ψi(x),  wi~N(0,1) In practice, the eigen-decomposition is truncated to a number N_eig of pairs.

In this work, we discretize the manifold S⊂ℝ3 using a triangulated mesh and solve Equation (5) using finite element shape functions. As such, we can obtain the stiffness matrix A and mass matrix M: (8)Aij=Ae=1nel∫B∇Nj·∇Ni dx,Mij=Ae=1nel∫BNjNi dx, where A represents the assembly of the local element matrices, and N are the finite element shape functions. Then, we solve the eigenvalue problem: (9)Av=λMv In practice, to compute the kernel in Equation (6) we use a portion of all the resulting eigenpairs, starting from the smallest eigenvalues. We also use the corresponding eigenvectors as the eigenfunctions with f(x_i) = v_i, where i is the node index at location x_i. Given that the eigenvalue problem is solved only once as a pre-processing step, this methodology provides an efficient way to compute the kernel and the prior in a manifold.

We finalize our Bayesian model description by prescribing the prior distributions for the kernel parameters. We assume the following distributions for the parameters θ = {η, ℓ}, (10)η~HalfNormal(σ=10000) (11)ℓ~Gamma(α=1,β=1). The posterior distribution over the model parameters θ = {η, ℓ} cannot be described analytically, and thus we must resort to approximate inference techniques to calibrate this Bayesian model on the available data. To this end, we use the NO-U-Turn sampler (NUTS) (Hoffman and Gelman, 2014), which is a type of Hamiltonian Monte Carlo algorithm, as implemented in NumPyro (Phan et al., 2019). We use one chain, and set the target accept probability to 0.9. The first 500 samples are used to adjust the step size of the sampler, and are later discarded. We use the subsequent 500 samples to statistically estimate the parameters θ.

Once we have completed the inference, we can make predictions y^* at new locations x^* in three steps. First, we compute the predictive posterior distribution of the latent function f*(x*)~N(μ(x*),Σ(x*)), which by construction follows a multi-variate normal distribution, with a mean μ and covariance Σ obtained by conditioning on the available training data (Rasmussen and Williams, 2006): (12)μ(x*)=k(x*,X)K-1f (13)Σ(x*)=k(x*,x*)-k(x*,X)K-1k(X,x*), where the covariance matrix K ∈ ℝ^{N × N} results from evaluating the kernel function k(·, ·; θ) at the locations of the input training data X and f = f(X), respectively. We then proceed by sampling μ, Σ using model parameters drawn from the estimated posterior distributions of θ and f. This will result in a number of random variables f^* that are independent and normally distributed, which we can be used to compute statistical averages as (14)f^*~N(μ^,Σ^),  μ^=1Ns∑i=1Nsμi,  Σ^=1Ns∑i=1NsΣi, where N_s is the number of samples considered for θ and f. We finally pass f^* through the logistic sigmoid function σ to obtain a distribution of class probabilities y^*.

In this work, we will assume that we have 2 information sources of different fidelity. We will call the high fidelity, computationally expensive, and hard to acquire information source with the subscript H and the inexpensive, faster to compute, low fidelity source with the subscript L. Now, our data set comes from these two sources D={(xLi,yLi)i=1NL,(xHi,yHi)i=1NH}={X,y}. We will postulate two latent functions f_H and f_L, respectively, that are related through an auto-regressive prior (Kennedy and O'Hagan, 2000). (15)fH(x)=ρfL(x)+δ(x) Under this model structure, the high fidelity function is expressed as a combination of the low fidelity function scaled by ρ, corrected with another latent function δ(x) that explains the difference between the different levels of fidelity. Following (Kennedy and O'Hagan, 2000), we assume Gaussian process priors on these latent functions. (16)fL~GP(0,k(x,x′;θL)) (17)δ~GP(0,k(x,x′;θH)). The vectors θ_L and θ_H contain the kernel hyper-parameters of this multi-fidelity Gaussian processes model. The choice of the auto-regressive model leads to a joint prior distribution over the latent functions that can be expressed as (Kennedy and O'Hagan, 2000). (18)f=[fLfH]~N([00],[KLLKLHKLHKHH]), with (19)KLL=       kL(XL,XL′;θL)KLH=ρkL(XL,XH′;θL)KHH=ρ2kL(XH,XH′;θL)+kH(XH,XH′;θH). The global covariance matrix K of this multi-fidelity Gaussian process model has a block structure corresponding to the different levels of fidelity, where K_HH and K_LL model the spatial correlation of the data observed in each fidelity level, and K_LH models the cross-correlation between the two levels of fidelity. We also have kernel parameters for the different levels of fidelity. We again use the Matérn as described in Section 2.4, which results in parameters θ_H = (η_H, ℓ_H), and θ_L = (η_L, ℓ_L). For these parameters and the scaling factor ρ, we consider the following prior distributions (20)ηH,ηL~HalfNormal(σ=10000) (21)ℓH,ℓL~Gamma(α=2,β=2) (22)ρ~Normal(μ=0,σ=10). We can perform inference and prediction for this model in the same way as for the single fidelity classifier, as detailed in Section 2.5. In particular, we can use Equations (12) and (13) with the entire covariance matrix K to obtain the conditional mean and covariance of f^*.

Here, we take advantage of the uncertainty predictions that are inherent to Gaussian processes and are absent in other types of classifiers, such as nearest neighbor. Specifically, at each active learning iteration, we train the classifier, and select the next point that should be included in our training data-set by solving the following optimization problem (Kapoor et al., 2007; Sahli Costabal et al., 2019): (23)xnew=arg minx∈Xcand|μ^(x)|Σ^(x), where X_cand represents a set of candidate locations that can be acquired. In our case, we use all the nodes in the mesh as candidates, except the ones at the boundaries which have artificially high variance. This active learning criterion presents a good balance between exploitation (sampling near the classification decision boundary) and exploration (discovering new inducible regions). It can be seen as promoting the selection of points that tend to be located near the decision boundary (σ(μ = 0) = 0.5), or points in regions with high uncertainty (as reflected by the posterior variance Σ). We keep adding points via this sequential active learning procedure until we have reached the desired number of samples.

We first create a synthetic example to test the performance of the proposed classifier. We study the length scale of different random fields that could represent the potential inducibility maps that we want to approximate in this study. In particular, we use a mesh based on the anatomy of the mid-layer described in Section 2.1. Here, we represent the left and right atria with 3,298 nodes and 6,335 triangles. First, we normalize the geometry by the largest standard deviation of one of its coordinates. In this way, we can use the same prior distributions regardless the particular geometry. Then, we generate Gaussian random fields on the atrial manifold with zero mean and the Matérn covariance kernel, as detailed in Equation (6). We use 1,000 eigenpairs to construct a computable kernel function approximation with ν = 3/2 and η = 1. We consider different length scales to simulate inducibility regions and assess the performance of the classifier: ℓ = {0.2, 0.4, 0.6, 0.8, 1.0}. Finally, we pass the resulting random field through the sigmoid function σ to obtain values between zero and one, which we round to the nearest integer to create discrete labels.

Examples of the resulting random fields can be seen in Figure 4, left column. We compare three different classifiers. First, as a baseline benchmark, we create a nearest neighbor classifier. Here, the prediction of an unknown point is based on the label of the closest data point. Since we are working with a manifold, we use the geodesic distance to find the closest point, which we compute using the heat method (Crane et al., 2013). As a data-set, we use a fixed design spread through the manifold surface. To select the locations, we first randomly pick a node in the mesh, and then we add the node that is further away from the initial node using the geodesic distance. Then, we iterate, finding the point that is further away from all the nodes already included in the data-set, until we reached the desired data-set size. The second classifier that we consider is a Gaussian process classifier, as described in the previous sections, that is trained on the same fixed experimental design. The final classifier is also a Gaussian process classifier, which we train with the first 20 samples of the fixed experimental design, and then we apply the proposed active learning procedure. For all the Gaussian process classifiers in these experiments, we set the number of eigenfunctions used to N_eig = 1,000.

In these examples, we test the performance of the three different classifiers, using between 20 and 100 samples, and 10 different random fields for each of the 5 length scales selected. To take into account potential imbalances of classes in the examples generated, we use the balanced accuracy score. This metric is defined as the arithmetic mean of the sensitivity and specificity as (24)balanced accuracy=12(# of predicted positives# of real positives+# of predicted negatives# of real negatives). In contrast to conventional accuracy, this metric will reflect if a classifier is predominately predicting one class due to the higher proportion of samples present in the data-set.

The results of this assessment are summarized in Figures 4, 5. We first observe in Figure 4, left column, that the complexity of the classification regions increases as the correlation length scale is reduced. In the same figure, we show the different classifiers trained with 100 samples. It is visually possible to note that the accuracy of the classifiers degrades as the length scale of the ground-truth classification surface is decreased. For the length scale ℓ = 0.2, some regions are not captured by the classifiers. We also note that the Gaussian process classification boundaries tend to be smoother than the nearest neighbor classifier. These differences are quantified in Figure 5. We first compare the improvements in accuracy between the nearest neighbor classifier and the Gaussian process classifier with a fixed design in the top row. These two methods are trained with identical data, and we observe that for most cases and number of samples, the Gaussian process classifier is more accurate than the nearest neighbor classifier. The accuracy improvements at 100 samples range on average from 0.8% at ℓ = 0.2 to 2.4% at ℓ = 0.8. Then, we compare the nearest neighbor classifier with the Gaussian process classier trained with active learning. These two classifiers only share the first 20 points of data. Then the active learning classifier judiciously selects the remaining samples attempting to maximize accuracy. We observe that the accuracy improvements are more pronounced with the active learning for ℓ = 0.4 − 1.0. The average improvements at 100 samples range from 3.0% at ℓ = 0.4 to 6.2% at ℓ = 0.8. For ℓ = 0.2, we see an average decrease in accuracy of 1.0% at 100 samples. In the last row of Figure 5, we see the average accuracies for the three classifiers at 100 samples, reflecting the improvements in accuracy obtained by the Gaussian process classifiers already described. We see that all classifiers tend to decrease their accuracy as the length scale decreases, which coincides with the increased complexity of classification boundaries for lower length scales seen in Figure 4. This detriment in performance becomes more pronounced between ℓ = 0.4 and ℓ = 0.2. This change corresponds with the average geodesic distance between points in the fixed design data-set, which is equal to 0.39. This metric is shown as a dashed vertical line bottom row plot of Figure 5. Classification regions with a characteristic size smaller than this value could be ignored by the classifiers, which is what we observe in the top row of Figure 5. In these cases, the uncertainty estimates used for active learning might be inadequate, leading to a worse performance compared to the longer length scale cases. Overall, we see that Gaussian process classifiers and active learning provide advantages in accuracy when compared to the baseline nearest neighbor classifier.

We examine the inducibility of the 9 models described in Section 2.1, specifically 3 different fibrotic patterns and 3 ablation strategies: no ablation, PVI, and PVI+BOX. For each model, we create a training set and test set, both containing 100 samples, using a fixed design, as described in Section 3.1 and shown in Figure 6. We run the model using each of these points as a pacing site and check whether AF was induced or not. For the training set, we also run the low fidelity model, obtaining 100 samples. In total we run 1800 high fidelity simulations and 900 low fidelity simulations. We test three different classifiers for both cases: a nearest neighbor classifier described in Section 3.1, a single-fidelity Gaussian process classifier described in Section 2.3, and a multi-fidelity Gaussian process classifier described in Section 2.6 with 100 low fidelity samples. We train the classifiers with different amounts of data from the training set, ranging from 20 to 100 points. For each level of data, we evaluate the performance of the classifier computing the balanced accuracy in the 100 samples of the test set.

The results of this numerical experiment are summarized in Figures 6–8 and Table 1. First, we note that training and predicting with the Gaussian process classifiers only takes a negligible fraction of the cost of high fidelity model, less than 5 min on a laptop. In Figure 6, we show the resulting classifiers trained with the same 50 and 100 high fidelity samples and also the low fidelity classifier trained with 100 samples. We see that the multi-fidelity classifier at 50 and 100 samples shares some features with the low fidelity classifier that are not present in the other two classifiers. Nonetheless, the multi-fidelity classifier is learning from the high fidelity data, as its balanced accuracy increases as the number of samples increases, as seen in Figure 8. We observe that the differences in accuracy tend to collapse as more data is available, showing small differences when 100 samples are provided to the classifiers. The multi-fidelity classifier has the biggest advantage in the small data regime, when it is trained with between 20 and 70 high fidelity samples. Perhaps surprisingly, we see that the low fidelity classifier is always more accurate than the single-fidelity classifiers trained with 20 samples. The cost of training the low fidelity classifier is approximately equivalent to the cost of acquiring 6.25 high fidelity samples, which makes it a cost-effective alternative to estimate the inducibility with limited budget. Along the same line, we compare the accuracies of the different classifiers for the different cases when the training with the equivalent cost of 40 high fidelity simulations in Figure 7B. This is the number of simulations that has been used in clinical studies to optimize the ablation treatment (Boyle et al., 2019). We observe that by using the multi-fidelity classifier we gain, on average, 5.4% points of accuracy comparing to the single-fidelity classifier and 5.7% comparing to the nearest neighbor classifier. Only in one case there was a decrease in accuracy when using the multi-fidelity classifier, but of only 0.45% points of accuracy.

We analyze the agreement between the low and high fidelity models by looking at training sets for all cases in Figure 7A. Overall, we find the low and high fidelity agree in 81.7% of the simulation. However, we see that the low fidelity model is biased towards predicting no AF when the high fidelity model is predicting AF. This is confirmed in every case, as can be seen in Table 1, where low fidelity inducibility is always lower than the high fidelity inducibility. A possible explanation is that the low fidelity model, being based on a coarser discretization of the atrial model, has fewer fine-grained features (fibrosis, anatomy, wall thickness) that might favor AF. It is also worth noting that we adapted the conduction velocity in the low fidelity model by increasing it to the level of the high fidelity one, a change that is potentially antiarrhythmic but that increased the correlation between the models and hence the overall performance of the multi-fidelity classifier. We also found that in the case of 50% fibrosis, the low fidelity model tends to predict proportionally more occurrence of AF when the high fidelity model is not predicting AF.

Finally, we see in Table 1 that the ablation strategies applied are decreasing the inducibility in all cases, both for the train, test and low fidelity sets. We see that pulmonary vein isolation has more impact on the inducibility than the subsequent box ablation for all cases, both in the train and the test set.