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

Gheorghe Craciun Department of Mathematics and Department of Biomolecular
Chemistry, University of Wisconsin-Madison; partially supported by NSF grants
1412643 and 1816238. Jiaxin Jin Department of Mathematics, University of
Wisconsin-Madison Polly Y. Yu Department of Mathematics, University of
Wisconsin-Madison; partially supported by a NSERC PGS-D award and NSF grant
1412643.

###### 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.

##  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

$\displaystyle\frac{d{\boldsymbol{x}}}{dt}\,\,=\\!\\!\sum_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}\in
G}k_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}}{\boldsymbol{x}}^{{\boldsymbol{y}}}({\boldsymbol{y}}^{\prime}-{\boldsymbol{y}}),$
(1)

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

$\displaystyle{\boldsymbol{y}}_{1}\rightarrow{\boldsymbol{z}}_{1}$
$\displaystyle=\begin{pmatrix}1\\\
0\end{pmatrix}\rightarrow\begin{pmatrix}2\\\
0\end{pmatrix},\qquad{\boldsymbol{y}}_{2}\rightarrow{\boldsymbol{z}}_{2}=\begin{pmatrix}1\\\
1\end{pmatrix}\rightarrow\begin{pmatrix}0\\\
2\end{pmatrix},\qquad{\boldsymbol{y}}_{3}\rightarrow{\boldsymbol{z}}_{3}=\begin{pmatrix}0\\\
1\end{pmatrix}\rightarrow\begin{pmatrix}0\\\ 0\end{pmatrix}.$

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.

XY$\bullet$$\bullet$$\bullet$ Figure 1: A reaction network $G$ in
$\mathbb{R}^{2}$ consisting of 3 reactions and 6 vertices. Under mass-action
kinetics, this network gives rise to the classical Lotka-Volterra model for
population dynamics.

######  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

$\displaystyle{\boldsymbol{y}}_{1}=\begin{pmatrix}0\\\
0\end{pmatrix},\quad{\boldsymbol{y}}_{2}=\begin{pmatrix}0\\\
2\end{pmatrix},\quad{\boldsymbol{y}}_{3}=\begin{pmatrix}3\\\
2\end{pmatrix},\quad\text{and}\quad{\boldsymbol{y}}_{4}=\begin{pmatrix}3\\\
0\end{pmatrix}.$

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

$\displaystyle{\boldsymbol{y}}_{5}=\begin{pmatrix}1\\\
1\end{pmatrix}\quad\text{ and }\quad{\boldsymbol{y}}_{6}=\begin{pmatrix}2\\\
1\end{pmatrix}.$

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.

XY$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$${\boldsymbol{y}}_{1}$${\boldsymbol{y}}_{2}$${\boldsymbol{y}}_{3}$${\boldsymbol{y}}_{4}$${\boldsymbol{y}}_{5}$${\boldsymbol{y}}_{6}$

XY$\bullet$$\bullet$$\bullet$$\bullet$

XY$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$

(a)(b)(c)

Figure 2: Examples of reaction networks (a) $G$, (b) $G^{\prime}$ and (c)
$G*$, with labels of vertices shown in (a). The dynamical systems generated by
the network (a) can also be generated by (b) or (c) for well-chosen rate
constants. Note that (b) and (c) are weakly reversible, and (b) is also
reversible.

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
${\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}$

$\displaystyle\frac{d{\boldsymbol{x}}}{dt}\,\,=\\!\\!\sum_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}\in
G}k_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}}{\boldsymbol{x}}^{{\boldsymbol{y}}}({\boldsymbol{y}}^{\prime}-{\boldsymbol{y}}),$
(2)

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

$\displaystyle\frac{d{\boldsymbol{x}}}{dt}$
$\displaystyle=k_{1}x\begin{pmatrix}1\\\
0\end{pmatrix}+k_{2}xy\begin{pmatrix}-1\\\
1\end{pmatrix}+k_{3}y\begin{pmatrix}0\\\
-1\end{pmatrix}=\begin{pmatrix}k_{1}x\hphantom{y}-k_{2}xy\\\ k_{2}xy-
k_{3}y\hphantom{x}\end{pmatrix}.$

This is the Lotka-Volterra population dynamics model.

