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.

Introduction

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 that models tissue-scale propagation of activation and cellular activation and recovery[12,13]. Models of cardiac electrical activation have become valuable research tools, but are also beginning to be used in the clinical setting to guide interventions in patients[14,15]. These applications require personalised models of both anatomy and electrophysiology to be constructed. Personalised anatomical models can be assembled from medical images, and statistical shape models enable the assessment of varying shape on electrical behaviour[16]. Calibration of EP models is difficult because of the limited measurements that can be made routinely in the clinical setting.

EP model parameters determine model behaviour and for a personalised model should be calibrated to reconstruct the heterogeneity in CV and ERP, as well as their dynamic behaviour, in the heart of a specific patient. Measurements of local activation time (LAT), which measures the time of arrival of the activation wavefront relative to the timing of a pacing stimulus, enable reconstruction of heterogeneous CV for pacing at a fixed rate[17,18]. Calibration to the dynamics of activation and recovery is more challenging. Both the quantity and type of data that can be recorded from patients are constrained by the clinical procedure, so it is difficult to determine spatial heterogeneity of repolarisation. An S1S2 pacing protocol can be used to measure restitution curves. The heart is paced for several beats at an initial pacing cycle length S1, followed by a stimulus with a shorter length S2. This protocol is repeated for different values of the S2 interval, and the shortest S2 that can elicit an activation indicates an upper bound for ERP at the stimulus site. While models can be calibrated to reconstruct CV(S2) restitution and ERP from LAT measurements with an S1S2 protocol[19,20], recent work raises doubts over whether model parameters can be identified uniquely from these types of measurement alone[21]. Biophysically detailed models of electrical activation have large numbers of parameters and many of these may be unidentifiable from restitution curve data[22]. There is a need for robust approaches that can interrogate cardiac tissue properties more thoroughly while at the same time minimising additional interventions.

In this paper, we present a novel method for probabilistic calibration of an electrophysiology simulator from spatially sparse measurements using a probabilistic model of electrophysiology parameters on a manifold representing the left atrium of the heart. We focus on estimating parameter fields that reconstruct heterogeneity in ERP. We chose to use a phenomenological EP model that captures the main features of cardiac activation and recovery. We determine two types of ERP measurements for calibrating EP parameters that determine excitability. EP parameters are modelled as latent Gaussian processes (GPs) on a manifold, and linked to observations via surrogate functions and a likelihood function designed for ERP measurements. We use Markov Chain Monte Carlo (MCMC) to obtain the posterior distribution of EP parameter fields across the atrium. We validate our method quantitatively by generating ground truths and calibrating to sparse data. The principles behind our method generalise to other measurement types, such as CV and APD restitution data, making our approach a step forward in the creation of digital twins capable of reproducing the complex dynamics of electrophysiology.

Results

Workflow

