Considering discrepancy when calibrating a mechanistic electrophysiology model

(a) Cardiac modelling

Cardiac models are typically a collection of mathematical functions governed by systems of ordinary and/or partial (when spatial dimensions are considered) differential equations, integrated using computational techniques, which produce responses that depend on the model inputs. Inputs can include model parameters, initial conditions, boundary conditions, and cellular, tissue or whole organ geometrical aspects. Inputs which have physiological meaning can sometimes be obtained by direct measurement, while others may need to be estimated via an indirect calibration procedure using experimental data. There are many examples of such cardiac models, at a variety of different scales, discussed in the papers of this special issue.

Mathematical modelling and computational simulation has been remarkably successful at providing insights into cardiac physiological mechanisms at cellular, tissue and whole organ scales [2–7]. In the majority of these quantitative efforts, models are derived based on simplified representations of complex biophysical systems and use in vitro and in vivo experimental data for calibration and validation purposes. Quantitative cardiac models have been a crucial tool for basic research for decades, and more recently have begun to transition into safety-critical clinical and pharmaceutical development applications [8–12]. The use of cardiac mathematical models in such applications will require high levels of credibility in the predictive model outputs, as well as an accurate quantification of the uncertainty in these predictions.

Parameters in cardiac models are often uncertain, mainly due to measurement uncertainty and/or natural physiological variability [13]. Thus, uncertainty quantification (UQ) methods are required to study uncertainty propagation in these models and help to establish confidence in model predictions. Parametric UQ is the process of determining the uncertainty in model inputs or parameters, and then estimating the resultant uncertainty in model outputs, thus testing the robustness of model predictions given our uncertainty in their inputs, and has been applied to a variety of cardiac models [14–19].

Another major source of uncertainty in modelling is uncertainty in the model structure, i.e. the form of the governing equations. There is always a difference between the imperfect model used to approximate reality, and reality itself; this difference is termed model discrepancy. Assessment of the robustness of model predictions given our uncertainty in the model structure, and methods to characterize model discrepancy, has received relatively little attention in this field (and mathematical/systems biology more generally). We have found only two published explicit treatments of discrepancy in cardiac electrophysiology models, in papers by Plumlee et al. [20,21]. In these studies, the assumption that ion channel rate equations follow an explicit form (such as that given, as we will see later, by equation (3.5)) was relaxed, and rates were allowed to be Gaussian processes (GPs) in voltage. A two-dimensional GP (in time and voltage) was then also added to the current prediction to represent discrepancy in current for a single step to any fixed voltage.

(b) Notation and terminology

Before discussing model discrepancy in detail, we introduce some notation and terminology. As the concepts introduced here are intended to be understood not just by a cardiac modelling audience, we provide a non-exhaustive list of terminology we have encountered in different fields to describe useful concepts relating to calibration and model discrepancy (and mathematical/computational modelling in general) in table 1.

Here, we delve into some of those concepts in more detail. Suppose a physiological system is modelled as y = f(θ, u), where f represents all governing equations used to model the system (also referred to as model form or model structure), θ is a vector of parameters characterizing the system, and u are known externally applied conditions or control variables applied in the particular experimental procedure. In a cardiac modelling context, these might represent a stimulus protocol, a drug concentration or the applied voltage protocol in a simulated voltage-clamp experiment. In general, θ = {θ_D, θ_C}, where values of θ_D are directly measured, and where values of θ_C are determined by calibration using the model f. Here, for simplicity of exposition, we assume θ_D is fixed (and known) and θ = θ_C.

We can distinguish between external conditions used for calibration, validation and prediction (that is, the application of the model, or context of use (CoU)), u_C, u_V, u_P, say. Suppose we have experimental data Y_C for calibration and Y_V for validation. A typical workflow, without UQ, is

—	Calibration: estimate $\hat{\mathit{\theta }}={\text{argmin}}_{\mathit{\theta }\in \mathbf{\Theta }}{d}_C(f(\mathit{\theta },u_C),Y_C)$, using some calibration distance function dC( · , · ) (e.g. a vector norm: $d_C(\mathbf{x},\mathbf{y})=||\mathbf{x}-\mathbf{y}||$), and some subset of parameter space Θ;
—	Validation: compare $y_V=f(\hat{\mathit{\theta }},u_V)$ against YV, either qualitatively or using a suitable validation distance $d_V(f(\hat{\mathit{\theta }},u_V),Y_V)$;
—	Context of use: compute $Y_P=f(\hat{\mathit{\theta }},u_P)$, or some quantity derived from this, to learn about the system or to make a model-based decision.