Given a mass-action system $G_{{\boldsymbol{k}}}$, Equation (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

$\displaystyle\sum_{{\boldsymbol{y}}_{1}\rightarrow{\boldsymbol{y}}_{2}\in
G}k_{{\boldsymbol{y}}_{1}\rightarrow{\boldsymbol{y}}_{2}}{\boldsymbol{x}}^{{\boldsymbol{y}}_{1}}({\boldsymbol{y}}_{2}-{\boldsymbol{y}}_{1})\,\,\,=\\!\\!\sum_{{\boldsymbol{y}}_{1}^{\prime}\rightarrow{\boldsymbol{y}}_{2}^{\prime}\in
G^{\prime}}k^{\prime}_{{\boldsymbol{y}}^{\prime}_{1}\rightarrow{\boldsymbol{y}}^{\prime}_{2}}{\boldsymbol{x}}^{{\boldsymbol{y}}^{\prime}_{1}}({\boldsymbol{y}}^{\prime}_{2}-{\boldsymbol{y}}^{\prime}_{1}),$
(3)

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 Equation (3), a necessary and sufficient condition for dynamical
equivalence is

$\displaystyle\sum_{{\boldsymbol{y}}_{0}\rightarrow{\boldsymbol{y}}\in
G}k_{{\boldsymbol{y}}_{0}\rightarrow{\boldsymbol{y}}}({\boldsymbol{y}}-{\boldsymbol{y}}_{0})\,\,\,=\\!\\!\sum_{{\boldsymbol{y}}_{0}\rightarrow{\boldsymbol{y}}^{\prime}\in
G^{\prime}}k^{\prime}_{{\boldsymbol{y}}_{0}\rightarrow{\boldsymbol{y}}^{\prime}}({\boldsymbol{y}}^{\prime}-{\boldsymbol{y}}_{0})$
(4)

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 Equation (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 same222The _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 2X -¿[$k$] X + Y is dynamically equivalent to the system in
Figure 3(b), consisting of the reactions 2X -¿[$k$] X + Y and 0
¡-[$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}$.

(a)XY$k$$\bullet$

(b)XY$k$$k^{\prime}$$k^{\prime}$$\bullet$$\bullet$

Figure 3: Two dynamically equivalent systems with different stoichiometric
subspaces. Trajectories are confined to the same affine invariant spaces
because their kinetic subspaces are the same.

######  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:

$\displaystyle k_{15}\begin{pmatrix}\,\,\,1\,\,\\\ 1\end{pmatrix}$
$\displaystyle=k^{\prime}_{12}\begin{pmatrix}\,\,\,0\,\,\\\
2\end{pmatrix}+k^{\prime}_{13}\begin{pmatrix}\,\,\,3\,\,\\\
2\end{pmatrix}+k^{\prime}_{14}\begin{pmatrix}\,\,\,3\,\,\\\
0\end{pmatrix},\quad$ $\displaystyle k_{25}\begin{pmatrix}1\\\
-1\end{pmatrix}$ $\displaystyle=k^{\prime}_{21}\begin{pmatrix}0\\\
-2\end{pmatrix}+k^{\prime}_{23}\begin{pmatrix}\,\,\,3\,\,\\\
0\end{pmatrix}+k^{\prime}_{24}\begin{pmatrix}3\\\ -2\end{pmatrix},$
$\displaystyle k_{36}\begin{pmatrix}-1\\\ -1\end{pmatrix}$
$\displaystyle=k^{\prime}_{31}\begin{pmatrix}-3\\\
-2\end{pmatrix}+k^{\prime}_{32}\begin{pmatrix}-3\\\
0\end{pmatrix}+k^{\prime}_{34}\begin{pmatrix}0\\\ -2\end{pmatrix},\quad$
$\displaystyle k_{46}\begin{pmatrix}-1\\\ 1\end{pmatrix}$
$\displaystyle=k^{\prime}_{41}\begin{pmatrix}-3\\\
0\end{pmatrix}+k^{\prime}_{42}\begin{pmatrix}-3\\\
2\end{pmatrix}+k^{\prime}_{43}\begin{pmatrix}\,\,\,0\,\,\\\ 2\end{pmatrix}.$

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

$\displaystyle\frac{d{\boldsymbol{x}}}{dt}\,\,=\\!\\!\sum_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}\in
G}k_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}}{\boldsymbol{x}}^{{\boldsymbol{y}}}({\boldsymbol{y}}^{\prime}-{\boldsymbol{y}}).$

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

$\displaystyle\frac{d{\boldsymbol{x}}}{dt}\,\,=\\!\\!\sum_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}\in
G}k_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}}{\boldsymbol{x}}_{0}^{{\boldsymbol{y}}}({\boldsymbol{y}}^{\prime}-{\boldsymbol{y}})={\boldsymbol{0}}.$
(5)

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

$\displaystyle
k_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}}{\boldsymbol{x}}_{0}^{{\boldsymbol{y}}}=k_{{\boldsymbol{y}}^{\prime}\rightarrow{\boldsymbol{y}}}{\boldsymbol{x}}_{0}^{{\boldsymbol{y}}^{\prime}}.$
(6)

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

$\displaystyle\sum_{{\boldsymbol{y}}_{0}\rightarrow{\boldsymbol{y}}^{\prime}\in
G}k_{{\boldsymbol{y}}_{0}\rightarrow{\boldsymbol{y}}^{\prime}}{\boldsymbol{x}}_{0}^{{\boldsymbol{y}}_{0}}\,\,=\\!\\!\sum_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}_{0}\in
G}k_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}_{0}}{\boldsymbol{x}}_{0}^{{\boldsymbol{y}}}.$
(7)

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

