Bayesian uncertainty analysis for complex systems biology models: emulation, global parameter searches and evaluation of gene functions

Fundamental challenges facing systems biology

Recent advances in genome sequencing techniques, a variety of ‘omic’ techniques and bioinformatic analyses, have led to an explosion of systems-wide biological data. Thus, identification of molecular components at the genome scale based on biological data has become possible. However, a major challenge in biology is to analyse and predict how functions in cells emerge from interactions between molecular components. Computational and mathematical modelling provide compelling tools to study the nonlinear dynamics of these complex interactions [1]. A particular example is kinetic modelling, in which the kinetics of each biological reaction are described in accordance with the corresponding biological process, and the properties of the whole system are described using differential equations: a common tool for analysing biological systems [2–4].

A critical problem found in the mathematical modelling of many complex biological systems, which is of particular severity in kinetic modelling, is that the models often contain large numbers of uncertain parameters (a common type being reaction rate parameters). In most cases, such kinetic parameters cannot be directly measured as experiments typically measure concentrations rather than rates. Even when such parameters can be measured ‘directly’, this is usually in experimental conditions that are significantly different from the cellular environment we wish to study. Therefore, we have to compare the mathematical model’s outputs with experimental observations, often in the form of measured concentrations and trends, and determine which values of the input or rate parameters will achieve an acceptable match between model and reality. This involves consideration of several sources of uncertainty including observation error, biological variability and the tolerance we place on the model’s accuracy, known as the model discrepancy. It is vital that we perform a global parameter search for all input parameter settings that achieve an acceptable match. This is because a single solution for the rate parameter values may suggest certain biological implications and give particular predictions for future experiments, both of which could be gravely misleading were we to explore the parameter space further and find several alternative solutions that give radically different implications and predictions. This is a mistake that is disturbingly common.

Unfortunately, performing global parameter searches over high dimensional spaces can be extremely challenging for several reasons, most notably: (a) the complex structure of the model and hence the complex way it imposes constraints on the parameters, (b) the substantial model evaluation time relative to the needs of the analysis, (c) the need for a careful assessment of an “acceptable match” that incorporates appropriately all the complexities and uncertainties of the comparison between the model and the real system, and (d) high dimensional spaces, being extremely large, require vast numbers of model evaluations to explore. For example, some spatial models of root development [5] require at least several minutes for a single evaluation. It is worth considering how large high dimensional spaces are: were we just to evaluate the model in question at the corners of the initial input space, in say 32 dimensions, we would require 2³²≃4.3 billion evaluations, which would take approximately 136 years if the model took 1 second per evaluation. However, global parameter searches are essential for any meaningful inference or prediction to be made about the biological system. Therefore this represents a fundamental challenge for the whole of systems biology. This article describes practical methodology to address this problem, based on Bayesian statistics methodology for the uncertainty analysis of complex models [6–9].

Bayesian emulation and uncertainty analysis

The issues surrounding the analysis of complex models under uncertainty, and specifically the global parameter search problem, are not unique to systems biology, and have been encountered in many different scientific disciplines. An area of Bayesian statistics has arisen to meet the demand of such analyses. This area, sometime referred to as the uncertainty analysis of computer models, centres around the construction of Bayesian emulators [6–9]. An emulator is a statistical construct that mimics the scientific model in question, providing predictions of the model outputs with associated uncertainty, at as yet unevaluated input parameter settings. The emulator is however, extremely fast to evaluate [10]. It provides insight into the model’s structure and, thanks to its speed, it can be used to help perform the global parameter search far more efficiently than approaches that just use the comparatively slow scientific model itself (for examples see [6, 8, 11–14]).

Many analyses and corresponding parameter searches still fail because an appropriate measure of an acceptable match between model and reality is not defined. This can lead to the use of badly behaved objective functions that do not properly capture the desired match criteria, and which are often harder to explore in high dimensions, due to increased numbers of ridges, spikes and local minima. The Bayesian emulation methodology we introduce naturally incorporates more detailed statistical models of the difference between the model outputs and the observed data, which allow the inclusion of important sources of uncertainty such as observational error and model discrepancy, the later being the upfront acknowledgement of the limitations of the current model. Various structures of increasing complexity are available for the representation of these uncertainties, depending on the requirements and importance of the study (see [6, 15–18] for examples and discussion).

