Prediction uncertainty assessment of a systems biology model requires a sample of the full probability distribution of its parameters

Model class

We consider models that are formulated in terms of differential equations and have the following form: ${\displaystyle \dot{x}\left(t\right)=f\left(x\left(t\right),\theta ,u\left(t\right)\right)}$ ${\displaystyle y\left(t\right)=g\left(x\left(t\right),\theta \right)}$ (1)${\displaystyle x\left(0\right)=h\left(\theta \right).}$ The dynamics of the internal state vector x is a function f of parameter vector θ∈ℝ^p, with initial condition x(0), and external input vector u(t). The output y∈ℝ^m is a positive function g of the internal state and θ. To obtain the output time series y(t), we have to numerically integrate $\dot{x}\left(t\right)=f\left(x\left(t\right),\theta ,u\left(t\right)\right)$, which is the time consuming part of the solution.

Estimation of the posterior distribution of the parameters

We adopt the Bayesian framework, in which a prior distribution of the parameters is combined with data to form the posterior distribution π(θ) of the parameters, according to Bayes’ theorem (2)${\displaystyle \pi \left(\theta \right)=p\left(\theta |y_d\right)=\frac{p\left(\theta \right)p\left(y_d|\theta \right)}{p\left(y_d\right)}.}$ Here p(θ) denotes the density of the prior distribution of the parameters based on prior knowledge and p(y_d|θ) the likelihood of the measured data y_d given θ. The denominator, p(y_d), is the marginal likelihood of the data y_d and acts as a θ-independent normalization constant which is not relevant in MCMC algorithms. We assume that the output y(t_i), i = 1, …, n contains uncorrelated Gaussian noise, with variance ${\sigma }_i^2$ per time point i, so the likelihood of the measured data y_d, given θ, is (3)${\displaystyle p\left(y_d|\theta \right)=\prod _{i=1}^n\frac{1}{\sqrt{2\pi \underset{i}{\overset{2}{\sigma }}}}exp\left(-\frac{{\left(y_i\left(\theta \right)-y_{d,i}\right)}^2}{2\underset{i}{\overset{2}{\sigma }}}\right),}$ with y_i(θ) the model output at time t_i given θ. Inserting (3) into (2) and taking the logarithm gives the log-posterior (4)${\displaystyle log\left(\pi \left(\theta \right)\right)=c-\frac{1}{2}{\chi }^2\left(\theta \right)+log\left(p\left(\theta \right)\right),}$ with c a constant independent of θ, and χ² a measure of the fitting error: (5)${\displaystyle {\chi }^2\left(\theta \right)=\sum _{i=1}^n\frac{{\left(y_i\left(\theta \right)-y_{d,i}\right)}^2}{\underset{i}{\overset{2}{\sigma }}}.}$ The penalized maximum likelihood parameter θ^PML maximizes (4). Draws from the posterior π(θ) are obtained by an MCMC algorithm, which, starting from a user-defined initial value of θ, generates a stochastic walk through the parameter space. In iteration k of the walk, a new candidate solution θ is proposed based on the current solution θ_k. We use symmetric proposal distributions centered at θ_k. In this case the proposed θ is accepted or rejected using the Metropolis acceptance probability min(1, r), where $r=\frac{\pi \left(\theta \right)}{\pi \left({\theta }_k\right)}$. So, the more probable a parameter vector is with respect to the data, the more probable it is to be accepted. In this paper the proposals for the walk are generated by the DE-MCz algorithm, which is an adaptive MCMC algorithm that uses multiple chains in parallel and exploits information from the past to generate proposals (ter Braak & Vrugt, 2008). As consecutive draws are dependent, it may be practical to reduce the dependence by thinning, that is, by storing every Kth draw (with K > 1). After a number of iterations (the burn-in period) the chain of points {θ_k} will be stationary distributed with a local density that represents the posterior π(θ). The burn-in iterations are discarded.

A quantifier of uncertainty in predicted systems dynamics

