Prior Predictive Checks for Mechanistic Models

Why Check Predictions Before Fitting?

A prior can look reasonable parameter by parameter and still imply unreasonable model behaviour. For example, a prior may assign plausible ranges to a production rate and an inhibition parameter, yet their combination may imply negative production at some doses. A prior predictive check asks what the model predicts before the observed data are used. It can reveal implausible concentrations, rates, and time courses while there is still an opportunity to revise the model or its prior assumptions.

For a mechanistic QSP model, a prediction usually begins with a parameter draw and a simulation of the model equations. The resulting trajectory is not itself a parameter draw: it is the consequence of the sampled parameters, initial conditions, inputs, and model structure. A useful check therefore examines both the sampled parameter values and their joint implications for the system.

Definition: Let \(\theta\) denote the model parameters, with prior density \(p(\theta)\). Let \(x(t;\theta)\) be the model state at time \(t\), obtained by solving the specified mechanistic model for parameter value \(\theta\). A prior predictive check simulates \(\theta\) from \(p(\theta)\), computes the associated state trajectories and any desired predictions, and assesses whether those predictions are scientifically plausible before conditioning on observed data.

For a deterministic ODE model, the trajectory is determined once the parameter values, initial conditions, and input schedule are fixed. If predictions also include a stochastic observation process, one can simulate a replicated observation conditional on each parameter draw. This tutorial focuses on the mechanistic trajectories and their concentrations and rates; observation-error assumptions should be kept distinct from the ODE’s biological implications.

Prior Predictive Distribution and Simulation

The prior predictive distribution combines uncertainty in the parameters with the model’s predictions conditional on those parameters. If \(Y\) is a predicted quantity, such as a concentration at a specified time, and \(p(Y\mid\theta)\) is its conditional distribution under the model, then

p(Y)=\int p(Y\mid\theta)p(\theta)\,d\theta.

For a deterministic prediction \(Y=f(\theta)\), \(p(Y\mid\theta)\) places all its probability at \(f(\theta)\). The integral then means that predictions vary because the parameter values vary. In practice, the integral is usually approximated by simulation: draw \(\theta^{(m)}\sim p(\theta)\), solve the model for each draw, and summarize or plot the resulting trajectories.

Specify the complete prior and model inputs.
Include the joint prior, units, initial conditions, dosing or perturbation schedule, and simulation time window. Preserve dependencies among parameters rather than sampling them independently if the prior specifies dependence.

Draw parameter sets.
For each draw, sample from the stated prior on the model’s parameter scale. Apply any parameter transformations consistently with the prior, as established in the earlier material on induced priors under parameter transformations.

Simulate and check.
Solve the mechanistic equations over the relevant times and inputs. Inspect parameter values, rates, states, and whole trajectories against scientific constraints and plausible ranges.

Trace problems to their source.
Determine whether an implausible result follows from the prior, the model structure, its inputs, or a numerical failure. Revise only after identifying which assumption generated the result.

The calculation relies on a usable prior distribution and a well-defined model solution for the parameter draws being examined. The initial conditions and inputs must be specified, and the solver must be sufficiently accurate for the features being assessed. If a prior places mass on parameter values for which the model is undefined or the ODE has no acceptable solution, that is itself an important warning, not a reason to silently discard those draws.

Worked Calculations With a Concentration Model

Consider a deliberately simple, hypothetical biomarker concentration \(C(t)\), measured in nM, under a constant drug concentration \(D\). Its dynamics are

\frac{dC}{dt} = k_{\mathrm{in}}\left(1-E_{\max}\frac{D}{IC_{50}+D}\right)-k_{\mathrm{out}}C.

Here \(k_{\mathrm{in}}\) has units nM/day, \(k_{\mathrm{out}}\) has units day\(^{-1}\), \(D\) and \(IC_{50}\) have the same concentration units, and \(E_{\max}\) is dimensionless. The bracketed term scales the baseline production rate. This model is used only to illustrate the checks; its assumptions include constant \(D\), a constant production term over the simulation, and first-order removal.

Worked Example: Checking One Draw and Its Trajectory

Take the hypothetical draw \(k_{\mathrm{in}}=12\) nM/day, \(k_{\mathrm{out}}=0.4\) day\(^{-1}\), \(E_{\max}=1.6\), \(IC_{50}=5\) nM, and \(D=5\) nM. Let \(C(0)=30\) nM. First check the net production rate.

