Data-Driven Uncertainty Quantification for Cardiac Electrophysiological Models: Impact of Physiological Variability on Action Potential and Spiral Wave Dynamics

Computational modeling of cardiac electrophysiology (EP) has recently transitioned from a scientific research tool to clinical applications. To ensure reliability of clinical or regulatory decisions made using cardiac EP models, it is vital to evaluate the uncertainty in model predictions. Model predictions are uncertain because there is typically substantial uncertainty in model input parameters, due to measurement error or natural variability. While there has been much recent uncertainty quantification (UQ) research for cardiac EP models, all previous work has been limited by either: (i) considering uncertainty in only a subset of the full set of parameters; and/or (ii) assigning arbitrary variation to parameters (e.g., ±10 or 50% around mean value) rather than basing the parameter uncertainty on experimental data. In our recent work we overcame the first limitation by performing UQ and sensitivity analysis using a novel canine action potential model, allowing all parameters to be uncertain, but with arbitrary variation. Here, we address the second limitation by extending our previous work to use data-driven estimates of parameter uncertainty. Overall, we estimated uncertainty due to population variability in all parameters in five currents active during repolarization: inward potassium rectifier, transient outward potassium, L-type calcium, rapidly and slowly activating delayed potassium rectifier; 25 parameters in total (all model parameters except fast sodium current parameters). A variety of methods was used to estimate the variability in these parameters. We then propagated the uncertainties through the model to determine their impact on predictions of action potential shape, action potential duration (APD) prolongation due to drug block, and spiral wave dynamics. Parameter uncertainty had a significant effect on model predictions, especially L-type calcium current parameters. Correlation between physiological parameters was determined to play a role in physiological realism of action potentials. Surprisingly, even model outputs that were relative differences, specifically drug-induced APD prolongation, were heavily impacted by the underlying uncertainty. This is the first data-driven end-to-end UQ analysis in cardiac EP accounting for uncertainty in the vast majority of parameters, including first in tissue, and demonstrates how future UQ could be used to ensure model-based decisions are robust to all underlying parameter uncertainties.

Computational modeling of cardiac electrophysiology (EP) has recently transitioned from a scientific research tool to clinical applications. To ensure reliability of clinical or regulatory decisions made using cardiac EP models, it is vital to evaluate the uncertainty in model predictions. Model predictions are uncertain because there is typically substantial uncertainty in model input parameters, due to measurement error or natural variability. While there has been much recent uncertainty quantification (UQ) research for cardiac EP models, all previous work has been limited by either: (i) considering uncertainty in only a subset of the full set of parameters; and/or (ii) assigning arbitrary variation to parameters (e.g., ±10 or 50% around mean value) rather than basing the parameter uncertainty on experimental data. In our recent work we overcame the first limitation by performing UQ and sensitivity analysis using a novel canine action potential model, allowing all parameters to be uncertain, but with arbitrary variation. Here, we address the second limitation by extending our previous work to use data-driven estimates of parameter uncertainty. Overall, we estimated uncertainty due to population variability in all parameters in five currents active during repolarization: inward potassium rectifier, transient outward potassium, L-type calcium, rapidly and slowly activating delayed potassium rectifier; 25 parameters in total (all model parameters except fast sodium current parameters). A variety of methods was used to estimate the variability in these parameters. We then propagated the uncertainties through the model to determine their impact on predictions of action potential shape, action potential duration (APD) prolongation due to drug block, and spiral wave dynamics. Parameter uncertainty had a significant effect on model predictions, especially L-type calcium current parameters. Correlation between physiological parameters was determined to play a role in physiological realism of action potentials. Surprisingly, even model outputs that were relative differences, specifically drug-induced APD prolongation, were heavily

