Understanding hormonal crosstalk in Arabidopsis root development via emulation and history matching

1.1 Use and understanding of scientific models in systems biology

One of the major challenges in biology is to understand how functions in cells emerge from molecular components. Computational and mathematical modelling is a key element in systems biology which enables the analysis of biological functions resulting from non-linear interactions of molecular components. The kinetics of each biological reaction can be systematically represented using a set of differential equations (Alves et al. 2006; Boogerd et al. 2007; Jamshidi and Palsson 2008; Liu et al. 2010; Smallbone et al. 2010). Due to the multitude of cell components and the complexity of molecular interactions, the kinetic models often involve large numbers of reaction rate parameters, that is parameters representing the rates at which reactions encapsulated by the model are occurring (Mobius and Laan 2015; Moore et al. 2015a, 2015c). Quantitative experimental measurements can be used to formulate the kinetic equations and learn about the associated rate parameters (Boogerd et al. 2007; Liu et al. 2010; Liu et al. 2013; Mobius and Laan 2015; Torres and Santos 2015). This in turn provides insight into the functions of the actual biological system.

An important question is therefore how much information about the kinetic equations and parameters can be obtained from an experimental measurement. Since a key aspect of experimental measurements in modern biological science is the study of the functions of specific genes, the answer to the above question is also important for understanding the role of each gene within the components of a biological system.

In plant developmental biology, a major challenge is to understand how plant growth is coordinated by interacting hormones and genes. Previously, a hormonal crosstalk network model – which describes how three hormones (auxin, ethylene and cytokinin) and the associated genes coordinate to regulate Arabidopsis root development – was constructed by iteratively combining regulatory relationships derived from experimental data and mathematical modelling (Liu et al. 2010, 2013; Moore et al. 2015a, 2015b, 2015c, 2017). However, for the mathematical model to best link with Arabidopsis root development, it is necessary to understand the parameter space of the model and identify all acceptable parameter combinations. Little is known about how acceptable parameter combinations of a model can be identified in light of specific experimental data. Therefore, this work explores how the acceptable parameter space of a complex model of hormonal crosstalk (Liu et al. 2013; Moore et al. 2015a, 2015b, 2015c, 2017) is assessed given a combination of quantitative experimental measurements and qualitative experimental trends by employing Bayesian history matching techniques (Craig et al. 1997; Cumming and Goldstein 2010; Gong and DiazDelaO 2017; Zhang et al. 2012). Additionally, we utilise these techniques to analyse how learning about the functions of a gene through particular relevant experiments can inform us about acceptable model parameter space.

1.2 Efficient analysis of scientific models

Complex systems biology models are frequently high dimensional and can take a substantial amount of time to evaluate (Moore et al. 2015c), thus comprehensive analysis of the entire input space, requiring vast numbers of model evaluations, may be unfeasible (Vernon et al. 2018). We are frequently interested, as is the case in this paper, in comparing the scientific model to observed data (usually measured with uncertainty), necessitating a possibly high dimensional parameter search. Our history matching approach aims to find the set of all possible combinations of input rate parameters which could have plausibly led to the observed data, given all sources of uncertainty involved with the model and the experimental data (Craig et al. 1997; Vernon et al. 2010a, 2018). This biologically relevant aim requires comprehensive exploration of the model’s behaviour over the whole input space, and therefore efficient techniques, such as emulation (Castelletti et al. 2012; Craig et al. 1997; Kennedy and O’Hagan 2001; Williamson et al. 2013), are required. An emulator mimics the biological model, but is substantially faster to evaluate, hence facilitating the large numbers of evaluations that are needed.

We are often keen to understand the contribution of particular sets of observations towards being able to answer critical scientific questions. Sequential incorporation of datasets into a history matching procedure, as presented in this article, is very natural and can allow us to attain such understanding. Comprehensive understanding and parameter searching of the hormonal crosstalk model for Arabidopsis root development (Liu et al. 2013), by sequentially history matching specific groups of experimental observations, is the focus of this paper.

2.1 Bayes linear emulation

In this section we review the process of constructing an emulator for a complex systems biology model. For more detail see Vernon et al. (2018). We represent the set of input rate parameters of the model as a vector x of length d, and the outputs of the model as vector f (x) of length q.

A Bayes linear emulator is a fast statistical approximation of the systems biology model built using a set of model runs, providing an expected value for the model output at a particular point x, along with a corresponding uncertainty estimate reflecting our beliefs about the uncertainty in the approximation (Goldstein 1999; Goldstein and Wooff 2007). The main advantage of emulation is its computational efficiency: often an emulator is several orders of magnitude faster to evaluate than the model it is mimicking. Emulation has been successfully applied across a variety of scientific disciplines such as climate science (Castelletti et al. 2012; Castruccio et al. 2014), cosmology (Bower et al. 2010; Heitmann et al. 2010), epidemiology (Farah et al. 2014), humanitarian relief (Overstall and Woods 2016), as well as systems biology (Vernon et al. 2018).

We index by i = 1,…,q the output components of the model. Each output component of the model f_i (x) can be represented in emulator form as presented by Vernon et al. (2010a):