It is worth noting that due to their speed, the use of emulators would greatly improve the efficiency of many forms of analysis that a modeller may wish to perform, e.g. for a fully Bayesian MCMC analysis [9, 19, 20] or for more direct global parameter searches such as [21]. However for high dimensional models, the particular strategy chosen for a parameter search is vital. Many approaches struggle due to being trapped in local minima (of which there may be many) or because they chase the scientifically spurious best match parameter setting. Here, we describe an efficient global parameter search method known as Bayesian history matching, which has proved very successful across a wide range of scientific disciplines including cosmology [6, 7, 15, 22–24], epidemiology [11, 25], oil reservoir modelling [8, 26–28], climate modelling [12], environmental science [16] and traffic modelling [29].

It utilises Bayesian emulators to reduce efficiently the input parameter space in iterations or waves, by identifying regions that are implausible as matches to the observed data, with the objective of identifying all acceptable input parameter settings. It is a careful approach that avoids many of the traps of common parameter search techniques.

Hormonal crosstalk network in Arabidopsis root development

Understanding how hormones and genes interact to coordinate plant growth is a major challenge in developmental biology. The activities of auxin, ethylene and cytokinin depend on the cellular context and exhibit either synergistic or antagonistic interactions. Previously, three of our authors developed a hormonal crosstalk network for a single Arabidopsis cell by iteratively combining modelling with experimental analysis [30]. Kinetic modelling was used to analyse how such a network regulates auxin concentration in the Arabidopsis root, by controlling the relative contribution of auxin influx, biosynthesis and efflux; and by integrating auxin, ethylene and cytokinin signalling [30]. Although some of the parameters in the model were based on experimental data, most parameters were chosen in an ad hoc way, by adjusting them to fit experimental data. Conditional on those somewhat ad hoc choices, it was shown that the hormonal crosstalk network quantitatively describes how the three hormones (auxin, ethylene, and cytokinin) interact via POLARIS peptide (PLS) [31, 32] to regulate plant root growth [30].

In this work we demonstrate the power of the Bayesian emulation methodology by applying it to the hormonal crosstalk network in Arabidopsis root development. Specifically, we explore the model’s 32-dimensional parameter space, and identify the set of all acceptable matches between model outputs and experimental data, taking into account major sources of uncertainty. This provides much insight into the model’s structure and the constraints imposed on the rate parameters by the current set of observed data. We apply the methodology to a second, competing model, and hence are able to investigate gene functions robustly. As an example, our analysis suggest that, in the context of the hormonal crosstalk network, POLARIS peptide (PLS) must have a role in positively regulating auxin biosynthesis.

The paper is organised as follows. In the “Methods” section we begin by defining a simple 1-dimensional toy model that we use to illustrate our definitions and to demonstrate the three main parts of the Bayesian methodology: linking the model to reality, Bayesian emulation, and history matching, before going on to compare the strengths and weaknesses of Bayesian history matching to more standard approaches. In the “Results” and “Discussion” section we describe in detail the application of this methodology to the full 32 dimensional Arabidopsis model, and discuss the relevant insights and biological implications obtained.

Simple 1-dimensional exponential example

We will temporarily assume the initial conditions are f ₀=f _t=0(x)=1. The system runs from t=0 to t=5 and we are at first interested in the value of f _t (x) at t=3.5. This mathematical model features an input or rate parameter x, which we wish to learn about. We do this using a measurement of the real biological system at t=3.5 which we denote z, which corresponds to, but is not the same as, the model output f _t=3.5(x). Note that usually, for models of realistic complexity, we would not have the analytic solution for f _t (x) given by Eq. (2). Instead we would resort to a numerical integration method to solve Eq. (1) that might require significant time for one model evaluation, ranging from less than a second to hours, days or even weeks, depending on the full complexity of the model [6, 12]. Such a computational cost, for a single evaluation of the model, means that a full global parameter search is computationally infeasible, especially when the model has many rate parameters and therefore a high dimensional input space, which may require vast numbers of evaluations to explore fully.

We typically begin the analysis by exploring the model’s behaviour for several different values of the unknown rate or input parameter x. Figure 1 (left panel), shows five evaluations of the model f _t (x) for different values of x between x=0.1 and 0.5, coloured red to purple respectively, with time on the x-axis. The measurement of the system is denoted z, and is represented as the black point in Fig. 1, with ± 3σ _e error bars representing observational error, defined precisely below. This measurement was made at t=3.5 shown as the vertical dashed line. The most important questions for the biologist at this point are: can the model match the observed data z at all, and if so, what is the entire set of input parameter choices that give rise to acceptable matches between model output and observed data? Figure 1 (right panel) represents this question as it now shows only f _t=3.5(x) but now represented purely as a function of the input parameter x on the x-axis, with the red to purple points consistent with those in the left panel. The observed data z is now represented as the solid black horizontal line, with the ± 3σ _e error bars as the horizontal black dashed lines. We see that there will be acceptable values of x approximately between 0.3 and 0.35.

