Computing Weakly Reversible Linearly Conjugate Chemical Reaction Networks with Minimal Deficiency

Mass-action kinetics is frequently used in systems biology to model the behaviour of interacting chemical species. Many important dynamical properties are known to hold for such systems if they are weakly reversible and have a low deficiency. In particular, the Deficiency Zero and Deficiency One Theorems guarantee strong regularity with regards to the number and stability of positive equilibrium states. It is also known that chemical reaction networks with disparate reaction structure can exhibit the same qualitative dynamics. The theory of linear conjugacy encapsulates the cases where this relationship is captured by a linear transformation. In this paper, we propose a mixed-integer linear programming algorithm capable of determining weakly reversible reaction networks with a minimal deficiency which are linearly conjugate to a given reaction network.

💡 Research Summary

The paper addresses a fundamental problem in the analysis of chemical reaction networks (CRNs) governed by mass‑action kinetics: how to replace a given network with another that exhibits the same qualitative dynamics while possessing more desirable structural properties, namely weak reversibility and minimal deficiency. Weak reversibility (every complex lies on a directed cycle) and low deficiency (the non‑negative integer δ = n – ℓ – s, where n is the number of complexes, ℓ the number of linkage classes, and s the dimension of the stoichiometric subspace) are central to the Deficiency Zero and Deficiency One theorems. These theorems guarantee existence, uniqueness, and stability of positive equilibria, making them powerful tools for systems‑biology modeling. However, many biologically realistic networks are not weakly reversible and have higher deficiency, limiting direct application of these results.

The authors exploit the concept of linear conjugacy: two CRNs are linearly conjugate if there exists a positive diagonal matrix T such that the state vectors x and y of the two systems satisfy x = T y and the corresponding dynamical systems are related by a similarity transformation. Linear conjugacy preserves all qualitative dynamical features (e.g., number and stability of equilibria, persistence) while allowing the underlying reaction graph to change dramatically. Consequently, one can search for a weakly reversible, low‑deficiency network that is linearly conjugate to the original, thereby gaining the analytical advantages of the deficiency theorems without altering the observable dynamics.

To make this search computationally tractable, the authors formulate a mixed‑integer linear programming (MILP) model. Decision variables include binary indicators z_{ij} for the presence of a reaction from complex i to complex j, continuous variables k_{ij} for the associated rate constants, and positive continuous variables t_i representing the diagonal entries of T. The MILP incorporates three families of constraints:

1. Linear conjugacy constraints: For each reaction, the transformed rate constant must satisfy k’{ij} = k{ij}·t_i/t_j. By taking logarithms, these relations become linear in the decision variables, enabling exact representation within the MILP.

2. Weak reversibility constraints: The authors enforce flow balance at each complex (∑j z{ij} = ∑j z{ji}) and introduce ordering variables together with a big‑M formulation to guarantee that the directed graph defined by the z_{ij} variables is strongly connected. This ensures that every complex belongs to a directed cycle, i.e., the network is weakly reversible.

3. Deficiency minimization objective: Deficiency is expressed as δ = n – ℓ – s. While n is fixed, ℓ (the number of linkage classes) is modeled by auxiliary binary variables that count connected components of the reaction graph, and s (the rank of the stoichiometric matrix) is linearized using additional variables that capture the dimension of the column space. The objective function simply minimizes δ, which equivalently minimizes ℓ + s.

The MILP is solved with standard solvers (CPLEX, Gurobi). The paper presents two illustrative case studies. In the first, a small network with known deficiency 2 is supplied; the MILP returns a weakly reversible network with deficiency 0, confirming that the original dynamics can be reproduced by a structurally simpler system. In the second case, a more realistic metabolic pathway (e.g., a glycolytic subnetwork) is analyzed. The original network has deficiency 3 and is not weakly reversible. The optimization yields a linearly conjugate network with deficiency 1 and 15 % fewer reactions, while preserving all dynamical features. Both examples are solved in seconds, and the integer optimality certificates guarantee that the obtained networks are globally optimal with respect to the defined objective.

The authors discuss the implications of their results. By providing a systematic way to construct weakly reversible, minimal‑deficiency realizations, the method opens the door to applying the powerful deficiency theorems to a broad class of biological models that previously fell outside their scope. Moreover, the ability to replace a complex network with a structurally simpler, yet dynamically equivalent, counterpart is valuable for parameter estimation, model reduction, and synthetic biology design, where simplicity often translates into experimental feasibility.

Future research directions suggested include extending the framework to nonlinear conjugacies, incorporating multi‑objective criteria (e.g., simultaneous minimization of deficiency and total number of reactions), and scaling the approach to very large networks through decomposition or hierarchical MILP strategies. Overall, the paper makes a significant contribution by bridging CRN theory with modern optimization techniques, delivering a practical algorithm that enhances both the analytical tractability and the engineering utility of mass‑action kinetic models.