where xAi represents the subset of active variables, that is the input components of x which are most influential for output f_i (x), g_ij are J_i known functions of xAi, and β_ij are the corresponding coefficients to the regression functions g_ij. ui(xAi) is a second-order weakly stationary stochastic process which captures residual variation in xAi. A priori, we assume that E[ui(xAi)]=0, along with the following covariance structure:

where SAi is the set of indices of the active inputs for output i. w_i (x) is a zero-mean “nugget”, or residual error, term with constant variance σwi2 over x and Cov[wi(x), wi(x′)]=0 for x ≠ x′. The nugget represents the effects of the remaining inactive input variables (Vernon et al. 2010a). We also make the assumption that:

for all i, j. There are several methods available for eliciting the parameters in Equation (2). We may have a priori beliefs as to sensible values of these parameters, perhaps as a result of expert knowledge of the general behaviour of the model, or from past computer experiments involving the same or a similar model. Alternatively, we may estimate the values of the parameters from the data, for example, using maximum likelihood estimates if we are happy with any implication of the specified prior distributions, or via predictive diagnostics such as Leave-One-Out-Cross-Validation (Andrianakis and Challenor 2011). Crucially, the resulting emulator should be assessed for adequacy using diagnostic measures, such as those discussed at the end of this section.

Suppose Di=(fi(x(1)) ,…, fi(x(n))) represents model output component i evaluated at n model runs performed at locations x⁽¹⁾,…., x⁽ⁿ⁾. The Bayes linear emulator output for simulator output component i at a new x is given by the Bayes linear update formulae (Goldstein 1999; Goldstein and Wooff 2007):

where the notation EDi[fi(x)] and VarDi[fi(x)] reflects the fact that we have adjusted our prior beliefs about f_i (x) by model runs D_i, and can be obtained for any point x using Equations (1) and (2). We note that in the literature it is common to assume normal and Gaussian process priors for β and u (x) in Equation (1) (Conti et al. 2009; Johnson et al. 2011; Kennedy and O’Hagan 2001), thus resulting in Gaussian process emulation (see Rasmussen and Williams (2006) for a general discussion of Gaussian processes). The resulting Bayesian update equations using Gaussian processes are practically similar to Equations (3) and (4) presented above, however, methodologically involve additional distributional assumptions. We would rather go as far as possible without making such distributional assumptions, as they may not always be valid, whilst still affecting the resulting inference. In this spirit, the Bayes linear framework is more similar to traditional kriging (Journel and Huijbregts 1978), noting that this term is now sometimes used to mean several different related approaches. Having said that, we note that kriging is derived from classical unbiased estimator arguments, whereas the Bayes linear paradigm follows from a foundational position, following DeFinetti (de Finetti 1974, 1975), that treats expectation as primitive and does not invoke concepts such as unbiasedness. The Bayes linear paradigm has been applied in a wide range of scenarios (Gosling et al. 2013; O’Hagan 1987; Whittle 1958): for example, in the context of the emulation of computer models, it has allowed tractable multilevel emulation due to multi-fidelity models, thus going beyond standard universal kriging (Craig et al. 1997; Cumming and Goldstein 2007).

Emulator design is the process of selecting the points in the input space x⁽¹⁾,…, x⁽ⁿ⁾ at which the simulator will be run in order to construct an emulator (Santner et al. 2003). A popular design choice in the computer model literature is the Maximin Latin Hypercube design (Currin et al. 1991; McKay et al. 1979), however, other options are also available (see, for example, Fisher (1937); Montgomery (2009)).

The quality of any emulator should be assessed using diagnostics to judge whether it is fit for purpose (Bastos and O’Hagan 2008). For example, we can calculate standardised prediction errors:

for a set of validation data. Large errors U_i (x) indicate conflict between simulator and emulator as a result of the emulator being overconfident in its predictions, and hence the emulator is not valid for inference. Systematically small errors indicate that the emulator is underconfident. As a rule of thumb, we expect most of the errors to lie between −3 and 3, appealing to Pukelsheim’s 3σ rule, which states that at least 95% of the probability mass of any unimodal continuous distribution will lie within ±3 standard deviations of the mean, regardless of asymmetry or skew (Pukelsheim 1994). To contrast with a Gaussian process emulator, these errors should theoretically follow a (standard) normal or a t-distribution, depending on the precise method of emulator update and inference (for example, approach to selecting the correlation length parameters). Poor diagnostics suggest that the emulator prior beliefs were misspecified, for example, as a result of incorrect prior specifications for the parameters β, σu2 and θ. Alternatively, it could be indication of an erratically behaved model that would require substantially more model runs in order to be emulated well.

2.1.1 Dimensional example