Methods to address these three fundamental issues are described in the next three sections.

Model discrepancy and linking the model to reality

Most systems biology models are developed to help explain and understand the behaviour of corresponding real world biological systems. An essential part of determining whether such models are adequate for this task is the process of comparing the model to experimental data. As a comparison of this kind involves several uncertainties that cannot be ignored, it is therefore vital to develop a clearly defined statistical model for the link between systems biology model f(x) and reality z. This allows for a meaningful definition of an ‘acceptable’ match between a model run and the observed data. Here we describe a simple yet extremely useful statistical model for the link between the biological model and reality, that has been successfully applied in a variety of scientific disciplines, for example climate, cosmology, oil reservoir modelling and epidemiology [6, 8, 11, 12].

where we represent the errors as additive, although more complex forms could be used. Note that z, y and e here represent vectors of random quantities, which will reduce to scalar random quantities if there is only one quantity of interest. A common specification [6] that we will employ here is to judge the errors to be independent from y, and unbiased with expectation E(e)=0 and, for the scalar case, \(\text {Var}(e) = \sigma _{e}^{2}\).

where the ε is the model discrepancy: a vector of uncertain quantities that represents directly the difference between the model and the real system. Again we treat y, f, x ^∗ and ε as vectors of random quantities. A simple and popular specification [6] would be to judge that ε is independent of f(x ^∗), x ^∗ and e, with E(ε)=0 and, in the scalar case, \(\text {Var}(\epsilon) = \sigma _{\epsilon }^{2}\). In a multivariate setting, where f(x) describes a vector of outputs (for example, with each output labelled by time t), the vector ε may have an intricate structure, possessing non-zero covariances between components of ε. This could capture the related deficiencies of the model across differing time points. Various structures of increasing complexity are available (for examples see [6, 8, 9]), along with methods for specification of their components [6, 16].

While the explicit inclusion of the model discrepancy term ε is unfamiliar, it is now standard practice in the statistical literature for complex models [7, 9, 17, 33]. Furthermore, any analysis performed without such a term is implicitly conditioned with the statement “given the model is a perfect representation of reality for some value of the inputs x”, a statement that is rarely true. The model discrepancy allows us to perform a richer analysis than before as we can now include any extra knowledge we have about the model’s deficiencies to improve our modelling of reality y, through the joint structure of ε (see for example [6, 16]). This is especially important in prediction: as a simple example, if our model undershot every auxin output we have measured so far, we may suspect that it will undershoot future measurements of auxin also, and may wish to build this into our prediction for future y [33].

We can specify probabilistic attributes of ε a priori, or learn about them by comparing to observed data. For direct specification, there are often various simple experiments that can be performed on the model itself to obtain assessments of σ _ε and other aspects if necessary. For example, often models are run from exact initial conditions, so performing a set of exploratory model evaluations with the initial conditions appropriately perturbed would provide a lower bound on σ _ε . See [16] where several such assessment methods are demonstrated, for more details.

1-dimensional example

Figure 2 (left panel) shows the case for the simple 1-dimensional exponential example model when both the observation error e and model discrepancy ε are ignored. The model f(x) is given by the purple line, while the observed data z is given by the horizontal black line. Here only one value of x can be viewed as acceptable (coloured green) while all others are unacceptable (red). This particular value of x is not unique in that if we were to perform the measurement again, due to measurement error we would get a different value for z and hence for x. More importantly, if the model had a second output, say corresponding to a different time, that also depended on the same input x, we would be extremely unlikely to be able to match both outputs to their measurements as we would have to obtain exact matches simultaneously for precisely the same value of x. Inferences and predictions about the biological system made from this case, using this value of x are not trustworthy.

Figure 2 (right panel) shows the far more realistic situation where we include both observation error e (the black dashed lines represent z±3σ _e ) and model discrepancy ε (the purple dashed lines show f(x)±3σ _ε ). As we have taken into account both major types of uncertainty there is now a range of acceptable values for x (green points) with borderline/unacceptable points in yellow/red. If we were to consider additional outputs of the model, we still have a chance to match them simultaneously to data for a subset of the currently acceptable points. If on the other hand we cannot find any acceptable points x even given the uncertainties represented by e and ε, then we can state that the model is inconsistent with the observed data and therefore most likely based on incorrect biological principles. Further, inclusion of the observation error and model discrepancy terms often aids a global parameter search as they tend to smooth the likelihood surface (or comparable objective function), making it both easier to explore while simultaneously more robust. They also help reduce the temptation to chase the often scientifically misleading global minimum, such as the lone green point in Fig. 2 (left panel), instead suggesting that the identification of a set of acceptable input points is the appropriate goal for such a search (see the green points in Fig. 2 (right panel)).

