Calibrating cardiac electrophysiology models using latent Gaussian processes on atrial manifolds

Models of electrical excitation and recovery in the heart have become increasingly detailed, but have yet to be used routinely in the clinical setting to guide personalized intervention in patients. One of the main challenges is calibrating models from the limited measurements that can be made in a patient during a standard clinical procedure. In this work, we propose a novel framework for the probabilistic calibration of electrophysiology parameters on the left atrium of the heart using local measurements of cardiac excitability. Parameter fields are represented as Gaussian processes on manifolds and are linked to measurements via surrogate functions that map from local parameter values to measurements. The posterior distribution of parameter fields is then obtained. We show that our method can recover parameter fields used to generate localised synthetic measurements of effective refractory period. Our methodology is applicable to other measurement types collected with clinical protocols, and more generally for calibration where model parameters vary over a manifold.

Mechanical contraction of the heart is initiated and synchronised by a travelling wave of electrical excitation and recovery that arises spontaneously in the natural pacemaker. The heart is made up of four chambers: the ventricles pump blood to the body and lungs, while the atria act as reservoirs and primers for the ventricles. A cardiac arrhythmia is a disturbance of regular heart rhythm resulting in a rapid, slow, or irregular rhythm. Atrial fibrillation (AF) is a common and increasingly prevalent cardiac arrhythmia 1 . AF can be sustained by re-entry, where electrical activation continually propagates into recovering tissue, creating a self-sustaining rotating wave 2 . Radio-frequency catheter ablation can be used to disrupt re-entrant circuits that act to sustain AF, but is not always effective 3 .
Two properties of cardiac tissue are important for the development of sustained re-entry, and these properties vary across atrial tissue. Conduction velocity (CV) describes the speed at which an activation wave spreads. The effective refractory period (ERP) is the minimum time interval between two successive stimuli that allows two activation waves to propagate and is related to action potential duration (APD), which is the interval between local activation (depolarization) and recovery (repolarization). Both CV and ERP decrease at shorter pacing intervals, and this dynamic behaviour and its spatial heterogeneity is important for determining the stability of re-entry 4,5 as well as the complex paths followed by electrical activation during AF 6 . Natural variability in the speed of the excitation wave and the dynamics of excitation and recovery exist both between individuals and within the heart of a single individual 7,8 . Cardiac tissue exhibits spatial heterogeneity with differences in ion channel conductances, gap junction distributions, and fibrotic remodelling across the heart 9 . These spatial heterogeneities in structural and functional properties lead to heterogeneity in ERP. The resulting dispersion in repolarisation properties is a mechanism for focal arrhythmia initiation 10 , and atrial fibrillation initiation through increasing vulnerability to re-entry 8,11 . Electrophysiology (EP) models describe how electrical activation diffuses through cardiac tissue. Local activation and recovery are represented by a set of differential equations describing a reaction-diffusion system