k_{\mathrm{in}}\left(1-E_{\max}\frac{D}{IC_{50}+D}\right) = 12\left(1-1.6\frac{5}{5+5}\right) = 12(1-0.8) = 2.4\ \text{nM/day}.

The net production is positive: \(1-0.8=0.2\), so it is \(0.2\times12=2.4\) nM/day. At equilibrium, production balances removal, giving \(C_{\mathrm{ss}}=2.4/0.4=6\) nM. The solution for a constant production rate is

C(t)=C_{\mathrm{ss}}+\bigl(C(0)-C_{\mathrm{ss}}\bigr)e^{-k_{\mathrm{out}}t} =6+24e^{-0.4t}\ \text{nM}.

At day 3, \(e^{-1.2}\approx0.3012\), so \(C(3)=6+24(0.3012)\approx13.23\) nM. This trajectory declines from 30 nM toward 6 nM and remains positive. The example checks not just the parameter values but their combined consequence for both the rate and the concentration.

Worked Example: Estimating the Prior Probability of Negative Production

Suppose the hypothetical prior for \(E_{\max}\) is Uniform(0, 2.4), with the other quantities and dose as above. At \(D=IC_{50}\), the production multiplier is \(1-E_{\max}/2\). Production is negative when \(E_{\max}>2\). The prior probability of that event is

P(E_{\max}>2)=\frac{2.4-2}{2.4-0}=\frac{0.4}{2.4}=\frac{1}{6}\approx0.1667.

Thus the prior assigns about 16.67% probability to a negative production rate at this dose. In a hypothetical simulation of 1,000 independent prior draws, suppose 171 produce a negative rate. The simulated proportion is \(171/1000=0.1710\), reasonably close to \(0.1667\). Its approximate Monte Carlo standard error is

\sqrt{\frac{(1/6)(5/6)}{1000}}\approx0.01179.

This is a check on the prior-model combination, not a hypothesis test about whether the prior is true. Whether negative production is scientifically impossible depends on the model interpretation. If the term represents literal production, it is a physical warning; if the equation is an empirical net-effect parameterization, its meaning requires separate justification.

For the particular draw \(E_{\max}=2.4\), the rate is \(12(1-2.4/2)=-2.4\) nM/day. With \(k_{\mathrm{out}}=0.4\) day\(^{-1}\), its equilibrium is \(-2.4/0.4=-6\) nM. Starting at 30 nM, the solution is \(C(t)=-6+36e^{-0.4t}\). It reaches zero when \(e^{-0.4t}=1/6\), or \(t=\log(6)/0.4\approx4.479\) days. The negative concentration thereafter is incompatible with a literal concentration state, exposing a consequential implication of this prior draw.

Worked Example: Screening Whole Trajectories

Use the same ODE, \(k_{\mathrm{in}}=12\) nM/day, \(k_{\mathrm{out}}=0.4\) day\(^{-1}\), \(D=IC_{50}=5\) nM, and \(C(0)=30\) nM. For four hypothetical prior draws \(E_{\max}=0.4, 1.0, 1.6,\) and \(2.4\), compute \(C(3)\) and \(C(6)\). Since \(D=IC_{50}\), the equilibrium for each draw is \(C_{\mathrm{ss}}=30(1-E_{\max}/2)\), and

C(t)=C_{\mathrm{ss}}+(30-C_{\mathrm{ss}})e^{-0.4t}.

\(E_{\max}\)	\(C_{\mathrm{ss}}\) (nM)	\(C(3)\) (nM)	\(C(6)\) (nM)
0.4	24	\(24+6e^{-1.2}=25.81\)	\(24+6e^{-2.4}=24.54\)
1.0	15	\(15+15e^{-1.2}=19.52\)	\(15+15e^{-2.4}=16.36\)
1.6	6	\(6+24e^{-1.2}=13.23\)	\(6+24e^{-2.4}=8.177\)
2.4	-6	\(-6+36e^{-1.2}=4.843\)	\(-6+36e^{-2.4}=-2.734\)

For example, \(e^{-2.4}\approx0.09072\), giving \(-6+36(0.09072)=-2.734\) nM for the last draw at day 6. This four-draw illustration is too small to estimate useful prior predictive quantiles, but it shows why inspecting complete trajectories matters: the problematic draw is still positive at day 3 and becomes negative later. In a real check, evaluate enough time points to detect transient or delayed failures, and compare trajectory ranges with scientifically justified bounds.

Interpreting Diagnostics and Failure Modes