Bayesian emulation of systems biology models

Many complex mathematical models have been developed and employed within the area of systems biology. Often these models have high dimensional input spaces in that they posses several input parameters, for example reaction rate parameters, that must be specified in order to evaluate the model. We represent the list of such inputs as the vector x, with individual inputs as x _k with k=1,…,d. The model may have any number of outputs, denoted as the vector f(x), with individual outputs as f _i (x) with i=1,…,q, the behaviour of which we want to investigate, possibly comparing some of these to observed data. For example, the index i may label the different times we are interested in, or the different chemical outputs of the model, or both. Most models are complex enough that they require numerical integration methods to solve, and hence take appreciable time to evaluate. This evaluation time can range anywhere from less than a second to minutes, hours or even days for highly sophisticated models: our approach is applicable in any of these cases, and adds more value as the dimensionality and evaluation time of the model increases.

A Bayesian emulator is a fast, approximate mimic of the full systems biology model. It gives insight into the structure of the model’s behaviour and can be used instead of the model in many complex calculations. The emulator gives a prediction of what the model’s output f(x) will be at a yet to be evaluated input point x, and additionally provides an associated uncertainty for that prediction (these are often expressed as posterior distributions, or simply expectations and variances in some cases). Critically an emulator is extremely fast to evaluate as it only requires a few matrix multiplications, and hence can be used to explore the input space more fully, as for example in a global parameter search.

with expectation zero and \(\text {Var} (w_{i}(x)) = \sigma ^{2}_{w_{i}}\), that represents the effects of the remaining inactive input variables [6].

The emulator, as given by Eq. (5), possesses various desirable features. The regression term, given by \(\sum _{j} \beta _{ij} g_{ij}(x_{A_{i}})\), is often chosen to represent say a third order polynomial in the active inputs. This would attempt to mimic the large scale global behaviour of the function f _i (x), and in many cases, will capture a large proportion of the model’s structure. (It is worth noting that reasonably accurate emulators can often be constructed just using regression models, for example using the lm() function in R. This can be a sensible first step, before one attempts the construction of a full emulator of the form given in Eq. (5).) The second term \(\phantom {\dot {i}}u_{i}(x_{A_{i}})\), the Gaussian process, mimics the local behaviour of f _i (x) and specifically its local deviations from the third order polynomial given by the regression terms. We can choose the list of active inputs \(x_{A_{i}}\) using various statistical techniques for example, classical linear model fitting criteria such as AIC or BIC [6]. A list of say p active inputs for a particular output f _i (x) means that we have reduced the input dimension from d to p dimensions, which can result in large efficiency gains. The small remaining effect of the inactive inputs is captured by the third term w _i (x) in Eq. (5).

We proceed by performing an initial set of carefully chosen model evaluations, often picked to be ‘space filling’, i.e. well spread out over the input space. For example we may use a maximin Latin hypercube design, an approximately orthogonal design which attempts to ensure there are no large holes in-between the run locations (see Fig. 3 and [34–36]). An n point Latin hypercube design is created by dividing the range of each input into n sub-intervals, and placing points to ensure there is only ever one point in each sub-interval (this can be done using the lhs() function in R [37]). Many such Latin hypercube designs are generated and the one with maximum minimum distance between points is chosen.

We then fit the emulator given by Eq. (5) to the set of model runs using our favourite statistical tools, guided by expert judgement. Specifically we would prefer a fully Bayesian approach if we required full probability distributions [9], and a Bayes Linear approach [38, 39], which we will describe below, if we required purely expectations, variances and covariances of f(x). We make certain pragmatic choices in the emulator construction process, for example, while we keep the regression coefficients β _ij uncertain, we may directly specify \(\sigma ^{2}_{u_{i}}\), \(\sigma ^{2}_{w_{i}}\) and θ _i a priori, or use suitable plugin estimates [6].

The emulators then provide an expectation and variance for the value of f(x) at an unexplored input point x. We can test the emulators using a series of diagnostics, for example checking their prediction accuracy over a new batch of runs [40]. See [10] for an introduction and [6, 7, 9] for detailed descriptions of emulator construction.

