Why Reproducibility Is Part of a Monte Carlo Analysis
A Monte Carlo PBPK analysis draws uncertain inputs, runs a model for each draw, and summarizes the resulting predictions. A reported result is more useful when another analyst can establish exactly which inputs produced it, inspect the corresponding model outputs, and rerun the computations. That requires more than saving a final mean or distribution summary: it requires recording how the simulated inputs were generated and preserving their links to model outputs.
The goal here is computational reproducibility: a specified workflow can be rerun to produce the same simulated inputs and, subject to the stated numerical tolerance, the same outputs. This is different from asking whether a simulation has run long enough to estimate a quantity precisely. Reproducibility makes a run traceable; it does not, by itself, establish that the Monte Carlo results are sufficiently accurate for their intended use.
What a Seed Does—and Does Not Do
A pseudorandom-number generator is a deterministic algorithm that produces a sequence of values designed to behave like random draws for the intended use. Its starting state is initialized, often using a seed. Given the same generator implementation, seed, and sequence of generator calls, the same sequence of values is generally produced. Sampling procedures then transform or combine those values to create model inputs.
Consequently, a seed is one part of a reproducibility record, not a complete one. The resulting inputs can change if the software uses a different generator, a changed sampling algorithm, a different call order, or a changed transformation. For example, adding a random draw for a diagnostic before drawing clearance can shift every subsequent value in a single shared stream. Changing software versions can also change how a distribution sampler converts generator values into draws.
A useful formal description is to write the simulated inputs as \(X_i=T(U_i;S)\) and the model output as \(Y_i=f(X_i;M)\). Here \(U_i\) is the generator output for simulation \(i\), \(S\) denotes the sampling specifications and transformations, and \(M\) denotes the model implementation and settings. Reproducing \(Y_i\) requires enough information to reproduce or recover \(U_i\), \(S\), and \(M\), as well as the model execution environment when numerical details matter.
For a deterministic model, the same input and same model execution should give the same result, apart from any explicitly allowed numerical tolerance. If the model itself contains stochastic elements, their random-number settings and their relation to the input simulation must also be recorded. Otherwise, a rerun can differ even when the sampled PBPK inputs are unchanged.
Specify the Entire Sampling Procedure
For every uncertain input, record its identity, units, marginal distribution, parameter values and parameterization, and any transformation used to obtain a physical value. For example, a lognormal distribution may be parameterized by the mean and standard deviation on the log scale or by moments on the original scale; writing only “lognormal” is ambiguous. For dependent inputs, record the dependence representation and how the joint samples were generated. If a bounded transformation, truncation rule, or rejection procedure is used, specify that as well.
Also record the sampling design and its settings: for example, ordinary Monte Carlo, Latin hypercube sampling, or stratification; the requested number of simulations; and any ordering or allocation rules. These approaches may target the same input distributions but do not generally produce the same set of simulated rows. Recording a seed without recording which sampling design used it is insufficient.
For nested sampling, keep the index structure explicit. A row should make clear which parameter-uncertainty draw and which individual-variability draw it represents. Do not flatten a nested design into an unlabelled list if doing so loses those relationships. A stable simulation identifier, plus any relevant outer- and inner-level identifiers, makes the design auditable.
A Practical Record for Each Run
A compact run manifest should identify what was run and how. At minimum, include the run date or release identifier, a unique run ID, the random-number generator and seed or stream settings, the sampling design and specifications, the model version, software and relevant package versions, and the numerical settings that can affect results. Include the input data or their stable references, the units and transformations, and the handling of failed or invalid simulations.
Keep the sampled inputs alongside outputs, not just in transient memory. Use one row per simulation (or a clearly defined record for each nested index combination), with the simulation ID, sampled values, outputs, units, and status. Preserve failed rows with an informative status and error description rather than silently dropping them. A failure may depend on a particular combination of inputs; removing that combination without recording it prevents a complete audit of what the workflow did.
Where possible, save a machine-readable manifest and data table, and retain a stable copy of the model code and configuration. A checksum or version-control identifier can help establish that a file or code version has not changed. A checksum is evidence about file identity, not a substitute for preserving the file itself. Likewise, recording a software version is useful, but an archived environment or container may be needed when exact long-term replay is important.
Record the distributions, parameterizations, dependence structure, transformations, sampling design, simulation count, model settings, and units before interpreting results.
Set the documented generator and seed or stream scheme, assign a unique run ID, and save the software and model identifiers.
Assign stable simulation IDs and retain the sampled inputs, nested indices where applicable, outputs, and success or failure status.
Rerun a small set or the full run from the saved specification, compare inputs and outputs under a stated equality rule, and document any discrepancy.
Worked Examples
Worked Example: Replaying a Small PBPK Input Set
Consider a hypothetical one-compartment calculation used only to illustrate record keeping. A 100 mg dose is given at time \(t=2\) h. Clearance \(CL\) is sampled uniformly from 4 to 6 L/h, and volume \(V\) uniformly from 20 to 30 L. For this simplified calculation, the concentration is \(C(t)=(D/V)\exp(-CLt/V)\). Suppose the run manifest records a particular generator, implementation, seed, sampling procedure, and a generated row with \(u_1=0.25\) and \(u_2=0.60\), where each \(u\) is a uniform draw on 0 to 1.
The predicted concentration is therefore
The saved record should link the simulation ID to both uniform draws, the resulting \(CL\) and \(V\), and the output \(C(2)\), with units and model version. The two uniform values here are hypothetical recorded draws; the example does not assert that a particular seed in a particular software package produces them. A real rerun must identify the generator and sampling implementation that produced its own draws.
Worked Example: Preserving Nested Simulation Indices
Suppose a hypothetical nested run uses two parameter draws and two individual draws per parameter draw. The model output is a concentration calculated by \(C=D/V\), with \(D=100\) mg. The parameter draw determines \(V\), and the individual draw determines a scaling factor \(r\), so the input volume is \(V_{\mathrm{ind}}=Vr\).
For parameter draw 1, let \(V=25\) L; for parameter draw 2, let \(V=30\) L. The two individual factors are \(r=0.8\) and \(r=1.2\). Retaining both indices gives four distinct records:
| Parameter ID | Individual ID | \(V\) (L) | \(r\) | \(V_{\mathrm{ind}}\) (L) | \(C\) (mg/L) |
|---|---|---|---|---|---|
| 1 | 1 | 25 | 0.8 | 20 | 5.00 |
| 1 | 2 | 25 | 1.2 | 30 | 3.33 |
| 2 | 1 | 30 | 0.8 | 24 | 4.17 |
| 2 | 2 | 30 | 1.2 | 36 | 2.78 |
For example, the first row has \(V_{\mathrm{ind}}=(25\ \mathrm{L})(0.8)=20\ \mathrm{L}\) and \(C=(100\ \mathrm{mg})/(20\ \mathrm{L})=5.00\ \mathrm{mg/L}\). The third has \(V_{\mathrm{ind}}=(30\ \mathrm{L})(0.8)=24\ \mathrm{L}\), giving \(C=4.17\ \mathrm{mg/L}\). Keeping both IDs makes it possible to identify which outer parameter draw and inner individual draw produced each result. A single row number without these keys would not express the nested structure.
Replay, Parallel Runs, and Diagnostics
Exact replay is easiest when the workflow is deterministic after the random inputs have been generated. One robust practice is to save the actual sampled inputs and then rerun the PBPK model using those stored inputs. This separates two checks: whether the sampling procedure can be replayed, and whether the model reproduces outputs from a fixed input set. It also provides a direct audit trail if the random-number software later changes.
Parallel execution needs particular care. If several workers are initialized with the same seed and generator, they may generate identical or overlapping streams. Instead, use a documented method for assigning independent streams or substreams, or a generator designed to provide reproducible parallel streams. Record the assignment of streams to workers or simulation IDs. If changing the number of workers changes which random values are assigned to simulations, exact row-by-row replay may require saving the generated inputs or using a stream scheme indexed independently of worker order.
Useful reproducibility checks include verifying that simulation IDs are unique, required inputs and outputs are present, units and dimensions agree with the specification, and the replayed sampled inputs match the retained inputs. For numerical outputs, define in advance whether equality must be bit-for-bit or within an absolute or relative tolerance. Floating-point calculations can differ slightly across hardware, libraries, or execution order. A tolerance should be justified by the use of the result; it should not be chosen after seeing an inconvenient discrepancy.
Investigate unexpected differences rather than simply rerunning until results look similar. Check whether the seed was reset at the correct point, the generator and package versions match, random calls occur in the same order, transformations and parameterizations are unchanged, parallel streams are assigned consistently, and model settings or input files have changed. Also confirm whether failed simulations were retained or filtered differently.
Worked Example: Finding a Pairing Error in an Audit
A hypothetical saved file contains two input-output records for \(C=D/V\), with \(D=100\) mg. The records are ID A: \(V=20\) L, \(C=5.00\) mg/L; and ID B: \(V=25\) L, \(C=4.00\) mg/L. During an audit, a rerun produces the same inputs, but a merged report lists ID A with 4.00 mg/L and ID B with 5.00 mg/L.
Recalculate by ID. For A, \(100/20=5.00\) mg/L; for B, \(100/25=4.00\) mg/L. The model outputs themselves are present, but the report has swapped their identifiers. Comparing only the set of output values would miss this error because both 5.00 and 4.00 appear in both files. Comparing the complete input-output pair by stable ID detects it. The correction is to repair the join or indexing logic and regenerate the report, not to change the model or sampling seed.
Common Pitfalls and Limits
A frequent mistake is to reason, “The seed is saved, so the analysis is reproducible.” A seed alone does not specify the generator, sampling algorithm, distributions, transformations, call order, model version, or parallel-stream assignment. The correction is to treat the seed as one field in a run manifest and to preserve the actual sampled inputs and outputs.
Other warning signs include different numbers of rows between the original and replay, duplicate simulation IDs, missing failure records, changed units, inputs that violate stated bounds, or outputs that cannot be traced back to a row. Such findings may indicate a data-handling or workflow change, even when summary results appear similar. A changed output need not mean the biological model is wrong; first establish whether the same input, model, and execution conditions were compared.
There are limits to what can be promised. A workflow may be reproducible only within a pinned software environment, and bitwise-identical results may not be feasible across hardware or numerical libraries. In that case, state the scope of reproducibility and a defensible numerical tolerance. Preserving inputs and outputs makes later review possible, but it does not prove that the distributions, model structure, or scientific assumptions are appropriate. Reproducibility supports verification; it is not a substitute for scientific evaluation.
Check Your Understanding
Answer each question, then compare with the short solutions.
- Why is recording a seed alone insufficient to guarantee the same sampled PBPK inputs?
- Name two fields that should be saved for each simulation in addition to its output.
- A nested run has parameter draw 3 and individual draw 8. Why might preserving both identifiers matter?
- Two reruns have the same output values but those values are assigned to different simulation IDs. Is the input-output record reproduced?
- What should an analyst specify when bitwise-identical model outputs cannot be guaranteed across computing environments?
Solutions. (1) The generator, implementation, sampling procedure, transformations, call order, and stream assignment can also change the draws. (2) For example, the sampled input values and a stable simulation ID; status and units are also important. (3) Both identifiers show which parameter and individual draws produced the result and preserve the nested design. (4) No. The values must be matched to the correct inputs and IDs, not merely compared as an unlabelled set. (5) State the scope of reproducibility and a justified numerical tolerance for comparing outputs.