An efficient characterization of complex-balanced, detailed-balanced, and weakly reversible systems

TL;DR

This paper provides an efficient method to identify when polynomial or power-law dynamical systems originate from reaction networks with properties like complex balancing, detailed balancing, or weak reversibility, which are linked to stability and thermodynamic equilibrium.

Contribution

It introduces a computationally efficient characterization technique for recognizing systems derived from reaction networks with special properties such as complex balancing and weak reversibility.

Findings

Characterization method for complex-balanced systems

Identification of detailed-balanced systems

Recognition of weakly reversible systems

Abstract

Very often, models in biology, chemistry, physics, and engineering are systems of polynomial or power-law ordinary differential equations, arising from a reaction network. Such dynamical systems can be generated by many different reaction networks. On the other hand, networks with special properties (such as reversibility or weak reversibility) are known or conjectured to give rise to dynamical systems that have special properties: existence of positive steady states, persistence, permanence, and (for well-chosen parameters) complex balancing or detailed balancing. These last two are related to thermodynamic equilibrium, and therefore the positive steady states are unique and stable. We describe a computationally efficient characterization of polynomial or power-law dynamical systems that can be obtained as complex-balanced, detailed-balanced, weakly reversible, and reversible mass-action systems.

Equations

$$\frac{d\mathbf{x}}{dt}=\sum _{\mathbf{y}\to {\mathbf{y}}^{\prime }\in G}{k}_{\mathbf{y}\to {\mathbf{y}}^{\prime }}{\mathbf{x}}^{\mathbf{y}}({\mathbf{y}}^{\prime }-\mathbf{y}),$$

$${\mathbf{y}}_1\to {\mathbf{z}}_1$$

$${\mathbf{y}}_1=\left(\begin{array}{c}0\\ 0\end{array}\right),{\mathbf{y}}_2=\left(\begin{array}{c}0\\ 2\end{array}\right),{\mathbf{y}}_3=\left(\begin{array}{c}3\\ 2\end{array}\right),\text{and}{\mathbf{y}}_4=\left(\begin{array}{c}3\\ 0\end{array}\right).$$

$${\mathbf{y}}_5=\left(\begin{array}{c}1\\ 1\end{array}\right)\text{ and }{\mathbf{y}}_6=\left(\begin{array}{c}2\\ 1\end{array}\right).$$

$$\frac{d\mathbf{x}}{dt}=\sum _{\mathbf{y}\to {\mathbf{y}}^{\prime }\in G}{k}_{\mathbf{y}\to {\mathbf{y}}^{\prime }}{\mathbf{x}}^{\mathbf{y}}({\mathbf{y}}^{\prime }-\mathbf{y}),$$

$$\frac{d\mathbf{x}}{dt}$$

$$\sum _{{\mathbf{y}}_1\to {\mathbf{y}}_2\in G}{k}_{{\mathbf{y}}_1\to {\mathbf{y}}_2}{\mathbf{x}}^{{\mathbf{y}}_1}({\mathbf{y}}_2-{\mathbf{y}}_1)=\sum _{{\mathbf{y}}_1^{\prime }\to {\mathbf{y}}_2^{\prime }\in G^{\prime }}{k}_{{\mathbf{y}}_1^{\prime }\to {\mathbf{y}}_2^{\prime }}^{\prime }{\mathbf{x}}^{{\mathbf{y}}_1^{\prime }}({\mathbf{y}}_2^{\prime }-{\mathbf{y}}_1^{\prime })$$

$$\sum _{{\mathbf{y}}_0\to \mathbf{y}\in G}{k}_{{\mathbf{y}}_0\to \mathbf{y}}(\mathbf{y}-{\mathbf{y}}_0)=\sum _{{\mathbf{y}}_0\to {\mathbf{y}}^{\prime }\in G^{\prime }}{k}_{{\mathbf{y}}_0\to {\mathbf{y}}^{\prime }}^{\prime }({\mathbf{y}}^{\prime }-{\mathbf{y}}_0)$$

$$k_{15}\left(\begin{array}{c}1\\ 1\end{array}\right)$$

$$k_{25}\left(\begin{array}{c}1\\ -1\end{array}\right)$$

$$k_{36}\left(\begin{array}{c}-1\\ -1\end{array}\right)$$

$$k_{46}\left(\begin{array}{c}-1\\ 1\end{array}\right)$$

$$\frac{d\mathbf{x}}{dt}=\sum _{\mathbf{y}\to {\mathbf{y}}^{\prime }\in G}{k}_{\mathbf{y}\to {\mathbf{y}}^{\prime }}{\mathbf{x}}^{\mathbf{y}}({\mathbf{y}}^{\prime }-\mathbf{y}).$$

$$\frac{d\mathbf{x}}{dt}=\sum _{\mathbf{y}\to {\mathbf{y}}^{\prime }\in G}{k}_{\mathbf{y}\to {\mathbf{y}}^{\prime }}{\mathbf{x}}_0^{\mathbf{y}}({\mathbf{y}}^{\prime }-\mathbf{y})=0.$$

$$k_{\mathbf{y}\to {\mathbf{y}}^{\prime }}{\mathbf{x}}_0^{\mathbf{y}}=k_{{\mathbf{y}}^{\prime }\to \mathbf{y}}{\mathbf{x}}_0^{{\mathbf{y}}^{\prime }}.$$

$$\sum _{{\mathbf{y}}_0\to {\mathbf{y}}^{\prime }\in G}{k}_{{\mathbf{y}}_0\to {\mathbf{y}}^{\prime }}{\mathbf{x}}_0^{{\mathbf{y}}_0}=\sum _{\mathbf{y}\to {\mathbf{y}}_0\in G}{k}_{\mathbf{y}\to {\mathbf{y}}_0}{\mathbf{x}}_0^{\mathbf{y}}.$$

$$\frac{d\mathbf{x}}{dt}=\sum _{\mathbf{y}\to {\mathbf{y}}^{\prime }\in G}{\nu }_{\mathbf{y}\to {\mathbf{y}}^{\prime }}(\mathbf{x})({\mathbf{y}}^{\prime }-\mathbf{y}),$$

$$\sum _{\mathbf{y}\to {\mathbf{y}}^{\prime }\in G}{J}_{\mathbf{y}\to {\mathbf{y}}^{\prime }}({\mathbf{y}}^{\prime }-\mathbf{y})=0.$$

$$J_{\mathbf{y}\to {\mathbf{y}}^{\prime }}=J_{{\mathbf{y}}^{\prime }\to \mathbf{y}}.$$

$$\sum _{{\mathbf{y}}_0\to {\mathbf{y}}^{\prime }\in G}{J}_{{\mathbf{y}}_0\to {\mathbf{y}}^{\prime }}=\sum _{\mathbf{y}\to {\mathbf{y}}_0\in G}{J}_{\mathbf{y}\to {\mathbf{y}}_0}.$$

$$\displaystyle\begin{array}[]{rrr}k^{\prime}_{\footnotesize\chemfig{0}\to\chemfig{Y}}=3,&\qquad k^{\prime}_{\footnotesize\chemfig{Y}\to\chemfig{X}+\chemfig{Y}}=1,&\qquad k^{\prime}_{\footnotesize\chemfig{X}+\chemfig{Y}\to\chemfig{0}}=1,\displaystyle\vphantom{\frac{}{}}\\ k^{\prime}_{\footnotesize\chemfig{Y}\to\chemfig{0}}=\frac{1}{2},&k^{\prime}_{\footnotesize\chemfig{X}+\chemfig{Y}\to 2\chemfig{X}}=\frac{5}{2},&k^{\prime}_{\footnotesize 2\chemfig{X}\to\chemfig{X}+\chemfig{Y}}=5,\displaystyle\vphantom{\frac{}{}}\end{array}$$

$$\sum _{{\mathbf{y}}_0\to \mathbf{y}\in G}{J}_{{\mathbf{y}}_0\to \mathbf{y}}(\mathbf{y}-{\mathbf{y}}_0)=\sum _{{\mathbf{y}}_0\to {\mathbf{y}}^{\prime }\in G^{\prime }}{J}_{{\mathbf{y}}_0\to {\mathbf{y}}^{\prime }}^{\prime }({\mathbf{y}}^{\prime }-{\mathbf{y}}_0).$$

$$\sum _{{\mathbf{y}}_0\to \mathbf{y}\in G}{J}_{{\mathbf{y}}_0\to \mathbf{y}}(\mathbf{y}-{\mathbf{y}}_0)=\sum _{{\mathbf{y}}_0\to \mathbf{y}\in G^{\prime }}{J}_{{\mathbf{y}}_0\to \mathbf{y}}^{\prime }(\mathbf{y}-{\mathbf{y}}_0)=\sum _{{\mathbf{y}}_0\to \mathbf{y}\in G^{\ast }}{J}_{{\mathbf{y}}_0\to \mathbf{y}}^{\ast }(\mathbf{y}-{\mathbf{y}}_0)$$

$$\sum _{{\mathbf{y}}_0\to \mathbf{y}\in G}{J}_{{\mathbf{y}}_0\to \mathbf{y}}(\mathbf{x})(\mathbf{y}-{\mathbf{y}}_0)-\sum _{{\mathbf{y}}_0\to {\mathbf{y}}^{\prime }\in G^{\prime }}{J}_{{\mathbf{y}}_0\to {\mathbf{y}}^{\prime }}^{\prime }(\mathbf{x})({\mathbf{y}}^{\prime }-{\mathbf{y}}_0)$$

$$=\sum _{{\mathbf{y}}_0\to \mathbf{y}\in G}{k}_{{\mathbf{y}}_0\to \mathbf{y}}{\mathbf{x}}^{{\mathbf{y}}_0}(\mathbf{y}-{\mathbf{y}}_0)-\sum _{{\mathbf{y}}_0\to {\mathbf{y}}^{\prime }\in G^{\prime }}{k}_{{\mathbf{y}}_0\to {\mathbf{y}}^{\prime }}^{\prime }{\mathbf{x}}^{{\mathbf{y}}_0}({\mathbf{y}}^{\prime }-{\mathbf{y}}_0)$$

$$=\frac{{\mathbf{x}}^{{\mathbf{y}}_0}}{{\mathbf{x}}_0^{{\mathbf{y}}_0}}\left(\sum _{{\mathbf{y}}_0\to \mathbf{y}\in G}{k}_{{\mathbf{y}}_0\to \mathbf{y}}{\mathbf{x}}_0^{{\mathbf{y}}_0}(\mathbf{y}-{\mathbf{y}}_0)-\sum _{{\mathbf{y}}_0\to {\mathbf{y}}^{\prime }\in G^{\prime }}{k}_{{\mathbf{y}}_0\to {\mathbf{y}}^{\prime }}^{\prime }{\mathbf{x}}_0^{{\mathbf{y}}_0}({\mathbf{y}}^{\prime }-{\mathbf{y}}_0)\right)$$

$$=\frac{{\mathbf{x}}^{{\mathbf{y}}_0}}{{\mathbf{x}}_0^{{\mathbf{y}}_0}}\left(\sum _{{\mathbf{y}}_0\to \mathbf{y}\in G}{J}_{{\mathbf{y}}_0\to \mathbf{y}}({\mathbf{x}}_0)(\mathbf{y}-{\mathbf{y}}_0)-\sum _{{\mathbf{y}}_0\to {\mathbf{y}}^{\prime }\in G^{\prime }}{J}_{{\mathbf{y}}_0\to {\mathbf{y}}^{\prime }}^{\prime }({\mathbf{x}}_0)({\mathbf{y}}^{\prime }-{\mathbf{y}}_0)\right)$$

$$=0.$$

$$k_{\mathbf{y}\to {\mathbf{y}}^{\prime }}^{\prime }=\frac{J_{\mathbf{y}\to {\mathbf{y}}^{\prime }}^{\prime }}{{\mathbf{x}}_0^{\mathbf{y}}}>0,$$

$$\sum _{{\mathbf{y}}_0\to {\mathbf{y}}^{\prime }\in G}{J}_{{\mathbf{y}}_0\to {\mathbf{y}}^{\prime }}({\mathbf{y}}^{\prime }-{\mathbf{y}}_0)=0,$$

