Estimating Rare-Event Probabilities from PBPK Simulations

Why Rare-Event Probabilities Need Special Care

A PBPK simulation may be used to estimate the probability that a predicted exposure exceeds a specified threshold. If the threshold is far into the upper tail of the predictive distribution, only a small fraction of simulated outcomes may exceed it. In that setting, reporting the fraction of exceedances without its Monte Carlo uncertainty can make a very imprecise estimate look definitive.

The target is a probability under a specified simulation setup: the chosen uncertain-input distributions and dependence structure, the PBPK model, the simulation conditions, and the definition of exceedance. It is not automatically the real-world frequency of harm. Direct Monte Carlo quantifies uncertainty from using a finite number of simulations to estimate that model-based probability; it does not by itself account for model misspecification or uncertainty in the chosen input distributions.

Definition: Let \(Y\) be a PBPK prediction from one draw of the specified input distribution, and let \(c\) be a fixed threshold. Define the exceedance indicator \(I_i=1\) if simulated prediction \(Y_i\) exceeds \(c\), and \(I_i=0\) otherwise. The rare-event probability is \(p=P(Y>c)\). From \(n\) independent simulations, let \(K=\sum_{i=1}^{n} I_i\) be the number of exceedances. The direct Monte Carlo estimate is \(\hat{p}=K/n\).

Assumptions and Binomial Uncertainty

The usual binomial calculation requires that each simulation corresponds to a draw from the target distribution and that the resulting event indicators are independent and have the same probability \(p\). Under these assumptions, each \(I_i\) is Bernoulli with probability \(p\), and \(K\sim\operatorname{Binomial}(n,p)\). The simulation outcome for one draw must be classified consistently: for example, the event might be \(Y>c\) or \(Y\ge c\), but the choice should be explicit and applied to every draw.

The estimate \(\hat{p}\) is unbiased, because \(E(\hat{p})=p\). Since the variance of a Bernoulli indicator is \(p(1-p)\), independence gives the variance of their average. Replacing the unknown \(p\) by its estimate yields the conventional estimated Monte Carlo standard error.

E(\hat{p})=p,\qquad \operatorname{Var}(\hat{p})=\frac{p(1-p)}{n},\qquad SE(\hat{p})=\sqrt{\frac{\hat{p}(1-\hat{p})}{n}}.

This standard error describes simulation-to-simulation variation in \(\hat{p}\) under the stated sampling assumptions. It is not the variability of individual PBPK predictions. For a small \(p\), the relative standard error is approximately \(1/\sqrt{np}\), provided \(np\) is not extremely small. Thus, the expected number of events, \(np\), is central: a very large total number of simulations can still give a noisy estimate if only a few events occur.

Key result: With independent, identically distributed direct Monte Carlo draws, the exceedance count is binomial and \(\hat{p}=K/n\) estimates the target probability. The standard error is useful when enough events occur, but a normal approximation based on it can be poor when \(K\) is small or zero.

Worked Example: Zero Exceedances Does Not Mean Zero Probability

Suppose 5,000 independent hypothetical PBPK simulations are run and none exceeds a specified exposure threshold. The estimated probability is \(\hat{p}=0/5{,}000=0\). The plug-in standard error is also zero, which is not evidence that the true model-based probability is zero; it is a failure of that plug-in estimate to express uncertainty when no events were observed.

A useful one-sided 95% exact binomial upper confidence bound for zero events is found by choosing \(p_U\) so that the probability of observing zero events at \(p_U\) is 0.05:

P(K=0\mid p_U)=(1-p_U)^n=0.05,\qquad p_U=1-0.05^{1/n}.

Substituting \(n=5{,}000\), \(p_U=1-0.05^{1/5000}\approx0.000599\). Thus, although no exceedances occurred, the data are compatible at the 95% one-sided confidence level with probabilities up to about 0.000599, or 0.0599% (about 6 per 10,000 simulations). This bound describes repeated-simulation coverage under the binomial model; it does not mean there is a 95% probability that \(p\) lies below this particular bound.