In this section we demonstrate emulation techniques on a simple one-dimensional example. We will suppose that we wish to emulate the simple function f(x)=0.1x+cos(x) in the range [0, 11π3], where we treat x as a rate parameter that we wish to learn about, and f (x) as a chemical concentration that we could measure. We assume an emulator of the form given by Equation (1) with covariance structure given by Equation (2). We assume a zero mean function, that is, gT(x)β=0, and that σu2=0.5, θ = 1.5 and σw2=0. We also specify a prior expectation E[f(x)]=0. Having specified our prior beliefs, we then use the update rules given by Equations (3) and (4) to obtain the adjusted expectation ED[f(x)] and variance VarD[f(x)] for f (x). The results of this emulation process are shown in the left panel of Figure 1. The blue lines represent the emulator expectation ED[f(x)] of the simulator output for the test points. The red lines represent the emulator mean ±3 emulator standard deviations, given as ED[f(x)]±3VarD[f(x)], these being bounds for a 95% credible interval, following Pukelsheim’s 3σ rule (Pukelsheim 1994). By comparison with the right panel of Figure 1, we can see that the emulator estimates the simulator output well, with some uncertainty. We note that we would not expect such large emulator uncertainty on such a smooth function as this, but have deliberately ensured that there is a large uncertainty for illustrative purposes, and in particular to highlight the effects of additional runs on reducing emulator uncertainty in the continuation of this example in Section 2.2.

Figure 1:

Left: Emulator expectation ED[f(x)] (blue) with emulator credible intervals ED[f(x)]±3VarD[f(x)] (red) for an emulator of f(x)=0.1x+cos(x) constructed using 8 training points. Right: Simulator function f(x)=0.1x+cos(x).

2.2 History matching