All quantities on the right hand side of Eqs. (8) and (9) can be calculated from Eqs. (5) and (6) combined with prior specifications for E(β _ij ), Var(β _ij ), \(\sigma ^{2}_{u_{i}}\), \(\sigma ^{2}_{w_{i}}\) and θ _i . Note that we could have used the entire collection of model outputs D={D ₁,D ₂,…,D _q } instead of just D _i in Eqs. (8) and (9), if we had specified a more complex, multivariate emulator [41].

\(\mathrm {E}_{D_{i}}(f_{i}(x))\) and \(\text {Var}_{D_{i}}(f_{i}(x))\) are used directly in the implausibility measures used for the global parameter searches described below.

1-dimensional example

We now demonstrate the construction of an emulator for the simple one dimensional exponential model. As there is only one output dimension, f(x) is now a scalar, so we drop the i index from Eqs. (5-9).

Figure 4 (left panel) shows output from such an emulator of the simple model defined by Eq. 1. We suppose that only n=5 runs of the model have been performed at the locations x ^(j)=0.1,0.2,0.3,0.4,0.5, which are shown as the purple points (these are the same as the five coloured points in Fig. 1). We therefore have the model output values

$$\begin{array}{*{20}l} D &= (f(x^{(1)}), f(x^{(2)}),\ldots,f(x^{(5)}))^{T} \\ &= (e^{0.1\times3.5}, e^{0.2\times3.5},\ldots,e^{0.5\times3.5})^{T} \end{array} $$

(10)

where again the output of interest has t=3.5.

We use a simplified form of the emulator given by Eq. (5), where we choose the polynomial terms β _j g _j (x _A ) to represent only a constant term β ₀. As we only have one input variable, there is no distinction between inactive and active variables so we also set w(x) to zero, and hence the emulator Eq. (5) reduces to

$$ f(x) \;\;=\;\; \beta_{0} + u(x) $$

(11)

For simplicity we treat the constant term β ₀ as known and hence set Var(β ₀)=0, and choose prior expectation E(β ₀)=β ₀=3.5, a value which we expect the function outputs to be approximately centred around. We specify the parameters in the covariance function for u(x) given by Eq. (6) to be σ _u =1.5 and θ=0.14 representing curves of moderate smoothness: this process will be discussed in more detail for the full Arabidopsis model.

All expectation, variance and covariance terms on the right hand side of Eqs. (8) and (9) can now be found using Eqs. (11), (6) and (10), for example,

$$\begin{array}{@{}rcl@{}} \mathrm{E}(f(x)) &=& \beta_{0} \end{array} $$

(12)

$$\begin{array}{@{}rcl@{}} \text{Var}(f(x)) &=& \sigma_{u}^{2} \end{array} $$

(13)

$$\begin{array}{@{}rcl@{}} \mathrm{E}(D) &=& (\beta_{0}, \ldots, \beta_{0})^{T} \end{array} $$

(14)

while Cov(f(x),D) is a row vector of length n with jth component

$$\begin{array}{@{}rcl@{}} \text{Cov}(f(x), D)_{j} &=& \text{Cov}(f(x), f(x^{(j)})) \\ &=& \sigma^{2}_{u} \text{exp}\left\{- \frac{\|x-x^{(j)}\|^{2}}{\theta^{2}}\right\} \end{array} $$

(15)

and similarly Var(D) is an n×n matrix with (j,k) element

$$\begin{array}{@{}rcl@{}} \text{Var}(D)_{jk} &=& \text{Cov}(f(x^{(j)}), f(x^{(k)})) \\ &=& \sigma^{2}_{u} \text{exp}\left\{- \frac{\|x^{(j)}-x^{(k)}\|^{2}}{\theta^{2}}\right\} \end{array} $$

(16)

We can now calculate the adjusted expectation and variance E _D (f(x)) and Var _D (f(x)) from Eqs. (8) and (9) respectively.

Figure 4 (left panel) shows E _D (f(x)) as a function of x as the blue line. We can see that it precisely interpolates the five known runs at outputs D, which is desirable as f(x) is a deterministic function. The blue line also gives a satisfactory estimate of the true function f(x)= exp(3.5x). The red pair of lines give the credible interval \(\mathrm {E}_{D}(f(x)) \pm 3 \sqrt {\text {Var}_{D}(f(x))}\) as a function of x. This defines a region between the lines that we believe is highly likely to contain the true function f(x). Another desirable feature of the emulator is that these credible intervals decrease to zero width at the five known run locations, as is appropriate for a deterministic function, as we precisely know the value of f(x) there (and because we have no inactive variables). Therefore when x is close to a known run we are more certain about the possible values of f(x), compared to when x is far from any such runs.