$\displaystyle\frac{d{\boldsymbol{x}}}{dt}\,\,=\\!\\!\sum_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}\in
G}\nu_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}}({\boldsymbol{x}})({\boldsymbol{y}}^{\prime}-{\boldsymbol{y}}),$

where
$\nu_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}}({\boldsymbol{x}})=k_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}}{\boldsymbol{x}}^{{\boldsymbol{y}}}$.
Once the non-linearity 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

$\displaystyle\sum_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}\in
G}J_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}}({\boldsymbol{y}}^{\prime}-{\boldsymbol{y}})={\boldsymbol{0}}.$
(8)

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

$\displaystyle
J_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{\prime}}=J_{{\boldsymbol{y}}^{\prime}\rightarrow{\boldsymbol{y}}}.$
(9)

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

$\displaystyle\sum_{{\boldsymbol{y}}_{0}\rightarrow{\boldsymbol{y}}^{\prime}\in
G}J_{{\boldsymbol{y}}_{0}\rightarrow{\boldsymbol{y}}^{\prime}}=\sum_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}_{0}\in
G}J_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}_{0}}.$ (10)

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

$\displaystyle\begin{array}[]{rrr}k^{\prime}_{\footnotesize\chemfig{0}\to\chemfig{Y}}=3,&\qquad
k^{\prime}_{\footnotesize\chemfig{Y}\to\chemfig{X}+\chemfig{Y}}=\frac{1}{2},&\qquad
k^{\prime}_{\footnotesize\chemfig{X}+\chemfig{Y}\to\chemfig{0}}=1,\displaystyle\vphantom{\frac{}{}}\\\
k^{\prime}_{\footnotesize\chemfig{Y}\to\chemfig{0}}=1,&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}$

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$.

XY$\bullet$$\bullet$$\bullet$$\bullet$312255 Figure 4: An example of a flux
system. The positive numbers on any 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 is complex-balanced.

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:

  1. 1.

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

$\displaystyle\sum_{{\boldsymbol{y}}_{0}\rightarrow{\boldsymbol{y}}\in
G}J_{{\boldsymbol{y}}_{0}\rightarrow{\boldsymbol{y}}}({\boldsymbol{y}}-{\boldsymbol{y}}_{0})\,\,=\\!\\!\sum_{{\boldsymbol{y}}_{0}\rightarrow{\boldsymbol{y}}^{\prime}\in{G^{\prime}}}J^{\prime}_{{\boldsymbol{y}}_{0}\rightarrow{\boldsymbol{y}}^{\prime}}({\boldsymbol{y}}^{\prime}-{\boldsymbol{y}}_{0}).$
(11)

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

$\displaystyle\sum_{{\boldsymbol{y}}_{0}\rightarrow{\boldsymbol{y}}\in
G}J_{{\boldsymbol{y}}_{0}\rightarrow{\boldsymbol{y}}}({\boldsymbol{y}}-{\boldsymbol{y}}_{0})\,\,=\\!\\!\sum_{{\boldsymbol{y}}_{0}\rightarrow{\boldsymbol{y}}\in
G^{\prime}}J^{\prime}_{{\boldsymbol{y}}_{0}\rightarrow{\boldsymbol{y}}}({\boldsymbol{y}}-{\boldsymbol{y}}_{0})\,\,=\\!\\!\sum_{{\boldsymbol{y}}_{0}\rightarrow{\boldsymbol{y}}\in
G^{*}}J^{*}_{{\boldsymbol{y}}_{0}\rightarrow{\boldsymbol{y}}}({\boldsymbol{y}}-{\boldsymbol{y}}_{0}),$

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:

  1. 1.

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

