Tutorials › Pharmacometric Workflows and Reproducible Modeling › Validating Simulated Outputs With Invariants

Simulation Pipelines · Tutorial 271 of 399

Validating Simulated Outputs With Invariants

Learn to turn model and data expectations into explicit checks that detect invalid simulation outputs without confusing numerical tolerance with scientific validity.

Advanced 12 min read

What You'll Learn

  • Define structural invariants for simulated profiles and derived metrics.
  • Choose and interpret absolute-plus-relative tolerances.
  • Check concentration bounds and dose-event jumps with worked calculations.
  • Verify a one-compartment mass-balance relationship.
  • Diagnose failed checks without silently altering simulation results.
  • Distinguish structural validation from comparison with reference outputs.

Why Validate Simulated Outputs With Invariants?

A simulation can complete without an error message and still produce outputs that violate basic expectations. A concentration profile might contain an impossible negative value, a dose event might be missing from the event record, or a derived exposure metric might be infinite. These checks do not establish that a model is scientifically appropriate. They establish whether the generated output is structurally consistent with specified assumptions and the simulation setup.

An invariant is an expectation that should hold for every valid output, or within a declared numerical tolerance. Some invariants come from the quantity being modeled, such as nonnegative concentration. Others follow from the dosing history or model equations, such as a bolus dose producing a specified change in amount. Still others are conditional: a concentration should decline during an interval with elimination and no input, but need not decline while a dose or infusion is active.

Definition: A validation invariant is a precisely stated property that a simulated output is required to satisfy, given the model, dosing history, units, and numerical method. A check passes when the property holds exactly or within a specified tolerance. A failed check identifies an inconsistency to investigate; by itself, it does not identify the cause.

A useful invariant states what is checked, the relevant output and units, the conditions under which the expectation applies, and how much numerical discrepancy is allowed. “The profile looks reasonable” is not reproducible. “For each subject and scheduled time, concentration is at least minus the declared tolerance in mg/L” is testable.

Exact Conditions and Numerical Tolerances

Some properties are exact by construction. For example, a simulation record may be required to contain one event for each scheduled dose. Other properties involve floating-point arithmetic or numerical integration. A solver may return a tiny negative value when the mathematical solution is nonnegative. Such a value should be assessed against a prespecified tolerance rather than judged by visual inspection.

Definition: For a computed quantity \(x\) and a reference value \(r\), an absolute-plus-relative tolerance check passes when \( |x-r| \leq \mathrm{atol}+\mathrm{rtol}|r| \). Here \(\mathrm{atol}\) is an absolute tolerance in the units of \(x\), and \(\mathrm{rtol}\) is a dimensionless relative tolerance. The values must be chosen for the quantity, scale, and numerical method before reviewing the output.

The absolute term matters near zero, where a relative difference can be unstable or undefined. The relative term permits proportionate error for large reference values. Tolerances should reflect the precision and accuracy required for the task; they are not adjustable allowances for an implausible result. For a lower-bound check \(x \geq L\), a common rule is to pass if \(x \geq L-\mathrm{tol}\), where the tolerance is stated in the units of \(x\).

A short justification explains why the combined tolerance is useful. If the expected value is \(r\), the absolute discrepancy is \(|x-r|\). Requiring it to be no more than \(\mathrm{atol}+\mathrm{rtol}|r|\) means that near zero the permitted discrepancy is approximately \(\mathrm{atol}\), while at large \(|r|\) the relative contribution grows with the scale. This does not derive a universal tolerance: numerical solvers, output precision, and the scientific purpose determine appropriate values.

For a one-sided bound, retain the sign of the difference when diagnosing a failure. If \(L=-0.001\) mg/L and \(x=-0.015\) mg/L, then \(x-L=-0.015-(-0.001)=-0.014\) mg/L. The signed difference is \(-0.014\) mg/L; its magnitude is \(0.014\) mg/L. Keeping these distinct makes the direction and size of the violation clear.

Worked Example: Checking a Concentration Bound

Suppose a hypothetical simulation produces concentrations in mg/L, for which the stated lower boundary is \(L=-0.001\) mg/L and the allowed discrepancy from that boundary is \(0.00050005\) mg/L. One simulated value is \(x=-0.015\) mg/L. Determine whether it passes.

The signed difference from the boundary is \(x-L=-0.015-(-0.001)=-0.014\) mg/L. Its absolute magnitude is \(|x-L|=0.014\) mg/L.