Since we want to compare the uncertainty in predictions of time courses over different time intervals and for different models, we need a quantifier that is independent from model specific issues, such as the number, and the dimensions of the variables. Also, biologists are generally interested in percentage difference rather than absolute difference, so the quantifier should be independent of the typical order of magnitude. To that end we introduce a measure for prediction uncertainty that satisfies these requirements, but we do not claim that this quantifier is unique. However, our present choice has the advantage that it allows a nice interpretation as will be explained below. We first define a measure for the prediction uncertainty of one component of the output and then average the outcomes over the components.

For a predicted output y_p(t, θ) on time interval [0, T], the quantifier Q represents the relative deviation of y_p(t, θ) from the penalized maximum likelihood prediction y_p(t, θ^PML), integrated over time: (6)${\displaystyle Q_i\left(\theta \right)=\frac{1}{T}{\int }_0^T{\left({log}_b\left(\frac{y_{p,i}\left(t,\theta \right)}{y_{p,i}\left(t,{\theta }^{PML}\right)}\right)\right)}^{2}dt,}$ with y_p,i the ith component of y_p, with i = 1, …, m, and $Q=\frac{1}{m}\sum _i{Q}_i$. Q integrates differences over time on a logarithmic scale with base b, and is therefore symmetric with respect to relative over- and underestimations (as opposed to, e.g., a sum of squared errors). The underlying assumption is that y is always positive, which holds true for all models describing concentration dynamics. When only a few time points are of biological interest, the integrand in (6) may be approximated by a summation. The base b characterizes the order of magnitude of the discrepancies that Q is sensitive to. For example, if $\frac{1}{b}<\frac{y_{p,i}}{y_{p,i}^{ML}}<b$ holds for all time points, this will result in Q < 1, while a few points outside this range will quickly result in Q > 1. In this paper we use b = 2, so only differences of a factor two or higher can result in Q > 1. The choice of b represents the maximum magnitude of deviations that a prediction is allowed to have, so it should in general be selected based on biological grounds.

If {θ} is the collection of sampled θ’s reliably representing density π(θ), then the density π(Q(θ)) is reliably represented by Q(θ) with θ∈{θ} (Robert, Marin & Rousseau, 2011). We denote the level α prediction uncertainty with Q_α, the α quantile of the distribution of Q. Q_α is therefore the deviation of a prediction relative to the penalized maximum likelihood prediction, at confidence level α (α-level deviation). We use α = 0.95 throughout.

Algorithm to estimate prediction uncertainty

The algorithm for estimating prediction uncertainty naturally falls apart in two sub-algorithms. Part I deals with estimating parameters by exploiting the prior knowledge, the model, and the data. This first part yields the sample of parameter values representing their posterior distribution. Part II is the focus of this paper, and performs the full computational uncertainty analysis by taking as input the sample of parameter values from the first part and by calculating the prediction for each member of this sample. Part I is computationally more intensive than Part II. Because prospective users of the model need to carry out only Part II, it is essential that the sample of parameter values obtained in Part I is stored and made available. In full:

• Part I: Estimation of the posterior parameter distribution

  1. To calculate the posterior π(θ) in Eq. (2), we use Eq. (3) for the likelihood, together with a log-uniform prior for the parameters, that is, a uniform prior for the logarithm of the parameters.

  2. A collection of solutions {θ} is generated by MCMC sampling; in this paper we use DE-MCz (ter Braak & Vrugt, 2008). This collection constitutes the sample approximating π(θ).

• Part II: Computational uncertainty analysis

  1. Q(θ) is computed for each θ∈{θ}. The α level prediction uncertainty Q_α is approximated by taking the largest Q after discarding the 100α% largest Q(θ) values.

  2. For visualization of the uncertainty in the predicted systems dynamics, for each θ∈{θ} the output y(t) is calculated. For each time point the $\frac{1}{2}\left(1-\alpha \right)100\%$ largest and smallest predicted values of y(t) are discarded, whereafter the minimum and maximum values are plotted, creating an envelope of predicted dynamics.

