Unit Tests for Equilibria and Limiting Cases

Why Test Equilibria and Limiting Cases?

An ordinary differential equation (ODE) function can run without errors and still encode the wrong biology or algebra. A sign error, a missing input term, or an incorrect state index may produce plausible trajectories under some parameter settings. Unit tests target simple situations in which the model must satisfy known properties, independently of a full simulation. Here the tests concern a one-state turnover model: its zero-input equilibrium, its behavior when a rate is zero or approaches zero, and derivatives of its steady state with respect to parameters.

The model is deliberately simple enough to derive its expected results exactly. The same testing logic can then be applied to components of larger mechanistic models. These tests check the right-hand side and steady-state calculations; they do not replace checks of initial conditions, conservation, or ODE integration.

Model, Equilibrium, and Test Tolerances

Let \(R(t)\) be a response amount or concentration, and let \(k_{\mathrm{in}}\) be a constant input or production rate. Let \(k_{\mathrm{out}}\) be a first-order elimination rate constant. The turnover model is

For example, if \(R\) is measured in mg/L and time in hours, \(k_{\mathrm{in}}\) has units mg/L/h and \(k_{\mathrm{out}}\) has units h\(^{-1}\). The right-hand side \(f\) therefore has units mg/L/h. The usual turnover assumptions here are that both rates are constant in time, the loss is proportional to the current \(R\), and the model is considered for \(k_{\mathrm{in}}\geq 0\), \(k_{\mathrm{out}}\geq 0\), and \(R\geq 0\).

A unit test compares a computed value with an expected value. Floating-point calculations should generally allow a small discrepancy rather than require exact equality. A useful rule is to pass a test when the absolute difference is no larger than an absolute tolerance plus a relative tolerance times the magnitude of the expected value:

The absolute tolerance protects comparisons near zero; the relative tolerance scales the allowance for larger values. Choose tolerances for the quantity being tested and state its units. A tolerance appropriate for \(R\) need not be appropriate for \(dR/dt\) or for a parameter derivative, since their units differ. Tolerances should be tight enough to expose meaningful errors, but not so tight that harmless roundoff causes failures.

Test the Zero-Input Equilibrium

In this model, zero input means setting \(k_{\mathrm{in}}=0\). For \(k_{\mathrm{out}}>0\), the expected equilibrium is \(R^*=0\). A direct test evaluates the right-hand side at that state and checks that it is zero within a tolerance. It can also evaluate the right-hand side at a positive state: with zero input, a positive response must decline, not grow. That second check is useful because a faulty implementation might return zero at \(R=0\) while still using the wrong sign for the loss term.

The distinction between an equilibrium test and a test at an arbitrary state matters. For \(k_{\mathrm{in}}>0\), substituting \(R=0\) gives \(f=k_{\mathrm{in}}>0\), not zero: the response rises from zero toward its positive equilibrium. Testing that \(R=0\) is an equilibrium for every parameter setting would therefore be incorrect.

Zero-Rate Limits and the Boundary at No Elimination

The equilibrium formula \(R^*=k_{\mathrm{in}}/k_{\mathrm{out}}\) requires \(k_{\mathrm{out}}>0\). It must not be evaluated by dividing by zero. At \(k_{\mathrm{out}}=0\), the model instead has \(dR/dt=k_{\mathrm{in}}\). If \(k_{\mathrm{in}}>0\), the response increases linearly and there is no equilibrium. If both rates are zero, the derivative is zero at every state, so every \(R\) is an equilibrium; there is no unique equilibrium selected by this model.

A limiting-case test asks whether the model behaves consistently as a rate approaches zero. For constant rates, the solution from an initial state \(R(0)=R_0\), when \(k_{\mathrm{out}}>0\), is