Results
Workflow. The computational model, or 'simulator' , that we seek to calibrate is composed of (i) a finite element mesh representing an atrial manifold x ∈ ; (ii) an electrophysiology model, that maps EP parameters θ l (x), l = 1, 2, . . . defined on the computational mesh to observable quantities; and (iii) a numerical solver for running EP simulations. Details on obtaining and processing a mesh for suitability in electrophysiology simulations are given in "Methods", including details on the example mesh used here. The EP model is the modified Mitchell-Schaeffer (mMS) model 23,24 , the parameters of which are effectively time-constants representing different phases of the action potential. We parameterize the mMS model with the following 5 parameters: CV max (x), τ in (x), τ out (x), τ open (x), APD max (x) . See "Methods" for details on the simulation model, parameterisation, and allowable parameter ranges, and details on the numerical implementation. We use the software open-CARP 25 to solve the mono-domain model for our simulations.
The main task in this work is to calibrate the simulator by inferring parameter fields θ l (x) from ERP measurements. Figure 1 represents our modeling workflow, which we summarize here. Our code is available in a Zenodo repository 26 .
Surrogate functions. The simulator can be used to map parameter fields to ERP fields ERP(x) = sim(θ l (x), l = 1, 2, . . .) . Given ERP observations at multiple locations on the atrial mesh, as well as an appropriate likelihood model for these observations, the simulator could be used in an MCMC setting in order to calibrate the parameter fields by obtaining samples from the posterior distribution for the EP fields. However, this is an extremely inefficient approach, since ERP depends only on local (rather than remote) tissue properties. We utilize a surrogate function (also called an "emulator") solution in which we learn the mapping from parameters to ERP. This surrogate function allows us to predict ERP at location x as ERP(x) = f (θ l (x), l = 1, 2, . . .) , bypassing the need to run the simulator directly for inference.
Gaussian process priors. The mesh has ≈ 10 5 vertices for which parameters need to be defined, but ERP measurements are restricted to a subset of these vertices, with number of observations on the order of 10 0 -10 1 . The electrophysiology parameter fields must be assumed to have low-rank structure, induced by spatial correlation, in order to make inferences about EP parameter values at locations other than ERP observation locations. This is achieved here by modeling the EP parameters using latent Gaussian process (GP) priors θ l (x) ∼ GP . We use Gaussian Process Manifold Interpolation (GPMI), a method we proposed for defining Gaussian process (GP) distributions on manifolds 27 . The approach uses solutions { k , φ k (x)} of the Laplacian (Laplace-Beltrami) eigenproblem on the mesh 28 .
Bayesian calibration. We perform probabilistic calibration with MCMC to obtain the posterior distribution of latent variables in the GPs. We utilize a likelihood function that we developed specifically for ERP measure-  Figure 2a shows sensitivity indices for two types of ERP:

