From Parameter Uncertainty to Mechanistic Quantities
A calibrated QSP model usually has many uncertain parameters, but the scientific question often concerns a quantity calculated from the model: how strongly a pathway is inhibited, how long a response stays above a target, or when a biomarker recovers. These quantities are functions of parameters and, when relevant, of a specified dosing regimen, population, and other scenario choices. Their posterior distributions summarize uncertainty about the mechanistic quantity after accounting for the observed data and the model.
The key calculation is to evaluate the quantity for each posterior draw, not to calculate it only at a single best-fitting parameter vector. Because the model may be nonlinear, a summary calculated at the posterior mean parameters need not equal the posterior mean of the derived quantity.
For example, if \(Q\) is pathway inhibition at a specified time, each posterior parameter draw produces a model-predicted treated response and a corresponding inhibition value. If \(Q\) is time above a target, each draw produces a trajectory from which the duration is calculated. The scenario must be stated: changing the dose, target, baseline definition, or patient profile can change the derived quantity even when the posterior parameter draws are unchanged.
Posterior Draws, Probabilities, and Quantiles
Suppose \(M\) usable posterior draws are available, written \(\theta^{(1)},\ldots,\theta^{(M)}\). For each draw, calculate \(Q^{(m)}=g(\theta^{(m)};c)\). The resulting values are draws from the posterior distribution of \(Q\), subject to the quality of the posterior simulation and the accuracy of the calculations.
Here, \(A\) is a set of values, \(I\{\cdot\}\) equals 1 when its condition is true and 0 otherwise, and \(\pi(\theta\mid y)\) denotes the posterior distribution of the parameters. The integral says that the probability of a derived quantity lying in a set is the posterior probability of parameter values that produce quantities in that set. With posterior draws, estimate this probability by the fraction of transformed draws in \(A\).
A quantile is another useful summary. Define the posterior \(q\)-quantile \(Q_q\) as the smallest value \(x\) for which \(P(Q\leq x\mid y,c)\geq q\). In a finite sample of draws, use the empirical inverse-CDF convention: sort the values as \(Q_{(1)}\leq\cdots\leq Q_{(M)}\), and take \(Q_q=Q_{(\lceil Mq\rceil)}\) for \(0<q\leq1\). This convention is stated explicitly because software packages can use different interpolation rules.
Equal-tail intervals are easy to calculate and interpret, but need not be the shortest interval containing the stated posterior probability. For a skewed or multimodal distribution they may also span a region with relatively little posterior density. A highest-posterior-density interval is an alternative when its definition and computation are appropriate, but for a multimodal distribution a single interval can conceal separated regions of high posterior probability. In such cases, show the distribution or report distinct regions and relevant probabilities as well as an interval.
Worked Examples
Worked Example: Posterior Inhibition of a Pathway
In a hypothetical calibration, suppose eight illustrative posterior draws are used to calculate inhibition of a pathway at 24 hours under one fixed dosing regimen. Define inhibition as \(I=1-R_{\mathrm{treated}}/R_{\mathrm{control}}\), where both responses are positive and refer to the same time and model scenario. The resulting draws, expressed as proportions, are \(0.42, 0.55, 0.48, 0.61, 0.51, 0.58, 0.46, 0.72\). These small numbers make the arithmetic transparent; they are not sufficient for a precise practical analysis.
Sorted, the draws are \(0.42, 0.46, 0.48, 0.51, 0.55, 0.58, 0.61, 0.72\). The empirical median is the fourth ordered draw because \(\lceil8(0.50)\rceil=4\), so the median inhibition is \(0.51\), or 51%. Four of eight draws are at or below \(0.51\), giving an empirical CDF value of \(4/8=0.50\).
The central 50% equal-tail interval uses \(Q_{0.25}=Q_{(2)}=0.46\) and \(Q_{0.75}=Q_{(6)}=0.58\), since \(\lceil8(0.25)\rceil=2\) and \(\lceil8(0.75)\rceil=6\). Thus this illustrative interval is 46% to 58% inhibition. The estimated posterior probability that inhibition exceeds 50% is \(5/8=0.625\), because five draws (\(0.55,0.61,0.51,0.58,0.72\)) are strictly greater than \(0.50\). In context, the draw-based estimate of the posterior probability of more than 50% inhibition at 24 hours is 0.625 under this model and regimen. The small draw count makes that estimate imprecise.
Worked Example: Time Above a Target
For a hypothetical drug-response analysis, define time above target as the total number of hours a model-predicted biomarker response exceeds a prespecified threshold during a 12-hour observation window. Suppose six illustrative posterior draws give durations \(0.8, 1.4, 1.9, 2.1, 2.8,\) and \(3.6\) hours. Assume the trajectory calculation uses the same threshold and time window for every draw, and that the solver resolves threshold crossings accurately.
The sorted durations are already in order. The empirical median is \(Q_{(3)}=1.9\) hours, since \(\lceil6(0.50)\rceil=3\); three of six values are at or below 1.9 hours, so the empirical CDF there is \(3/6=0.50\). The central 50% interval is \(Q_{(2)}\) to \(Q_{(5)}\), or 1.4 to 2.8 hours: \(\lceil6(0.25)\rceil=2\) and \(\lceil6(0.75)\rceil=5\).
For the posterior probability that time above target exceeds 2 hours, the strict inequality selects \(2.1, 2.8,\) and \(3.6\) hours. The estimate is \(3/6=0.50\). A duration equal to 2 hours, if present, would not count for this event; a question using “at least 2 hours” would instead use \(Q\geq2\). The threshold, interval, and strictness of the event are part of the quantity being summarized, not details to leave implicit.
Worked Example: Biomarker Recovery Time
Consider a hypothetical biomarker with baseline value \(B_0\). Define recovery time as the first time after the post-dose nadir at which the model-predicted biomarker reaches at least \(0.90B_0\). Suppose all six illustrative posterior trajectories cross this threshold within the simulation window, producing recovery times of \(18,20,21,23,27,\) and \(35\) hours. The same baseline definition and threshold are used for each draw.
The empirical median is the third ordered value, \(Q_{(3)}=21\) hours, because \(\lceil6(0.50)\rceil=3\). The central 50% interval is \(Q_{(2)}\) to \(Q_{(5)}\), or 20 to 27 hours. Four draws (18, 20, 21, and 23 hours) are below 24 hours, so the estimated posterior probability of recovery before 24 hours is \(4/6=0.667\). The result is a probability about the model-defined recovery time under the specified scenario, not a guarantee that an individual patient will recover by that time.
If a trajectory does not reach the threshold by the simulation horizon, its recovery time is not known to be infinite: it is only known to exceed that horizon, unless the model defines behavior beyond it. Do not discard such draws and summarize only those that crossed. Report the posterior probability of recovery by the horizon, or use a carefully defined censored-quantity summary appropriate to the question.
Why the Transformation Must Be Applied Draw by Draw
The posterior distribution of a derived quantity is generally not obtained by inserting a posterior mean parameter into the model. The reason is that nonlinear functions do not usually commute with averaging. A simple justification uses the inhibition function above. If the treated response \(R\) varies across posterior draws while control response is fixed at 1, then inhibition is \(g(R)=1-R\); in this linear case, \(g\) of the mean equals the mean of \(g\). But if inhibition is calculated from a ratio of uncertain responses, \(g(R_t,R_c)=1-R_t/R_c\), generally \[ 1-\frac{E(R_t\mid y)}{E(R_c\mid y)} \ne E\left(1-\frac{R_t}{R_c}\,\middle|\,y\right). \] The discrepancy can arise from nonlinearity and posterior dependence between the responses. For a trajectory-derived time or threshold crossing, the operation is more nonlinear still: small parameter changes can shift a crossing time substantially or remove the crossing altogether. Therefore, calculate the mechanistic quantity for each joint parameter draw and summarize those values.
Write down the mathematical definition, units, dosing regimen, population profile, time window, and any threshold or baseline convention.
Use parameter values from the same draw together. Combining separate marginal draws can destroy posterior dependence and produce quantities that do not reflect the fitted model.
Evaluate the response or trajectory and apply the definition to obtain one quantity per draw. Check units, numerical solver accuracy, and whether the quantity is defined for every draw.
Report a location summary, a credible interval, and probabilities tied to clinically or mechanistically meaningful thresholds. State the scenario and distinguish posterior uncertainty from Monte Carlo error.
Diagnostics, Alternatives, and Limitations
A posterior summary is only as trustworthy as the draws and the derived-quantity calculation. The earlier material on MCMC convergence, effective sample size, and multimodality applies to transformed quantities as well as to parameters. A chain that has not explored relevant posterior regions cannot be repaired by calculating a median or interval. Check convergence and mixing for the derived quantity itself: a parameter may appear adequately explored while a threshold time, tail probability, or rare event remains poorly estimated.
For an event probability estimated from \(M\) independent draws, the approximate Monte Carlo standard error is \(\sqrt{\hat p(1-\hat p)/M}\). For correlated MCMC draws, \(M\) is not the number of independent observations; autocorrelation reduces information, so use an effective sample size or an appropriate chain-based Monte Carlo error estimate. This error measures simulation precision in estimating a posterior summary. It is distinct from posterior uncertainty, represented by the spread of the transformed draws.
- Check the definition and units. A pathway inhibition expressed as a proportion is not a concentration; time above target is measured in time units. Confirm the same threshold, baseline, dosing scenario, and observation window are used across draws.
- Inspect unusual values and boundaries. A proportion may be constrained to \([0,1]\); a negative duration suggests a coding or definition error. Very long recovery times may indicate a genuine slow-recovery tail, an inadequate simulation horizon, or a numerical problem. Investigate rather than deleting inconvenient draws.
- Check numerical stability. Repeat selected trajectory calculations with tighter solver tolerances or a finer time grid. If a reported threshold crossing changes materially, the numerical calculation is not yet precise enough for the summary.
- Look for multiple regions or skewness. A posterior mean and standard deviation can obscure separated modes or a long tail. Plot transformed draws or their distribution, and report probabilities or quantiles that directly answer the question.
A plug-in estimate evaluates the quantity at a representative parameter value, such as the posterior mean or median. It is compact, but does not describe uncertainty and can be misleading for nonlinear functions. A mean and standard deviation summarize a distribution efficiently when it is reasonably symmetric and unimodal, but can be inadequate for bounded, skewed, or multimodal quantities. Quantiles and posterior event probabilities are often more interpretable for QSP outputs. None of these summaries addresses whether the model is biologically adequate; they are conditional on the specified model and data.
Keep posterior parameter uncertainty separate from observation or population variability. If \(g\) is calculated from a deterministic model trajectory for each parameter draw, its distribution describes uncertainty induced by the parameters (and any random effects explicitly included in the draw). Adding observation error or simulating new individual outcomes answers a predictive question and changes the quantity being summarized. State clearly which source of variation is included.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions below.
- Six sorted posterior draws for recovery time are 2, 3, 3, 4, 6, and 9 days. Under the empirical inverse-CDF convention, what is the median, and what fraction of draws are at or below it?
- From 200 posterior draws, 146 give pathway inhibition above 0.40. Estimate the posterior probability of inhibition above 0.40.
- Why should a time-above-target quantity be calculated for each posterior trajectory rather than from the trajectory at the posterior mean parameters?
- A recovery-time simulation ends at 72 hours, and some trajectories have not reached the recovery threshold by then. What should be reported rather than simply omitting those draws?
Solutions. (1) The median index is \(\lceil6(0.5)\rceil=3\), so the empirical 50th percentile is 3 days. Three of six draws are at or below 3 days, giving empirical probability \(3/6=0.50\); the repeated value means the empirical CDF jumps across 3, but its value at 3 is 0.50. (2) The estimate is \(146/200=0.730\), or 73.0%, of posterior draws. (3) The threshold duration is a nonlinear function of the trajectory and its crossing times; in general, calculating it after averaging parameters or trajectories does not give the posterior distribution of the duration. (4) Report the probability of recovery by 72 hours and the probability of not having recovered by 72 hours, or use a clearly defined method that accounts for times known only to exceed the horizon. Do not condition the summary silently on recovery having occurred within 72 hours.