Sensitivity and Uncertainty Analysis of Two Human Atrial Cardiac Cell Models Using Gaussian Process Emulators

Biophysically detailed cardiac cell models reconstruct the action potential and calcium dynamics of cardiac myocytes. They aim to capture the biophysics of current flow through ion channels, pumps, and exchangers in the cell membrane, and are highly detailed. However, the relationship between model parameters and model outputs is difficult to establish because the models are both complex and non-linear. The consequences of uncertainty and variability in model parameters are therefore difficult to determine without undertaking large numbers of model evaluations. The aim of the present study was to demonstrate how sensitivity and uncertainty analysis using Gaussian process emulators can be used for a systematic and quantitive analysis of biophysically detailed cardiac cell models. We selected the Courtemanche and Maleckar models of the human atrial action potential for analysis because these models describe a similar set of currents, with different formulations. In our approach Gaussian processes emulate the main features of the action potential and calcium transient. The emulators were trained with a set of design data comprising samples from parameter space and corresponding model outputs, initially obtained from 300 model evaluations. Variance based sensitivity indices were calculated using the emulators, and first order and total effect indices were calculated for each combination of parameter and output. The differences between the first order and total effect indices indicated that the effect of interactions between parameters was small. A second set of emulators were then trained using a new set of design data with a subset of the model parameters with a sensitivity index of more than 0.1 (10%). This second stage analysis enabled comparison of mechanisms in the two models. The second stage sensitivity indices enabled the relationship between the L-type Ca²⁺ current and the action potential plateau to be quantified in each model. Our quantitative analysis predicted that changes in maximum conductance of the ultra-rapid K⁺ channel I_Kur would have opposite effects on action potential duration in the two models, and this prediction was confirmed by additional simulations. This study has demonstrated that Gaussian process emulators are an effective tool for sensitivity and uncertainty analysis of biophysically detailed cardiac cell models.

1. Introduction

The cardiac action potential arises from the movement of ions through channels, pumps, and exchangers in the cell membrane. At any instant, the current carried by each ionic species depends on potential difference across the cell membrane, as well as ion concentrations and the dynamics of ion channel gating. The complex interplay of currents produces depolarization and then repolarization of the membrane, which then acts to trigger release of Ca²⁺, initiating mechanical contraction (Fink et al., 2011).

The first model of the action potential in a cardiac myocyte was developed over 50 years ago (Noble, 1962), and since then a series of more detailed models have been developed as experimental techniques and data have improved. The present generation of models provide detailed reconstructions of the cardiac action potential (Fink et al., 2011), and computational models of cardiac cells and tissue have become valuable research tools because they can encode biophysical mechanisms into a quantitative framework, and so can be used to test and construct hypotheses (Clayton et al., 2011).

Although these detailed models are capable of simulating the behavior of real cardiac myocytes, this veracity comes at the price of complexity. Models of the cardiac action potential typically comprise a system of coupled, stiff, and non-linear ordinary differential equations. There are many model parameters and boundary conditions, which we will refer to as model inputs from here onward. These model inputs can be derived from experimental data, using approaches based on those pioneered by Hodgkin and Huxley in squid giant axon (Hodgkin and Huxley, 1952). However, experimental data are subject to variability and error arising from both limitations of experimental methods as well as intrinsic variability in cardiac cells. Some of these inputs, such as binding affinities and reaction rate constants, can be considered to have fixed values because they have a physical basis. However, others, such as maximum ion channel conductance, depend on the ion channel density in the cell membrane as well as other factors that are variable. These quantities may therefore vary from one cell to another, and even from beat to beat in the same cell. These considerations underlie three specific problems. First, errors and variability in data are typically not taken into account when calibrating model inputs, and taking an average of experimental data can distort model behavior (Pathmanathan et al., 2015). Second, data from different sets of experiments can result in different models of the same cell type. For example, there are several models of the human atrial action potential, all based on human data, but which show different types of behavior (Cherry and Evans, 2008; Wilhelms et al., 2012). Finally, a further complication arises from the modular nature of cardiac cell models. The equations for a particular ion channel, pump, or exchanger are often re-used in different models and so the provenance of model inputs may be very difficult to establish (Niederer et al., 2009).

Addressing these problems requires tools and approaches that can quantify how model behaviors and outputs depend on model inputs. These include sensitivity analysis, which aims to quantify the change in model output resulting from changes in one or more model inputs, and uncertainty analysis, which considers how uncertain or variable model inputs defined by a distribution or range of values influence model outputs (Saltelli et al., 2019). However, the level of detail included in the present generation of cardiac cell models means that formal sensitivity and uncertainty analysis is difficult, and so a detailed examination of how model inputs influence model behavior often relies on large numbers of numerical simulations where a different set of inputs is used for each model run (Koivumäki et al., 2014). These datasets can be used for regression analysis, which enables the sensitivity of model outputs to changes in model inputs to be assessed (Sarkar et al., 2012). Another approach is to select a set of inputs, or population of models, that produce action potentials in the range of experimental observations (Britton et al., 2013; Sánchez et al., 2014; Muszkiewicz et al., 2015). A drawback of these approaches arises from the high dimensional input space for the models; a very large number of model evaluations is needed to investigate the input space thoroughly (Clarke et al., 2008), although recent work indicates that this challenge can be overcome by constructing new models that are designed for uncertainty analysis (Pathmanathan et al., 2019).

Methods for sensitivity and uncertainty analysis using probabilistic approaches have been developed and applied in areas including climate modeling (Lee et al., 2011, 2013), and are beginning to be used for cardiac and cardiovascular flow models (Eck et al., 2016; Mirams et al., 2016). One of these approaches is to use a Gaussian process (GP) as an emulator or surrogate of a model.

A GP is a flexible non-parametric regression tool widely used for machine learning, which effectively interpolates an output surface. The GP is initially trained using a set of design data composed of inputs and the desired output. Once trained, a GP can be evaluated very quickly to estimate an output for a new and unseen set of inputs, and so can emulate a cardiac cell model, or simulator (Chang et al., 2015). A GP can treat uncertainty and variability explicitly, providing a probabilistic (mean and variance) estimate of the output, and so is an ideal tool for sensitivity and uncertainty analysis. It can be trained on a relatively small number of simulator runs, and so offers advantages over other approaches, such as Monte Carlo methods that require large numbers of simulator runs (Sánchez et al., 2014; Melis et al., 2017).

The aim of this study was to extend our previous work on GP emulators of cardiac cell models (Chang et al., 2015; Johnstone et al., 2015; Coveney and Clayton, 2018) to demonstrate that this approach can be used on a larger scale for a systematic sensitivity and uncertainty analysis of two biophysically detailed models of the human atrial action potential. Our objectives were (i) to undertake a comprehensive sensitivity analysis of both models, taking into account first order effects and interactions; (ii) to compare sensitivity indices calculated using emulators with those calculated using a regression technique; (iii) to identify a subset of model inputs that have the greatest influence on model outputs; and (iv) to use these inputs to compare the behavior of the two models quantitatively.

2. Materials and Methods

2.1. Human Atrial Cell Models

Several models of the human atrial action potential have been developed and are reviewed in detail elsewhere (Cherry and Evans, 2008; Grandi et al., 2011; Wilhelms et al., 2012; Sánchez et al., 2014). We selected two models for the present study, both based on data from human atrial cells. The first model was the Courtemanche model (Courtemanche et al., 1998). The second model was an extension of the model developed by Nygren et al. (1998), with modifications to the I_Kur and I_to currents as well as the movement of Na⁺ (Maleckar et al., 2009) which we refer to as the Maleckar model.

We chose these models because both represent the action potential of human atrial cells, and so have clinical relevance. They have been used for tissue and whole-organ scale simulations of atrial fibrillation (McDowell et al., 2011; Krogh-Madsen et al., 2012; Colman et al., 2013; Tobón et al., 2013). Furthermore, both models have a comparable set of ion channels, pumps, and exchangers, but have different representations of intracellular Ca²⁺ handling and different action potential shapes (Cherry and Evans, 2008).

2.2. Model Inputs and Outputs