Intervals and the Problem of Few Events

A common quick interval is \(\hat{p}\pm1.96SE(\hat{p})\). It relies on a normal approximation to the sampling distribution of \(\hat{p}\), which is unreliable when the count of events is small. It can even produce a negative lower endpoint or a zero-width interval when \(K=0\). For small counts, a Wilson score interval is generally more defensible than this unadjusted Wald interval. It remains an approximation; an exact binomial interval is another option when exact coverage is required.

For a two-sided 95% Wilson interval, let \(z=1.96\). The lower and upper limits are given by the following expression with the minus and plus signs, respectively. The interval stays within 0 and 1 and does not collapse to zero width when no events are observed.

\frac{\hat{p}+\frac{z^2}{2n}\ \pm\ z\sqrt{\frac{\hat{p}(1-\hat{p})}{n}+\frac{z^2}{4n^2}}} {1+\frac{z^2}{n}}.

Worked Example: Twelve Exceedances in Twenty Thousand Simulations

In a hypothetical run of \(n=20{,}000\) independent simulations, \(K=12\) predicted exposures exceed the prespecified threshold. The probability estimate is

\hat{p}=\frac{12}{20{,}000}=0.000600, \qquad SE(\hat{p})= \sqrt{\frac{0.000600(1-0.000600)}{20{,}000}} \approx0.000173.

The estimated probability is 0.0600%, or 6.00 exceedances per 10,000 simulations. The standard error is about 28.9% of the estimate, so the estimate is not very precise relative to its size. For comparison, the naive normal interval would be approximately \(0.000600\pm1.96(0.000173)\), or 0.000261 to 0.000939. With only 12 events, that interval's normal approximation warrants caution.

Using the Wilson formula, the denominator is \(1+1.96^2/20{,}000\approx1.000192\). The center is \((0.000600+1.96^2/(2\times20{,}000))/1.000192\approx0.000696\), and the half-width is approximately 0.000353. The resulting 95% Wilson interval is approximately 0.000343 to 0.001049, or 0.0343% to 0.1049%. It expresses substantial uncertainty relative to the estimated probability. The interval concerns the probability under this simulation setup, not a range containing 95% of individual exposure predictions.

Planning and Diagnosing a Rare-Event Simulation

Simulation size should be considered in terms of the number of expected events, not just the total number of runs. For small \(p\), the relative standard error is approximately \(1/\sqrt{np}\). A desired relative standard error \(r\) therefore calls for roughly \(np\ge1/r^2\) expected events. This is a planning approximation, not a guarantee: the realized count varies, and confidence-interval width should be checked after simulation.

Worked Example: Planning for a Small Probability

Suppose a preliminary estimate suggests an exceedance probability near \(10^{-5}\). A team wants the relative Monte Carlo standard error to be no more than 20% as a planning target. Using the small-probability approximation, require \(1/\sqrt{np}\le0.20\), so \(np\ge25\). With \(p\approx10^{-5}\), this implies

n\ge\frac{25}{10^{-5}}=2{,}500{,}000.

At this size, the expected event count is 25, and the approximate relative standard error is \(1/\sqrt{25}=0.20\). This is a planning calculation based on an uncertain preliminary probability. If the goal were instead a roughly 20% relative half-width for a normal 95% interval, the approximate requirement would be \(1.96/\sqrt{np}\le0.20\), or \(np\ge(1.96/0.20)^2\approx96.0\), corresponding to about 9.60 million runs. Either target should be followed by inspection of the observed event count and an appropriate interval.

A practical analysis should check more than the final fraction. Verify that every intended draw produced a valid prediction, that failures were not silently removed in a way that changes the target population, and that the threshold and units match the scientific question. Inspect whether the observed count is large enough for the chosen uncertainty summary. Repeatability of results across larger runs or independent random seeds can reveal instability, although agreement between runs does not prove that the model or input distributions are appropriate.

