^{1}U.S. Food and Drug Administration, Center for Devices and Radiological Health, Silver Spring, MD, United States^{2}Department of Experimental Cardiology, Masonic Medical Research Institute, Utica, NY, United States^{3}School of Physics, Georgia Institute of Technology, Atlanta, GA, United States

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.

## 1. 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 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 (Mirams et al., 2020); 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 extra-cellular 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 (Pathmanathan et al., 2015), 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 (Pathmanathan et al., 2019), 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 (Pathmanathan et al., 2019) 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.

## 2. Methods

### 2.1. 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\mu $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*), *r*_{∞}(*V*), *s*, *d*_{∞}(*V*), *f*, *x*_{r}, *y*_{∞}(*V*), and *x*_{s} are probability of gates being open [activation gates: *m*, *r*_{∞}, *d*_{∞}, *x*_{r}, *x*_{s}; inactivation gates: *h*, *z*_{∞}(*V*), *s*, *f*, *y*_{∞}(*V*)], and *E*_{Na}, *E*_{K}, *E*_{Ca} are the Nernst potentials for sodium, potassium, and calcium, respectively. Per Equations (2)–(7), gates *z*, *r*, *d*, and *y* are not state variables (not governed by an ODE); these gates are taken to instantaneously reach steady state. Gating variables, *m, h, s, f, x*_{r}, *x*_{s} have dynamics governed by the ODE

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:

where τ_{h0}, δ_{h}, ${\tau}_{m}^{*}$, ${\tau}_{s}^{*}$, ${\tau}_{f}^{*}$, ${\tau}_{xr}^{*}$, ${\tau}_{xs}^{*}$ are positive constants. 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} (Pathmanathan et al., 2019)] is available upon request.

### 2.2. 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, ${\tau}_{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.

#### 2.2.1. *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 steady-state 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}_{\text{K}}^{\text{exp}}$ represents unknown effective *E*_{K} in the experiments, which is dependent on the unknown intra and extracellular potassium ion concentrations. ${E}_{z}^{\text{exp}}$ represents the effective *E*_{z} in the experiments, and actual *E*_{z} was computed by shifting ${E}_{\text{z}}^{\text{exp}}$ to correspond to *E*_{K} = −85 mV, i.e., *E*_{z} = ${E}_{\text{z}}^{\text{exp}}-({E}_{\text{K}}^{\text{exp}}+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 “two-stage 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.

#### 2.2.2. *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 ${\tau}_{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 ${\tau}_{s}^{*}$, we first defined ${\tau}_{s}^{*}$ more precisely [than in our previous work (Pathmanathan et al., 2019)] 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 ${\tau}_{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, ${\tau}_{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.

#### 2.2.3. *I*_{CaL} Parameter Uncertainty

The *L*-type inward calcium current *I*_{CaL} has six parameters: maximum conductance *g*_{CaL}, half-activation voltage *E*_{d}, activation slope *k*_{d}, half-inactivation voltage *E*_{f}, inactivation slope *k*_{f}, and inactivation time constant ${\tau}_{f}^{*}$.