Although this expression contains division by \(k_{\mathrm{out}}\), its limit as \(k_{\mathrm{out}}\) tends to zero is finite. Since \(1-e^{-k_{\mathrm{out}}t}\) is approximately \(k_{\mathrm{out}}t\) for small \(k_{\mathrm{out}}t\), the second term tends to \(k_{\mathrm{in}}t\), and the first tends to \(R_0\). Thus the zero-elimination limiting trajectory is \(R_0+k_{\mathrm{in}}t\), in agreement with the ODE when \(k_{\mathrm{out}}=0\). This limiting result is a useful check on the model and on code that evaluates the solution formula; the subtraction in \(1-e^{-k_{\mathrm{out}}t}\) can also lose numerical precision when its argument is very small.

A complementary zero-input limit is \(k_{\mathrm{in}}\to 0\) with \(k_{\mathrm{out}}>0\). The equilibrium \(k_{\mathrm{in}}/k_{\mathrm{out}}\) then tends to zero. The trajectory tends toward \(R_0e^{-k_{\mathrm{out}}t}\), which decays to the zero-input equilibrium. Testing this limit helps expose code that mistakenly leaves a nonzero production term when the input parameter is set to zero.

Derive and Test Steady-State Derivatives

Steady-state derivatives describe how the equilibrium changes when a parameter changes. They are not derivatives of the transient state with respect to time. For \(k_{\mathrm{out}}>0\), start from the equilibrium condition \(k_{\mathrm{in}}-k_{\mathrm{out}}R^*=0\). Solving gives \(R^*=k_{\mathrm{in}}/k_{\mathrm{out}}\). Differentiating this expression with respect to each parameter gives the expected sensitivities:

The first derivative is positive when \(k_{\mathrm{out}}>0\): increasing the input raises the equilibrium. The second derivative is nonpositive for \(k_{\mathrm{in}}\geq 0\) and \(k_{\mathrm{out}}>0\); it is strictly negative when \(k_{\mathrm{in}}>0\), and equals zero when \(k_{\mathrm{in}}=0\). In particular, increasing the elimination rate does not change the zero equilibrium when there is no input. These sign and boundary checks should be part of the tests, not assumed from intuition alone.

The same derivatives can be obtained by differentiating the equilibrium condition. For example, differentiating \(k_{\mathrm{in}}-k_{\mathrm{out}}R^*=0\) with respect to \(k_{\mathrm{out}}\), while holding \(k_{\mathrm{in}}\) fixed, gives \(-R^*-k_{\mathrm{out}}(\partial R^*/\partial k_{\mathrm{out}})=0\). Therefore \(\partial R^*/\partial k_{\mathrm{out}}=-R^*/k_{\mathrm{out}}=-k_{\mathrm{in}}/k_{\mathrm{out}}^2\). This provides an independent algebraic check of the quotient-rule result.

Designing Tests and Diagnosing Failures

A compact test suite should cover distinct properties rather than repeat the same parameter setting. The following sequence is appropriate for this model:

Failures can help localize the bug. A nonzero derivative at the zero-input equilibrium suggests a residual input, a wrong state value, or an indexing error. A positive derivative at a positive state with no input suggests the sign of the loss term may be reversed. A finite steady-state value returned at \(k_{\mathrm{out}}=0\) with positive \(k_{\mathrm{in}}\) suggests an invalid division or an inappropriate special case. A derivative with the wrong units or sign can indicate that a parameter was confused with a state or that the derivative formula was implemented incorrectly.

A unit test of the right-hand side is different from integrating the ODE and comparing a trajectory with an analytic solution. The earlier material on testing ODE functions against analytic solutions addresses trajectory comparisons; equilibrium and limiting-case tests check targeted properties, often with less dependence on solver settings. Neither establishes that every possible parameter combination is correct. Tests should use multiple parameter values, including boundary cases, and should be supplemented by inspection of units and model equations.

A further limitation is that these checks use a one-state, constant-rate model with a known equilibrium. In a coupled model, an equilibrium may require solving a system of equations, and a zero-rate boundary may produce multiple equilibria or no equilibrium. The principles remain useful, but expected results must be derived for the model actually being tested. Do not apply \(k_{\mathrm{in}}/k_{\mathrm{out}}\) to a model whose input or loss terms have a different structure.