Sensitivity analysis and surrogate functions.
an ERP measurement for S1S2 with S1 600 ms, denoted here as ERP S2 , and another type of ERP measurement for S1S2S3 pacing for S1 600 ms S2 300 ms, denoted here as ERP S3 . The S1S2S3 protocol, consisting of N S1 beats, 1 S2 beat, and 1 S3 beat, is introduced in this paper. We have determined that these measurements can be Figure 1. Workflow for using ERP observations to calibrate model parameters. A simplified 'surrogate simulation' consisting of a strip of tissue paced from one end is used to determine ERP values for a design of parameters. Surrogate functions are fit to these data, allowing rapid prediction of ERP from model parameters.
Electrophysiology parameter fields are modelled as Gaussian processes on the atrial manifold, using a reducedrank formulation relying on eigenfunctions and eigenvalues of the Laplace-Beltrami operator on the atrial mesh. ERP measurements using an S1S2 (or S1S2S3) protocol measure whether successful activation occurs for S2 (or S3) times t 1 < t 2 < t 3 . . . etc. Using a likelihood function designed for ERP measurements and the surrogate functions, the likelihood can be evaluated given hyperparameters for the GP models. This allows for probabilistic calibration by obtaining the posterior distribution of hyperparameters using MCMC. www.nature.com/scientificreports/ used to calibrate EP parameters sufficiently to reproduce not only these ERP measurements, but also the time required for the action potential to reach various levels of repolarization recovery (e.g. APD 20 and APD 90 , the time required for 20% and 90% recovery). It is a key finding that the S1S2S3 protocol can be used (alongside the standard S1S2 protocol) to disentangle the contributions of separate parts of the action potential to the value of ERP, without needing to measure the action potential directly. The sensitivity indices in Fig. 2a show that these ERP measurements are mainly determined by τ out and APD max , which approximately correspond to the duration of the repolarization and plateau phases of the action potential respectively. Calibration of other parameters, which determine some aspects of the shape of restitution curves but which do not strongly impact ERP, require both CV and APD restitution curve data from an S1S2 protocol 21 . For this reason, we have determined to use ERP to calibrate θ 1 ≡ τ out and θ 2 ≡ APD max . Figure 2b shows contour plots of the surrogate functions for ERP. A discontinuity occurs in the ERP S3 surface for parameters combinations resulting in ERP S2 > 285 ms, so data for ERP S2 > 280 ms were discarded before fitting this function. Note that the majority of clinical ERP S2 measurements fall in the range 170-270 ms, so even 280 ms could be considered as an upper limit 29 . We used a left atrial mesh generated from a scan of an individual performed at St Thomas' Hospital (see "Methods" for details). We created ground truth parameter fields for τ out and APD max in order to verify our calibration approach. We used 10 measurement locations, placed at random using a maximin design that excluded sites close to the mesh boundaries. The resolution of the S1S2 and S1S2S3 protocol was set to 10 ms. We used 24 eigenfunctions for representing each of the two parameter fields in Eq. (5), which we found to be sufficient to capture spatial variation while allowing good posterior sampling. For MCMC we used 5000 iterations, using 8 chains, discarding the first 50% of the samples as 'burn-in' , and randomly thinning the remaining samples by a factor of 100 to give 200 posterior samples. Figure 3 shows the true parameter fields, and the posterior mean and standard deviation of the calibrated parameter fields. Figure 4 shows the true ERP fields, the posterior mean and standard deviation of the ERP fields (calculated from ERP samples, which are calculated from the parameter field posterior samples), and the Independent Standard Errors (ISE) of ERP (the absolute difference between true and posterior mean, divided by the posterior standard deviation). Measurement locations are shown as spheres in Figs. 3 and 4, colored by the corresponding values at each location. Figure 5 shows the APD simulation results from the atrial simulator using the ground truth parameter fields and the posterior mean of the calibrated parameter fields.
The prediction of EP parameter fields τ out and APD max and ERP fields ERP S2 and ERP S3 captures the ground truth extremely well. Predictions on the pulmonary veins, which are effectively regions of extrapolation, deviate from the ground truth more than other regions on the main body of the atrium. These deviations are on the order of the S2 and S3 resolution, and the posterior variance is higher in these regions. Uncertainty increases with distance from the measurement locations. The ISE scores show that the distribution of ERP predicted by the model covers the ground truth well as nearly all values are less than 3. The ISE for ERP S2 on the left atrial appendage is above 3, which may be caused by a combination of high ground-truth values for τ out (which are not effectively probed by the measurements) and insufficient basis functions to capture high spatial variation in this region of the mesh. APD from the full atrial simulator using the posterior mean of the parameters (the maximum a posteriori estimate could have been used instead) matches the ground truth values very closely, demonstrating that the action potential has been calibrated well using only ERP measurements.
We also performed quantitative validation across a broad range of designs. Figure 6 shows these validation results, for different configurations of the S1S2(S3) pacing protocol (number of ERP observations, resolution of S2 and S3 intervals) and different heterogeneity for ERP, controlled by different correlation lengthscales for www.nature.com/scientificreports/ generated APD max and τ out ground-truth fields. A unit of kernel lengthscale is approximately 3.2 mm for this mesh; see "Methods" for details. Prediction values of ERP are based on the maximum a posteriori estimate of the parameters, and here we use 32 eigenfunctions per EP parameter field in order to better model fields with more rapid spatial variation. Root Mean Squared Error (RMSE) is reduced with increasing lengthscale (less ERP heterogeneity), decreasing S2 and S3 resolution (more precise measurements), and increasing number of observations. We note that our likelihood function introduces a small amount of bias, discussed below, which for S2 and S3 resolution 10ms causes RMSE to increase slightly from 20 to 40 observations. Overall, the quantitative validation suggests that little is gained above 20 observation locations.

