Variable elimination in post-translational modification reaction networks with mass-action kinetics

We define a subclass of Chemical Reaction Networks called Post-Translational Modification systems. Important biological examples of such systems include MAPK cascades and two-component systems which are well-studied experimentally as well as theoretically. The steady states of such a system are solutions to a system of polynomial equations with as many variables as equations. Even for small systems the task of finding the solutions is daunting. We develop a mathematical framework based on the notion of a cut, which provides a linear elimination procedure to reduce the number of variables in the system. The steady states are parameterized algebraically by a set of “core” variables, and the non-negative steady states correspond to non-negative values of the core variables. Further, minimal cuts are the connected components in the species graph and provide conservation laws. A criterion for when a set of independent conservation laws can be derived from cuts is given.

💡 Research Summary

The paper addresses the challenging problem of determining steady‑state solutions for post‑translational modification (PTM) systems, a subclass of chemical reaction networks that model many important signaling pathways such as MAPK cascades and two‑component systems. In a PTM system the species are divided into substrates S and intermediate complexes Y, and only four types of reactions are allowed: (i) formation of a complex from two substrates (Rₐ), (ii) dissociation of a complex back to its two substrates (R_b), (iii) conversion between complexes (R_c), and (iv) direct conversion between substrates (R_d). With mass‑action kinetics the dynamics are described by a polynomial ODE system that is quadratic in the substrate concentrations and linear in the complex concentrations. The steady‑state equations are obtained by setting all time‑derivatives to zero, yielding a system of as many polynomial equations as unknowns, which is generally intractable even for modest network sizes.

The authors introduce the notion of a “cut” – a subset S_α ⊂ S that satisfies two conditions: (a) no two substrates in S_α interact directly, and (b) the set is closed under the 1‑link relation (if a substrate in S_α is 1‑linked to another substrate or to a complex, the linked element must also belong to S_α). Minimal cuts correspond precisely to the connected components of a graph built from the species and 1‑link edges; this graph is called a non‑interacting graph. By selecting a cut, the authors separate the species into “core” variables S \ S_α and the cut species S_α. Because the steady‑state equations for the complexes Y are linear, they can be solved explicitly in terms of the core variables, yielding rational expressions for every Y_k. Substituting these expressions back eliminates all Y‑variables and reduces the original polynomial system to a smaller system that involves only the core substrates. Thus, the steady‑state variety is parametrized algebraically by the core variables; any non‑negative choice of core variables produces a non‑negative steady state for the whole network.

A second major contribution is the systematic derivation of conservation laws from cuts. The stoichiometric subspace Γ is defined by the reaction vectors, and its orthogonal complement Γ^⊥ contains vectors that give linear combinations of species whose total amount is invariant. The authors prove that each minimal cut generates a vector in Γ^⊥, and therefore each non‑interacting graph yields an independent conservation law of the form
 ω_l = Σ_{S_i∈S_l} S_i + Σ_{Y_k∈Y_l} Y_k = constant.
These conserved quantities further reduce the dimensionality of the steady‑state problem because the constants are fixed by the initial total amounts. When the set of conservation laws obtained from all minimal cuts is independent, the total reduction equals the number of cuts, and the remaining degrees of freedom are exactly the core variables.

The theoretical framework is illustrated with several examples. A small network with five substrates and three complexes is worked out in detail: the authors eliminate Y₃ directly, solve the linear subsystem for Y₁ and Y₂, and identify the conserved quantity S₅ + Y₃. They then discuss how the same procedure scales to larger, biologically relevant systems such as multi‑layer MAPK cascades, phosphorelay two‑component systems, and networks that include self‑interactions. In each case, appropriate cuts are identified, the complexes are expressed as rational functions of the core substrates, and the resulting conservation laws are derived from the connected components of the associated non‑interacting graphs.

The paper also acknowledges situations where a cut may not exist—for instance, when every substrate is linked to every other substrate—so that the non‑interacting graph would have to contain all species, violating the cut definition. In such cases the method does not directly yield independent conservation laws, and additional algebraic or graph‑theoretic techniques would be required.

Importantly, the authors extend the earlier work of Thomson and Gunawardena (2017), which handled only single‑step phosphorylation cycles, to a much broader class of PTM networks, including cascades, phosphorelays, and self‑interacting reactions. By treating kinetic parameters as symbolic positive constants, the approach avoids the need for numerical parameter values, making it especially valuable for biological systems where rate constants are often unknown or only roughly estimated.

In summary, the paper provides a rigorous, graph‑based variable‑elimination scheme for PTM reaction networks under mass‑action kinetics. By defining cuts, expressing intermediate complexes as rational functions of core substrates, and extracting conservation laws from the structure of the species graph, the authors achieve a substantial reduction in the number of variables needed to describe steady states. This framework offers a powerful analytical tool for studying multistationarity, parameter dependence, and the qualitative behavior of complex signaling networks without resorting to exhaustive numerical simulations.