The calibration stage has many different names (table 1).

In practice, there a number of reasons why we cannot infer parameter values with certainty. The most commonly considered situation is when the link between the data and the model output is stochastic, e.g. because of measurement error on Y_C or because of model discrepancy. Computing the uncertainty about θ based on noisy data Y_C is referred to as ‘inverse UQ’, and requires a statistical model of the experimental data to be specified. For example, when considering measurement error, a common choice is to assume independent identically distributed zero-mean Gaussian errors on all data points, in which case (neglecting model discrepancy; see later) our model for the data is

with $\mathit{ϵ}={({ϵ}_1,{ϵ}_2,\dots )}^{\mathrm{\top }}$, where ${ϵ}_i\sim \mathcal{N}(0,{\sigma }^2)$. There are many different approaches to solving inverse UQ problems (e.g. [22,23]), most of which are based on inferring probability distributions to describe the relative likelihood that each different parameter set is consistent with the available data. Though a number of different methods to solve inverse UQ problems have been applied in cardiac electrophysiology [13], the most common is a Bayesian approach, which combines prior information about the parameters, π(θ), with the probability of observing the data given each parameter π(Y_C|θ) (referred to as the likelihood of θ), to find a posterior distribution over the parameters For an introduction to Bayesian methods, see [24,25]. For the i.i.d. Gaussian error model (equation (1.1)), the likelihood is given by where $||x|{|}_2^2=\sum _i{x}_i^2$, and n is the number of data points.

Another potential source of uncertainty about θ can occur when the parameter varies across the (or a) population. Estimating population variability in θ requires multiple Y_C recordings, $\{Y_C^{(1)},Y_C^{(2)},\dots \}$. Multilevel or hierarchical models can then be used: we assume the parameters for population i are drawn from some distribution θ⁽ⁱ⁾ ∼ π(θ|ψ), and infer the population parameters ψ, see [26].

Once uncertainty in θ (given the data) has been determined, the impact of this uncertainty on validation simulations Y_V or CoU simulations Y_P can be computed by propagating the uncertainty through the model f in the validation/CoU simulations, e.g.

This is referred to as ‘uncertainty propagation’ or the ‘posterior predictive distribution’. Uncertainty in the prediction of Y_V helps provide a more informed comparison to the observed validation data (especially if experimental error in Y_V is also accounted for). Uncertainty in Y_P enables a more informed model-based decision-making process.

(c) Model discrepancy

UQ as outlined above does not account for the fact that the model is always an imperfect representation of reality, due to limited understanding of the true data-generating mechanism and perhaps also any premeditated abstraction of the system. The model discrepancy is the difference between the model and the true data-generating mechanism, and its existence has implications for model selection, calibration and validation, and CoU simulations.

For calibration, the existence of model discrepancy can change the meaning of the estimated parameters. If we fail to account for the model discrepancy in our inference, our parameter estimates, instead of being physically meaningful quantities, will have their meaning intimately tied to the model used to estimate them (we end up estimating ‘pseudo-true’ values; see §2c). The estimated parameter values depend on the chosen model form, and the uncertainty estimates obtained during inverse parameter UQ tell us nothing about where the true value is (only how confident we are about the pseudo-true values). In other words, there is no guarantee the obtained θ will match true physiological values of any parameters that have a clear physiological meaning.

We can try to restore meaning to the estimated parameters by including a term to represent the model discrepancy in our models. Validation, in particular, provides an opportunity for us to identify possible model discrepancy. In many cases, validation, rather than being considered as an activity for confirming a ‘model is correct’, is better considered as a method for estimating the model discrepancy. To maximize the likelihood that the validation can discern model discrepancy, the validation data should ideally be ‘far’ from the calibration data, and as close to the CoU as possible.

(a) Model calibration

We calibrate the model using a train of five APs stimulated under a 1 Hz pacing protocol as the calibration data. Before attempting to do this fitting exercise, the appropriately sceptical reader might ask whether we are attempting to do something sensible. Will we get back information on all the parameters we want, or will we just find one good fit to the data among many equally plausible ones, indicating non-identifiability of the parameters?