$\displaystyle\hphantom{11}\,\sum_{{\boldsymbol{y}}_{0}\rightarrow{\boldsymbol{y}}\in
G}J_{{\boldsymbol{y}}_{0}\rightarrow{\boldsymbol{y}}}({\boldsymbol{x}})({\boldsymbol{y}}-{\boldsymbol{y}}_{0})\,\,-\\!\\!\sum_{{\boldsymbol{y}}_{0}\rightarrow{\boldsymbol{y}}^{\prime}\in
G^{\prime}}J^{\prime}_{{\boldsymbol{y}}_{0}\rightarrow{\boldsymbol{y}}^{\prime}}({\boldsymbol{x}})({\boldsymbol{y}}^{\prime}-{\boldsymbol{y}}_{0})$
$\displaystyle=\sum_{{\boldsymbol{y}}_{0}\rightarrow{\boldsymbol{y}}\in
G}k_{{\boldsymbol{y}}_{0}\rightarrow{\boldsymbol{y}}}{\boldsymbol{x}}^{{\boldsymbol{y}}_{0}}({\boldsymbol{y}}-{\boldsymbol{y}}_{0})\,\,-\\!\\!\sum_{{\boldsymbol{y}}_{0}\rightarrow{\boldsymbol{y}}^{\prime}\in
G^{\prime}}k^{\prime}_{{\boldsymbol{y}}_{0}\rightarrow{\boldsymbol{y}}^{\prime}}{\boldsymbol{x}}^{{\boldsymbol{y}}_{0}}({\boldsymbol{y}}^{\prime}-{\boldsymbol{y}}_{0})$
$\displaystyle=\frac{{\boldsymbol{x}}^{{\boldsymbol{y}}_{0}}}{{\boldsymbol{x}}_{0}^{{\boldsymbol{y}}_{0}}}\left(\sum_{{\boldsymbol{y}}_{0}\rightarrow{\boldsymbol{y}}\in
G}k_{{\boldsymbol{y}}_{0}\rightarrow{\boldsymbol{y}}}{\boldsymbol{x}}_{0}^{{\boldsymbol{y}}_{0}}({\boldsymbol{y}}-{\boldsymbol{y}}_{0})\,\,-\\!\\!\sum_{{\boldsymbol{y}}_{0}\rightarrow{\boldsymbol{y}}^{\prime}\in
G^{\prime}}k^{\prime}_{{\boldsymbol{y}}_{0}\rightarrow{\boldsymbol{y}}^{\prime}}{\boldsymbol{x}}_{0}^{{\boldsymbol{y}}_{0}}({\boldsymbol{y}}^{\prime}-{\boldsymbol{y}}_{0})\right)$
$\displaystyle=\frac{{\boldsymbol{x}}^{{\boldsymbol{y}}_{0}}}{{\boldsymbol{x}}_{0}^{{\boldsymbol{y}}_{0}}}\left(\sum_{{\boldsymbol{y}}_{0}\rightarrow{\boldsymbol{y}}\in
G}J_{{\boldsymbol{y}}_{0}\rightarrow{\boldsymbol{y}}}({\boldsymbol{x}}_{0})({\boldsymbol{y}}-{\boldsymbol{y}}_{0})\,\,-\\!\\!\sum_{{\boldsymbol{y}}_{0}\rightarrow{\boldsymbol{y}}^{\prime}\in
G^{\prime}}J^{\prime}_{{\boldsymbol{y}}_{0}\rightarrow{\boldsymbol{y}}^{\prime}}({\boldsymbol{x}}_{0})({\boldsymbol{y}}^{\prime}-{\boldsymbol{y}}_{0})\right)$
$\displaystyle={\boldsymbol{0}}.$

∎

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

$\displaystyle
k^{\prime}_{{\boldsymbol{y}}\to{\boldsymbol{y}}^{\prime}}=\frac{J^{\prime}_{{\boldsymbol{y}}\to{\boldsymbol{y}}^{\prime}}}{{\boldsymbol{x}}_{0}^{{\boldsymbol{y}}}}>0,$

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].

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 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

$\displaystyle\sum_{{\boldsymbol{y}}_{0}\to{\boldsymbol{y}}^{\prime}\in
G}J_{{\boldsymbol{y}}_{0}\to{\boldsymbol{y}}^{\prime}}({\boldsymbol{y}}^{\prime}-{\boldsymbol{y}}_{0})={\boldsymbol{0}},$
(12)

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_ appear555I.e.,
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
${\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

$\displaystyle
P_{(G,{\boldsymbol{J}})}({\boldsymbol{y}}^{*})\,\,=\\!\\!\sum_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{*}\in
G}J_{{\boldsymbol{y}}\rightarrow{\boldsymbol{y}}^{*}}\,\,-\sum_{{\boldsymbol{y}}^{*}\rightarrow{\boldsymbol{y}}^{\prime}\in
G}J_{{\boldsymbol{y}}^{*}\rightarrow{\boldsymbol{y}}^{\prime}}.$ (13)

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.

$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$${\boldsymbol{y}}^{*}$${\boldsymbol{y}}_{1}$${\boldsymbol{y}}_{3}$${\boldsymbol{y}}_{2}$${\boldsymbol{y}}_{4}$${\boldsymbol{z}}$

$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$${\boldsymbol{y}}^{*}$${\boldsymbol{y}}_{1}$${\boldsymbol{y}}_{3}$${\boldsymbol{y}}_{2}$${\boldsymbol{y}}_{4}$${\boldsymbol{z}}$

(a) (b)

Figure 5: Illustrating the idea behind Lemma 4.3 in $\mathbb{R}^{3}$. (a)
Assume that ${\boldsymbol{y}}^{*}$ is a virtual source in the flux system
$(G,{\boldsymbol{J}})$. In (b) is an equivalent flux system
$(G^{\prime},{\boldsymbol{J}}^{\prime})$, obtained by redirecting fluxes from
${\boldsymbol{z}}\to{\boldsymbol{y}}^{*}\to{\boldsymbol{y}}_{j}$ as fluxes
from ${\boldsymbol{z}}\to{\boldsymbol{y}}_{j}$.