$$P_{(G,\mathbf{J})}({\mathbf{y}}^{\ast })=\sum _{\mathbf{y}\to {\mathbf{y}}^{\ast }\in G}{J}_{\mathbf{y}\to {\mathbf{y}}^{\ast }}-\sum _{{\mathbf{y}}^{\ast }\to {\mathbf{y}}^{\prime }\in G}{J}_{{\mathbf{y}}^{\ast }\to {\mathbf{y}}^{\prime }}.$$

$$\sum _{j=1}^M{J}_{\mathbf{z}\to {\mathbf{y}}_j}^{\prime }\left({\mathbf{y}}_j-\mathbf{z}\right)$$

$$-P_{(G^{\prime },{\mathbf{J}}^{\prime })}(\mathbf{z})$$

$$J_{{\mathbf{y}}_i\to {\mathbf{y}}_k}^{\prime }=J_{{\mathbf{y}}_i\to {\mathbf{y}}_k}+J_{{\mathbf{y}}^{\ast }\to {\mathbf{y}}_k}\frac{J_{{\mathbf{y}}_i\to {\mathbf{y}}^{\ast }}}{{\sum }_{{\mathbf{y}}_j\to {\mathbf{y}}^{\ast }\in G}{J}_{{\mathbf{y}}_j\to {\mathbf{y}}^{\ast }}}$$

Full text

1 Introduction

Many mathematical models in biology, chemistry, physics, and engineering are obtained from nonlinear interactions between several species or populations, such as (bio)chemical reactions in a cell or a chemical reactor, population dynamics in an ecosystem, or kinetic interactions in a gas or solution [GunaNts, HJ72, CasianGAC1, TDS, ctf06, CraciunBinh, Banaji_2013, Fein72, FeinLectNts, Feinberg_1987, Savageau_Voit, Sontag1, angeli_tutorial]. Very often, these models are generated by a graph of interactions according to specific kinetic rules; mass-action kinetics for reaction network models is one such example [RevCY].

If the graph underlying the mass-action system in a given reaction network has some special properties, then the associated dynamical system is known (or conjectured) to have certain dynamical properties. For example, dynamical systems generated by reversible reaction networks are known to have at least one positive steady state within each linear invariant subspace [Boros_2018]. Moreover, these models are known to be persistent and permanent if the number of species is small and are conjectured to have these properties for any number of species [CasianGAC2, CasianGAC1]. The same situation occurs for weakly reversible reaction networks, i.e., for networks where each reaction is part of a cycle (see Figure 2(b) and (c) for examples of such networks). For descriptions of other important classes of networks, see [ABCJ_2018].

Moreover, after some restrictions on the parameter values, weakly reversible networks give rise to complex-balanced systems, which are known to have a unique locally stable steady state within each linear invariant subspace. This steady state is known to be globally stable under some additional assumptions [CasianGAC1, DaveGAC, geomGAC, CasianGAC2] and is actually conjectured to be globally stable even without these assumptions [CasianGAC1, CraciunGAC]. If a reaction network is a complex-balanced system under mass-action kinetics, then other relevant models, ranging from continuous-time Markov chain models [Anderson_Craciun_Kurtz] to reaction-diffusion models [Desvillettes_2017, method_of_lines] and delay differential equation models [delayCB], are also stable in some sense.

It turns out that the same dynamical system can be generated by a multitude of reaction networks [cpIdentifiability, Sze2011, WRAlgDense, WRAlgSparse, HarsToth1987]. Therefore, if a system is generated by a network that does not enjoy a specific graphical property (e.g., not weakly reversible), we can ask whether the same system may be generated by a weakly reversible network. Others have asked this question before and formulated algorithms for a given number of complexes [PolylTimeAlgWR, Sze2011, WRAlgDense, WRAlgSparse] and applied the results to designing control systems [Sontag1, DEMASControl]. In order to determine whether a given system is generated by a weakly reversible or complex-balanced system, one would have to determine if it can be done using $n$ number of complexes for all $n\geq 1$.

In this paper we develop a theory of dynamical equivalence between mass-action systems (or more generally, polynomial or power-law dynamical systems) and weakly reversible and complex-balanced systems. Our results allow us to reformulate this dynamical equivalence problem as a linear feasibility problem whose dimension depends only on the size of the original system.

In order to describe our main results, we need to introduce some definitions and notations (these notions will be described in further detail in Section 2). For our purposes here, a reaction network is an oriented graph $G=(V_{G},E_{G})$ with vertex set $V_{G}$ and edge set $E_{G}$ such that $V_{G}\subseteq\mathbb{R}^{n}$. If ${\boldsymbol{y}}$, ${\boldsymbol{y}}^{\prime}\in V_{G}$ and $({\boldsymbol{y}},{\boldsymbol{y}}^{\prime})$ is an edge in $E_{G}\subseteq V_{G}\times V_{G}$, then we write ${\boldsymbol{y}}\to{\boldsymbol{y}}^{\prime}\in G$. With these notations, a dynamical system generated by $G$ (according to mass-action kinetics) is a system of ordinary differential equations on $\mathbb{R}_{>0}^{n}$ given by

[TABLE]

where ${\boldsymbol{x}}\in\mathbb{R}_{>0}^{n}$, ${\boldsymbol{x}}^{{\boldsymbol{y}}}=x_{1}^{y_{1}}x_{2}^{y_{2}}\cdots x_{n}^{y_{n}}$, and $k_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}}>0$ for all ${\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}\in G$. We will denote the dynamical system (1) by $G_{{\boldsymbol{k}}}$, where ${\boldsymbol{k}}$ is the vector of parameters $k_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}}$ for all ${\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}\in G$.

One of our main results is the following.

Theorem.

A mass-action system $G_{{\boldsymbol{k}}}$ is dynamically equivalent to some complex-balanced mass-action system if and only if it is dynamically equivalent to a complex-balanced mass-action system $G^{\prime}_{{\boldsymbol{k}}^{\prime}}$ that only uses the vertices of $G$, i.e., with $V_{G^{\prime}}\subseteq V_{G}$.

This theorem is useful not only for finding complex-balanced realizations of mass-action systems but also because for the first time, it gives us a computationally feasible way to decide if such realizations exist, as we only need to check if they exist for graphs $G^{\prime}$ that have $V_{G^{\prime}}\subseteq V_{G}$.

We will see in Section 4 that we can restrict the set $V_{G^{\prime}}$ even more: without loss of generality we can assume that it is contained in the set of “source vertices” of $G$. We have also obtained similar results for other important classes of mass-action systems: detailed-balanced, weakly reversible and reversible systems. Moreover, our results are shown for flux systems, which allows for other types of kinetics beside mass-action kinetics (Section 3).

Reaction networks and mass-action systems, along with all other relevant terms, are defined in Section 2. We view a reaction network as a directed graph embedded in Euclidean space. In Section 3, we define fluxes on a reaction network and relate them back to mass-action systems. Section 4 contains our main results for complex-balanced realizations, weakly reversible realizations, detailed-balanced realizations and reversible realizations. We make a brief comment on the implication of our results on the network’s deficiency. Finally, we present the relevant feasibility problems in Section 5.

2 Reaction Networks and Mass-Action Systems

Chemical reaction networks appear at the intersection of biology, biochemistry, chemistry, engineering, and mathematics. Different notations are used in the literature; here we explain the notations used throughout this paper. Introductions to chemical reaction network theory can be found in [RevCY, GunaNts, FeinLectNts].

Definition 2.1**.**

A reaction network (or simply a network) is a directed graph $G=(V_{G},E_{G})$ embedded in Euclidean space, with no self-loops, i.e., $V_{G}\subseteq\mathbb{R}^{n}$ and $E_{G}\subseteq V_{G}\times V_{G}$ and $({\boldsymbol{y}},{\boldsymbol{y}})\not\in E_{G}$ for any ${\boldsymbol{y}}\in V_{G}$.

When there is no ambiguity, we simply write $G=(V,E)$.

**Remark. **Vertices are points in $\mathbb{R}^{n}$, so an edge $e\in E$ can be regarded as a bona fide vector in $\mathbb{R}^{n}$. We denote an edge $e=({\boldsymbol{y}},{\boldsymbol{y}}^{\prime})$ as ${\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}$, which is associated to a reaction vector ${\boldsymbol{y}}^{\prime}-{\boldsymbol{y}}\in\mathbb{R}^{n}$. We also write ${\boldsymbol{y}}\to{\boldsymbol{y}}^{\prime}\in G$ instead of ${\boldsymbol{y}}\to{\boldsymbol{y}}^{\prime}\in E$.

The dimension $n$ of the ambient Euclidean space is the number of chemical species involved in the reaction network $G$. An edge in the set $E$ is called a reaction. A vertex in $V$ is also known as a reaction complex. The source vertex of a reaction ${\boldsymbol{y}}\to{\boldsymbol{y}}^{\prime}$ is the vertex ${\boldsymbol{y}}$, while ${\boldsymbol{y}}^{\prime}$ is the product vertex. Let $V_{s}\subseteq V$ denote the set of source vertices, i.e., the set of vertices that is the source of some reaction.

The vector space spanned by the reaction vectors is the stoichiometric subspace $S=\operatorname{span}_{\mathbb{R}}\{{\boldsymbol{y}}^{\prime}-{\boldsymbol{y}}\colon{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}\in G\}$. For any positive vector ${\boldsymbol{x}}_{0}\in\mathbb{R}_{>0}^{n}$, the affine polytope $({\boldsymbol{x}}_{0}+S)_{>}=({\boldsymbol{x}}_{0}+S)\cap\mathbb{R}_{>0}$ is known as the stoichiometric compatibility class of ${\boldsymbol{x}}_{0}$. A reaction network $G$ is reversible if ${\boldsymbol{y}}^{\prime}\rightarrow{\boldsymbol{y}}\in G$ whenever ${\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}\in G$; for simplicity, we denote such a pair of reactions by ${\boldsymbol{y}}\rightleftharpoons{\boldsymbol{y}}^{\prime}$. It is weakly reversible if every connected component of $G$ is strongly connected, i.e., every reaction ${\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}\in G$ is part of an oriented cycle.

Example 2.2**.**

Figure 1 shows a reaction network $G$ in $\mathbb{R}^{2}$ with $6$ vertices and $3$ reactions. The reactions are

[TABLE]

The stoichiometric subspace, which is the linear span of the reaction vectors, is $\mathbb{R}^{2}$. In particular, any stoichiometric compatibility class is all of $\mathbb{R}_{>0}^{2}$. The reaction network $G$ is neither reversible nor weakly reversible.

Example 2.3**.**

Three more examples of reaction networks are presented in Figure 2. The reaction networks (a) $G$, (b) $G^{\prime}$, and (c) $G^{*}$ share the vertices

[TABLE]

The reaction networks $G$, $G^{*}$ have two additional vertices

[TABLE]

The set of four reactions of $G$ is $E_{G}=\{{\boldsymbol{y}}_{1}\rightarrow{\boldsymbol{y}}_{5},\,{\boldsymbol{y}}_{2}\rightarrow{\boldsymbol{y}}_{5},\,{\boldsymbol{y}}_{3}\rightarrow{\boldsymbol{y}}_{6},\,{\boldsymbol{y}}_{4}\rightarrow{\boldsymbol{y}}_{6}\}$. The set of reactions of $G^{\prime}$ is $E_{G^{\prime}}=\{{\boldsymbol{y}}_{1}\rightleftharpoons{\boldsymbol{y}}_{2},\,{\boldsymbol{y}}_{2}\rightleftharpoons{\boldsymbol{y}}_{3},\,{\boldsymbol{y}}_{3}\rightleftharpoons{\boldsymbol{y}}_{4},\,{\boldsymbol{y}}_{4}\rightleftharpoons{\boldsymbol{y}}_{1},\,{\boldsymbol{y}}_{1}\rightleftharpoons{\boldsymbol{y}}_{3},\,{\boldsymbol{y}}_{2}\rightleftharpoons{\boldsymbol{y}}_{4}\}$. The set of reactions of $G^{*}$ is $E_{G^{*}}=\{{\boldsymbol{y}}_{1}\rightleftharpoons{\boldsymbol{y}}_{5}\rightleftharpoons{\boldsymbol{y}}_{2},\,{\boldsymbol{y}}_{3}\rightleftharpoons{\boldsymbol{y}}_{6}\rightleftharpoons{\boldsymbol{y}}_{4},\,{\boldsymbol{y}}_{5}\rightleftharpoons{\boldsymbol{y}}_{6},\,{\boldsymbol{y}}_{5}\rightarrow{\boldsymbol{y}}_{3},\,{\boldsymbol{y}}_{5}\rightarrow{\boldsymbol{y}}_{4}\}$. The networks $G^{\prime}$ and $G^{*}$ are weakly reversible, and $G^{\prime}$ is also reversible. The stoichiometric subspace is $S=\mathbb{R}^{2}$ for all three networks.

