From Simulated Profiles to Time-Specific Uncertainty
A PBPK simulation can produce a concentration at every time for each sampled set of inputs. The resulting collection of profiles describes uncertainty in the predicted concentration-time trajectory under the chosen simulation design. A useful summary is a band plotted around the trajectory: at each time, its lower and upper edges show quantiles of the simulated concentrations at that time.
The key qualification is at that time. A pointwise band summarizes a separate marginal distribution at each fixed time. It does not, by itself, say that a specified proportion of complete simulated profiles stays inside the two edges for the whole observation period. That distinction matters whenever a figure is read as evidence about an entire trajectory rather than about concentrations at particular times.
The interpretation depends on what was sampled. If uncertain parameters were sampled, the band reflects that parameter uncertainty as represented by their sampling distributions. If the design also samples between-subject variability, the distribution additionally reflects that variability. The band includes only sources represented in the simulation. It does not automatically include measurement error, residual variability, or future observations unless those are explicitly part of the simulated quantity.
For continuous distributions, the probability that \(C(t)\) lies between the two central quantiles is \(1-\alpha\), at a fixed time \(t\). With discrete or otherwise tied output, coverage can be greater than \(1-\alpha\), depending on the quantile convention. Neither statement is a probability about the unknown model being correct. The interval describes uncertainty in the simulated prediction under the stated model and input assumptions.
Constructing the Empirical Band
In practice, the predictive distribution is usually represented by a finite set of simulated concentrations. At each selected time, sort the \(M\) outputs from smallest to largest. One transparent empirical convention is the nearest-rank quantile: the estimated \(p\)-quantile is the observation with rank \(\lceil pM\rceil\) in that sorted list. Other quantile interpolation conventions are also used, so a report should identify the convention when it could affect the result.
For a central 90% band, set \(\alpha=0.10\) and take the 5th and 95th percentiles at every time. The samples are sorted separately at each time. The simulation index attached to a lower or upper quantile can differ across times: the profile supplying the lower edge at one time need not supply the lower edge at another.
Worked Example: Building a Central 90% Band
Suppose 100 hypothetical PBPK simulations produce concentrations at 1, 4, and 12 hours. The table gives the 5th, 50th, and 95th empirical percentiles at each time, in mg/L, using nearest ranks. Thus the 5th percentile is the 5th sorted value and the 95th percentile is the 95th sorted value.
| Time (h) | 5th percentile (mg/L) | 50th percentile (mg/L) | 95th percentile (mg/L) |
|---|---|---|---|
| 1 | 2.1 | 3.0 | 4.6 |
| 4 | 5.8 | 9.2 | 15.1 |
| 12 | 1.4 | 3.1 | 8.7 |
At 4 hours, the central 90% pointwise interval is \([5.8,15.1]\) mg/L. Under the simulated predictive distribution, it describes the middle 90% of concentrations at that fixed time, subject to empirical quantile approximation. At 12 hours, the corresponding interval is \([1.4,8.7]\) mg/L. Plotting the 5th and 95th percentile values across time gives the two edges of the pointwise band; the median values can be plotted as a separate summary curve.
The width at 4 hours is \(15.1-5.8=9.3\) mg/L, while the width at 12 hours is \(8.7-1.4=7.3\) mg/L. These are time-specific spreads in the simulated predictions. They are not standard errors of the median, nor do they alone measure Monte Carlo error in the estimated percentiles.
Why Pointwise Coverage Is Not Whole-Curve Coverage
A concentration-time profile is a collection of values that are dependent across time. A pointwise interval only imposes a coverage statement on one member of that collection at a time. Even if every time-specific interval has 90% coverage, the probability that one complete simulated profile lies inside every interval can be less than 90%.
The reason can be seen directly from the definition of a quantile. At a fixed time \(t\), for a continuous distribution, the lower quantile leaves probability \(\alpha/2\) below it and the upper quantile leaves probability \(\alpha/2\) above it. Therefore the probability between them is \(1-\alpha\). This calculation concerns the random variable \(C(t)\) alone; it does not calculate the probability of satisfying the interval conditions at multiple times simultaneously.
The distinction also follows from elementary probability bounds. If \(A_j\) is the event that the concentration is inside its pointwise interval at time \(t_j\), and each of two times has 90% coverage, then \(P(A_1)=P(A_2)=0.9\). For independent coverage events, \(P(A_1\cap A_2)=0.9\times0.9=0.81\). Without independence, the joint probability is not determined by the two marginal probabilities; it can be as low as 0.8 or as high as 0.9. For many time points, the dependence structure and number of points further affect the chance of covering all of them.
Worked Example: Interpreting Coverage at Two Times
Consider a hypothetical model with central 90% pointwise intervals at 2 and 10 hours. Suppose the events of falling inside the two intervals are independent, and each has probability 0.90 under the predictive distribution.
The probability of being inside both intervals is \(P(A_1\cap A_2)=P(A_1)P(A_2)=0.90\times0.90=0.81\). Thus, in this hypothetical case, 81% of complete two-time outcomes are inside both intervals, even though each individual time point has 90% coverage.
The independence assumption is only for this calculation; concentrations across time in PBPK profiles are generally dependent. Without that assumption, the marginal 90% statements alone are insufficient to give the joint coverage. The example shows why a pointwise 90% label must not be read as a 90% guarantee for an entire profile.
A band can still be useful for displaying how marginal prediction uncertainty changes over time. Its shape can reveal times when simulated concentrations are more dispersed or when uncertainty is asymmetric. But its edges are assembled from separate time-specific distributions, not necessarily from two physically realizable profiles. In particular, a lower-edge point at one time and a lower-edge point at another time may come from different simulations.
Assumptions, Diagnostics, and Alternatives
The calculation is straightforward, but interpretation relies on several conditions. The PBPK model must be evaluated at the requested times using consistent units and a consistent definition of concentration. The sampled inputs and any dependence among them must represent the intended predictive target. Simulations must be valid model outputs rather than failed runs silently treated as ordinary values. Finally, the Monte Carlo sample must be large enough for the quantiles of interest to be estimated with useful precision.
Diagnostics should examine both the profiles and the quantile calculations. Check that all simulations use the same time grid or that interpolation to a common grid is justified; verify that concentrations and time units are consistent; and inspect sorted outputs at selected times for extreme values, ties, or unexpected multimodality. Plot individual simulated profiles as well as the quantile edges when feasible. A suspicious kink or crossing in quantile curves can be a genuine feature of changing marginal distributions, but it can also signal inadequate time resolution, inconsistent interpolation, or too few simulations. It is not automatically a model error.
Quantile estimates can be particularly unstable in the tails when \(M\) is small. Increasing the number of simulations may help, but the number alone does not certify precision. Repeated batches or other Monte Carlo diagnostics can reveal whether estimated quantiles fluctuate materially. A narrow-looking band is not evidence of precise quantile estimation: it describes the spread of model outputs, whereas Monte Carlo error describes uncertainty in the numerical estimate of that spread.
A pointwise band is preferable to reporting isolated prediction intervals when the goal is to display marginal uncertainty over a concentration-time profile. A table of intervals at prespecified times can communicate the same time-specific information more compactly. A mean plus or minus a standard deviation is a different summary: it is not generally a central prediction interval, especially for skewed or multimodal output distributions. A simultaneous band for an entire curve targets a different coverage question and requires a method designed for joint coverage; pointwise quantiles do not become simultaneous merely by being drawn together.
Worked Example: A One-Sided Pointwise Limit
Suppose 100 hypothetical simulations give concentrations at 8 hours, with the 95th sorted value equal to 7.2 mg/L. Using the nearest-rank convention, the empirical 95th percentile is \(7.2\) mg/L because \(\lceil0.95\times100\rceil=95\).
A one-sided pointwise upper limit at 8 hours is therefore \(7.2\) mg/L. Under the simulated predictive distribution, about 95% of concentrations at 8 hours are at or below this limit, with the empirical approximation and quantile convention understood. This statement is useful if the question concerns a concentration threshold at 8 hours. It does not assert that 95% of complete profiles remain below the corresponding upper quantile curve at every time.
The lower and upper edges of a central 90% interval would instead use the 5th and 95th percentiles, respectively. The 95th percentile here supplies the upper edge of that central interval as well, but the one-sided interpretation does not require or imply an accompanying 5th-percentile lower edge.
Common Pitfalls
- Calling a pointwise band a whole-profile guarantee. The wrong reasoning is that a 90% band contains 90% of entire curves because it is drawn across time. The correction is that 90% refers to coverage at each fixed time; joint coverage across times is a separate probability.
- Assuming each band edge is one simulated profile. The lower quantile at each time is calculated from that time's sorted outputs. Different simulations can determine the edge at different times.
- Leaving the uncertainty target unstated. Quantiles reflect the inputs and variability represented in the simulation. State whether the target includes parameter uncertainty, subject variability, or other sampled components, and do not imply omitted sources are included.
- Treating the band width as Monte Carlo error. The band summarizes predictive spread. Monte Carlo error concerns how accurately a finite simulation estimates that spread or its quantiles.
- Reporting percentiles without units or a quantile convention. State the concentration and time units, coverage level, and, when relevant, how empirical quantiles were calculated.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions below.
- With 200 simulations, what nearest-rank observations estimate the 2.5th and 97.5th percentiles at a fixed time?
- A central 90% pointwise interval at 6 hours is \([1.2,4.8]\) mg/L. State its interpretation and one conclusion it does not support.
- If two independent time-specific coverage events each have probability 0.95, what is the probability that both occur?
- Why might the lower edge of a pointwise band fail to trace any single simulated PBPK profile?
- Name two features of the simulation or output that should be checked before interpreting a pointwise band.
Solutions. (1) The ranks are \(\lceil0.025\times200\rceil=5\) and \(\lceil0.975\times200\rceil=195\), so the 5th and 195th sorted concentrations estimate the quantiles. (2) Under the specified predictive distribution, the interval describes the middle 90% of concentrations at 6 hours, approximately for a finite simulation. It does not establish 90% coverage of complete concentration-time profiles. (3) Independence gives \(0.95\times0.95=0.9025\), or 90.25%. (4) Each time is sorted separately, so different simulations can supply the lower quantile at different times. (5) Suitable checks include whether the sampled inputs represent the intended uncertainty target, whether concentration and time units are consistent, whether simulations failed, whether outputs at different times are comparable, and whether the Monte Carlo sample is adequate for the chosen quantiles.