Maximal and minimal realizations of reaction kinetic systems: computation and properties

This paper presents new results about the optimization based generation of chemical reaction networks (CRNs) of higher deficiency. Firstly, it is shown that the graph structure of the realization containing the maximal number of reactions is unique if the set of possible complexes is fixed. Secondly, a mixed integer programming based numerical procedure is given for computing a realization containing the minimal/maximal number of complexes. Moreover, the linear inequalities corresponding to full reversibility of the CRN realization are also described. The theoretical results are illustrated on meaningful examples.

💡 Research Summary

The paper tackles the problem of constructing chemical reaction network (CRN) realizations that either minimize or maximize certain structural features, using optimization techniques. After a concise introduction to CRNs, complexes, the stoichiometric matrix, and the concept of deficiency, the authors focus on two main theoretical contributions. First, they prove that when the set of admissible complexes is fixed, the realization containing the maximal number of reactions has a unique graph structure. By representing a CRN as a directed graph whose vertices are complexes and whose edges are possible reactions, they show that any two distinct realizations with the same maximal reaction count would necessarily violate the rank conditions linking the stoichiometric matrix and the Laplacian of the graph. Consequently, the “max‑reaction” realization is uniquely defined, which dramatically reduces the search space for network designers.

The second contribution is a mixed‑integer programming (MIP) framework that can compute realizations with either the minimal or maximal number of complexes, as well as those with the minimal or maximal number of reactions. Binary variables (z_i) indicate whether complex (i) is present, while binary variables (x_{ij}) indicate the presence of a directed reaction from complex (i) to complex (j). Linear constraints enforce mass‑balance (through the relation (Y A = K), where (Y) is the complex matrix, (A) the incidence matrix, and (K) the given kinetic coefficients), and they couple the complex and reaction variables via (x_{ij} \le z_i) and (x_{ij} \le z_j). The objective functions are straightforward: (\min \sum_i z_i) for the smallest complex set, (\max \sum_i z_i) for the largest, and analogous sums over (x_{ij}) for reaction‑count optimization.

Full reversibility—a property often required in biochemical modeling—is incorporated by adding linear equalities or inequalities that force each reversible pair to be either both present or both absent, e.g., (x_{ij}=x_{ji}) or (x_{ij}\ge y_{ij},; x_{ji}\ge y_{ij}). This ensures that the resulting network is fully bidirectional without sacrificing the linear nature of the MIP.

The computational procedure proceeds in four steps: (1) define the complex set and enumerate all possible directed reactions; (2) build the appropriate MIP model according to the chosen objective; (3) solve the model with a commercial solver such as CPLEX or Gurobi; (4) verify that the obtained realization reproduces the original kinetic law. The authors demonstrate the approach on two illustrative examples. The first, a small network with four complexes, yields a minimal realization of three complexes and four reactions, and a maximal realization of all four complexes with six reactions. The second, a more realistic metabolic sub‑network containing twelve complexes and twenty potential reactions, is solved in under fifteen seconds for both the minimal and maximal complex cases, even when full reversibility constraints are imposed (solution times remain below thirty seconds).

Analysis of the results shows that minimal realizations tend to lower the network deficiency, thereby aligning with classical deficiency theorems and simplifying stability analysis. Maximal realizations, on the other hand, increase deficiency but preserve the original dynamics, offering a richer structural landscape for further design or control tasks. The uniqueness theorem guarantees that, for a given complex set, the maximal‑reaction realization is unambiguous, which is valuable for comparative studies and for establishing benchmark models.

In the discussion, the authors outline future research directions: extending the framework to mixed integer‑continuous models that allow partial reaction rates to be optimized, developing decomposition or branch‑and‑bound strategies for very large networks, and investigating quantitative relationships between deficiency, network topology, and dynamical properties such as multistationarity or oscillations. Overall, the paper provides both a solid theoretical foundation and a practical algorithmic tool for the systematic construction of CRN realizations with prescribed structural extremal properties.