There are two stages to UQ, which we will refer to as uncertainty characterization (UC) and uncertainty propagation (UP). Uncertainty characterization is the quantification of the uncertainty in model inputs. Here, “inputs” is a broad term for any quantity in the model whose value is based on real-world data. This includes model parameters, boundary and initial conditions, loading conditions, underlying geometry and material properties. The most common reasons for uncertainty in a model input are measurement uncertainty (e.g., the inability to measure a quantity exactly) and natural variability (e.g., variability in a physiological quantity across individuals in a population). The aim of UC is to determine probability distributions describing each of the inputs. This is generally a data-driven task that can be especially difficult for complex models with large numbers of parameters, where even estimating mean values can be challenging. The second stage, uncertainty propagation, involves propagating the input uncertainty through the model to derive the resultant uncertainty in important model outputs. This can be a computationally-challenging task, since it typically requires large numbers of simulations to be run. For a detailed introduction to UQ see Smith (2013).

There appear to be different interpretations of sensitivity analysis by different communities, and unfortunately no consistent terminology for distinguishing between them. In one interpretation, SA uses the distributions for model inputs (identified in the UC stage as discussed above), and involves calculating how the uncertainty in the model output can be apportioned to the uncertainties in inputs (Saltelli et al., 2008). In this interpretation, SA and UP are complementary activities: first UC is performed, then UP calculates the output uncertainty and SA identifies which inputs are responsible for that output uncertainty. This is global sensitivity analysis (GSA), since the entire range of permissible parameter values is considered. Local sensitivity analysis (LSA), on the other hand, focuses on how model outputs are affected when parameters are perturbed from a nominal set of values. GSA using empirically-derived input distributions provides a fundamentally different measure of sensitivity compared to LSA, since it uses information derived from experiment that is not used in LSA³. For a detailed introduction to SA see Saltelli et al. (2008).

To consider what rigorous fully comprehensive UQ might entail for a whole-heart EP model, first consider the components of such models as illustrated in Figure 1. Whole-heart models usually include a cellular action potential (AP) model (“cell model”), which are typically sets of ordinary differential equations (ODEs), and comprised of multiple sub-models of e.g., ion channels and sub-cellular processes. Notable human cell models proposed in recent years include the ten Tusscher model [TP06] (ten Tusscher and Panfilov, 2006) (system of 19 ODEs), the O'Hara & Rudy model [OR11] (O'Hara et al., 2011) (system of 41 ODEs) and the Grandi model (Grandi et al., 2010) (38 ODEs). These cell models are coupled to partial differential equations (PDEs) governing electrical wave propagation, usually either the monodomain or bidomain equations (Clayton et al., 2011). These are solved over a computational mesh representing the heart anatomy. To understand what a comprehensive UQ analysis might entail for a whole-heart model, consider all the inputs listed in Figure 1, using the broad definition of input as any quantity whose value is derived from real world data. All of these inputs are uncertain to some degree, and therefore a fully comprehensive UQ analysis, loosely equivalent to analyses done in other fields, would require characterization of uncertainty or variability in all of these inputs. This includes all the physiological parameters within the cell model (often numbering in the hundreds, for example the OR11 model has over 250 parameters), the geometry of the heart (for which there may be significant variability in shape and size across the human population), and the fiber-sheet orientation (often highly uncertain due to difficulty in measurement in vivo). Characterizing the uncertainty in all of these quantities is an immense experimental challenge. Moreover, a full uncertainty characterization also requires correlations between these quantities to be identified across the human population. For example, there is evidence for correlation of maximal conductances of I_Na and I_Kr (Milstein et al., 2012) and half-activation and half-inactivation voltages I_Na (Clerx, 2017). The inability or experimental difficulty in measuring such correlations is one of the biggest factors prohibiting comprehensive UQ with cardiac models. Correctly accounting for correlation may be key to successful UQ analysis, because, as implied in Figure 1, neglecting correlations may lead to overly wide ranging output distributions (through oversampling of non-physiological regions of parameter space).

