Algebraic Elimination Using Conserved Totals

Why Use a Conserved Total?

A reaction network may track several species even though their amounts cannot vary independently. If a particular weighted sum of species is conserved, one species can sometimes be recovered from the others and the known total. This removes a state variable without approximating the reaction kinetics or assuming that any reaction is fast.

The reduction is exact only on the set of states satisfying the conserved total. The eliminated species must therefore be reconstructed at every reduced state, and the reduced solver must stay in a domain where that reconstruction is physically feasible. These conditions distinguish algebraic elimination from an approximation to the dynamics.

Conserved Totals and the Reduced Equations

Definition: For a reaction-network ODE \(\dot{x}=S v(x)\), a vector \(c\) defines a conserved total if \(c^{\mathsf T}S=0\). Along any solution, \(c^{\mathsf T}x=T\) is then constant, where \(T=c^{\mathsf T}x(0)\). Here \(x\) is the vector of species amounts or concentrations, \(S\) is the stoichiometric matrix, and \(v(x)\) is the vector of reaction rates.

For a single conserved total, write \(c^{\mathsf T}x=\sum_{i=1}^{m}c_i x_i=T\). If \(c_m\ne 0\), designate \(x_m\) as the dependent species and solve the total for it:

x_m=h(y)=\frac{T-\sum_{i=1}^{m-1}c_i x_i}{c_m}, \qquad y=(x_1,\ldots,x_{m-1})^{\mathsf T}.

Let \(f(x)=S v(x)\), so that \(\dot{x}=f(x)\). Substitute \(x_m=h(y)\) into the first \(m-1\) equations. The reduced state \(y\) then satisfies:

\dot{y}=g(y)= \begin{pmatrix} f_1(y,h(y))\\ \vdots\\ f_{m-1}(y,h(y)) \end{pmatrix}.

This reduction presumes that the initial condition satisfies the stated total, that the total is genuinely conserved by the model, and that the same amount or concentration convention is used throughout. In particular, if concentrations are used, the conservation relation must be valid for the compartment volumes and transfers represented in the model. The dependent species can be chosen only when its coefficient in the total is nonzero. If several totals are available, the resulting linear constraints may allow several species to be eliminated; each elimination must leave a solvable, independent set of constraints.

Why the Elimination Is Exact

Along a solution of the full ODE, differentiating the total gives:

\frac{d}{dt}(c^{\mathsf T}x) =c^{\mathsf T}\dot{x} =c^{\mathsf T}S v(x) =0.

Thus, a solution starting on \(c^{\mathsf T}x=T\) remains on that set. Its first \(m-1\) components satisfy the reduced equations by direct substitution. Conversely, for a reduced solution \(y(t)\), reconstruct \(x_m=h(y(t))\). Differentiating the conservation relation gives:

c_m\dot{x}_m=-\sum_{i=1}^{m-1}c_i\dot{x}_i.

Since \(c^{\mathsf T}f(x)=0\), the full vector field obeys \(c_m f_m(x)=-\sum_{i=1}^{m-1}c_i f_i(x)\). The reduced equations set \(\dot{x}_i=f_i(x)\) for \(i<m\), so the reconstructed derivative satisfies \(\dot{x}_m=f_m(x)\). This proves that the reconstructed trajectory solves the omitted equation too. The argument requires \(c_m\ne0\) and the same conserved total throughout; it does not require a fast-reaction or steady-state assumption.

Worked Example: Reversible Conversion Between Two Species

Consider the hypothetical reaction \(A\rightleftharpoons B\), with forward rate constant \(k_f=0.3\ \mathrm{min}^{-1}\), reverse rate constant \(k_r=0.1\ \mathrm{min}^{-1}\), and initial amounts \(A(0)=2\ \mathrm{nM}\), \(B(0)=8\ \mathrm{nM}\). The mass-action equations are \(\dot{A}=-k_f A+k_r B\) and \(\dot{B}=k_f A-k_r B\). Adding them gives \(\frac{d}{dt}(A+B)=0\), so \(T=10\ \mathrm{nM}\).

