Why Test Relationships Between Runs?
A complex QSP model may have no exact solution against which to compare every simulated trajectory. That does not mean its implementation cannot be checked. Often, the model’s structure implies that two carefully related simulation runs must also have related outputs. A test of that implication is a metamorphic test.
For example, a linear drug-disposition subsystem may imply that doubling the dose doubles the drug-induced concentration at each time point. A change from one unit of a state variable to another should change its numerical values by a conversion factor, while leaving the represented physical trajectory unchanged. Neither test requires knowing the full correct trajectory in advance.
Here, an input can include the dose, regimen, parameters, initial conditions, units, or virtual individual. An output can be a state trajectory, a derived quantity, or a selected endpoint. The essential requirement is that the relation follows from stated model structure and assumptions—not merely from an expectation that the results “look reasonable.”
Specify the Relation Before Running the Test
A useful test has four explicit parts: the model conditions under which the relation should hold, the transformation applied, the output expected from that transformation, and the numerical tolerance. For a deterministic simulator, write the intended relationship as \(F(Tx)=M(F(x))\), where \(M\) describes how the output should change. For a numerical solver, replace exact equality by a tolerance rule.
The index \(i\) can identify an output, time point, or state component. The absolute tolerance \(\mathrm{atol}\) protects comparisons near zero; the relative tolerance \(\mathrm{rtol}\) scales the allowance to the expected magnitude. Tolerances should reflect numerical precision and the test’s purpose. A tolerance so wide that a meaningful implementation error passes makes the test uninformative.
These are implementation checks, not statistical hypothesis tests: there is no \(H_0\), \(H_a\), p-value, or inference about a population. For stochastic simulators, run-to-run variability requires additional treatment; the examples here concern deterministic relationships.
Dose Proportionality in a Linear Subsystem
Dose proportionality is a valuable test, but it is conditional—not a universal property of QSP models. Suppose the drug amount or concentration states \(x(t)\) follow a linear system with dose-proportional input \(u_D(t)=D u_1(t)\), and the initial drug-related state is zero. For fixed parameters, write
If the dose is multiplied by a positive factor \(c\), the input becomes \(c u_D(t)\). The state \(c x(t)\) then satisfies the transformed model, because its derivative is \(c A x(t)+c B u_D(t)=A[cx(t)]+B[c u_D(t)]\), and its initial state is still zero. Uniqueness of the solution to this linear initial-value problem therefore gives \(x_{cD}(t)=c x_D(t)\). A linear output \(z(t)=C x(t)\) has the same scaling relation.
This reasoning requires linear dynamics, unchanged parameters and regimen shape, and initial conditions that scale appropriately. A nonzero baseline should be removed before testing proportionality if the claim concerns the drug-induced increment. Saturable clearance, nonlinear binding, feedback, or nonlinear downstream pharmacology can break proportionality. In particular, a proportional concentration relation does not imply a proportional clinical response.
Worked Example: Checking Dose Proportionality
Consider a hypothetical one-compartment model with an instantaneous dose \(D\), volume \(V=5\) L, and first-order elimination rate \(k=0.2\ \mathrm{h}^{-1}\). Its concentration is \(C_D(t)=(D/V)e^{-kt}\). The model assumptions are fixed \(V\) and \(k\), linear elimination, and no background drug concentration. Compare doses of 10 mg and 20 mg.
The dose ratio is \(c=20/10=2\), so the expected relation is \(C_{20}(t)=2C_{10}(t)\). At \(t=1\) h, for example,
The same calculation at selected times gives the following hypothetical simulator outputs, rounded for display. The test should use unrounded values where available.
| Time (h) | 10 mg output (mg/L) | 20 mg output (mg/L) | Twice 10 mg output (mg/L) |
|---|---|---|---|
| 0 | 2.0000 | 4.0000 | 4.0000 |
| 1 | 1.6375 | 3.2749 | 3.2750 |
| 2 | 1.3406 | 2.6813 | 2.6812 |
| 4 | 0.8987 | 1.7973 | 1.7974 |
The small displayed differences arise from rounding and are consistent with the expected relationship. For example, at 1 h the difference between the 20 mg output and twice the 10 mg output is about \(0.0001\) mg/L. Whether it passes depends on a tolerance chosen before examining the result. A large or systematic discrepancy would prompt investigation of dose handling, model equations, or solver configuration; it would not by itself identify which is wrong.
Invariance Under State-Unit Rescaling
A second useful relation checks whether changing the units used internally for a state preserves its physical meaning. Let \(x\) be a state vector and let \(y=Sx\), where \(S\) is an invertible diagonal matrix of unit-conversion factors. If the original model is \(dx/dt=f(t,x,\theta)\), then the transformed model must use
This follows from the chain rule: \(dy/dt=S(dx/dt)\), and the original state is \(x=S^{-1}y\). Parameters, initial conditions, inputs, and observation calculations must also be expressed consistently. The test compares physical quantities after converting outputs to a common unit; it does not demand identical numerical state values in different units.
Worked Example: Rescaling a Concentration State
A hypothetical model has one concentration state \(x\) stored in mg/L and governed by \(dx/dt=-0.5x\ \mathrm{h}^{-1}\), with \(x(0)=2\) mg/L. A second implementation stores the same state as \(y\) in ng/mL. Since \(1\) mg/L equals \(1000\) ng/mL, the conversion is \(y=1000x\).
The transformed initial condition is \(y(0)=2000\) ng/mL. Its derivative is \(dy/dt=1000(dx/dt)=-0.5y\), with time still measured in hours. The physical trajectory is \(x(t)=2e^{-0.5t}\) mg/L, so the expected rescaled trajectory is \(y(t)=2000e^{-0.5t}\) ng/mL.
At \(t=1\) h, \(x(1)=2e^{-0.5}=1.21306\) mg/L and the converted value is \(1000(1.21306)=1213.06\) ng/mL. At \(t=2\) h, \(x(2)=2e^{-1}=0.735759\) mg/L and the converted value is \(735.759\) ng/mL. Suppose the rescaled implementation reports 1213.1 and 735.8 ng/mL at those times. Converting those values back gives 1.2131 and 0.7358 mg/L, respectively. The differences from the first implementation are \(0.00004\) and \(0.000041\) mg/L, consistent with the displayed rounding.
A test of this kind can expose an unconverted initial condition or a rate/input term still expressed in the old units. Passing it supports consistency of the transformation; it does not establish that the underlying model is biologically correct.
Superposition as a Relationship Between Runs
Metamorphic tests need not compare only scaled versions of one input. For a linear system with zero initial state and two input components, the response to both inputs together equals the sum of the separate responses. This relation can check a complex numerical implementation even when no exact trajectory is known. It is invalid if the relevant subsystem is nonlinear, or if interactions between the inputs are part of the model.
Worked Example: Testing Superposition of Two Inputs
In a hypothetical linear subsystem, run the simulator with input A alone, input B alone, and both inputs together. At each time point, the output is a two-component vector in arbitrary effect units. The model assumptions are linear dynamics, a zero initial state, and unchanged parameters across the three runs.
| Time | A alone | B alone | Expected sum | Both inputs together |
|---|---|---|---|---|
| 1 h | (1.2000, 0.4000) | (0.3000, 0.5000) | (1.5000, 0.9000) | (1.5002, 0.8999) |
| 2 h | (0.8000, 0.7000) | (0.6000, 0.2000) | (1.4000, 0.9000) | (1.3998, 0.9001) |
At 1 h, componentwise addition gives \((1.2000+0.3000,\ 0.4000+0.5000)=(1.5000,\ 0.9000)\). The combined-run differences are \(0.0002\) and \(-0.0001\). At 2 h, the expected sum is \((1.4000,\ 0.9000)\), and the differences are \(-0.0002\) and \(0.0001\). With a predeclared absolute tolerance of \(0.001\) effect units and relative tolerance \(10^{-4}\), all four component comparisons pass: even the largest absolute difference, \(0.0002\), is less than the absolute allowance alone. The result is consistent with superposition at these outputs and times, not proof that every state or parameter is correct.
Practical Design, Diagnostics, and Limits
A reliable metamorphic test makes its assumptions inspectable. Record the input transformation, the expected output mapping, the compared outputs and times, and the tolerances. Use the same model version and settings except for the intended transformation. If testing state-unit rescaling, transform every affected quantity—not just the state label. If testing a dose relation, hold the regimen shape and relevant covariates fixed.
When a test fails, first check whether the stated relation actually applies to the model and the selected output. Then check transformations and units, initial conditions, event handling, output mapping, and numerical tolerances. A failure may indicate a coding defect, a violated assumption, or an incorrect test specification. Conversely, a pass is limited: two runs can share the same defect and still obey the tested relation.
- Check across the relevant domain. A relation that passes at one time or dose may fail elsewhere. Select several informative times, states, and input settings without turning the test into an unbounded search.
- Inspect the size and pattern of residuals. Small isolated discrepancies may be numerical; consistent scaling errors or growing discrepancies may indicate a conversion or input-handling problem.
- Do not confuse numerical tolerance with scientific uncertainty. Tolerances address computational comparison. They do not represent uncertainty in model parameters, patients, or clinical outcomes.
- Use complementary checks. Earlier material on analytical reference solutions, invariants, limiting cases, and dimensional assertions offers different checks. Those can be more direct when an exact solution or conserved quantity is available.
Metamorphic tests are especially useful when exact references are unavailable, but they do not replace comparison with analytical solutions where those exist, nor do they establish predictive validity. Their strength depends on choosing a relation that is both mathematically justified and sensitive to plausible implementation errors. They also do not localize a failure automatically: a failed relation tells the developer where to investigate, not which line of code is wrong.
Common Pitfall
Correction: Dose proportionality follows only under conditions such as linear dose-proportional dynamics and a suitably scaled initial state. Nonlinear clearance, feedback, or response functions can make the relationship false even when implemented correctly. State the conditions and derive the expected relation before treating a failure as evidence of a defect.
A related mistake is to compare raw numerical states after changing units. If one run reports mg/L and another reports ng/mL, their numbers should differ by a factor of 1000 for the same concentration. Compare after conversion to a common unit, using a tolerance appropriate to that unit and magnitude.
Check Your Understanding
Answer each question, then compare with the solutions below.
- A linear drug state starts at zero and has dose-proportional input. What output relation is expected when dose changes from \(D\) to \(3D\), and name one condition that could invalidate it.
- A state is converted from mg/L to ng/mL. What numerical conversion factor is required, and what should be compared between the two runs?
- For a superposition test, the separate outputs are 0.7 and 1.1 effect units, while the combined output is 1.8004. What is the absolute discrepancy from the expected value?
- Why does a metamorphic test pass not prove that a QSP model is biologically correct?
Solutions. (1) The expected output is multiplied by 3, provided the linear dynamics and relevant initial conditions support that scaling; saturable elimination is one condition that can invalidate it. (2) Multiply the mg/L numerical value by 1000 to express it in ng/mL, then compare values in a common unit. (3) The expected combined output is \(0.7+1.1=1.8\), so the absolute discrepancy is \(|1.8004-1.8|=0.0004\) effect units. (4) The test checks one specified relationship between runs. A model can satisfy that relationship while containing a shared implementation error or while being inconsistent with biological evidence.