Specify the event.
State the PBPK output, threshold, units, and whether equality counts as exceedance.

Check the sampling design.
Confirm that the direct Monte Carlo draws are independent draws from the target joint input distribution and that each simulation is handled consistently.

Count and estimate.
Record \(K\) and \(n\), calculate \(\hat{p}=K/n\), and report the event count as well as the probability.

Quantify and diagnose uncertainty.
Use a binomial interval suited to the event count; treat zero or very few events as limited information, not proof of negligible risk.

Alternatives, Limitations, and Common Pitfalls

Direct Monte Carlo is attractive because the estimate is simply the observed event fraction and the binomial uncertainty calculation is transparent. Its main limitation for rare events is computational: many simulations may be needed to observe enough exceedances. More specialized approaches, such as importance sampling or splitting methods, can focus computation on the tail. They are not ordinary unweighted binomial estimates: their probability estimates and uncertainty calculations must account for the altered sampling scheme and any weights. They also require checks that the rare-event region is represented correctly. For a straightforward baseline estimate, direct Monte Carlo provides a useful reference when feasible.

The binomial uncertainty calculation does not apply unchanged if event indicators are dependent, have different probabilities, or arise from a sampling design such as stratified or Latin hypercube sampling. Such designs can be useful for uncertainty propagation, but their dependence structure means that treating the total exceedance count as an ordinary binomial count may give an incorrect standard error. Similarly, repeated time points from one simulated individual are not automatically independent replicates. Define the event at the intended simulation-unit level and use an uncertainty method appropriate to the actual design.

Mistake: “No events were observed, so the probability is zero.” The estimate is zero, but finite simulations cannot rule out a small positive probability. Report an upper bound or state that the simulation did not resolve probabilities below its effective precision.
Mistake: reporting only a percentage. “0.06%” hides whether this came from 12 events in 20,000 simulations or a much larger count in many more simulations. Report \(K\), \(n\), the estimate, and an interval or standard error.
Mistake: trusting a symmetric normal interval for a tiny estimate. It may be poorly calibrated or extend below zero. Use a Wilson or exact binomial interval when event counts are small.
Mistake: interpreting Monte Carlo precision as scientific certainty. A narrow interval only indicates precise computation of the probability under the stated PBPK model and input assumptions. It does not establish that those assumptions accurately represent clinical reality.

In short, rare-event reporting should make the event definition and simulation denominator visible, quantify finite-simulation uncertainty, and acknowledge when the observed event count is too small for a stable estimate. A zero or tiny estimate is often a reason to plan more simulations or consider a validated rare-event method, not a reason to claim that exceedance cannot occur.

Check Your Understanding

Answer each question using the binomial model where its assumptions are stated.

In 10,000 independent simulations, 5 outcomes exceed the threshold. What is the direct Monte Carlo estimate?
Why is the plug-in standard error zero when no events are observed, and why should that not be interpreted as certainty?
For an estimated probability near \(2\times10^{-4}\), approximately how many independent simulations are needed for 100 expected events?
Why might the ordinary binomial standard error be inappropriate for an exceedance count from a dependent or stratified sampling design?

Solutions. (1) \(\hat{p}=5/10{,}000=0.000500\), or 0.0500%. (2) Substituting \(\hat{p}=0\) into \(\sqrt{\hat{p}(1-\hat{p})/n}\) gives zero, but this plug-in formula cannot show uncertainty at a boundary estimate; a one-sided upper bound remains positive. (3) Set \(np\approx100\): \(n\approx100/(2\times10^{-4})=500{,}000\) simulations. (4) The binomial variance assumes independent, identically distributed indicators with common probability. Dependence or unequal sampling probabilities can change the variance, so a design-appropriate uncertainty calculation is needed.