Given a polynomial ODE system, can you find a chemical reaction network that generates it? Yes — Hárs and Tóth showed in 1979 [1] that a polynomial ODE $\dot{x} = P(x)$ is realizable by a mass-action CRN if and only if it has the “lack of negative cross-effects” property (the negative terms in $\dot{x}_k$ must all contain $x_k$ as a factor), and they gave an explicit construction. Finding a CRN is easy.

But finding a minimal CRN — one with the fewest reactions or the fewest complexes — is a different problem entirely. When is that problem trivial, and when does it require real mathematics?

The Canonic Mechanism

The first sum collects all production — each reaction $y(q) \to y(p)$ at rate $k_{pq}\, x^{y(q)}$ increases species $m$ by $y^m(p)$ copies. The second sum collects all consumption — the same reaction decreases species $m$ by $y^m(q)$ copies. The positive and negative contributions can be lumped together, and when $y^m(q) = 0$ the consumption term vanishes entirely (species $m$ does not appear in complex $q$, so it is neither produced nor consumed by reactions from that complex).

This leads to the Hárs-Tóth necessary condition: in any kinetic differential equation, each negative monomial in $\dot{x}_m$ must contain $x_m$ as a factor (the “lack of negative cross-effects”). If a monomial appears with a negative coefficient in $\dot{x}_m$ but does not involve $x_m$, then it cannot come from degradation of species $m$ and no CRN can produce it.

The Hárs-Tóth construction (the canonic mechanism) is the converse: given any polynomial ODE satisfying this condition, write down one reaction per monomial term. A positive term $c \cdot x^{\alpha}$ in $\dot{x}_m$ means species $m$ is being produced at rate $c \cdot x^{\alpha}$, so write a reaction with source complex $\alpha$ whose product complex adds one copy of species $m$: $\alpha \to \alpha + e_m$ with rate constant $c$. A negative term $-c \cdot x^{\alpha}$ in $\dot{x}_m$ means species $m$ is being consumed, so write a reaction that removes one copy: $\alpha \to \alpha - e_m$ with rate constant $c$.

Consider the following four-species ODE from Hárs and Tóth’s original paper:

The canonic mechanism reads off one reaction per monomial term. For instance, the term $-2\alpha\, x^2 v$ in $\dot{x}$ means species $X$ is consumed at rate $2\alpha\, x^2 v$: the source complex is $2X + V$, and the reaction removes one $X$, giving $2X + V \to X + V$ with rate $2\alpha$. The term $+3\alpha\, x^2 v$ in $\dot{y}$ means species $Y$ is produced at rate $3\alpha\, x^2 v$: same source complex $2X + V$, and the reaction adds one $Y$, giving $2X + V \to 2X + Y + V$ with rate $3\alpha$. And so on for every term.

The result is nine reactions organized into three star-shaped trees, one for each source complex:

The canonic mechanism for the four-species ODE. Each source complex fans out into separate reactions for each species it affects. Nine reactions, twelve complexes.

Nine reactions, twelve complexes. The construction is mechanical — no thinking required.

A Minimal Realization

But as Hárs and Tóth themselves observed, the same ODE is generated by just three reactions forming a cycle:

A minimal realization: $V + 2X \xrightarrow{\alpha} 3Y + 2V \xrightarrow{\beta} 4Z \xrightarrow{\gamma} V + 2X$. Three reactions, three complexes.

Under mass-action kinetics, the reaction $V + 2X \xrightarrow{\alpha} 3Y + 2V$ contributes a stoichiometric change of $(-2, +3, 0, +1)$ at rate $\alpha\, x^2 v$. The reaction $3Y + 2V \xrightarrow{\beta} 4Z$ contributes $(0, -3, +4, -2)$ at rate $\beta\, y^3 v^2$. The reaction $4Z \xrightarrow{\gamma} V + 2X$ contributes $(+2, 0, -4, +1)$ at rate $\gamma\, z^4$. Summing contributions to each species gives the original ODE.

Three reactions instead of nine. Three complexes instead of twelve. What happened?

The canonic construction treats every positive and negative contribution independently. It never notices that “losing two $X$ molecules” and “gaining three $Y$ molecules” and “gaining one $V$ molecule” can all be the same event. The minimal realization merges them: a single reaction $V + 2X \to 3Y + 2V$ simultaneously accounts for contributions to $\dot{x}$, $\dot{y}$, and $\dot{v}$.

Why the Trivial Case Is Trivial

The canonic mechanism works by dumping all the positive contributions and all the negative contributions into separate reactions. There is no constraint on which positive terms match with which negative terms — every monomial gets its own reaction, and the stoichiometric bookkeeping is handled by construction.