Eliminate \(B\) using \(B=10-A\). The reduced equation is:

\dot{A}=-0.3A+0.1(10-A)=1-0.4A \quad \mathrm{nM\,min}^{-1}.

At the initial state, \(\dot{A}=1-0.4(2)=0.2\ \mathrm{nM\,min}^{-1}\). The original equation gives the same rate: \(-0.3(2)+0.1(8)=-0.6+0.8=0.2\ \mathrm{nM\,min}^{-1}\). The solution is \(A(t)=2.5-0.5e^{-0.4t}\ \mathrm{nM}\), with \(t\) in minutes, and reconstruction gives \(B(t)=10-A(t)\). For example, at \(t=5\ \mathrm{min}\), \(A\approx2.432\ \mathrm{nM}\) and \(B\approx7.568\ \mathrm{nM}\); both are nonnegative and sum to \(10\ \mathrm{nM}\).

Feasibility and Nonnegative Reconstruction

Algebraic elimination alone does not guarantee that reconstructed species are nonnegative. For a total with nonnegative coefficients, the reduced state must satisfy the constraint that makes \(h(y)\ge0\), as well as \(x_i\ge0\) for the retained species. These inequalities define the feasible part of the conserved set. For instance, in the two-species example, \(0\le A\le10\ \mathrm{nM}\) is required to have both \(A\ge0\) and \(B=10-A\ge0\).

For reaction systems with nonnegative reaction rates that vanish when a required reactant is absent, the full nonnegative state space is forward-invariant: a reaction cannot continue consuming a species once that species has reached zero. If the initial state is nonnegative and satisfies the conserved total, the exact full trajectory remains nonnegative, so its reduced representation and reconstruction do as well. Numerical integration is a separate matter: a solver may step slightly outside the feasible set, especially near a boundary. Check both the conservation residual and every retained and reconstructed species over time.

Worked Example: Two Conserved Totals in a Reversible Reaction

Consider the hypothetical reaction \(X+Y\rightleftharpoons Z\), with forward rate \(k_fXY\) and reverse rate \(k_rZ\). Let \(k_f=0.02\ \mathrm{nM}^{-1}\mathrm{min}^{-1}\), \(k_r=0.1\ \mathrm{min}^{-1}\), and initial concentrations \(X=3\ \mathrm{nM}\), \(Y=8\ \mathrm{nM}\), \(Z=2\ \mathrm{nM}\). The stoichiometric change is \((-1,-1,+1)^{\mathsf T}\), so \(X+Z\) and \(Y+Z\) are conserved. Their totals are \(5\ \mathrm{nM}\) and \(10\ \mathrm{nM}\), respectively.

Use \(Z\) as the single reduced state. The totals reconstruct the other species as \(X=5-Z\) and \(Y=10-Z\), with concentrations in \(\mathrm{nM}\). The reduced equation is:

\dot{Z}=0.02(5-Z)(10-Z)-0.1Z \quad \mathrm{nM\,min}^{-1}.

At \(Z=2\ \mathrm{nM}\), the reconstructed values are \(X=3\ \mathrm{nM}\) and \(Y=8\ \mathrm{nM}\). Substitution gives:

\dot{Z}=(0.02)(3)(8)-(0.1)(2) =0.48-0.20 =0.28\ \mathrm{nM\,min}^{-1}.

The feasible interval is \(0\le Z\le5\ \mathrm{nM}\): \(Z\) must be nonnegative, and \(X=5-Z\) must also be nonnegative. The other reconstructed species is nonnegative throughout this interval because \(Y=10-Z\ge5\ \mathrm{nM}\). Thus this reduced state and its feasible interval preserve nonnegative concentrations for all three species.

Worked Example: A Weighted Conserved Total

For a hypothetical reaction \(2A\rightarrow B\), let the reaction rate be \(kA^2\), with \(k=0.05\ \mathrm{nM}^{-1}\mathrm{min}^{-1}\). The equations are \(\dot{A}=-2kA^2\) and \(\dot{B}=kA^2\). The weighted total \(A+2B\) is conserved because \(\dot{A}+2\dot{B}=-2kA^2+2kA^2=0\). Suppose \(A(0)=4\ \mathrm{nM}\) and \(B(0)=3\ \mathrm{nM}\), giving \(T=4+2(3)=10\ \mathrm{nM}\).

