From Input Distributions to PBPK Predictions
A PBPK model maps a set of physiological, physicochemical, and drug-specific inputs to predicted quantities such as concentration at a specified time, area under the concentration-time curve, or clearance. When inputs are uncertain, one prediction from a single parameter set does not show how that uncertainty affects the output. Sampling repeatedly from the specified input distributions and running the model for every sampled set provides a numerical approximation to the output distribution.
The aim is not to create arbitrary variation around a preferred parameter set. Each input distribution should represent the uncertainty or variability that the analysis intends to propagate, and its units and support must match the parameter. Earlier material on measurement error, biological variability, and knowledge uncertainty distinguishes these sources; here, the sampling procedure treats the supplied distributions as given. Unless a joint distribution specifies otherwise, the examples below use independent draws as an explicit simplifying assumption.
The model, output definition, input distributions, and sampling design should be fixed before interpreting the results. For example, “peak concentration” is incomplete unless its dose, route, population, and time horizon are specified. If the model is deterministic, identical parameter sets and settings should give identical outputs, apart from documented numerical tolerances. Randomness in the sampling process is not the same as randomness in the PBPK model.
Sampling a Parameter Set
A basic Monte Carlo procedure draws each parameter according to its specified distribution, assembles the draws into a vector, and evaluates the model. For an independent-input design, each parameter is sampled independently across inputs and simulation replicates. That assumption is consequential: if the inputs are dependent, independently drawing their marginal distributions does not reproduce their joint distribution. The appropriate dependent sampling construction is a separate issue, so independence should never be presented as an innocuous computational default.
Name the model output, dose and simulation conditions, and whether the intended inputs represent parameter uncertainty, population variability, or another stated source.
For each input, record its distribution, parameterization, units, and allowable values. A positive parameter needs a distribution that cannot generate negative values; a fraction must remain between zero and one.
Draw \(N\) parameter sets under the stated sampling design and run the PBPK model once for each set. Retain the input set, output, and run status together.
Verify bounds and run completion, then summarize the output distribution and assess whether Monte Carlo error is small enough for the intended use.
Distribution-based transformations are often useful. For a positive parameter with a lognormal specification, one may draw \(Z\sim N(0,1)\) and set \(\theta=\exp(\mu_{\log}+\sigma_{\log}Z)\), where \(\mu_{\log}\) and \(\sigma_{\log}\) are the mean and standard deviation on the log scale. For a uniform input on \([a,b]\), draw \(U\sim\mathrm{Uniform}(0,1)\) and set \(\theta=a+(b-a)U\). These transformations preserve the specified support. They do not establish that the distribution itself is scientifically appropriate; that choice must come from the input evidence and modeling rationale.
Worked Example: Drawing a Positive Intrinsic Clearance
Suppose a hypothetical PBPK analysis specifies intrinsic clearance as lognormal with median \(12\ \mu\mathrm{L}/\mathrm{min}/10^6\) cells and geometric standard deviation \(1.5\). For this parameterization, \(\mu_{\log}=\ln(12)\) and \(\sigma_{\log}=\ln(1.5)\). Assume one standard normal draw is \(Z=0.5\).
The transformed draw is
This draw is positive, as required by the lognormal distribution. A different standard normal draw gives a different clearance value, and the resulting values are passed into the model one parameter set at a time. This calculation shows how to draw from the stated distribution; it does not imply that a lognormal distribution is automatically suitable for every clearance.
Worked Example: Sampling a Bounded Fraction
Suppose, hypothetically, that a fraction \(F\) is specified as uniform from \(0.20\) to \(0.60\), and one draw from \(\mathrm{Uniform}(0,1)\) is \(U=0.25\). The inverse transformation gives
The sampled fraction \(0.30\) lies within the specified interval. More generally, because \(0\leq U\leq1\), \(a+(b-a)U\) lies between \(a\) and \(b\) whenever \(a\leq b\). The bounded support is therefore respected by construction, subject to checking the implementation and any subsequent transformations.
What the Sample Summaries Mean
For \(N\) successful runs, the simulated mean output is \(\bar{Y}=\sum_{i=1}^{N}Y_i/N\). The empirical quantiles describe the spread of outputs induced by the sampled input distributions and the model. For example, the 2.5th and 97.5th percentiles are an approximate central 95% interval under that input distribution and model. They are not automatically a 95% confidence interval for a population parameter, nor do they include uncertainty that was omitted from the specified inputs or model.
Monte Carlo sampling also introduces finite-sample error: another random batch of the same size would generally produce slightly different summaries. For independent, identically distributed outputs with finite variance, the sample mean has standard error estimated by \(s_Y/\sqrt{N}\), where \(s_Y\) is the sample standard deviation of the outputs. This estimates simulation error in the estimated mean; it is not the spread of the PBPK output distribution itself.
Worked Example: Assessing Monte Carlo Error and Run Count
In a hypothetical pilot run, five parameter sets produce predicted peak concentrations of \(5,\ 7,\ 6,\ 8,\) and \(4\ \mathrm{mg/L}\). Their mean is
The squared deviations from the mean sum to \(1+1+0+4+4=10\), so the sample standard deviation and estimated standard error are
This pilot is far too small to establish a precise output distribution, but it can illustrate a run-count calculation. If a larger pilot estimates \(s_Y=1.8\ \mathrm{mg/L}\), and the chosen target is an estimated standard error no greater than \(0.10\ \mathrm{mg/L}\) for the output mean, the approximate requirement is
Choosing \(N=400\) would give an estimated standard error of about \(1.8/\sqrt{400}=0.090\ \mathrm{mg/L}\), if the pilot spread is representative. This calculation targets precision for the mean only. It does not guarantee equally precise tail quantiles, and it depends on the independent-run and finite-variance assumptions used in the justification.
Bounds, Failures, and Reproducibility
A sampling run should retain enough information to diagnose unexpected results. At minimum, check that sampled values are finite, have the intended units, and lie within their specified supports. Check model outputs for missing values, numerical failures, and implausible extremes. A failed run should not silently disappear: report its frequency and investigate whether failure reflects a coding problem, a numerical limitation, or an input region where the model is undefined. Dropping failed runs without accounting for them can change the output distribution being summarized.
Bounds checks are not the same as a full biological plausibility assessment. A value can satisfy its own marginal range while forming an implausible combination with other inputs. Earlier material on plausibility checks addresses such combinations; this tutorial’s sampling checks establish whether draws obey the specified input distributions and whether the model evaluated them successfully.
For reproducibility, record the distribution definitions and parameter values, input units, sampling design, sample size, model version, dose and simulation settings, random-number generator and seed, and any transformations or rounding. A seed alone is insufficient if the generator, software, draw order, or model implementation changes. Saving the generated parameter sets and outputs provides a direct audit trail.
A common alternative to simple Monte Carlo is Latin hypercube sampling, which stratifies each marginal distribution so the sampled values cover its range more evenly for a fixed run count. It can improve coverage for some summaries, but it is not a guarantee of better precision for every PBPK output and does not by itself encode dependence among parameters. The earlier material on the Delta Method offers a local approximation based on sensitivities; direct sampling instead evaluates the model at many input sets and can represent distribution shape more fully, at greater computational cost. Neither approach repairs misspecified input distributions or an inadequate model.
Common Pitfalls and Practical Limits
- Confusing output spread with Monte Carlo error. The standard deviation \(s_Y\) describes variation across propagated inputs; \(s_Y/\sqrt{N}\) estimates simulation error in the sample mean.
- Calling an empirical percentile interval a confidence interval. It summarizes the modeled output distribution under the stated inputs; its interpretation depends on what those distributions represent.
- Ignoring support or units. A bounded fraction, positive clearance, and volume cannot be sampled interchangeably. Verify that the distribution and transformation match both the physical parameter and its units.
- Assuming independent draws without saying so. Independent marginal sampling is valid only for an analysis that specifies or deliberately assumes independence. It can misrepresent the input distribution when dependence matters.
- Stopping when a run completes. Successful execution does not establish adequate sample size. Compare Monte Carlo error with the precision needed for the intended summary, and inspect tail summaries separately if they drive a decision.
Sampling is a numerical propagation method, not evidence that the chosen distributions are correct. Results remain conditional on the input distributions, model structure, implementation, and simulation conditions. A larger run can reduce Monte Carlo error, but it cannot remove uncertainty that was not represented or correct structural model error.
Check Your Understanding
Answer each question, then compare with the solutions.
- A model output has \(s_Y=2.4\ \mathrm{mg/L}\) across independent runs. Approximately how many runs are needed for an estimated standard error of \(0.12\ \mathrm{mg/L}\) for its mean?
- For a uniform input on \([0.1,0.5]\), what value results from \(U=0.75\)? Does it satisfy the specified bounds?
- In one sentence, distinguish the output standard deviation from the estimated standard error of the sample mean.
- Why does recording a random-number seed alone not guarantee reproducibility?
Solutions. (1) \(N\geq(2.4/0.12)^2=20^2=400\) runs. (2) \(0.1+(0.5-0.1)(0.75)=0.40\), which is within \([0.1,0.5]\). (3) The standard deviation describes propagated output spread, whereas the estimated standard error \(s_Y/\sqrt{N}\) describes simulation error in the estimated mean under independent sampling. (4) Reproduction also depends on the generator, software and model versions, draw order, input distributions, and simulation settings.