Why Invariants Are Useful Verification Tests
A mechanistic model can produce smooth, plausible-looking trajectories even when its implementation contains an error. A reaction term may be missing from one state equation; a transfer may be added to one compartment but not removed from another; or a bound species may be counted incorrectly in a total. These defects can be difficult to notice from plots alone. Conservation and invariant tests turn known structural relationships into checks on the executable model.
The tests are verification tools: they ask whether the code respects properties implied by the equations it is supposed to implement. They do not establish that the equations adequately represent biology or that a model is suitable for a particular prediction claim. The earlier material on analytical reference solutions provides another way to check an implementation when a tractable reference solution is available. Invariant tests are especially useful when the full trajectory has no convenient closed-form solution.
From Model Equations to an Invariant
Let \(x(t)\) be a vector of model states, such as amounts of chemical species or amounts in compartments. For a reaction or transfer system, write the state equations as
Here \(S\) is the stoichiometric matrix: each column describes how one reaction or transfer changes the states. The vector \(v(x,t)\) contains the corresponding rates. A row entry of \(S\) is negative when that process consumes a state and positive when it produces one.
The matrix condition gives a short derivation. For a fixed vector \(c\),
Therefore \(c^\mathsf{T}x(t)=c^\mathsf{T}x(0)\), provided the equations have the stated form and no unaccounted external source or sink changes the total. The result does not depend on the particular rate law: the rates may be nonlinear or time-varying, while the stoichiometric cancellation still holds.
A conserved total is only one kind of invariant. For a nonlinear candidate \(g(x)\), the corresponding condition is \(\nabla g(x)^\mathsf{T}f(x,t)=0\), where \(f(x,t)=dx/dt\), throughout the relevant state domain. An expected monotonic quantity is different: its derivative has a known sign, such as \(dg/dt\leq 0\), rather than being identically zero. Such sign checks can reveal errors that preserve a total.
Worked Tests for Binding and Conversion
Worked Example: Checking a Reversible Binding Model
Consider hypothetical amounts \(A\), \(B\), and \(C\), where \(A+B\rightleftharpoons C\). The forward rate is \(k_{\mathrm{on}}AB\), and the reverse rate is \(k_{\mathrm{off}}C\). Assume there are no other reactions, sources, or sinks, and all three states are expressed in the same amount units. The equations are
Binding consumes one unit of free \(A\) and one of free \(B\) while producing one unit of \(C\). Consequently \(A+C\) and \(B+C\) should each be conserved. For example,
Suppose the hypothetical initial amounts are \(A=10\), \(B=8\), and \(C=0\) nmol. At a saved simulation time, the executable returns \(A=4\), \(B=2\), and \(C=6\) nmol. The initial totals are \(A+C=10+0=10\) nmol and \(B+C=8+0=8\) nmol. The returned totals are \(4+6=10\) nmol and \(2+6=8\) nmol, so both checks pass at this time.
If instead the executable returns \(B=1.7\) nmol with the other values unchanged, then \(B+C=1.7+6=7.7\) nmol. The discrepancy is \(7.7-8=-0.3\) nmol. That failure suggests a defect in the implemented bookkeeping or a missing process, but it does not by itself identify the faulty line of code. A further check can evaluate the equations’ implied derivatives, not just their saved output values.
Worked Example: Combining Conservation With Monotonicity
Suppose a hypothetical closed system converts \(X\) irreversibly to \(Y\) at rate \(kX\), with \(k\geq 0\). There are no other reactions or external fluxes. Then
Adding the equations gives \(d(X+Y)/dt=0\), so \(X+Y\) is conserved. Under the stated condition \(X\geq 0\), the directional checks are \(dX/dt=-kX\leq 0\) and \(dY/dt=kX\geq 0\). Thus \(X\) should not increase and \(Y\) should not decrease along the trajectory.
For hypothetical outputs in mg, let the initial state be \(X=12\), \(Y=0\). At time 1, suppose \(X=8\), \(Y=4\); at time 2, suppose \(X=5.5\), \(Y=6.5\). The totals are \(8+4=12\) mg and \(5.5+6.5=12\) mg, so conservation holds at both times. But \(X\) increased from 5.0 to 5.5 mg if an additional saved output at time 1.5 reports \(X=5.0\) mg. That increase violates the expected monotonicity, despite the conserved total. If the same output reports \(Y=7.0\) mg, the total becomes \(5.5+7.0=12.5\) mg as well, giving a separate conservation failure.
The monotonicity test is conditional on the specified model and domain. A source of \(X\), a reverse conversion, or a negative rate would invalidate the stated directional expectation. Before treating an increase as an implementation error, check that the executable actually represents the closed, irreversible conversion model assumed here.
Balances When the System Has External Fluxes
A total need not remain constant when the model includes dosing, elimination, degradation, or another external flux. In that case, test a balance that accounts for the flux rather than demanding conservation of the internal states alone. For example, transfer between two compartments cancels from their sum, while input and elimination change it in known directions.
Worked Example: Auditing Transfer, Input, and Elimination
Consider hypothetical drug amounts \(A\) in a central compartment and \(B\) in a peripheral compartment, both in mg. Drug transfers from \(A\) to \(B\) at rate \(k_{12}A\), returns at rate \(k_{21}B\), enters through an input rate \(u(t)\), and is eliminated from \(A\) at rate \(k_eA\). Assume cumulative input \(D_{\mathrm{in}}\) and cumulative elimination \(E\) are tracked in mg, with derivatives \(dD_{\mathrm{in}}/dt=u(t)\) and \(dE/dt=k_eA\). The equations for the compartments are
When the compartment equations are added, the internal transfers cancel:
Therefore the balance quantity \(A+B-D_{\mathrm{in}}+E\) is invariant. Indeed, its derivative is \(u(t)-k_eA-u(t)+k_eA=0\).
Suppose the initial amounts and cumulative quantities are all zero. At a hypothetical time, let \(A=5.1\) mg, \(B=2.5\) mg, \(D_{\mathrm{in}}=10.0\) mg, and \(E=2.4\) mg. The balance is \(5.1+2.5-10.0+2.4=0.0\) mg, as expected. At a later time, suppose \(A=3.4\) mg, \(B=2.5\) mg, \(D_{\mathrm{in}}=10.0\) mg, and \(E=4.1\) mg. The balance is \(3.4+2.5-10.0+4.1=0.0\) mg again.
If the later peripheral amount were instead \(2.8\) mg, the balance would be \(3.4+2.8-10.0+4.1=0.3\) mg. The nonzero residual flags inconsistency with the stated input, elimination, and transfer bookkeeping. The residual alone cannot distinguish a wrong transfer term from a wrong cumulative-elimination calculation, so the component equations and recorded fluxes should then be inspected.
Implementing the Checks and Diagnosing Failures
A useful test starts from the model’s declared states, reactions, and external fluxes—not from a pattern noticed in one output plot. Derive the expected total or directional relationship, calculate it at the initial state and at multiple saved times, and record the residual. For a conserved quantity \(I(t)\), a simple residual is \(r(t)=I(t)-I(0)\). For a balance with known inputs and outputs, construct the full balance quantity before calculating its residual.
Floating-point arithmetic and numerical integration mean that a computed residual need not be exactly zero. Set a tolerance appropriate to the solver, state magnitudes, and units. One practical approach is to compare the absolute residual with an absolute tolerance plus a relative tolerance times a representative scale for the conserved total. Report both the residual and the scale; a residual of \(0.01\) mg has a different interpretation for a total of \(0.1\) mg than for a total of \(10^6\) mg. A tolerance should be tight enough to catch meaningful implementation defects and should not be enlarged merely to make a failing test pass.
Check directional conditions over the trajectory as well as at selected endpoints. Sparse output can miss a brief reversal, and a monotonicity check between saved points does not prove the condition held at every internal solver step. Conversely, tiny numerical fluctuations near a boundary can create an apparent sign violation. Diagnose such a result by examining the scale and duration of the excursion, solver behavior, and the actual equations; do not silently discard it.
Several complementary approaches serve different purposes. An analytical reference solution can compare the computed trajectory with a known answer when one is available. Invariant and balance tests can apply even when no complete reference trajectory is available, because they check structural relationships. Visual inspection may reveal grossly implausible behavior, but it is less systematic and can miss compensating errors. These approaches complement one another; none alone establishes that the implementation is correct.
| Check | What it tests | Important limitation |
|---|---|---|
| Conserved total or balance | Whether specified contributions and fluxes add up consistently | Different errors can cancel and leave the total unchanged |
| Monotonicity | Whether a quantity moves only in its model-implied direction | Valid only under the assumptions that imply that direction |
| Analytical reference comparison | Agreement with a known solution for a selected case | A reference may not be available for the full model |
Common Pitfalls and Limits
- Testing the wrong total. For binding, free and bound forms both contribute to the relevant moiety total. Derive the weights from the reaction bookkeeping rather than assuming every state contributes equally.
- Ignoring external fluxes. A changing amount is not a conservation failure if input or elimination is present and the balance accounts for it.
- Checking only one quantity. A total can remain correct while individual states violate their expected directions. Use both conservation and monotonicity when each follows from the equations.
- Checking only a final output. Errors can appear temporarily and later cancel. Evaluate multiple times and, where practical, inspect derivatives or solver outputs within the interval.
- Treating a pass as proof of correctness. A test can miss errors that preserve the same invariant, and an incorrectly specified invariant can produce misleading results. Independent checks and careful model-to-code traceability remain necessary.
These tests are strongest when the claimed property is derived directly from the intended equations and implemented in a check that is as independent as practical from the code being tested. A pass supports a narrow conclusion: the tested trajectories satisfy the checked relationships within the selected tolerance. It does not show that every equation, parameter, or mechanistic assumption is correct.
Check Your Understanding
Work through each question, then compare your reasoning with the solutions below.
- For \(P\rightleftharpoons Q\) in a closed system, state a conserved total and justify it from the state equations.
- A model of \(M\to N\) reports \(M+N=20\) at every saved time, but \(M\) rises between two times. What additional property may have failed, and under what assumptions?
- In a two-compartment model with internal transfer and elimination from the first compartment, why does the sum of compartment amounts generally fail to be invariant?
- A balance residual is \(0.02\) mg for a total near \(0.05\) mg. Why is the residual’s absolute value alone insufficient to declare the test passed?
Solutions. 1. If \(dP/dt=-r_f+r_r\) and \(dQ/dt=r_f-r_r\), then \(d(P+Q)/dt=0\), so \(P+Q\) is conserved. 2. The monotonicity condition for \(M\) may have failed: for an irreversible conversion with no source and nonnegative rate, \(dM/dt\leq0\), so \(M\) should not rise. A source or reverse reaction would change that conclusion. 3. Internal transfer cancels, but elimination removes material from the two-compartment system, so the sum changes unless elimination is included in an augmented balance. 4. The residual is large relative to the total: \(0.02/0.05=0.4\), or 40%. Whether it passes depends on a justified tolerance and scale, not the absolute residual alone.