To address these questions, we first look at inferring the parameters of the original Model T (as well as inferring the noise model parameter, σ). We use equation (1.1) with Gaussian noise giving the likelihood in equation (1.2), together with a uniform prior distribution from 0.1 × to 10 × the original parameters of Model T. We take two different approaches to calibration. Firstly, we find a point estimate using a global optimization algorithm [35] to find the optimal model parameters (with no estimate of uncertainty). Secondly, we approximate the full posterior distribution using Markov chain Monte Carlo (MCMC). All inference is done using an open-source Python package, PINTS [36], and simulations are performed in Myokit [37].

The results are shown in electronic supplementary material, figure S1. This exercise results in a narrow plausible distribution of parameters very close to the ones that generated the data, and we conclude that the model parameters are identifiable with the given data. Additionally, electronic supplementary material, figure S1 shows that when using samples of these distributions to make predictions, all of the forward simulations are very closely grouped around the synthetic data for the I_Kr block CoU.

We now attempt the fitting exercise using Model F (i.e. the misspecified model). The fitted model prediction (using the maximum a posteriori (MAP) parameter estimate) is shown in figure 2a. The agreement between the calibrated model output and the synthetic data would be considered excellent if these were real experimental data. Therefore, it is tempting to conclude that this calibrated model gives accurate predictions, and that the model discrepancy is minor. But can we trust the predictive power of the model in other scenarios based solely on the result we see in figure 2a?

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

(b) Discrepant model predictions

Interestingly, the calibrated Model F gives very accurate predictions for the 2 Hz pacing validation protocol (data that are not used to estimate the parameters), as shown in figure 2b. Such rate-adaptation predictions are used commonly as validation evidence for AP models. At this stage, we may be increasingly tempted to conclude that we have a good model of this system’s electrophysiology.

But if one now uses the model to predict the effect of drug-induced I_Kr block, the catastrophic results are shown in the bottom right panel of figure 2. The calibrated Model F fails to repolarize, completely missing the true I_Kr block response of a modest AP duration prolongation. This example highlights the need for thorough validation and the CoU-dependence of model validation, but also the difficulty in choosing appropriate validation experiments.

We can also quantify the uncertainty in parameter estimates and predictions while continuing to ignore the discrepancy in Model F’s kinetics. Again, we use equation (1.2) together with a uniform prior to derive the posterior distribution of the parameters. The marginals of the posterior distribution, estimated by MCMC, and the point estimates obtained by optimization are shown in figure 3. The posterior distribution is very narrow (note the scale), which suggests that we can be confident about the parameter values. The resulting posterior predictions, shown in figure 2c, give a very narrow bound. By ignoring model discrepancy we have become highly (and wrongly) certain that the catastrophically bad predictions are correct.

It is worth noting that all of the issues above arise from the fact that the model discrepancy was ignored during calibration. In the scenario of no model discrepancy, i.e. when fitting Model T to the data, none of the issues above occurred, as shown in electronic supplementary material, figure S1.

To conclude our motivation of this paper, we can see that neglecting discrepancy in the model’s equations is dangerous and can lead to false confidence in predictions for a new context of use. We discuss methods that have been suggested to remedy this in §3.

(c) A statistical explanation

To understand what happens when we fit an incorrect model to data, let us first consider the well-specified situation where the data generating process (DGP) has probability density function (pdf) g(y), and for which we have data y_i ∼ g( · ) for i = 1, …, n. Then suppose we are considering the models $\mathcal{P}=\{p(y\mid \mathit{\theta }):\mathit{\theta }\in \mathbf{\Theta }\}$, i.e. a collection of pdfs parameterized by unknown parameter θ. If the DGP g is in $\mathcal{P}$, i.e. we have a well-specified model so that for some θ₀ ∈ Θ, we have g( · ) = p( · |θ₀), then asymptotically, as we collect more data (and under suitable conditions [38]), the maximum-likelihood estimator converges to the true value θ₀ almost surely