INTRODUCTION
Computational modeling of cardiac electrophysiology has long been used to generate hypotheses for experimental verification and investigate mechanisms underlying normal and pathological cardiac physiology. In recent years computational cardiac electrophysiological modeling has transitioned to clinical and regulatory applications. These include clinical trials that evaluate the ability of personalized whole-heart models to guide ablation therapy 1 , and the Comprehensive in vitro Pro-arrhythmia Assay (CiPA) program, which uses a computational model of the action potential (AP) as part of a framework for assessing drug cardiotoxicity (Strauss et al., 2018). These advances have coincided with a dramatic rise in the use of physics-based and physiological modeling in drug/device regulatory submissions supporting safety and efficacy/effectiveness of medical products, or in software within medical devices (Faris and Shuren, 2017;Morrison et al., 2018). In parallel, the medical devices community and other healthcare communities have focused attention on methods and best practices for ensuring the reliability of computational modeling approaches (ASME V&V 40, 2018).
Understanding the uncertainty in model predictions is acknowledged as a critical component of model credibility assessment (Oberkampf et al., 2004;National Research Council, 2012;ASME V&V 40, 2018). Just as measurement error is key to interpreting experimental measurements, supplementing computational model predictions with an estimate of uncertainty vastly improves the ability to make informed decisions. Uncertainty quantification (UQ) is the science of characterizing uncertainties in computational models. While this includes uncertainty in the model form, i.e., uncertainty in the best equations to model a system, the most common form of UQ involves two stages. First, characterizing the uncertainty in the model inputs-all the real-world measured quantities that are used in the model, such as model parameters, boundary conditions, and initial conditions. Second, "propagating" these uncertainties through the model to obtain the uncertainty in model outputs of interest. We refer to these two stages as uncertainty characterization and uncertainty propagation, respectively. With UQ inputs and outputs are generally represented using probability distributions, rather than taking fixed values.
For physiological models, the main reasons for input uncertainty are experimental measurement error and natural physiological variability. Both may be significant-typically orders of magnitude greater than measurement error and natural 1 https://clinicaltrials.gov/ct2/show/NCT03536052 variability within engineering systems. Physiological models are typically non-linear, so the impact of input uncertainty on model outputs is difficult to predict without simulation. Both these challenges apply to cardiac electrophysiological modeling. Accordingly, there has been a lot of recent interest in UQ within the cardiac modeling community, and many groups are overcoming the challenges. For a detailed discussion on challenges and previous work see our previous discussion in Pathmanathan et al. (2019) (introduction), and a recent special edition devoted to cardiac model UQ ; also see the Discussion section in the present paper which will compare the approach used herein with other prominent approaches.
Cardiac action potential models are typically sets of ordinary differential equations (ODEs). Several hundred models for different species have been developed (Noble et al., 2012); it is not uncommon for a model to have dozens of equations governing dynamic behavior of transmembrane voltage, membrane gating variables and intra-cellular and extracellular ionic concentrations, altogether involving hundreds of parameters. For this reason, all previous cardiac AP model UQ has been limited in one or both of the following: (i) only considering uncertainty in a restricted subset of the parameters, for example considering ionic current maximum conductances as uncertain but keeping parameters for gating dynamics fixed; and (ii) using arbitrarily chosen variation in the inputs, for example letting a parameter vary in a range, such as ±10 or 50% around its mean value, rather than estimating true uncertainty in the input due to measurement error or physiological variability. In our first work on this topic , we estimated uncertainty due to physiological variability in two parameters characterizing steady-state inactivation of the fast sodium current, and propagated this through the model to assess its impact on the action potential. This was an end-to-end UQ analysis, but only considered one component (steady-steady inactivation) of just one ion channel, i.e., was heavily limited by (i)-uncertainty in the vast majority of parameters was neglected. In our previous paper , we developed a novel canine action potential model which included representations of six major currents (fast sodium, inward rectifier, transient outward, L-type calcium, rapidly and slowly activating delayed rectifier), but had just seven state variables and 33 parameters (excluding environmental parameters). This allowed us to overcome limitation (i) and perform a comprehensive analysis where we accounted for uncertainty in all 33 parameters, including all conductances, steady state activation/inactivation parameters, and time constant parameters. This provided a wealth of information on the general robustness of the model. However, it was still limited by (ii), because we prescribed arbitrary input variation.
In this paper we expand upon our previous work  by using experimental data to estimate the parameter uncertainty, and compute the impact of the parameter uncertainty on the action potential and spiral wave dynamics. Using our canine model developed in Pathmanathan et al. (2019), uncertainty due to population variability was estimated for all parameters in the five currents active during repolarization (inward rectifier, transient outward, L-type calcium, rapidly and slowly activating delayed rectifier currents); 25 parameters in all. We did not attempt to estimate uncertainty in fast sodium current parameters due to the known challenges in obtaining high quality experimental data on this current under physiological conditions which makes even estimating average values of parameters challenging (Berecki et al., 2010), let alone estimating population variability. We used a variety of methods for estimating population variability. Where relevant data was available, we estimated correlation between parameters across the population. Overall, we generated probabilistic representation of all 25 repolarization parameters in the model. We then propagated this uncertainty through the action potential model, to compute the impact on: (a) the action potential shape; (b) action potential duration (APD) prolongation following drug-induced partial block of the I Kr current; and (c) spiral wave dynamics. We also investigated the role parameter correlation plays in governing AP shape. To our knowledge this is the first time such a data-driven end-to-end uncertainty quantification has been performed for more than a handful of AP model parameters, let alone for the vast majority, and the first such analysis in tissue. Our results provide for the first time information on the expected uncertainty in cardiac AP and spiral wave models given experimentally-derived parameter uncertainty in all repolarization parameters, and demonstrates how future UQ could be used to ensure cardiac model-based decisions are robust to all underlying parameter uncertainties.

Cell Model
The canine cell model of Pathmanathan et al. (2019) was used. In this model transmembrane voltage is governed by where V is transmembrane voltage, C m = 1µF cm −2 is the specific membrane capacitance per unit area, and I Na , I K1 , I to , I CaL , I Kr , and I Ks are ionic currents (respectively: rapid sodium, inward rectifier, transient outward, L-type calcium, rapidly and slowly activating delayed rectifier), and I stim is a stimulus current. The currents are formulated as where the g X are ion channel maximal conductances, m, h, z ∞ (V), where Y ∞ is the steady-state function and τ Y is the time "constant." Steady state functions for all gates are sigmoidal with − for activation gates and + for inactivation gates, where E Y is the half-activation/inactivation voltage for that gating variable and k Y > 0 controls the slope of the sigmoid. For time constants, voltage dependence was not modeled except for the h-gate: Overall, this cell model has 33 parameters, excluding E Na , E K , and E Ca . In Pathmanathan et al. (2019), "nominal" values of each of these parameters were defined, and the impact of arbitrarily chosen variation was analyzed. In this present paper, we studied the impact of data-driven uncertainty in all I K1 , I to , I CaL , I Kr , and I Ks parameters (25 parameters in total). For the other parameters (eight I Na parameters, and E Na , E K , and E Ca ) nominal values from Pathmanathan et al. (2019) were used. A CellML version of this model [using mean values of repolarization parameters (see Table 1 in section 3) and nominal values of I Na parameters and E Na , E K , and E Ca ] is available upon request.

Uncertainty Characterization
For the 25 repolarization parameters, our aim was to derive a first-estimate of population variability based on available data. Since this the first time such a UQ analysis has been performed on this scale for cardiac models, we consider relatively crude estimates acceptable for some parameters, as long as they are derived from experimental data. We assumed the type of distribution (normal or log-normal) and estimated distribution scale parameters (e.g., means and variances). Parameters which are physiologically constrained to be positive (conductances g X , slopes k X , time constants, τ * X ) were assumed to follow log-normal distributions. The other parameters (half activation/inactivation voltages E X ) were assumed to follow normal distributions. Different types of data was available for the different repolarization currents, and accordingly a variety of different methods was used to characterize the parameter uncertainty. These are discussed in the following sections.

