Why a simulated balance needs its own check
A numerical trajectory can look smooth and biologically plausible while failing to preserve a quantity that the model says is conserved. Conversely, a small change in a conserved total does not automatically mean that the model is wrong: it may reflect numerical error, floating-point rounding, interpolation at reported times, or an inconsistency between the balance being checked and the quantities included in it. A conservation residual makes the discrepancy explicit and gives it a scale and a sign.
The earlier material on conserved moieties established how a closed reaction system can have constant linear combinations of its state amounts. Here the focus is the residual calculated from a simulated trajectory, and how that residual behaves when solver settings change. The residual is a diagnostic, not a replacement for checking that the model equations themselves express the intended balance.
For several conserved quantities, collect their vectors as columns of a matrix \(L\). The residual vector is \(r(t)=L^\mathsf{T}x(t)-c\), where \(c=L^\mathsf{T}x(0)\). Each component should be inspected; a small aggregate score can conceal a larger error in one balance. Use amounts in the balance. If the model state is a concentration in a changing-volume compartment, convert it to amount using the relevant volume before applying an amount balance.
What the residual says about the equations and the solver
For a closed reaction system, write the differential equation as \(\dot{x}=Sv(x)\), where \(S\) is the stoichiometric matrix and \(v(x)\) is the vector of reaction rates. A conserved vector satisfies \(l^\mathsf{T}S=0\). Therefore $$ \frac{d}{dt}\left(l^\mathsf{T}x(t)\right) =l^\mathsf{T}\dot{x}(t) =l^\mathsf{T}Sv(x(t)) =0. $$ With a valid initial condition, the exact solution has \(l^\mathsf{T}x(t)=c\) at every time. This is the balance the numerical result is being tested against.
The distinction between equation-level conservation and trajectory-level residual is essential. If the implemented right-hand side satisfies \(l^\mathsf{T}f(x)=0\) for every relevant state \(x\), the continuous model conserves the quantity. A numerical solver then approximates that model. Depending on the integration method and implementation, its reported states may still show small discrepancies from roundoff, approximate internal solves, or interpolation. Some common Runge–Kutta updates preserve linear invariants particularly well when their stages and arithmetic satisfy the balance, so a large residual is not an inevitable consequence of using a numerical solver.
If instead \(l^\mathsf{T}f(x)\ne 0\), the continuous equations themselves change the quantity. Tightening a solver tolerance cannot repair that structural balance violation. For an open system, the quantity need not be constant: compare its change with the specified cumulative inputs and outputs.
The integral may itself be evaluated numerically. Thus, a nonzero open-system residual can include error from both state integration and flux integration. Report how the flux integral was obtained when the expected balance is small enough that this numerical error matters.
Worked Example: Closed-System Residual at Reported Times
Suppose a hypothetical closed binding system has two state amounts, free ligand \(L\) and bound ligand \(RL\), both measured in nmol. The conserved ligand total is \(L+RL\), and the initial state is \(L(0)=9.00\) nmol and \(RL(0)=3.00\) nmol. At two later output times, a solver reports the following values.
| Time (h) | \(L\) (nmol) | \(RL\) (nmol) | \(L+RL\) (nmol) |
|---|---|---|---|
| 0 | 9.00 | 3.00 | 12.00 |
| 1 | 7.42 | 4.59 | 12.01 |
| 2 | 6.18 | 5.83 | 12.01 |
The initial conserved amount is \(c=9.00+3.00=12.00\) nmol. At 1 h, $$ r(1)=7.42+4.59-12.00=0.01\ \text{nmol}. $$ At 2 h, \(r(2)=6.18+5.83-12.00=0.01\) nmol. Using the initial total as the scale, the relative residual at each later time is \(0.01/12.00=0.0008333\), or about \(0.0833\%\). The residual is positive: the reported states contain that much more total ligand than the initial balance.
This calculation establishes the size of the observed discrepancy, but not its cause. To assess numerical drift, repeat the same model and initial state with tighter error tolerances and/or a smaller maximum step, then compare residuals on the same time grid. If the residual falls substantially and the trajectory otherwise converges, that supports a numerical-error explanation. It is evidence, not proof that the equations and implementation are correct.
Interpreting solver refinements
A solver’s relative and absolute error tolerances govern its error-control procedure; they are not direct upper bounds on a conservation residual. In particular, state variables with different magnitudes may receive different effective error control, and the conserved quantity may combine several state values. A residual can therefore be larger or smaller than either tolerance. Compare residuals empirically across runs rather than declaring a balance valid solely because it is below the tolerance setting.
A useful comparison holds the model, parameters, initial conditions, and output times fixed while changing one numerical control at a time. For example, first tighten tolerances with the same maximum step, then reduce the maximum step while keeping the tolerances fixed. Record the maximum absolute residual over the time window and, where useful, the maximum scaled residual. If both settings change at once, a reduced residual still indicates improvement, but it is harder to identify which change mattered.
Residuals need not decrease monotonically at every individual output time as settings change. Adaptive solvers can choose different sequences of steps; rounding and interpolation can also affect reported values. The stronger question is whether the trajectory and balance residual stabilize under reasonable refinement. A persistent, non-convergent, or unexpectedly large residual calls for checking both the balance definition and the model’s right-hand side.
Worked Example: Comparing Residuals Across Solver Settings
For the same hypothetical closed system above, suppose two runs use the same initial state and output interval. The solver reports these largest absolute residuals over 0–2 h:
| Run | Relative tolerance | Maximum step (h) | Maximum \(|r|\) (nmol) |
|---|---|---|---|
| A | \(10^{-4}\) | 0.50 | 0.0120 |
| B | \(10^{-7}\) | 0.10 | 0.00018 |
For a reference total of 12.00 nmol, Run A has maximum scaled residual \(0.0120/12.00=0.00100\), or \(0.100\%\). Run B has \(0.00018/12.00=0.0000150\), or \(0.00150\%\). The absolute residual fell by a factor of \(0.0120/0.00018=66.7\). In this hypothetical comparison, the much smaller residual under the refined settings is consistent with solver-related drift.
The conclusion is limited: both tolerance and maximum step changed, so this comparison does not isolate their individual effects. Nor does the reduction prove that every state or every balance is correct. A follow-up can change one setting at a time and inspect residuals for each conserved quantity. The tolerance values themselves are not the measured conservation error; the residual values are.
Checking whether the equations conserve the quantity
A direct structural test evaluates \(l^\mathsf{T}f(x)\), where \(f(x)\) is the implemented right-hand side. If the intended quantity is conserved, this expression should be zero for all relevant states, not merely at the initial condition. For a closed reaction network \(f(x)=Sv(x)\), the check is \(l^\mathsf{T}S=0\). For an open model, the right-hand side should instead match the intended external balance.
This test separates a property of the equations from an error accumulated by integration. It does not by itself establish that the chosen balance represents the intended biology; that requires a correct definition of the conserved quantity and system boundary. But if the derivative of the quantity is nonzero when it should be zero, a solver refinement is not the remedy.
Worked Example: A Nonzero Derivative Balance
Consider three amounts \(A,B,C\) in nmol and two reaction rates \(v_1,v_2\) in nmol/h. The intended reactions are \(A\to B\) at rate \(v_1\) and \(B\to C\) at rate \(v_2\), so the intended total \(A+B+C\) is conserved. Suppose an implemented right-hand side is $$ \dot{A}=-v_1,\qquad \dot{B}=v_1-v_2,\qquad \dot{C}=0. $$ Then $$ \frac{d}{dt}(A+B+C)=-v_1+(v_1-v_2)+0=-v_2. $$
If \(v_2=0.20\) nmol/h over a 2 h interval, the implemented equations reduce the total by \(0.20\times2=0.40\) nmol. From an initial total of 10.0 nmol, the equations predict 9.60 nmol at 2 h. The expected conservation residual is therefore \(-0.40\) nmol, not a small integration discrepancy. A numerical solver may accurately integrate this nonconserving right-hand side, and tighter tolerances would make the simulated result approach that balance violation more closely.
This derivative check identifies a structural inconsistency in the equations relative to the intended total. The next task is to inspect the model definition and implementation to locate its cause; the residual alone does not identify which state or term should be changed.
Diagnostics, alternatives, and limitations
A practical diagnostic sequence is to define the balance and its units first, compute the residual at all relevant output times, and then compare the residual with a clearly stated scale. Next, test whether the right-hand side satisfies the expected derivative balance. If it does, perform solver refinements and examine whether residuals and trajectories stabilize. If it does not, investigate the equations before interpreting solver settings.
Several approaches answer related but different questions. Inspecting each state trajectory can reveal implausible values but may not reveal a small balance discrepancy. Checking \(l^\mathsf{T}f(x)\) tests the equations locally but does not measure the error in a particular numerical trajectory. Refining tolerances and step sizes tests numerical sensitivity, but cannot validate the biological balance definition. Projection of simulated states onto a balance constraint can force a zero residual, but it changes the computed states and can conceal a flawed right-hand side; it should not replace diagnosis.
For relative residuals, the denominator must be chosen and reported. A scale such as the initial conserved amount is convenient when it is nonzero and representative. When that amount is zero or very small, use a meaningful positive amount scale, such as a stated characteristic amount, and report the absolute residual as well. In a system with several balances, scale each component separately rather than combining quantities with incompatible units.
Residuals are diagnostics, not universal pass/fail thresholds. Their acceptable size depends on the model’s scales, the required accuracy, and how the simulation will be used. A small residual does not prove that the model is biologically correct, and a large residual should not be attributed to the solver without checking the equations and balance definition.
Check Your Understanding
Answer each question using the balance definition and units, then compare with the solutions.
- A closed system begins with \(A=4\) nmol and \(B=6\) nmol. A simulated state has \(A=2.5\) nmol and \(B=7.6\) nmol. Calculate the absolute and signed conservation residual.
- An open system has initial amount 8 nmol, cumulative input 3 nmol, cumulative output 1.5 nmol, and simulated final amount 9.4 nmol. Calculate the open-system residual.
- For an intended closed balance \(M=A+B\), the implemented equations give \(\dot{A}=-kA\) and \(\dot{B}=kA\). Is \(M\) structurally conserved? Show the derivative check.
- A trajectory has a persistent residual after tolerances are tightened, and the implemented equations give \(dM/dt\ne0\) for a quantity intended to be conserved. Is further tolerance tightening the primary correction?
Solutions. 1. The initial total is 10 nmol and the simulated total is \(2.5+7.6=10.1\) nmol. Thus \(r=10.1-10=0.1\) nmol and \(|r|=0.1\) nmol. 2. The expected final amount is \(8+3-1.5=9.5\) nmol, so \(r_{\mathrm{open}}=9.4-9.5=-0.1\) nmol. 3. \(dM/dt=\dot{A}+\dot{B}=-kA+kA=0\), so the equations structurally conserve \(M\). 4. No. The nonzero derivative indicates that the equations do not satisfy the intended balance; investigate that structural discrepancy before treating solver refinement as the solution.