A structural difference between Q and the visualization is that Q integrates the deviation over the total function y(t), while the visualization displays prediction uncertainty per time point and per component of y(t). An alternative to definition (6) is to define Q as the deviation of log_b(y_p(t, θ)) from the average log-prediction, thus replacing log_b(y_p(θ^PML)) with $\overline{{log}_b\left(y_{p,i}\left(\theta \right)\right)}$.

Illustrative example

To demonstrate what kind of effects can be expected when studying prediction uncertainty of models with parameter uncertainty, we use a model example that represents a highly simplified case of self-regulation of gene transcription: (7)${\displaystyle \dot{y}\left(t\right)=\frac{10y{\left(t\right)}^2}{{\theta }_1+y{\left(t\right)}^2}-{\theta }_{2}y\left(t\right)+y\left(t\right)u\left(t\right).}$ Here, y(t) represents protein concentration and $\dot{y}\left(t\right)$ its rate of change. Changes in y(t) stem from self-regulation (the first term in the equation, which is a Hill function Alon, 2006), decay (the second term), and an input function (the third term) given by u(t) = sin(t). The input represents some external signal, e.g., a triggering mechanism (van Mourik et al., 2010). The unknown parameters involved in self-regulation and decay are θ₁ and θ₂, respectively. They are to be estimated from time series data, for which we used simulated data at ten time points (Fig. 2A). We assumed independent Gaussian noise on each measurement, with a mean of 10% of the y value, and a log-uniform prior for the parameter distribution. We estimated θ₁ and θ₂ by generating 1000 draws from their distribution using MCMC sampling. This yielded a sample of 1000 draws. From this sample we derived a 95% credible region of θ₁ and θ₂ (Fig. 2B) and also point-wise 95% credible intervals (envelopes) for three sets of predictions (Fig. 2C). The prediction sets differ in starting value y(0) or the input function u(t). The envelopes for prediction (Fig. 2C) are obtained by inserting each of the 1000 draws into the (perturbed) system, calculating the prediction on a fine grid of time points by evaluating the (modified) differential equation (7) and summarizing the 2.5 and 97.5 percentiles of the predictions at each time point. Details are given in the Methods section.

The data give much more information about θ₂ than about θ₁ (Fig. 2B, note the different scales on the axes), that is, the dynamics of the time series is sensitive to changes in θ₂, but not sensitive to changes in θ₁; the model thus shows sloppiness. The reason is not hard to understand: since we have used a high initial value, most of the y(t) curve has values much bigger than θ₁ = 10, so that the parameter θ₁ has little influence on the dynamics (see Eq. (7)). In more complex models, credible regions may have far more complicated shapes than in this illustrative example.

The prediction sets Prediction 1 and Prediction 2 (Fig. 2C) are obtained by starting the system for each draw of the sample at a low initial concentration (y(0) = 1) and at a very high concentration (y(0) = 10⁴), respectively. The uncertainty in prediction increases in time for Prediction 1 and is much higher than that for Prediction 2. The predicted dynamics is thus highly sensitive to parameter uncertainty for Prediction 1, but fairly insensitive for Prediction 2. Based on Predictions 1 and 2 one might guess that the uncertainty is related to the magnitude of y(t), but this guess is wrong as the third prediction set shows. In this set, we predict the time course starting at the same low value as in Prediction 1 but now the input amplitude is increased 7-fold (u(t) = 7sin(t)). This drives the concentration dynamics into an oscillating motion which can be predicted rather precisely (Fig. 2C), since now the input u(t) dominates the dynamics.

We also estimated prediction uncertainty via a local sensitivity analysis, namely linearized covariance analysis (LCA, see the Supplemental Information for details). LCA extrapolates parameter- and prediction uncertainty from a second- and first order Taylor expansion, respectively, and does not take into account the higher order terms that in this case constrain the credible region and the predictions, but note that linearization can also lead to predictions that appear more precise than they really are. Consequently, LCA dramatically overestimates the uncertainties of Predictions 1 and 3 (Fig. 2D).