I K1 Parameter Uncertainty
The inward rectifier current I K1 has three parameters: maximum conductance g K1 , half-inactivation voltage E z , and slope of steadystate inactivation at half-inactivation voltage k z . Previously, in Pathmanathan et al. (2019), we determined values for the three I K1 parameters by fitting to averaged voltage clamp recordings from canine epicardial cells (n = 12). Here, E exp K represents unknown effective E K in the experiments, which is dependent on the unknown intra and extracellular potassium ion concentrations. E exp z represents the effective E z in the experiments, and actual E z was computed by shifting E exp z to correspond to E K = −85 mV, i.e., E z = E exp z − (E exp K + 85). In essence, E K and E z were both fit and E z shifted to correspond to E K = −85 mV.
In the present paper, to calibrate with uncertainty, we fitted to recordings from the n = 12 individual cells. A "twostage approach" was used. Stage 1 was to apply the above process using each cell's I-V curve (rather than the averaged data), to obtain parameters for each cell. Stage 2 was to fit a multivariate probability distribution to the set of (g K1 , E z , k z ) values. Allowing for possible correlation between parameters across the population, we assumed a multivariate normal distribution and computed the mean vector and covariance matrix from the 12 samples.

I to Parameter Uncertainty
The transient outward current I to has six parameters: maximum conductance g to , half-activation voltage E r , activation slope k r , half-inactivation voltage E s , inactivation slope k s , and inactivation time constant τ * s . Previously, in Pathmanathan et al. (2019), nominal values for g to , E r and k r were determined by fitting to averaged voltage clamp recordings from canine epicardial cells (n = 16). E K was not identifiable from this dataset so was set to be −85 mV. In the present paper, to calibrate with uncertainty, the two stage approach was used as described in section 2.2.1, which again allows for possible correlation between these parameters. Inactivation parameters, E s and k s were determined by fitting F(E s , k s , α) = (1/(1 + exp((V − E s )/k s )) + α)/(1 + α) to voltage clamp recordings from canine epicardial cells (n = 14). The function is a sigmoid that decreases from 1 to α, and was chosen because the data exhibited an experimental artifact where peak current did not drop to zero at higher voltage. α was fit for each cell (to obtain more accurate E s , k s ) but not used in the model. The two-stage approach was again used.
To assess variability in τ * s , we first defined τ * s more precisely [than in our previous work ] to be the average value of voltage-dependent τ s (V) for voltages in the range 10-50 mV. This is consistent with the value of τ * s used in Pathmanathan et al. (2019), and is necessary to be able to meaningfully ask what the population variability in this quantity is (see Discussion in Pathmanathan et al., 2019). Using this definition, τ * s was computed for eight canine epicardial cells (raw data behind Figure 2 in Cordeiro et al., 2012). A log-normal distribution was assumed, and scale parameters µ, σ estimated from the eight samples.

I CaL Parameter Uncertainty
The L-type inward calcium current I CaL has six parameters:  (2012) reports means and standard errors for each of these values, so here we estimated uncertainty due to population variability by simply calculating standard deviations from the standard errors. E d and E f were assumed to be normally distributed with these means and standard deviations. k d and k f were assumed to be log-normally distributed, with parameters µ and σ 2 that were computed by inverting the relationship between (mean,variance) and (µ, σ ) for log-normal random variables. This provides relatively crude (compared to methods used for I K1 and I to ) but data-driven estimates of uncertainty for these four parameters, although this approach provides no indication of any correlations between them.
In Pathmanathan et al. (2019), the nominal value of g CaL was determined by finding the value of g CaL such that (5) with f = 1 and E d , k d at nominal values matched the experimental peak current of the I-V curve in Figure 3 (EPI) of Iyer et al. (2012) (using digitized data). Here, we first estimated experimental peak current standard deviation (n = 17 cells) from digitized peak current standard error in the same figure. We then determined (µ, σ ) values for log-normally distributed g CaL such that this variability in g CaL , together with the above variability in E d and k d , when propagated through (5), gave rise to the experimental mean and standard deviation of peak current. This approach was only possible because the variability in E d and k d was not sufficient to explain the experimental variability in peak current. Specifically, peak current was estimated to be −3.90 ± 1.40 mS/cm 2 (mean ± standard deviation) from Iyer et al. (2012) (Figure 3, EPI). With constant g CaL = 0.115 (nominal value from Pathmanathan et al., 2019), propagating the above uncertainty in E d and k d through (5) leads to simulated peak current of −4.12 ± 0.41, less variability than in the experiments. Hence it is possible to match the experimental mean and standard deviation by introducing variability in g CaL .
In Pathmanathan et al. (2019), the nominal value of τ * f was chosen based on Figure 8 of Xiao et al. (2006). Here, we first re-defined τ * f as the average value of τ f (V) over voltages between 0 and 40 mV (to enable us to meaningfully ask what is the population variability in τ * f ). We then digitized the male and female means and standard errors in Figure 8 of Xiao et al. (2006), computed standard errors from standard errors, computed grouped (both sexes) means and standard deviation at each voltage, and then computed the mean and standard deviation of τ * f using the above definition.

I Kr Parameter Uncertainty
The rapidly activating delayed rectifier current I Kr has six parameters: maximum conductance g Kr , half-activation voltage E xr , activation slope k xr , activation time constant τ * xr , halfinactivation voltage E y , and inactivation slope k y .
In Pathmanathan et al. (2019) E xr , k xr , E y , and E y were taken from Table 1 of Berecki et al. (2005). Here we used standard errors reported in that table to estimate uncertainty in these parameters, using the same method as described above for I CaL parameters in section 2.2.3.
The parameters g Kr , τ * xr , together with I Ks conductance g Ks (below), are fundamentally different to all other parameters in this model. In Pathmanathan et al. (2019) all other parameters had nominal values derived directly from voltage clamp data or data in the literature. These three parameters however were calibrated using the full action potential model, to match experimental restitution data. Therefore, an appropriate method for quantifying the variability in these parameters is not obvious. In principle, we could estimate population variability in the restitution curve, and determine probability distributions for g Kr , τ * xr , and g Ks that give rise to that variability in restitution. However, there are two problems with such an approach. First, it is inconsistent with our aim of studying the impact of parameter uncertainty on action potential characteristics, since AP characteristics would be used to determine the parameter uncertainty. More importantly, to be fully rigorous, the quantified uncertainty in all the other parameters should be accounted for (in the same way uncertainty in E d and k d was included when calibrating g CaL to variability in peak I CaL in section 2.2.3). But we already know from the results in Pathmanathan et al. (2019) that moderate levels of uncertainty lead to very wide ranging APDs and action potentials that do not repolarize, and therefore the uncertainty in simulated restitution curve due to variability in the other parameters, will be greater than the experimentally observed variability in restitution.
Therefore, we used the following approach. In Table 1 (see below),σ LN is defined as the ratio of standard deviation to mean, for log-normally distributed parameters. Noting that the other conductances and time constants haveσ LN values of 28, 39, 13, 35, 31, and 22%, we set mean values for these three parameters to be their nominal value (from Pathmanathan et al., 2019) and chose standard deviations corresponding toσ LN = 40%, a conservatively wide-ranging distribution. Thus, these parameters have distributions that are not arbitrarily chosen, but only based on information regarding other parameters.