or equivalently $p(\cdot \mid {\hat{\mathit{\theta }}}_n)$ converges to g( · ). Similarly, for a Bayesian analysis (again under suitable conditions [39]), the posterior will converge to a Gaussian distribution centred around the true value θ₀, with variance that shrinks to zero at the asymptotically optimal rate (given by the Cramér–Rao lower bound), i.e. where y_1:n = (y₁, …, y_n), and $\mathcal{I}({\mathit{\theta }}_0)$ is the Fisher information matrix for the true parameter value θ₀.

However, when our model is misspecified, i.e. $g\notin \mathcal{P}$ (there is no θ ∈ Θ for which g( · ) = f( · |θ)), if we do inference for θ ignoring the discrepancy, then we usually still get asymptotic convergence of the maximum-likelihood estimator and Bayesian posterior [40,41]. However, instead of converging to a true value (which does not exist), we converge to the pseudo-true value

where $\text{KL}(g||p)=\int g(x)\log (g(x)/p(x))\mathrm{d}x$ is the Kullback–Leibler divergence from p to g (a measure of the difference between two distributions). In other words, we converge upon the model, $p(\cdot \mid {\mathit{\theta }}^{\ast })$, which is closest to the DGP as measured by the Kullback–Leibler divergence (figure 4).

• Download figure
• Open in new tab
• Download PowerPoint

Perhaps more importantly from a UQ perspective, as well as getting a point estimate that converges to the wrong value, we still usually get asymptotic concentration at rate 1/n, i.e. the posterior variance shrinks to zero. That is, we have found model parameters that are wrong, and yet we are certain about this wrong value. The way to think about this is that the Bayesian approach is not quantifying our uncertainty about a meaningful physical parameter θ₀, but instead it gives our uncertainty about the pseudo-true value ${\mathit{\theta }}^{\ast }$. Consequently, we can not expect our calibrated predictions

to perform well, as we saw in the AP example above.

This leaves us with two options. We can either extend our model class $\mathcal{P}$ in the hope that we can find a class of models that incorporates the DGP (and which is still sufficiently simple that we can hope to learn the true model from the data), or we can change our inferential approach.

(a) Ion channel model example

We now illustrate the difficulty of accounting for model discrepancy in a tutorial example. We demonstrate that it can be hard to determine the appropriate information to include in δ, and that different functional forms for δ can lead to different parameter estimates.

We consider three structurally different models: Models A, B and C. We take Model C as the ground truth model in this particular example, and use it to perform synthetic voltage-clamp experiments and generate synthetic data. The goal is to use Models A and B to explain the generated synthetic data, assuming we have no knowledge about the ground truth Model C. This tutorial aims to demonstrate the importance of considering model discrepancy, jointly with model selection, to represent given data with unknown true DGP.

We use the hERG channel current as an example, and use three different model structures (shown in figure 5). Model A is a variant of the traditional Hodgkin–Huxley model, described in Beattie et al. [43]; Model B is used in Oehmen et al. [44]; and Model C is adapted from Di Veroli et al. [45].

• Download figure
• Open in new tab
• Download PowerPoint

All three ion channel models can be expressed using a Markov model representation. For a model with a state vector, x = (x₁, x₂, …)^T, then in each case x evolves according to

where M is the Markov matrix describing the transition rates between states. Markov models are linear coupled ordinary differential equations (ODEs) with respect to time, t, and states, x. Typically, the components in the Markov matrix, M, are nonlinear functions of voltage, V(t), which in these voltage-clamp experiments is an externally prescribed function of time known as the ‘voltage-clamp protocol’ (i.e. u_C in equation (1.1)). The observable, the macroscopic ionic current, I, measured under V(t), is where g is the maximum conductance, E is the reversal potential and $\mathcal{O}$ is the sum of all ‘open states’ in the model.

Take Model B as an example. Its state vector, x, and Markov matrix, M, can be written as

where x_i is the probability a gate is in state i (or equivalently, the proportion of gates which are in state i), with $\sum x_i=1$. The parameters k_i,j represent the transition rates from state x_i to state x_j. Note that for all our models, there is just one open state so that $\mathcal{O}=\text{O}$. For all three models, each transition rate, k_i,j, is voltage dependent and takes the form with two parameters (A_i,j, B_i,j) to be inferred. This yields a total of 12 parameters for Model B which we denote as {p₁, …, p₁₂}, together with the maximum conductance, g, to be found. Similarly for Model A, it has eight parameters {p₁, …, p₈} together with g, to be inferred.

(b) Synthetic experiments