Discussion
In this paper, we have developed a workflow for calibrating an electrophysiology simulator from sparse measurements of excitability. This was done by representing the spatially varying parameter fields as Gaussian processes on a manifold, and linking these parameters to excitability observations through non-linear surrogate functions (emulators). Using a likelihood function for ERP observations, we performed probabilistic calibration to obtain the posterior distribution of the EP parameter fields. Both visual and quantitative comparison demonstrates that this workflow can successfully calibrate a simulator to ERP to a high level of accuracy. The nature of ERP observations, in which only the interval containing ERP is observed (and the possible brackets around this interval are fixed by the S1S2(S3) protocol), is that the ability to learn more by adding www.nature.com/scientificreports/ observations is strongly limited above a certain point. Figure 6 demonstrates that this limit is reached faster for smaller S2 and S3 resolution. Our likelihood function does introduce a very small amount of bias, since the true likelihood should be constant in the pacing interval, but our approximation decreases on approaching the interval edges. A simple solution would be to pad the ERP observation brackets, which would remove the bias but reduce the precision. Without the assumption that measurements at locations give information about quantities at nearby locations, i.e. spatial correlation, inference about tissue properties beyond measurement sites would not be possible and atrial tissue would need to be sampled everywhere. Such regularization might make it difficult to capture discontinuous changes in tissue properties, although it would be difficult to measure such abrupt changes in tissue behaviour using sparse measurements. It may be possible to utilize other personal data (e.g. scans) or prior information (e.g. a database of clinical measurements) to assist with calibration.
The latent Gaussian process model serves two purposes. Firstly, a run of the electrophysiology model requires specification of parameters at all points on the mesh, and the Gaussian process enables this specification via interpolation between measurement locations. Secondly, we assume that parameter values at neighbouring locations on the mesh are likely to be similar, which means that we need to do joint inference for the parameters at  www.nature.com/scientificreports/ the measurement locations, rather than inferring parameters at each measurement location independently. In developing our method, we first attempted such an independent inference approach, in which parameters are calibrated at each measurement location independently and then interpolated over the manifold using GPMI, but we were not able to obtain satisfactory results. Our current workflow easily allows more complex spatial modeling using multiple latent GPs per EP parameter field, each with independent covariance kernels and hyperparameters that can be freely given suitable priors. It also provides the benefit of being able to constrain the posterior distribution by directly manipulating the posterior samples based on a priori knowledge, such that parameter values (or the tissue properties depending on these parameters) should fall within a certain physiological range. Our proposed workflow for calibration is suitable for other types of data. We have previously shown that Gaussian processes can be used as surrogate functions for CV, APD, and ERP restitution curves 21 . Observations from these restitution curves at different locations over the atrium could be included in calibration simply by including additional contributions to the likelihood function and using "Restitution Curve Emulators" to map from EP parameters to the corresponding restitution curves. Our approach here solves the problem of representing the EP parameter fields on a manifold so as to make probabilistic calibration to sparse measurements into a tractable problem. This allows for propagating uncertainty from measurements through to an ensemble of calibrated models.

Methods
Electrophysiology model. The modified Mitchell-Schaeffer (mMS) cell model 23,24 for mono-domain tissue simulations with isotropic diffusion is expressed in the following equations: where V m is a normalised membrane voltage, h is a gating parameter that controls recovery, and J stim is an externally applied stimulus. The 4 cell model parameters τ = (τ in , τ close , τ out , τ open ) are time-constants that approximately characterize stages of the action potential sequence, and D is conductivity. We fixed the excitation threshold V gate to 0.1. As in 21 , we reparameterized the model as follows: In this new parameter space, weighted combinations of valid parameters are also valid parameters, which means that spatial interpolation of valid parameters will produce valid parameters. We refer to these transformed parameters simply as 'parameters' . The valid ranges of these parameters are set as CV max 0.1-1. Atrial mesh. To generate the mesh for the simulator, the left atrial blood pool was segmented from a contrast enhanced magnetic resonance angiogram scan performed at St Thomas' Hospital 30 . This segmentation was meshed using a marching cubes algorithm in CEMRGApp 31 , and the resulting surface was remeshed to a regular edge length of 0.3mm using mmgtools software 32 , corresponding to around 110,000 vertices, which is sufficient for simulation with the MMS model. This mesh can be found here 33 , and is also included with our code 26 .