I Ks Parameter Uncertainty
The slowly activating delayed rectifier current I Ks has four parameters: maximum conductance g Ks , half-activation voltage E xs , activation slope k xs , and activation time constant τ * xs . In Pathmanathan et al. (2019) nominal values of E xs , k xs , and τ * xs uncertainty were taken from Liu and Antzelevitch (1995) (text); here we set mean values as the nominal values and computed standard deviations from reported standard errors using values in Liu and Antzelevitch (1995). For g Ks , see section 2.2.4.

Comparing Parameter Uncertainty
Once uncertainty in model parameters has been quantified, and before the impact of that uncertainty on model outputs was investigated, we wished to directly compare the uncertainty in parameters. To do so, we introduced two relative measures of uncertainty. For half activation/inactivation voltages (normallydistributed), we definedσ N as standard deviation divided by a representative range of 100mV. For all other parameters (all constrained to be positive, all log-normally distributed), we definedσ LN as standard deviation divided by mean. Bothσ N and σ LN are analogous toσ in Pathmanathan et al. (2019).

Single Cell and Tissue Simulations
Single cell simulations were run using Chaste, a general purpose package for physiological simulations (Mirams et al., 2013). ODEs were solved in Chaste using the CVODE adaptive ODE solver (Hindmarsh et al., 2005). Initial conditions were dependent on parameter values (see Pathmanathan et al., 2019). Electrical activity in tissue was simulated using the monodomain equations coupled to (2)- (7), where χ = 1, 400cm −1 is the surfacearea-to-volume ratio, and σ = 1.4 mS/cm is the bulk conductivity. Tissue simulations were performed with a GPU-based finite difference solver implemented in WebGL 2.0 using the Abubu.js library (Kaboudian et al., 2019a,b). The WebGL 2.0 implementation of cardiac models provide significant speedup and better performance compared to other GPU implementations, such as NVIDIA CUDA or OpenACC (Kaboudian et al., 2019c). Simulations were performed on a 16 by 16 cm square domain, with initial conditions that were the same as single cell simulations except a gradient in V was prescribed in the x-direction (linear, −85 mV to −5 mV) and a gradient in h was prescribed in the y-direction (linear, ranging from h = 0.1 to 0.6). With these initial conditions a spiral wave spontaneously forms at t = 0, consistently across the majority of the parameter values used. We quantified the impact of the parameter uncertainty on two quantities of interest: (i) number of phase singularities [PSs; centers of reentrant waves, computed as described previously Galappaththige et al., 2019)] at times t = 1, 000, 1, 500, and 2, 000 ms; and (ii) for the node at the center of the domain, average cycle length for all full beats in the window t = 1, 000 to t = 2, 000 ms. (Note that this is an average over time for a single parameter set, not an average over parameter space).

Uncertainty Propagation, Sensitivity Analysis, Non-behavioral Analysis
The same analytic methods as in our previous paper were used; see Pathmanathan et al. (2019) for a full description. Briefly, simple Monte Carlo sampling was used for uncertainty propagation. For single cell simulations we used N = 100, 000 and for tissue simulations we used N = 1, 000 due to increased computational cost. All simulation results are used for computing histograms, but when plotting representative simulations we plotted 1,000 random selected action potential traces for single cell simulations and 10 for tissue, and 100 randomly selected snapshots of activity in tissue. For global sensitivity analysis (determining how much the uncertainty in an output can be apportioned to the underlying uncertainties in each input), variance-based sensitivity indices were computed using the Saltelli sampling method and the SALib library (Herman and Usher, 2017) (N = 10, 000 per uncertain parameter). Non-behavioral analysis was used to identify which parameters were responsible for different classes of behavior observed (e.g., determining which parameters where highly influential in AP repolarization failure, when this occurred). Cumulative distribution functions (CDFs) of "parameter value given behavior occurred" and "parameter value given behavior did not occur" were computed and the Kolmogorov-Smirnov test was used to test if the CDFs were statistically different, which indicates that the parameter is influential in whether the behavior occurs or not. Highly influential parameters were defined as those for which the test statistic was >0.2. Figure 1 plots the individualized fits of (8) to I K1 voltage clamp data recorded from 12 cells. Each color represents a different cell. Each point in parameter space represents the fitted values for that cell. Also plotted as ellipses are 90% confidence regions for each fitted parameter, using the variance-covariance matrixσ (J T T) −1 , whereσ is the estimated residual variance and J is the Jacobian matrix with entries J j = ∂f ∂p j . For most cells the calibration uncertainty (the ellipses) is small, but somewhat surprisingly for a few cells the calibration uncertainty was similar in size to the population variability across the cells (range of points). Note that confidence regions are presented here for illustration and not included in subsequent analysis. See discussion of limitations in section 4. See Supplementary Figures 1, 2, for the corresponding results for I to activation and inactivation parameters. Some correlation was observed between I K1 parameters, between I to activation parameters, and to a lesser extent between I to inactivation parameters (see trends in fitted parameter values in Figure 1C, Supplementary Figures 1, 2).