Even if the uncertainty in all the model inputs could somehow be characterized, uncertainty propagation raises two further challenges. This first is the computational expense of running a massive number of whole heart simulations to propagate the input uncertainty through the model. Whole heart simulations are computationally expensive and generally require high performance computing resources; running millions (or more) of simulations to obtain estimates of model outputs is likely to be prohibitively expensive. Second, as the inputs are varied, it is very likely that some of the simulations will fail in some way or display a variety of behaviors (e.g., repolarization failure). Suppose this occurs during a formal UQ analysis for which the aim is to demonstrate that a model-derived conclusion (e.g., a clinical decision made using the model) is robust to the underlying uncertainties. If model failure occurs, the conclusion/decision is not robust to the input uncertainties. The second challenge is that with such complex models and such large parameter dimensions, it may be unclear how to resolve this issue when it occurs.

It is debatable as to whether anything similar to the above hypothetical example will ever be attained, but it is important to recognize that it is loosely equivalent to the level of UQ rigor used in other fields where modeling is used in safety critical decision making (Oberkampf and Roy, 2010). Currently, nothing that in any way resembles the above has been performed. Previous UQ/GSA research on cardiac cell models has typically involved a relatively small number of parameters in comparison with the total number of parameters in the model. Examples of traditional UQ analyses (i.e., empirically-derived input distributions, propagated through a model) include (Pathmanathan et al., 2015), in which UQ was performed in two parameters determining steady state I_Na inactivation, and a recent publication for the CiPA project (Chang et al., 2017), in which UQ was used to determine the impact of uncertainty in drug binding and drug ionic current block on drug risk stratification. While this provided important insight on the robustness of CiPA model predictions, it only considered uncertainty in drug related parameters, not uncertainty in the parameters of the cell model. Various alternative modeling frameworks have been devised for handling uncertainty, especially physiological variability. The population of models approach (Britton et al., 2013; Muszkiewicz et al., 2016; Passini et al., 2017) (essentially) derives distributions for certain parameters by calibration to AP recordings using acceptance/rejection criteria, and then performs UP by running simulations using the derived population. These papers have generally accounted for variability in conductances, with vast majority of the remaining parameters in the model held fixed. Others that have focused on conductance uncertainty include (Johnstone et al., 2016), who performed Bayesian calibration to AP data to calibrate conductances with uncertainty representing calibration error, and (Chang et al., 2015), who used a Gaussian process emulator for efficient UP and GSA after prescribing conductance uncertainty. A series of works (Sobie, 2009; Sarkar and Sobie, 2011; Sarkar et al., 2012) study the effect of variability by performing UP and GSA using multivariate regression, after introducing variability in a large number of parameters, for example up to 40 parameters including conductances, time constant scaling factors and steady-state voltage offsets in Sobie (2009), although even this was a minority of the total number of parameters in the cell model used. Hu et al. (2018) use the polynomial chaos method for efficient UP and GSA. They integrated variability in all parameters of two ionic currents (I_Kto and I_Kur) of a mouse cell model with 14 currents, but uncertainty in parameters in the other 12 currents was not considered. Numerous others have integrated uncertainty in cardiac models and performed UC, UP and/or GSA, for example (Geneser et al., 2006; Sadrieh et al., 2014; Krogh-Madsen et al., 2017; Costabal et al., 2019); see also the review of Ni et al. (2018). We believe that in all cases the number of parameters varied was significantly less than the total number of parameters in the model.

It is important to appreciate that the motivation behind much of the previous work on cardiac model variability was to understand cardiac (patho-)physiology and potentially develop new therapeutic targets or methods. Our ultimate motivation is different: we wish to understand the feasibility of comprehensive UQ for cardiac-model-based tools. There are two types of model input in such tools: fixed inputs (quantities that take the same value every time the tool is used), and variable inputs (quantities that take different values when the tool is used, i.e., the inputs into the tool itself). For example, for the CiPA computational tool, which predicts pro-arrhythmic risk of a drug, all parameters within the cell model are fixed inputs, whereas drug binding parameters and drug ionic current block are variable inputs. For a hypothetical clinical tool that uses patient imaging data to generate a patient-specific heart model and provide diagnostic/therapeutic guidance to the clinician, the fixed inputs include all parameters in the underlying cell model. The imaging data are variable inputs (Recall that we are using ‘input’ as a generic term for any quantity derived from real-world data). The relatively low number of variable inputs in cardiac-model-based tools suggest UQ limited to just the variable inputs will often be feasible (see for example Chang et al., 2017). We wish to understand if it is possible to demonstrate that clinical decisions made using cardiac-model-based tools are robust to all underlying uncertainties—both in the variable and the fixed inputs—the highest possible level of rigor.