We let Model C be the ground truth DGP and simulate data from it (using parameter values estimated from real room temperature data by Beattie et al. [43], where g = 204 nS). We add i.i.d. Gaussian noise with zero mean and standard deviation σ = 25 pA to the simulated data. We generate data under three different voltage-clamp protocols, V(t). These are a sinusoidal protocol (see top plot in figure 6) and an AP series protocol from Beattie et al. [43] (see electronic supplementary material, figure S9), and the staircase protocol from Lei et al. [26,46] (figure 6b).

• Download figure
• Open in new tab
• Download PowerPoint

(c) Standard calibration ignoring model discrepancy

To calibrate the model (without considering any model discrepancy), we assume a statistical model of the form of equation (1.1), which has the same observation noise model as our synthetic data. The likelihood of model parameter θ, having observed the data y = y_1:n, is given by equation (1.3).

We use the sinusoidal protocol (figure 6a) as the calibration protocol; the AP series protocol (electronic supplementary material, top, figure S9) and the staircase protocol (electronic supplementary material, figure S9) are used as validation data. We use a global optimization algorithm [35] to fit the model parameters using their maximum-likelihood estimates. All inference is done using PINTS [36].

The fitting results of Models A and B are shown in figure 6. Using different starting points in the optimization gives almost exactly the same parameter sets each time. Although both models fit the calibration data reasonably well, neither matches perfectly, due to model discrepancy. While the exact forms of the model discrepancy differs between the two models, both models notably fail to reproduce the correct form of the current decay following the step to −120 mVshortly after 2000 ms.

The validation predictions for the staircase protocol are also shown in figure 6. Unlike in the sinusoidal protocol, where Model A generally gives a better prediction than Model B, in the staircase protocol, it is more evident that the model discrepancy traits are different for each model. For example, Model B appears to give slightly better predictions of the current during the first 10 000 ms, whereas after this point Model A begins to give better predictions.

(d) Calibration with model discrepancy

We now consider an approach that allows us to incorporate model discrepancy when doing parameter inference and making predictions. We adapt the method proposed in [42] and instead of assuming independent errors in equation (1.1), which corresponds to assuming a diagonal covariance matrix for the vector of errors ϵ, we consider an additive discrepancy model of the form given by equation (3.1), giving a correlated (non-diagonal) error structure. We consider three different choices for the discrepancy δ(v_C), and jointly infer θ and δ. Note that we allow for a different choice of input v_C, compared to the input of model f, u_C.

First, we model δ as a sparse-GP [47,48], for which we adapted the implementation in PyMC3 [49] using Theano [50]. The radial basis function was used for the results presented here; we also tried two other GP covariance functions (the exponential covariance function and the Matérn 3/2 covariance function) in electronic supplementary material, §S7(c), where we found the impact of the choice of covariance functions in this problem is not as sensitive as the formulation of the discrepancy models. We explore two possibilities: choosing v_C to be either (i) t (time); or (ii) O,V (the open probability, $\mathcal{O}$ in equation (3.3), and the voltage, V). In fitting the model, we estimate hyperparameters associated with the GP covariance function, and condition the model on the observed discrepancies. For the GP(O, V) model, this means that we assume that the discrepancy is a function O and V so that we use the observed combinations of (O, V, δ) to predict future discrepancies; in the GP(t) model, it means that we assume the discrepancy process is always similarly distributed in time (which will not be a sensible assumption in most situations). Full details are provided in the electronic supplementary material, §S2.

As a third approach, we model discrepancy δ and the white noise error ϵ, as an autoregressive-moving-average (ARMA) model of order p, q [51]. If e_t = δ_t(v_c) + ϵ_t is the residual at time t, then an ARMA(p, q) model for e_t is

where ${\nu }_t\sim \mathcal{N}(0,{\tau }^2)$, and φ₁, …, φ_p and ζ₁, …, ζ_q are, respectively, the coefficients of the autoregressive and moving-average part of the model. We used the StatsModels [52] implementation, and assumed p = q = 2 throughout. Note that when using the ARMA model, we do not condition on the observed discrepancy sequence (so the mean of the ARMA process remains zero, unlike in the GP approaches), but only use it to correlate the discrepancy structure in time. In general, there is an interesting connection between GPs discretely sampled regularly in time, and autoregressive processes [47], but here we treat the ARMA process differently to how we use GP discrepancies, and use the data only to estimate the ARMA parameters, not to condition the process upon the observed temporal structure, i.e. we use the ARMA process as a simple approach for introducing correlation into the residuals to better account for the discrepancy, not to correct the discrepancy (as is done with the GP). The motivation is that if the mechanistic model is correct, the residuals should be uncorrelated, but for misspecified models they will typically be correlated. For further details, please refer to electronic supplementary material, §S3.