Uncertainty in Parameters
The final probability distributions that were derived for each parameter are provided in the Appendix, and summarized in Table 1. The non-zero off-diagonal elements of the covariance matrices for I K1 and I to parameters (Appendix) indicate the observed correlation. Looking first atσ N , it is remarkable that all half-activation/inactivation values except E xs have similar levels of uncertainty: each hasσ N in the 4-8% range. This is despite a variety of methods and data sources being used to obtain these values. Uncertainty in E xs is much larger (σ N = 34%). This value was derived directly from the large standard error for E xs reported in Liu and Antzelevitch (1995), and potentially dominated by measurement error rather than true variability. Although the experimental data on I Ks activation in canine epicardial cells are sparse, a more recent study (Obreztchikova et al., 2006) exhibited variation less than (Liu and Antzelevitch, 1995) and more comparable with the other parameters. Uncertainty in the five conductances g X was in theσ LN = 28-40% range (two chosen to be 40%). Different data and methods were used to determine the variability of g K1 , g to , and g CaL ; it is interesting that the results are similar in magnitude. There was less uncertainty in slopes k X , all but one around 15%. Time constant uncertainty was in the 20-40% range. All values are considerably higher than the 1, 3, and 5% we prescribed in our previous work . Figure 2 illustrates the impact of the quantified parameter uncertainty on the predicted action potential. Figure 2A plots 1,000 sample APs, with all 25 parameters sampled from the probability distributions provided in the Appendix. Unsurprisingly, given the wide range of APs observed in our previous work with justσ = 5%, a very wide range of APs was observed. Figures 2B-F plot sample APs with all parameters in each current uncertain but all other parameters fixed at mean values, with APD histograms (N = 100, 000, requiring about 45s of computational time on a desktop using 20 CPUs) and sensitivity indices in insets. Sensitivity indices were not computed for the I K1 and I to results because for these currents there is parameter correlation, and methodology for calculating sensitivity indices given correlated parameters is not well-established (Most, 2012). The uncertainty in I CaL parameters has the greatest impact on the action potential, with early-after depolarizations (EADs) and repolarization failure occurring in some cases. No such behavior is observed for any other current. Uncertainty in I Kr and I Ks significantly affects APD. For I Ks the sensitivity analysis shows that this is almost entirely due to the uncertainty in E xs , which as discussed in section 3.1 is most likely dominated by measurement error. Surprisingly, the distribution of APD in the I Ks results is bi-modal. To investigate why, we plotted APD as a function of E xs , both with other parameters held fixed and with other parameters varying (see Supplementary Figure 3). The figure shows that E xs takes values in the range −75 to 150 mV (perhaps wider-ranging than that which occurs in reality as discussed above). There is a non-linear relationship between APD and E xs , but at greater values of E xs (e.g., >75 mV) the relationship is nearly flat (APD always near 450 ms), which is expected because at these values of E xs I Ks will not activate and therefore changes in E xs will not impact APD. This explains the clustering of APD values near this value, hence the second peak in the histogram in Figure 2F.

Impact of Parameter Uncertainty on Action Potential
Uncertainty in I K1 parameters has the least impact on AP shape. Uncertainty in I to parameters has significant effect on phase 1 of the AP, but with large probability (p = 0.97) does not heavily affect APD. There is a small probability (p = 0.03) of AP shortening (APD < 200 ms).
Non-behavioral analysis (see section 2.4) was performed to determine parameters responsible for different behaviors. I CaL APs (N = 100, 000) were generated and categorized as being either: normal; exhibiting loss of dome morphology; exhibiting EADs; or repolarization failure. I to APs were categorized as either being normal or exhibiting reduced AP, delineated based on whether the AP exhibited spike-and-dome morphology or not (see Figure 2C). The results are provided in Table 2. Table 2 shows that variability in E d and E f (steady state halfactivation and half-inactivation voltages for the d and f gates in the I CaL current) are influential in determining EADs or repolarization failure. Previously in Pathmanathan et al. (2019), we came to the same conclusion using arbitrary variation, and we hypothesized that these parameters may be correlated in reality. The probability distributions derived for E d and E f (see Appendix) assume no correlation, because we had no data allowing us to infer any correlation (section 2.2.3). In other words, the covariance matrix for (E d , E f ) is diagonal. To assess if missing parameter correlation could be responsible for the wide range of APs observed under I CaL uncertainty, we introduced correlation between E d and E f at a range of different levels, and computed the probability of loss of dome morphology, EADs, or repolarization failure. Different values of correlation coefficient r were chosen, spanning r = 0 (no correlation, i.e., above results) to r = 1 (perfect correlation). Off-diagonal elements of the (E d , E f ) covariance matrix were set to covariance values corresponding to the chosen r. The results are plotted in Figure 3. All behaviors became less likely with increasing E d -E f correlation, and the probability of EADs or repolarization failure approaches zero at perfect correlation. EADs and repolarization failure are likely related to the I CaL window current, which is directly impacted by the values of E d and E f . This raises the question of whether correlation between E d and E f impacts the magnitude of the window current, or location, or both. We computed the distribution of window current magnitude and location, for different values of r, and determined that correlation does not affect I CaL window current location but it does affect window current magnitude (see Supplementary Figure 4). Specifically, correlation between E d and E f reduces the probability of the window current having larger magnitude, and through this reduces the probability of EADs.

Impact of Parameter Correlation
For I K1 and I to , our probability distributions include parameter correlation (see covariance matrices in Appendix). To further investigate the importance of parameter correlation, we removed these correlations (i.e., made the covariance matrices diagonal), and re-computed APs. Results are shown in Figure 4A for I K1 and Supplementary Figure 5 for I to . Loss of correlation led to slightly wider-ranging APs; interestingly loss of I K1 parameter correlation led to some non-physiological APs (see arrow in Figure 4A). To understand why, we visualized the parameters values that led to the non-physiological APs in parameter space (see Figure 4B). The non-physiological APs were associated with the region of parameter space where E z and k z both took low values, and additionally g K1 was in the lower half of its set of values. In the original probability distribution, before correlation was removed, E z and k z were somewhat negatively correlated, which meant there was very small probability of sampling the region where E z and k z both took low values. When correlation was removed, the probability of sampling this region became non-negligible. E z and k z both taking low values, together with smaller g K1 , corresponds to a small or near-zero I K1 current which gives rise to the non-physiological APs.
Overall, both correlation experiments (Figures 3, 4) suggest parameter correlation could play an important role in ensuring simulated APs are physiological when parameter variability is modeled.