For this goal, the current state of the art for cardiac cell model UQ/GSA is far from an ideal situation in which empirically-derived uncertainty in all cell model parameters is known and can be propagated through the model. Here we take a concrete step toward this goal by developing a novel canine cardiac cell model containing just 36 parameters to describe six major ionic currents using Hodgkin-Huxley formulations. This cell model was developed with UQ in mind, using carefully chosen simplifications to keep the number of state variables and parameters low; the approach was motivated by the philosophy that a simple(r) model with UQ may be more useful than a complex model without UQ (National Research Council, 2012). We then perform cell model UP and GSA accounting for uncertainty in all cell model parameters (excluding three environmental parameters). For this paper input distributions will be prescribed rather than empirically-derived, a common approach in cardiac modeling (Chang et al., 2015; Hu et al., 2018) (also see Table 1 of Ni et al., 2018 for a review), here taken as an initial step. In a future paper, we plan to perform full UQ with this cell model using empirically-derived input distributions—see section Discussion. As far as we are aware, the present paper represents the first time UP and GSA has been performed accounting for variation in all conductance and gating kinetics parameters in a cell model. We demonstrate the general feasibility of, and elucidate some of the challenges to, the computational side to comprehensive cell model UQ. We will demonstrate that our model is robust to (sufficiently small) interacting uncertainty in all parameters—the first time this has been shown for any cardiac cell model of moderate complexity, we believe. We further study model robustness by increasing parameter uncertainty, and explore when and how new model behaviors (such as repolarization failure or loss of spike-and-dome morphology) occur, and use Monte Carlo filtering to identify which parameters are responsible. We conclude with a discussion of the implications of these results for comprehensive UQ of whole heart models, in the context of the challenges laid out above.

A novel canine cell model was developed. Transmembrane voltage was modeled using an ODE with six ionic currents

where V is transmembrane voltage, Cm=1μF cm⁻² 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). I_stim is a square wave pulse used to initialize activity. Employing a Hodgkin-Huxley formulation for each current, each current is the product of a maximal conductance, voltage-dependent gating variable and a driving force. These were modeled as:

where the g_X are ion channel maximal conductances, m, h, z, r, s, d, f, x_r, y, and x_s are gating variables (activation gates: m, r, d, x_r, x_s; inactivation gates: h, z, s, f, y), and E_Na, E_K, E_Ca are the Nernst potentials for sodium, potassium and calcium, respectively. Each gating variable, Y say, has dynamics governed by the ODE

where Y_∞(V) is the steady state activation/inactivation function and τ_Y(V) is the voltage-dependent time constant. For steady state activation/inactivation functions, we chose sigmoid functions (ten Tusscher et al., 2004):

(− sign for activation gating variables, + sign for inactivation gating variables), where E_Y is the half-activation/inactivation voltage for that gating variable and k_Y > 0 controls the slope of the sigmoid. We did not model the full voltage dependence of all time constants that are typically included in other models. Instead, we chose

where τ_h0, δ_h, τm*, τs*, τf*, τxr*, τxs* are positive constants. This means that gating variables z, r, d, y are assumed to instantaneously reach their steady state values. Therefore, we can replace (2) to (7) with

Moreover, the gating variables which instantaneously reach steady state are not state variables for the system of ODEs, limiting the total number of state variables. Overall, the model has seven state variables, V, m, h, s, f, x_r, x_s, and is therefore a system of seven ODEs, and has 36 parameters—both significantly less than other cell models. The parameters are listed in Table 1.

‘Nominal’ parameter values were derived from either literature data or new fits to data from voltage clamp experiments on isolated canine cardiac myocytes. Table 1 provides an overview of the provenance of the parameter values. Except for three I_Na and four I_Kr parameters, all parameters were derived from canine data taken from adult canine epicardial myocytes under physiological conditions (i.e., temperature of 36-37 C, physiological extra- and intracellular ion concentrations).

Seventeen parameter values were taken directly from the literature. One (g_CaL) was derived by fitting to digitized data. Eleven parameters (E_m, k_m; I_K1 and I_to parameters) were fit to raw canine voltage clamp data. Four parameters, g_Na, g_Kr, τxr* and g_Ks, were not derived from data on their respective currents. Instead, g_Na was determined by simulating propagation in tissue (see below) to compute 1D conduction velocity, and chosen so that 1D conduction velocity was equal to 60 cm/s, to match longitudinal conduction velocity measurements in canine epicardial tissue (Kadish et al., 1986). Note that g_Na was calibrated in this way before the values of g_Kr, τxr* and g_Ks were finalized, but conduction velocity was found to be essentially independent of those values.