For all methods, i.i.d. noise, GP(t), GP(O, V) and ARMA(2, 2), we infer the posterior distribution of the parameters (equation (1.2)), where the priors are specified in electronic supplementary material, §S4. We use an adaptive covariance MCMC method in PINTS [32,36] to sample from the posterior distributions. The trace plots of the samples are shown in electronic supplementary material, §S7. The inferred (marginal) posterior distributions for Model A are shown in figure 7, and they are used to generate the posterior predictive distributions shown in figure 8. Electronic supplementary material, figure S16 shows the same plots for Model B. Note that the choice of the discrepancy model can shift the posterior distribution significantly, both in terms of its location and spread. In particular, the ARMA(2, 2) model gives a much wider posterior than the other discrepancy models.

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

Figure 8 shows the posterior predictive distributions of Model A with the calibration protocol using the four discrepancy models (electronic supplementary material, figure S17 for Model B), i.e.predicting the data used in training. The top panel shows the sinusoidal voltage protocol, and the panels underneath are calibrated model predictions with i.i.d. noise (blue), GP(t) (orange), GP(O, V) (green) and ARMA(2, 2) (red). The calibration data are shown in grey. Visually, we can see that the two GP models, GP(t) (orange) and GP(O, V) (green), fit the data with high accuracy; later we will see one of them is overfitting, while the other is not. The ARMA(2, 2) model (red) increases the width of the posterior (compared to i.i.d. noise), but its posterior mean prediction does not follow the data as closely as the two GP models.

Table 2 shows the root mean square errors (RMSEs) of the posterior mean predictions for all of the models, and is coloured so that yellow highlights the best performing model and red the worst. The first row of the table shows the results for the calibration (sine wave) protocol, and it is clear that the GP(t) and GP(O, V) models give the best RMSE values for the calibration data. Note that the RMSE only assesses the accuracy of the point estimate (given by the posterior mean). Table S1 in the electronic supplementary material gives the posterior predictive log-likelihoods; the log-likelihood is a proper scoring rule [53] which assesses the entire predictive distribution, not just the mean prediction. The ARMA(2, 2) and GP(O, V) models achieve the highest log-likelihood scores on the calibration data (best all round predictions when accounting for uncertainty).

Figure 9 shows the prediction results for the staircase validation protocol for Model A (electronic supplementary material, figure S18 for Model B) using different discrepancy models, with the same layout as figure 8. Similar figures for the AP protocol predictions are shown in electronic supplementary material, figures S9 (Model A) and S19 (Model B). The GP(t) discrepancy model is conditioned to give the same temporal discrepancy pattern as in the calibration protocal, and is unable to change its predicted discrepancy in any way for the validation protocol; i.e. the GP(t) discrepancy predicts as if it were still under the sinusoidal protocol. Thus, there is some residual from the calibration protocol shown in the GP(t) (orange) prediction for the staircase protocol, e.g. see ‘wobbly’ current at approxiametely 7000 ms as pointed at by the blue arrow.

• Download figure
• Open in new tab
• Download PowerPoint

For Model A, it is interesting to see that the RMSE of the point prediction (the posterior mean) in table 2 (top) is best for the i.i.d. noise model with the GP(O, V) model only a little worse. Note that the GP(O, V) model is able to capture and accurately predict the tail current after the two activation steps, as indicated by the red arrows in figure 9—a visible area of model mismatch in our calibration without model discrepancy. The uncertainty quantification in the predictions is poor for all of the discrepancy models, but from electronic supplementary material, table S1 we can see that when we assess the uncertainty in the prediction, the i.i.d. noise model is the worst performing model (as for intervals where the prediction is wrong, each error is equally surprising, whereas in correlated models, the first error in any interval makes subsequent errors more probable). The unstructured ARMA(2, 2) and GP(O, V) models score highest for their uncertainty quantification.