######  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 Equation (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

$\displaystyle\sum_{j=1}^{M}J^{\prime}_{{\boldsymbol{z}}\rightarrow{\boldsymbol{y}}_{j}}\left({\boldsymbol{y}}_{j}-{\boldsymbol{z}}\right)$
$\displaystyle=\sum_{j=1}^{M}J_{{\boldsymbol{y}}^{*}\rightarrow{\boldsymbol{y}}_{j}}\left({\boldsymbol{w}}_{j}-{\boldsymbol{w}}_{0}\right)=\underbrace{\sum_{j=1}^{M}J_{{\boldsymbol{y}}^{*}\rightarrow{\boldsymbol{y}}_{j}}{\boldsymbol{w}}_{j}}_{=\,\,{\boldsymbol{0}}}\,-\,\,{\boldsymbol{w}}_{0}\sum_{j=1}^{M}J_{{\boldsymbol{y}}^{*}\rightarrow{\boldsymbol{y}}_{j}}=-J_{{\boldsymbol{z}}\rightarrow{\boldsymbol{y}}^{*}}{\boldsymbol{w}}_{0},$

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,

$\displaystyle-P_{(G^{\prime},{\boldsymbol{J}}^{\prime})}({\boldsymbol{z}})$
$\displaystyle=\sum_{j=1}^{M}J^{\prime}_{{\boldsymbol{z}}\rightarrow{\boldsymbol{y}}_{j}}=\sum_{j=1}^{M}J_{{\boldsymbol{y}}^{*}\rightarrow{\boldsymbol{y}}_{j}}=J_{{\boldsymbol{z}}\rightarrow{\boldsymbol{y}}^{*}}=-P_{(G,{\boldsymbol{J}})}({\boldsymbol{z}}).$

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,

$\displaystyle
J^{\prime}_{{\boldsymbol{y}}_{i}\rightarrow{\boldsymbol{y}}_{k}}=J_{{\boldsymbol{y}}_{i}\rightarrow{\boldsymbol{y}}_{k}}+J_{{\boldsymbol{y}}^{*}\rightarrow{\boldsymbol{y}}_{k}}\frac{J_{{\boldsymbol{y}}_{i}\rightarrow{\boldsymbol{y}}^{*}}}{\sum_{{\boldsymbol{y}}_{j}\to{\boldsymbol{y}}^{*}\in
G}J_{{\boldsymbol{y}}_{j}\rightarrow{\boldsymbol{y}}^{*}}}$ (14)

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

$\displaystyle\sum_{j=1}^{M}J_{{\boldsymbol{y}}^{*}\rightarrow{\boldsymbol{y}}_{j}}=\sum_{i=1}^{N}J_{{\boldsymbol{z}}_{i}\rightarrow{\boldsymbol{y}}^{*}}.$

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$,

$\displaystyle
J^{\prime}_{{\boldsymbol{z}}_{1}\rightarrow{\boldsymbol{y}}_{j}}-J_{{\boldsymbol{z}}_{1}\rightarrow{\boldsymbol{y}}_{j}}$
$\displaystyle=\theta
J_{{\boldsymbol{y}}^{*}\rightarrow{\boldsymbol{y}}_{j}},$ $\displaystyle
J^{\prime}_{{\boldsymbol{y}}^{*}\rightarrow{\boldsymbol{y}}_{j}}-J_{{\boldsymbol{y}}^{*}\rightarrow{\boldsymbol{y}}_{j}}$
$\displaystyle=-\theta
J_{{\boldsymbol{y}}^{*}\rightarrow{\boldsymbol{y}}_{j}},$

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

$\displaystyle\hphantom{=\,\,}\left(\text{Final flow from
${\boldsymbol{z}}_{1}$}\right)-\left(\text{Initial flow from
${\boldsymbol{z}}_{1}$}\right)$
$\displaystyle=\sum_{j=1}^{M}J^{\prime}_{{\boldsymbol{z}}_{1}\to{\boldsymbol{y}}_{j}}({\boldsymbol{y}}_{j}-{\boldsymbol{z}}_{1})-J_{{\boldsymbol{z}}_{1}\rightarrow{\boldsymbol{y}}^{*}}({\boldsymbol{y}}^{*}-{\boldsymbol{z}}_{1})$
$\displaystyle=\theta\underbrace{\sum_{j=1}^{M}J_{{\boldsymbol{y}}^{*}\rightarrow{\boldsymbol{y}}_{j}}({\boldsymbol{y}}_{j}-{\boldsymbol{y}}^{*})}_{=\,\,{\boldsymbol{0}}}+\theta\underbrace{\sum_{j=1}^{M}J_{{\boldsymbol{y}}^{*}\rightarrow{\boldsymbol{y}}_{j}}}_{=\,\,\sum_{j=1}^{N}J_{{\boldsymbol{z}}_{i}\rightarrow{\boldsymbol{y}}^{*}}}({\boldsymbol{y}}^{*}-{\boldsymbol{z}}_{1})-J_{{\boldsymbol{z}}_{1}\rightarrow{\boldsymbol{y}}^{*}}({\boldsymbol{y}}^{*}-{\boldsymbol{z}}_{1})$
$\displaystyle={\boldsymbol{0}}.$

