Why State Indices Need Explicit Names
An ODE model is usually written in terms of named biological quantities, while a numerical integrator receives a vector of numbers. Translating between those representations is not merely a coding detail: if the value for one species is read as another, or a derivative is returned in the wrong position, the program can integrate a different model from the one expressed by the equations. A fixed, documented mapping makes the translation inspectable and testable.
This tutorial focuses on the mapping between biological states and vector indices. It does not change the biological equations. The task is to ensure that each named state occupies exactly one position and that the derivative returned at that position belongs to that same state.
The mapping is one-to-one: no two state names share an index, and every index from 1 through \(m\) is assigned to a state. It must remain fixed wherever the state vector and its derivative are read or written. Mathematical notation commonly uses indices starting at 1. A programming language may start array indices at 0; in that case, the implementation needs a consistent translation, such as storing mathematical index \(i\) at array position \(i-1\). The biological mapping should still be documented in one unambiguous convention.
Specify the Mapping Before Implementing the Equations
A useful specification lists each biological name, its vector index, the meaning of its value, and the meaning and units of its derivative. Consider a two-state reversible transfer model. Let \(A\) and \(B\) be amounts in two compartments, measured in nmol. Material transfers from \(A\) to \(B\) at rate \(k_{AB}A\), and from \(B\) to \(A\) at rate \(k_{BA}B\). Assume a well-mixed amount in each compartment, constant positive rate coefficients \(k_{AB}\) and \(k_{BA}\) in \(\mathrm{h}^{-1}\), and no production, degradation, or other inputs and outputs. These assumptions give the named equations
The transfer terms appear with opposite signs because each transfer removes amount from one state and adds the same amount to the other. One valid mapping is \(A\mapsto 1\), \(B\mapsto 2\). Under this mapping, the state vector is \(x=(A,B)^\mathsf{T}\), and the derivative vector must be ordered as \((dA/dt,dB/dt)^\mathsf{T}\).
| Biological state | Vector index | Stored value | Derivative returned at that index |
|---|---|---|---|
| \(A\) | 1 | Amount of \(A\), nmol | \(-k_{AB}A+k_{BA}B\), nmol/h |
| \(B\) | 2 | Amount of \(B\), nmol | \(k_{AB}A-k_{BA}B\), nmol/h |
The mapping does not have to follow the order in which states appear in the equations. For instance, \(B\mapsto 1\), \(A\mapsto 2\) is also valid. What matters is that both the input state vector and output derivative vector follow the same declared order. A named mapping table is therefore part of the model specification, not an informal memory aid.
Worked Example: Calculate the Derivative Vector in the Declared Order
Use \(A\mapsto 1\), \(B\mapsto 2\), with \(A=8\) nmol, \(B=2\) nmol, \(k_{AB}=0.30\ \mathrm{h}^{-1}\), and \(k_{BA}=0.10\ \mathrm{h}^{-1}\). Find the derivative vector.
First calculate each directed transfer rate. From \(A\) to \(B\), the rate is \(k_{AB}A=(0.30\ \mathrm{h}^{-1})(8\ \mathrm{nmol})=2.4\ \mathrm{nmol/h}\). From \(B\) to \(A\), it is \(k_{BA}B=(0.10\ \mathrm{h}^{-1})(2\ \mathrm{nmol})=0.2\ \mathrm{nmol/h}\).
Thus \(dA/dt=-2.4+0.2=-2.2\ \mathrm{nmol/h}\), while \(dB/dt=2.4-0.2=2.2\ \mathrm{nmol/h}\). Because the declared order is \(A,B\), the state vector and returned derivative are
The first output is the derivative of \(A\), not simply “the first equation evaluated”; it is the derivative of whichever state the mapping assigns to index 1. Here those descriptions agree because index 1 is \(A\).
Why a Permuted Order Can Still Be Correct
A different ordering changes the arrangement of values, not the biological model, provided the equations are rearranged consistently. Let \(z=(A,B)^\mathsf{T}\) be the state vector in the order used to write the named equations, and let \(g(z)\) be their right-hand side. If the stored vector instead uses \(x=(B,A)^\mathsf{T}\), then \(x=Pz\), where
Since \(P^{-1}=P\), the derivative in the \(B,A\) order is \(P g(P^{-1}x)\). This follows by differentiating \(x=Pz\): \(dx/dt=P(dz/dt)=Pg(z)\), and substituting \(z=P^{-1}x\). The swap of values and the swap of derivative components must happen together.
For the reversible transfer model, \(x=(B,A)^\mathsf{T}\) therefore requires the first output to be \(dB/dt\) and the second to be \(dA/dt\). Returning the equations in \(A,B\) order while storing states in \(B,A\) order is not a harmless formatting difference: it assigns each computed rate to the wrong biological state.
Worked Example: Check a Reversed State Order
Use the same amounts and rate coefficients as above, but declare \(B\mapsto 1\) and \(A\mapsto 2\). Write the state and derivative vectors in that order.
The biological derivatives have not changed: \(dA/dt=-2.2\ \mathrm{nmol/h}\) and \(dB/dt=2.2\ \mathrm{nmol/h}\). The stored state is now \(x=(B,A)^\mathsf{T}=(2,8)^\mathsf{T}\ \mathrm{nmol}\). The output must be \(F=(dB/dt,dA/dt)^\mathsf{T}\), so
This is the correct output for the reversed mapping. Returning \((-2.2,2.2)^\mathsf{T}\) would be incorrect: it would tell the numerical method that \(B\) is decreasing and \(A\) is increasing, contrary to the rates calculated from the named equations at these state values.
A Direct Test for Swaps and Omissions
A dependable check starts with the name-to-index table and applies it to a deliberately simple state. For each state, calculate its derivative from the named equation, then place that result at the index assigned to that state. Separately, make a term ledger: every biological process that changes a state should appear in that state’s equation with the appropriate contribution. For a transfer between two states, the corresponding loss and gain should both be present.
For the closed reversible model, adding the equations gives a useful balance check:
Thus the total amount \(A+B\) is conserved under the stated assumptions. This check can reveal an inconsistent transfer pair, but it cannot by itself prove that the indices are correct. A swap of both state values and derivative outputs preserves the same total. A component-by-component comparison against the named equations is still required.
Likewise, a balance check can fail to detect a process omitted from both equations. The test must match the expected balance for the model as specified and should be paired with checks of individual state derivatives. A useful test state has different values in the two components: if \(A=B\), swapping their stored values leaves the input vector unchanged and can conceal a mapping error.
Worked Example: Detect a Missing Transfer Contribution
Extend the two-state model by adding a constant input \(u\) to \(A\) and an output from \(A\) at rate \(k_{\mathrm{out}}A\). Assume the same reversible transfers as before, with \(u=1.0\ \mathrm{nmol/h}\), \(k_{\mathrm{out}}=0.05\ \mathrm{h}^{-1}\), and the values from the earlier calculation. The complete equations are
For \(A=8\) nmol and \(B=2\) nmol, the \(A\) equation gives \(1.0-2.4+0.2-0.4=-1.6\ \mathrm{nmol/h}\). The \(B\) equation gives \(2.4-0.2=2.2\ \mathrm{nmol/h}\). With \(A\mapsto1\), \(B\mapsto2\), the expected derivative vector is \((-1.6,2.2)^\mathsf{T}\ \mathrm{nmol/h}\).
Suppose an implementation accidentally omits the reverse-transfer contribution \(+k_{BA}B\) from the \(A\) equation. It would return \((-1.8,2.2)^\mathsf{T}\ \mathrm{nmol/h}\), differing from the expected first component by \(-0.2\ \mathrm{nmol/h}\). The total balance implied by the complete equations is
The incorrect output sums to \(-1.8+2.2=0.4\ \mathrm{nmol/h}\), so it also fails this balance check. The discrepancy identifies a missing contribution, while checking each component against its named equation establishes which component is wrong.
Practical Checks and Representation Choices
Before integrating a model, verify the mapping at the interface between named states and vectors. These checks are specific to state indexing; they do not replace checks of the biological assumptions or equations.
- Check coverage and uniqueness. Confirm that every declared state appears once in the mapping and that no index is repeated or left unassigned.
- Check both directions. For each state name, identify its index; for each index, identify its state name. These two lookups should be inverses.
- Use unequal test values. Evaluate the right-hand side at distinct state values so exchanging components changes the expected result.
- Check each returned component by name. Compare the value at index \(\pi(s)\) with the derivative of state \(s\) calculated directly from its equation.
- Check paired contributions and balances. Where the model specifies a transfer or known total balance, verify that the relevant component terms agree with it.
An indexed vector is compact and works naturally with numerical integration routines, but its biological meaning depends on an external mapping. A name-based structure can make individual states easier to recognize and reduce accidental use of the wrong position, but it may require a conversion to the ordered vector expected by a numerical routine. A mapping table is useful with either representation: it documents the intended order and provides a reference for tests. Automatically generating both the state vector and derivative vector from a single declared ordering can reduce opportunities for inconsistent manual edits, but the generated order still needs to be checked against the biological names.
A common mistake is to assume that a derivative vector is correct because it has the right length, or because its components sum to an expected total. Neither condition establishes that each component belongs to the intended state. The correction is to test the name-to-index correspondence component by component, using unequal state values and the explicit mapping. A second pitfall is changing the order of the state vector in one place but not updating the derivative order; treat the mapping as fixed and update every dependent representation together.
Check Your Understanding
Try each problem, then compare your reasoning with the solutions below.
- For \(A\mapsto1\), \(B\mapsto2\), what biological derivative must be returned at vector index 2?
- In the reversible model, let \(A=4\) nmol, \(B=1\) nmol, \(k_{AB}=0.2\ \mathrm{h}^{-1}\), and \(k_{BA}=0.1\ \mathrm{h}^{-1}\). Calculate the derivative vector for the order \(A,B\).
- For the same values as question 2, write the state and derivative vectors if the declared order is \(B,A\).
- Why is testing the model only at \(A=B\) a weak way to detect a swap of the two state values?
- For the input-and-output model, what is the expected rate of change of \(A+B\) when \(u=1.0\ \mathrm{nmol/h}\), \(k_{\mathrm{out}}=0.05\ \mathrm{h}^{-1}\), and \(A=8\) nmol?
Solutions. 1. Index 2 corresponds to \(B\), so the returned component must be \(dB/dt\). 2. The transfers are \(k_{AB}A=(0.2)(4)=0.8\ \mathrm{nmol/h}\) and \(k_{BA}B=(0.1)(1)=0.1\ \mathrm{nmol/h}\). Therefore \(dA/dt=-0.8+0.1=-0.7\ \mathrm{nmol/h}\) and \(dB/dt=0.8-0.1=0.7\ \mathrm{nmol/h}\), giving \((-0.7,0.7)^\mathsf{T}\ \mathrm{nmol/h}\). 3. The state vector is \((B,A)^\mathsf{T}=(1,4)^\mathsf{T}\ \mathrm{nmol}\). The derivative vector follows the same order, so it is \((dB/dt,dA/dt)^\mathsf{T}=(0.7,-0.7)^\mathsf{T}\ \mathrm{nmol/h}\). 4. When \(A=B\), exchanging the two values leaves the state vector unchanged, so a test may not expose a state-reading swap. Use unequal values and check each output against its named equation. 5. The transfer terms cancel in the total balance, leaving \(d(A+B)/dt=u-k_{\mathrm{out}}A=1.0-(0.05)(8)=0.6\ \mathrm{nmol/h}\).