For Model B, the GP(O, V) discrepancy model gives the best overall predictions for both the staircase and the AP protocols, although when we examine the contributions of the mechanistic and discrepancy models, we see that an element of non-identifiability between them has arisen (electronic supplementary material, §S7b). In terms of the posterior predictive log-likelihood, electronic supplementary material table S1 (bottom) again highlights that the ARMA(2, 2) and GP(O, V) models tend to be better than the i.i.d. noise and GP(t) models.

Electronic supplementary material, figures S10, S11 and S12 show the model discrepancy for Model A for the sine wave protocol, AP protocol and staircase protocol, respectively; electronic supplementary material, figures S20, S21 and S22 show the same plots for Model B. Electronic supplementary material, figures S21 and S22, in particular, highlight that the GP(t) model has, by design, learnt nothing of relevance about model discrepancy for extrapolation under an independent validation protocol (in which V(t), and indeed the range of t differs from that of the training protocol). Furthermore, the discrepancy model is based only on information extending to 8000 ms (the duration of the training protocol), after which the credible interval resorts to the width of the GP prior variance. By contrast, the GP(O, V) model learns, independently of t, the discrepancy under combinations of (O, V) present in the training data (such as the activation step to 40 mV followed by a step to −120 mV), which is why it is able to better predict the tail current after the two activation steps. Finally, the ARMA(2, 2) model has zero mean with similar 95% credible intervals to the i.i.d. noise model, but has correlated errors and so scores better in terms of the posterior predictive log-likelihood. The ion channel (ODE) model-only predictions for the sine wave protocol, AP protocol and staircase protocol are shown in electronic supplementary material, figures S13, S14 and S15 for Model A and figures S23, S24 and S25 for Model B.

For a given dataset, the RMSE and log-likelihood values in table 2 and electronic supplementary material table S1 are comparable across models. Note that Model A is more accurate than Model B on all datasets and with all discrepancy models. With Model A, none of the discrepancy models are able to improve the mean predictions over the i.i.d. noise model performance, but the GP(O, V) comes close (in RMSE) while being able to capture some of the nonlinear dynamics that Model A misses, as discussed above. With Model B, the GP(O, V) model gives the best mean predictions (as measured by the RMSE). The GP(t) model achieves a better score on the calibration data, but by over-fitting the data. The ARMA(2, 2) model consistently gives the best posterior predictive log-likelihood values for Models A and B, as it gives a wider posterior distribution compared to other methods (figure 7). Over-confident predictions are heavily penalized by the log-likelihood, which explains the large differences observed in these values.

To conclude, we have used two different incorrect model structures (Models A, B) to fit synthetic data generated from a third model (Model C). We considered both ignoring and incorporating discrepancy when calibrating the model. Calibrating with discrepancy improved predictions notably for Model B, but not for Model A. Although our problem was a time-dependent (ODE) system, constructing the discrepancy model as a pure time-series based function is not necessarily useful in predicting unseen situations; we found the GP(O, V) model performed best at correcting the point prediction from the models.

(a) Recommendations

Very rarely do computational studies use more than one model to test the robustness of their predictions to the model form. We should bear in mind that all models are approximations and so when we are comparing to real data, all models have discrepancy. Here, we have seen, using synthetic data from an assumed true data-generating model, how fragile the calibration process can be for models with discrepancy, and how this discrepancy manifests itself in predictions of unseen situations. Synthetic data studies, simulating data from different parameter sets and different model structures, allow the modeller to test how well the inverse problem can be solved and how robust predictions from the resulting models are [54]. We strongly recommend performing such studies to learn more about the chosen, and alternative, models, as well as the effects of the model choice on parameter calibration and subsequent predictions. To develop our field further, it will be important to document the model-fitting process, and to make datasets and infrastructure available to perform and reproduce these fits with different models [55].

(b) Open questions and future work

The apparent similarity of the AP models we looked at (summarized in figure 1) is a challenge for model calibration. A number of papers have emphasized that more information can be gained to improve parameter identifiability with careful choice of experimental measurements, in particular by using membrane resistance [30,34], or other protocols promoting more information-rich dynamics [31,32] and some of these measurements may be more robust to discrepancy than others.