A reaction network $G$ is associated to a dynamical system, by assuming that each reaction ${\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}$ proceeds according to a rate function $\nu_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}}({\boldsymbol{x}})$, where ${\boldsymbol{x}}\in\mathbb{R}_{>0}^{n}$ is the vector of concentrations of the chemical species in the system. One of the most extensively studied kinetic systems is mass-action kinetics, where $\nu_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}}({\boldsymbol{x}})$ is a monomial whose exponent vector is ${\boldsymbol{y}}$.

Definition 2.4**.**

Let $G=(V,E)$ be a reaction network, and let ${\boldsymbol{k}}=(k_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}})_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}\in G}\in\mathbb{R}_{>0}^{E}$ be a vector of rate constants. We call the weighted directed graph $G_{{\boldsymbol{k}}}$ a mass-action system, whose associated dynamical system is the system on $\mathbb{R}_{>0}^{n}$

[TABLE]

where ${\boldsymbol{x}}^{{\boldsymbol{y}}}=x_{1}^{y_{1}}x_{2}^{y_{2}}\cdots x_{n}^{y_{n}}$. By convention, ${\boldsymbol{x}}^{{\boldsymbol{0}}}=1$.

It is convenient to refer to $k_{{\boldsymbol{y}}\to{\boldsymbol{y}}^{\prime}}$ even when ${\boldsymbol{y}}\to{\boldsymbol{y}}^{\prime}\not\in G$, in which case we mean $k_{{\boldsymbol{y}}\to{\boldsymbol{y}}^{\prime}}=0$. We adopt the convention that the empty sum is ${\boldsymbol{0}}$, i.e., $\sum_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}\in\varnothing}k_{{\boldsymbol{y}}\to{\boldsymbol{y}}^{\prime}}({\boldsymbol{y}}^{\prime}-{\boldsymbol{y}})={\boldsymbol{0}}$.

Example 2.5**.**

We revisit Example 2.2 under the assumption of mass-action kinetics. The dynamical system associated to this reaction network $G=(V,E)$ for an arbitrary vector of rate constants ${\boldsymbol{k}}=(k_{j})_{{\boldsymbol{y}}_{j}\rightarrow{\boldsymbol{z}}_{j}\in G}\in\mathbb{R}_{>0}^{E}$ is

[TABLE]

This is the Lotka–Volterra population dynamics model.

Given a mass-action system $G_{{\boldsymbol{k}}}$, (2) uniquely defines its associated dynamical system; however, many different reaction networks can give rise to the same dynamical system under mass-action kinetics. It has been known for a long time that if a reaction network has some special properties (e.g., reversible, weakly reversible, deficiency zero), then the mass-action system is known to have certain dynamical properties (e.g., existence of positive steady state, local and global stability). Therefore, given a mass-action system, we are interested in networks with richer structural properties that give rise to same dynamical systems. If two mass-action systems give rise to the same associated dynamical systems, we say they are dynamically equivalent [Sze2011, WRAlgDense, WRAlgSparse, cpIdentifiability].

Definition 2.6**.**

Two mass-action systems $G_{{\boldsymbol{k}}}$ and $G^{\prime}_{{\boldsymbol{k}}^{\prime}}$ are dynamically equivalent if

[TABLE]

for all ${\boldsymbol{x}}\in\mathbb{R}_{>0}^{n}$. We say that $G^{\prime}_{{\boldsymbol{k}}^{\prime}}$ is another realization of $G_{{\boldsymbol{k}}}$.

**Remark. **From (3), a necessary and sufficient condition for dynamical equivalence is

[TABLE]

for all ${\boldsymbol{y}}_{0}\in V_{G}\cup V_{G^{\prime}}$.111It is possible that either ${\boldsymbol{y}}_{0}\not\in V_{G}$ or ${\boldsymbol{y}}_{0}\not\in V_{G^{\prime}}$. Then one side of (4) is an empty sum, which by convention is ${\boldsymbol{0}}$.

Note that in the associated dynamical system of a mass-action system, $\frac{d{\boldsymbol{x}}}{dt}$ belongs to the stoichiometric subspace $S$. Moreover, $\mathbb{R}_{>0}^{n}$ is forward invariant under mass-action kinetics, i.e., if ${\boldsymbol{x}}(0)\in\mathbb{R}_{>0}^{n}$, then ${\boldsymbol{x}}(t)\in\mathbb{R}_{>0}^{n}$ for all $t\geq 0$ [FeinLectNts]. Consequently, the trajectory ${\boldsymbol{x}}(t)$ is confined to the stoichiometric compatibility class $({\boldsymbol{x}}(0)+S)_{>}$ for all $t\geq 0$.

**Remark. **The stoichiometric subspaces for dynamically equivalent systems can in principle be different. However, the kinetic subspaces for the two systems must be the same.222The kinetic subspace of a dynamical system $\frac{d{{\boldsymbol{x}}}}{dt}={{\boldsymbol{f}}({\boldsymbol{x}})}$ on a domain $\Omega$ is the linear subspace generated by $\{{{\boldsymbol{f}}({\boldsymbol{x}})\colon{\boldsymbol{x}}}\in\Omega\}$ [FH77]. For a mass-action system, the kinetic subspace is a subset of the stoichiometric subspace $S$. For example, the system in Figure 3(a), made of the reaction \ce2X -¿[$k$] X + Y, is dynamically equivalent to the system in Figure 3(b), consisting of the reactions \ce2X -¿[$k$] X + Y and \ce0 ¡-[$k^{\prime}$] Y -¿[$k^{\prime}$] 2Y. By definition, the two systems have different stoichiometric subspaces. However, in these systems, the trajectory starting at ${\boldsymbol{x}}_{0}\in\mathbb{R}_{>0}^{n}$ is confined to the affine space ${\boldsymbol{x}}_{0}+\mathbb{R}(-1,1)^{T}$ because their kinetic subspace is $\mathbb{R}(-1,1)^{T}$.

Example 2.7**.**

For the networks in Figure 2, let $k_{ij}>0$ be the rate constant on the reaction ${\boldsymbol{y}}_{i}\rightarrow{\boldsymbol{y}}_{j}\in G$; let $k^{\prime}_{ij}$ be the rate constant on the reaction ${\boldsymbol{y}}_{i}\rightarrow{\boldsymbol{y}}_{j}\in{G^{\prime}}$. Suppose $k_{ij}$ and $k_{pq}^{\prime}$ satisfy the following equations:

[TABLE]

Then $G_{{\boldsymbol{k}}}$ and $G^{\prime}_{{\boldsymbol{k}}^{\prime}}$ are dynamically equivalent. The linear constraints on the rate constants arise from vector decomposition of the reaction vectors starting at the source vertices of $G$ and $G^{\prime}$.

In fact, if ${\boldsymbol{k}}$, ${\boldsymbol{k}}^{\prime}$, and ${\boldsymbol{k}}^{*}$, where ${\boldsymbol{k}}^{*}$ is a vector of rate constants for $G^{*}$, satisfy some linear relations, the three mass-action systems $G_{{\boldsymbol{k}}}$, $G^{\prime}_{{\boldsymbol{k}}^{\prime}}$ and $G^{*}_{{\boldsymbol{k}}^{*}}$ are dynamically equivalent.

Mass-action systems give rise to very diverse dynamics. For example, weakly reversible deficiency zero mass-action systems have exactly one locally asymptotically stable steady state (within the same stoichiometric compatibility class). Yet there are other mass-action systems that have periodic orbits or limit cycles [osc1, osc2, osc3] and others that admit multiple steady states (within the same stoichiometric compatibility class) [cf05, cf06, Banaji_Craciun_2009], and even chaotic dynamics [RevCY, chaosphys]. We refer the reader to [RevCY, FeinLectNts, GunaNts, angeli_tutorial] for an introduction to mass-action systems. In this paper, we focus on several kinds of steady states of mass-action systems.

Definition 2.8**.**

Let $G_{{\boldsymbol{k}}}$ be a mass-action system with the associated dynamical system

[TABLE]

A state ${\boldsymbol{x}}_{0}\in\mathbb{R}_{>0}^{n}$ is a positive steady state if

[TABLE]

A positive steady state ${\boldsymbol{x}}_{0}\in\mathbb{R}_{>0}^{n}$ is detailed-balanced if for every ${\boldsymbol{y}}\rightleftharpoons{\boldsymbol{y}}^{\prime}\in G$, we have

[TABLE]

A positive steady state ${\boldsymbol{x}}_{0}\in\mathbb{R}_{>0}^{n}$ is complex-balanced if for every vertex ${\boldsymbol{y}}_{0}\in V_{G}$, we have

[TABLE]

Intuitively, detailed balancing is when fluxes across every pair of reversible reactions are balanced; this is intimately related to the notion of microreversibility or dynamical equilibrium in physical chemistry [Boltzmann_1887, Boltzmann_1896]. Complex balancing is when fluxes through every vertex (i.e., reaction complex) is balanced.

3 Fluxes on Reaction Networks

Most dynamical systems associated to reaction networks are nonlinear [HJ72, genMAS18, Savageau_Voit]. While nonlinear dynamical systems are generally difficult to study, the analysis of reaction networks is sometimes facilitated by the linear constraints arising from the network structure and stoichiometry.

To illustrate what we mean, consider mass-action kinetics. The (generally nonlinear) dynamical system under mass-action kinetics has the form

[TABLE]

where $\nu_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}}({\boldsymbol{x}})=k_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}}{\boldsymbol{x}}^{{\boldsymbol{y}}}$. Once the nonlinearity is hidden inside the reaction rate function $\nu_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}}({\boldsymbol{x}})$, the linear structure remaining becomes apparent.

Enumerate the set of reactions, $E=\{{\boldsymbol{y}}_{j}\rightarrow{\boldsymbol{y}}_{j}^{\prime}\}_{j=1}^{|E|}$, and let ${\boldsymbol{\nu}}({\boldsymbol{x}})=(\nu_{{\boldsymbol{y}}_{j}\rightarrow{\boldsymbol{y}}_{j}^{\prime}}({\boldsymbol{x}}))_{j=1}^{|E|}$ be a vector consisting of the reaction rate functions. Define the stoichiometric matrix $N\in\mathbb{R}^{n\times|E|}$ as the matrix whose $j$th column is the $j$th reaction vector ${\boldsymbol{y}}^{\prime}_{j}-{\boldsymbol{y}}_{j}$. Then the dynamical system above can be written succinctly as $\frac{d{\boldsymbol{x}}}{dt}=N{\boldsymbol{\nu}}({\boldsymbol{x}})$.

In order to deal with the underlying linear structure, we do not keep track of the concentrations that give rise to ${\boldsymbol{\nu}}({\boldsymbol{x}})$ but leave it as a vector of unknowns. For this reason, we denote the value ${\boldsymbol{\nu}}({\boldsymbol{x}})$ simply as ${\boldsymbol{J}}$ and call it a flux vector.

Definition 3.1**.**

A flux vector ${\boldsymbol{J}}=(J_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}})_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}\in G}\in\mathbb{R}_{>0}^{E}$ on a reaction network $G=(V,E)$ is a vector of positive numbers. The number $J_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}}$ is called the flux of the reaction ${\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}$, and the pair $(G,{\boldsymbol{J}})$ is called a flux system.333 The word “system” in “flux system” is in the sense of a system of linear equations, rather than a dynamical system.

As with the rate constants, it may be convenient to refer to $J_{{\boldsymbol{y}}\to{\boldsymbol{y}}^{\prime}}$ even when ${\boldsymbol{y}}\to{\boldsymbol{y}}^{\prime}\not\in G$, in which case $J_{{\boldsymbol{y}}\to{\boldsymbol{y}}^{\prime}}=0$.