Impact of Parameter Uncertainty on Drug-Induced Action Potential Duration Prolongation
Often, cardiac models are not used to predict absolute quantities. Many times users are only interested in relative differences, for example changes in simulation results under drug block or when a proposed therapy is applied. While we observed large variability in predicted APs in Figure 2, we hypothesized that relative differences would be much more robust (less variable) to the underlying parameter uncertainty. To test this hypothesis, we computed the uncertainty in APD prolongation to partial block of the I Kr current. Figure 5A illustrates the control action potential (no block) and the action potential under 50% block of I Kr , without including impact of parameter uncertainty. APD prolongation was 70.3 ms. Figure 5B then demonstrates the impact of I K1 variability (only). The control uncertainty is the same as shown in Figure 2B. Surprisingly, APD, a relative quantity, is heavily affected by underlying uncertainty, and its total uncertainty is similar in magnitude to the control uncertainty. This is because the APDs under partial I Kr block (red traces/histogram) are very wide ranging. Figure 6 repeats this analysis for the full range of I Kr block (0-100%), and with variability in all five currents. Parameter uncertainty always has a significant affect on results, especially at greater levels of block.

Impact of Parameter Uncertainty on Spiral Wave Dynamics
Finally, we investigated the impact of the parameter uncertainty on reentrant spiral wave simulations in tissue (2D, 16 by 16 cm square domain). The results of the baseline simulations are shown in Figure 7, which uses mean values of all parameters. Two seconds of activity were simulated, taking about 40 s of computation using a desktop with a GeForce GTX 960 graphics card with 1,024 CUDA cores (typical utilization: 50-70%). We then considered each current in turn, randomly sampling N = 1, 000 parameters for each, with parameters in other currents kept at mean values. Results are presented in Figures 8-10. Figure 8 presents 100 snapshots of the activity at t = 1, 000 ms, for 100 of the 1,000 randomly chosen I to parameter sets. Figure 9 is the equivalent figure for I CaL parameter uncertainty. Corresponding figures for I K1 , I Kr , and I Ks , and for all currents at t = 2, 000 ms, are provided in the Supplementary Material. Figure 10 presents the distributions of number of phase singularities (PSs) and average cycle length for the node at the center of the domain. This figure also includes 10 random traces for the center node, as a simple visual indication of the output variability. First, it is interesting that across all the 5,000 simulations, it was possible to induce reentrant activity lasting at least 1,000 ms in about 80% of the simulations, all using the same induction protocol. This observation is based on the number of simulations with one or more phase singularities at t = 1, 000 ms, given that zero phase singularities implies no activity (blue subplots in Figures 8, 9) or activity about to die out. This demonstrates that spiral wave inducibility is reasonably robust to the underlying parameter uncertainty, perhaps more so than might have been expected.
Considering first I K1 and I to , uncertainty in these parameters had very little impact on inducibility, and limited impact on average cycle length. Figure 8 shows that most of the I to simulations had a similar voltage profile after 1,000 ms as the baseline simulation, though there was a small probability of low wavelength circular spirals. We observed in Figure 2C that some random I to parameters led to shortened action potentialsit is these parameters which lead to the circular spirals in the 2D simulations. For both I K1 and I to the distributions of average cycle length were similar to the single cell APD distributions (Figures 2B,C), though for I K1 less skewed and for I to more wide-ranging. I CaL uncertainty had the greatest impact on the tissue results, with near 50% of simulations not inducing or dying out by 1,000 ms, and with a bimodal wide-ranging distribution of average cycle length for the remaining simulations. As with I to , parameters corresponding to shortened action potentials led to circular low-wavelength spirals. Interestingly, highly complex behavior occurred more often under I CaL uncertainty than for the other currents: 3% of I CaL simulations had >10 PSs at t = 2, 000 ms, in contrast to <1% for the other currents (results not plotted). Overall, I CaL parameter uncertainty appears to strongly affect the type of reentrant behavior; the results again emphasize the importance of this current on model results and the need to reduce or at least better characterize uncertainty in its parameters. Finally, I Kr and I Ks had a moderate impact on the tissue simulations: uncertainty in these parameters had less of an effect than I CaL parameters but more of an effect than I K1 and I to parameters, for example activity had died out after 2 s in about half I Kr and I Ks simulations.

