Empirical quantification of predictive uncertainty due to model discrepancy by training with an ensemble of experimental designs: an application to ion channel kinetics

TL;DR

This work proposes to characterise uncertainty owing to model discrepancy with an ensemble of parameter sets, each of which results from training to data from a different protocol, that will help select more suitable ion channel models for future studies and will be widely applicable to a range of biological modelling problems.

Abstract

When mathematical biology models are used to make quantitative predictions for clinical or industrial use, it is important that these predictions come with a reliable estimate of their accuracy (uncertainty quantification). Because models of complex biological systems are always large simplifications, model discrepancy arises - where a mathematical model fails to recapitulate the true data generating process. This presents a particular challenge for making accurate predictions, and especially for making accurate estimates of uncertainty in these predictions. Experimentalists and modellers must choose which experimental procedures (protocols) are used to produce data to train their models. We propose to characterise uncertainty owing to model discrepancy with an ensemble of parameter sets, each of which results from training to data from a different protocol. The variability in predictions from this ensemble provides an empirical estimate of predictive uncertainty owing to model discrepancy, even for unseen protocols. We use the example of electrophysiology experiments, which are used to investigate the kinetics of the hERG potassium ion channel. Here, 'information-rich' protocols allow mathematical models to be trained using numerous short experiments performed on the same cell. Typically, assuming independent observational errors and training a model to an individual experiment results in parameter estimates with very little dependence on observational noise. Moreover, parameter sets arising from the same model applied to different experiments often conflict - indicative of model discrepancy. Our methods will help select more suitable mathematical models of hERG for future studies, and will be widely applicable to a range of biological modelling problems.

Figures (10)

• Figure 1: Under model discrepancy, parameter estimates depend on the design used for training. Panels to the left of the dotted line correspond to designs containing $n=11$ observations at times ($T_1, \ldots, T_4$ as shown in panel a). Panels on the right show designs with $n=101$ observations, ($T_1', \ldots, T_4'$ as shown in panel b) (a) and (b): representative datasets generated by the DGP shown with the solid black line (Equation \ref{['eqn:toy_deterministic']}) with points indicating observations (sampled using Equation \ref{['eqn:toy_observations']}) and the fitted discrepant model (Equation \ref{['eqn:simple_example_discrepant_y']}), with calibrated ${\boldsymbol{\theta}}$) (grey dashed lines). (c) and (d): the parameter estimates for each design, each fitted to one of ten repeats of the DGP. (e) and (f): 99% Bayesian credible regions obtained using MCMC, a uniform prior and a single repeat of the DGP. (g) and (h): predictions using the discrepant model fitted using a single repeat of each protocol (using the estimates shown in e and f), showing the true DGP (black), discrepant model predictions (red), and the difference between predictions (grey).
• Figure 2: The structural differences between the two Markov model structures used in this paper for synthetic data generation and model training. (a): the four-state Beattie model used in both Case I and Case II. (b): the five-state Wang model used only for Case II. When a channel is in the open/conducting (O) state (green) current is able to flow. Whereas, when the model is in the other closed (C) or inactivated (I) states, no current can flow. The arrows adjacent to each model structure indicate the direction in which rates increase as the voltage increases.
• Figure 3: Left: a range of different input voltage-clamp protocols (forcing functions) used in this study. Right: corresponding synthetic output data IKr simulated using the Beattie model with noise added as described in Section \ref{['sec:computational_methods']}. Here, we generate and plot data observed at a 10 kHz sampling rate. Training protocols (all protocols except $d_0$) were tested for numerical identifiability fink_markov_2009: inverse problems performed on synthetic data with repeatedly sampled random noise yielded parameter estimates with little variability.
• Figure 4: The set of predictions (Equation \ref{['eqn:set_of_predictions']}) shown for parameter estimates obtained by training with different values of $\lambda$ to synthetic data under $d_1, \ldots, d_5$ (using the Beattie model). The synthetically generated data used for model validation are shown in grey and and the spread of the predictions is highlighted in yellow. (a): the voltage trace for $d_0$. (b): The set of predictions with $\lambda = \frac{1}{4}$. (c): the set of predictions with $\lambda = 1$, that is, under the assumption of the correct maximal conductance ($g$). (d): the set of predictions with $\lambda = 4$. N.B. the 'angular' nature of the current is not a plotting artefact, but reflects the fact the voltage clamp (a) is constructed from a series of linear ramps for compatibility with automated voltage clamp machines.
• Figure 5: Discrepancy in parameter estimates and subsequent currents when a non-discrepant model is fitted to synthetic data, with all parameters free except the maximal conductance, $g$, which is scaled by some factor $\lambda$, ($g=\lambda g^*$), where $g^*$ is the true value. (a): estimates of $\theta_1$ and $\theta_2$ obtained by training with different protocols for 10 repeats of the DGP. The lines (linearly interpolated using 17 values for $\lambda \in [ \frac{1}{4}, 4]$) show how the estimates from each protocol improve as $\lambda\rightarrow 1$. (b): $d_0$ voltage protocol. (c): the spread of predictions of $I_\text{Kr}$ under the $d_0$ protocol using the parameter estimates in Column a. (d): a heatmap showing the predictive error obtained by training and validating for each pair of protocols. Here c corresponds with the top row of each heatmap, as indicated.