The Courtemanche and Maleckar cell models include components that represent membrane electrophysiology as well as intracellular Ca²⁺ storage, uptake, and release. We chose to concentrate on inputs that control the maximum current density carried by ion channels, pumps, and exchangers in the cell membrane as well as those that control the rate and magnitude of uptake and release of intracellular Ca²⁺. We also selected the cell capacitance C_m, and the extracellular concentrations ${\left[N{a}^{+}\right]}_o$, ${\left[K^{+}\right]}_o$, and ${\left[C{a}^{2+}\right]}_o$ as additional inputs. The inputs examined in this study are listed in Table 1, where the central values given are the default for each model. We used notation for model inputs and currents as given in the original manuscripts (Courtemanche et al., 1998; Maleckar et al., 2009), and as used in the CellML implementations (http://cellml.org). For example we refer to the HERG current as I_Kr with maximum conductance G_Kr, and the Kv1.5 current as I_Kur with maximum conductance G_Kur.

TABLE 1

Table 1. Range of inputs used for design data in each cell model.

The rationale for this choice was that each of the selected inputs can be considered uncertain (i.e., not a physical constant), yet has a biophysical interpretation. Maximum conductances of ion channels, pumps, and exchangers depend on protein expression, and so could be expected to vary within an individual cell at different times as well as from cell to cell. Cell size and capacitance vary from cell to cell. This natural variability can be considered to be aleatoric uncertainty, which is irreducible (Mirams et al., 2016). On the other hand, the kinetics of transmembrane currents are related to ion channel biophysics, and so could be considered epistemic uncertainty, which can in principle be reduced.

Our analysis proceeded in two stages. In Stage 1 the influence of all of the inputs listed in Table 1 was examined for a fixed pacing cycle length of 1,000 ms. In Stage 2, a subset of the inputs was selected on the basis of their Stage 1 sensitivity index (see below), and a new set of emulators was built using this subset as inputs. As for Stage 1, the simulations in Stage 2 were paced at a fixed pacing cycle length of 1,000 ms for 39 beats, followed by an S2 stimulus. The diastolic interval (DI) of the S2 beat in the S1–S2 pacing sequence was used as an additional input.

Cardiac cell models produce an output that is a time series of states. Of these, membrane voltage V_m and intracellular Ca²⁺ concentration ${\left[C{a}^{2+}\right]}_i$ describe the time course of action potentials and Ca²⁺. To investigate how cell model inputs influence action potential shape, we selected nine features of the action potential that quantify its shape, based on biomarkers used in related work (Britton et al., 2013; Sánchez et al., 2014) as well as the minimum and maximum ${\left[C{a}^{2+}\right]}_i$. These eleven outputs are shown in Figure 1 and are listed below.

• dV_m/dt_max—Maximum slope of the action potential upstroke.

• V_max—Peak voltage of the action potential.

• V₂₀, V₄₀, V₆₀, and V₈₀—Membrane voltage measured at 20, 40, 60, and 80% of APD₉₀.

• APD₅₀ and APD₉₀—Action potential duration at 50 and 90% of repolarization.

• RestV_m—Resting membrane potential, calculated as the average membrane voltage over a 10 ms period, 100 ms prior to the action potential upstroke.

• $C{a}_{min}^{2+}$ and $C{a}_{max}^{2+}$—Minimum and maximum intracellular ${\left[C{a}^{2+}\right]}_i$.

FIGURE 1

Figure 1. Action potential biomarkers. Nine action potential biomarkers were used as outputs to characterize each simulator run. (A) Courtemanche model, (B) Maleckar model. Biomarker labels are a: dV/dt_max, b: V_max, c–f: V₂₀, V₄₀, V₆₀, and V₈₀, g: APD₅₀, h: APD₉₀, i: RestV_m. $C{a}_{max}^{2+}$ and $C{a}_{min}^{2+}$ not shown.

2.3. Model Implementation

Both cell models were implemented in Matlab (Mathworks Inc.), using Matlab code automatically generated from the CellML repository (http://cellml.org). The models were solved using the ode15s time adaptive solver for stiff systems of ODEs with the relative and absolute tolerance both set to 10⁻⁶, and the maximum time step set to 0.5 ms.

To ensure that both cell models remained stable over the range of inputs used to train the GP emulators, we made some small modifications. Previous studies (e.g., Wilhelms et al., 2012) have identified an instability in the Courtemanche model that arises from a gradual drift in the intracellular concentrations [Na⁺]_i and ${\left[K^{+}\right]}_i$. We therefore fixed ${\left[N{a}^{+}\right]}_i$ and ${\left[K^{+}\right]}_i$ at their default initial values of 11.17 and 139.00 mM, respectively in the Courtemanche model implementation. In the Maleckar model, we fixed the I_{K, Ach} current to zero.

For each run, action potentials were initiated by an inward current of 2,000 pA delivered for 2 ms in the Courtemanche model, and 750/C_m pA/pF for 6 ms in the Maleckar model. In Stage 1, each run comprised 40 action potentials at a cycle length of 1,000 ms, with a check to ensure that APD₉₀ had reached steady state. In Stage 2, each run was composed of 39 S1 action potentials at a cycle length of 1000 ms and a final S2 stimulus delivered at an S1S2 interval determined by the APD₉₀ of the final S1 beat, plus an offset of 10 ms, plus a diastolic interval (DI) with a range of 50–450 ms sampled using a Latin hyper-cube design with the other selected inputs as described below.

2.4. Gaussian Process Emulators

Our overall approach is described in detail in a previous paper (Chang et al., 2015). We treat the cardiac cell models as simulators which produce a vector of model outputs y_s as a function of a vector of model inputs (parameters) x such that y_s = f_s(x). An emulator is then a statistical model of the simulator, sometimes known as a meta-model, a surrogate model, or a response surface model. The emulator approximates the model as y_e = f_e(x), where the emulator output approximates the simulator output y_e ≈ y_s for a given input x.

In the present study we specified the emulator as a GP, where the GP hyperparameters are optimized using a set of simulator runs called design data. When the GP has been trained, the posterior prediction y_e at an input x^* can be evaluated, which is a probability density with an expectation and a variance. The variance for the prediction y_e expresses uncertainty in the prediction of the simulator behavior at x^* (Oakley and O'Hagan, 2002, 2004).

2.5. Simulator Runs for Emulator Design Data

For Stage 1 we generated design data from 300 runs of each cell model implemented as described above. For each simulator run, a different set of inputs was generated using an optimized Latin hyper-cube design. We chose a range of variation that was a trade off between examining the effect of co-variation in model inputs while minimizing the number of model runs resulting in unphysiological behavior. We selected inputs within a range of ±50% of the central values (model defaults) given in Table 1 (i.e., from central value × 0.5 to central value × 1.5), except for G_K1, C_m, and extracellular ion concentrations. To reduce the incidence of unphysiological behavior (see below), these inputs had ranges of ±25, ±25, and ± 10%, respectively. Each input was scaled to vary between [0…1] over the designated range. A set of output biomarkers was obtained from the final action potential. A further set of 150 model runs were then used for emulator validation (see below). For Stage 2, a second set of design data was produced from 200 simulator runs of each model, with a reduced set of inputs sampled from a Latin hyper-cube as described in the results section, and other inputs set to their central value. A further set of 100 simulator runs were used for emulator validation.

Outputs from the model runs used for Stage 1 design data in the Courtemanche model and Maleckar model are shown in Figure 2. A wide variation in action potential shapes and durations were elicited by varying the model inputs.

FIGURE 2

Figure 2. Design data. Action potentials and calcium transients produced by Latin hyper-cube sampling as described in the main text, and running each model with a pacing cycle length of 1,000 ms. Traces shown in gray were used for emulator design data, and those shown in blue were excluded (see text for details). (A,B) Courtemanche model. (C,D) Maleckar model.

Model runs were removed from the design data if there was unphysiological behavior in the model outputs evidenced by pacemaking activity, a resting potential > −60 mV, APD₉₀ >600 ms, or if APD₉₀ of the 39th and 40th beats differed by more than 5% indicating alternans. Using these criteria, 5 out of 300 Courtemanche model runs in Stage 1 were removed, because of a long APD or APD alternans. In the Maleckar model 14 out of 300 model runs in Stage 1 were removed, all with pacemaking activity or failure to repolarize. In Figure 2 the removed runs are highlighted in blue. In Stage 2 no model runs were removed for the Courtemenache model, and 5 model runs were removed for the Maleckar model.

2.6. GP Emulator Training

Our approach to training and using GP emulators is described in full detail elsewhere (Chang et al., 2015; Johnstone et al., 2015). Mathematical details including the expression for the posterior prediction of the emulator are provided in Supporting Information as well as in Kennedy and O'Hagan (2001); Oakley and O'Hagan (2002, 2004), and the Python code used in this study for emulator training, validation, sensitivity and uncertainty analysis is available from https://github.com/samcoveney/maGPy. Briefly, each emulator was composed of a mean function and a zero mean GP,

$\begin{array}{ll}{f}_e(\text{x})=m(\text{x})+g(\text{x});& (1)\end{array}$

with a linear mean,

$\begin{array}{ll}m(\text{x})=h{(\text{x})}^T\beta ={\beta }_0+{\beta }_1{x}_1+\dots +{\beta }_P{x}_P,& (2)\end{array}$

and a zero-mean GP,

$\begin{array}{ll}g(\text{x})~{GP}(0,{\sigma }^{2}c(\text{x},{\text{x}}^{\prime })),& (3)\end{array}$

where the covariance has a Radial Basis Function form,

$\begin{array}{ll}c(\text{x},{\text{x}}^{\prime })=exp\left[-{\displaystyle \sum _{p=1}^P}\frac{{\left(x_p-x_p^{\prime }\right)}^2}{{\delta }_p^2}\right].& (4)\end{array}$

In these expressions x = (x₁, x₂, …, x_P) are the P inputs (parameters). The emulator hyperparameters β, δ, and σ define the emulator and need to be fitted to design data as the emulator is trained. β is a vector of length P + 1, δ a vector of length P, and σ² is a scalar. An optimized value for δ was obtained by maximum log-likelihood fitting to model inputs and outputs in the design data, assuming weak prior information on β and σ² (Kennedy and O'Hagan, 2001), and with a fixed nugget of 10⁻⁷ (Andrianakis and Challenor, 2012). Mathematical details for the training and fitting procedure are given in the Supporting Information. To avoid the training process becoming trapped in a local maximum, we repeated each fit ten times, each with a different set of randomly chosen initial values for the hyperparameters. The fit with the greatest log-likelihood was then selected. We produced a separate emulator for each of the outputs shown in Figure 1.

2.7. Sensitivity and Uncertainty Analysis

Sensitivity and uncertainty analysis can be seen as distinct but related topics; where variance based sensitivity analysis identifies the contribution of variance in each input to variance in each output, and uncertainty analysis concentrates on estimating the uncertainty in model outputs (Saltelli et al., 2019). A GP can be evaluated for uncertain inputs x where each input can be either a fixed value or a probability density expressed as a mean and a variance. A GP emulator can therefore be used for variance-based sensitivity analysis.

We calculated a first order sensitivity index for each combination of input and output (Oakley and O'Hagan, 2004). For each input w, the first order sensitivity index describes how much the output variance would be reduced if x_w is fixed, while all other inputs are uncertain and are described by a mean and variance. The first order index is expressed as the ratio of variance in the emulator output when x_w is fixed to variance in the output when all inputs are considered uncertain.

$\begin{array}{ll}{S}_w=\frac{Var[E(f_e(\text{x}|x_w))]}{Var[f_e(\text{x})]}& (5)\end{array}$

To capture the effect of interactions between the inputs, a total effect index can be calculated. This describes the reduction in output variance when x_w is uncertain and all other inputs are fixed, denoted as x_~w. It is also expressed as a proportion of the output variance when all inputs are considered uncertain.

$\begin{array}{ll}{S}_{Tw}=\frac{Var[f_e(\text{x})]-Var[f_e(\text{x}|x_{~w})]}{Var[f_e(\text{x})]}& (6)\end{array}$

The difference between S_Tw and S_w is then the contribution of all interactions between x_w and x_~w to the variance in the output. These quantities were calculated using expressions given in the Supporting Information and described in Oakley and O'Hagan (2004).

To calculate these indices, each uncertain input was assigned a mean of 0.5 in normalized units defined by the input ranges given in Table 1. Uncertain inputs that varied ±50% were then assigned a variance of 0.02, G_K1 and C_m were assigned a variance of 0.04, and the extracellular ionic concentrations were assigned a variance of 0.1.

2.8. Main Effects

The first order and total effect indices are both expressed as a ratio of variances, with a value in the range [0…1], and do not indicate whether an increase in output per change in input is positive or negative. To provide a sign to each index, we also calculated the main effect of each input on each output. The main effect for a single input x_w is the emulator output averaged over uncertain inputs, when x_w has a fixed value. Main effects for each combination of input and output were calculated using the procedure described in the Supporting Information, where inputs x_w were assigned fixed values in the range [0, 0.01…1]. The sign of the sensitivity indices was determined from the gradient of the main effect close to the central value of each input.

2.9. Comparison With Regression-Based Sensitivity Indices

Several recent studies have calculated sensitivity indices based on partial least squares (PLS) regression (Sobie, 2009; Sarkar and Sobie, 2010). In this approach, each model output is assumed to be a weighted sum of inputs. Thus, the model is described by the linear relationship

$\begin{array}{ll}\text{y}=\text{x}\text{B},& (7)\end{array}$

where y = (y₁, y₂, …, y_M) is a vector of M outputs, x = (x₁, x₂, …, x_P) a vector of P inputs, and B a P × M matrix of regression coefficients. An estimate of the matrix B, B_PLS, can be found by PLS regression on a set of design data obtained from N model runs, that generates an N × M matrix of inputs Y and an N × P matrix of outputs X. Each element of X, x_{i, j} is regularized by subtracting the mean $\overline{x_j}$ and dividing by the standard deviation of x_j, and each element of Y is regularized in the same way. B_PLS is found by minimizing the difference ||Ŷ−Ŷ||, where Ŷ = XB_PLS.

The matrix B_PLS can be interpreted as a matrix of sensitivity indices, provided the linear model holds. Each element of B_PLS, b_{i, j} describes how changing input x_i results in a corresponding change in output y_j. In both cases the change is relative to the mean value, and is a fraction of its standard deviation.

For comparison with variance based sensitivity indices, we calculated PLS sensitivity indices from the Stage 1 design data used to train the GP emulators. The input and output matrices X and Y were constructed from regularized design data inputs and outputs. The regression matrix B_PLS was then calculated using the Matlab function mvregress.

2.10. Emulator Validation

Each emulator was validated against an independent set of 150 simulator runs for Stage 1 and 100 simulator runs for Stage 2. For each output, we calculated the mean average predicted error (MAPE) and the median individual standard error (ISE) for each validation run. The MAPE was given by

$\begin{array}{ll}MAPE=\frac{100\%}{N}\frac{1}{\overline{y_s}}{\displaystyle \sum _{n=1}^N}|y_s^n-y_e^n|,& (8)\end{array}$

where N was the number of validation runs, $y_s^n$ simulator output for run n, and $y_e^n$ the posterior mean emulator output for run n. We used the mean of the simulator output $\overline{y_s}$ as a denominator instead of $\left|y_s^n\right|$ to avoid bias associated with small values of $\left|y_s^n\right|$.

$\begin{array}{ll}ISE=\frac{|y_s^n-y_e^n|}{\sqrt{Var(n,n)}},& (9)\end{array}$

where $y_s^n$ was the simulator output for run n, $y_e^n$ the posterior mean emulator output for run n, and Var(n, n) the posterior emulator variance for run n.

For most stage Stage 1 and Stage 2 emulators the MAPE was <10% and the median ISE was <1.0. Most of the differences between the output from the emulator and the output of the simulator for a given set of inputs were small, and so the fit of the emulators was considered acceptable. A table of MAPE and ISE is provided in the Supporting Information.

3. Results

3.1. Stage 1 Sensitivity Indices

The first order sensitivity indices for both cell models are shown in Figure 3. Each row of the figure corresponds to one of the model outputs, and each column represents a model input. We allocated a sign to each sensitivity index based on the slope of the main effect. Main effects are described in more detail below and illustrated in Figure 6. The sum of the absolute values of the sensitivity indices for each output is given to the right of each grid. Since first order sensitivity indices are a ratio, a sum close to one indicates that almost all of the output variance is accounted for by the variance on each input. Smaller values for this sum, such as those for APD₅₀ and APD₉₀, can indicate interactions among the inputs.

FIGURE 3

Figure 3. Stage 1 first order sensitivity indices. The first order sensitivity index for each input and output is given, with a sign based on the gradient of the mean effect. Sensitivity indices <0.05 not shown, to assist visualization. The numbers at the right hand side of the table indicate the sum of the absolute values of sensitivity indices along each row. (A) Courtemanche model. (B) Maleckar model.

The total effect indices are shown in Figure 4. For each combination of input and output, the difference between the total effect index and the first order index reflects interactions with the other inputs. The sum of these differences across all of the inputs is shown at the right hand side of Figure 4. In most cases the first order and total effect indices were similar, and the sum of differences was small (≤0.12) indicating that the effect of interactions was also small. However, the sum of differences was larger for V₈₀ in the Courtemanche model, as well as for APD₅₀ and APD₉₀ in both models. We conclude that in both models there are some interactions between inputs and that these interactions have an effect on repolarization. However, comparison of the sensitivity indices in Figures 3, 4 shows that these interactions appear to be distributed among all of the inputs.

FIGURE 4

Figure 4. Stage 1 total effect indices. The total effect sensitivity index for each input and output is given, indices <0.01 not shown. Numbers on right hand side are the sum across each row of the differences between the total effect index and the absolute value of the first order index. (A) Courtemanche model. (B) Maleckar model.

For comparison, Figure 5 shows sensitivity indices obtained by PLS regression on the emulator design data, and a comparison with variance based first order indices. The comparison plots show broad agreement, with the first order indices $S_i\approx B_i^2$. This relationship arises from the different definitions of S_i and B_i based on variance and standard deviation, respectively (Saltelli et al., 2012).

FIGURE 5

Figure 5. Multivariate regression sensitivity indices. The sensitivity index based on multivariate regression for each input and output is shown for the Courtemanche model (A) and Maleckar model (B). In the lower panel, these sensitivity indices are compared with those obtained from GP emulators for the Courtemanche model (C), and the Maleckar model (D).

Overall these sensitivity indices show the contribution of uncertainty in each input to uncertainty in each output. Thus, the main contributors to uncertainty in dV/dt_max are the Na⁺ channel maximum conductances G_Na and P_Na, and the membrane capacitance C_m. The sign of the sensitivity indices show that these act in opposite directions, as would be expected from the role played by the Na⁺ current in depolarization: increasing Na⁺ current acts to increase dV/dt_max, whereas increasing C_m acts to decrease dV/dt_max. The bigger influence of C_m in the Maleckar model arises because the stimulus current density in this model scales with 1/C_m; a larger C_m results in a smaller stimulus current, which in turn produces a smaller dV/dt_max and a smaller V_max.

In both models, these sensitivity indices can be interpreted to show that increased outward currents (for example arising from increased G_Kur) act to decrease both the voltage of the action potential plateau and action potential duration, whereas increased inward currents resulting from increased G_GaL and G_bCa have an opposite effect. This behavior is broadly what would be expected, and confirms that the sensitivity indices quantify model behavior, and can reasonably be extended to understand more complex relationships between inputs and outputs, as well as comparing behavior of the two models. In the Maleckar model I_NaKMax has the opposite effect to G_bCa. In the Courtemanche model, $C{a}_{max}^{2+}$ and $C{a}_{min}^{2+}$ are influenced by G_CaL, G_bCa, I_NaCaMax and inputs that control Ca²⁺ handling I_rel and K_up, whereas in the Maleckar model G_CaL and K_NaCa have a negligible effect, and this reflects the different formulation of Ca²⁺ handling in the two models (Cherry and Evans, 2008).

The sensitivity analysis shows that several inputs influence APD₅₀, APD₉₀, RestV_m, and $C{a}_{max}^{2+}$. Figure 6 shows the main effect of each input on these outputs, for each cell model. The main effect shows the expected value of the output as each input is fixed and varied in turn across the normalized range 0…1 corresponding to the input ranges given in Table 1, while all other inputs are considered uncertain. The residual variance arising from the uncertain inputs accounts for the fact that the main effects do not converge exactly for an input value of 0.5.

FIGURE 6

Figure 6. Stage 1 selected main effects plots. Main effects for APD₅₀, APD₉₀, RestV_m, and $C{a}_{max}^{2+}$. (A) Courtemanche model. (B) Maleckar model.

Some of the main effects are comparable between the two models, for example increasing G_bCa acts to increase APD₉₀, RestV_m, and $C{a}_{max}^{2+}$. Several of the effects are non-linear, for example the main effect of G_CaL on APD₅₀ and APD₉₀. However, the overall picture is complex, and it is hard to compare the different models. In order to simplify the analysis, we selected a subset of inputs for Stage 2 of the analysis based on their sensitivity indices as described below.

3.2. Stage 2 Sensitivity Analysis

For Stage 2, we concentrated on inputs that strongly influenced action potential shape and duration, with first order sensitivity index of more than 0.1. To simplify the analysis further, we excluded G_Na and P_Na as these inputs mainly influence action potential upstroke and amplitude. We also excluded extracellular concentrations, since these are tightly controlled in normal physiological conditions, and we excluded the inputs directly involved in the storage, uptake and release of intracellular Ca²⁺ because we sought to concentrate on action potential shape and duration. The Stage 2 inputs selected for the Courtemanche model were G_K1, G_to, G_KurMult, G_CaL, G_bCa, I_NaKMax, I_NaCaMax, and I_PCaMax, and for the Maleckar model G_K1, G_t, G_Kur, G_CaL, G_bCa, I_NaKMax, and C_m. All other inputs were assigned their central value from Table 1. In addition, the DI of the final action potential was introduced as an additional input to explore the dynamic behavior of the model, and this final action potential was the one that was analyzed.

The sensitivity indices for Stage 2 are shown in Figures 7, 8, and the main effects for APD₅₀, APD₉₀, RestV_m, and $C{a}_{max}^{2+}$ in Figure 9.

FIGURE 7

Figure 7. Stage 2 first order sensitivity indices. The first order sensitivity index for each input and output is given, with a sign based on the gradient of the mean effect. Sensitivity indices <0.05 not shown. The numbers at the right hand side of the table indicate the sum of sensitivity indices along each row. (A) Courtemanche model. (B) Maleckar model.

FIGURE 8

Figure 8. Stage 2 total effect sensitivity indices. The total effect index for each input and output is given, indices <0.01 not shown. (A) Courtemanche model. (B) Maleckar model.

FIGURE 9

Figure 9. Stage 2 selected main effects plots. Stage 2 main effects for APD₅₀, APD₉₀, RestV_m, and $C{a}_{max}^{2+}$. (A) Courtemanche model. (B) Maleckar model.

The first order and total sensitivity indices for Stage 2 were similar to those obtained in the Stage 1 analysis, which is consistent with the idea that the inputs fixed for Stage 2 had only a small effect on the outputs. DI was included as an additional input for Stage 2. In the Courtemanche model both first order and total effect indices for DI were larger than for the Maleckar model, indicating that DI has a greater influence on both the shape and duration of the action potential in the Courtemanche model compared to the Maleckar model.

The main effects plots show opposite effects of inward and outward currents on action potential duration; however APD₅₀ as a proportion of APD₉₀ in the Maleckar model was considerably shorter than in the Courtemanche model as a result of the different action potential shape, and so both the sensitivity indices and main effects for APD₅₀ may not be easily comparable. This observation is reflected in the larger main effect of DI on APD₅₀ in the Courtemanche model. The main effects plots also show several non-linear relationships, for example the main effect of G_Kur was non-linear in each of the outputs shown in Figure 9.

An additional observation from Figure 9 was that the main effect of G_Kur on APD₉₀ in the Maleckar model was larger, and in an opposite direction, to the main effect in the Courtemanche model. This is a prediction from the emulator analysis, and we undertook additional simulations of the cell models to verify that the emulator prediction was correct. In these additional simulations G_Kur was decreased to 50% and increased to 150% of its default value in both models, and each model was run for 40 beats paced at a cycle length of 1,000 ms as in the Stage 1 analysis. The outcome of these simulations is shown in Figure 10, which shows simulated action potentials for each model, together with the principal inward and outward currents that act during the action potential plateau.

FIGURE 10

Figure 10. Effect of magnitude of I_Kur on inward and outward currents. Simulated action potential, outward current I_Kur, inward current I_{Ca, L}, and outward current I_Kr. In each case the final beat of 40 is shown, with a pacing cycle length of 1,000 ms. (A) Courtemanche model. (B) Maleckar model.

The top panels of Figure 10 show that the prediction from the emulator main effects is confirmed in these simulations: decreasing G_Kur results in a longer APD in the Maleckar model and a slightly shorter APD in the Courtemanche model, and vice-versa. Based on the simulation results shown in Figure 10, we would speculate that the mechanism by which G_Kur influences APD is the secondary effect of changes in action potential plateau voltage on the balance of inward and outward currents during the plateau and repolarization of the action potential, and this is consistent with earlier mechanistic studies (Greenstein et al., 2000).

Changing G_Kur influenced the voltage of the action potential plateau in both models, but had a different effect on the timing of repolarization. The time course of I_to in both models was similar, and was not strongly influenced by G_Kur and so is not shown. A decrease in G_Kur reduced the outward current during the initial part of the action potential plateau. This resulted in an increased voltage during the plateau, and a larger inward I_{Ca, L}, which is voltage-dependent and acts to prolong the action potential plateau. In turn, the increased plateau voltage resulted in greater activation of the outward current I_Kr, which is larger in the Courtemanche model compared to the Maleckar model. Thus, in the Courtemanche model, a decrease in G_Kur resulted in increased I_Kr, with little change in action potential duration. In the Maleckar model I_Kr is much smaller, and so the increased plateau voltage did not result in increased outward current during repolarization, and so action potential duration was prolonged.

3.3. APD Restitution

The Stage 2 analysis included diastolic interval as an input, which enabled us use the emulators to examine how different inputs affect APD restitution. In Figure 11 we have plotted a surface showing the expected value of APD₉₀, colored by the 95% credible interval (see Supporting Information). In each of these plots the emulators were evaluated with all inputs assigned fixed values with no uncertainty, and so the 95% credible intervals reflect only uncertainty in the emulator predictions, with no uncertainty arising from uncertainty in the inputs. We assigned a value of 0.5 in normalized units to all of the inputs, except for DI and another inputs that were varied in each plot; these were assigned fixed values between 0 … 1.

FIGURE 11

Figure 11. APD₉₀ as a function of DI and other inputs. Each plot shows a surface of predicted mean APD₉₀, colored by the 95% credible interval. (A) Courtemanche model. (B) Maleckar model.

Overall, APD₉₀ restitution was flatter in the Maleckar model compared to the Courtemanche model, which is consistent with other studies (Cherry and Evans, 2008; Wilhelms et al., 2012). The overall effect of G_K1, G_CaL, and G_bCa on APD₉₀ was similar in each model, with I_NaKMax in the Courtemanche model having a similar effect to G_Kur in the Maleckar model.

However the shape of the APD₉₀ restitution was modulated to some extent by the other inputs shown. In the Courtemanche model, decreasing G_bCa resulted in steepening of APD₉₀ restitution, with a marked decrease in APD₉₀ for short DI as G_bCa was changed from 150 to 50% of its central value. The Maleckar model also showed a marked decrease in APD₉₀ as G_bCa was reduced, but this change was seen across the full range of DI.

3.4. Uncertainty Analysis

Uncertainty analysis aims to quantify how uncertain or variable model inputs affect uncertainty in a model output. For a cardiac cell model, uncertainty analysis enables the effects of hypothesized natural variability in model inputs to be quantified. Using the GP emulators for each cell model, we assessed how the variance in model outputs changed as variance in all of the model inputs was increased, and the results for APD₉₀ are shown in Figure 12.

FIGURE 12

Figure 12. Uncertainty in APD₉₀ and V_max. Coefficient of variation (standard deviation divided by mean) in APD₉₀ for Stage 1 (A) and Stage 2 plotted against imposed uncertainty on all inputs (B) for each model, with different mean values of 0.4, 0.5, and 0.6 in normalized units.

Figure 12A shows how the uncertainty in APD₉₀ (expressed as coefficient of variation) changed as all uncertainty in all the Stage 1 inputs (expressed as standard deviation) increased, and Figure 12B shows the effect of increasing uncertainty in all the Stage 2 inputs. Uncertainty in APD₉₀ increased monotonically with uncertainty in the inputs, for both Stage 1 and Stage 2 inputs, and in both cell models. Changing the mean value of the inputs had a small effect on output uncertainty, larger mean values reduced output uncertainty.

Reducing the number of uncertain inputs for Stage 2 acted to reduce the output uncertainty for the Maleckar model, but increased output uncertainty slightly in the Courtemanche model. This rather unexpected finding may be due to increased interactions among the inputs, or a more problematic emulator fit. However, in the Maleckar model, the total effects for APD₉₀ in Figure 4B are higher than the total effects for the Courtemanche model, indicating a greater degree of interaction in the Maleckar than Courtemanche models. The smaller number of inputs for Stage 2 could then result in less uncertainty arising from interactions among uncertain inputs for the Maleckar model, and a consequent reduction in uncertainty. However, fixing some inputs for the Stage 2 analysis may have had increased interactions in the Courtemanche model (Figure 8A), leading to an increase in uncertainty in predicted APD₉₀.

4. Discussion

The main findings of this study are:

• We have demonstrated the use of GP emulators for a systematic and quantitative uncertainty and sensitivity analysis of biophysically detailed cardiac cell models.

• The first order and total effect sensitivity indices obtained using GP emulators showed both relationships that would be expected, such as the dependence of action potential upstroke on G_Na, as well as those that were more unexpected, such as the almost negligible effect of G_Ks on APD. Most of the variance in model outputs was accounted for by first order sensitivities, so the interactions between the inputs considered in the present study are small.

• There was broad agreement between first order sensitivity indices obtained using GP emulators and those obtained using PLS regression.

• A subset of model inputs was identified, and behaviors predicted by the emulators were confirmed by model simulations.

We discuss these findings in more detail below, and then highlight limitations, challenges, and future directions.

4.1. Gaussian Process Emulation of Biophysical Models

As models of cardiac cell and tissue electrophysiology become more widely used, it is becoming increasingly important to understand how different components of the models influence model behavior, and especially how these different components interact. Biophysically detailed cardiac cell models are complex, with many interacting parts. Some of these model components may be inherited from earlier models and experiments (Niederer et al., 2009), and the process by which model parameters are fitted is also fragile when there are uncertainties associated with experimental data (Pathmanathan et al., 2015). The development and evaluation of tools for sensitivity and uncertainty analysis of cardiac models is therefore an important and growing area (Mirams et al., 2016), but much remains to be done.

The ability to evaluate emulators cheaply can be valuable for model calibration, where thorough exploration of a high dimensional input space is required (Coveney and Clayton, 2018), or if the GP emulator is used to construct outputs distributions with Monte Carlo methods. A key benefit of a GP emulator approach is the explicit handling of uncertainty. Under the assumption that inputs and outputs have Gaussian distributions, the variance on the emulator output can be calculated directly, given variances on the inputs. This enables the direct calculation of sensitivity indices, as well as enabling a systematic investigation of the way that output uncertainties depend on uncertainties in the inputs.

In this study, we have demonstrated how GP emulators can be used to systematically and quantitatively analyse biophysically detailed cardiac cell models. A GP emulator must be trained on design data. The number of simulator runs required to for design data remains an open question, and will depend on the complexity of the model output surfaces. A typical rule of thumb is to use ten times the number of inputs, and based on our previous experience (Chang et al., 2015) we opted for 300 runs for Stage 1 and 200 runs for Stage 2. We also trained the Stage 1 emulators on design data sets composed of 200 and 400 simulator runs, and obtained similar sensitivity indices to those presented here. However, further work is required to develop methods to determine the number of simulator runs needed as well as suitable metrics to determine emulator quality. Emulators should be trained on design data that fill the input space evenly, and we chose to use Latin hypercube sampling in this study (McKay et al., 1979). Other methods, such as orthogonal sampling (Bingham et al., 2009) may provide a better sampling strategy.

4.2. Sensitivity and Uncertainty Analysis

Recent studies have pioneered the use sensitivity indices obtained by partial least squares (PLS) regression of simulator outputs on simulator inputs, which allows a calculation of sensitivity indices (Sobie, 2009; Sarkar et al., 2012; Koivumäki et al., 2014). This approach is straightforward to implement, and we have found that it gives sensitivity indices that agree well with the square root of the first order index obtained using the GP approach (Figure 5), and the reason for this appears to be that the PLS and GP indices are based on variance and standard deviation, respectively. The overall agreement indicates that both approaches can yield similar first order sensitivity indices, although the GP easily enables calculation of interaction effects as well as first order indices. Other approaches for uncertainty and sensitivity analysis based on generalized polynomial chaos expansion have also been developed and used for analysis of cardiovascular system models (Eck et al., 2016). These approaches also enable calculation of sensitivity indices, but the relative merits of these different approaches are only beginning to be explored (Johnston et al., 2018). Both GP emulators and polynomial chaos expansions enable the explicit treatment of uncertainties, and so offer advantages for more comprehensive model analysis.

4.3. Human Atrial Cell Models

Models of the human atrial action potential are a subject of research interest and clinical relevance because heterogeneity in action potential shape and duration in different parts of the atria is associated with vulnerability to atrial fibrillation (Varela et al., 2016), and persistence of atrial fibrillation is associated with remodeling of the atrial action potential (Krogh-Madsen and Christini, 2012; Colman et al., 2013). There are several different models of human atrial myocytes, each with different properties (Cherry and Evans, 2008; Cherry et al., 2008; Wilhelms et al., 2012). More recent models are extensively based on data from human myocytes (Grandi et al., 2011). Most analyses of these models have focussed on the mechanisms that change action potential duration because a reduced APD increases vulnerability to atrial fibrillation, and APD can be modulated pharmacologically (Nygren et al., 1998; Zhang et al., 2005; Cherry and Evans, 2008; Koivumäki et al., 2014; Sánchez et al., 2014). These previous studies have highlighted the importance of I_CaL, as well as I_K1 and I_Kur, in regulating APD. The present study adds to our understanding of the Courtemanche and Maleckar models by providing a more comprehensive view of how the model parameters affect the shape and duration of the action potential, as well as the maximum and minimum of the Ca²⁺ transient.

Increased inward current tends to increase amplitude of the upstroke and plateau as well as increasing APD, and increased outward current tends to have the opposite effect. The present study has highlighted how action potential duration and shape depends on the net flow of charge across the cell membrane. It is well-accepted that net current flow is finely balanced so small changes in inward and/or outward currents can influence action potential shape and duration. Our approach enables these effects to be investigated in a quantitative way. The two models examined in this study represent repolarization using a different balance of currents, and this difference is seen in the main affects plot for APD90 in Figure 9 and the consequences can be seen in Figure 10. The relative magnitudes of I_CaL and I_Kr in the two models are different, with much smaller I_CaL and I_Kr in the Maleckar model. This may explain why G_bCa exerts a stronger influence over $C{a}_{min}^{2+}$ and $C{a}_{max}^{2+}$ in the Maleckar model compared to the Courtemanche model.

Overall, the difference between the first order sensitivity indices (Figure 3) and the total effect indices (Figure 4) was small. The sum of these differences for each output (right hand column in Figures 4, 8) indicates more interactions in the Courtemanche model than in the Maleckar model, and that these interactions tend to affect the plateau of the action potential. These observations mean that overall the interactions between the inputs in these models do not have a strong effect on the outputs, and so we can conclude that the inputs examined in this study tend to act independently. This is a potentially important feature of the models, which could be exploited for model calibration as well as for examining mechanisms of remodeling and pharmacological action. However, it remains to be seen whether this independence is a feature of real cardiac myocytes.

4.4. Limitations, Challenges, and Future Directions

The use of emulators to probe detailed biophysical models is at an early stage, and so there are several limitations and challenges associated with the present study.

4.4.1. Choice of Inputs

We concentrated on the effect of maximum conductances in the present study, as this reduced the complexity of the analysis. The rationale for this approach was an assumption that kinetic parameters are determined by biophysics, and so less prone to variation than the expression of ion channels, pumps, and exchangers. However, a detailed sensitivity analysis of I_Kr dynamics in the Courtemanche model showed that these kinetic parameters influence APD (Chang et al., 2017), and other studies have highlighted difficulties in calibrating ion channel dynamics using traditional approaches as well as showing that different formulations can have an important effect on the magnitude and time course of an ion channel current (Beattie et al., 2018). In the present study our focus was on the action potential rather than Ca²⁺ handling, and a detailed sensitivity analysis of the mechanisms of Ca²⁺ storage, release, and uptake in each model would be a valuable extension to the work presented here. So far a fully comprehensive analysis has only been done for specially constructed models (Pathmanathan et al., 2019). Nevertheless, a complete sensitivity analysis of biophysically detailed models, possibly using a hierarchical approach, remains an important challenge.

4.4.2. Simulator Instability

One of the issues with a complete sensitivity analysis, highlighted in Pathmanathan et al. (2019), is that parts of the simulator input space may generate implausible behaviors. For a cardiac cell model these behaviors might be a numerical instability, spontaneous beats, or failure to repolarize. In the present study we removed simulator runs from the design data where model behavior was implausible, or where the simulator runs produced action potential alternans. We considered this to be a pragmatic approach. However, it is clearly an area for improvement because the location of these regions of input space conveys information about the model, and approaches where these locations are encoded explicitly show promise (Ghosh et al., 2018). This consideration is especially important if the inputs have a greater range of variation than those considered in the present study, to represent the effects of cellular remodeling, pathological changes, or drug action. Extending the analysis presented here to a greater range of inputs is an important future direction.

4.4.3. Choice of Outputs

We selected a range of action potential features for our model outputs, these were based on measures used to describe experimental action potentials and aim to capture the main features of the action potential shape and duration. Our main focus was on the action potential. We included the maximum and minimum intracellular Ca²⁺ concentration as additional outputs, but did not consider the duration of the Ca²⁺ transient. We would not consider our choice of outputs to be definitive, and there may be better choices. A principal component analysis of the design data used in the present study showed that 95% of the output variance was accounted for by the first 6 principal components for the Courtemanche model and the first 5 for the Maleckar model. Parameterizing the action potential and Ca²⁺ transient so that they are described by a minimal set of features is likely to be important not only for model analysis but also for model calibration (Coveney and Clayton, 2018). Emulators that emulate time-dependent outputs have been developed, but are not yet widely used (Conti et al., 2009; Conti and O'Hagan, 2010), but could be a promising tool for extending work in this area.

4.4.4. Future Directions

The aim of this study was to demonstrate the utility of sensitivity and uncertainty analysis with GP emulators. A detailed and hypothesis-driven mechanistic study was beyond the scope of the present work, but will be a valuable next step. Extending the use of emulators from models of cardiac cells to models of cardiac tissue is an important future direction, and initial studies are promising (Lawson et al., 2018). At present, tissue calculations are computationally expensive, especially for personalized meshes. However, the need to evaluate uncertainty in model predictions for use in the clinical setting requires computationally efficient approaches, and we anticipate exciting developments in this area.

Data Availability Statement

The datasets generated for this study can be found in Clayton and Coveney (2019), Figshare Dataset—https://figshare.shef.ac.uk/articles/Design_data_for_sensitivity_analysis_of_human_atrial_myocyte_models/11309756.

Author Contributions

RC conceived the study and undertook the analysis. SC wrote all the software used for sensitivity and uncertainty analysis. Manuscript preparation and review was shared by both the authors.

Funding

This work was funded by the UK Engineering and Physical Sciences Research Council through grants EP/P010741/1 and EP/K037145/1.

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.

Acknowledgments

We were grateful to Dr. Eugene Chang for helpful discussions in the early stages of preparing this manuscript. This manuscript has been released as a Pre-Print at https://www.biorxiv.org/content/10.1101/818047v1 (Coveney and Clayton, 2019).

Supplementary Material

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

Mathematical Details: Details of the procedures used to fit and evaluate the Gaussian process emulators, and to calculate main effects and sensitivity indices.

References

Andrianakis, I., and Challenor, P. G. (2012). The effect of the nugget on Gaussian process emulators of computer models. Comput. Stat. Data Anal. 56, 4215–4228. doi: 10.1016/j.csda.2012.04.020

CrossRef Full Text | Google Scholar

Beattie, K. A., Hill, A. P., Bardenet, R., Cui, Y., Vandenberg, J. I., Gavaghan, D. J., et al. (2018). Sinusoidal voltage protocols for rapid characterisation of ion channel kinetics. J. Physiol. 596, 1813–1828. doi: 10.1113/JP275733

PubMed Abstract | CrossRef Full Text | Google Scholar

Bingham, D., Sitter, R. R., and Tang, B. (2009). Orthogonal and nearly orthogonal designs for computer experiments. Biometrika 96, 51–65. doi: 10.1093/biomet/asn057

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Chang, E. T. Y., Coveney, S., and Clayton, R. H. (2017). Variance based sensitivity analysis of IKr in a model of the human atrial action potential using Gaussian process emulators. Funct. Imaging Modell. Heart 10263, 249–259. doi: 10.1007/978-3-319-59448-4_24

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

Cherry, E. M., and Evans, S. J. (2008). Properties of two human atrial cell models in tissue: restitution, memory, propagation, and reentry. J. Theor. Biol. 254, 674–690. doi: 10.1016/j.jtbi.2008.06.030

PubMed Abstract | CrossRef Full Text | Google Scholar

Cherry, E. M., Hastings, H. M., and Evans, S. J. (2008). Dynamics of human atrial cell models: restitution, memory, and intracellular calcium dynamics in single cells. Prog. Biophys. Mol. Biol. 98, 24–37. doi: 10.1016/j.pbiomolbio.2008.05.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Clarke, R., Ressom, H. W., Wang, A., Xuan, J., and Liu, M. C. (2008). The properties of high-dimensional data spaces: implications for exploring gene and protein expression data. Nat. Rev. Cancer 8, 37–49. doi: 10.1038/nrc2294

PubMed Abstract | CrossRef Full Text | Google Scholar

Clayton, R., and Coveney, S. (2019). Design Data for Sensitivity Analysis of Human Atrial Myocyte Models, Figshare Dataset. Available online at: https://figshare.shef.ac.uk/articles/Design_data_for_sensitivity_analysis_of_human_atrial_myocyte_models/11309756 (accessed December 13, 2019).

Google Scholar

Clayton, R. H., Bernus, O., Cherry, E. M., Dierckx, H., Fenton, F. H., Mirabella, L., et al. (2011). Models of cardiac tissue electrophysiology: progress, challenges and open questions. Prog. Biophys. Mol. Biol. 104, 22–48. doi: 10.1016/j.pbiomolbio.2010.05.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Colman, M. A., Aslanidi, O. V., Kharche, S., Boyett, M. R., Garratt, C. J., Hancox, J. C., et al. (2013). Pro-arrhythmogenic effects of atrial fibrillation induced electrical remodelling–insights from 3D virtual human atria. J. Physiol. 17, 4249–4272. doi: 10.1113/jphysiol.2013.254987

CrossRef Full Text | Google Scholar

Conti, S., Gosling, J. P., Oakley, J. E., and O'Hagan, A. (2009). Gaussian process emulation of dynamic computer codes. Biometrika 96, 663–676. doi: 10.1093/biomet/asp028

CrossRef Full Text | Google Scholar

Conti, S., and O'Hagan, A. (2010). Bayesian emulation of complex multi-output and dynamic computer models. J. Stat. Plan. Infer. 140, 640–651. doi: 10.1016/j.jspi.2009.08.006

CrossRef Full Text | Google Scholar

Courtemanche, M., Ramirez, R. J., and Nattel, S. (1998). Ionic mechanisms underlying human atrial action potential properties: insights from a mathematical model. Am. J. Physiol. 275, H301–H321. doi: 10.1152/ajpheart.1998.275.1.H301

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Coveney, S., and Clayton, R. H. (2019). Sensitivity and uncertainty analysis of two human atrial cardiac cell models using Gaussian process emulators. bioRxiv 818047. doi: 10.1101/818047

CrossRef Full Text | Google Scholar

Eck, V. G., Donders, W. P., and Sturdy, J. (2016). Advances in 0D and 1D models for circulation. A guide to uncertainty quantification and sensitivity analysis for cardiovascular applications. Int. J. Numer. Methods Biomed. Eng. 32:e02755. doi: 10.1002/cnm.2755

CrossRef Full Text | Google Scholar

Fink, M., Niederer, S. A., Cherry, E. M., Fenton, F. H., Koivumaki, J. T., Seemann, G., et al. (2011). Cardiac cell modelling: observations from the heart of the cardiac physiome project. Prog. Biophys. Mol. Biol. 104, 2–21. doi: 10.1016/j.pbiomolbio.2010.03.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Ghosh, S., Gavaghan, D. J., and Mirams, G. R. (2018). Gaussian process emulation for discontinuous response surfaces with applications for cardiac electrophysiology models. arXiv 1805.10020.

Google Scholar

Grandi, E., Pandit, S. V., Voigt, N., Workman, A. J., Dobrev, D., Jalife, J., et al. (2011). Human atrial action potential and Ca²⁺ model: sinus rhythm and chronic atrial fibrillation. Circ. Res. 109, 1055–1066. doi: 10.1161/CIRCRESAHA.111.253955

PubMed Abstract | CrossRef Full Text | Google Scholar

Greenstein, J. L., Wu, R., Po, S., Tomaselli, G. F., and Winslow, R. L. (2000). Role of the calcium-independent transient outward current I(to1) in shaping action potential morphology and duration. Circ. Res. 87, 1026–1033. doi: 10.1161/01.RES.87.11.1026

PubMed Abstract | CrossRef Full Text | Google Scholar

Hodgkin, A., and Huxley, A. (1952). A quantitative description of membrane current and its application to conduction and excitation in nerve. J. Physiol. 117, 500–544. doi: 10.1113/jphysiol.1952.sp004764

PubMed Abstract | CrossRef Full Text | Google Scholar

Johnston, B. M., Coveney, S., Chang, E. T., Johnston, P. R., and Clayton, R. H. (2018). Quantifying the effect of uncertainty in input parameters in a simplified bidomain model of partial thickness ischaemia. Med. Biol. Eng. Comput. 56, 761–780. doi: 10.1007/s11517-017-1714-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Johnstone, R. H., Chang, E. T. Y., Bardenet, R., de Boer, T. P., Gavaghan, D. J., Pathmanathan, P., et al. (2015). Uncertainty and variability in models of the cardiac action potential: can we build trustworthy models? J. Mol. Cell. Cardiol. 96, 49–62. doi: 10.1016/j.yjmcc.2015.11.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Kennedy, M. C., and O'Hagan, A. (2001). Bayesian calibration of computer models. J. R. Stat. Soc. B Stat. Methodol. 63, 425–464. doi: 10.1111/1467-9868.00294

CrossRef Full Text | Google Scholar

Koivumäki, J. T., Seemann, G., Maleckar, M. M., and Tavi, P. (2014). In silico screening of the key cellular remodeling targets in chronic atrial fibrillation. PLoS Comput. Biol. 10:e1003620. doi: 10.1371/journal.pcbi.1003620

PubMed Abstract | CrossRef Full Text | Google Scholar

Krogh-Madsen, T., Abbott, G. W., and Christini, D. J. (2012). Effects of electrical and structural remodeling on atrial fibrillation maintenance: a simulation study. PLoS Comput. Biol. 8:e1002390. doi: 10.1371/journal.pcbi.1002390

PubMed Abstract | CrossRef Full Text | Google Scholar

Krogh-Madsen, T., and Christini, D. J. (2012). Nonlinear dynamics in cardiology. Annu. Rev. Biomed. Eng. 14, 179–203. doi: 10.1146/annurev-bioeng-071811-150106

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, L. A., Carslaw, K. S., Pringle, K. J., Mann, G. W., and Spracklen, D. V. (2011). Emulation of a complex global aerosol model to quantify sensitivity to uncertain parameters. Atmos. Chem. Phys. 11, 12253–12273. doi: 10.5194/acp-11-12253-2011

CrossRef Full Text | Google Scholar

Lee, L. A., Pringle, K. J., Reddington, C. L., Mann, G. W., Stier, P., Spracklen, D. V., et al. (2013). The magnitude and causes of uncertainty in global model simulations of cloud condensation nuclei. Atmos. Chem. Phys. 13, 8879–8914. doi: 10.5194/acp-13-8879-2013

CrossRef Full Text | Google Scholar

Maleckar, M. M., Greenstein, J. L., Giles, W. R., and Trayanova, N. a. (2009). K⁺ current changes account for the rate dependence of the action potential in the human atrial myocyte. Am. J. Physiol. Heart Circ. Physiol. 297, H1398–H1410. doi: 10.1152/ajpheart.00411.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

McDowell, K. S., Arevalo, H. J., Maleckar, M. M., and Trayanova, N. A. (2011). Susceptibility to arrhythmia in the infarcted heart depends on myofibroblast density. Biophys. J. 101, 1307–1315. doi: 10.1016/j.bpj.2011.08.009

PubMed Abstract | CrossRef Full Text | Google Scholar

McKay, M. D., Beckman, R. J., and Conover, W. J. (1979). Comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics 21, 239–245. doi: 10.1080/00401706.1979.10489755

CrossRef Full Text | Google Scholar

Melis, A., Clayton, R., and Marzo, A. (2017). Bayesian sensitivity analysis of a 1D vascular model with Gaussian process emulators. Int. J. Numer. Methods Biomed. Eng. 33:e2882. doi: 10.1002/cnm.2882

PubMed Abstract | CrossRef Full Text | Google Scholar

Mirams, G. R., Pathmanathan, P., Gray, R. A., Challenor, P., and Clayton, R. (2016). White paper: uncertainty and variability in computational and mathematical models of cardiac physiology. J. Physiol. 594, 1–33. doi: 10.1113/JP271671

CrossRef Full Text | Google Scholar

Muszkiewicz, A., Britton, O. J., Gemmell, P., Passini, E., Carlos, S., Zhou, X., et al. (2015). Variability in cardiac electrophysiology: using experimentally- calibrated populations of models to move beyond the single virtual physiological human paradigm. Prog. Biophys. Mol. Biol. 120, 115–127. doi: 10.1016/j.pbiomolbio.2015.12.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Niederer, S. A., Fink, M., Noble, D., and Smith, N. P. (2009). A meta-analysis of cardiac electrophysiology computational models. Exp. Physiol. 94, 486–495. doi: 10.1113/expphysiol.2008.044610

PubMed Abstract | CrossRef Full Text | Google Scholar

Noble, D. (1962). A modification of the Hodgkin-Huxley equations applicable to Purkinje fibre action and pacemaker potentials. J. Physiol. 160, 317–352. doi: 10.1113/jphysiol.1962.sp006849

CrossRef Full Text | Google Scholar

Nygren, A., Fiset, C., Firek, L., Clark, J. W., Lindblad, D. S., Clark, R. B., et al. (1998). Mathematical model of an adult human atrial cell the role of K currents in repolarization. Circ. Res. 82, 63–81. doi: 10.1161/01.RES.82.1.63

PubMed Abstract | CrossRef Full Text | Google Scholar

Oakley, J., and O'Hagan, A. (2002). Bayesian inference for the uncertainty distribution of computer model outputs. Biometrika 89, 769–784. doi: 10.1093/biomet/89.4.769

CrossRef Full Text | Google Scholar

Oakley, J. E., and O'Hagan, A. (2004). Probabilistic sensitivity analysis of complex models: a Bayesian approach. J. R. Stat. Soc. B Stat. Methodol. 66, 751–769. doi: 10.1111/j.1467-9868.2004.05304.x

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Saltelli, A., Aleksankina, K., Becker, W., Fennell, P., Ferretti, F., Holst, N., et al. (2019). Why so many published sensitivity analyses are false: A systematic review of sensitivity analysis practices. Environ. Modell. Softw. 114, 29–39. doi: 10.1016/j.envsoft.2019.01.012

CrossRef Full Text | Google Scholar

Saltelli, A., Ratto, M., Tarantola, S., and Campolongo, F. (2012). Update 1 of: sensitivity analysis for chemical models. Chem. Rev. 112, PR1–PR21. doi: 10.1021/cr200301u

PubMed Abstract | CrossRef Full Text | Google Scholar

Sánchez, C., Bueno-Orovio, A., Wettwer, E., Loose, S., Simon, J., Ravens, U., et al. (2014). Inter-subject variability in human atrial action potential in sinus rhythm versus chronic atrial fibrillation. PLoS ONE 9:e105897. doi: 10.1371/journal.pone.0105897

PubMed Abstract | CrossRef Full Text | Google Scholar

Sarkar, A. X., Christini, D. J., and Sobie, E. A. (2012). Exploiting mathematical models to illuminate electrophysiological variability between individuals. J. Physiol. 590, 2555–2567. doi: 10.1113/jphysiol.2011.223313

PubMed Abstract | CrossRef Full Text | Google Scholar

Sarkar, A. X., and Sobie, E. A. (2010). Regression analysis for constraining free parameters in electrophysiological models of cardiac cells. PLoS Comput. Biol. 6:e1000914. doi: 10.1371/journal.pcbi.1000914

PubMed Abstract | CrossRef Full Text | Google Scholar

Sobie, E. A. (2009). Parameter sensitivity analysis in electrophysiological models using multivariable regression. Biophys. J. 96, 1264–1274. doi: 10.1016/j.bpj.2008.10.056

PubMed Abstract | CrossRef Full Text | Google Scholar

Tobón, C., Ruiz-Villa, C. A, Heidenreich, E., Romero, L., Hornero, F., and Saiz, J. (2013). A three-dimensional human atrial model with fiber orientation. Electrograms and arrhythmic activation patterns relationship. PLoS ONE 8:e50883. doi: 10.1371/journal.pone.0050883

PubMed Abstract | CrossRef Full Text | Google Scholar

Varela, M., Colman, M. A., Hancox, J. C., and Aslanidi, O. V. (2016). Atrial heterogeneity generates re-entrant substrate during atrial fibrillation and anti-arrhythmic drug action: mechanistic insights from canine atrial models. PLoS Comput. Biol. 12:e1005245. doi: 10.1371/journal.pcbi.1005245

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilhelms, M., Hettmann, H., Maleckar, M. M., Koivumäki, J. T., Dössel, O., and Seemann, G. (2012). Benchmarking electrophysiological models of human atrial myocytes. Front. Physiol. 3:487. doi: 10.3389/fphys.2012.00487

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, H., Garratt, C. J., Zhu, J., and Holden, A. V. (2005). Role of up-regulation of IK1 in action potential shortening associated with atrial fibrillation in humans. Cardiovasc. Res. 66, 493–502. doi: 10.1016/j.cardiores.2005.01.020

PubMed Abstract | CrossRef Full Text | Google Scholar