This idea of fluxes on a reaction network may be familiar to anyone who has worked with stoichiometric network analysis or flux balance analysis. One form of the analysis is to solve the linear equation $N{\boldsymbol{J}}={\boldsymbol{0}}$, where the unknown vector ${\boldsymbol{J}}$ has nonnegative coordinates [FBARev, FBARev2]. Since we are interested in relating network structure with dynamics, if ${\boldsymbol{y}}\to{\boldsymbol{y}}^{\prime}\in G$, we impose that $J_{{\boldsymbol{y}}\to{\boldsymbol{y}}^{\prime}}>0$. Also if ${\boldsymbol{y}}\rightleftharpoons{\boldsymbol{y}}^{\prime}$ is a reversible reaction in $G$, then $J_{{\boldsymbol{y}}\to{\boldsymbol{y}}^{\prime}}$ and $J_{{\boldsymbol{y}}^{\prime}\to{\boldsymbol{y}}}$ are two positive components of the vector ${\boldsymbol{J}}$. A solution ${\boldsymbol{J}}>{\boldsymbol{0}}$ of the equation $N{\boldsymbol{J}}={\boldsymbol{0}}$ corresponds to a positive steady state if ${\boldsymbol{J}}={\boldsymbol{\nu}}({\boldsymbol{x}}_{0})$ for some ${\boldsymbol{x}}_{0}\in\mathbb{R}_{>0}^{n}$. We define the flux analogues of positive steady state, detailed-balanced steady state, and complex-balanced steady state.

Definition 3.2**.**

A steady state flux on a network $G=(V,E)$ is a flux vector ${\boldsymbol{J}}\in\mathbb{R}_{>0}^{E}$ satisfying

[TABLE]

A flux ${\boldsymbol{J}}\in\mathbb{R}_{>0}^{E}$ is said to be detailed-balanced if for every ${\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}\in G$, we have

[TABLE]

A flux ${\boldsymbol{J}}\in\mathbb{R}_{>0}^{E}$ is said to be complex-balanced if for every ${\boldsymbol{y}}_{0}\in V$, we have

[TABLE]

A steady state flux is a positive vector ${\boldsymbol{J}}$ in $\ker N$, where the stoichiometric matrix $N$ has the reaction vectors as its columns. As a shorthand, we refer to the flux system $(G,{\boldsymbol{J}})$ as detailed-balanced if ${\boldsymbol{J}}$ is a detailed-balanced flux on $G$. Similarly defined is a complex-balanced flux system on $G$. It will be clear from context whether a complex-balanced system refers to a mass-action system or a flux system.

Example 3.3**.**

An example of a flux system $(G,{\boldsymbol{J}})$ is shown in Figure 4. The positive number labelled on each edge ${\boldsymbol{y}}\to{\boldsymbol{y}}^{\prime}$ is the flux $J_{{\boldsymbol{y}}\to{\boldsymbol{y}}^{\prime}}$ of that reaction.

Note that this flux system could have risen from a mass-action system. For example, suppose the numbers labelled on the edges are taken to be rate constants ${\boldsymbol{k}}$, and the state of the system is ${\boldsymbol{x}}={\boldsymbol{1}}$. Then $(G,{\boldsymbol{J}})$ would be the flux system based off of the mass-action system $G_{{\boldsymbol{k}}}$.

There is no unique mass-action system that gives rise to a fixed flux system. For example, on the reaction network shown in Figure 4, suppose that the rate constants are taken to be

[TABLE]

and that the state of the system is ${\boldsymbol{x}}_{0}=(1,2)^{T}$; then it can be shown that $(G,{\boldsymbol{J}})$ is the flux system of the mass-action system $G_{{\boldsymbol{k}}^{\prime}}$ at the state ${\boldsymbol{x}}_{0}$.

This flux system $(G,{\boldsymbol{J}})$ is complex-balanced. For example, at the vertex $(0,1)$ corresponding to $\chemfig{Y}$, there is one reaction going into it with flux value $3$, and there are two reactions leaving this vertex, with sum of fluxes being $2+1=3$.

Whenever a flux vector arises from mass-action kinetics, i.e., $J_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}}=k_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}}{\boldsymbol{x}}^{{\boldsymbol{y}}}$, classical results for mass-action systems carry over to flux systems, as summarized in the following two lemmas.

Lemma 3.4**.**

Let $G_{\boldsymbol{k}}$ be a mass-action system, and fix ${\boldsymbol{x}}\in\mathbb{R}_{>0}^{n}$. For each edge ${\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}\in G$, define $J_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}}=k_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}}{\boldsymbol{x}}^{{\boldsymbol{y}}}$, so that ${\boldsymbol{J}}=(J_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}})_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}\in G}$ is a flux vector on the network $G$. The following hold:

The flux vector ${\boldsymbol{J}}$ is a steady state flux on $G$ if and only if ${\boldsymbol{x}}$ is a positive steady state of $G_{\boldsymbol{k}}$. 2. 2.

The flux vector ${\boldsymbol{J}}$ is detailed-balanced if and only if ${\boldsymbol{x}}$ is a detailed-balanced steady state for $G_{\boldsymbol{k}}$. 3. 3.

The flux vector ${\boldsymbol{J}}$ is complex-balanced if and only if ${\boldsymbol{x}}$ is a complex-balanced steady state for $G_{\boldsymbol{k}}$.

Lemma 3.5**.**

If $G$ admits a detailed-balanced flux, then $G$ is reversible; if $G$ admits a complex-balanced flux, then $G$ is weakly reversible. If a flux is detailed-balanced on $G$, then it is also complex-balanced; if a flux is complex-balanced, then it is also a steady state flux.

Proof.

Let ${\boldsymbol{J}}$ be a flux vector on a network $G$ — either detailed-balanced or complex-balanced or merely a steady state flux. On $G$, define a mass-action system $G_{\boldsymbol{k}}$ with rate constants $k_{{\boldsymbol{y}}\to{\boldsymbol{y}}^{\prime}}=J_{{\boldsymbol{y}}\to{\boldsymbol{y}}^{\prime}}$ for each ${\boldsymbol{y}}\to{\boldsymbol{y}}^{\prime}\in G$. Then ${\boldsymbol{x}}_{0}=(1,\cdots,1)^{T}$ is a (detailed-balanced or complex-balanced or positive) steady state.

Lemma 3.5 follows from classical results on mass-action systems [GunaNts, FeinLectNts, Feinberg_1987, Fein72, Horn72, Horn_1974]. ∎

As we have seen in the previous section, some mass-action systems are dynamically equivalent; similarly there are flux equivalent systems. We define an equivalence relation for flux systems in $\mathbb{R}^{n}$.

Definition 3.6**.**

Two flux systems $(G,{\boldsymbol{J}})$ and $(G^{\prime},{\boldsymbol{J}}^{\prime})$ are flux equivalent if for every vertex ${\boldsymbol{y}}_{0}\in V_{G}\cup V_{G^{\prime}}$,444As before, we adopt the convention that the empty sum is ${\boldsymbol{0}}$, i.e., $\sum_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}\in\varnothing}J_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}}({\boldsymbol{y}}^{\prime}-{\boldsymbol{y}})={\boldsymbol{0}}$. we have

[TABLE]

We denote equivalent flux systems by $(G,{\boldsymbol{J}})\sim(G^{\prime},{\boldsymbol{J}}^{\prime})$ and say that $(G^{\prime},{\boldsymbol{J}}^{\prime})$ is a realization of $(G,{\boldsymbol{J}})$.

Lemma 3.7**.**

Flux equivalence is an equivalence relation.

Proof.

That flux equivalence is symmetric and reflexive is clear. Suppose $(G,{\boldsymbol{J}})\sim(G^{\prime},{\boldsymbol{J}}^{\prime})$ and $(G^{\prime},{\boldsymbol{J}}^{\prime})\sim(G^{*},{\boldsymbol{J}}^{*})$. Transitivity follows from

[TABLE]

for any ${\boldsymbol{y}}_{0}\in V_{G}\cup V_{G^{\prime}}\cup V_{G^{*}}$. Note that if ${\boldsymbol{y}}_{0}\not\in V_{G^{\prime}}$, then the sums above are all ${\boldsymbol{0}}$. ∎

Suppose a flux vector arises from a mass-action system; one expects the notion of dynamical equivalence to line up with that of flux equivalence.

Proposition 3.8**.**

Let $G_{\boldsymbol{k}}$, $G^{\prime}_{{\boldsymbol{k}}^{\prime}}$ be mass-action systems, and fix ${\boldsymbol{x}}\in\mathbb{R}_{>0}^{n}$. For each edge ${\boldsymbol{y}}\to{\boldsymbol{y}}^{\prime}\in G$, let $J_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}}=k_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}}{\boldsymbol{x}}^{{\boldsymbol{y}}}$, so that ${\boldsymbol{J}}({\boldsymbol{x}})=(J_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}})_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}\in G}$ is a flux vector on $G$. Similarly, define the flux vector ${\boldsymbol{J}}^{\prime}({\boldsymbol{x}})=(J^{\prime}_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}})_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}\in G^{\prime}}$ on $G^{\prime}$, where $J^{\prime}_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}}=k^{\prime}_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}}{\boldsymbol{x}}^{{\boldsymbol{y}}}$. Then the following are equivalent:

The mass-action systems $G_{\boldsymbol{k}}$ and $G^{\prime}_{{\boldsymbol{k}}^{\prime}}$ are dynamically equivalent. 2. 2.

The flux systems $(G,{\boldsymbol{J}}({\boldsymbol{x}}))$, $(G^{\prime},{\boldsymbol{J}}^{\prime}({\boldsymbol{x}}))$ are flux equivalent for all ${\boldsymbol{x}}\in\mathbb{R}_{>0}^{n}$. 3. 3.

The flux systems $(G,{\boldsymbol{J}}({\boldsymbol{x}}))$, $(G^{\prime},{\boldsymbol{J}}^{\prime}({\boldsymbol{x}}))$ are flux equivalent for some ${\boldsymbol{x}}\in\mathbb{R}_{>0}^{n}$.

Proof.

It is clear that statements 1 and 2 are equivalent, and that statement 2 implies statement 3. Showing the implication of statement 1 from statement 3 will complete the proof. Let ${\boldsymbol{x}}_{0}\in\mathbb{R}_{>0}^{n}$ be a vector such that $(G,{\boldsymbol{J}}({\boldsymbol{x}}_{0}))\sim(G^{\prime},{\boldsymbol{J}}^{\prime}({\boldsymbol{x}}_{0}))$. For any ${\boldsymbol{y}}_{0}\in V_{G}\cup V_{G^{\prime}}$ and arbitrary ${\boldsymbol{x}}\in\mathbb{R}_{>0}^{n}$, we have

[TABLE]

∎

**Remark. **The proof above holds for kinetics other than mass-action type. For each (source) vertex ${\boldsymbol{y}}\in V_{G}\cup V_{G^{\prime}}$, define a rate function $\nu_{{\boldsymbol{y}}}\colon\mathbb{R}_{>0}^{n}\to\mathbb{R}_{>0}$. Then the above proposition holds when the flux vectors are defined to be $J_{{\boldsymbol{y}}\to{\boldsymbol{y}}^{\prime}}=k_{{\boldsymbol{y}}\to{\boldsymbol{y}}^{\prime}}\nu_{{\boldsymbol{y}}}({\boldsymbol{x}})$ for each ${\boldsymbol{y}}\to{\boldsymbol{y}}^{\prime}\in G$, and $J^{\prime}_{{\boldsymbol{y}}\to{\boldsymbol{y}}^{\prime}}=k^{\prime}_{{\boldsymbol{y}}\to{\boldsymbol{y}}^{\prime}}\nu_{{\boldsymbol{y}}}({\boldsymbol{x}})$ for each ${\boldsymbol{y}}\to{\boldsymbol{y}}^{\prime}\in G^{\prime}$.

