Deriving Conserved Moieties from the Stoichiometric Matrix

Why Look for Conserved Moieties?

A reaction network can change the amounts of individual species while preserving the total amount of some underlying material. For example, conversion of one species into another changes both species amounts, but may preserve their sum. In a QSP model, identifying such conserved combinations provides a structural check on the model and clarifies which quantities the reactions can change.

The stoichiometric matrix gives a systematic way to find these combinations. The key is to look for weights on the species that make the weighted sum unchanged by every reaction. These weights are found in the matrix’s left null space.

Definition: For a network with \(n\) species and \(m\) reactions, the stoichiometric matrix \(S\) is an \(n \times m\) matrix. Its entry \(S_{ij}\) is the net change in the amount of species \(i\) when reaction \(j\) occurs once, using negative entries for consumed species and positive entries for produced species. If \(x\) is the vector of species amounts and \(v(x)\) is the vector of reaction rates, the closed-network balance equations are \(\frac{dx}{dt}=S v(x)\). A vector \(c\) is a conserved-moiety vector if \(c^T S=0\). It defines the conserved quantity \(c^T x\).

Here \(x\) contains amounts, not concentrations. Each component of \(c\) weights the amount of its corresponding species. For example, if one molecule of species \(D\) contains two units of a moiety and one molecule of \(A\) contains one, the corresponding conserved amount is \(A+2D\). All amounts in a particular weighted sum must use a consistent amount basis.

The closed-network assumption matters: the balance equation includes only the reactions represented in \(S\), with no processes that add or remove the tracked species from the network. The reaction rates may depend nonlinearly on the state; the conservation result does not require mass-action kinetics or constant rates.

From the Matrix to an Invariant

Let \(c\) be a fixed vector of species weights. Along a model trajectory, the rate of change of the weighted amount is obtained by differentiating \(c^T x\). Since \(c\) is fixed, the chain rule gives:

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

If \(c^T S=0\), this derivative is zero at every state, regardless of the values of the reaction rates. Therefore \(c^T x\) stays constant in time. The same condition can be understood one reaction at a time: the change in the weighted amount caused by reaction \(j\) is the dot product of \(c\) with column \(j\) of \(S\). It is zero exactly when that column contributes zero to \(c^T S\).

c^T S=0 \quad\Longrightarrow\quad c^T x(t)=c^T x(0).

This condition must hold for every column of \(S\), not just for the net change observed at one particular state. A cancellation that happens only because of a special rate or initial condition does not establish a stoichiometric conserved moiety.

Worked Examples

Worked Example: Reversible Interconversion

Consider two species, \(A\) and \(B\), that interconvert through \(A\rightarrow B\) and \(B\rightarrow A\). The species order is \((A,B)\), and the reaction order is forward then reverse. The stoichiometric matrix is:

S= \begin{bmatrix} -1 & 1\\ 1 & -1 \end{bmatrix}.

Write \(c=(c_A,c_B)^T\) and solve \(c^T S=0\). Multiplication gives:

c^T S = \begin{bmatrix}c_A&c_B\end{bmatrix} \begin{bmatrix}-1&1\\1&-1\end{bmatrix} = \begin{bmatrix}-c_A+c_B&c_A-c_B\end{bmatrix}.

Both entries are zero when \(c_A=c_B\). Choosing \(c_A=c_B=1\) gives the invariant \(A+B\). Suppose the initial amounts are \(A=12\) nmol and \(B=3\) nmol. Then the conserved total is \(12+3=15\) nmol. For any forward and reverse rates \(v_1\) and \(v_2\), the species-amount rates are \(\frac{dA}{dt}=-v_1+v_2\) and \(\frac{dB}{dt}=v_1-v_2\). Their sum is zero, so the total remains 15 nmol.

Worked Example: Association of Two Distinct Species

Consider \(A+B\rightleftharpoons C\), with species order \((A,B,C)\) and reaction order forward then reverse. The matrix is:

S= \begin{bmatrix} -1&1\\ -1&1\\ 1&-1 \end{bmatrix}.

For \(c=(c_A,c_B,c_C)^T\), either column gives the same condition:

-c_A-c_B+c_C=0, \qquad \text{so} \qquad c_C=c_A+c_B.

There are two independent choices, for example \(c=(1,0,1)^T\) and \(c=(0,1,1)^T\). They give the invariants \(A+C\) and \(B+C\). The first counts the amount of \(A\)-type material, whether free as \(A\) or present in \(C\); the second does the corresponding accounting for \(B\)-type material.

Take hypothetical initial amounts \(A=8\) nmol, \(B=5\) nmol, and \(C=2\) nmol. The invariant totals are \(A+C=10\) nmol and \(B+C=7\) nmol. If the forward and reverse rates at some instant are \(v_1=0.6\) nmol/min and \(v_2=0.2\) nmol/min, then:

\frac{dx}{dt} =S \begin{bmatrix}0.6\\0.2\end{bmatrix} = \begin{bmatrix}-0.4\\-0.4\\0.4\end{bmatrix} \text{ nmol/min}.

Thus \(\frac{d}{dt}(A+C)=-0.4+0.4=0\) nmol/min and \(\frac{d}{dt}(B+C)=-0.4+0.4=0\) nmol/min. Both totals are unchanged, even though each species amount is changing.

Worked Example: Dimerization and Stoichiometric Coefficients

Now consider \(2A\rightleftharpoons D\), with species order \((A,D)\) and forward and reverse reaction columns. The matrix is:

S= \begin{bmatrix} -2&2\\ 1&-1 \end{bmatrix}.

For \(c=(c_A,c_D)^T\), the left-null-space condition is:

c^T S = \begin{bmatrix}-2c_A+c_D&2c_A-c_D\end{bmatrix} =0.

Therefore \(c_D=2c_A\), and choosing \(c_A=1\) gives \(c=(1,2)^T\). The invariant is \(A+2D\), not \(A+D\): one unit of \(D\) represents two units of the material in \(A\).

Suppose \(A=10\) μmol and \(D=3\) μmol initially. The conserved amount is \(10+2(3)=16\) μmol of \(A\)-equivalent material. If the forward and reverse rates are \(0.4\) μmol of reaction events per minute and \(0.1\) μmol of reaction events per minute, the net reaction rate is \(0.3\) μmol of events per minute. Consequently, \(\frac{dA}{dt}=-2(0.3)=-0.6\) μmol/min and \(\frac{dD}{dt}=0.3\) μmol/min. The weighted rate is \(-0.6+2(0.3)=0\) μmol/min, as required.

Independent Invariants and Matrix Rank

The equations \(c^T S=0\) are homogeneous linear equations in the species weights. Solving them yields a basis for the left null space. Any linear combination of basis vectors is also a conserved-moiety vector, so the basis is a compact description of all independent stoichiometric invariants.

If \(S\) has \(n\) rows, its left null space has dimension \(n-\operatorname{rank}(S)\). Therefore, the number of linearly independent stoichiometric invariants is \(n-\operatorname{rank}(S)\). In the association example, \(n=3\) and the two reaction columns are negatives of each other, so \(\operatorname{rank}(S)=1\). There are \(3-1=2\) independent invariants, matching the two totals found directly.

The particular basis is not unique. Multiplying an invariant by a nonzero constant, or taking a linear combination of two invariants, produces another valid invariant. For biological interpretation, a convenient basis often uses simple nonnegative integer weights. A mathematical null-space vector is not automatically a biologically meaningful moiety: its weights must also make sense for the species being represented.

Derivation: By the rank-nullity theorem, the nullity of \(S^T\) is \(n-\operatorname{rank}(S^T)\). Since \(\operatorname{rank}(S^T)=\operatorname{rank}(S)\), the dimension of the set of vectors satisfying \(S^T c=0\), equivalently \(c^T S=0\), is \(n-\operatorname{rank}(S)\). Each independent vector in this space supplies one independent conserved weighted sum.