Figure 4 (right panel) shows an emulator as applied to a more complex 1-dimensional function. Here the true function is

$$ f(x) \;\;=\;\; 3 \,x \sin\left(\frac{5\pi(x-0.1)}{ 0.4}\right) $$

(17)

which has been simulated at only 10 input points evenly spread between x ⁽¹⁾=0.1 and x ⁽¹⁰⁾=0.5. Here the prior emulator specifications were as in the previous example, but with E(β ₀)=0, σ _u =0.6 and θ=0.06 allowing for functions with more curvature, centred around zero. As before the blue and red lines show E _D (f(x)) and \(\mathrm {E}_{D}(f(x)) \pm 3 \sqrt {\text {Var}_{D}(f(x))}\) as functions of x. The true function f(x) is given by the solid black line and it can be seen that it lies within the credible region for all x, only getting close to the boundary for x>0.5. This demonstrates the power of the emulation process: with the use of only 10 points the emulator accurately mimics a reasonably complex function with five turning points. We will demonstrate the effectiveness of emulators in higher dimensions for the main Arabidopsis model example.

History matching: an efficient global parameter search

Bayesian emulation is very useful in a variety of situations. As emulators are extremely fast to evaluate, they can replace the original model in any larger calculation, for example when designing future experiments. They can also provide much structural insight into the behaviour of the model. One of the most important applications of emulation is to the problem of performing a global parameter search. In this section we describe a powerful iterative global search method known as history matching, which has been successfully employed in a variety of scientific disciplines including galaxy formation [6, 7, 15, 22, 24], epidemiology [11, 25, 42, 43], oil reservoir modelling [8, 26–28], climate modelling [12], environmental science [16, 44] and traffic modelling [29]. Many of these applications involved models with substantial runtime, for which the process of emulation is vital.

History matching is designed to answer these questions. It proceeds iteratively and employs implausibility measures to determine parts of the input space that can be discarded from further investigation.

The numerator of Eq. (18) gives the distance between the emulator expectation \(\mathrm {E}_{D_{i}}(f_{i}(x))\) and the observation z _i , while the denominator standardises this quantity by all the relevant uncertainties regarding this distance: the emulator variance \(\text {Var}_{D_{i}}(f_{i}(x))\), the model discrepancy variance Var(ε _i ) and the observation error variance Var(e _i ). This structure is a direct consequence of Eqs. (3) and (4). A large value of I _i (x) for a particular x implies that we would be unlikely to obtain an acceptable match between f _i (x) and z _i were we to run the model there. Hence we can discard the input x from the parameter search if I _i (x)>c, for some cutoff c, and refer to such an input as implausible. We may choose the cutoff c by appealing to Pukelsheim’s 3-sigma rule [45], which is the powerful result that states that for any continuous, unimodal distribution, 95% of its probability must lie within ± 3σ, regardless of asymmetry, skew, or heavy tails, which suggests that a choice of c=3 could be deemed reasonable [6].

where Q represents the collection of all outputs, or some important subset of them (often we will only emulate a small subset of outputs in early iterations). A more robust approach would be to consider the second or third maximum implausibility, hence allowing for some inaccuracy of the emulators [6]. Also, multivariate implausibility measures are available (see [6] for details), but these require a more detailed prior specification, for example this requires covariances between different components of e and ε. Note that a low value of the implausibility I _M (x) does not imply that the input point x is ‘good’ or ‘plausible’ as it still may lead to a poor fit to outputs that have not been included in Q yet. Also, low implausibility at x may occur because of a high emulator variance \(\text {Var}_{D_{i}}(f_{i}(x))\) which once resolved following further runs of the model, may then lead to a high implausibility at x. Hence we refer to low implausibility inputs x as “non-implausible", consistent with the literature [6–8, 11, 24, 25, 28].

We then analyse the form of the non-implausible volume \(\mathcal {X}\), the behaviour of model evaluations from different locations within it and the corresponding biological implications. If the entire input space has been deemed implausible in step 6 b), this may imply that the model is inconsistent with the observed data, with respect to the specified uncertainties. This could be because the biological principles that underlie the model’s structure are incorrect, and hence remodelling is required. Or that we may have underestimated the observation errors, the model discrepancy or even the emulator uncertainty, although emulator diagnostics [40] combined with the choice of fairly conservative cutoffs should make the latter unlikely. Note that concluding that a far larger model discrepancy is needed is essentially stating that the model is highly inaccurate, and may, for example be judged unfit for purpose.

1-dimensional example