In the following proposition, we reduce a nonlinear problem about mass-action systems to a linear problem about flux systems. Instead of showing that a mass-action system is dynamically equivalent to a complex-balanced (or detailed-balanced) system, it suffices to show that an appropriately defined flux system is flux equivalent to a complex-balanced (or detailed-balanced) system.

Proposition 3.9**.**

Let $G_{{\boldsymbol{k}}}$ be a mass-action system, and let ${\boldsymbol{x}}_{0}\in\mathbb{R}_{>0}^{n}$. For each edge ${\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}\in G$, define $J_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}}=k_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}}{\boldsymbol{x}}_{0}^{{\boldsymbol{y}}}$, so that ${\boldsymbol{J}}=(J_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}})_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}\in G}$ is a flux vector on the network $G$. Suppose $(G,{\boldsymbol{J}})$ is flux equivalent to $(G^{\prime},{\boldsymbol{J}}^{\prime})$, where ${\boldsymbol{J}}^{\prime}$ is complex-balanced; then $G_{\boldsymbol{k}}$ is dynamically equivalent to a mass-action system $G^{\prime}_{{\boldsymbol{k}}^{\prime}}$, where ${\boldsymbol{x}}_{0}$ is a complex-balanced steady state for $G^{\prime}_{{\boldsymbol{k}}^{\prime}}$. Similarly, if $(G,{\boldsymbol{J}})$ is flux equivalent to a detailed-balanced flux system $(G^{\prime},{\boldsymbol{J}}^{\prime})$, then $G_{\boldsymbol{k}}$ is dynamically equivalent to a mass-action system $G^{\prime}_{{\boldsymbol{k}}^{\prime}}$, where ${\boldsymbol{x}}_{0}$ is a detailed-balanced steady state for $G^{\prime}_{{\boldsymbol{k}}^{\prime}}$.

Proof.

For each edge ${\boldsymbol{y}}\to{\boldsymbol{y}}^{\prime}\in G^{\prime}$, define its rate constant to be

[TABLE]

so that $G^{\prime}_{{\boldsymbol{k}}^{\prime}}$ is a mass-action system. By Proposition 3.8, the mass-action systems $G_{\boldsymbol{k}}$ and $G^{\prime}_{{\boldsymbol{k}}^{\prime}}$ are dynamically equivalent, and by Lemma 3.4, ${\boldsymbol{x}}_{0}$ is a complex-balanced steady state if ${\boldsymbol{J}}^{\prime}$ is a complex-balanced flux on $G^{\prime}$, and if ${\boldsymbol{J}}^{\prime}$ is detailed-balanced on $G^{\prime}$, then ${\boldsymbol{x}}_{0}$ is a detailed-balanced steady state. ∎

4 Complex balancing without additional vertices

The identification of possible network structures associated to a biochemical system, say, from experimental data, is closely related to identifying key players in the system (e.g., enzymes in metabolic networks, genes in genetic networks). While the general nonuniqueness implies that network identification may often be impossible, it may still be desirable to compute equivalent systems — whether that be dynamical equivalence or flux equivalence — in order to conclude that the system has better properties than first suspected, e.g., weak reversibility or complex balance. This problem is not new [Sze2011, cpIdentifiability].

In recent years, the engineering community has utilized properties of mass-action systems in novel ways to designing and analyzing control systems [angeli_tutorial, Sontag1, KineticFeedbackDesign, DEMASControl]. For example, the controllers can be added in such a way that the resulting system is a complex-balanced mass-action system; from this, one can conclude that the control system has a unique positive steady state and local stability [KineticFeedbackDesign, DEMASControl]. Moreover, very general results have been obtained on the stability of complex-balanced systems with delay [delayCB].

Thus, there is strong incentive for developing effective computational methods to find structurally better dynamically equivalent systems. One approach uses linear programming, but an objective function must be chosen. To reduce the search space, one can decide to search for a realization with the maximal and minimal number of edges [WRAlgDense, WRAlgSparse]. Nonetheless, the set of vertices to be included in the reaction network must be chosen ahead of time.

In the examples of Figure 3, the mass-action systems systems are dynamically equivalent, but one uses an additional source vertex, whose weighted vectors sum to zero. Intuition may say that additional vertices can only improve the chance to find a network with desirable properties, as additional parameters provide extra degrees of freedom. Even if that is the case, the question of computability arises. Even if by adding new vertices to the network, one can produce an equivalent complex-balanced system, there is no a priori bound on the number of new vertices needed. One cannot realistically add new vertices ad infinitum.

Fortunately, we prove that no additional vertices are needed in order to check if a given system admits complex-balanced realizations. Thus, to check whether or not a network can admit a complex-balanced realization becomes a finite calculation, one that can be done by searching through the admissible domain as done in linear programming. Although the motivation came from mass-action systems, we prove our results in the more general setting of flux systems.

Our approach is to show that any such additional vertices in the network can be removed without changing the properties desired, namely, complex-balanced or weak reversibility. Such additional vertices will be called virtual sources.

Definition 4.1**.**

A vertex ${\boldsymbol{y}}_{0}\in V_{s}$ is a virtual source of the flux system $(G,{\boldsymbol{J}})$ if

[TABLE]

where the sum is over all edges with ${\boldsymbol{y}}_{0}$ as its source.

If the flux system $(G,{\boldsymbol{J}})$ arises from a mass-action system, then ${\boldsymbol{y}}_{0}\in V_{s}$ is a virtual source if and only if the monomial ${\boldsymbol{x}}^{{\boldsymbol{y}}_{0}}$ does not appear555That is the monomial does not appear after simplifying. on the right-hand side of the associated dynamical system (2). For example, if we consider fluxes that arise from mass-action kinetics in the network in Figure 3(b), the vertex $\chemfig{Y}$ is a virtual source.

In this section, we prove that if a flux vector on a weakly reversible reaction network is complex-balanced and has a virtual source ${\boldsymbol{y}}^{*}$, then there is an equivalent complex-balanced flux system that does not involve ${\boldsymbol{y}}^{*}$ at all. In short, virtual sources ${\boldsymbol{y}}^{*}$ are not needed for complex balancing.

Just as an arbitrary concentration vector ${\boldsymbol{x}}\in\mathbb{R}_{>0}^{n}$ may not be a complex-balanced steady state for a weakly reversible mass-action system, so we may want to speak of fluxes that are not complex-balanced. To keep track of how far a flux vector is from being complex-balanced, we define the potential at a vertex to be the difference between incoming and outgoing fluxes.

Definition 4.2**.**

Let $G=(V,E)$ be a reaction network, and let ${\boldsymbol{J}}\in\mathbb{R}_{>0}^{E}$ be a flux vector on $G$. The potential at a vertex ${\boldsymbol{y}}^{*}\in V$ is the scalar quantity

[TABLE]

**Remark. **The flux vector ${\boldsymbol{J}}$ is complex-balanced on $G$ if and only if $P_{(G,{\boldsymbol{J}})}({\boldsymbol{y}})=0$ for all ${\boldsymbol{y}}\in V_{s}$.

By an abuse of notation, if ${\boldsymbol{y}}^{*}\not\in G$, we still refer to the potential $P_{(G,{\boldsymbol{J}})}({\boldsymbol{y}}^{*})$ by setting it to be $P_{(G,{\boldsymbol{J}})}({\boldsymbol{y}}^{*})=0$.

In showing that virtual sources are not needed for complex balancing, the idea is to redirect the fluxes flowing into a virtual source ${\boldsymbol{y}}^{*}$ to other vertices while maintaining flux equivalence. If we are doing nothing more than redirecting flow of fluxes, the potential at every vertex does not change; therefore, we preserve complex balancing for the resulting flux system. This type of construction appeared first in [KineticFeedbackDesign] to show that new monomials were not necessary in feedback design.

We have to simultaneously keep track of the potential at each vertex and flux equivalence. We illustrate the key idea of Lemma 4.3 in Figure 5.

Lemma 4.3**.**

Consider a reaction network $G$ consisting of the reactions ${\boldsymbol{z}}\to{\boldsymbol{y}}^{*}$ and ${\boldsymbol{y}}^{*}\rightarrow{\boldsymbol{y}}_{j}$ for $j=1,2,\dots,M$. Suppose ${\boldsymbol{y}}^{*}$ is a virtual source for a flux system $(G,{\boldsymbol{J}})$ and its potential is $P_{(G,{\boldsymbol{J}})}({\boldsymbol{y}}^{*})=0$. Then there exists a flux equivalent system $(G^{\prime},{\boldsymbol{J}}^{\prime})$ such that ${\boldsymbol{y}}^{*}\not\in V_{G^{\prime}}$, and the potential at each vertex is preserved, i.e., $P_{(G,{\boldsymbol{J}})}({\boldsymbol{y}}_{j})=P_{(G^{\prime},{\boldsymbol{J}}^{\prime})}({\boldsymbol{y}}_{j})$ for $1\leq j\leq M$ and $P_{(G,{\boldsymbol{J}})}({\boldsymbol{z}})=P_{(G^{\prime},{\boldsymbol{J}}^{\prime})}({\boldsymbol{z}})$.

The flux system $(G^{\prime},{\boldsymbol{J}}^{\prime})$ can be obtained constructively: remove the edges ${\boldsymbol{z}}\rightarrow{\boldsymbol{y}}^{*}$ and ${\boldsymbol{y}}^{*}\rightarrow{\boldsymbol{y}}_{j}$, and add the edges ${\boldsymbol{z}}\rightarrow{\boldsymbol{y}}_{j}$ with fluxes $J^{\prime}_{{\boldsymbol{z}}\rightarrow{\boldsymbol{y}}_{j}}=J_{{\boldsymbol{y}}^{*}\rightarrow{\boldsymbol{y}}_{j}}$.

Proof.

For $j=1,2,\dots,M$, let ${\boldsymbol{w}}_{j}={\boldsymbol{y}}_{j}-{\boldsymbol{y}}^{*}$ and ${\boldsymbol{w}}_{0}={\boldsymbol{z}}-{\boldsymbol{y}}^{*}$ denote the reaction vectors. First, remove the edges ${\boldsymbol{y}}^{*}\rightarrow{\boldsymbol{y}}_{j}$ coming out of ${\boldsymbol{y}}^{*}$. Because ${\boldsymbol{y}}^{*}$ is a virtual source, $\sum_{j=1}^{M}J_{{\boldsymbol{y}}^{*}\rightarrow{\boldsymbol{y}}_{j}}{\boldsymbol{w}}_{j}={\boldsymbol{0}}$, so the resulting flux system is still equivalent to the original. Note that in this new flux system, only ${\boldsymbol{z}}$ is a source vertex.

Next, we redirect the reaction ${\boldsymbol{z}}\rightarrow{\boldsymbol{y}}^{*}$. Instead of the reaction ${\boldsymbol{z}}\rightarrow{\boldsymbol{y}}^{*}$ with flux $J_{{\boldsymbol{z}}\rightarrow{\boldsymbol{y}}^{*}}$, we have $M$ reactions ${\boldsymbol{z}}\rightarrow{\boldsymbol{y}}_{j}$ with fluxes $J^{\prime}_{{\boldsymbol{z}}\rightarrow{\boldsymbol{y}}_{j}}=J_{{\boldsymbol{y}}^{*}\rightarrow{\boldsymbol{y}}_{j}}$. Let $(G^{\prime},{\boldsymbol{J}}^{\prime})$ denote this newest flux system.

Recall that flux equivalence means (11) holds at each vertex of $G$ and $G^{\prime}$. Here we only need to look at the vertex ${\boldsymbol{z}}$ to show that $(G^{\prime},{\boldsymbol{J}}^{\prime})\sim(G,{\boldsymbol{J}})$. Note that ${\boldsymbol{y}}_{j}-{\boldsymbol{z}}={\boldsymbol{w}}_{j}-{\boldsymbol{w}}_{0}$. From $P_{(G,{\boldsymbol{J}})}({\boldsymbol{y}}^{*})=0$, we also have $\sum_{j=1}^{M}J_{{\boldsymbol{y}}^{*}\rightarrow{\boldsymbol{y}}_{j}}=J_{{\boldsymbol{z}}\rightarrow{\boldsymbol{y}}^{*}}$. Thus, the weighted sum of vectors coming out of ${\boldsymbol{z}}$ is

[TABLE]