Finally, g_Kr, τxr*, and g_Ks were simultaneously determined by fitting simulated APD95 restitution data to data taken from Volders et al. (2003) (data digitized from Figures 5A,B, average of two sets of control experiments), while constraining the ratio of peak I_Kr to peak I_Ks under 1Hz pacing matched experimental measurements (0.46 μA/cm² and 0.11 μA/cm²; obtained by digitizing Figure 3 of Jost et al., 2013). It was determined that these three parameters are identifiable (at least locally) given this protocol. The parameters were not identifiable from restitution data if this constraint was not included, or if τxs* was also included in the fit.

The values of the initial conditions were derived from the parameter values. Initial V was set to be E_K, and initial values of the gating variables were set to be the steady state value of those gating variables at initial V.

Single cell and tissue simulations were run using Chaste, a general purpose package for computationally demanding physiological simulations. Chaste has been extensively tested and demonstrated to solve the governing equations of cardiac electrophysiology highly accurately (Niederer et al., 2011; Pathmanathan et al., 2012; Pathmanathan and Gray, 2014). For single cell simulations, the ODEs were solved in Chaste using an adaptive ODE solver. Simulation of electrical propagation through tissue was modeled using the monodomain formulation:

coupled to Equations (8)–(13), where χ = 1, 400cm⁻¹ is the surface-area-to-volume ratio, and σ = 1.4 mS/cm is the bulk conductivity. Tissue simulations were performed using the cardiac solver in Chaste, which uses the finite element method to solve the PDE (Pathmanathan et al., 2010).

UC involves determining probability distributions for each parameter, ideally based on experimental data and/or subject-matter expertise. The distributions should cover the parameters' uncertainty range, given the sources of uncertainty being accounted for (population variability, measurement uncertainty, etc). Determining empirally-derived probability distributions each of the parameters in Table 1 is an arguably feasible but difficult task. As discussed in section 1, here we prescribe input distributions. Specifically, we assumed that all parameters are independent and either normally-distributed (for parameters without physiological constraints) or log-normally distributed (parameters physiologically constrained to be positive). A ‘hyper-parameter’, σ^, was introduced, that controls the uncertainty across all parameters. By varying σ^, we can control the total amount of parameter uncertainty, and evaluate robustness of the model as uncertainty increases.

For parameters that are half-activation or half-inactivation voltages (E_m, E_h, E_z, E_r, E_d, E_f, E_xr, E_y, E_xs), we chose all parameters to be normally distributed with mean equal to the nominal value (Table 1), and standard deviation proportional to (with proportionality constant σ^) a reference range of 100 mV:

where p_nom is the nominal value and reference R = 100 mV. For conductances (g values), half-(in/)activation slopes (k values), time constants (τ values), and scaling factor δ_h, all of which are necessarily positive, we set:

Given the properties of the log normal distribution, p then has mean value p_nom and variance σ^2pnom2+O(σ^4), that is, standard deviation proportional to the nominal value, with proportionality constant, σ^.

As an illustration of the parameter ranges for a given σ^, consider parameters g_Na, E_m, g_CaL, and E_d, which are the maximal conductances and half-activation voltages of I_Na and I_CaL. These have nominal values of 12 mS/μF, -52mV, 0.115 mS/μF, -0.7 mV respectively (Table 1). When σ^=5%, the parameters become uncertain with the following 95% confidence intervals: g_Na: [10.9,13.2] mS/μF, E_m: [-61.8,-42.2] mV, g_CaL: [0.104,0.127] mS/μF, E_r: [-10.5,9.1] mV. Note that if we had chosen variances to be proportional to the absolute mean value for all parameters, as done elsewhere, E_r would be 95% certain to lie within [-0.77,-0.63] mV – much more tightly constrained than activation voltages not near zero. The use of a reference range R ensures inactivation/activation voltages have similar variability (for a given σ^).

Parameters E_K, E_Na, E_Ca were fixed as we consider them environmental parameters that can be controlled by simulating different extracellular ionic concentrations. Cm, χ and σ (monodomain conductivity) were also held fixed, for simplicity.