We now demonstrate the history matching process as applied to the simple 1-dimensional exponential example. Figure 5 (left panel) shows the emulator expectation and credible intervals as in Fig. 4, however now the observation z plus observed error has been included as the horizontal black solid and dashed lines respectively. Here we have set the model discrepancy to zero (σ _ε =0) and reduced the size of the observation errors σ _e for clarity. Also given are the implausibilities I(x) as represented by the colours on the x=axis: red, yellow and green for high (I(x)>3.5), borderline (3.5<I(x)<3) and low (I(x)<3) implausibility respectively.

The non-implausible space \(\mathcal {X}_{1}\) at wave 1 is the full initial range of the rate parameter x, which is 0.075<x<5.25. If we impose cutoffs of I(x)<3 then this defines the wave 2 non-implausible space \(\mathcal {X}_{2}\) as shown by the green region of the x-axis in Fig. 5 (left panel).

We then perform the second wave by designing a set of two more runs over \(\mathcal {X}_{2}\), reconstructing the emulator over this region, and recalculating the implausibility measure I(x). The results of this second wave are shown in Fig. 5 (right panel). It can be seen that the emulator is now highly accurate over the \(\mathcal {X}_{2}\) region and that the non-implausible region in green has been further reduced. As the emulator is now far more accurate than the corresponding observation error, we may stop the analysis with this wave as \(\mathcal {X}_{3} \simeq \mathcal {X}\), implying that further runs will do little to reduce the non-implausible region further. Note that providing we have enough runs in each wave, we would often create new emulators at each wave, defined only over the current green non-implausible region [6], instead of updating emulators from previous waves, as in Fig. 5. See the Additional file 1 for R code to reproduce the example model output, discrepancy, emulation and history matching plots of Figs. 1, 2, 4, and 5 respectively.

History matching and Bayesian MCMC

Here we discuss the standard form of a full Bayesian analysis, and compare it to the above history matching approach, highlighting the relative strengths and weaknesses of each method.

History matching attempts to answer efficiently some of the most important questions that a modeller may have, identifying if the model can match the data, and where in input space such acceptable matches can be found. It requires only a limited specification related to the key uncertain quantities, in terms of means, variances and covariances. A fully Bayesian approach goes further [46], and delivers a posterior distribution across all uncertain quantities, which has the benefit of providing probabilistic answers to most scientific questions e.g. in this context it gives the posterior distribution of the location of the true input x ^∗. However, it requires a more detailed prior specification of joint probability distributions across all these quantities, and critically, it also assumes the existence of a single true x ^∗ (and the accuracy of the statistical model that defines it). This may not be judged appropriate, for example, in a situation where two different x ^∗s were found that gave good matches to two different subsets of the observed data. In more detail, a fully Bayesian specification in the context of Eqs. (3) and (4) requires the multivariate distributions π(z|y),π(y|f(x ^∗)),π(f(x ^(j))|x ^(j)) for a collection of inputs x ^(j), and a prior distribution π(x ^∗) over the true input x ^∗. Meaningful specifications of this form can be difficult to make, for example, often familiar distributional forms are assumed such as the multivariate normal distribution, but such choices are often made for mathematical convenience or computational tractability. Really, choices of this kind demand a careful justification, without which results such as the posterior for x ^∗ rapidly lose meaning. For example, the fully Bayesian approach will, after substantial calculation, return a posterior for x ^∗, which may be quite narrow, even if the model cannot fit the observed data at all. History matching however may quickly discover this mismatch after a few waves, making further analysis unnecessary.

The second drawback of the fully Bayesian approach is that it is often hard to perform the necessary calculations, and therefore various numerical schemes are required, the most popular being MCMC [47]. While MCMC has enjoyed much success, issues still remain over the convergence of an MCMC algorithm for even modest dimensional problems [48]. Often the likelihood may be highly multimodal, and therefore vast numbers of model evaluations are usually required to reach convergence, making MCMC prohibitively expensive for models of even moderate evaluation time. In contrast, the calculations for history matching are relatively fast and simple.

A third issue is that of robustness: small changes in the full Bayesian specification, especially involving the likelihood, can lead to substantial changes in the posterior. These sensitivities can go unnoticed and can be hard to analyse [49–52], but will call into question the resulting scientific conclusions.