Summary of Results
In general data-driven, comprehensive (in the sense of considering all parameters) UQ in cardiac modeling has been considered infeasible, due to the complexity of cardiac cellular models. All previous work in this domain has either limited the uncertain parameters to a subset (usually minority subset) of the total number of parameters, and/or (usually and) prescribed the parameter uncertainty rather than estimating it from experimental data. The aim of this work was to quantify the uncertainty, due to population variability, in repolarization parameters of a canine cardiac action potential model, and to understand the subsequent uncertainty in model predictions. Model outputs that were studied included the action potential, relative effects on action potential duration of drug-induced current block, and "biomarkers" of reentrant behavior in tissue simulations. Our results demonstrate for the first time the feasibility of comprehensive data-driven cardiac UQ (though we only accounted for uncertainty in repolarization parameters; I Na parameters were left fixed).
A variety of approaches were used to characterize uncertainty in the repolarization model parameters: (i) fitting parameters to data derived from a number of individual cells; (ii) back-calculating standard deviations from reported parameter standard errors; (iii) fitting parameters with uncertainty, to observable data (e.g., peak current) given or after backcalculating uncertainty in those observables; and (iv) specifying a conservative estimate of uncertainty based on uncertainty in related parameters for other currents. Interestingly, the different methods (i)-(iii) provided similar estimates of uncertainty for related parameters (i.e., steady-state activation/inactivation voltages; time constants; conductances). The same is true of method (iv) by definition. The ensuing parameter uncertainty was substantial, and greater than that typically prescribed when arbitrary variation is used [including our previous work ]. When the parameter uncertainty was propagated through the AP model, there was  substantial uncertainty in predicted action potentials, which was not surprising given our previous results in Pathmanathan et al. (2019) where even moderate prescribed parameter uncertainty led to significant output uncertainty. I K1 parameter uncertainty had least impact on the AP; followed by I Kr . For both these currents, AP shape was relatively robust (invariant) to the underlying uncertainty and no anomalous APs arose. I Ks parameter uncertainty led to surprisingly large output uncertainty, although this was caused by large uncertainty in one parameter that was most likely dominated by experimental measurement error rather than genuine population variability (see section 3.1). When considering uncertainty in I to parameters, there was a large probability of little impact on AP shape, and a small probability of loss of spike and dome morphology and shortened APs. Uncertainty in I CaL parameters strongly affected APs, including loss of spike and dome morphology and shortened APs, early-after-depolarizations and repolarization failure ( Figure 2D). These results suggest that when experimental resources are limited, focused experiments to better characterize I CaL parameter variability may be the best use of resources.
We expected that relative differences, such as change in APD under drug-induced current block, would be more robust to the underlying uncertainty, even if absolute values of model outputs, such as APD, exhibited large uncertainty. Surprisingly, this was not the case; the former exhibited similar magnitudes of uncertainty as the latter. While this is an unfortunate result, analyses such as that performed in Figure 6 could be used to identify level of block at which response is relatively robust to underlying uncertainty. For example, 50% appears a reasonable choice based on Figure 6. Additionally, Figure 6 reveals the currents (I CaL and I to ) that it may be important to focus experimental investigation when attempting to minimize output uncertainty.
We had limited data available for building parameter correlation into our statistical model; covariance was included for I K1 and some I to parameters only. However, we investigated the importance of correlation, by removing correlation from I K1 or I to , and by artificially introducing it to two I CaL parameters. We observed a clear impact of the correlation: I K1 correlation was responsible for the model exhibiting realistic AP shape, and the introduction of hypothesized I CaL correlation reduced the probability of anomalous APs occurring. This strongly supports the expectation that parameter correlation exists in reality. Moreover, the results suggest correlation will FIGURE 6 | Impact of parameter variability on predictions of APD prolongation under block of I Kr . (A) APD prolongation as a function of I Kr block, with no variability accounted for. (B-F) Corresponding results considered variability in specific current, with parameters in other currents held fixed. Black lines are means, error bars are ±1 standard deviation, histograms show distribution of ( APD | APD computable). In some cases APD was not computable ("failure"), because either control or block AP did not repolarize. p F represents probability of such failures. All histograms use 20 bins. often need to be accounted for in cardiac electrophysiological modeling if population variability is modeled. Note that we only considered intra-current parameter correlation, but there is also evidence of inter-current parameter correlation, especially in maximal conductances due to co-expression of currents. For example, Milstein et al. (2012) shows potential correlation between g Na and g K1 in rat, Rees et al. (2018) presents evidence for correlation between g CaL and potassium current conductances in mouse, and recently Ballouz et al. (2019) analyzed RNA-Seq datasets to identify correlation between g CaL and g Kr in human. Fully characterizing such relationships will likely be one of the biggest challenges in cardiac AP model UQ.
Note that our ability to draw direct conclusions about human AP models based on these results is limited. While canine is frequently used as an experimental model for human, there are known differences which limit the ability to extrapolate even average behaviors from canine to human, for example related to I to characteristics (Akar et al., 2004) or smaller I K1 , I Kr , and I Ks in human vs. dog which results in a lower repolarization reserve (Jost et al., 2013). Species differences in regards to variability are even more uncertain, and it is difficult to estimate the extent to which the conclusions in this work (importance of I CaL variability on action potential, potential importance of parameter correlation for physiological behavior, and uncertainty in relative model outputs) will apply to human action potential models. Further research is needed to verify if the findings apply to human models, but the methods and workflow presented in this paper provide a FIGURE 7 | Baseline spiral wave simulation using mean values of all parameters. A rotating spiral wave was induced on a square two-dimensional sheet by setting gradients in voltage and h at t = 0. (Top) Progression of the spiral wave at selected times (color represents transmembrane voltage). Number of phase singularities in snapshots: 40 ms: one; 350 ms: one; 670 ms: one; 1,000 ms: three; 1,500 ms: one; 2,000 ms: three. (Bottom) Transmembrane voltage as a function of time for the node at the center of the tissue.
pathway to begin such research. Our results do suggest that appealing to the fact that a quantity of interest is only a relative difference is not a particularly strong justification for not performing UQ (including for human models), in the absence of other evidence.

Comparison With Other Work
It is worth comparing our approach with related cardiac UQ efforts. In this sub-section we focus on the method used to to characterize parameter variability; see Pathmanathan et al. (2019) for a comparison with other publications in regards to the number of parameters varied. As stated earlier, it is common to introduce arbitrarily-chosen uncertainty in parameters when studying impact of parameter uncertainty on cardiac AP model predictions (e.g., Sadrieh et al., 2014;Chang et al., 2015;Hu et al., 2018;Ballouz et al., 2019). We are not aware of other analyses that propagate entirely data-driven (not arbitrarilychosen) independently-derived parameter uncertainty though a cardiac AP model, apart from our initial work on this subject for I Na inactivation , and articles that consider uncertainty in non-endogenous parameters [e.g., drug-response parameters (Chang et al., 2017;Costabal et al., 2019)]. There are also methods that use the full action potential model for parameter uncertainty estimation, such as traditional Bayesian calibration methods (Houston et al., 2020), historymatching approaches (Coveney and Clayton, 2018), and the recent population of models (POM) approach (Britton et al., 2013;Gong et al., 2020), all of which are related (Whittaker et al., 2020). Bayesian calibration methods compute a posterior distribution for parameter values. The interpretation of this distribution depends on the type of data the model is fit to. If the model is fit to single values of an output quantity (e.g., APD90), FIGURE 8 | Snapshots of activity at t = 1, 000 ms for 100 random parameters, under I to parameter uncertainty. the posterior distribution represents calibration uncertainty (similar to the ellipses in Figure 1C) (see e.g., Houston et al., 2020). If the model is fit to a distribution of the output quantity, based on information on how the quantity varies across a population, the posterior distribution represents parameter variability (and calibration uncertainty). The POM approach is similar, and involves specifying bounds on quantities of interest that are considered physiologically reasonable, and retaining points in parameter space, from a pre-selected region, for which the model output falls within the pre-specified range. This has recently been extended to specifying output distributions instead of just bounds, and computing a multivariate parameter distribution (Lawson et al., 2018), which could be considered analogous to the parameter distributions we derived (Appendix).
The key difference between the POM approach and our approach used here is that we aimed to avoid using the action potential model to quantify parameter uncertainty-all parameter uncertainty estimates were based on independent data or comparison to related parameters (section 2.2.4). This was because our focus was model evaluation, rather than model application; specifically we wanted to understand the impact of parameter uncertainty on the model, which would not be possible had we used to the model to compute parameter uncertainty. The POM approach does not test the model, and in fact the reliability of the parameter distributions obtained is founded upon the reliability of the model. Our approach is especially a "bottomup" forward UQ analysis, whereas the Bayesian calibration/POM approach is a backward UQ which has the potential for inferring true parameter distributions (especially for parameters which are difficult to measure experimentally) and correlations between parameters (especially inter-current correlations), but is heavily reliant on confidence in the model equations. Ultimately, both approaches have value in different stages of the model development and application workflow. One possibility is using the bottom-up forward UQ approach in iterative model development and experimental resource allocation, with the resultant probability distributions used to inform the region of parameter space to be considered in a subsequent backward UQ stage.