With these constructed distributions, we can study how the model behaves as we transition from σ^=0 (nominal model, no parameter uncertainty) to σ^=1% (small amount of uncertainty in all parameters (except equilibrium potentials), covering a small parameter-dependent range) to σ^=5% (larger uncertainty in all parameters, covering a larger parameter-dependent range). We re-iterate though that these distributions are prescribed rather than empircally-derived. Still, assessing the impact of uncertainty using prescribed distributions is common practice (Chang et al., 2015; Hu et al., 2018) and provides important insight about the behavior of computational model. Future work will focus on more accurate parameter uncertainty estimates derived from experimental data. In the discussion, section 4, we will describe how results in this paper provide information on allocating resources wisely when performing experiments for improved UC estimates.

Various AP characteristics or quantities of interest (QOIs) were analyzed. For the upstroke, these included: Threshold, the minimum stimulus current needed to induce depolarization (current applied for 0.5 ms); MaxUpstrokeVelocity, the maximum rate of change of voltage during depolarization; TimeOfMaxUpstrokeVelocity, the time of this maximum; and action potential amplitude (APA), defined as the difference between maximum voltage during upstroke and resting membrane potential. Plateau and repolarization QOIs analyzed included: NotchMin, the voltage attained at the local minimum during the action potential notch (see Figure 2) (if a notch is observed), NotchMax, the voltage attained at the local maximum at the end of the action potential notch (if a notch is observed), and action potential duration (APD), the time from activation for transmembrane potential to go below -70 mV. We measured APD using a fixed threshold (-70 mV) rather than APD90 (the time for 90% repolarization from maximum voltage), because APD90 is defined using APA and APD might then be determined to be sensitive to whichever parameters strongly influence APA.

For 1D simulations using a 1 cm long strand of tissue, QOIs included the above AP characteristics at the point 0.75 cm away from the stimulus site, as well as conduction velocity (CV).