In synthetic data, fitting the model used to generate the data will recover the same parameter set from any different protocol (where there is sufficient information to identify the parameters). But in the presence of discrepancy, fitting the same model to data from different protocols/experiments will result in different parameter sets, as the models make the best possible compromise (as shown schematically in figure 4). This phenomenon may be an interesting way to approach and quantify model discrepancy.

If the difference between imperfect model predictions represented the difference between models and reality then this may also provide a way to estimate discrepancy. For instance, the largest difference between the ion channel Model A and B predictions in the staircase protocol was at the point in time that both of them showed largest discrepancy (figure 6). Some form of Bayesian model averaging [56], using variance-between-models to represent discrepancy, may be instructive if the models are close enough to each other and reality, but can be misleading if the ensemble of models is not statistically exchangeable with the DGP [57,58] or if there is some systematic error (bias) due to experimental artefacts [59].

In time-structured problems, rather than adding a discrepancy to the final simulated trajectory, as we have done here, we can instead change the dynamics of the model directly. It may be easier to add a discrepancy term to the differential equations to address misspecification, than it is to correct their solution, but the downside is that this makes inference of the discrepancy computationally challenging. One such approach is to convert the ODE to a stochastic differential equation [60,61], i.e. replace $\mathrm{d}\mathbf{x}/\mathrm{d}t=f_{\theta }(\mathbf{x},t)$ by $\mathrm{d}\mathbf{x}=f_{\theta }(\mathbf{x},t)\mathrm{d}t+{\mathrm{\Sigma }}^{1/2}\mathrm{d}{W}_t$ where W_t is a Brownian motion with covariance matrix Σ. This turns the deterministic ODE into a stochastic model and can improve the UQ, but cannot capture any structure missing from the dynamics. We can go further and attempt to modify the underlying model equations, by changing the ODE system to

where again δ(x) is an empirical term to be learnt from the data. For example, this has been tried with a discretized version of the equations using a parametric model for δ [62], with GPs [63], nonlinear autoregressive exogenous (NARX) models [64] and deep neural networks [65]. Computation of posterior distributions for these models is generally challenging, but is being made easier by the development of automatic-differentiation software, which allows derivative information to be used in MCMC samplers, or in variational approaches to inference (e.g. [66,67]).

Ultimately, modelling our way out of trouble, by expanding the model class, may prove impossible given the quantity of data available in many cases. Instead, we may want to modify our inferential approach to allow the best judgements possible about the parameters given the limitation of the model and data. Approaches such as approximate Bayesian computation (ABC) [68] and history-matching [69,70] change the focus from learning a statistical model within a Bayesian setting, to instead only requiring that the simulation gets within a certain distance of the data. This change, from a fully specified statistical model for δ to instead only giving an upper bound for δ, is a conservative inferential approach where the aim is not to find the best parameter values, but instead rule out only obviously implausible values [71,72].

For example, in the AP model from §2, instead of taking a Bayesian approach with an i.i.d. Gaussian noise model, we can instead merely try to find parameter values that get us within some distance of the calibration data (for details, see electronic supplementary material, figure S2). In the electronic supplementary material, we describe a simple approach, based on the methods presented in [73], where we find 1079 candidate parameter sets that give a reasonable match to the calibration data. When we use these parameters to predict the 2 Hz validation data, and the 75% I_Kr block CoU data, we get a wide range of predictions that incorporate the truth (electronic supplementary material, figure S3)—for a small subset of 70 out of 1079, we get good predictions and not the catastrophic prediction shown in figure 2. By acknowledging the existence of model discrepancy, the use of wider error bounds (or higher-temperature likelihood functions) during the fitting process may avoid fitting parameters overly-precisely. However, we have no way of knowing which subset of remaining parameter space is more plausible (if any) without doing these further experiments; testing the model as close as possible to the desired context of use helps us spot such spurious behaviour.

This paper has focused on the ion channel and AP models of cardiac electrophysiology. There is an audit of where uncertainty appears in cardiac modelling and simulation in this issue [74]. The audit highlights many other areas where discrepancy may occur: in assumptions homogenizing the subcellular scale to the models we have here; or at the tissue and organ scales in terms of spatial heterogeneity, cell coupling or mechanical models for tissue contraction and fluid-solid interaction. All of these areas need attention if we are to prevent model discrepancy producing misleading scientific conclusions or clinical predictions.