A prior predictive check is not a search for trajectories that resemble an observed data set; the data have not yet been used. Instead, ask whether the prior-model combination generates outputs that make sense under the intended biology and experimental context. Useful summaries include the fraction of simulations violating a hard constraint, the ranges of concentrations and rates, and plots of many trajectories over time. When constraints are soft rather than absolute, describe how often predictions enter regions judged implausible and why those regions matter.

Check units and scales. A rate in day\(^{-1}\) must multiply time in days in the ODE. Convert dose and concentration units consistently before interpreting saturable terms.
Inspect joint effects. A parameter value can be plausible on its own but problematic in combination with another parameter, an initial condition, or a dose.
Check more than endpoints. A final concentration can appear acceptable even when the trajectory has an extreme peak, crosses zero, or changes implausibly quickly.
Separate scientific failure from numerical failure. Repeat suspicious solves with tighter tolerances or a suitable solver. If the result changes materially with numerical settings, diagnose the computation before interpreting the trajectory as a model implication.
Do not silently remove failed draws. Record how many draws fail and why. Frequent failures may indicate that the prior includes invalid regions or that the model is not well-defined across its intended parameter support.

The diagnostic thresholds require scientific justification: a concentration bound should reflect the modeled state and context, not be selected merely because it makes a plot look tidy. When a result is implausible, possible responses include revising the prior support or dependence, changing the mechanistic formulation, or correcting a unit or input specification. A prior predictive check identifies the issue; it does not decide which response is biologically appropriate.

How This Check Differs From Alternatives

Reviewing marginal prior ranges is simpler, but it misses interactions and the nonlinear transformations created by the model. Prior predictive simulation carries each joint draw through the equations, so it can reveal implausible combinations and time courses directly. Conversely, a plausible parameter table does not guarantee plausible predictions.

Checks made after fitting answer a different question: they assess predictions conditional on information from the observed data and the fitted model. Such checks cannot replace a prior predictive check, because data-informed predictions may conceal a prior that was excessively broad or scientifically incoherent before fitting. Nor does a prior predictive check establish that the model is correct or that its parameters will be identifiable. Its scope is narrower: it assesses the prior’s implications under the specified model and inputs.

A practical limitation is computational cost. Large ODE systems, stiff dynamics, or many prior draws may make simulation slow. A smaller pilot can locate obvious problems, but it may miss rare failures; increase the number of draws when estimating tail probabilities, and report simulation counts and uncertainty for estimated proportions. Prior predictive summaries also depend on the chosen scenario. A prior might yield plausible baseline trajectories but implausible responses under a treatment dose, so check scientifically relevant inputs and time horizons.

Common Mistake

Common mistake: “Every parameter prior has a biologically plausible range, so the model’s prior predictions must also be plausible.” This does not follow: nonlinear equations combine parameters, initial conditions, and inputs. For example, \(E_{\max}=2.4\) is within Uniform(0, 2.4), but at \(D=IC_{50}\) it makes the production multiplier \(1-2.4/2=-0.2\). The correction is to simulate the joint prior through the model and inspect rates and full trajectories, not to rely on marginal parameter ranges alone.

Check Your Understanding

Work through the questions, then compare with the solutions below.

In one sentence, distinguish a prior parameter draw from a prior predictive trajectory.
For \(E_{\max}=1.0\), \(D=IC_{50}\), and \(k_{\mathrm{in}}=12\) nM/day, calculate the net production rate.
Under \(E_{\max}\sim\mathrm{Uniform}(0,2.4)\), what proportion of prior draws imply negative production at \(D=IC_{50}\)?
Why can checking only a day-3 concentration miss the failure in the \(E_{\max}=2.4\) trajectory?
Name one way to distinguish a numerical solver problem from an implausible model implication.

Solutions. (1) A parameter draw is a sample of model inputs from the prior; its trajectory is the model solution generated from that draw and the specified initial conditions and inputs. (2) The multiplier is \(1-1/2=0.5\), so production is \(12(0.5)=6\) nM/day. (3) Negative production occurs for \(E_{\max}>2\), so the proportion is \((2.4-2)/2.4=1/6\), or about 16.67%. (4) That draw gives \(C(3)\approx4.843\) nM, still positive, but \(C(6)\approx-2.734\) nM; the crossing occurs after day 3. (5) Repeat the solve with tighter tolerances or another appropriate solver and compare the trajectories; a material change suggests numerical error needs investigation before biological interpretation.