History matching concerns the problem of finding the set of inputs to a model for which the corresponding model outputs give acceptable matches to observed historical data, given our state of uncertainty about the model itself and the measurements. History matching has been successfully applied across many scientific disciplines including oil reservoir modelling (Craig et al. 1996, 1997; Cumming and Goldstein 2009, 2010; Oliver and Chen 2011), engineering (Gardner et al. 2018; Gong and DiazDelaO 2017), epidemiology (Andrianakis et al. 2015, 2017a, 2017b; McCreesh et al. 2017), climate modelling (Williamson et al. 2013) and systems biology (Vernon et al. 2018). Here we provide a brief summary of the history matching procedure (see Vernon et al. (2010a, 2018 for more details).

We need a general structure to describe the link between a complex model and the corresponding physical system. We use the direct simulator approach, otherwise known as the best input approach (Goldstein and Rougier 2006), where we posit that there exists a value x^⋆ such that f (x^⋆) best represents the real biological system, which we denote as y (Goldstein and Rougier 2006, 2009). We then formally link the ith output of the model to the ith real system value y_i via

and link the experimental observation z_i corresponding to output i to the real system via

where we assume fi(x⋆)╨ϵi╨ei, with a╨b indicating that random variables a and b are uncorrelated (Goldstein and Wooff 2007). Here, fi(x⋆) is the model run at best input x^⋆, ϵ_i is a random variable which reflects our uncertainty due to discrepancy between the model run at the best possible input combination setting and the real world (Arendt et al. 2012; Brynjarsdottir and O’Hagan 2014; Goldstein et al. 2013; Kennedy and O’Hagan 2001), and e_i is a random variable which incorporates our beliefs about the error between each desired real world quantity and its corresponding observed measurement. We assume E[ϵi]=E[ei]=0, Var[ϵi]=σϵi2 and Var[ei]=σei2. The connection between system, observation and model given by (6) and (7) is simple but well-used (Andrianakis et al. 2015; Craig et al. 1997; Goldstein and Rougier 2006), and judged sufficient for our purposes. For discussion of more advanced approaches see Goldstein and Rougier (2009).

We then aim to find the set X⋆ of all input combinations x that are consistent with Equations (6) and (7), that is those that will provide acceptable matches between model output and data. Note that classifying points in this way can lead to X⋆ being empty, an informative conclusion which would contradict the posited existence of x^⋆ in Equation (6), and imply that the model may not be fit for purpose. To analyse whether a point x∈X⋆ it is practical to use implausibility measures for each output i, as given, for example, in Craig et al. (1996, 1997; Vernon et al. (2010a):

If I_i (x) is large this suggests that we would be unlikely to obtain an acceptable match between model output and observed data were we to run the model at x. This is after taking into account all the uncertainties associated with the model and the measurements. We develop a combined implausibility measure over multiple outputs such as IM(x)=maxiIi(x), I2M(x)=maxi({Ii(x)}∖IM(x)) and I3M(x)=maxi({Ii(x)}∖{IM(x), I2M(x)}) (Vernon et al. 2010a), where S₁∖S₂ is general set notation for in set S₁ but not set S₂. We class x as implausible if the values of these measures lie above suitable cutoff thresholds (Craig et al. 1997; Vernon et al. 2010a).

History matching using emulators proceeds as a series of iterations, called waves, discarding regions of the input parameter space at each wave. At the kth wave emulators are constructed for a selection of well-behaved outputs Q_k over the non-implausible space remaining after wave k − 1. These emulators are used to assess implausibility over this space where points with sufficiently large values are discarded to leave a smaller set Xk remaining (Vernon et al. 2010a, 2018).

The history matching algorithm is as follows:

1. Let X0 be the initial domain space of interest and set k = 1.

2. Generate a design for a set of runs over the non-implausible space Xk−1, for example using a maximin Latin hypercube with rejection (Vernon et al. 2010a).

3. Check to see if there are new, informative outputs that can now be emulated accurately and add them to the previous set Q_k − 1 to define Q_k.

4. Use the design of runs to construct new, more accurate emulators defined only over Xk−1 for each output in Q_k.

5. Calculate implausibility measures over Xk−1 for each of the outputs in Q_k.

6. Discard points in Xk−1 with I (x) > c to define a smaller non-implausible region Xk.

7. If the current non-implausible space Xk is sufficiently small, go on to step 8. Otherwise repeat the algorithm from step 2 for wave k + 1. The non-implausible space is sufficiently small if it is empty or if the emulator variances are small in comparison to the other sources of uncertainty (σϵ2 and σe2), since in this case more accurate emulators would do little to reduce the non-implausible space further.

8. Generate a large number of acceptable runs from Xk, sampled according to scientific goal.

It should be the case that X⋆⊆Xk⊆Xk−1 for all k, where X⋆={x: maxiIi(x)<c} for some threshold c, where each I_i (x) is calculated using expression (8) with E[fi(x)]=fi(x) and Var[fi(x)]=0, that is were we to know the simulator output everywhere. The choice of cutoff c = 3 is frequently chosen, motivated by Pukelsheim’s 3-sigma rule (Pukelsheim 1994), which in this case implies that P(Ii(x)<3|x=x⋆)>0.95 for any unimodal continuous distribution for the combined error term ϵi+ei+(fi(x)−E[fi(x)]). This iterative procedure is powerful as it quickly discards large regions of the input space as implausible based on a small number of well behaved (and hence easy to emulate) outputs. In later waves, outputs that were initially hard to emulate, possibly due to their erratic behaviour in uninteresting parts of the input space, become easier to emulate over the much reduced space Xk. Careful consideration of the initial non-implausible space X0 is important. It should be large enough such that no potentially scientifically interesting input combinations are excluded. A more in-depth discussion of the benefits of this history matching approach, especially in problems requiring the use of emulators, may be found in Vernon et al. (2018).

Note that whilst in this article we take EDi[fi(x)] and VarDi[fi(x)] to be adjusted beliefs resulting from a Bayes linear emulator, such as those discussed in Section 2.1, these quantities can be substituted with the mean and variance resulting from the corresponding estimates of a Gaussian process emulator if preferred (see, for example, Hamdi et al. (2017) and Gardner et al. (2018)). In this case, history matching can proceed practically similarly. The main difference lies in the choice of implausibility cutoff threshold c, which will be lower in comparison to that used in a Bayes linear framework whilst ensuring the same upper bound on the probabilities of a false rejection (say around 2 as opposed to 3 above). This lower threshold can be utilised by appealing to the features of the spread of probability mass of the predictions assumed by making use of Gaussian processes, as opposed to general rules such as Pukelsheim’s 3σ rule discussed above. A brief comparison of Bayes linear emulators and Gaussian process emulators was presented in Section 2.1. For further discussion and a comparison of using Bayes linear emulators and Gaussian process emulators within a history match, see Vernon et al. (2010c).

2.2.1 Dimensional example continued

Figure 2 (left panel) shows the emulator expectation and conservative bounds for the 95% credible intervals, as given by Figure 1 (left panel), however, now an observation z = −0.3 along with observed error is included as solid and dashed lines respectively. In this example, we let model discrepancy be 0, and set the measurement standard error σe=0.05. Along the bottom of the figure is the implausibility I (x) for each x value represented by colour: red for large implausibility values, orange and yellow for borderline implausibility, and green for low implausibility (I (x) < 3) (Pukelsheim 1994).

Figure 2:

Left panel: Emulators for the simple 1D example f(x)=0.1x+cos(x) as given by the left panel of Figure 1. The blue line represents the emulator’s updated expectation ED[f(x)] and the pair of red lines give the credible interval ED[f(x)]±3VarD[f(x)]. Observation z along with observed error are shown as horizontal black solid and dashed lines respectively. The implausibilities I (x) are represented by the colours on the x-axis, with red representing high implausibility, orange and yellow representing borderline implausibility, and green representing low implausibility (I (x) < 3). Right panel: The wave 2 emulator for the same function, now including three additional wave 2 runs.

X0 is the full initial range of x, that is 0≤x≤11π3. X1 is as shown by the green regions in Figure 2 (left panel). Wave 2, shown in Figure 2 (right panel), involves designing a set of three more runs over X1, constructing another emulator over this region and calculating implausibility measures for each x∈X1. This second emulator is considerably more accurate than the observed error, thus X2≈X⋆, so the analysis can be stopped at this point as extra runs would do little to further reduce the non-implausible space.

2.3 Sequential history matching of observations

Part of the novelty of our history matching approach involves dividing experimental observations into subsets and sequentially performing a history match on the model using each group of observations. Much scientific insight can be gained from performing a history match, however, using all output components simultaneously can mask which experiments are informative for certain aspects of the scientific system.

Breaking the data down into subsets and sequentially adding them to the history match is a novel approach which allows for further scientific insight. Most prominently, it not only allows inferences to be made about the system quantities associated with the model input parameters, but also provides insight into the links between quantities associated with both the input and output. Note that adding model outputs sequentially in this way is different from bringing outputs in sequentially due to emulator capability (step 3 of the algorithm) (Vernon et al. 2010a). We will explore this in detail for the Arabidopsis model.

2.4 History matching and Bayesian inference

In this section we discuss how history matching is informative for aiding the understanding of physical systems described by a computer model, and can also be used in conjunction with alternative methods of analysis, such as a general Bayesian analysis, Bayesian optimisation and Approximate Bayesian Computation (ABC). In doing so we will compare some of the similar and differing features between history matching and the standard form of a general Bayesian analysis.

History matching is a computationally efficient and practical approach to identifying if a model is consistent with observed data, and, if so, utilising the key uncertainties within the problem to identify where in the input space acceptable matches lie (Craig et al. 1997). History matching is applicable for situations where the observed data is provided as either qualitative experimental trends or quantitative experimental measurements (or a combination of the two such as in the Arabidopsis model application in the following sections), attempting to answer some of the main questions that a modeller may have. In contrast, a general Bayesian framework typically requires full probabilistic specification of all uncertain quantities, providing a theoretically coherent method to obtain probabilistic answers to scientific questions. For example, in the context of a direct simulator, as given by Equations (6) and (7), a general Bayesian analysis will provide a posterior distribution for the location of the true best input x^⋆, whereas a history match provides a set, which may be empty, of points that could not implausibly be x^⋆ under the posited assumption that x^⋆ exists. History matching therefore has the benefit of avoiding the challenging task of making a full joint distributional specification over all uncertain quantities, to which resulting analyses may be sensitive.

Regardless of how prior distributions have been specified, performing the necessary calculations for a full Bayesian analysis is hard and typically computationally intensive, for example by making use of time-consuming numerical schemes such as Markov Chain Monte Carlo (MCMC) (Brooks et al. 2011). Many model evaluations are required to thoroughly explore the multi-modal likelihoods over the entire input space. Emulators can facilitate these large numbers of model evaluations at the cost of uncertainty (Higdon et al. 2008; Kennedy and O’Hagan 2001). However, since the likelihood function is constructed from all outputs of interest, we need to be able to emulate with sufficient accuracy all such outputs, including their possibly complex joint behaviour. Erratically behaved outputs may lead to emulators with large uncertainty or emulators which fail diagnostics. The consequential likelihood, and hence posterior, may be sensitive to the misspecification of these emulators.

Another related issue is that the posterior distribution of a general Bayesian analysis may be concentrated over a very small subspace of the initial input domain of interest X0. This is particularly true for models with relatively high-dimensional parameter spaces, and for which the initial ranges of interest for the parameters (that is, those defining X0) are large, such as for the Arabidopsis model discussed in the following sections. In such cases, a sufficiently accurate emulator over the whole input space will still require far too many model evaluations, hence alternative approaches to performing analysis, often sequential in nature and also involving the use of emulators, must be used. A variety of approaches are proposed in the literature. One group of such approaches is to use more sophisticated iterative MCMC algorithms, such as Population MCMC (Mohamed et al. 2012). Another aim in such a setting is Bayesian optimisation, using measures such as Expected Improvement and Expected Quantile Improvement to find points in the input space which result in minimised distance between model output and observed data (Forrester 2010; Picheny et al. 2013). Alternatively, there are several approaches drawing on the popular field of ABC (Smith and Gelfand 1992; Wilkinson 2013), such as ABC-SMC (Sequential Monte Carlo, Toni et al. (2009)), Bayes linear adjustment based ABC (Nott et al. 2014) and Rare Event ABC (Prangle et al. 2018). Whilst such approaches can effectively handle intractable likelihoods, finding very small regions in relatively high-dimensional spaces can still be challenging.

History matching is designed to efficiently cut out the uninteresting parts of the input space, thus allowing more accurate emulators to be constructed over the region of interest X⋆, where the vast majority of the mass of the posterior distribution should lie. A more detailed analysis, whether this be general Bayesian or an alternative such as ABC, can then be performed within this much smaller volume of input space, where the probabilistic specifications can now be considered more carefully. For example, Wilkinson (2014) performs a history match as a precursor to ABC. Although in that paper the likelihood function is the subject of the history match, rather than the model outputs themselves, it nevertheless highlights the benefits of carrying out a history match prior to performing a more detailed analysis. Directly emulating the likelihood may have computational advantages by reducing the emulated output dimension, which in some cases may reduce the number of required emulators substantially (although crucially probably not the number of model runs). However, the likelihood may exhibit complex behaviour which is difficult to emulate, particularly if there are several erratically behaved outputs. For discussions comparing History Matching and ABC, see, for example, Holden et al. (2018) and McKinley et al. (2018). For further information about how history matching can fit into a Bayesian paradigm, we refer the interested reader to Vernon et al. (2018), and also Wang et al. (2018), where history matching is used to help elicit reasonable prior choices from expert-elicited information.

To conclude this section, we suggest that history matching can be seen as both a useful precursor (rather than a competitor) to a more detailed analysis, and also as a form of analysis in its own right, particularly for modellers who do not wish to make the more detailed specifications required for a general Bayesian analysis. This latter situation is the motivation for this work, and we will discuss at length the insights that can be gained from having sequentially history matched the Arabidopsis model. Whilst the discussion around the results relates to the Arabidopsis model in particular, the history matching measures used within the discussion are generic and can be utilised in the context of history matching in many different fields of application.

3.1 Model Structure

3.1.1 Description and network

The model represents the hormonal crosstalk of auxin, ethylene and cytokinin of Arabidopsis root development as a set of 18 differential equations, given in Table 1, which must be solved numerically. The model takes an input vector of 45 rate parameters (k₁, k_1a, k₂,…) and produces an output vector of 18 chemical concentrations ([Auxin], [X], [PLSp],…). Note that, for simplicity, we refer to all components of the model, including hormones, proteins and mRNA, as “chemicals”. Experiments accumulated over many years have established certain relationships between some of the 18 concentrations. For example, either manipulation of the PLS gene or exogenous application of IAA (a form of auxin), cytokinin or ACC (ethylene precursor) affects model outputs [Auxin], [CK], [ET] and [PIN]. We use initial conditions for the model, given in Table 2, that are consistent with Liu et al. (2010, 2013.

Table 1:

Arabidopsis model differential equations.

d[Auxin]dt=k1a1+[X]k1+k2+k2a[ET]1+[CK]k2b[PLSp]k2c+[PLSp]+VIAA[IAA]KmIAA+[IAA]−(k3+k3a[PIN1pm]k3auxin+[Auxin])[Auxin]

d[X]dt=k16−k16a[CTR1∗]−k17[X]

d[PLSp]dt=k8[PLSm]−k9[PLSp]

d[Ra]dt=−k4[Auxin][Ra]+k5[Ra*]

d[Ra*]dt=k4[Auxin][Ra]−k5[Ra*]

d[CK]dt=k18a1+[Auxin]k18−k19[CK]+VCK[cytokinin]KmCK+[cytokinin]

d[ET]dt=k12+k12a[Auxin][CK]−k13[ET]+VACC[ACC]KmACC+[ACC]

d[PLSm]dt=k6[Ra*]1+[ET]k6a−k7[PLSm]

d[Re]dt=k11[Re*][ET]−(k10+k10a[PLSp])[Re]

d[Re*]dt=−k11[Re*][ET]+(k10+k10a[PLSp])[Re]

d[CTR1]dt=−k14[Re*][CTR1]+k15[CTR1*]

d[CTR1*]dt=k14[Re*][CTR1]−k15[CTR1*]

d[PIN1m]dt=k20ak20b+[CK][X][Auxin]k20c+[Auxin]−k1v21[PIN1m]

d[PIN1pi]dt=k22a[PIN1m]−k1v23[PIN1pi]−k1v24[PIN1pi]+k25a[PIN1pm]1+[Auxin]k25b

d[PIN1pm]dt=k1v24[PIN1pi]−k25a[PIN1pm]1+[Auxin]k25b

d[IAA]dt=0

d[cytokinin]dt=0

d[ACC]dt=0

Table 2:

The list of 18 original model outputs, along with their initial conditions. The values of 0 or 1 for IAA, cytokinin and ACC correspond to no feeding or feeding of Auxin, Cytokinin or Ethylene respectively. See Liu et al. (2010) and Liu et al. (2013) for details.

Output	Initial concentration	Output	Initial concentration
Auxin	0.1	Re∗	0.3
X	0.1	CTR1	0
PLSp	0.1	CTR1∗	0.3
Ra	0	PIN1m	0
Ra∗	1	PIN1pi	0
CK	0.1	PIN1pm	0
ET	0.1	IAA	0 or 1
PLSm	0.1	Cytokinin	0 or 1
Re	0	ACC	0 or 1

The hormonal crosstalk network of auxin, cytokinin and ethylene for Arabidopsis root development is shown in Figure 3. The auxin, cytokinin and ethylene signalling modules correspond to the model of Liu et al. (2010). The PIN functioning module is the additional interaction of the PIN proteins introduced in Liu et al. (2013). Solid arrows represent conversions whereas dotted arrows represent regulations. The v_i represent reactions in the biological system and link to the rate parameters k_i on the right hand side of the equations in Table 1. For full details of the model see Liu et al. (2013).

Figure 3:

The Arabidopsis model network for the interaction of PIN, PLS and hormonal crosstalk. The auxin, cytokinin and ethylene signalling modules correspond to the model of Liu et al. (2010). The PIN functioning module is the additional interaction of the PIN proteins introduced in Liu et al. (2013). Solid arrows represent conversions whereas dotted arrows represent regulations. The v_i represent reactions in the biological system and link to the rate parameters k_i on the right hand side of the differential equations in Table 1.

3.2 Eliciting the necessary information for history matching

To perform a history match, we need to understand how real-world observations relate to model outputs, thus aiding the specification of observed values z_i, model discrepancy terms σϵi and measurement error terms σei. History matching is a versatile technique which can deal with observations of varying quality, such as we have for the Arabidopsis model.

3.3 Input ratio ranges

The initial ranges of values for the 31 input parameters were chosen based on those in the literature (Liu et al. 2010) and further analysis of the model (Liu et al. 2013), and are shown in Table 3. Many of the input ranges were chosen to cover an order of magnitude either side of the single satisfactory input parameter setting found in Liu et al. (2010). Some parameters of particular interest were subsequently increased to allow a wider exploration of the input parameter space. This gave us a large initial input space X0 which was thought to be suitable for our purposes. The logged ratio ranges were all converted to the range [−1, 1] prior to analysis.

We now apply the technique of sequential history matching using Bayes linear emulation to the Arabidopsis model (Liu et al. 2013). Analysis of the results, after history matching to each of Datasets A, B and C, will involve consideration of the following:

• The volume reduction of the non-implausible input space (Vernon et al. 2018).

• Input plots of the non-implausible space (Vernon et al. 2018).

• The variance resolution of individual inputs and groups of inputs.

• Output plots of the non-implausible space (Vernon et al. 2018).

• The degree to which each output was informative for learning about each input.

3.4 Insights from initial simulator runs

A wave 1 set of 2000 training runs were designed using a maximin Latin hypercube design over the initial input space X0. Figure 4 shows the wave 1 output runs f_i (x) for all 32 outputs considered. The targets for the history match, as given by the intervals zi±3σci and the ranges in Table 4, are shown as vertical error bars. Black error bars represent Dataset A outputs, blue error bars represent Dataset B outputs and red error bars represent Dataset C outputs. Note that the measurements for Datasets B and C are shown for illustrative purposes, and in general may have been obtained at a later point in time to when the wave 1 model runs were performed. In this case, however, the model would still produce a value for all output components (both those with and without corresponding physical measurements).

Figure 4:

2000 wave 1 output runs f_i (x) for all 32 outputs considered. The targets for the history match, as given by the intervals zi±3σci and the ranges in Table 4, are shown as vertical error bars. Black error bars represent Dataset A outputs, blue error bars represent Dataset B outputs and red error bars represent Dataset C outputs. The horizontal black line at zero represents zero trend.

Figure 4 gives substantial insight into the general behaviour of the model over the initial input space X0, for example informing us about model outputs that can take extreme values, for example, f_c_Auxin, f_c_ET and f_af_c_PLSm. More importantly, the runs also inform us as to the class of possible observed data sets that the model could have matched, and hence gives insight into the model’s flexibility. There exist outputs with constrained ranges. In particular, many outputs seem to be constrained to being either positive or negative, for example, the logged trend for pls_CK must be positive and that of PLSox_CK must be negative. If such constrained outputs, which are consequences of the biological structure of the model, are found to be consistent with observations, this provides (partial) evidence for the model’s validity. Conversely, we may be concerned about an overly flexible model that was capable of reproducing any combination of positive or negative observed data values for outputs in this dataset. Specifically, we should doubt claims that such a model has been validated by comparison to this data, as it would have inevitably matched any possible data values and hence arguably may not contain much inherent biological structure at all.

There are some outputs for which the majority of the wave 1 runs already go through the corresponding error bars, for example PLSox_Auxin and PLSox_ET. This is an indication that these outputs did not help much to constrain the input space. Despite this, none of the wave 1 runs pass through all of the target intervals of the outputs in Dataset A simultaneously, thus already suggesting that the volume of the final non-implausible space would be small or indeed zero.

3.5 History matching the model

We outline the general decisions required to perform the history match. Several packages are available that perform standard Gaussian Process emulation, possibly with Automatic Relevance Determination (Neal 1997; Williams and Rasmussen 1996), for example the BACCO (Hankin 2005) and GPfit (MacDonald et al. 2015) packages in R (R Core Team) or GPy (since 2012) for Python, which may be used as an alternative to the emulators we describe here. Emulators accurate enough to reduce the size of the non-implausible space at each wave to some degree are sufficient for our purposes. When constructing emulators, we decided to put more detail into the mean function, but incorporate more complicated structures for the residual process at each wave, thus sequentially increasing the complexity of the emulators at each wave. We provide a summary of the choices made in the history match at each wave in Table 5, including the dataset history matched to (column 2), the number of design runs (column 3), the implausibility cut-off thresholds (columns 4–6) and the emulation strategy (column 7), each of which is discussed in more detail below.

The amount of space that was cut out after each wave is shown in Table 6. We let V(Xi) represent the volume of the non-implausible space after wave i, as judged by the emulators, and V(XG) represent the volume of the space with acceptable matches to the observed data in Dataset G, as judged using actual model runs (hence without emulator error). Then columns 2 and 3 give the proportion of the previous wave and initial non-implausible spaces respectively still classed as non-implausible, and columns 5 and 6 give the proportion of the wave i and initial non-implausible spaces giving rise to actual acceptable matches to the data in Dataset G. The proportion of space cut out at each wave is influential for deciding the number of waves and emulator technique at each wave. In addition, Table 6 presents the radical space reduction obtained by performing the history match. A proportion of 6.1 × 10⁻⁷ of the original space was still considered non-implausible after history matching to Dataset A. A proportion of only 8.5 × 10⁻¹⁰ of the original space was still considered non-implausible after history matching to Datasets A and B, thus the 5 trends in Dataset B, for exogenous application of ACC, facilitated an additional reduction of 3 orders of magnitude. After all experimental observations had been matched to, the non-implausible space had been reduced to a proportion of 7.2 × 10⁻¹² of the original space, thus the 5 trends in Dataset C, for measurement of POLARIS gene expression, refocused the set by another 2 orders of magnitude. Such small proportions of the original space being classed as non-implausible means that acceptable runs within these spaces would likely be missed by more ad-hoc parameter searching methods of analysis.

Linear model emulators with uncorrelated residual processes were used in the initial waves since they are very cheap to evaluate, substantially more so even than emulators involving a correlated residual process, which may only be slightly more accurate (Andrianakis et al. 2017a). For these emulators, we estimated the value of σui2 to be the estimated variance parameter from the linear model fit. As the amount of space being classed as implausible at each wave started to drop, we introduced emulators with a Gaussian correlation residual process. There are various methods in the literature for assessing correlation length parameters, as explained in Section 2.1. Some of the methods in the literature for picking the correlation lengths θ and variance parameter σui2 tend to be computationally intensive and their result highly sensitive to the sample of simulator runs (Andrianakis and Challenor 2009, 2011). The choice was therefore made, at wave 5, to use a single correlation length parameter value of θ = 2 for all input-output combinations, and to fit σui2 using the corresponding linear model fit, these choices being checked using emulator diagnostics. The motivation for this choice of correlation length parameter was made by appealing to the heuristic argument made in Vernon et al. (2010a) that the regression residuals may be derived from a polynomial of order one higher than the fitted polynomials, the alteration in the chosen value taking into account the higher dimensionality of the input space.

At wave 6 the complexity of the residual process was increased still further by splitting the active inputs xAi for each output emulator into five groups based on similar strength of effect, and using maximum likelihood to fit the same correlation length to all inputs in each group, along with the variance parameter σui2. This extension to the literature of fitting several different correlation length parameter values strikes a balance between the stability of the maximum likelihood process (which can become very challenging were we to include 31 separate correlation lengths for each of the 31 input components) and the overall complexity of the residual process. It should be noted that maximum likelihood makes use of probabilistic distributions to make an assessment of the correlation length parameters. Whilst it can be argued that such an approach lies outside of the Bayes linear paradigm in which we presented the construction of our emulators, it provides a useful tool for calculating adequate emulators which satisfy diagnostics. At wave 8 we introduced the Dataset B outputs by first using linear model emulators for the new outputs only, and then using emulators with residual correlation processes for all outputs. In waves 12 and 13 we incorporated emulators with residual processes for the Dataset C outputs.

The number of design points per wave was largely kept constant at 2000. 2000 was deemed a suitable number of runs per wave as it meant that the matrix calculations involved in the emulator were reasonable, whilst allowing adequate coverage of the non-implausible input space with simulator runs. At waves 11 and 13, 3500 design runs were used to build more accurate emulators.

In terms of design, a maximin Latin hypercube (Currin et al. 1991; McKay et al. 1979) was deemed sufficient for our needs as we required a simple and efficient space-filling design. The speed of the simulator meant that more structured and tactical designs were unnecessary for our requirements. At wave 1 we constructed a Latin hypercube of size 2000 to build emulators for each of the outputs. At waves 2–7 we first built a large maximin Latin hypercube design containing a large number of points over the smallest hyper-rectangle enclosing the non-implausible set. We then used all previous wave emulators and implausibility measures to evaluate the implausibility of all the proposed points in the design (Vernon et al. 2010a). Any points that did not satisfy the implausibility cut-offs were discarded from further analysis. If a single Latin hypercube was not sufficient to generate enough design points, multiple Latin hypercubes were taken in turn and the remaining points in each were taken to be the design. From wave 8 onwards an alternative sampling scheme was necessary to generate approximately uniform points from the non-implausible sets, since generating points using Latin hypercubes became infeasible due to the size of the non-implausible space. There are several alternative ways presented in the literature to approximately sample uniformly distributed points over the non-implausible space (Andrianakis et al. 2015, 2017a; Gong and DiazDelaO 2017; Williamson and Vernon 2013). We used a simple Metropolis-Hastings MCMC algorithm (Brooks et al. 2011), which provided adequate coverage of the non-implausible space.

At each wave we performed diagnostics on the mean function linear model, the emulator and the implausibility criteria using 200 points in the non-implausible space. The diagnostic test for implausibility which we used was the one described in Vernon et al. (2018). It compares the data implausibility cut-off criteria IMdata(x), that is the implausibility evaluated at a known diagnostic run, against the chosen implausibility cut-off criteria for the emulator outputs.

Many waves were necessary to complete this history matching procedure due to the complex structure of the Arabidopsis model. We see from the first few waves that linear model emulators are sufficient for learning a great deal about the input parameter space, but that including full emulators with correlated residuals at later waves can be useful. In addition, emulators have greatly increased the efficiency of our analysis over simply running the model. To make a comparison, the model itself takes 30 s to run on a standard laptop computer (not very slow, but slow enough to cause problems) in comparison to the tens of millions of runs per second for the linear model emulators and tens of thousands of runs per second for the more complex later-wave emulators. At the end of the procedure we obtained many runs satisfying each of the datasets A, B and C. We now go on to describe the results of the parameter search using various graphical representation techniques and discuss their biological implications.

4.1 Gaining insight into specific learning objectives from history matching results