The computational model, or ‘simulator’, that we seek to calibrate is composed of (i) a finite element mesh representing an atrial manifold ; (ii) an electrophysiology model, that maps EP parameters 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: . See “Methods” for details on the simulation model, parameterisation, and allowable parameter ranges, and details on the numerical implementation. We use the software openCARP[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 from ERP measurements. Figure 1 represents our modeling workflow, which we summarize here. Our code is available in a Zenodo repository[26].

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 reduced-rank 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 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.

Surrogate functions

The simulator can be used to map parameter fields to ERP fields . 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 as , bypassing the need to run the simulator directly for inference.

Gaussian process priors

The mesh has 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 –. 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 . We use Gaussian Process Manifold Interpolation (GPMI), a method we proposed for defining Gaussian process (GP) distributions on manifolds[27]. The approach uses solutions 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 measurements, which accounts for how an S1S2 pacing protocol to determine ERP effectively measures the S2 interval in which ERP lies, rather than measuring ERP directly.

Sensitivity analysis and surrogate functions

Figure 2a shows sensitivity indices for two types of ERP: an ERP measurement for S1S2 with S1 600 ms, denoted here as , and another type of ERP measurement for S1S2S3 pacing for S1 600 ms S2 300 ms, denoted here as . 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 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. and , 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.

Figure 2

(a) Sensitivity analysis of ERP for the mMS electrophysiology model for all parameters. and have the largest sensitivity indices for both types of ERP. (b) Surrogate functions for predicting ERP from and . A discontinuity in occurs at  ms, so this region is shown in white.

The sensitivity indices in Fig. 2a show that these ERP measurements are mainly determined by and , 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 and . Figure 2b shows contour plots of the surrogate functions for ERP. A discontinuity occurs in the surface for parameters combinations resulting in  ms, so data for  ms were discarded before fitting this function. Note that the majority of clinical measurements fall in the range 170–270 ms, so even 280 ms could be considered as an upper limit[29].

Synthetic experiments

To test our methodology, we ran synthetic experiments as detailed in “Methods”. 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 and 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.

Figure 3

Ground truth and predicted (posterior mean and standard deviation) electrophysiology parameter fields. The spheres show the location of ERP measurements. All units are milliseconds.

Figure 4

Ground truth and predicted (posterior mean and standard deviation) ERP fields, and independent standard error (ISE) plots comparing predicted and ground truth ERP. The spheres show the location of ERP measurements. All units are milliseconds.

Figure 5

Ground truth and predicted APD values from simulation. All units are milliseconds.

The prediction of EP parameter fields and and ERP fields and 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 on the left atrial appendage is above 3, which may be caused by a combination of high ground-truth values for (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 generated and 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.

Figure 6

Quantitative validation for ERP prediction, for different lengthscales (1 unit kernel lengthscale  mm), number of observations, and S2 and S3 resolution. The model prediction is from the maximum a posteriori estimate. 32 basis functions are used for both EP parameter fields. For each combination of lengthscale, observations, and S2 and S3 resolution, 45 different samples of EP parameter fields are generated, and 5 different observation designs are used for each sample. The RMSE scores of ERP for each validation run is calculated over all mesh vertices. The error bars are 1 standard deviation.

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 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 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 is a normalised membrane voltage, h is a gating parameter that controls recovery, and is an externally applied stimulus. The 4 cell model parameters are time-constants that approximately characterize stages of the action potential sequence, and D is conductivity. We fixed the excitation threshold 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 0.1–1.5 m/s, 0.01–0.30 ms, 1–30 ms, 65–215 ms, 120–270 ms.

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 , the ERP value under an S1S2 protocol for S1 600 ms, and , 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. , as a function of a single model input, e.g. . 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

Surrogate functions

The map from EP parameters (inputs) to high dimensional tissue responses (outputs), such as restitution curves, has been modelled previously using Gaussian processes[21]. Here, cubic polynomials in both and were fit to corresponding values of and , generated from an optimized Latin hyper-cube design of 100 values of , and , keeping  ms and  ms in order to produce ERP and APD values in a range observed in human atrial tissue[29]. was varied for robustness, but has negligible effects on ERP, as confirmed by negligible fitting residuals. There is a discontinuity in for parameter values producing  ms, so data for  ms were discarded before fitting these functions. We refer to these polynomial fits as ‘surrogate functions’ for and , denoted as and respectively, as they allow for determining ERP without running simulations.

Gaussian process priors

We model the EP parameters fields, , as spatially correlated random fields defined on the atrial manifold, i.e. . We use Gaussian Process Manifold Interpolation (GPMI), a method we proposed for defining Gaussian process distributions on manifolds[27]. The approach uses solutions 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 can then be represented using the following probabilistic model, which uses the K smallest eigenvalue solutions to the Laplacian eigenproblem:with hyperparameters mean , amplitude , and lengthscale , for and . 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 must be given a Gaussian prior in order for this model to approximate a Gaussian process. The function is the spectral density corresponding to the choice of covariance kernel, with the square root of the eigenvalue 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 over the atrium, it is possible to calibrate the parameter fields and by obtaining the posterior distribution of the hyperparameters in Eq. (5). For convenience, we collect the hyperparameters into the vector . Defining the ERP measurements as , then we can write the Bayesian inference problem as:where and represent observations of and 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 lies between two S2 values and can be expressed in the following way (see Fig. 1 for a graphical representation):

Observations can be linked to the hyperparameters via the GP fields defined by Eq. (5), which determine the EP parameters at positions on the atrial mesh, and by the surrogate functions, which map these EP parameter values to ERP values:where and 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 , the true likelihood is given by a truncated uniform, or ‘top-hat’, distribution. In other words, is equal 1 if produces ERP in the specified interval, and 0 otherwise. The surrogate functions can be used to predict ERP from the EP parameters, which are determined by the GP fields 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 centered on each sub-interval with standard deviation equal to the sub-interval width . We choose , such that . For an observation (and similarly for ) 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 log-likelihood, 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 with (where is constant), which assists with defining priors for . We used the following priors on the hyperparameters: , , and . 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 and , and is set to values explained in Results and below. The generated samples for and are then scaled and offset into the full allowable parameters ranges. The same operation is performed for , which is needed for atrial simulations.

Certain combinations of and 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 . 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 where 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 and 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 field was generated for use in this simulation, and and 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 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 reached 0.7 on upstroke, and recovery (repolarization) was measured when fell to 0.8 (), 0.7, (), 0.5 (), and 0.1 (). APD values are the time between activation and recovery. Simulations results for and are given in Results.