and $(G^{\prime},{\boldsymbol{J}}^{\prime})\sim(G,{\boldsymbol{J}})$.

Finally, we prove that the potentials are unchanged. Trivially $P_{(G,{\boldsymbol{J}})}({\boldsymbol{y}}^{*})=P_{(G^{\prime},{\boldsymbol{J}}^{\prime})}({\boldsymbol{y}}^{*})=0$. Also $P_{(G,{\boldsymbol{J}})}({\boldsymbol{y}}_{j})=J_{{\boldsymbol{y}}^{*}\rightarrow{\boldsymbol{y}}_{j}}=J^{\prime}_{{\boldsymbol{z}}\rightarrow{\boldsymbol{y}}_{j}}=P_{(G^{\prime},{\boldsymbol{J}}^{\prime})}({\boldsymbol{y}}_{j})$ for $j=1,2,\dots,M$. Last but not least,

[TABLE]

We have shown that the resulting flux system $(G^{\prime},{\boldsymbol{J}}^{\prime})$ is flux equivalent to the original flux system $(G,{\boldsymbol{J}})$, and the potential at each vertex is preserved. ∎

**Remark. **In Lemma 4.3, the source vertex ${\boldsymbol{z}}$ may not be distinct from ${\boldsymbol{y}}_{j}$.

We now arrive at our main technical theorem (Theorem 4.4), a generalization of Lemma 4.3. Here, the virtual source ${\boldsymbol{y}}^{*}$ may have multiple reactions coming into it and coming out of it. The proof will be an induction on the number of edges flowing into ${\boldsymbol{y}}^{*}$. At each step, we redirect a fraction of the fluxes flowing through ${\boldsymbol{y}}^{*}$ from one incoming edge.

Theorem 4.4**.**

Let $(G,{\boldsymbol{J}})$ be a complex-balanced flux system on reaction network $G=(V,E)$. Suppose that ${\boldsymbol{y}}^{*}\in V$ is a virtual source. Then there exists an equivalent complex-balanced flux system $(G^{\prime},{\boldsymbol{J}}^{\prime})$ with $V_{G^{\prime}}=V\setminus\{{\boldsymbol{y}}^{*}\}$. Moreover,

[TABLE]

for any ${\boldsymbol{y}}_{i}$ such that ${\boldsymbol{y}}_{i}\to{\boldsymbol{y}}^{*}\in G$ and any ${\boldsymbol{y}}_{k}$ such that ${\boldsymbol{y}}^{*}\to{\boldsymbol{y}}_{k}\in G$, and $J^{\prime}_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}}=J_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}}$ for all other edges ${\boldsymbol{y}}\to{\boldsymbol{y}}^{\prime}$.

Proof.

Let $N$ be the number of reactions with ${\boldsymbol{y}}^{*}$ as target, i.e., $N=|\{{\boldsymbol{z}}\rightarrow{\boldsymbol{y}}^{*}\in G\}|$. Enumerate the sources as ${\boldsymbol{z}}_{1},{\boldsymbol{z}}_{2},\dots,{\boldsymbol{z}}_{N}$. Let $M$ be the number of reactions with ${\boldsymbol{y}}^{*}$ as sources, i.e., $M=|\{{\boldsymbol{y}}^{*}\rightarrow{\boldsymbol{y}}\in G\}|$. Enumerate the targets as ${\boldsymbol{y}}_{1},{\boldsymbol{y}}_{2},\dots,{\boldsymbol{y}}_{M}$. Since ${\boldsymbol{y}}^{*}$ is a virtual source, it is in the relative interior of the convex hull of the targets ${\boldsymbol{y}}_{j}$. From complex balancing, we have $P_{(G,{\boldsymbol{J}})}({\boldsymbol{y}}^{*})=0$, or

[TABLE]

Let $\theta=\frac{J_{{\boldsymbol{z}}_{1}\rightarrow{\boldsymbol{y}}^{*}}}{\sum_{i}J_{{\boldsymbol{z}}_{i}\rightarrow{\boldsymbol{y}}^{*}}}$ be the fraction of flux to be redirected from ${\boldsymbol{z}}_{1}\rightarrow{\boldsymbol{y}}^{*}$. We apply the construction described in Lemma 4.3 to the incoming edge ${\boldsymbol{z}}_{1}\rightarrow{\boldsymbol{y}}^{*}$, and the outgoing edges ${\boldsymbol{y}}^{*}\rightarrow{\boldsymbol{y}}_{j}$ for $j=1,2,\dots,M$. Let $(G^{\prime},{\boldsymbol{J}}^{\prime})$ denote the flux system after the diversion. More precisely, $J^{\prime}_{{\boldsymbol{z}}_{1}\rightarrow{\boldsymbol{y}}^{*}}=0$,

[TABLE]

and the fluxes on all other edges unchanged from ${\boldsymbol{J}}$.

Checking for flux equivalence at ${\boldsymbol{z}}_{1}$ before and after the diversion, we see that

[TABLE]

At all other vertices, the net flux is unchanged.

In terms of potentials, at ${\boldsymbol{z}}_{1}$, we have

[TABLE]

At each ${\boldsymbol{y}}_{j}$:

[TABLE]

At ${\boldsymbol{y}}^{*}$:

[TABLE]

The new flux system $(G^{\prime},{\boldsymbol{J}}^{\prime})$ after diverting the flux from ${\boldsymbol{z}}_{1}\rightarrow{\boldsymbol{y}}^{*}$ is still complex-balanced, as the potential is unchanged from those of $(G,{\boldsymbol{J}})$. Moreover, $(G^{\prime},{\boldsymbol{J}}^{\prime})$ and $(G,{\boldsymbol{J}})$ are flux equivalent. In addition, at ${\boldsymbol{y}}^{*}$, we have

[TABLE]

i.e, ${\boldsymbol{y}}^{*}$ is a virtual source for $(G^{\prime},{\boldsymbol{J}}^{\prime})$.

Thus we have recovered all the hypotheses stated in the theorem. The only difference between $(G,{\boldsymbol{J}})$ and $(G^{\prime},{\boldsymbol{J}}^{\prime})$ is that $G^{\prime}$ contains $N-1=|\left\{{\boldsymbol{z}}\rightarrow{\boldsymbol{y}}^{*}\in G^{\prime}\right\}|$ reactions with ${\boldsymbol{y}}^{*}$ as target vertex. By induction on the number $|\left\{{\boldsymbol{z}}\rightarrow{\boldsymbol{y}}^{*}\in G^{\prime}\right\}|$, there exists a flux system $(G^{*},{\boldsymbol{J}}^{*})$ that is flux equivalent to $(G,{\boldsymbol{J}})$, and for which ${\boldsymbol{J}}^{*}$ is a complex-balanced flux on $G^{*}$. Finally, because $P_{(G^{*},{\boldsymbol{J}}^{*})}({\boldsymbol{y}}^{*})=0$, but there are no incoming reactions to ${\boldsymbol{y}}^{*}$, it follows that there are no outgoing reactions from ${\boldsymbol{y}}^{*}$, i.e., ${\boldsymbol{y}}^{*}\not\in V_{G^{*}}$. ∎

When does a flux system (or a reaction network) admit a complex-balanced realization? Theorem 4.4 implies that virtual sources do not need to be considered. Theorem 4.5 below is the basis behind several relevant numerical methods in Section 5 for determining if a flux system (or a reaction network) is equivalent to complex-balanced.

Theorem 4.5**.**

Let $(G,{\boldsymbol{J}})$ be a flux system, and $V_{G,s}$ its set of source vertices. Then $(G,{\boldsymbol{J}})$ is flux equivalent to some complex-balanced flux system if and only if $(G,{\boldsymbol{J}})$ is flux equivalent to some complex-balanced flux system $(G^{\prime},{\boldsymbol{J}}^{\prime})$ where $V_{G^{\prime}}\subseteq V_{G,s}$.

Proof.

One direction is trivial. To prove the other direction, suppose $(G,{\boldsymbol{J}})$ is a flux system that is flux equivalent to some complex-balanced flux system $(\widetilde{G},\widetilde{{\boldsymbol{J}}})$. If ${\boldsymbol{y}}^{*}\in V_{\widetilde{G}}\setminus V_{G,s}$, the set $\{{\boldsymbol{y}}^{*}\to{\boldsymbol{y}}\in G\}$ is empty; flux equivalent demands that

[TABLE]

Theorem 4.4 implies we can maintain flux equivalence and complex balance even after dropping the vertex ${\boldsymbol{y}}^{*}$ from $V_{\widetilde{G}}$. Repeating this process for all vertices not in $V_{G,s}$ ultimately implies that there is a complex-balanced flux system $(G^{\prime},{\boldsymbol{J}}^{\prime})$ such that $(G^{\prime},{\boldsymbol{J}}^{\prime})\sim(G,{\boldsymbol{J}})$ and in addition $V_{G^{\prime}}\subseteq V_{G,s}$. ∎

Theorem 4.6**.**

Let $G$ be a reaction network, and $V_{G,s}$ its set of source vertices. Then the following are equivalent:

1. (i)

There exists a flux vector ${\boldsymbol{J}}$ such that $(G,{\boldsymbol{J}})$ is flux equivalent to some complex-balanced flux system. 2. (ii)

There exists a flux vector ${\boldsymbol{J}}$ such that $(G,{\boldsymbol{J}})$ is flux equivalent to some complex-balanced flux system $(G^{\prime},{\boldsymbol{J}}^{\prime})$, where $V_{G^{\prime}}\subseteq V_{G,s}$.

Proof.

The proof follows immediately from Theorem 4.5. ∎

Theorem 4.7**.**

A mass-action system $G_{{\boldsymbol{k}}}$ is dynamically equivalent to some complex-balanced system if and only if it is dynamically equivalent to a complex-balanced system $G^{\prime}_{{\boldsymbol{k}}^{\prime}}$ that only uses the source vertices, i.e., $V_{G^{\prime}}\subseteq V_{G,s}$.

Proof.

This theorem follows from Proposition 3.8 and Theorem 4.5. Suppose $G_{\boldsymbol{k}}$ is dynamically equivalent to some complex-balanced mass-action system $\widetilde{G}_{\widetilde{{\boldsymbol{k}}}}$. Define the appropriate fluxes ${\boldsymbol{J}}$ on $G$ and $\widetilde{{\boldsymbol{J}}}$ on $\widetilde{G}$; by Proposition 3.8, the two flux systems are flux equivalent. Theorem 4.5 holds if and only if $(G,{\boldsymbol{J}})$ is flux equivalent to some complex-balanced flux system $(G^{\prime},{\boldsymbol{J}}^{\prime})$ where $V_{G^{\prime}}\subseteq V_{G,s}$. Define the appropriate mass-action system $G^{\prime}_{{\boldsymbol{k}}^{\prime}}$ (see Proposition 3.9); we have one direction of this theorem. The other direction is trivially true. ∎

All of our theorems thus far have been concerned with flux systems; in the case of mass-action systems, implicit in everything is the existence of a complex-balanced steady state. However, the idea of redirecting fluxes can be adapted to show the surprising result that weak reversibility can be accomplished (if at all) with no extra vertices.

Theorem 4.8**.**

A mass-action system $G_{{\boldsymbol{k}}}$ is dynamically equivalent to some weakly reversible mass-action system if and only if it is dynamically equivalent to a weakly reversible mass-action system $G^{\prime}_{{\boldsymbol{k}}^{\prime}}$ that only uses its source vertices, i.e., $V_{G^{\prime}}\subseteq V_{G,s}$.

Proof.

Without loss of generality, we may suppose that $G_{{\boldsymbol{k}}}$ is a weakly reversible mass-action system for which there exists a virtual source ${\boldsymbol{y}}^{*}$. As in Theorem 4.4, we remove the vertex ${\boldsymbol{y}}^{*}$ by redirecting the reactions flowing through it. Since $G$ is weakly reversible, there exists some vertex ${\boldsymbol{z}}$ such that ${\boldsymbol{z}}\to{\boldsymbol{y}}^{*}\in G$. As before, we will try to replace pairs of reactions ${\boldsymbol{z}}\to{\boldsymbol{y}}^{*}$ and ${\boldsymbol{y}}^{*}\to{\boldsymbol{y}}$ with ${\boldsymbol{z}}\to{\boldsymbol{y}}$.