At all other vertices, the net flux is unchanged.

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

$\displaystyle
P_{(G^{\prime},{\boldsymbol{J}}^{\prime})}({\boldsymbol{z}}_{1})-P_{G}({\boldsymbol{z}}_{1})$
$\displaystyle=-\sum_{j=1}^{M}J^{\prime}_{{\boldsymbol{z}}_{1}\rightarrow{\boldsymbol{y}}_{j}}+J_{{\boldsymbol{z}}_{1}\rightarrow{\boldsymbol{y}}^{*}}=-\theta\sum_{i=1}^{N}J_{{\boldsymbol{z}}_{i}\rightarrow{\boldsymbol{y}}^{*}}+J_{{\boldsymbol{z}}_{1}\rightarrow{\boldsymbol{y}}^{*}}=0.$

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

$\displaystyle
P_{(G^{\prime},{\boldsymbol{J}}^{\prime})}({\boldsymbol{y}}_{j})-P_{(G,{\boldsymbol{J}})}({\boldsymbol{y}}_{j})$
$\displaystyle=\left(J^{\prime}_{{\boldsymbol{z}}_{1}\rightarrow{\boldsymbol{y}}_{j}}+J^{\prime}_{{\boldsymbol{y}}^{*}\rightarrow{\boldsymbol{y}}_{j}}\right)-\left(J_{{\boldsymbol{z}}_{1}\rightarrow{\boldsymbol{y}}_{j}}+J_{{\boldsymbol{y}}^{*}\rightarrow{\boldsymbol{y}}_{j}}\right)=0.$

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

$\displaystyle
P_{(G^{\prime},{\boldsymbol{J}}^{\prime})}({\boldsymbol{y}}^{*})-P_{(G,{\boldsymbol{J}})}({\boldsymbol{y}}^{*})=-J_{{\boldsymbol{z}}_{1}\to{\boldsymbol{y}}^{*}}+\theta\sum_{j=1}^{M}J_{{\boldsymbol{y}}^{*}\to{\boldsymbol{y}}_{j}}=0$

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

$\displaystyle\sum_{{\boldsymbol{y}}^{*}\rightarrow{\boldsymbol{y}}\in
G^{\prime}}J^{\prime}_{{\boldsymbol{y}}^{*}\rightarrow{\boldsymbol{y}}}({\boldsymbol{y}}-{\boldsymbol{y}}^{*})$
$\displaystyle=(1+\theta)\\!\\!\sum_{{\boldsymbol{y}}^{*}\rightarrow{\boldsymbol{y}}\in
G^{\prime}}J_{{\boldsymbol{y}}^{*}\rightarrow{\boldsymbol{y}}}({\boldsymbol{y}}-{\boldsymbol{y}}^{*})={\boldsymbol{0}},$

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

$\displaystyle{\boldsymbol{0}}=\sum_{{\boldsymbol{y}}^{*}\to{\boldsymbol{y}}\in{\widetilde{G}}}\widetilde{J}_{{\boldsymbol{y}}^{*}\to{\boldsymbol{y}}}({\boldsymbol{y}}-{\boldsymbol{y}}^{*}).$

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.

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

$\displaystyle
k^{\prime}_{{\boldsymbol{z}}_{j}\to{\boldsymbol{y}}_{i}}=k_{{\boldsymbol{z}}_{j}\to{\boldsymbol{y}}_{i}}+\alpha_{j}\frac{\beta_{i}}{\sum\beta_{s}},$

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

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

$\displaystyle\sum_{i=1}^{M}\beta_{i}{\boldsymbol{y}}_{i}=\sum_{i=1}^{M}\beta_{i}{\boldsymbol{y}}^{*}.$

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}$:

$\displaystyle\sum_{i=1}^{M}(k^{\prime}_{{\boldsymbol{z}}_{1}\to{\boldsymbol{y}}_{i}}-k_{{\boldsymbol{z}}_{1}\to{\boldsymbol{y}}_{i}})({\boldsymbol{y}}_{i}-{\boldsymbol{z}}_{1})$
$\displaystyle=\sum_{i=1}^{M}\alpha_{1}\frac{\beta_{i}}{\sum\beta_{s}}({\boldsymbol{y}}_{i}-{\boldsymbol{z}}_{1})$
$\displaystyle=\frac{\alpha_{1}}{\sum\beta_{s}}\left(\sum_{i=1}^{M}\beta_{i}{\boldsymbol{y}}^{*}-\sum_{i=1}^{M}\beta_{i}{\boldsymbol{z}}_{1}\right)$
$\displaystyle=\alpha_{1}({\boldsymbol{y}}^{*}-{\boldsymbol{z}}_{1}),$

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,