Due to these issues, we support the fully Bayesian approach, but only for cases where such detailed calculations are warranted, say for a well tested, accurate biological model, which possesses well understood model deficiencies, and which is to be combined with data that have a trusted observation error structure, and critically, where full probabilistic results are deemed essential. If, instead, the main concern is to check whether the model can fit the data; to see what regions of input parameter space give acceptable matches, and for this to be used for further model development, then a history match may be the more appropriate analysis. Even if one wishes to forecast the results of future experiments [33, 53] and to make subsequent decisions, history matching can be sufficient, as one can either re-weight appropriately the samples generated in the final wave, as was done in [11], or as stated in step 7 of the history matching algorithm, use more sophisticated sampling in the non-implausible region depending on the scientific question.

Now the use of emulators can of course facilitate the large number of model evaluations required for Bayesian MCMC algorithms, admittedly at the expense of increased uncertainty (see for example [9, 54]). However, a more serious problem is then encountered. The likelihood function, which is a core component of Bayesian calculations, is constructed from all outputs of interest (and therefore attempts to describe both the ‘good’ and ‘bad’ inputs simultaneously). Hence we need to be able to emulate with sufficient accuracy all such outputs, including their possibly complex joint behaviour. This may be extremely challenging as often, especially in early waves, there may exist several erratically behaved outputs that are extremely difficult to emulate, which will dramatically fail emulator diagnostics [40]. Unfortunately, the likelihood and hence the posterior may be highly sensitive to such poorly constructed emulators. Therefore, from a purely practical perspective, employing the full Bayesian paradigm using inadequate emulators constructed at wave 1, may be unwise [6, 15].

History matching as a precursor to MCMC

If one does wish to perform a fully Bayesian analysis on a well tested biological model, we would usually recommend performing a history match first [15]. This would identify the non-implausible region \(\mathcal {X}\) which should contain the vast majority of the posterior. Then MCMC or an equivalent Bayesian calculation (such as importance sampling) can be performed within \(\mathcal {X}\), using the accurate emulators defined in the final wave (for an example of this see [11]).

This is because if the model is of modest to high dimension, the posterior may often only occupy a tiny fraction of the original input space \(\mathcal {X}_{1}\). Unless the model is very fast to evaluate, we would need to use emulators to overcome the MCMC convergence issues, but performing enough model runs to construct sufficiently accurate emulators over the whole of \(\mathcal {X}_{1}\) would be extremely inefficient. Iterative history matching naturally provides the accurate emulators that we would need, defined only over \(\mathcal {X}\), which should contain the posterior.

Note that the history match cuts out space efficiently, based on small numbers of easy to emulate outputs in early waves, and designs appropriate runs for the next wave that have good coverage properties for the current non-implausible region. Alternative iterative strategies, such as using MCMC at each wave to generate samples from the current posterior (which includes high emulator uncertainty), for use as model runs in the next wave, may be highly inefficient and could run into a series of difficulties. Such strategies would not fully exploit the smoothness of the model output, and may tend to cluster points together around the current posterior mode (hence squandering runs, for smooth deterministic models), leading to poor coverage of the full non-implausible space. Such strategies may also be highly misled by inaccurate emulators and the subsequent posterior sensitivity combined with multimodal likelihood issues, leading to clustered designs in the wrong parts of the input space, and in some cases a lack of convergence.

To conclude, if one really desires a fully Bayesian analysis, then performing a history match first can greatly improve efficiency. In this way we view history matching not as a straight competitor to alternative approaches, but instead a complimentary technique, which has many benefits in its own right. See [6] and the extended discussion in [15] for more details of this argument.

History matching and ABC

Another Bayesian technique that has been developed more recently, and compared (somewhat cautiously) to history matching, by for example [55], is that of Approximate Bayesian Computation or ABC [56]. While the two approaches seem to share some superficial similarities, they are fundamentally different in their goal and the principled way each approach is set up and implemented. For example, ABC attempts to approximate full Bayesian inference and hence to obtain an approximate Bayesian posterior distribution (critically, using a tolerance that tends to zero). History matching is not an inference procedure, as it is simply attempting to rule out all the input space that is clearly inconsistent with the data given the model discrepancy and observation error (which are meaningful ‘tolerances’ that critically will never tend to zero). It is worth noting that if one attempts to specify a meaningful minimum size for the tolerance in ABC, one is arguably not really employing ABC anymore, but is instead just back using Bayesian inference (as shown by [56]) in the form of a sampling-resampling algorithm (as described for example by [57]). History matching does not attempt to probabilize the remaining input space in any way, which can result in increased efficiency of the parameter search. We have directly compared and contrasted the two approaches in [58], where we demonstrated that a powerful version of ABC failed to find any part of parameter space that matched the observed data for a 22-dimensional stochastic epidemiology model, while history matching found the correct part of input space and many good matches using approximately half the number of runs that the (failed) ABC-SMC approach used.