Enumerate the set $\{{\boldsymbol{y}}^{*}\to{\boldsymbol{y}}\in G\}$ as $\{{\boldsymbol{y}}^{*}\to{\boldsymbol{y}}_{i}\}_{i=1}^{M}$, and enumerate the set $\{{\boldsymbol{z}}\to{\boldsymbol{y}}^{*}\in G\}$ as $\{{\boldsymbol{z}}_{j}\to{\boldsymbol{y}}^{*}\}_{j=1}^{N}$. For simplicity, let $\alpha_{j}=k_{{\boldsymbol{z}}_{j}\to{\boldsymbol{y}}^{*}}$, and let $\beta_{i}=k_{{\boldsymbol{y}}^{*}\to{\boldsymbol{y}}_{i}}$. Informally speaking, in place of the reactions ${\boldsymbol{z}}_{j}\to{\boldsymbol{y}}^{*}$ and ${\boldsymbol{y}}^{*}\to{\boldsymbol{y}}_{i}$, we shall have the reaction ${\boldsymbol{z}}_{j}\to{\boldsymbol{y}}_{i}$ with rate constant $k^{\prime}_{{\boldsymbol{z}}_{j}\to{\boldsymbol{y}}_{i}}=\alpha_{j}\frac{\beta_{i}}{\sum\beta_{s}}$. More precisely, let $G^{\prime}$ be the graph after deleting the vertex ${\boldsymbol{y}}^{*}$ and its adjacent edges from $G$, and (if needed) the edges ${\boldsymbol{z}}_{j}\to{\boldsymbol{y}}_{i}$ added for all $i=1,2,\dots,M$ and $j=1,2,\dots,N$. On $G^{\prime}$, take the rate constants to be $k^{\prime}_{{\boldsymbol{z}}_{j}\to{\boldsymbol{y}}^{*}}=k^{\prime}_{{\boldsymbol{y}}^{*}\to{\boldsymbol{y}}_{i}}=0$ and

[TABLE]

and all other rate constants same as in $G_{\boldsymbol{k}}$.

The assumption that ${\boldsymbol{y}}^{*}$ is a virtual source can be written as

[TABLE]

Now to check for dynamical equivalence at ${\boldsymbol{z}}_{1}$, we consider the differences due to the reactions ${\boldsymbol{z}}_{1}\to{\boldsymbol{y}}_{i}$:

[TABLE]

which is the contribution from the reaction ${\boldsymbol{z}}_{1}\to{\boldsymbol{y}}^{*}$. Since other reactions were untouched, we have dynamical equivalence at ${\boldsymbol{z}}_{1}$. There is nothing special about $j=1$; the same holds for all source vertices ${\boldsymbol{z}}_{2},{\boldsymbol{z}}_{3},\dots,{\boldsymbol{z}}_{N}$.

Finally, given any cycle ${\boldsymbol{v}}_{1}\to{\boldsymbol{v}}_{2}\to\cdots\to{\boldsymbol{v}}_{\ell}\to{\boldsymbol{v}}_{1}$ in $G^{\prime}$, whenever an edge ${\boldsymbol{z}}_{j}\to{\boldsymbol{y}}_{i}$ appears in the cycle, replace it with two edges ${\boldsymbol{z}}_{j}\to{\boldsymbol{y}}^{*}\to{\boldsymbol{y}}_{i}$, and obtain a cycle in $G$. Therefore, $G^{\prime}$ is still weakly reversible. ∎

We extend the above results (Theorems 4.4-4.8) to detailed-balanced fluxes and/or reversible networks. We summarize these results in the following theorems.

Theorem 4.9**.**

Let $(G,{\boldsymbol{J}})$ be a detailed-balanced flux system on a reaction network $G=(V,E)$. Suppose that ${\boldsymbol{y}}^{*}\in V$ is a virtual source. Then there exists an equivalent detailed-balanced flux system $(G^{\prime},{\boldsymbol{J}}^{\prime})$ with $V_{G^{\prime}}=V\setminus\{{\boldsymbol{y}}^{*}\}$. Moreover,

[TABLE]

for any ${\boldsymbol{y}}_{i},{\boldsymbol{y}}_{k}$ connected to ${\boldsymbol{y}}^{*}$ in $(G,{\boldsymbol{J}})$. Let other fluxes remain unchanged from $(G,{\boldsymbol{J}})$. In particular, $(G,{\boldsymbol{J}})$ is flux equivalent to some detailed-balanced flux system if and only if $(G,{\boldsymbol{J}})$ is flux equivalent to some detailed-balanced flux system $(G^{\prime},{\boldsymbol{J}}^{\prime})$ where $V_{G^{\prime}}\subseteq V_{G,s}$.

Proof.

As in Theorem 4.4, we divert fluxes away from ${\boldsymbol{y}}^{*}$. We only need to check detail balancing. Consider any two vertices ${\boldsymbol{y}}_{i}\neq{\boldsymbol{y}}_{k}$ where ${\boldsymbol{y}}_{i}\rightleftharpoons{\boldsymbol{y}}^{*},{\boldsymbol{y}}_{k}\rightleftharpoons{\boldsymbol{y}}^{*}\in G$. Using the fact that the flux system was originally detailed-balanced, i.e., $J_{{\boldsymbol{y}}\to{\boldsymbol{y}}^{\prime}}=J_{{\boldsymbol{y}}^{\prime}\to{\boldsymbol{y}}}$, we obtain

[TABLE]

For any other pairs of reversible reaction, detail balancing is inherited from $(G,{\boldsymbol{J}})$. In other words, $(G^{\prime},{\boldsymbol{J}}^{\prime})$ is detailed-balanced. ∎

Theorem 4.10**.**

A mass-action system $G_{{\boldsymbol{k}}}$ is dynamically equivalent to some reversible system if and only if it is dynamically equivalent to a reversible system $G^{\prime}_{{\boldsymbol{k}}^{\prime}}$ that only uses its source vertices, i.e., $V_{G^{\prime}}\subseteq V_{G,s}$.

Proof.

We assume that $G_{\boldsymbol{k}}$ is reversible and has a virtual source ${\boldsymbol{y}}^{*}\in V_{G}$. We will replace the reactions $\{{\boldsymbol{y}}^{*}\rightleftharpoons{\boldsymbol{y}}_{i}\in G\}$ by modifying/adding the reactions $\{{\boldsymbol{y}}_{i}\rightleftharpoons{\boldsymbol{y}}_{k}\colon{\boldsymbol{y}}_{i}\rightleftharpoons{\boldsymbol{y}}^{*},{\boldsymbol{y}}_{k}\rightleftharpoons{\boldsymbol{y}}^{*}\in G\}$. For any ${\boldsymbol{y}}_{i}$, ${\boldsymbol{y}}_{j}$ such that ${\boldsymbol{y}}_{i}\rightleftharpoons{\boldsymbol{y}}^{*}$, ${\boldsymbol{y}}_{k}\rightleftharpoons{\boldsymbol{y}}^{*}\in G$, let $k^{\prime}_{{\boldsymbol{y}}_{i}\to{\boldsymbol{y}}^{*}}=k^{\prime}_{{\boldsymbol{y}}^{*}\to{\boldsymbol{y}}_{i}}=0$ and

[TABLE]

Similar to Theorem 4.8, it can be shown that $G_{{\boldsymbol{k}}}$ and $G^{\prime}_{{\boldsymbol{k}}^{\prime}}$ are dynamically equivalent. Moreover, by symmetry of construction, $G^{\prime}$ is reversible. ∎

Note that related results have been obtained recently for the problem of kinetic feedback design involving complex-balanced and weakly reversible systems [KineticFeedbackDesign]. Here, for the problem of dynamical equivalence, we show that a given system admits a dynamically equivalent system that is complex-balanced (or weakly reversible, or detailed-balanced, or reversible) if and only if such a system exists using only the complexes that are already present in the original system.

4.1 Connection to deficiency theory

Within the reaction network theory literature, deficiency is a well-known quantity defined for a network $G$. Equipped with mass-action kinetics, networks with low deficiency are known to enjoy special dynamical properties under mass-action kinetics. For example, the famous deficiency zero theorem says that a weakly reversible deficiency zero network is complex-balanced for any choices of rate constants [HJ72, Feinberg_1987]. As we have introduced, complex-balanced systems enjoy properties such as uniqueness and stability of steady states, existence of a Lyapunov function, and the steady states admit a monomial parametrization [GunaNts, FeinLectNts, RevCY, HJ72, Fein95]. Despite the strong implications, deficiency has a relatively simple definition.

Definition 4.11**.**

Let $G=(V_{G},E_{G})$ be a reaction network with $\ell_{G}$ connected components. Suppose the dimension of the stoichiometric subspace $S$ is $s=\dim S$; then the deficiency of the network $G$ is the nonnegative integer

[TABLE]

It can be shown that $\delta_{G}=\dim(\operatorname{ker}Y\cap\operatorname{im}I_{G})$, where $Y$ is the stoichiometric matrix, with the vertices as its columns, and $I_{G}$ the incidence matrix of $G$ [Joh14]. It follows that $\delta_{G}$ is a nonnegative integer. When the network is weakly reversible, we also have $\delta_{G}=\dim(\operatorname{ker}Y\cap\operatorname{im}A_{{\boldsymbol{k}}})$, where $-A_{{\boldsymbol{k}}}^{T}$ is the Laplacian of the weighted graph $G_{{\boldsymbol{k}}}$ [GunaNts, FeinLectNts].

Deficiency continues to play an important role in the analysis of reaction networks and mass-action systems. In our procedure for removing virtual vertices, deficiency always decreases. This is similar to a result obtained in [KineticFeedbackDesign], where the removal of additional monomials that function as controls in a feedback system also leads to a decrease in deficiency.

Theorem 4.12**.**

Let $G_{{\boldsymbol{k}}}$ be a weakly reversible mass-action system with deficiency $\delta_{G}$. Suppose it has a virtual source ${\boldsymbol{y}}^{*}$. Let $G^{\prime}_{{\boldsymbol{k}}^{\prime}}$ be the weakly reversible mass-action system as produced in Theorem 4.8, dynamically equivalent to $G_{{\boldsymbol{k}}}$ with $V_{G^{\prime}}=V_{G}\setminus\{{\boldsymbol{y}}^{*}\}$. Then the deficiency of $G^{\prime}_{{\boldsymbol{k}}^{\prime}}$ is $\delta_{G^{\prime}}=\delta_{G}-1$.

Proof.

In the proof of Theorem 4.8, we replaced the reactions ${\boldsymbol{z}}\to{\boldsymbol{y}}^{*}$ and ${\boldsymbol{y}}^{*}\to{\boldsymbol{y}}$ with the reaction ${\boldsymbol{z}}\to{\boldsymbol{y}}$ by choosing appropriate rate constants. It is clear that $|V_{G^{\prime}}|=|V_{G}|-1$, and the number of linkage classes stays the same. We claim that the stoichiometric subspace $S$ remains unchanged. Thus, the drop in deficiency is due to the removal of the vertex ${\boldsymbol{y}}^{*}$, and $\delta_{G^{\prime}}=\delta_{G}-1$.

First enumerate the reactions coming out of ${\boldsymbol{y}}^{*}$ as ${\boldsymbol{y}}^{*}\to{\boldsymbol{y}}_{j}$, and enumerate the reactions going into ${\boldsymbol{y}}^{*}$ as ${\boldsymbol{z}}_{i}\to{\boldsymbol{y}}^{*}$. Let $S_{0}$ be the span of the reaction vectors “untouched” by our procedure, more precisely,

[TABLE]

Let $S_{G}$ be the stoichiometric subspace of $G$, in particular,

[TABLE]

and $S_{G^{\prime}}$ be the stoichiometric subspace of $G^{\prime}$, where

[TABLE]