Because \(0.014\) mg/L is greater than the declared tolerance \(0.00050005\) mg/L, this value fails the bound check. The negative signed difference shows that it is below the boundary; it is not merely a tolerated small numerical undershoot. The value and its subject and time should be retained in the diagnostic record for investigation.

Checks for Profiles, Events, and Derived Metrics

An invariant should be tied to the specific structure of the simulation. For concentrations, a nonnegative lower bound may be appropriate when the model represents a physical concentration. A tiny negative solver artifact may be tolerated if the rule permits it, but a large negative value should be flagged. Do not automatically replace every negative value with zero: clipping can conceal solver instability, a unit conversion error, or a bug in the model or output transformation.

Dose-event checks link the dosing history to the simulated state. For a one-compartment model with an instantaneous intravenous bolus of dose \(D\) and volume \(V\), concentration immediately after the event should exceed concentration immediately before it by \(D/V\), provided there are no simultaneous state changes. The event time, dose, units, and pre- and post-event values must be aligned. A check at a nearby time instead of the event time can be misleading, because elimination or other processes may also change the state.

Worked Example: Verifying an Intravenous Bolus Event

In a hypothetical one-compartment simulation, a \(100\) mg bolus is given at 8 h and the apparent volume is \(20\) L. The event record contains one bolus at 8 h. The model output immediately before the event is \(1.50\) mg/L and immediately after is \(6.495\) mg/L. Use a tolerance of \(0.010\) mg/L for the jump check.

The expected concentration increment is \(D/V=100\text{ mg}/20\text{ L}=5.00\) mg/L. Therefore the expected post-event concentration is \(1.50+5.00=6.50\) mg/L.

The observed post-event value differs from expectation by \(6.495-6.50=-0.005\) mg/L, with absolute discrepancy \(0.005\) mg/L. Since \(0.005 \leq 0.010\), the jump check passes. The event count and timestamp also match the dosing record in this example. If the event were absent or recorded at a different time, the numerical jump alone would not establish correct dose-event alignment.

Derived metrics need their own checks. A metric such as area under the concentration-time curve should be finite, expressed in the expected units, and compatible with its input profile. Finiteness is necessary but not sufficient: a finite value can still be wrong because times were misordered or units were inconsistent. Where the metric is calculated from sampled points, checks can also verify that the time sequence is ordered and that the calculation used the intended interval.

Worked Example: Checking a Finite Exposure Metric

A hypothetical profile has concentrations \(4\), \(2\), and \(1\) mg/L at times \(0\), \(1\), and \(2\) h. Calculate the trapezoidal area from 0 to 2 h and state two structural checks.

For adjacent points, the trapezoidal area is the time difference multiplied by the average of the endpoint concentrations. Thus, \( \mathrm{AUC}=1(4+2)/2+1(2+1)/2=3+1.5=4.5 \) mg·h/L.

The value \(4.5\) mg·h/L is finite and has the expected concentration-times-time units. The times are ordered \(0<1<2\) h, and all three concentrations are nonnegative. These checks do not establish that the profile is scientifically plausible or that the selected interval is appropriate; they confirm only the stated structural properties.

Model-Specific Conservation and Monotonicity

Conservation checks compare changes in a system state with known inputs and outputs. Monotonicity checks assert that a quantity moves in one direction under conditions that make that direction necessary. Neither should be applied as a generic rule to every profile. For example, concentration may rise during infusion or after a dose, and a conserved quantity may not exist in a model with clearance or degradation unless those losses are included in the accounting.

For a one-compartment amount \(A(t)\) with an instantaneous dose \(D\) and first-order elimination rate constant \(k\), the amount after the dose is \(A(0)=D\), and during a no-input interval it follows \(A(t)=D e^{-kt}\). The eliminated amount by time \(t\) is \(D-A(t)\). Adding the amount remaining to the amount eliminated gives \(A(t)+[D-A(t)]=D\), so the dose is accounted for. This identity assumes the specified one-compartment model, no other input or loss pathway, and consistent units.

Worked Example: Checking One-Compartment Mass Balance

Consider a hypothetical \(100\) mg bolus in a one-compartment model with \(k=0.100\) h\(^{-1}\). At 10 h, the simulated amount is \(36.8\) mg and the model-derived cumulative eliminated amount is \(63.2\) mg. Check the mass balance and the expected direction of change.