Eliminating \(B\) gives \(B=(10-A)/2\), while the reduced equation is \(\dot{A}=-2(0.05)A^2\). The rate at the initial state is \(\dot{A}=-2(0.05)(4^2)=-1.6\ \mathrm{nM\,min}^{-1}\). For \(A\ge0\), the reconstruction is nonnegative exactly when \(A\le10\ \mathrm{nM}\), so the feasible interval is \(0\le A\le10\ \mathrm{nM}\).

The reduced solution is \(A(t)=A(0)/(1+2kA(0)t)\). At \(t=2\ \mathrm{min}\), \(A=4/[1+2(0.05)(4)(2)]=4/1.8\approx2.222\ \mathrm{nM}\), and \(B=(10-2.222)/2\approx3.889\ \mathrm{nM}\). Directly, \(kA^2\approx0.05(2.222)^2\approx0.2469\ \mathrm{nM\,min}^{-1}\), which is also \(-\dot{A}/2\). The weighted total is \(2.222+2(3.889)=10.000\ \mathrm{nM}\), up to rounding.

Diagnostics, Alternatives, and Limitations

A practical check is to reconstruct the eliminated species at every solver output and monitor the total residual \(c^{\mathsf T}x-T\). Also evaluate the omitted full equation at reconstructed states: its right-hand side should agree with the derivative implied by the conservation relation. A substantial mismatch can indicate a coding error in the substitution, an incorrect stoichiometric total, or inconsistent units. Negative reconstructed amounts indicate either that the reduced state has left the feasible domain or that the initial totals or model equations were specified inconsistently.

The main alternatives serve different purposes. Keeping the full ODE avoids choosing dependent coordinates, but retains a redundant state and may permit small numerical drift in a conserved total. Algebraic elimination removes that redundancy exactly, but makes feasible-domain constraints explicit and can make equations more complicated. A timescale-based reduction instead approximates dynamics under additional assumptions; conservation-based elimination does not. It is also limited to totals that are actually conserved: if the model includes an input, loss, or transfer that changes the proposed total, the constant-total reduction is invalid unless the modeled system is expanded or the changing total is handled explicitly.

Common pitfall: A conserved total does not make every choice of retained state physically valid. For example, \(B=T-A\) is negative if \(A>T\). The correction is to derive and enforce the feasible domain from nonnegativity of every reconstructed and retained species, rather than treating the algebraic formula as valid for arbitrary solver states.

Key takeaway: A linear conserved total provides an exact algebraic expression for a dependent species. Substitute that expression into the retained equations, verify the omitted equation through the conservation identity, and restrict the reduced state to values that reconstruct nonnegative species.

Check Your Understanding

Answer each question, then check the solutions below.

For \(A\rightleftharpoons B\) with \(A+B=T\), if \(A\) is retained, what is the reconstructed \(B\), and what condition makes it nonnegative?
In \(X+Y\rightleftharpoons Z\), why can the three species be represented by one reduced state when both \(X+Z\) and \(Y+Z\) are fixed?
For \(2A\rightarrow B\) with rate \(kA^2\), verify that \(A+2B\) is conserved.
Name two checks that can reveal an error or failure in a numerically integrated reduced model.

Solutions. (1) \(B=T-A\), which is nonnegative when \(A\le T\); also require \(A\ge0\) if both species are to be nonnegative. (2) The two independent totals reconstruct \(X=T_{XZ}-Z\) and \(Y=T_{YZ}-Z\), leaving only \(Z\) free, subject to nonnegativity constraints. (3) Differentiating gives \(\frac{d}{dt}(A+2B)=-2kA^2+2kA^2=0\). (4) Monitor the conserved-total residual and check that all reconstructed species remain nonnegative; comparing the omitted equation with the derivative implied by the total is another useful check.