Clearly, $S_{G^{\prime}}\subseteq S_{G}$, since ${\boldsymbol{y}}_{j}-{\boldsymbol{z}}_{i}=({\boldsymbol{y}}_{j}-{\boldsymbol{y}}^{*})+({\boldsymbol{y}}^{*}-{\boldsymbol{z}}_{i})\in S_{G}$. Moreover, because $G$ is weakly reversible, the edge ${\boldsymbol{y}}^{*}\to{\boldsymbol{y}}_{j}$ is a part of a cycle; therefore, $S_{G}=\operatorname{span}_{\mathbb{R}}\{S_{0},\,{\boldsymbol{y}}^{*}-{\boldsymbol{z}}_{i}\}_{i}$. Finally, we note that ${\boldsymbol{y}}^{*}$ is in the convex hull of the vertices ${\boldsymbol{y}}_{j}$, and thus ${\boldsymbol{y}}^{*}-{\boldsymbol{z}}_{i}\in\operatorname{span}_{\mathbb{R}}\{{\boldsymbol{y}}_{j}-{\boldsymbol{z}}_{i}\}_{j}$, which implies $S_{G}\subseteq S_{G^{\prime}}$. In other words, $S_{G}=S_{G^{\prime}}$ and $\delta_{G^{\prime}}=\delta_{G}-1$. ∎

5 Numerical methods

In this section, we characterize when a flux system or a mass-action system is equivalent to a complex-balanced system. We also describe a method to determine when a mass-action system is dynamically equivalent to a complex-balanced or weakly reversible system.

5.1 Flux equivalence to complex-balancing

Is a steady state flux system $(G,{\boldsymbol{J}})$ flux equivalent to a complex-balanced one? The answer lies in the following linear feasibility problem for an unknown vector ${\boldsymbol{J}}^{\prime}$. Enumerate the set of source vertices in $G$ as $\{{\boldsymbol{y}}_{1},{\boldsymbol{y}}_{2},\dots,{\boldsymbol{y}}_{N}\}$. Search for ${\boldsymbol{J}}^{\prime}=(J^{\prime}_{{\boldsymbol{y}}_{i}\to{\boldsymbol{y}}_{j}})_{i\neq j}\in\mathbb{R}^{N^{2}-N}$ satisfying

[TABLE]

If such a flux vector ${\boldsymbol{J}}^{\prime}$ exists, then $(G,{\boldsymbol{J}})$ is flux equivalent to a complex-balanced system. If no such flux vector ${\boldsymbol{J}}^{\prime}$ exists, then $(G,{\boldsymbol{J}})$ is not flux equivalent to a complex-balanced system.

Equation (17) is the flux equivalence condition, while (18a) ensures that the new flux system is complex-balanced. Equation (17) alone checks for flux equivalence between any two given systems $(G,{\boldsymbol{J}})$ and $(G^{\prime},{\boldsymbol{J}}^{\prime})$.

Example 5.1**.**

We return to the network $G$ in Figure 2(a) and Example 2.7. The network has $6$ vertices, $4$ of which are sources, and $4$ reactions. At the moment, we consider a flux system on the graph $G$ and ask, for what flux ${\boldsymbol{J}}$ is the flux system $(G,{\boldsymbol{J}})$ equivalent to a complex-balanced one? One can show that (17)-(19) hold if and only if

[TABLE]

A chosen flux ${\boldsymbol{J}}$ that satisfies (20) is flux equivalent to a complex-balanced system, whose network is a subgraph of $G^{\prime}$ of Figure 2(b). The details of this characterization will be in an upcoming paper [STN].

**Remark. **The setup for the detailed-balanced case is defined analogously. We keep (17) and (19) and include the equation

[TABLE]

5.2 Dynamical equivalence to complex balancing

We considered above a set of equalities and inequalities necessary and sufficient for a flux system to be equivalent to a complex-balanced one. If the flux system arises from mass-action kinetics, we can write down an analogous system of equalities and inequalities necessary and sufficient for dynamical equivalence to a complex-balanced system.

Consider a mass-action system $G_{{\boldsymbol{k}}}$, whose vertices are points in $\mathbb{R}^{n}$, and enumerate the set of source vertices in $G$ as $\{{\boldsymbol{y}}_{1},{\boldsymbol{y}}_{2},\dots,{\boldsymbol{y}}_{N}\}$. We set up a nonlinear feasibility problem for unknowns ${\boldsymbol{k}}^{\prime}$ and ${\boldsymbol{x}}$. Search for vectors ${\boldsymbol{k}}^{\prime}=(k^{\prime}_{{\boldsymbol{y}}_{i}\to{\boldsymbol{y}}_{j}})_{i\neq j}\in\mathbb{R}^{N^{2}-N}$ and ${\boldsymbol{x}}\in\mathbb{R}^{n}$ satisfying

[TABLE]

If such ${\boldsymbol{k}}^{\prime}$ and ${\boldsymbol{x}}$ exist, then $G_{{\boldsymbol{k}}}$ is dynamically equivalent to a complex-balanced system with ${\boldsymbol{x}}$ a complex-balanced steady state. If no such rate constants and steady state exist, then $G_{{\boldsymbol{k}}}$ is not dynamically equivalent to a complex-balanced system.

Equation (21) enforces dynamical equivalence. Equations (22) and (24) imply that ${\boldsymbol{x}}$ is a positive complex-balanced steady state for an equivalent mass-action system; hence ${\boldsymbol{x}}$ is a positive steady state of $G_{{\boldsymbol{k}}}$. Note that in the inequality (23), some $k^{\prime}_{{\boldsymbol{y}}_{i}\to{\boldsymbol{y}}_{j}}$ can be zero, which implies that ${\boldsymbol{y}}_{i}\to{\boldsymbol{y}}_{j}$ is not a reaction in the equivalent network.

Equations (21)-(24) generally form a nonlinear problem. Despite that, for networks with additional structure, one may be able to extract more information about the rate constants. One such example is the network $G$ in Figure 2(a). For this network we can completely characterize the parameter values for which the associated mass-action system has a complex-balanced realization.

Example 5.2**.**

Consider a mass-action system on the network $G$ of Figure 2(a) and Example 2.7, with rate constants

[TABLE]

By a calculation, (21)-(24) hold if and only if

[TABLE]

Again, a complex-balanced realization is a subgraph of $G^{\prime}$ in Figure 2(b). More precisely, it is the reversible square with one pair of reversible diagonal (either ${\boldsymbol{y}}_{1}\rightleftharpoons{\boldsymbol{y}}_{3}$ or ${\boldsymbol{y}}_{2}\rightleftharpoons{\boldsymbol{y}}_{4}$); which diagonal is needed depends on the magnitudes of $k_{1}k_{3}$ and $k_{2}k_{4}$. The details of this characterization can be found in an upcoming paper [STN].

The complex-balanced realization described (the subgraph of $G^{\prime}$ in Figure 2(b)) has deficiency $\delta_{G^{\prime}}=1$. It is known that if its eight rate constants lie in a toric ideal of codimension $\delta_{G^{\prime}}=1$, then the mass-action system is complex-balanced [TDS]. While these eight rate constants are related to $k_{1},\,k_{2},\,k_{3}$, and $k_{4}$ by several linear equations, we found one explicit condition (25) for when the mass-action system $G_{\boldsymbol{k}}$ of Figure 2(a) is dynamically equivalent to a complex-balanced system.

Finally, note that the network of Example 5.2 gives rise to systems that are equivalent to complex-balanced for certain choices of rate constants, but not for other choices of rate constants. In a follow-up paper we will show that an entire class of networks give rise to systems that are equivalent to complex-balanced for all choice of rate constants. More precisely, we will prove that systems generated by single-target networks that have their (unique) target vertex in the strict relative interior of the convex hull of its source vertices are dynamically equivalent to detailed-balanced mass-action systems for any choice of rate constants [STN].

5.3 Existence of a weakly reversible realization for a mass-action system

While complex-balanced mass-action systems are weakly reversible, not all weakly reversible mass-action systems are complex-balanced. There has been much work on determining when a weakly reversible mass-action system is complex-balanced or not. Nonetheless, weakly reversible mass-action systems always have at least one positive steady state within each stoichiometric compatibility class [Boros_2018] and are conjectured to be persistent, and even permanent [CasianGAC1].

We present a simple nonlinear feasibility problem to determine when a mass-action system is dynamically equivalent to a weakly reversible one. Recall that a mass-action system is weakly reversible if and only if it is complex-balanced for some choice of rate constants. We introduce a scaling factor $\alpha_{{\boldsymbol{y}}_{i}\to{\boldsymbol{y}}_{j}}$ in order to decouple the dynamical equivalence condition from the complex-balanced condition.

Consider a mass-action system $G_{{\boldsymbol{k}}}$, whose vertices are points in $\mathbb{R}^{n}$, and enumerate the set of source vertices in $G$ as $\{{\boldsymbol{y}}_{1},{\boldsymbol{y}}_{2},\dots,{\boldsymbol{y}}_{N}\}$. We set up a nonlinear feasibility problem for unknown rate constants ${\boldsymbol{k}}^{\prime}$ and a scaling factor ${\boldsymbol{\alpha}}$. Search for vectors ${\boldsymbol{k}}^{\prime}=(k_{{\boldsymbol{y}}_{i}\to{\boldsymbol{y}}_{j}})_{i\neq j}$ and ${\boldsymbol{\alpha}}=(\alpha_{{\boldsymbol{y}}_{i}\to{\boldsymbol{y}}_{j}})_{i\neq j}\in\mathbb{R}^{N^{2}-N}$ satisfying.

[TABLE]

If such ${\boldsymbol{k}}^{\prime}$ and ${\boldsymbol{\alpha}}$ exist, then $G_{{\boldsymbol{k}}}$ is dynamically equivalent to a weakly reversible mass-action system. If no solution exists, then $G_{\boldsymbol{k}}$ is not dynamically equivalent to a weakly reversible system.

Equation (26) enforces dynamical equivalence. Equation (27) can be regarded as a complex balancing condition that uses a different set of rate constants $\alpha_{{\boldsymbol{y}}_{i}\to{\boldsymbol{y}}_{j}}k^{\prime}_{{\boldsymbol{y}}_{i}\to{\boldsymbol{y}}_{j}}$. Since $\alpha_{{\boldsymbol{y}}_{i}\to{\boldsymbol{y}}_{j}}k^{\prime}_{{\boldsymbol{y}}_{i}\to{\boldsymbol{y}}_{j}}\neq 0$ if and only if $k^{\prime}_{{\boldsymbol{y}}_{i}\to{\boldsymbol{y}}_{j}}\neq 0$, we preserve the graph structure of $G^{\prime}_{{\boldsymbol{k}}^{\prime}}$. It is well-known that a reaction network is weakly reversible if and only if it is complex-balanced for some choice of rate constants [TDS]. The scaling factor ${\boldsymbol{\alpha}}$ frees the rate constants from the dynamical equivalence constraint.

Note that while (26)-(29) are simple to describe, more sophisticated, computationally efficient methods have been developed [PolylTimeAlgWR, WRAlgDense]. Weak reversibility is a condition of the underlying directed graph. Ultimately one is imposing conditions on the incidence matrix or the Kirchhoff matrix of the network. Algorithms to find weakly reversible realization for a fixed vertex set have been proposed initially using mixed-integer linear programming [WRAlgDense, WRAlgSparse] and later by a polynomial time algorithm based on linear programming [PolylTimeAlgWR]. However, as with previous work on complex-balanced realizations, one must fix the set of vertices to be used in the computation. According to Theorem 4.8, it suffices to find an equivalent network using the existing source vertices. Therefore, the mixed-integer linear programming algorithms proposed in [WRAlgDense, WRAlgSparse] and the polynomial time algorithm in [PolylTimeAlgWR] can be used in conjunction with Theorem 4.8 to completely characterize whether or not a mass-action system $G_{\boldsymbol{k}}$ is dynamically equivalent to a weakly reversible one.

6 Conclusion

If we are looking for a complex-balanced realization of a given polynomial (or power-law) dynamical system, there exists no a priori limit on the number of vertices in the objective network. Moreover, there are no a priori choices for the locations of the vertices. Here we prove that a solution exists if and only if the objective network can be constructed by using only the vertices that are already present in the original system (i.e., the exponents of the monomial terms present in the original system). We also prove that the same is true for detailed-balanced, reversible and weakly reversible systems.

7 Acknowledgment

This work was partially supported by NSF grants DMS-1412643 and DMS-1816238. The work of the third author was partially supported by an NSERC PGS-D award.