$\displaystyle
J^{\prime}_{{\boldsymbol{y}}_{i}\to{\boldsymbol{y}}_{k}}=J_{{\boldsymbol{y}}_{i}\to{\boldsymbol{y}}_{k}}+J_{{\boldsymbol{y}}^{*}\to{\boldsymbol{y}}_{k}}\frac{J_{{\boldsymbol{y}}_{i}\to{\boldsymbol{y}}^{*}}}{\sum_{{\boldsymbol{y}}_{j}\to{\boldsymbol{y}}^{*}\in
G}J_{{\boldsymbol{y}}_{j}\to{\boldsymbol{y}}^{*}}}$ (15)

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

$\displaystyle J^{\prime}_{{\boldsymbol{y}}_{i}\to{\boldsymbol{y}}_{k}}$
$\displaystyle=J_{{\boldsymbol{y}}_{i}\to{\boldsymbol{y}}_{k}}+J_{{\boldsymbol{y}}^{*}\to{\boldsymbol{y}}_{k}}\frac{J_{{\boldsymbol{y}}_{i}\to{\boldsymbol{y}}^{*}}}{\sum_{{\boldsymbol{y}}_{j}\to{\boldsymbol{y}}^{*}\in
G}J_{{\boldsymbol{y}}_{j}\to{\boldsymbol{y}}^{*}}}$
$\displaystyle=J_{{\boldsymbol{y}}_{k}\to{\boldsymbol{y}}_{i}}+J_{{\boldsymbol{y}}^{*}\to{\boldsymbol{y}}_{i}}\frac{J_{{\boldsymbol{y}}_{k}\to{\boldsymbol{y}}^{*}}}{\sum_{{\boldsymbol{y}}_{j}\to{\boldsymbol{y}}^{*}\in
G}J_{{\boldsymbol{y}}_{j}\to{\boldsymbol{y}}^{*}}}$
$\displaystyle=J^{\prime}_{{\boldsymbol{y}}_{k}\to{\boldsymbol{y}}_{i}}.$

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

$\displaystyle
k^{\prime}_{{\boldsymbol{y}}_{i}\to{\boldsymbol{y}}_{j}}=k_{{\boldsymbol{y}}_{i}\to{\boldsymbol{y}}_{j}}+k_{{\boldsymbol{y}}_{i}\to{\boldsymbol{y}}^{*}}\left(\frac{k_{{\boldsymbol{y}}^{*}\to{\boldsymbol{y}}_{j}}}{\sum
k_{{\boldsymbol{y}}^{*}\to{\boldsymbol{y}}_{s}}}\right).$

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. ∎

###  4.1 Connection to deficiency theory

Within the Chemical 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

$\displaystyle\delta_{G}=|V_{G}|-\ell_{G}-s.$ (16)

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
[DEMASControl], where the removal of additional monomials that function as
controls in a feedback system also lead 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 Theorems 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,

$\displaystyle
S_{0}=\operatorname{span}_{\mathbb{R}}\\{{\boldsymbol{y}}\to{\boldsymbol{y}}^{\prime}\in
G\colon{\boldsymbol{y}}\neq{\boldsymbol{y}}^{*}\text{ or
}{\boldsymbol{y}}^{\prime}\neq{\boldsymbol{y}}^{*}\\}.$

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

$\displaystyle
S_{G}=\operatorname{span}_{\mathbb{R}}\\{S_{0},\,{\boldsymbol{y}}_{j}-{\boldsymbol{y}}^{*},\,{\boldsymbol{y}}^{*}-{\boldsymbol{z}}_{i}\\}_{i,j},$

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

$\displaystyle
S_{G^{\prime}}=\operatorname{span}_{\mathbb{R}}\\{S_{0},\,{\boldsymbol{y}}_{j}-{\boldsymbol{z}}_{i}\\}_{i,j}.$

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, so 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 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:

$\displaystyle\sum_{j\neq
i}J^{\prime}_{{\boldsymbol{y}}_{i}\to{\boldsymbol{y}}_{j}}({\boldsymbol{y}}_{j}-{\boldsymbol{y}}_{i})$
$\displaystyle=\sum_{j\neq
i}J_{{\boldsymbol{y}}_{i}\to{\boldsymbol{y}}_{j}}({\boldsymbol{y}}_{j}-{\boldsymbol{y}}_{i}),\qquad\text{
for }i=1,2,\dots,N$ (17) $\displaystyle\sum_{j\neq
i}J^{\prime}_{{\boldsymbol{y}}_{i}\to{\boldsymbol{y}}_{j}}$
$\displaystyle=\sum_{j\neq
i}J^{\prime}_{{\boldsymbol{y}}_{j}\to{\boldsymbol{y}}_{i}},\qquad\text{ for
}i=1,2,\dots,N$ (18a) $\displaystyle{\boldsymbol{J}}^{\prime}$
$\displaystyle\geq 0.$ (19)

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 Equation (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 Equations (17)-(19) hold if and only
if

$\displaystyle J_{1}=J_{3},\quad J_{2}=J_{4},\quad\text{and
}\quad\frac{1}{5}\leq\frac{J_{1}}{J_{2}}\leq 5.$ (20)

A chosen flux ${\boldsymbol{J}}$ that satisfies Equation (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 set up for the detailed-balanced case is defined analogously. We
keep Equation (17) and (19), and include the equation

$\displaystyle
J^{\prime}_{{\boldsymbol{y}}_{i}\to{\boldsymbol{y}}_{j}}=J^{\prime}_{{\boldsymbol{y}}_{j}\to{\boldsymbol{y}}_{i}},\qquad\text{
for }1\leq i\neq j\leq N.$ (18b)

###  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:

$\displaystyle\sum_{j\neq
i}k^{\prime}_{{\boldsymbol{y}}_{i}\to{\boldsymbol{y}}_{j}}({\boldsymbol{y}}_{j}-{\boldsymbol{y}}_{i})$
$\displaystyle=\sum_{j\neq
i}k_{{\boldsymbol{y}}_{i}\to{\boldsymbol{y}}_{j}}({\boldsymbol{y}}_{j}-{\boldsymbol{y}}_{i}),\qquad\text{
for }i=1,2,\dots,N$ (21) $\displaystyle\sum_{j\neq
i}k^{\prime}_{{\boldsymbol{y}}_{i}\to{\boldsymbol{y}}_{j}}{\boldsymbol{x}}^{{\boldsymbol{y}}_{i}}$
$\displaystyle=\sum_{j\neq
i}k^{\prime}_{{\boldsymbol{y}}_{j}\to{\boldsymbol{y}}_{i}}{\boldsymbol{x}}^{{\boldsymbol{y}}_{j}},\qquad\text{
for }i=1,2,\dots,N$ (22) $\displaystyle{\boldsymbol{k}}^{\prime}$
$\displaystyle\geq 0,$ (23) $\displaystyle{\boldsymbol{x}}$ $\displaystyle>0.$
(24)

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

$\displaystyle k_{{\boldsymbol{y}}_{1}\to{\boldsymbol{y}}_{5}}$
$\displaystyle=k_{1},\quad
k_{{\boldsymbol{y}}_{2}\to{\boldsymbol{y}}_{5}}=k_{2},\quad
k_{{\boldsymbol{y}}_{3}\to{\boldsymbol{y}}_{6}}=k_{3},\quad\text{and}\quad
k_{{\boldsymbol{y}}_{4}\to{\boldsymbol{y}}_{6}}=k_{4}.$

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

$\displaystyle\frac{1}{25}\leq\frac{k_{1}k_{3}}{k_{2}k_{4}}\leq 25.$ (25)

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 in
Equation (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:

$\displaystyle\sum_{j\neq
i}k^{\prime}_{{\boldsymbol{y}}_{i}\to{\boldsymbol{y}}_{j}}({\boldsymbol{y}}_{j}-{\boldsymbol{y}}_{i})$
$\displaystyle=\sum_{j\neq
i}k_{{\boldsymbol{y}}_{i}\to{\boldsymbol{y}}_{j}}({\boldsymbol{y}}_{j}-{\boldsymbol{y}}_{i}),\qquad\text{
for }i=1,2,\dots,N$ (26) $\displaystyle\sum_{j\neq
i}\alpha_{{\boldsymbol{y}}_{i}\to{\boldsymbol{y}}_{j}}k^{\prime}_{{\boldsymbol{y}}_{i}\to{\boldsymbol{y}}_{j}}$
$\displaystyle=\sum_{j\neq
i}\alpha_{{\boldsymbol{y}}_{j}\to{\boldsymbol{y}}_{i}}k^{\prime}_{{\boldsymbol{y}}_{j}\to{\boldsymbol{y}}_{i}},\qquad\text{
for }i=1,2,\dots,N$ (27) $\displaystyle{\boldsymbol{k}}^{\prime}$
$\displaystyle\geq 0,$ (28) $\displaystyle{\boldsymbol{\alpha}}$
$\displaystyle>0.$ (29)

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 Equations (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
intially 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 conjuction 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 Acknowledgments

This work was supported in part by the National Science Foundation under
grants DMS–1412643 and DMS–1816238, and National Science and Engineering
Research Council of Canada.

## References