From Uncertain Inputs to an Exposure Distribution
A PBPK simulation maps physiological and drug-related inputs to predicted concentrations and other outputs. When those inputs are uncertain, one simulation at a single selected parameter vector gives only one conditional prediction. Forward Monte Carlo propagation instead samples plausible input vectors, runs the PBPK model for each vector, and summarizes the resulting output values. The central task is not just to run the model repeatedly: it is to preserve the intended joint input distribution and apply the same defined exposure metric to every run.
For example, a question might concern plasma concentration at a specified time after a dose, in a specified population or scenario. The result of propagation is a distribution of that metric under the specified input uncertainty. It is not automatically a distribution that includes between-subject variability, and it does not include model-structure uncertainty unless those sources have been represented explicitly.
A common model representation is a system of differential equations, such as \(d\mathbf{x}/dt=\mathbf{f}(t,\mathbf{x};\boldsymbol{\theta})\), with specified initial conditions and dosing or other inputs. The metric must be defined before summarizing results: for instance, concentration in a named compartment at a fixed time, with concentration units and time stated. Different choices of metric answer different questions, even when they use the same simulated trajectory.
Conditions for a Meaningful Propagation
Propagation answers a well-defined question only when the input distribution, model scenario, and output metric are all explicit. The joint distribution should represent the uncertainty being studied, including the dependence among inputs. A sampled vector must be treated as a coherent set: its entries belong together under that joint distribution. Drawing each marginal independently is appropriate only if independence is part of the specified joint distribution.
- Input definition: each draw is within the intended domain of the PBPK model, and the marginal distributions and dependence structure have a defensible interpretation.
- Fixed scenario: dose, route, schedule, population definition, and other non-sampled settings are held as specified. Any setting intended to vary must be represented in the sampled inputs.
- Consistent simulation: the same model equations, initial conditions, dosing rules, and numerical standards are applied across draws. The numerical solver must adequately resolve the model over the time range needed for the metric.
- Defined output: the metric, its units, and how it is extracted from a trajectory are stated in advance and applied consistently.
If these conditions hold, differences among the \(Y_i\) reflect the modeled consequences of input uncertainty, subject to simulation error. They do not prove that the model is biologically correct. A narrow output distribution can arise because the sampled inputs have little influence on the chosen metric, but it can also arise because relevant uncertainty was omitted or represented too narrowly.
Why the Output Sample Represents a Distribution
The mathematical basis is the induced, or pushforward, distribution of the metric. For any threshold \(y\), its cumulative distribution function is
Here, \(\mathbf{1}\{\cdot\}\) equals 1 when the condition in braces is true and 0 otherwise. If the \(\boldsymbol{\theta}_i\) are independent draws from the specified joint distribution, then each \(Y_i\) has the induced distribution of \(Y\). For example, the proportion of outputs at or below \(y\) estimates \(F_Y(y)\). If the output has a finite mean, the sample average estimates that mean:
Under independent sampling, \(E(\bar{Y})=E(Y)\): the estimate is unbiased, provided the expectation exists and each simulation returns the intended metric. Its random error decreases as more draws are used, though the speed and adequacy of that decrease depend on the output distribution. This result justifies averaging simulated outputs; it does not imply that a small collection of runs describes every tail feature accurately.
Worked Example: Summarizing a Fixed-Time Concentration
Suppose a hypothetical PBPK analysis samples six joint parameter vectors from a specified uncertainty distribution. Each simulation uses the same dose and scenario, and the metric is plasma concentration at 4.0 h after dosing, in mg/L. The following are hypothetical solver outputs; the six rows are six complete joint draws, not separately recombined marginal values.
| Draw | Clearance input (L/h) | Fraction unbound | Concentration at 4.0 h (mg/L) |
|---|---|---|---|
| 1 | 8.0 | 0.080 | 2.1 |
| 2 | 7.5 | 0.090 | 2.4 |
| 3 | 7.0 | 0.100 | 2.8 |
| 4 | 6.5 | 0.110 | 3.0 |
| 5 | 6.0 | 0.120 | 3.7 |
| 6 | 5.5 | 0.130 | 4.0 |
The Monte Carlo estimate of the mean concentration is
The sample standard deviation describes the spread of these six outputs around their sample mean:
Thus, in this hypothetical computation, the estimated mean is 3.00 mg/L and the output sample standard deviation is 0.735 mg/L. The minimum and maximum simulated values are 2.1 and 4.0 mg/L. Six draws are enough to illustrate the calculations, not to establish a reliable description of the tails or a stable final summary.
Joint Draws Matter When the Model Is Nonlinear
The same marginal distributions can produce different output distributions when their dependence differs. This matters because PBPK outputs are generally nonlinear functions of several inputs. Replacing a joint sample by separately shuffled or independently drawn columns changes which parameter combinations are sent through the model; it therefore changes the propagated question.
Worked Example: Preserving Dependence Between Inputs
Consider a deliberately simplified hypothetical exposure mapping, used only to make the arithmetic transparent:
where \(F\) is dimensionless and \(CL\) is clearance in L/h, so \(Y\) is in mg/L. Suppose the marginal sampled values are \(F=0.4,0.6,0.8\) and \(CL=2,3,4\) L/h. One possible positively paired joint sample is \((0.4,2),(0.6,3),(0.8,4)\). The output values are
Now pair the same marginal values in the opposite order: \((0.4,4),(0.6,3),(0.8,2)\). The outputs become
Both samples have the same marginal values, and both have mean output \( (1+2+4)/3=2.33\) mg/L in the second pairing versus \(2.00\) mg/L in the first. The first has no output spread; the second has sample standard deviation
The pairings represent different dependence structures and therefore different joint distributions. In an actual PBPK workflow, the model replaces this illustrative mapping, but the principle is the same: sample and pass forward complete vectors drawn from the intended joint distribution. The marginal distributions alone do not determine the distribution of a nonlinear output.
Run, Extract, and Check Each Simulation
Record the scenario, joint input distribution, and precise exposure metric, including compartment, time or other defining feature, and units.
Generate each vector using the specified joint sampling design. Retain each vector as a unit and record its draw identifier.
Apply consistent initial conditions, doses, equations, and numerical settings. Record solver status and any warnings alongside the parameter vector.
Apply the same extraction rule to each trajectory. Check units, numerical plausibility, and model-specific quantities such as mass balance where applicable.
Report the number of completed valid runs, the metric definition, and suitable summaries of the simulated distribution, in context and with units.
Worked Example: Investigating a Failed Simulation
In a hypothetical analysis, the metric is plasma concentration at 6 h after dosing. Five joint parameter vectors are sampled. Four runs complete initially, while the fifth reports a solver failure. The failed run must not be silently omitted: its parameter vector and solver message should be inspected, because a failure could indicate a numerical problem or a draw outside the model's valid domain.
Suppose inspection shows that the fifth draw is within the specified input support, and a rerun with a numerically adequate solver setting completes. The five valid hypothetical outputs are 1.8, 2.0, 2.4, 3.1, and 8.2 mg/L. The mean is
The 8.2 mg/L result is high relative to the other four, but it is not a reason by itself to delete that run. If it is a valid consequence of the sampled inputs, removing it changes the target distribution. The mean of the four lower outputs alone would be \((1.8+2.0+2.4+3.1)/4=2.325\) mg/L, which answers a different, selectively truncated question. If the run cannot be made to complete, the failure should be reported and investigated; the remaining successful runs should not be presented as though they were an unfiltered sample from the intended output distribution.
Choosing Summaries and Recognizing Failure Modes
The output sample can be summarized in several ways, but each summary describes a different feature. The mean and sample standard deviation describe the center and spread on the original scale when those summaries are useful. The range reports only the smallest and largest simulated values and is strongly dependent on the number of draws. Empirical quantiles can describe locations in the simulated distribution, but their accuracy depends on having enough informative output draws. A summary should never replace a clear statement of the metric and the distribution from which it was computed.
Several checks help reveal when propagation is going wrong:
- Input checks: confirm that sampled values respect support constraints and that dependence has not been lost through independent resampling or accidental column shuffling.
- Simulation checks: review failed runs, warnings, implausible discontinuities, and numerical behavior near extreme inputs. Repeat calculations with suitable solver settings when numerical adequacy is in doubt.
- Output checks: verify the metric extraction, time units, concentration units, and compartment labels. A trajectory can be valid while the post-processing code selects the wrong time or state variable.
- Distribution checks: inspect the output values for skew, multiple clusters, or a few influential results. These patterns may be genuine consequences of the joint inputs, but may also reveal coding errors or an input distribution inconsistent with the intended use.
A common mistake is to run the model at the mean input vector and treat that result as the mean output. In general, for a nonlinear model, \(h\{E(\boldsymbol{\theta})\}\ne E\{h(\boldsymbol{\theta})\}\). The correction is to propagate the input draws through the model first, then summarize the resulting outputs. Another mistake is to discard a large valid output as an “outlier” without a prespecified scientific or computational reason. Investigate it; do not change the target distribution merely because the result is inconvenient.
How Forward Propagation Differs from Alternatives
| Approach | What it provides | Main limitation for this question |
|---|---|---|
| Single nominal-parameter simulation | A prediction conditional on one selected input vector | Does not show the output distribution induced by uncertain inputs |
| Local sensitivity analysis | How outputs respond to small input changes near a selected point | Does not by itself represent a full joint input distribution or its nonlinear output consequences |
| Forward sampling through the PBPK model | A simulated sample from the output distribution under the specified input distribution | Requires many model evaluations and inherits limitations in the input distributions and model |
Monte Carlo, Latin hypercube, and stratified designs are alternative ways to select input draws, as covered in the earlier material on sampling. Whatever design is used, the simulation step is the same: each complete draw is run through the PBPK model, and its output is recorded. Forward propagation is limited by the assumed input distributions, the PBPK model and scenario, numerical accuracy, and finite computational effort. It does not repair uncertain or misspecified inputs, nor does it establish that the model captures every relevant biological process.
Check Your Understanding
Use the definitions and workflow above to answer these questions.
- Three hypothetical simulation outputs for a fixed-time concentration are 2, 3, and 7 mg/L. Calculate their Monte Carlo mean.
- Why can separately drawing each uncertain input from its marginal distribution give the wrong propagated distribution?
- A run produces a high concentration but has no solver warning and its inputs are within their specified ranges. Should it automatically be removed? Explain.
- Why is the model prediction at the mean parameter vector generally not a substitute for the mean of propagated predictions?
Solutions. (1) \(\bar{Y}=(2+3+7)/3=4.00\) mg/L. (2) Independent marginal draws impose independence; if the intended inputs are dependent, the resulting combinations and thus the output distribution differ from those under the specified joint distribution. (3) No. A valid high output can be part of the intended distribution; investigate its inputs and simulation, but do not delete it solely because it is high. (4) The model mapping is generally nonlinear, so \(h\{E(\boldsymbol{\theta})\}\) need not equal \(E\{h(\boldsymbol{\theta})\}\); simulate the sampled vectors first and summarize their outputs.