Variance-based global sensitivity analysis (GSA) was performed by computing main and total Sobol sensitivity indices (Sobol', 1990). Sobol sensitivity indices provide fractional measures of the effect of the each parameter's uncertainty on the resultant variance of the model output. GSA approaches are far more computationally expensive than commonly-used local SA approaches, such as computing partial derivatives of the output with respect to each input, or one-at-a-time variation of each parameter. However, local SA approaches only provide information on sensitivity about the base point, resulting in incomplete or misleading information for highly nonlinear systems. GSA on the other hand explores the entire parameter space, using the distributions specified for the model inputs (section 2.4), and provides a more complete understanding of output sensitivity to inputs, including input interactions. Below is a brief introduction to this method; for a more complete introduction we refer the reader to Saltelli et al. (2008).

For a quantity of interest q = q(p), where p_i is the i-th parameter distributed as described in section 2.4, the expectation and variance of q are given by E(q)=∫q(p)dp and V(q)=∫q2(p)dp−(E(q))2. The i-th first-order Sobol sensitivity index is defined as

where V_i(q) = V(E(q|p_i)) (inner expectation over all parameters except p_i; outer variance over p_i only). The i-th first order Sobol index can be interpreted as the fraction of the output variance that can be attributed to the variance of parameter p_i, not accounting for interaction with other parameters. Similarly, second-order indices measuring the fraction of the output variance that can be attributed to interaction of parameters p_i and p_j (only) are defined as S_ij = V(E(q|p_i, p_j)) − S_i − S_j, and so on for higher-order indices. The total sensitivity index is given by S_Ti = S_i + Σ_jS_ij + Σ_{j, k}S_ijk + …. The total sensitivity index measures the total contribution of the parameter p_i to the variance of q, including through interactions with other parameters.

First-order and total Sobol indices were calculated using the python library SALib (Herman and Usher, 2017). The number of points sampled per dimension was chosen so that error estimates of the sensitivity indices, as provided by SALib, was small (see later figures).

For QOIs that are not computationally cheap to calculate (such as conduction velocity), computation of Sobol indices can be expensive since the total number of simulations, equal to the total number of points sampled, can be hundreds or thousands multiplied by the number of parameters. For such QOIs, the method of elementary effects, also referred to as Morris Screening (Morris, 1991), was used to identify parameters that were not influential (that is, the output has very low sensitivity to that parameter); those parameters were then excluded from the Sobol indices calculation. Morris Screening is a computationally cheap method (requiring only a few hundred simulations) that generally has a low false-positive rate when used to screen out non-influential parameters (Saltelli et al., 2008).

Simple Monte Carlo sampling was performed using the parameter distributions [Equations (14), (15)]; histograms and statistics were computed from the resultant model outputs. Monte Carlo approaches are straightforward to implement but can be very slow to converge. In all cases the number of samples, N, was chosen such that the output mean (standard deviation) were converged to three (two) significant figures, when estimated using either N or N/2 samples.

For some regions in parameter space, the model may transition from normal behavior to a different type of dynamics. For example, the action potential may fail to repolarize. If this occurs, probability distributions and sensitivity indices are not computable for derived quantities such as APD. There are several possibilities to consider when the model displays different classes of behavior:

Are the observed behaviors representative of what occurs in reality? One approach to addressing this question is to identify the mechanism underlying the behavior in the model, and confirming if that mechanism is physiologically reasonable. An alternative approach is direct statistical validation. For example, suppose a cell model fails to repolarize with small probability p, when parameter uncertainty representing population variability is included. It can then be asked if a random sample of isolated cardiomyocytes would also exhibit repolarization failure under identical stimuli with the same probability p.

If it is not believed that the observed behaviors are representative of reality, then one or both of the following should be considered:

2. Error in model form (structure of equations): one of the assumptions underlying the chosen governing equations may be being violated, at least in some regions of parameter space

Determining which path to take is difficult, and ideally requires empirically-determined input distributions, and therefore not in scope for the present paper. However, whichever of the above apply, it is very useful to identify which parameters are responsible for observed behavior. ‘Regionalized sensitivity analyses’ (Saltelli et al., 2008) are a group of methods that can be used to analyze model results exhibiting different types of behavior. They have been used in other fields but, as far as we are aware, has never been applied to cardiac models before. Here we use the Monte-Carlo filtering method (Saltelli et al., 2008), as follows: M points in parameter space, p₁, …, p_M, are randomly sampled, and the model solved using each. The sampled points are then split into two sets, those for which the behavior was observed, and those for which it did not. Each parameter p_i is then considered in turn. If p_i plays a role in determining whether the behavior occurs, the marginal cumulative distributions functions (CDFs) of (p_i|Behavior occurred) and (p_i|Behavior did not occur) will differ. The Kolmogorov-Smirnov test was be used to test whether the two CDFs were statistically different or not (α = 0.01). The Smirnov test statistic, D_stat (maximum difference between the two CDFs) was used to measure how influential a parameter was. Influential parameters were categorized as highly influential if the test statistic D_stat > 0.2.

We described in section 1.2 how the sheer number of quantities in cardiac models that are derived from physiological data is one of the greatest difficulties for cardiac UQ. Especially problematic are the number of parameters in modern cell models, for example, the OR11 model (O'Hara et al., 2011) has more than 250 numeric quantities in its governing equations that were derived from data in some way. Therefore, we developed a novel 36 parameter model, where each parameter has a clear physiological interpretation, for which it is arguably feasible to introduce empirically-derived uncertainty in all parameters. The model has various simplifications such as instantaneous gating and no intracellular ionic concentrations, which means it is not suitable for applications involving the mechanisms of excitation-contraction coupling, ischemia or reproducing the exact details of ionic current traces during voltage clamp experiments. However, it is suitable for investigating the impact of simultaneously variation in all parameters, and could serve as a starting point for development of incrementally more complex models, for example for developing a ‘minimally-complex’ model for a chosen application, for which comprehensive UQ remains possible. In general, there is a trade-off between less complex models which will likely exhibit less physiologically accurate behavior vs. more complex models which are less transparent and have greater numbers of uncertain inputs (Huberts et al., 2018). For cardiac models the crucial question is: can simpler models be as predictive as complex models for clinically-relevant applications? If so, simpler models have the advantage that fully testing robustness to parameter uncertainty is possible. This question has not been explored in any great depth but deserves greater attention—most recent research has been focused on whether the currently-available complex models are predictive for clinical applications.

Even with simple(r) models, characterizing the true uncertainty (whether population variability or measurement uncertainty) in all the parameters is a difficult experimental task (even ignoring the correlation question, discussed below). Sensitivity analysis and regionalized sensitivity analysis following UC with prescribed input distributions, as performed here, can provide guidance for allocating resources wisely. For example, the results in Table 2 and Figure 5 suggest it may be an inefficient use of resources to run experiments to estimate the population variability in g_Na, τm*, g_K1, g_to, k_r, E_s, k_s, k_d, k_f, k_xr, τxr*, g_Ks, k_xs, or τxr*, since none of these parameters appear in Table 2 or influence any of the QOIs in Figure 5. However, once improved uncertainty estimates have been obtained for the other parameters, the sensitivity analysis should be repeated, using crude but wide-ranging uncertainty in the above parameters, to confirm that they are still uninfluential. Also, and very importantly, when testing a model for a specific model application, the sensitivity analysis should be performed for the QOI to be used in decision-making, not just generic QOIs. Accordingly, the above results on parameter sensitivity may not extrapolate to spiral wave dynamics.

There are several methods that can be used to characterize the uncertainty due to population variability in cell model parameters, and the choice depends on the various factors, including the type of parameter, method used to derive the nominal parameter value (Table 1) and the availability of data from individual cells. For parameters for which voltage clamp data is available from individual cells, probability distributions can be derived by fitting parameters to individualized data, and then fitting probability distributions to the resultant parameters, similar to what has already been done for E_h and k_h in Pathmanathan et al. (2015). This is not possible for parameters for which only averaged data is available; for those a reasonable first approximation may be to use similar magnitude uncertainty as for corresponding parameters in other currents. Parameters that are calibrated using the full AP model, as opposed to derived from voltage clamp data or the literature (four parameters—see Table 1) may need to be re-calibrated with uncertainty, perhaps using Bayesian methods Johnstone et al., 2016). When suitable data is not available, a common approach is uncertainty elicitation using subject matter expertise (Morris et al., 2014). Obtaining/estimating uncertainty in time constant parameters will likely pose the greatest challenge, for two reasons. First, experimentally measuring time constants is more difficult than steady-state behavior. Second, time constants have known voltage-dependence, but in our model we approximated them as being independent of voltage (except for τ_h). This raises the question of what exactly is meant by population variability in this quantity. One solution is to define the parameter more precisely, for example as the average time constant over a pre-specified voltage range. This is a then well-defined quantity for which it is meaningful to ask what the population variability is.

Many of the parameters in cardiac cell models might be correlated across the species population. For example, Clerx describes how half-activation and half-inactivation voltages for I_Na appear to be correlated in human (Milstein et al., 2012; Clerx, 2017) shows potential correlation between g_Na and g_K1. The second major challenge mentioned in section 1 is the fact that identifying such correlations is experimentally very difficult, if not impossible. Again, model results can be used to guide experiments. First, correlations involving uninfluential parameters (see previous subsection) do not need to be considered. For influential parameters, distributions of model outputs assuming different input correlation structures could be computed, to determine parameters for which introducing correlation makes a difference. Alternatively, parameters which give rise to different types of model behavior could be assessed (as was presented in Figure 8), which may provide insight into which parameters to investigate experimentally. For example, we could hypothesize from the results in Figure 8 that the half-activation (E_d) and half-inactivation (E_f) voltages for I_CaL are positively correlated in reality (for canine), since larger values of |E_d − E_f| were associated with different types of non-physiological action potentials. This hypothesis is based on assumptions about the accuracy of the model form, and therefore needs to be verified experimentally.

We did not run any 2D or 3D simulations for this paper, nor was computational cost a focus of this paper. In fact, two of the methods used, Sobol sensitivity analysis and simple Monte Carlo sampling, are fairly computationally demanding even for single-cell analyses (although initial Morris Screening (section 2.6) reduces the computational cost for computing the sensitivity indices). Therefore, our results provide limited insight into the challenge of running large numbers of whole heart simulations to perform UQ for whole-heart simulation-based outputs. We expect that emulators of whole-heart models will need to be developed to overcome this challenge; see ongoing research (Chang et al., 2015; Ghosh et al., 2018; Lawson et al., 2018). Note that sensitivity indices were very similar for QOIs in 0D and 1D (Figure 5), as were output distributions (results not presented), suggesting that uninfluential parameters in single-cell simulations may be uninfluential in tissue simulations. It is not clear that this conclusion will hold for simulations involving arrhythmias however.

Even if comprehensive UC and UP could be performed with cardiac models, it is likely that a range of behaviors will be observed, and many common model outputs will not be computable for many of the sampled points in parameter space, for example APD when repolarization failure occurs. If this occurs, scientific conclusions or clinical decisions based on that model output are not robust to the input uncertainty. The final challenge described in section 1 was the fact that in situations such as this, it may not be clear how to proceed. In this paper we have observed how a range of AP dynamics occur as uncertainty in parameters is increased, and demonstrated how Monte Carlo filtering can be used to identify exactly which parameters are influential for each behavior (see section 3). This provides a pathway to investigate the root cause for the observed behaviors. In section 2.8 we discussed how one possibility is that the model is accurately reproducing different dynamics that occurs in reality. If this is not the case, the conclusion should be that the model failure occurred, in which case there are two options. Either (i) the model is an inaccurate representation of reality, at least in some regions of parameter space; and/or (ii) the uncertainty representations of some of the parameters are inappropriate in some way (e.g., do not represent true variability or do not account for correlations that occur in reality). For example, the results from section 3 suggest experimentally investigating whether I_CaL half-activation and half-inactivation voltages are correlated, because larger values of |E_f − E_d| were associated with both oscillatory APs and with loss of dome. Section 3 also showed how uncertainty in I_Na half-activation and half-inactivation voltages E_m and E_h, and in I_K1 half-inactivation voltages E_z, was responsible for non-physiological low voltage oscillatory behavior, suggesting focused investigation into those aspects of the model equations or refinement of uncertainty ranges for those parameters.

We end with a discussion of implications of our work for clinical applications of cardiac models. We focus on cardiac models as predictive tools, such as the in silico model used in the CiPA program for assessing drug cardiotoxicity (Li et al., 2018; Strauss et al., 2018). As discussed in section 1.3, for such tools, inputs can be categorized as either fixed (taking the same value every time the tool is used) or variable (the converse). For example, for the CiPA computational tool, all parameters within the cell model are fixed inputs, whereas drug binding parameters and drug ionic current block are variable inputs. It is undeniably important to perform UQ in variable inputs. Typically, uncertainty in those values with be due to measurement uncertainty, and if the tool output is not robust to this measurement uncertainty then the tool is not reliable. The importance of performing UQ in the fixed inputs is more debatable. For fixed inputs, the greatest source of uncertainty will often be due to population variability. Here, UQ in the fixed inputs can provide confidence that clinical decisions derived from the tool are robust to that underlying variability, especially as clinical trial results are extrapolated to a broader patient population. As an illustrative example, consider a hypothetical tool that has two inputs, patient-specific height (variable input) and average patient weight (fixed input). The clinical decision made using the tool certainly needs to be robust to any measurement error in patient height. Since any given patient may not be average weight, UQ to demonstrate that the clinical decision is robust to the uncertainty in weight, for the intended patient population, would provide additional confidence in the tool, and potentially reduce the need for a clinical trial cohort to fully cover the range of weights in the intended patient population. This is a simplistic example but the same ideas apply for variable and fixed inputs in patient-specific clinical tools based on personalized cardiac models (Gray and Pathmanathan, 2018).

The results in this paper represent a step toward cardiac model-based tools that are demonstrably robust to underlying uncertainties, for both fixed and variable inputs, but there remain many challenges to be overcome. For one, we used prescribed input uncertainty in this paper, the next step is to investigate AP variability with empirically-derived input distributions. If (when) overly large AP variability is observed, this will need to be addressed by one or more of: (i) refining the model governing equations; (ii) refining parameter uncertainty estimates; (iii) determining which parameters are correlated; and/or (iv) identifying parameters that should potentially be personalized in clinical applications because their uncertainty due to population variability leads to too much output uncertainty. This will likely result in several iterations of the following: model form refinement; parameter uncertainty characterization; robustness analysis as performed here; and analysis of different behaviors observed (section 2.8). This will need to be performed initially for simple pacing protocols (as performed here) and then arrhythmia-inducing protocols. Our model is canine, but the same approach can used to develop a human model, although parameterization will be more challenging. There are various open questions about how best to integrate uncertainty into tissue simulations (Ni et al., 2018). Finally, simplified models will be need to be shown to be predictive for clinical applications, and incrementally improved if there are not. This could be performed iteratively with UQ in each step as additional complexity is added. An incremental bottom-up approach may also allow formal methods for accounting for model form uncertainty (structural uncertainty) to be used (Mirams et al., 2016). Overall, although there is a long way to go, we believe this approach—development and iterative refinement of simplified models for which UQ in all parameters is possible—represents a complementary pathway for developing models for clinical cardiac applications, in comparison to the traditional approach where established high complexity models are used in clinical applications and UQ in all parameters is not possible. We expect that such research will reveal more information about the consequences of physiological variability in cardiac models, and on the general credibility of cardiac (and other physiological) models.

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.