This is the key observation: when there are no constraints on how terms pair up, the realization problem is trivial. You read the ODE, you write down the reactions. No algorithm needed.

When Constraints Make It Interesting

The minimal realization problem becomes genuinely hard when there are structural constraints on the CRN. Here are three examples:

1. Minimizing reactions or complexes. Given a fixed set of candidate complexes, the set of all CRN realizations of a given ODE forms a polyhedral cone (the rate constants are nonneg variables satisfying linear stoichiometric constraints). Finding the realization with the fewest reactions is a mixed-integer linear program [2]. Finding the realization with the fewest complexes adds a second layer of integer variables.

2. Weak reversibility. A CRN is weakly reversible if every connected component of the reaction graph is strongly connected — for every reaction $A \to B$, there is a path back from $B$ to $A$. Weakly reversible CRNs with deficiency zero are guaranteed to have a unique globally stable equilibrium (by the Deficiency Zero Theorem [3]). Craciun and Jin [4] gave a combinatorial algorithm for finding such realizations when they exist. The minimal realization above happens to be weakly reversible (it is a cycle), which is a nice bonus — but not all ODEs admit such clean realizations.

3. No self-production. In population protocol constructions [5], a species must not catalyze its own production. That is, if $X_k$ appears as a reactant, it should not also appear (with higher multiplicity) as a product in the same reaction. This “forbidden pairing” constraint means you cannot simply lump all positive and negative terms together — some matchings are prohibited. It is precisely this constraint that makes the PP-to-NAP construction in [5] require careful monomial splitting and a matching argument.

The Network Flow Perspective

There is a natural network flow interpretation. Complexes are nodes. Reactions are directed edges with nonneg capacities (rate constants). The mass-action constraint says the total stoichiometric effect of all reactions must equal the target polynomial vector field. This is a flow balance condition: the net species production/consumption at each complex node must sum to the right thing.

Minimizing the number of reactions is then a minimum-cost network flow problem: use the fewest edges to achieve the required flow balance. Minimizing the number of complexes means minimizing the number of nodes in the flow network.

From this viewpoint, the canonic mechanism is the trivial flow: one source-to-sink path for each demand unit. The minimal realization consolidates flows, sharing nodes and edges. The Hárs-Tóth example shows this consolidation at its most dramatic: three star-shaped trees collapse into a single triangle.

A Matching Theory Question

Here is a question that, as far as the literature indicates, has not been asked:

Consider a bipartite graph where one side is the set of monomials appearing in the ODE (the “demands”) and the other side is a candidate set of complexes (the “suppliers”). An edge connects a monomial to a complex if that complex can serve as a source or product in a reaction accounting for that monomial’s stoichiometric contribution. A valid realization requires covering all monomials.

Hall’s condition says: a complete matching exists if and only if every subset $S$ of monomials has at least $|S|$ neighboring complexes. When this condition fails, it gives a lower bound on the number of complexes needed.

The gap between this combinatorial lower bound and the number of complexes actually used is reminiscent of CRN deficiency in Feinberg’s sense ($\delta = |C| - \ell - s$, where $C$ is the complex set, $\ell$ the number of linkage classes, and $s$ the rank of the stoichiometric subspace). Whether there is a direct, clean connection between Hall-type deficiency and Feinberg deficiency remains open.

The moral: the CRN realization problem is only interesting when there are constraints. Without them, the canonic construction solves everything in one pass. The real mathematics begins when you ask for fewer reactions, fewer complexes, weak reversibility, no self-production, or other structural properties. Then matching theory, network flow, and integer programming all enter the picture — and the simple one-reaction-per-monomial approach is no longer good enough.

鸳鸯绣了从教看，莫把金针度与人。

Show the embroidered mandarin ducks, but do not hand out the golden needle.

— 元好问《论诗三十首》

References

1. V. Hárs, J. Tóth, “On the Inverse Problem of Reaction Kinetics,” in Qualitative Theory of Differential Equations, Colloq. Math. Soc. János Bolyai, vol. 30, 1981. [Recompiled PDF]
2. G. Szederkényi, “Computing sparse and dense realizations of reaction kinetic systems,” J. Math. Chem. 47, 2010; G. Szederkényi, K. M. Hangos, T. Péni, “Maximal and minimal realizations of reaction kinetic systems,” MATCH 65(2), 2011.
3. M. Feinberg, Foundations of Chemical Reaction Network Theory, Applied Mathematical Sciences vol. 202, Springer, 2019.
4. G. Craciun, W. Jin, “An algorithm for finding weakly reversible deficiency zero realizations of polynomial dynamical systems,” J. Nonlinear Sci., 2025.
5. X. Huang, T. H. Klinge, J. I. Lathrop, “Computing real numbers with large-population protocols,” DNA 28, 2022.