Analytical Reference Solutions for ODE Verification

Why Analytical References Matter

A QSP component can produce smooth, plausible trajectories and still implement the wrong equation. Verification asks whether the executable code solves the mathematical model that was specified. When a component has an exact solution, that solution provides a reference independent of the numerical solver: evaluate both at the same times and with the same parameters, inputs, units, and initial conditions, then quantify their difference.

Analytical references are especially useful for small components embedded in larger models. A turnover state, a binding process under fixed ligand concentration, or a linear compartment subsystem can be tested in isolation before interactions make discrepancies harder to locate. Agreement does not establish that the equations adequately represent biology or that the model is suitable for a particular prediction; it establishes a narrower claim about implementation of the tested equations under the tested conditions.

First-Order Turnover Under Constant Input

Consider a state \(R(t)\), such as the amount or concentration of a turnover-controlled quantity, governed by constant input and first-order loss:

Assume \(k_{\mathrm{in}}\) and \(k_{\mathrm{out}}\) are constant, \(k_{\mathrm{out}}>0\), and there are no time-varying inputs, events, or additional coupled states. If \(R\) is measured in units \(U\), then \(k_{\mathrm{in}}\) has units \(U/\text{time}\) and \(k_{\mathrm{out}}\) has units \(1/\text{time}\). The steady state is \(R_{\mathrm{ss}}=k_{\mathrm{in}}/k_{\mathrm{out}}\). Subtracting this steady state from the differential equation gives \(d(R-R_{\mathrm{ss}})/dt=-k_{\mathrm{out}}(R-R_{\mathrm{ss}})\), whose solution is exponential decay of the deviation from steady state:

Reversible Binding With Fixed Free Ligand

A useful binding reference is available when free ligand concentration \(L\) is held constant and receptor total concentration \(R_{\mathrm{tot}}\) is fixed. Let \(C(t)\) be the concentration of complex. For association rate constant \(k_{\mathrm{on}}\) and dissociation rate constant \(k_{\mathrm{off}}\), the model is

Here \(k_{\mathrm{on}}\) has units concentration\(^{-1}\) time\(^{-1}\), and \(k_{\mathrm{off}}\) has units time\(^{-1}\). Because \(L\) is fixed, the equation is linear in \(C\). Define \(k_{\mathrm{obs}}=k_{\mathrm{on}}L+k_{\mathrm{off}}\). Setting the derivative to zero gives \(C_{\mathrm{ss}}=k_{\mathrm{on}}LR_{\mathrm{tot}}/k_{\mathrm{obs}}\). The difference \(C-C_{\mathrm{ss}}\) obeys \(d(C-C_{\mathrm{ss}})/dt=-k_{\mathrm{obs}}(C-C_{\mathrm{ss}})\), giving the exact trajectory below.

This is a pseudo-first-order reference, not a general solution for a closed system in which free ligand is depleted as complex forms. Using it for a closed, mass-balanced binding system would compare the implementation with the wrong equations.

Two-Compartment Linear Model

Linear compartment systems can also have exact reference trajectories. Consider a hypothetical IV bolus into the central compartment, with central amount \(A_1\) and peripheral amount \(A_2\):

The rate constants \(k_{10}, k_{12}, k_{21}\) have units time\(^{-1}\). The initial conditions for a central bolus of dose \(D\) are \(A_1(0)=D\) and \(A_2(0)=0\). With constant rate constants, this is a linear system with constant coefficients. Its modes are exponentials associated with the eigenvalues of its coefficient matrix. For \(k_{10}=0.40\), \(k_{12}=0.30\), and \(k_{21}=0.20\) h\(^{-1}\), the characteristic equation is

Its roots are \(-0.10\) and \(-0.80\) h\(^{-1}\). Substitution into the differential equations, together with the initial conditions, gives the following exact amount trajectories for any bolus dose \(D\):

The expressions satisfy the initial conditions: \(A_1(0)=D\) and \(A_2(0)=0\). They also satisfy the equations on differentiation. For example, differentiating \(A_1\) gives \(-0.10(D/7)e^{-0.10t}-0.80(6D/7)e^{-0.80t}\); inserting the displayed \(A_1,A_2\) into \(-(0.70)A_1+0.20A_2\) gives the same expression. The second equation is verified by differentiating \(A_2\) and substituting the same trajectories. This provides a direct check of the derived reference.

How to Make the Comparison

A verification comparison should be reproducible and localize errors rather than merely report that two curves differ. Choose time points that include the initial time, several interior times, and later times at which the dynamics have changed appreciably. Evaluate the exact expression at precisely the output times from the executable model. For each state and time, calculate an error such as

Here \(y_{\mathrm{sim}}\) and \(y_{\mathrm{ref}}\) denote simulated and reference values in the same units; \(\mathrm{atol}\) is an absolute tolerance in those units and \(\mathrm{rtol}\) is a dimensionless relative tolerance. A scaled error no greater than one satisfies this chosen pointwise tolerance. Near a zero reference value, the absolute tolerance prevents an unstable relative-error ratio. Tolerances should reflect the required numerical accuracy and the precision of the implementation, not be selected after seeing the result.

Diagnostics, Alternatives, and Limitations

The pattern of error often points to the source. A nonzero error at time zero suggests an initial-condition, state-mapping, or output-alignment problem. A discrepancy that grows steadily can indicate a wrong rate constant, a unit conversion error, or an implementation defect. A correct early trajectory but wrong long-time value may indicate a steady-state or input error. A mismatch in only one state of a coupled system may indicate an indexing or output-mapping mistake. Conversely, small differences consistent with the chosen tolerance may arise from floating-point arithmetic and the numerical solver's approximation.

The analytical reference is the strongest direct comparison when the same mathematical assumptions yield a closed-form solution. A separately implemented high-accuracy numerical solution can serve as a useful cross-check when the equation is more complicated, but it is not an exact reference and could share a coding error with the implementation under test. Visual overlays are useful for seeing broad trajectory shape, but can conceal small, systematic discrepancies. Analytical checks also complement, rather than replace, code review and broader testing: they cover only the equations and conditions represented by the selected test cases.

A frequent mistake is to compare a simulation with a familiar formula without checking whether the assumptions match. For example, applying the fixed-free-ligand binding solution to a closed ligand–receptor system, or comparing a concentration output with a reference amount, creates a discrepancy even if each calculation is internally correct. The correction is to write down the exact implemented equations and state definitions first, verify units and conditions, and only then derive or select the reference. Another pitfall is to regard a close match at one time as proof of correctness; the comparison should cover multiple states and times, including the initial condition and the relevant time course.

These examples are deliberately simple. Their role is not to establish biological validity, assess qualification evidence, or test all behaviors of a complete QSP model. They provide transparent, quantitative checks of selected ODE components under specified conditions. More complex components may not admit a convenient closed form, and even a correct component-level match does not guarantee that interactions, event handling, or the full executable model are correct.