How to Check a Candidate Invariant

Fix the species and reaction order.
Record which species each row represents and which reaction each column represents. Interpret every sign and stoichiometric coefficient using that ordering.

Propose or solve for the weights.
For a proposed \(c\), compute \(c^T S\). To find all candidates, solve the linear system \(S^T c=0\).

Verify every reaction column.
Check that the weighted sum of stoichiometric changes is zero for each reaction separately, then interpret \(c^T x\) using the species amounts.

Check the model’s scope.
Confirm that the tracked network is closed with respect to the reactions represented. A valid invariant for the listed reaction columns need not remain constant if the model includes additional processes outside those columns.

Diagnostics, Alternatives, and Common Mistakes

A useful algebraic diagnostic is to calculate \(c^T S\) and inspect every entry. A nonzero entry means that at least one reaction changes the proposed weighted total. If a result seems biologically implausible, first check for a species-order mismatch, a missing stoichiometric coefficient, or a sign error in the matrix. These errors can produce a mathematically correct null-space calculation for the wrong matrix.

Finding invariants from the matrix is more systematic than guessing conserved totals by inspection, particularly when a network has many species and reactions. Directly tracking a known material can still be a useful independent check: assign each species the number of units it contains and verify that each reaction preserves the total. Simulation provides another check: calculate \(c^T x(t)\) at several times and see whether it stays constant. Simulation alone, however, does not establish the structural result; a numerical trajectory might hide a matrix error or reflect a special rate pattern.

Using the right null space instead of the left. The equation \(S z=0\) concerns combinations of reaction columns, where \(z\) has one entry per reaction. Conserved species weights instead satisfy \(c^T S=0\), where \(c\) has one entry per species.
Assuming every species is individually conserved. Individual species amounts usually change. It is a weighted combination that remains constant when the reaction changes cancel under those weights.
Ignoring stoichiometric coefficients. For \(2A\rightleftharpoons D\), the conserved combination is \(A+2D\); \(A+D\) changes during the reaction.
Calling a special trajectory a structural invariant. A quantity that happens to remain unchanged for one parameter choice or initial condition is not necessarily conserved by the stoichiometry. The defining test is \(c^T S=0\).

The method identifies invariants implied by the specified reaction stoichiometry. It does not by itself establish that the species list or reaction network fully represents the biology. Its result also depends on the closed-network balance equation: a conserved sum for the listed reactions may change if the model is expanded to include processes that alter the tracked amounts.

Check Your Understanding

Try each problem, then compare your reasoning with the solutions.

For \(A\rightarrow B\) with species order \((A,B)\), write \(S\) and give one conserved-moiety vector.
For \(A+B\rightarrow C\), use species order \((A,B,C)\). Find two independent conserved combinations.
For \(2A\rightarrow D\), suppose \(A=6\) μmol and \(D=2\) μmol. What is the conserved amount in \(A\)-equivalent units?
What does \(S z=0\) describe, and why is it not the test for a conserved species combination?

Solutions. (1) \(S=\begin{bmatrix}-1\\1\end{bmatrix}\). The vector \(c=(1,1)^T\) satisfies \(c^T S=-1+1=0\), so \(A+B\) is conserved. (2) The single reaction column is \((-1,-1,1)^T\). The condition is \(-c_A-c_B+c_C=0\), so \(c_C=c_A+c_B\). Choosing \(c=(1,0,1)^T\) and \(c=(0,1,1)^T\) gives \(A+C\) and \(B+C\). (3) The invariant is \(A+2D\), which equals \(6+2(2)=10\) μmol of \(A\)-equivalent material. (4) \(S z=0\) describes combinations of reaction columns whose net stoichiometric change is zero; \(z\) indexes reactions. A conserved species weighting instead has one weight per species and must satisfy \(c^T S=0\).