Sensitivity analysis.
To determine ERP S2 , the ERP value under an S1S2 protocol for S1 600 ms, and ERP S3 , the ERP value under an S1S2S3 protocol for S1 600 ms and S2 300 ms, we utilized a surrogate simulation: a strip of tissue with homogeneous parameters, paced from one end with the corresponding protocol, with activation measured in the strip centre 21 . The strip simulation is set up to match the atrial simulation as closely as possible (space and time discretization, cell model time-step subdivision, numerical integration, etc). We obtain simulation results with an optimized Latin hyper-cube design of 500 parameter combinations in the parameter range explained above.
Variance-based sensitivity analysis was performed by fitting a General Additive Model (GAM) to model outputs, e.g. ERP S2 , as a function of a single model input, e.g. APD max . The expectation of the GAM is then a line through a point-cloud of input-output pairs. The variance of this line (evaluated at the inputs) divided by the variance of the point-cloud gives an approximate sensitivity index of the input on that output 34,35 . This method can be repeated for all inputs and all outputs. We implement GAMs using the LinearGAM function with 10 splines from the Python module PYGAM 36 . The sensitivity index of output y for input x can then be calculated as puts), such as restitution curves, has been modelled previously using Gaussian processes 21 . Here, cubic polynomials in both θ 1 = τ out and θ 2 = APD max were fit to corresponding values of ERP S2 and ERP S3 , generated from an optimized Latin hyper-cube design of 100 values of CV max , τ out and APD max , keeping τ in = 0.05 ms and τ open = 120 ms in order to produce ERP and APD values in a range observed in human atrial tissue 29 . CV max was varied for robustness, but has negligible effects on ERP, as confirmed by negligible fitting residuals. There is a discontinuity in ERP S3 for parameter values producing ERP S2 ≈ 285 ms, so data for ERP S2 > 280 ms were discarded before fitting these functions. We refer to these polynomial fits as 'surrogate functions' for ERP S2 and ERP S3 , denoted as f 1 (θ 1 , θ 2 ) and f 2 (θ 1 , θ 2 ) respectively, as they allow for determining ERP without running simulations.
Gaussian process priors. We model the EP parameters fields, θ l (x) , as spatially correlated random fields defined on the atrial manifold, i.e. x ∈ . We use Gaussian Process Manifold Interpolation (GPMI), a method we proposed for defining Gaussian process distributions on manifolds 27 . The approach uses solutions { k , φ k (x)} of the Laplacian (Laplace-Beltrami) eigenproblem on the mesh 28 . Using GPMI allows us to represent fields on the atrium using a coordinate system that uses these eigenfunctions as a basis, enabling us to calibrate parameter fields on any given atrial manifold. The prior for each parameter field θ l (x) can then be represented using the following probabilistic model, which uses the K smallest eigenvalue solutions to the Laplacian eigenproblem: with hyperparameters mean m l , amplitude α l , and lengthscale ρ l , for l = 1, 2 and k = 1, . . . , K . The lengthscales determine the distance over which values are correlated, with larger lengthscale corresponding to smoother parameter fields. The units of the lengthscale hyperparameters are determined by the spatial units of the mesh on which the eigenproblem is solved. See below for details. The hyperparameter vector η l ∈ R K must be given a Gaussian prior in order for this model to approximate a Gaussian process. The function S √ k , ρ l is the spectral density corresponding to the choice of covariance kernel, with the square root of the eigenvalue √ k being the 'frequency' argument to this function. In this work, we use the spectral density for the radial basis function (exponentiated quadratic) kernel, but other stationary kernels could be used.
It is possible and tractable to perform inference directly on the simulation mesh by solving the Laplacian eigenproblem on this mesh. But for convenience, a lower resolution mesh of 5000 vertices was used, with vertices that are a subset of the higher resolution mesh vertices. The lower resolution mesh is produced using a simulated annealing algorithm to optimally choose a subset of 5000 nodes before meshing these new nodes to form a new surface. The routines for this are found in the quLATi package 37 , and method details are given in our previous work 27 . This lower resolution mesh allows for calculation of eigenfunctions with fewer computational resources, less data storage for eigenfunctions, and is convenient for plotting. Values can be easily transferred to the simulation mesh via interpolation (we use the 'interpolate' function in the Python software Pyvista 38 ). However, this 'two-mesh' approach is entirely optional. The units of lengthscale parameter ρ in Eq. (5) can be empirically related to geodesic distance by drawing many GP samples for a given kernel function on the mesh (using a lengthscale that allows for correlations to approach zero for some pairs of vertices since the mesh is finite), calculating the correlation between these samples at many pairs of vertices, and fitting (via least squares) the kernel function to correlation as a function of geodesic distance between the pairs of vertices. For the mesh in this work, 1 unit of kernel lengthscale corresponds to approximately 3.2 mm.
Bayesian calibration. Given ERP measurements at different locations x i over the atrium, it is possible to calibrate the parameter fields θ 1 (x) ≡ τ out (x) and θ 2 (x) ≡ APD max (x) by obtaining the posterior distribution of the hyperparameters in Eq. (5). For convenience, we collect the hyperparameters into the vector ψ := (m 1 , m 2 , α 1 , α 2 , ρ 1 , ρ 2 , η 1 , η 2 ) . Defining the ERP measurements as y , then we can write the Bayesian inference problem as: where y 1 (x i ) and y 2 (x i ) represent observations of ERP S2 and ERP S3 respectively. We assume that both types of ERP are measured at each location, but this is not a requirement as terms can just be replaced with 1 if the corresponding measurement is not performed. Clinically, S1S2 protocols are performed by decreasing S2 by S2 until successful activation does not occur on the S2 beat. Therefore, observations of ERP are only observations of an interval in which ERP lies. The observation that each ERP value at a measurement location x i lies between two S2 values t s and t s+1 can be expressed in the following way (see Fig. 1 for a graphical representation): where ψ 1 and ψ 2 represent partitions of ψ for each EP parameter field.
An S1S2 pacing protocol to determine ERP effectively measures the S2 interval in which ERP lies. Defining the lower bound of this interval by I and the interval width by S2 , the true likelihood is given by a truncated uniform, or 'top-hat' , distribution. In other words, p(ERP ∈ [I, I + �S2] | ψ) is equal 1 if ψ produces ERP in the specified interval, and 0 otherwise. The surrogate functions f m (θ 1 , θ 2 ) can be used to predict ERP from the EP parameters, which are determined by the GP fields θ l (x) depending on the hyperparameters ψ.
However, it is more convenient to work with an approximation to this top-hat distribution, which we previously derived for use with ERP measurements 21 . This top-hat likelihood can be approximated by dividing the interval into N sub-intervals, with a normal distribution N (c i , s) centered on each sub-interval c i = I + (i − 1/2)�S2/N with standard deviation equal to the sub-interval width s = �S2/N . We choose N = S2 , such that s = 1 . For an observation y 1 (x * ) := ERP S2 ∈ [I, I + �S2] (and similarly for ERP S3 ) the likelihood can be approximated as: The shape of this likelihood function is a top-hat with smoothed sides and no discontinuities, such that the likelihood is approximately constant in the interval but rapidly falls to zero near the interval edges. This approximate top-hat distribution has infinite support, allowing gradient-based MCMC to be performed. For the loglikelihood, the readily available logsumexp function is used to prevent numerical underflow (this function is available in STAN 39,40 , which is used for MCMC for this work). Note that this approximate top-hat distribution integrates to 1, rather than having approximately constant value 1 (in the interval), but constant factors do not matter for MCMC so we retain this form for simplicity. Also note that with a small adjustment this likelihood can be used with a Gaussian process surrogate function that predicts mean and variance 21 .
We use STAN (via PyStan) 39,40 to perform Hamiltonian MCMC, which yields samples from the posterior distribution of the parameters ψ . See Results for details. We can use these samples to calculate samples of the EP fields over the entire atrium using Eq. (5), from which ERP samples can be calculated using the surrogate functions. We modify Eq. (5) by replacing α l with α l /|� 1 | (where � 1 (x) is constant), which assists with defining priors for α l . We used the following priors on the hyperparameters: ρ l ∼ InvGamma(1.01, 20) , m l ∼ Uniform(−∞, +∞) , and α l ∼ InvGamma (1,5) . We found that these priors consistently allow for recoverability of both EP parameter fields. Eq. (6) gives the prior for the remaining hyperparameters.
Parameter samples. To generate 'ground truth' parameter fields, we draw samples from a Gaussian process defined by Eq. (5) with Matern 5/2 spectral density function using 256 eigenfunctions. We set parameters m = 0 and α = 1 , and ρ is set to values explained in Results and below. The generated samples for τ out and APD max are then scaled and offset into the full allowable parameters ranges. The same operation is performed for CV max , which is needed for atrial simulations.
Certain combinations of τ out and APD max correspond to regions of parameter space that produce unrealistic ERP. This is handled for both 'ground truth' samples and posterior samples by identifying mesh nodes where the parameters produce values ERP S2 > 280 ms . The parameter values at these nodes are then replaced by a weighted average of parameter values at other nodes with acceptable ERP values, weighted by 1/d 4 BH where d BH is biharmonic distance. Biharmonic distance, calculated from the Laplacian eigenvalues and eigenfunctions, is significantly cheaper to calculate than geodesic distance, and avoids topological issues from using Euclidean distance 41 to interpolate values on a manifold. This procedure allows to constrain parameter samples efficiently and effectively, and is far simpler than attempting to encode such constraints into MCMC. Synthetic experiment. We created ground truth parameter fields for τ out and APD max in order to verify our calibration approach (see above for details). ERP values were calculated using the surrogate functions. For ERP measurements, we generate a design of measurement locations using an optimized 'maximin' hypercube design, excluding mesh sites within 0.6 cm of mesh boundaries as potential sites for these measurements since clinical measurements are unlikely to sample these regions. The resolution of the S1S2 and S1S2S3 protocols are set to values specified in Results. A lengthscale of 20 was used for the example shown in Figs. 3, 4 and 5.
APD values were obtained using the atrial simulator (simulations of the mono-domain equation with the mMS model using the software openCARP). Tissue was paced for 8 beats from near the coronary sinus, and depolarization and repolarization were measured on the final beat. A spatially varying CV max field was generated for use in this simulation, and τ in and τ open were fixed as described above. Simulations were run either for ground truth parameter fields or for predicted parameter fields resulting from calibration. Note that the parameters described in this manuscript were transformed back into the original parameters for the mMS model for running www.nature.com/scientificreports/ simulations in openCARP. The diffusion time-step was 0.1 ms and the ionic current time-step was 0.02 ms. The mMS action potential is normalized to have minimum 0 and maximum 1. Activation (depolarization) was measured when V m reached 0.7 on upstroke, and recovery (repolarization) was measured when V m fell to 0.8 ( APD 20 ), 0.7, ( APD 30 ), 0.5 ( APD 50 ), and 0.1 ( APD 90 ). APD values are the time between activation and recovery.
Simulations results for APD 20 and APD 90 are given in Results.