The analytical amount is \(100e^{-0.100(10)}=100e^{-1}\approx36.79\) mg, consistent with the simulated \(36.8\) mg to the reported precision. The balance is \(36.8+63.2=100.0\) mg, matching the 100 mg dose. The balance residual, defined here as remaining plus eliminated minus dose, is \(100.0-100=0.0\) mg.

With no input after the bolus and \(k>0\), \(e^{-kt}\) decreases as time increases, so amount and concentration should decline. For instance, the exact amount at 5 h is \(100e^{-0.5}\approx60.65\) mg, greater than the amount at 10 h, approximately \(36.79\) mg. A nondecreasing profile over this no-input interval would warrant investigation. This check would not be valid across a subsequent dose.

A Practical Validation Sequence

Validation is most useful when checks are explicit, localizable, and run before summaries obscure individual failures. A useful sequence is:

1
State the conditions.
Record which model, regimen, time interval, units, and output variable the invariant applies to.
2
Define the calculation and tolerance.
Specify the reference value or bound, the residual or discrepancy, and any absolute or relative tolerance before inspecting results.
3
Run checks at the replicate level.
Retain the replicate, subject, time, event, observed value, expected value, and failure magnitude needed to locate a problem.
4
Investigate failures before summarizing.
Check event alignment, units, state transitions, solver behavior, and metric inputs; do not silently discard or repair a failure.

Finite-value checks are also essential before calculating summaries, but they should not be confused with bounds or model-specific invariants. The earlier material on failed and nonfinite simulation replicates addresses how to handle failed replicates; here the focus is on checking whether outputs that can be evaluated obey structural expectations. A failed invariant should be reported with its condition and tolerance, not hidden by a population average.

Diagnostics, Alternatives, and Limitations

Repeated failures at the same dose time suggest checking event alignment or dose units. Failures that appear only at later times may point to solver instability, extrapolation, or a transformation error. A sudden negative concentration isolated to one replicate can indicate a replicate-specific parameter or state problem. A mass-balance discrepancy that grows over time may reflect a missing term in the accounting, such as an input or loss pathway. These patterns help prioritize investigation but do not prove a particular cause.

ApproachWhat it can establishMain limitation
Structural invariant checksWhether outputs satisfy declared constraints or identitiesOnly checks properties that have been specified correctly
Comparison with a reference outputWhether a run agrees with a known result under matching settingsA shared error in the reference and new run can go undetected
Visual profile reviewCan reveal broad patterns, discontinuities, or suspicious shapesSmall or localized violations can be missed; judgments may vary

These approaches complement one another. An invariant can pass even when the model implementation is wrong in a way that preserves the invariant, and a reference comparison can pass if both runs share the same defect. Conversely, a reference comparison may fail because settings differ even though each output obeys its own structural checks.

The most common mistake is treating an invariant as a universal biological law. “Concentration must decrease with time” is wrong when doses or infusion inputs occur. The correction is to state the condition: for a specified no-input interval in a model with positive elimination and no other process that increases concentration, the expected direction is decreasing. Another mistake is to call any tiny violation acceptable after seeing it. Choose the tolerance in advance, report the residual, and investigate values outside it. Tolerance is a numerical rule, not a declaration that a model is valid.

Key takeaway: A useful invariant is conditional, measurable, and documented with units and a prespecified tolerance. It is a targeted test of structural consistency, not a substitute for scientific model assessment.

Check Your Understanding

Answer each question, then compare with the solutions.

  1. A lower bound is \(0\) mg/L, the tolerance is \(0.002\) mg/L, and a simulated concentration is \(-0.001\) mg/L. Does it pass the stated bound check? What does that result mean?
  2. A 50 mg bolus is given into a one-compartment model with volume 10 L. If the pre-dose concentration is 2 mg/L, what post-dose concentration is expected, assuming no simultaneous state change?
  3. Why should a monotonic-decline check be restricted to a no-input interval?
  4. A mass-balance check has a residual of \(0.20\) mg when its declared tolerance is \(0.05\) mg. What should be reported, and what should not be done automatically?

Solutions. (1) It passes because \(-0.001 \geq -0.002\) mg/L. This means it meets the declared numerical bound, not that the underlying model is scientifically valid. (2) The increment is \(50/10=5\) mg/L, so the expected post-dose concentration is \(2+5=7\) mg/L. (3) Input can increase the state, so decline is not required while a dose or infusion is active. (4) Report the \(0.20\) mg residual and the \(0.05\) mg tolerance as a failure, with enough context to locate it. Do not silently alter the state or remove the replicate without investigation.

Clinical Biostats