In Pathmanathan et al. (2019), nominal values for *E*_{d}, *k*_{d}, *E*_{f}, and *k*_{f}, were taken directly from Iyer et al. (2012) (Table 2, normal, EPI). The same table in Iyer et al. (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 ${\tau}_{f}^{*}$ was chosen based on Figure 8 of Xiao et al. (2006). Here, we first re-defined ${\tau}_{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 ${\tau}_{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 ${\tau}_{f}^{*}$ using the above definition.

#### 2.2.4. *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 ${\tau}_{xr}^{*}$, half-inactivation 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}, ${\tau}_{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}, ${\tau}_{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), ${\widehat{\sigma}}_{\text{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 ${\widehat{\sigma}}_{\text{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 ${\widehat{\sigma}}_{\text{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.

#### 2.2.5. *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 ${\tau}_{xs}^{*}$. In Pathmanathan et al. (2019) nominal values of *E*_{xs}, *k*_{xs}, and ${\tau}_{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.

#### 2.2.6. 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 (normally-distributed), we defined ${\widehat{\sigma}}_{\text{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 ${\widehat{\sigma}}_{\text{LN}}$ as standard deviation divided by mean. Both ${\widehat{\sigma}}_{\text{N}}$ and ${\widehat{\sigma}}_{\text{LN}}$ are analogous to $\widehat{\sigma}$ in Pathmanathan et al. (2019).

### 2.3. 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 surface-area-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 (Pathmanathan and Gray, 2015; 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).

### 2.4. 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.

## 3. Results

### 3.1. Uncertainty in Parameters

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 $\widehat{\sigma}{({J}^{T}T)}^{-1}$, where $\widehat{\sigma}$ is the estimated residual variance and *J* is the Jacobian matrix with entries ${J}_{j}=\frac{\partial f}{\partial {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).

**Figure 1**. Individualized fits of *I*_{K1} parameters. **(A)** Experiment voltage clamp data (stars) and fitted model (line) for *n* = 12 cells. **(B)** Corresponding residuals (difference between experiment and fitted model). **(C)** Cell-specific parameters (solid diamonds) and cell-specific 90% confidence regions (ellipses) for the four fitted parameters, *g*_{K1}, ${E}_{z}^{exp}$, *k*_{z}, and ${E}_{K}^{exp}$, for each of the cells. See text for definitions of 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 ${\widehat{\sigma}}_{\text{N}}$, it is remarkable that all half-activation/inactivation values except *E*_{xs} have similar levels of uncertainty: each has ${\widehat{\sigma}}_{\text{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 $({\widehat{\sigma}}_{\text{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 ${\widehat{\sigma}}_{\text{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 (Pathmanathan et al., 2019).

### 3.2. Impact of Parameter Uncertainty on Action Potential

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 $\widehat{\sigma}=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.

**Figure 2**. Impact on action potential of uncertainty (due to population variability, derived from experimental data) in all **(A)** or individual currents **(B–F)**. Each figure displays 1,000 sample action potentials, together with APD histograms (*N* = 100, 000). Pie charts show first-order Sobol sensitivity indices.

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**. Results from non-behavioral analysis to determine which parameters are highly influential in different behaviors that occur in simulated action potentials under *I*_{CaL} and *I*_{to} parameter uncertainty.

### 3.3. Impact of Parameter Correlation

Table 2 shows that variability in *E*_{d} and *E*_{f} (steady state half-activation 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.

**Figure 3**. Effect of *introducing* correlation between half-activation voltage *E*_{d} and half-inactivation voltage *E*_{f} on *I*_{CaL} variability results. Probabilities of repolarization failure, early-after depolarizations (EADs) and loss of dome morphology all decreased with increasing 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.

**Figure 4**. Effect of neglecting inferred correlations in *I*_{K1} parameters. **(A)** 1,000 sample action potentials with histograms of APD (*N* = 100, 000) in inset. Red arrow indicates some non-physiological behavior that arises upon neglecting the correlation between *I*_{K1} parameters. **(B)** Corresponding sampled parameters *g*_{K1}, *E*_{z}, *k*_{z}. Blue dots: sampling from original distribution which included parameter correlation. Red dots *and* black diamonds: samples from distribution where correlation was removed. Black diamonds: parameter values which led to non-physiological action potentials.

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.

### 3.4. 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.

**Figure 5**. Impact of parameter variability on predictions of APD prolongation under 50% block of *I*_{Kr}. **(A)** APD prolongation under 50% *I*_{Kr} block, with no parameter variability. **(B)** corresponding action potentials (1,000 samples) when *I*_{K1} variability is included, with histograms of the control APs, blocked APs, and ΔAPD in inset.

**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.

### 3.5. 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.

**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.

**Figure 8**. Snapshots of activity at *t* = 1, 000 ms for 100 random parameters, under *I*_{to} parameter uncertainty.

**Figure 9**. Snapshots of activity at *t* = 1, 000 ms for 100 random parameters, under *I*_{CaL} parameter uncertainty.

**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.

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 potentials—it 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.

## 4. Discussion

### 4.1. 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 back-calculating 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 (Pathmanathan et al., 2019)]. 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 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 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.

### 4.2. 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 arbitrarily-chosen) independently-derived parameter uncertainty though a cardiac AP model, apart from our initial work on this subject for *I*_{Na} inactivation (Pathmanathan et al., 2015), 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), history-matching 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), 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 “bottom-up” 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.

### 4.3. 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 back-calculating 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 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.

### 4.4. 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.

### 4.5. 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 computationally-demanding 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.

## Author Contributions

PP devised the project, implemented the simulation and analysis code, ran the simulations and analysis, and wrote the paper. SG implemented and ran the analysis code, and reviewed the paper. JC provided the experimental data and reviewed the paper. AK and FF provided the simulation code and reviewed the paper. RG developed the model and provided the feedback on all aspects of the project. All authors contributed to the article and approved the submitted version.

## Funding

AK and FF are supported by funding from the National Institutes of Health (Award 1R01HL143450-01) and National Science Foundation (Award NSF1446675.Q1). This study was supported by the Free and Accepted Masons of New York, Florida, Massachusetts, Connecticut, Maryland, Wisconsin, Washington, and Rhode Island (to JC).

## Conflict of Interest

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

## Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2020.585400/full#supplementary-material

## Footnote

## References

Akar, F. G., Wu, R. C., Deschenes, I., Armoundas, A. A., Piacentino, III, V., Houser, S. R., et al. (2004). Phenotypic differences in transient outward K^{+} current of human and canine ventricular myocytes: insights into molecular composition of ventricular Ito. *Am. J. Physiol. Heart Circ. Physiol*. 286, H602–H609. doi: 10.1152/ajpheart.00673.2003

ASME V& 40 (2018). *Assessing Credibility of Computational Models Through Verification and Validation: Application to Medical Devices*. Washington, DC: Standard, American Society of Mechanical Engineers.

Ballouz, S., Mangala, M. M., Perry, M. D., Heitmann, S., Gillis, J. A., Hill, A. P., et al. (2019). Co-expression of calcium channels and delayed rectifier potassium channels protects the heart from proarrhythmic events. *bioRxiv* 659821. doi: 10.1101/659821

Berecki, G., Wilders, R., De Jonge, B., Van Ginneken, A. C., and Verkerk, A. O. (2010). Re-evaluation of the action potential upstroke velocity as a measure of the Na^{+} current in cardiac myocytes at physiological conditions. *PLoS ONE* 5:e0015772. doi: 10.1371/journal.pone.0015772

Berecki, G., Zegers, J. G., Verkerk, A. O., Bhuiyan, Z. A., de Jonge, B., Veldkamp, M. W., et al. (2005). Herg channel (dys) function revealed by dynamic action potential clamp technique. *Biophys. J*. 88, 566–578. doi: 10.1529/biophysj.104.047290

Britton, O. J., Bueno-Orovio, A., Van Ammel, K., Lu, H. R., Towart, R., Gallacher, D. J., et al. (2013). Experimentally calibrated population of models predicts and explains intersubject variability in cardiac cellular electrophysiology. *Proc. Natl. Acad. Sci. U.S.A*. 110, E2098–E2105. doi: 10.1073/pnas.1304382110

Chang, E. T., Strong, M., and Clayton, R. H. (2015). Bayesian sensitivity analysis of a cardiac cell model using a gaussian process emulator. *PLoS ONE* 10:e0130252. doi: 10.1371/journal.pone.0130252

Chang, K. C., Dutta, S., Mirams, G. R., Beattie, K. A., Sheng, J., Tran, P. N., et al. (2017). Uncertainty quantification reveals the importance of data variability and experimental design considerations for *in silico* proarrhythmia risk assessment. *Front. Physiol*. 8:917. doi: 10.3389/fphys.2017.00917

Cordeiro, J. M., Calloe, K., Moise, N. S., Kornreich, B., Giannandrea, D., Di Diego, J. M., et al. (2012). Physiological consequences of transient outward K^{+} current activation during heart failure in the canine left ventricle. *J. Mol. Cell. Cardiol*. 52, 1291–1298. doi: 10.1016/j.yjmcc.2012.03.001

Costabal, F. S., Matsuno, K., Yao, J., Perdikaris, P., and Kuhl, E. (2019). Machine learning in drug development: characterizing the effect of 30 drugs on the QT interval using Gaussian process regression, sensitivity analysis, and uncertainty quantification. *Comput. Methods Appl. Mech. Eng*. 348, 313–333. doi: 10.1016/j.cma.2019.01.033

Coveney, S., and Clayton, R. H. (2018). Fitting two human atrial cell models to experimental data using Bayesian history matching. *Prog. Biophys. Mol. Biol*. 139, 43–58. doi: 10.1016/j.pbiomolbio.2018.08.001

Faris, O., and Shuren, J. (2017). An FDA viewpoint on unique considerations for medical-device clinical trials. *N. Engl. J. Med*. 376, 1350–1357. doi: 10.1056/NEJMra1512592

Galappaththige, S. K., Pathmanathan, P., Bishop, M., and Gray, R. (2019). Effect of heart structure on ventricular fibrillation in the rabbit: a simulation study. *Front. Physiol*. 10:564. doi: 10.3389/fphys.2019.00564

Gong, J. Q., Susilo, M. E., Sher, A., Musante, C. J., and Sobie, E. A. (2020). Quantitative analysis of variability in an integrated model of human ventricular electrophysiology and β-adrenergic signaling. *J. Mol. Cell. Cardiol*. 143, 96–106. doi: 10.1016/j.yjmcc.2020.04.009

Gray, R. A., and Pathmanathan, P. (2018). Patient-specific cardiovascular computational modeling: diversity of personalization and challenges. *J. Cardiovasc. Transl. Res*. 11, 80–88. doi: 10.1007/s12265-018-9792-2

Herman, J., and Usher, W. (2017). Salib: an open-source python library for sensitivity analysis. *J. Open Source Softw*. 2:97. doi: 10.21105/joss.00097

Hindmarsh, A. C., Brown, P. N., Grant, K. E., Lee, S. L., Serban, R., Shumaker, D. E., et al. (2005). Sundials: Suite of nonlinear and differential/algebraic equation solvers. *ACM Trans. Math. Softw*. 31, 363–396. doi: 10.1145/1089014.1089020

Houston, C., Marchand, B., Engelbert, L., and Cantwell, C. (2020). Reducing complexity and unidentifiability when modelling human atrial cells. *Philos. Trans. R. Soc. A* 378:20190339. doi: 10.1098/rsta.2019.0339

Hu, Z., Du, D., and Du, Y. (2018). Generalized polynomial chaos-based uncertainty quantification and propagation in multi-scale modeling of cardiac electrophysiology. *Comput. Biol. Med*. 102, 57–74. doi: 10.1016/j.compbiomed.2018.09.006

Huberts, W., Heinen, S. G., Zonnebeld, N., van den Heuvel, D. A., de Vries, J.-P. P., Tordoir, J. H., et al. (2018). What is needed to make cardiovascular models suitable for clinical decision support? A viewpoint paper. *J. Comput. Sci*. 24, 68–84. doi: 10.1016/j.jocs.2017.07.006

Iyer, V., Heller, V., and Armoundas, A. A. (2012). Altered spatial calcium regulation enhances electrical heterogeneity in the failing canine left ventricle: implications for electrical instability. *J. Appl. Physiol*. 112, 944–955. doi: 10.1152/japplphysiol.00609.2011

Jost, N., Virág, L., Comtois, P., Ördög, B., Szuts, V., Seprényi, G., et al. (2013). Ionic mechanisms limiting cardiac repolarization reserve in humans compared to dogs. *J. Physiol*. 591, 4189–4206. doi: 10.1113/jphysiol.2013.261198

Kaboudian, A., Cherry, E. M., and Fenton, F. H. (2019a). Large-scale interactive numerical experiments of chaos, solitons and fractals in real time via GPU in a web browser. *Chaos Solit. Fract*. 121, 6–29. doi: 10.1016/j.chaos.2019.01.005

Kaboudian, A., Cherry, E. M., and Fenton, F. H. (2019b). Real-time interactive simulations of large-scale systems on personal computers and cell phones: Toward patient-specific heart modeling and other applications. *Sci. Adv*. 5:eaav6019. doi: 10.1126/sciadv.aav6019

Kaboudian, A., Velasco-Perez, H. A., Iravanian, S., Shiferaw, Y., Cherry, E. M., and Fenton, F. H. (2019c). *A Comprehensive Comparison of GPU Implementations of Cardiac Electrophysiology Models*. Cham: Springer International Publishing.

Lawson, B. A., Burrage, K., Burrage, P., Drovandi, C. C., and Bueno-Orovio, A. (2018). Slow recovery of excitability increases ventricular fibrillation risk as identified by emulation. *Front. Physiol*. 9:1114. doi: 10.3389/fphys.2018.01114

Lei, C. L., Clerx, M., Whittaker, D. G., Gavaghan, D. J., De Boer, T. P., and Mirams, G. R. (2020a). Accounting for variability in ion current recordings using a mathematical model of artefacts in voltage-clamp experiments. *Philos. Trans. R. Soc. A* 378:20190348. doi: 10.1098/rsta.2019.0348

Lei, C. L., Ghosh, S., Whittaker, D. G., Aboelkassem, Y., Beattie, K. A., Cantwell, C. D., et al. (2020b). Considering discrepancy when calibrating a mechanistic electrophysiology model. *Philos. Trans. R. Soc. A* 378:20190349. doi: 10.1098/rsta.2019.0349

Liu, D.-W., and Antzelevitch, C. (1995). Characteristics of the delayed rectifier current (IKr and IKs) in canine ventricular epicardial, midmyocardial, and endocardial myocytes: a weaker IKs contributes to the longer action potential of the m cell. *Circ. Res*. 76, 351–365. doi: 10.1161/01.RES.76.3.351

Milstein, M. L., Musa, H., Balbuena, D. P., Anumonwo, J. M., Auerbach, D. S., Furspan, P. B., et al. (2012). Dynamic reciprocity of sodium and potassium channel expression in a macromolecular complex controls cardiac excitability and arrhythmia. *Proc. Natl. Acad. Sci. U.S.A*. 109, E2134–E2143. doi: 10.1073/pnas.1109370109

Mirams, G., Arthurs, C., Bernabeu, M., Bordas, R., Cooper, J., Corrias, A., et al. (2013). Chaste: an open source C++ library for computational physiology and biology. *PLoS Comput. Biol*. 9:e1002970. doi: 10.1371/journal.pcbi.1002970

Mirams, G., Niederer, S., and Clayton, R. (2020). The fickle heart: uncertainty quantification in cardiac and cardiovascular modelling and simulation. *Philos. Trans. R. Soc. A* 378:20200119. doi: 10.1098/rsta.2020.0119

Morrison, T. M., Pathmanathan, P., Adwan, M., and Margerrison, E. (2018). Advancing regulatory science with computational modeling for medical devices at the FDA's Office of Science and Engineering Laboratories. *Front. Med*. 5:241. doi: 10.3389/fmed.2018.00241

Most, T. (2012). “Variance-based sensitivity analysis in the presence of correlated input variables,” in *Proceedings of 5th International Conference on Reliable Engineering Computing (REC)* (Brno).

National Research Council (2012). *Assessing the Reliability of Complex Models: Mathematical and Statistical Foundations of Verification, Validation, and Uncertainty Quantification*. Washington, DC: The National Academies Press.

Noble, D., Garny, A., and Noble, P. (2012). How the Hodgkin-Huxley equations inspired the Cardiac Physiome project. *J. Physiol*. 590, 2613–2628. doi: 10.1113/jphysiol.2011.224238

Oberkampf, W., Trucano, T., and Hirsch, C. (2004). Verification, validation, and predictive capability in computational engineering and physics. *Appl. Mech. Rev*. 57, 345–384. doi: 10.1115/1.1767847

Obreztchikova, M. N., Patberg, K. W., Plotnikov, A. N., Ozgen, N., Shlapakova, I. N., Rybin, A. V., et al. (2006). IKr contributes to the altered ventricular repolarization that determines long-term cardiac memory. *Cardiovasc. Res*. 71, 88–96. doi: 10.1016/j.cardiores.2006.02.028

Pathmanathan, P., Cordeiro, J. M., and Gray, R. A. (2019). Comprehensive uncertainty quantification and sensitivity analysis for cardiac action potential models. *Front. Physiol*. 10:721. doi: 10.3389/fphys.2019.00721

Pathmanathan, P., and Gray, R. A. (2015). Filament dynamics during simulated ventricular fibrillation in a high-resolution rabbit heart. *BioMed Res. Int*. 2015:720575. doi: 10.1155/2015/720575

Pathmanathan, P., Shotwell, M., Gavaghan, D., Cordeiro, J., and Gray, R. (2015). Uncertainty quantification of fast sodium current steady-state inactivation for multi-scale models of cardiac electrophysiology. *Prog. Biophys. Mol. Biol*. 117, 4–18. doi: 10.1016/j.pbiomolbio.2015.01.008

Rees, C. M., Yang, J.-H., Santolini, M., Lusis, A. J., Weiss, J. N., and Karma, A. (2018). The Ca^{2+} transient as a feedback sensor controlling cardiomyocyte ionic conductances in mouse populations. *Elife* 7:e36717. doi: 10.7554/eLife.36717

Sadrieh, A., Domanski, L., Pitt-Francis, J., Mann, S. A., Hodkinson, E. C., Ng, C.-A., et al. (2014). Multiscale cardiac modelling reveals the origins of notched T waves in long QT syndrome type 2. *Nat. Commun*. 5:5069. doi: 10.1038/ncomms6069

Strauss, D. G., Gintant, G., Li, Z., Wu, W., Blinova, K., Vicente, J., et al. (2018). Comprehensive *in vitro* Proarrhythmia Assay (CiPA) update from a Cardiac Safety Research Consortium/Health and Environmental Sciences Institute/FDA Meeting. *Ther. Innov. Regul. Sci*. 53, 519–525. doi: 10.1177/2168479018795117

Whittaker, D. G., Clerx, M., Lei, C. L., Christini, D. J., and Mirams, G. R. (2020). Calibration of ionic and cellular cardiac electrophysiology models. *Wiley Interdiscipl. Rev. Syst. Biol. Med*. 12:e1482. doi: 10.1002/wsbm.1482

Xiao, L., Zhang, L., Han, W., Wang, Z., and Nattel, S. (2006). Sex-based transmural differences in cardiac repolarization and ionic-current properties in canine left ventricles. *Am. J. Physiol. Heart Circ. Physiol*. 291, H570–H580. doi: 10.1152/ajpheart.01288.2005

## 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 ${\tau}_{X}^{*}$ have units ms.

For *I*_{K1}:

For *I*_{to}:

For *I*_{CaL}:

For *I*_{Kr}:

For *I*_{Ks}:

Keywords: uncertainty quantification, sensitivity analysis, variability, electrophysiology, correlation

Citation: Pathmanathan P, Galappaththige SK, Cordeiro JM, Kaboudian A, Fenton FH and Gray RA (2020) Data-Driven Uncertainty Quantification for Cardiac Electrophysiological Models: Impact of Physiological Variability on Action Potential and Spiral Wave Dynamics. *Front. Physiol.* 11:585400. doi: 10.3389/fphys.2020.585400

Received: 20 July 2020; Accepted: 20 October 2020;

Published: 19 November 2020.

Edited by:

Patrick M. Boyle, University of Washington, United StatesReviewed by:

Elisa Passini, University of Oxford, United KingdomTrine Krogh-Madsen, Weill Cornell Medicine, Cornell University, United States

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

*Correspondence: Pras Pathmanathan, pras.pathmanathan@fda.hhs.gov