This simple example illustrates that prediction uncertainty depends on the type of prediction and is hard to foresee intuitively. We cannot evaluate the prediction uncertainty on the basis of some simple guidelines after only inspecting the credible region of the parameters. Because of the nonlinearities in the system, the credible region is untruthful for evaluating predictions and their uncertainty (Savage, 2012). These results suggest that a full computational uncertainty analysis, as performed here (Fig. 2C), has always to be done in order to estimate the consequences of parameter uncertainty for prediction uncertainty. The key element in this analysis is the sample representing the probability distribution of the parameters.

It is convenient to quantify prediction uncertainty, that is, the variability among predicted time courses. We do so on the basis of the quantity Q defined in Eq. (7) in the Methods section, which expresses the deviation between two time courses, one being calculated for an arbitrary set of parameter values and the other using the penalized maximum likelihood estimator of the parameters. In our simulation studies, the latter is replaced by the true parameter values. The quantity Q is calculated for each draw of the parameters. The higher the Q value of a prediction, the higher its prediction uncertainty is. Each prediction set yields 1000 values of Q, which are summarized by their 95% percentile, indicated by Q_0.95. Q_0.95 values smaller than 1 indicate tight predictions. For example, the Q_0.95 values for Predictions 1, 2 and 3 in Fig. 2C are 18, 0.02 and 1.5, respectively. Even for this simple model the uncertainty among different predictions thus varies by a factor 900 which is an order of 2.9 on logarithmic scale.

Scanning a variety of systems biology models

Next, we study prediction uncertainty for six models from the BioModels database (Li et al., 2010). We selected these models to represent some cross section of systems biology models; The systems they represent greatly differ in numbers of variables, parameters, and the types of equations. The parameters in these models represent process rates. The models include circadian rhythm, metabolism, and signalling (Laub & Loomis, 1998; Levchenko, Bruck & Sternberg, 2000; Tyson, 1991; Goldbeter, 1991; Vilar et al., 2002; Leloup & Goldbeter, 1999), see the Supplemental Information for details. For each model, a data set is created in the same manner as described for the simple model above, with parameter values and initial conditions as provided in the database. Each data set was analyzed using an MCMC sampling algorithm yielding a sample of 1000 draws representing the distribution of the parameters. This number turned out to be large enough for our prediction purposes as checked in the Supplemental Information. Subsequently, different sets of predictions were made. For each prediction set, the system was perturbed and the time course of the perturbed system was predicted on the basis of each of the 1000 draws of parameter values, whereafter prediction uncertainty was quantified on the basis of Q_0.95. Perturbations consisted of artificially increasing and decreasing each parameter in turn by a factor of 100 (different factors yielded similar results, see Supplemental Information). For a model with p parameters, we so obtain 2p different perturbations, prediction sets, and Q_0.95 values.

For the six models the prediction uncertainty is displayed in terms of the Q_0.95 quantifiers (Fig. 3). For the collective fits we used a data set with 10 data points per variable, with an assumed Gaussian noise of σ_i = 0.1y_d,i per data point i. (Fig. 3A). The sizes and locations of the Q_0.95 distribution ranges considerably differ per model. For example, the Tyson model gives mostly precise predictions, and the Laub model gives mostly uncertain predictions. Next, the accuracy of the data was increased artificially by decreasing the noise with a factor 10, and observed the effect on the Q_0.95 distribution ranges (Fig. 3B). With more accurate data, the Q_0.95 values are reduced, but not in a very systematic way. For the Vilar model, for example, all predictions have become precise, except for one. More importantly, for all six models the Q_0.95 values scatter enormously, in ranges that span 2 to 4 orders of magnitude. There is no clear connection between the number of parameters and the maximum Q_0.95 values or the size of the ranges. By consequence, prediction uncertainty is not a characteristic of a calibrated model, but must be evaluated on a case-by-case basis.