Limitations
There are various limitations to the methods applied in this paper. First, we used relatively crude methods to estimate population variability for several parameters, such as backcalculating standard deviation from reported standard errors, and we assumed the distribution shape (normal or log-normal) rather than inferring from the data, or treating the uncertainty as epistemic and using interval analysis (Oberkampf et al., 2004). We consider these reasonable choices given the scarcity of relevant data, because this paper is the first time such an analysis has been attempted and there were dozens of parameters that needed to be considered. It is plausible however that many of distributions we provided in the Appendix are wider-ranging than in reality, and there is considerable potential for refinement with focused experimentation and more sophisticated analysis. Moreover, we neglected the impact of experimental bias and measurement error. Experimental variability may be responsible for some or much of the observed variability (Lei et al., 2020a). A second limitation is the simple two-step method used to estimate the population mean and covariance matrices, for currents (I K1 and I to ) where cell-wise data was available. This approach does not account for calibration uncertainty (ellipses in Figure 1C. A better approach is to use a hierarchical method, such as non-linear mixed effects which would essentially down-weight cells for which there is substantial calibration uncertainty (large ellipses) when inferring the population mean and variance. Another limitation is the fact that we considered each current in isolation when propagating the uncertainty through the model. This is at odds with the global sensitivity analysis paradigm in which the full parameter space given the specified probability distributions should be explored. We could have performed a global sensitivity analysis of the results of Figure 2A (where all parameters were varied), and then assessed the contributions of each of the different currents to the range of APs observed (non-behavioral analysis would also have been needed, since repolarization failure occurred; see Pathmanathan et al., 2019). However, we believe it is more intuitive and useful in an initial analysis (and therefore appropriate for this paper) to consider FIGURE 10 | Impact of parameter uncertainty on the spiral wave results. Left column of figures: 10 sample traces of transmembrane voltage for the center node, using random parameters for the current under consideration. Middle column: distribution of number of phase singularities (0, 1, >1) at t = 1, 000, 1, 500, and 2, 000 ms (labeled "1s" etc), given the parameter uncertainty. Right column: distribution of average cycle length in the t = 1, 000 to t = 2, 000 ms window, conditional on at least 1 cycle in this window. each of the currents in isolation. Finally, it should be noted that the simplicity of the model impacts the type of variability observed, for example there was little possibility of significant variability in resting membrane potential because the model did not include any exchanger or background currents and because we did not allow environmental parameters (E K , E Na , E Ca ) to vary.

Outlook
The ultimate goal of this line of research is the development of AP models with parameter uncertainty quantified in the form of probability distributions, rather than parameters taking fixed values, ideally together with the model and parameter distributions being such that the resultant output uncertainty matches that observed in reality-a very strong level of validation.
Achieving this admittedly ambitious goal will require refinement of simplified models, such as that used here, and more accurate characterizations of parameter uncertainty (see limitations above; our initial estimates of parameter uncertainty are quite crude). It is also likely that "model discrepancy"-the difference between model equations and reality due from limitations in (or incorrect) model form-will need to be characterized using formal model discrepancy techniques (Lei et al., 2020b), since simplified models for which comprehensive UQ is feasible generally exhibit greater discrepancy (Huberts et al., 2018). It is worth noting that these approaches are relevant to patient-specific models too. With rare exceptions, patient-specific models contain parameters that are not personalized to patient data, in addition to the personalized parameters (Gray and Pathmanathan, 2018). The non-personalized parameters are usually fixed to population average values, but in principle it is better, and more aligned to the UQ paradigm, to characterize them using probability distributions that account for population variability (ideally, probability distributions conditional upon the available patient data, such as sex or age or covariates that are personalized). Statistical representations of population variability will be useful even for parameters that will be personalized to patients, because those distributions can serve as informative priors in Bayesian calibration.

Conclusion
Overall, our results provide for the first time information on the expected uncertainty in cardiac AP and spiral wave models, given experimentally-derived parameter uncertainty in all repolarization parameters. We have demonstrated that robustness of model outputs to parameter uncertainty should not be assumed, even for model outputs that are relative differences. Our results begin to reveal the role of correlation between parameters in the cardiac action potential. We have demonstrated that it is feasible to perform data-driven forward uncertainty quantification in cell model parameters considering uncertainty in the vast majority of parameters, including assessment of the impact of uncertainty in computationallydemanding tissue simulations. We believe that use of the methods and approaches presented here will lead to improved models and ultimately more reliable decision-making, when cardiac models are used in clinical, regulatory, or product development applications.

DATA AVAILABILITY STATEMENT
All data used in model generation is already published (see text), except that used for parameterizing the I K1 model (see Table 1). This data is available upon result.

DISCLOSURE
The mention of commercial products, their sources, or their use in connection with the material reported herein is not to be construed as either an actual or implied endorsement of such products by the Department of Health and Human Services.

APPENDIX
In the below conductances g X have units mS/cm 2 , half (in)activation voltages E X have units mV, slopes k X have units mV, and time constants τ * X have units ms. For I K1 :