#  An algorithm for finding weakly reversible deficiency zero realizations of
polynomial dynamical systems

Gheorghe Craciun  Department of Mathematics, University of Wisconsin–Madison
Department of Biomolecular Chemistry, University of Wisconsin–Madison  Jiaxin
Jin  Department of Mathematics, The Ohio State University  Polly Y. Yu  NSF-
Simons Center for Mathematical and Statistical Analysis of Biology, Harvard
University

###### Abstract

Systems of differential equations with polynomial right-hand sides are very
common in applications. On the other hand, their mathematical analysis is very
challenging in general, due to the possibility of complex dynamics: multiple
basins of attraction, oscillations, and even chaotic dynamics. Even if we
restrict our attention to mass-action systems, all of these complex dynamical
behaviours are still possible. On the other hand, if a polynomial dynamical
system has a weakly reversible deficiency zero (WR0) realization, then its
dynamics is known to be remarkably simple: oscillations and chaotic dynamics
are ruled out and, up to linear conservation laws, there exists a single
positive steady state, which is asymptotically stable. Here we describe an
algorithm for finding WR0 realizations of polynomial dynamical systems,
whenever such realizations exist.

##  1 Introduction

By a _polynomial dynamical system_ we mean a system of ODEs with polynomial
right-hand side, of the form

$\begin{split}\frac{dx_{1}}{dt}&=p_{1}(x_{1},...,x_{n}),\\\
\frac{dx_{2}}{dt}&=p_{2}(x_{1},...,x_{n}),\\\ &\qquad\quad\vdots\\\
\frac{dx_{n}}{dt}&=p_{n}(x_{1},...,x_{n}),\\\ \end{split}$ (1)

where $p_{i}(x_{1},\ldots,x_{n})\in\mathbb{R}[x_{1},\ldots,x_{n}]$. In
general, such systems are very difficult to analyze due to nonlinearities and
feedback that may give rise to bifurcations, multiple basins of attraction,
oscillations, and even chaotic dynamics. The second part of Hilbert’s 16th
problem (about the number of limit cycles of polynomial dynamical systems in
the plane) is still essentially unsolved, even for quadratic polynomials
[Ilyashenko2002]. Even the simplest object associated to (1), its steady state
set, is central to real algebraic geometry.

In terms of applications, polynomial dynamical systems often show up in, for
example, chemistry, biology, and population dynamics. In these models, the
variable $x_{i}$ typically represents concentration, population, or another
quantity that is strictly positive, so the domain of (1) is restricted to the
positive orthant. For example, in an infectious disease model, an infectious
individual might covert a susceptible individual; this would contribute a
“${}+bxy$” term to $\frac{dx}{dt}$, where $x$ is the population of susceptible
individuals, $y$ the infectious population, and $b>0$ a parameter measuring
the contact rate. Collecting all contributing terms results in an interaction
network. An active area of research is to relate the structure of the
interaction network to the dynamics generated by it [Johnston2014,
PerezmillanDickensteinShiuConradi2012, ThomsonGunawardena2009, JoshiShiu2015,
JoshiShiu2017, BanajiPantea2016, BanajiPantea2018, ShinarFeinberg2010,
MinchevaRoussel2007].

Conversely, one may start with (1) from experimental data, with little or no
information on the generating interaction network. One may try to elucidate
the underlying interaction network; however, without additional assumptions, a
polynomial dynamical system is not uniquely generated by one interaction
network, but infinitely many [CraciunPantea2008]. This lack of identifiability
of the underlying network can actually be leveraged to analyze the dynamics:
if a network with certain properties can be found to generate (1), then we may
be able to immediately infer its dynamical behaviour.

A class of systems whose dynamics is very well understood is the family of
complex-balanced systems [HornJackson1972], which are also called toric
dynamical systems [CraciunDickensteinShiuSturmfels2009]. They can never
exhibit oscillations or chaotic dynamics, and, up to linear conservation laws,
there exists a single positive steady state, which is locally asymptotically
stable [HornJackson1972]. Moreover, this steady state is conjectured to be a
global attractor [Horn1974_GAC].

Not only are the dynamical properties of complex-balanced systems well
understood, but also the network and parameter structures that characterize
them [horn1972necessary]. While in general, there are algebraic conditions on
the parameters necessary for complex-balancing, the exception to this rule is
the case of weakly reversible and deficiency zero (WR0) networks – these
systems are complex-balanced for any choices of parameters, in a sense that
will be made clear below. This fact is very important in applications, because
the exact values of the coefficients in the polynomial right-hand sides of
these dynamical systems are often very difficult to estimate accurately in
practice.

In this paper we describe an efficient algorithm for determining whether a
given polynomial dynamical system admits a WR0 realization, and for finding
such a realization whenever it exists (see Algorithm 1). Our algorithm does
not require solving the differential equation (1), nor does it require solving
for its steady state set. Instead, the algorithm, making use of the geometric
and log-linear structure of WR0 networks, requires as its inputs only the
monomials and the matrix of coefficients. If a WR0 realization exists, in
Theorem 3.12 we provide a bijection between the positive steady state set of
(1) and the solution to a system of linear equations.

The paper is organized as follows. In Section 2 we introduce interaction
networks as embedded in $\mathbb{R}^{n}$ and formalize their relations to
polynomial dynamical systems; we also introduce complex-balanced systems, WR0
networks, and other relevant notions and results. In Section 3.1 we describe
our algorithm for finding a WR0 realization of a given polynomial dynamical
system, whose steady state set is studied in Section 3.2. Our algorithm
applies to the case of where the coefficients in the polynomials are
unspecified; we consider such systems in Section 3.3.

##  2 Background

Throughout this work, we denote by $\mathbb{R}_{\geq}^{n}$ and
$\mathbb{R}_{>}^{n}$ the sets of vectors with non-negative and positive
entries respectively. Similarly, $\mathbb{Z}_{\geq}^{n}$ is the set of vectors
with non-negative integer components. Vectors are typically denoted
${\boldsymbol{x}}$, ${\boldsymbol{y}}$, or ${\boldsymbol{w}}$. We denote by
$\dot{{\boldsymbol{x}}}$ the time-derivative $\frac{d{\boldsymbol{x}}}{dt}$.
For any ${\boldsymbol{x}}\in\mathbb{R}_{>}^{n}$ and
${\boldsymbol{y}}\in\mathbb{R}^{n}$, define the operation
${\boldsymbol{x}}^{{\boldsymbol{y}}}=x_{1}^{y_{1}}x_{2}^{y_{2}}\cdots
x_{n}^{y_{n}}$. If
$\mathbf{Y}=\begin{pmatrix}{\boldsymbol{y}}_{1}&{\boldsymbol{y}}_{2}&\dots&{\boldsymbol{y}}_{n}\end{pmatrix}$,
then
${\boldsymbol{x}}^{\mathbf{Y}}=({\boldsymbol{x}}^{{\boldsymbol{y}}_{1}},{\boldsymbol{x}}^{{\boldsymbol{y}}_{2}},\ldots,{\boldsymbol{x}}^{{\boldsymbol{y}}_{n}})^{\top}$.
The support of a vector ${\boldsymbol{x}}\in\mathbb{R}^{n}$ is the set of
indices $\mathrm{supp}({\boldsymbol{x}})=\\{i\colon x_{i}\neq 0\\}$.

###  2.1 Dynamical systems and Euclidean embedded graphs

In this section, we introduce the Euclidean embedded graph (E-graph), a
directed graph in $\mathbb{R}^{n}$, and explain how a system of differential
equations with polynomial right-hand side (a polynomial dynamical system) is
defined by it.

######  Definition 2.1.

A _Euclidean embedded graph_ (_E-graph_) in $\mathbb{R}^{n}$ is a directed
graph $(V,E)$, where $V$ is a finite subset of $\mathbb{R}_{\geq}^{n}$, and
there are neither self-loops nor isolated vertices. Denote by $V_{s}$ the set
of source vertices.

Let
$V=\\{{\boldsymbol{y}}_{1},{\boldsymbol{y}}_{2},\ldots,{\boldsymbol{y}}_{m}\\}$.
An edge $({\boldsymbol{y}}_{i},{\boldsymbol{y}}_{j})$, or $(i,j)\in E$, is
also denoted ${\boldsymbol{y}}_{i}\to{\boldsymbol{y}}_{j}$. Since vertices are
points in $\mathbb{R}^{n}$, an edge can be regarded as a bona fide vector
between vertices. An _edge vector_ ${\boldsymbol{y}}_{j}-{\boldsymbol{y}}_{i}$
is associated to the edge ${\boldsymbol{y}}_{i}\to{\boldsymbol{y}}_{j}$.

For the purpose of using E-graphs to study polynomial dynamical systems, we
assume $V_{s}\subset\mathbb{Z}_{\geq}^{n}$, even though most results stated in
this paper hold for $V\subset\mathbb{R}_{\geq}^{n}$.

The set of vertices $V$ of $(V,E)$ is partitioned by its connected components,
which we identify by the subset of vertices that belong to that connected
component. If every connected component is strongly connected, i.e., every
edge is part of a cycle, then $(V,E)$ is said to be _weakly reversible_.

Two geometric properties of the E-graph will become important to our analysis
of polynomial dynamical system. The first is a notion of affine independence
within each connected component; the second is a notion of linear independence
between connected components.

######  Definition 2.2.

An E-graph $(V,E)$ has _affinely independent connected components_ if the
vertices in each connected component are affinely independent, i.e., if
$\\{{\boldsymbol{y}}_{0},{\boldsymbol{y}}_{1},\ldots,{\boldsymbol{y}}_{r}\\}\subseteq
V$ is a connected component, then the set
$\\{{\boldsymbol{y}}_{j}-{\boldsymbol{y}}_{0}\colon j=1,2,\ldots,r\\}$ is
linearly independent.

######  Definition 2.3.

Let $(V,E)$ be an E-graph. For any $U\subseteq V$, the _associated linear
subspace of $U$_ is
$S(U)=\operatorname{span}\\{{\boldsymbol{y}}_{j}-{\boldsymbol{y}}_{i}\colon{\boldsymbol{y}}_{i},\,{\boldsymbol{y}}_{j}\in
U\\}$. The _associated linear space_ 111In reaction network theory literature,
the associated linear space is called the _stoichiometric subspace_ of the
network. of $(V,E)$ is

$\displaystyle
S=\operatorname{span}\\{{\boldsymbol{y}}_{j}-{\boldsymbol{y}}_{i}\colon{\boldsymbol{y}}_{i}\to{\boldsymbol{y}}_{j}\in
E\\}.$

If $U$ defines a connected component of $(V,E)$, then $S(U)\subseteq S$.
Indeed, if $V_{1}$, $V_{2},\ldots,V_{\ell}$ are the connected components, then
$S=S(V_{1})+S(V_{2})+\cdots+S(V_{\ell})$.

Thus far, we have defined an E-graph, and introduced several objects and
properties associated to it. We now turn our attention to how such a graph is
canonically associated to dynamics, by assigning a positive weight to each
edge.

######  Definition 2.4.

Let $(V,E)$ be an E-graph. For each
${\boldsymbol{y}}_{i}\to{\boldsymbol{y}}_{j}\in E$, let $\kappa_{ij}>0$ be its
weight, and let ${\boldsymbol{\kappa}}=(\kappa_{ij})\in\mathbb{R}_{>}^{E}$.
The _associated dynamical system_ on $\mathbb{R}_{>}^{n}$ of the weighted
E-graph $(V,E,{\boldsymbol{\kappa}})$ is

$\displaystyle\frac{d{\boldsymbol{x}}}{dt}=\sum_{(i,j)\in
E}\kappa_{ij}{\boldsymbol{x}}^{{\boldsymbol{y}}_{i}}({\boldsymbol{y}}_{j}-{\boldsymbol{y}}_{i}).$
(2)

It is sometimes convenient to refer to $\kappa_{ij}$ even though
${\boldsymbol{y}}_{i}\to{\boldsymbol{y}}_{j}$ may not be an edge in the
network. In such cases, set $\kappa_{ij}=0$.

######  Remark 2.5.

We defined the domain of (2) to be $\mathbb{R}_{>}^{n}$. Systems of ODEs with
polynomial right-hand side do not in general leave $\mathbb{R}_{>}^{n}$
forward-invariant, but if we assume $V\subset\mathbb{Z}_{\geq}^{n}$, the
positive orthant $\mathbb{R}_{>}^{n}$ is indeed forward-invariant under (2)
[Sontag2001].

It is clear that the right-hand side of (2) lies in the associated linear
space $S$, so any solution to (2) is confined to a translate of $S$. By the
above remark, any solution to (2) where $V\subset\mathbb{Z}_{\geq}^{n}$ with
initial condition ${\boldsymbol{x}}_{0}\in\mathbb{R}_{>}^{n}$ is confined to
$({\boldsymbol{x}}_{0}+S)\cap\mathbb{R}_{>}^{n}$, which is called the
_invariant polyhedron of ${\boldsymbol{x}}_{0}$_.

$1$$\bullet$$\bullet$$3$$2$$\bullet$$\bullet$$3$$5$ (a)

$\bullet$$\bullet$$\bullet$$\bullet$$3$$5$$3$$1$$3$$3$$5$$1$$1$ (b)

$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$6$$10$$6$$4$ (c)

Figure 1: Weighted E-graphs from Example 2.6.

######  Example 2.6.

We illustrate the notions and notations defined above. Figure 1 shows three
examples of weighted E-graphs. The graphs in Figures 1(a) and 1(b) are weakly
reversible, but that in Figure 1(c) is not. The graph in Figure 1(a) has two
connected components, each of which is affinely independent; however, those in
Figures 1(b) and 1(c) do _not_ have affinely independent connected components.

The associated dynamical system of Figure 1(a) is

$\displaystyle\frac{d}{dt}\begin{pmatrix}x_{1}\\\ x_{2}\end{pmatrix}$
$\displaystyle=3\begin{pmatrix}2\\\
2\end{pmatrix}+5x_{1}^{2}\begin{pmatrix}[r]-2\\\
2\end{pmatrix}+2x_{1}^{2}x_{2}^{2}\begin{pmatrix}-2\\\
-2\end{pmatrix}+3x_{2}^{2}\begin{pmatrix}[r]2\\\ -2\end{pmatrix}$
$\displaystyle=\begin{pmatrix}6-10x_{1}^{2}-4x_{1}^{2}x_{2}^{2}+6x_{2}^{2}\\\\[5.0pt]
6+10x_{1}^{2}-4x_{1}^{2}x_{2}^{2}-6x_{2}^{2}\end{pmatrix}.$ (3)

The source vertices play the role of exponents in the monomials, thus the set
of source vertices $V_{s}$ determines the monomials in the associated
dynamical system.

It so happens that the weighted E-graphs in Figures 1(b) and 1(c) also have
(2.6) as their associated dynamical systems. We say that the three weighted
E-graphs in Figure 1 are _dynamically equivalent_ , and the weighted graphs
are _realizations_ of the dynamical system (2.6); we define these terms
precisely in Definition 2.9. This example demonstrates that while a weighted
E-graph is associated to a unique dynamical system, the converse is not true;
there is in general infinitely many realizations of a given polynomial
dynamical system [CraciunPantea2008]. This work is concerned with finding a
realization that guarantee certain algebraic and stability properties.

Another way to study the vector field generated by (2) is to use a linear
combination of some fixed vectors, one for each monomial, with the
coefficients given by the strength of the monomials at that point. We give a
name to those fixed vectors.

######  Definition 2.7.

Let $(V,E,{\boldsymbol{\kappa}})$ be a weighted E-graph, and
${\boldsymbol{y}}_{i}\in V_{s}$. The _net direction vector from
${\boldsymbol{y}}_{i}$_ is

$\displaystyle{\boldsymbol{w}}_{i}=\sum_{{\boldsymbol{y}}_{j}\in
V}\kappa_{ij}({\boldsymbol{y}}_{j}-{\boldsymbol{y}}_{i}).$

The _matrix of net direction vectors_ of $(V,E,{\boldsymbol{\kappa}})$ is

$\displaystyle\mathbf{W}=\begin{pmatrix}\vphantom{\displaystyle\sum}{\boldsymbol{w}}_{1}&{\boldsymbol{w}}_{2}&\cdots&{\boldsymbol{w}}_{m}\end{pmatrix}.$

For convenience, we may refer to the net direction vector even if
${\boldsymbol{y}}_{i}\not\in V_{s}$; in this case, let the net direction
vector be zero. Such a net direction vector will _not_ show up as a column of
$\mathbf{W}$.

The matrix $\mathbf{W}$ from Definition 2.7 is also well defined when we start
not with a weighted E-graph, but with a fixed polynomial dynamical system of
the form

$\displaystyle\frac{d{\boldsymbol{x}}}{dt}=\sum_{i=1}^{m}{\boldsymbol{x}}^{{\boldsymbol{y}}_{i}}{\boldsymbol{w}}_{i}.$
(4)

Note that any polynomial dynamical systems can be uniquely written as such,
for some ${\boldsymbol{y}}_{1}$,
${\boldsymbol{y}}_{2},\ldots,{\boldsymbol{y}}_{m}\in\mathbb{Z}_{\geq}^{n}$
distinct, and ${\boldsymbol{w}}_{1}$,
${\boldsymbol{w}}_{2},\ldots,{\boldsymbol{w}}_{m}\in\mathbb{R}^{n}$ non-zero.

######  Definition 2.8.

Consider the polynomial dynamical system (4). The _matrix of source vertices
$\mathbf{Y}_{s}$_ and the _matrix of net direction vectors $\mathbf{W}$_ of
(4) are

$\mathbf{Y}_{s}=\begin{pmatrix}\vphantom{\displaystyle\sum}{\boldsymbol{y}}_{1}&{\boldsymbol{y}}_{2}&\cdots&{\boldsymbol{y}}_{m}\end{pmatrix}\quad\text{and}\quad\mathbf{W}=\begin{pmatrix}\vphantom{\displaystyle\sum}{\boldsymbol{w}}_{1}&{\boldsymbol{w}}_{2}&\cdots&{\boldsymbol{w}}_{m}\end{pmatrix}.$

Clearly, $\dot{{\boldsymbol{x}}}=\mathbf{W}{\boldsymbol{x}}^{\mathbf{Y}_{s}}$.

Thus far, we start with a weighted E-graph $(V,E,{\boldsymbol{\kappa}})$, and
from it, define a dynamical system. The goal of the present work is the
converse direction: start with a polynomial dynamical system, find some
$(V,E,{\boldsymbol{\kappa}})$, ideally with certain properties, that gives
rise to such dynamics. For example, (4) is generated by the graph
${\boldsymbol{y}}_{i}\xrightarrow{\,\,1\,\,}{\boldsymbol{y}}_{i}+{\boldsymbol{w}}_{i}$,
for $i=1,2,\ldots,m$. As Example 2.6 illustrates, there are in general many
weighted E-graphs that can generate the same dynamics.

######  Definition 2.9.

A _realization_ of a polynomial dynamical system
$\dot{{\boldsymbol{x}}}={\boldsymbol{f}}({\boldsymbol{x}})$ is a weighted
E-graph $(V,E,{\boldsymbol{\kappa}})$ whose associated dynamical system is
precisely $\dot{{\boldsymbol{x}}}={\boldsymbol{f}}({\boldsymbol{x}})$. Two
realizations of $\dot{{\boldsymbol{x}}}={\boldsymbol{f}}({\boldsymbol{x}})$
are said to be _dynamically equivalent_.

######  Lemma 2.10 ([CraciunJinYu2019]*Definition 2.6).

The weighted E-graphs $(V,E,{\boldsymbol{\kappa}})$ and
$(V^{\prime},E^{\prime},{\boldsymbol{\kappa}}^{\prime})$ are dynamically
equivalent if and only if the net direction vector from ${\boldsymbol{y}}_{i}$
in $(V,E,{\boldsymbol{\kappa}})$ coincides with that in
$(V^{\prime},E^{\prime},{\boldsymbol{\kappa}}^{\prime})$, for all
${\boldsymbol{y}}_{i}\in V_{s}\cup V^{\prime}_{s}$.

###### Proof.

This follows from the linear independence of monomials as functions on
$\mathbb{R}_{>}^{n}$. ∎

###  2.2 Complex-balanced systems and WR0 systems

General polynomial dynamical systems can display a wide range of dynamical
behaviours, ranging from stable or unstable steady states, limit cycles, and
even chaos. In this work, we are interested in the family of _complex-balanced
systems_ , which enjoy various algebraic and stability properties.

######  Definition 2.11.

Let $(V,E,{\boldsymbol{\kappa}})$ be a weighted E-graph in $\mathbb{R}^{n}$,
and let $\dot{{\boldsymbol{x}}}={\boldsymbol{f}}({\boldsymbol{x}})$ be its
associated dynamical system. A state
${\boldsymbol{x}}^{*}\in\mathbb{R}_{>}^{n}$ is said to be a _positive steady
state_ if ${\boldsymbol{f}}({\boldsymbol{x}}^{*})={\boldsymbol{0}}$. Let
$V_{>}({\boldsymbol{f}})$ be the set of positive steady states. A state
${\boldsymbol{x}}^{*}>{\boldsymbol{0}}$ is a _complex-balanced steady state_
if at every ${\boldsymbol{y}}_{i}\in V$, we have

$\displaystyle\sum_{(i,j)\in
E}\kappa_{ij}({\boldsymbol{x}}^{*})^{{\boldsymbol{y}}_{i}}=\sum_{(j,i)\in
E}\kappa_{ji}({\boldsymbol{x}}^{*})^{{\boldsymbol{y}}_{j}}.$

The equations above can be interpreted as balancing the fluxes flowing across
the vertex ${\boldsymbol{y}}_{i}$. If a weighted E-graph
$(V,E,{\boldsymbol{\kappa}})$ admits one complex-balanced steady state, then
every positive steady state is complex-balanced [HornJackson1972]; such a
$(V,E,{\boldsymbol{\kappa}})$ is called a _complex-balanced system_.

These systems first arose from the study of chemical systems under mass-action
kinetics, as a generalization of thermodynamic equilibrium. The following
theorem lists some of the most important results about complex-balanced
systems. For more details, see [YuCraciun2018_review, GunaNts, Feinberg1987].

######  Theorem 2.12 ([HornJackson1972]).

Let $(V,E,{\boldsymbol{\kappa}})$ be a complex-balanced system, with steady
state ${\boldsymbol{x}}^{*}\in\mathbb{R}_{>}^{n}$, and associated linear space
$S$. Then the following are true:

  1. (i)

All positive steady states are complex-balanced, and there is exactly one
steady state within each invariant polyhedron.

  2. (ii)

Any complex-balanced steady state ${\boldsymbol{x}}$ satisfies
$\ln{\boldsymbol{x}}-\ln{\boldsymbol{x}}^{*}\in S^{\perp}$.

  3. (iii)

The function

$\displaystyle L({\boldsymbol{x}})=\sum_{i=1}^{n}x_{i}(\ln x_{i}-\ln
x^{*}_{i}-1),$

defined on $\mathbb{R}_{>}^{n}$, is a strict Lyapunov function within each
invariant polyhedron $({\boldsymbol{x}}_{0}+S)\cap\mathbb{R}_{>}^{n}$, with a
global minimum at the corresponding complex-balanced steady state.

  4. (iv)

Every complex-balanced steady state is asymptotically stable with respect to
its invariant polyhedron.

Beside these properties, complex-balanced systems enjoy other remarkable
algebraic and dynamical properties. For example, the set of positive steady
states $V_{>}({\boldsymbol{f}})$ admits a monomial parametrization
[CraciunDickensteinShiuSturmfels2009, PerezmillanDickensteinShiuConradi2012].
Each positive steady state ${\boldsymbol{x}}^{*}$ is in fact linearly stable
with respect to its invariant polyhedron [Johnston_note,
BorosMuellerRegensburger2020]. Complex-balanced systems are also conjectured
to be _persistent_ and _permanent_ 222Roughly speaking, persistence is the
property that starting in $\mathbb{R}_{>}^{n}$, the solution is always bounded
away from the boundary of $\mathbb{R}_{>}^{n}$, and permanence occurs when
solutions always converge to a compact subset of the invariant polyhedron.
[CraciunNazarovPantea2013_GAC]. Moreover, the unique steady state is
conjectured to be _globally_ stable within its invariant polyhedron
[Horn1974_GAC]. The _Persistence_ and _Permanence Conjectures_ have been
proved in several cases, such as when there is only one connected component
[Anderson2011_GAC, BorosHofbauer2019_GAC], or the ambient state space is
$\mathbb{R}^{2}$ [CraciunNazarovPantea2013_GAC], or the E-graph is strongly
endotactic [GopalkrishnanMillerShiu2014], or the associated linear space $S$
is of dimension two and all trajectories are bounded [Pantea2012_GAC]. The
_Global Attractor Conjecture_ has also been proved if there is only one
connected component [Anderson2011_GAC, BorosHofbauer2019_GAC], or the E-graph
is strongly endotactic [GopalkrishnanMillerShiu2014], or the ambient state
space is $\mathbb{R}^{3}$ [CraciunNazarovPantea2013_GAC], or when the
associated linear space $S$ is of dimension at most three [Pantea2012_GAC].

Besides dynamical stability, complex-balanced systems are characterized graph-
theoretically and algebraically. Horn proved in [Horn1972] that
$(V,E,{\boldsymbol{\kappa}})$ is complex-balanced if and only if $(V,E)$ is
weakly reversible and ${\boldsymbol{\kappa}}$ satisfies some algebraic
equations, the number of which is measured by a non-negative integer called
the _deficiency_ of $(V,E)$.

######  Definition 2.13.

Let $(V,E)$ be an E-graph with $\ell$ connected components, and let $S$ be its
associated linear space. The _deficiency_ of $(V,E)$ is the integer
$\delta=|V|-\ell-\dim S$.

The notion of deficiency can also be applied to the connected components.
Suppose $V_{1}$, $V_{2},\ldots,V_{\ell}$ are the connected components of
$(V,E)$. The _deficiency of a connected component $V_{p}$_ is
$\delta_{p}=|V_{p}|-1-\dim S(V_{p})$. It is easy to see that

$\displaystyle\delta\geq\sum_{p=1}^{\ell}\delta_{p},$

with equality if and only if $S(V_{1})$, $S(V_{2}),\ldots,S(V_{\ell})$ are
linearly independent. If $\delta=0$, then necessarily $\delta_{p}=0$ for all
$p$.

If $(V,E)$ is weakly reversible and $\delta=0$, then the associated dynamical
system is always complex-balanced, regardless of the choice of
${\boldsymbol{\kappa}}$. This result is known as the Deficiency Zero Theorem
[Feinberg1972, horn1972necessary]. The deficiency is a property of the
E-graph, not of the associated dynamical system, yet in the case of deficiency
zero, it has strong implications on the dynamics. The goal of this paper is to
search for weakly reversible and deficiency zero (WR0) realizations for
polynomial dynamical systems, which are automatically complex-balanced, and
therefore obey the properties listed in Theorem 2.12.

Deficiency also has a _geometric interpretation_ ; $\delta=0$ if and only if
$(V,E)$ has affinely independent connected components $S_{1}$,
$S_{2},\ldots,S_{\ell}$, and the subspaces $S(V_{1})$,
$S(V_{2}),\ldots,S(V_{\ell})$ are linearly independent
[CraciunJohnstonSzederkenyiTonelloTothYu2020, Theorem 9]. Later we make use of
this interpretation when searching for WR0 realizations.

The system (2) admits a matrix decomposition that aids in studying complex-
balanced steady states. For a weighted E-graph $(V,E,{\boldsymbol{\kappa}})$
where $|V|=m$, its associated dynamical system (2) can decomposed be as
$\dot{{\boldsymbol{x}}}=\mathbf{Y}\mathbf{A}_{{\boldsymbol{\kappa}}}{\boldsymbol{x}}^{\mathbf{Y}}$
[HornJackson1972], where

$\displaystyle\mathbf{Y}=\begin{pmatrix}\vphantom{\displaystyle\sum}{\boldsymbol{y}}_{1}&{\boldsymbol{y}}_{2}&\cdots&{\boldsymbol{y}}_{m}\end{pmatrix}$

is a matrix whose columns are the vertices (including both sources and
targets); ${\boldsymbol{x}}^{\mathbf{Y}}$ is the vector of monomials whose
$i$th component is ${\boldsymbol{x}}^{{\boldsymbol{y}}_{i}}$, and the
_Kirchoff matrix_

$\displaystyle[\mathbf{A}_{{\boldsymbol{\kappa}}}]_{ij}=\left\\{\begin{array}[]{cl}\kappa_{ji}&\text{if
}{\boldsymbol{y}}_{i}\to{\boldsymbol{y}}_{j}\in E\\\\[5.0pt]
\displaystyle-\sum_{r}\kappa_{jr}&\text{for }i=j\\\\[5.0pt]
0&\text{otherwise}\end{array}\right.$

is the negative transpose of the graph Laplacian of
$(V,E,{\boldsymbol{\kappa}})$. In general, the $i$th component of
$\mathbf{A}_{{\boldsymbol{\kappa}}}{\boldsymbol{x}}^{\mathbf{Y}}$

$\displaystyle[\mathbf{A}_{{\boldsymbol{\kappa}}}{\boldsymbol{x}}^{\mathbf{Y}}]_{i}=\sum_{(j,i)\in
E}\kappa_{ji}{\boldsymbol{x}}^{{\boldsymbol{y}}_{j}}-{\boldsymbol{x}}^{{\boldsymbol{y}}_{i}}\\!\\!\sum_{(i,j)\in
E}\kappa_{ij}$

measures the net flux passing through the $i$th vertex, so a complex-balanced
steady state ${\boldsymbol{x}}^{*}$ is a solution to the equation
$\mathbf{A}_{{\boldsymbol{\kappa}}}({\boldsymbol{x}}^{*})^{\mathbf{Y}}={\boldsymbol{0}}$.

A subgraph $(V_{0},E_{0})\subseteq(V,E)$ is a _terminal strongly connected
component_ if it is strongly connected, and there does not exist an edge in
$E$ from a vertex in $V_{0}$ to a vertex in $V\setminus V_{0}$. The kernel of
$\mathbf{A}_{{\boldsymbol{\kappa}}}$ is supported on the terminal strongly
connected components:

######  Theorem 2.14 ([FeinbergHorn1977]).

Let $\mathbf{A}_{{\boldsymbol{\kappa}}}$ be the Kirchoff matrix of
$(V,E,{\boldsymbol{\kappa}})$ with terminal strongly connected components
$V_{1}$, $V_{2},\ldots,V_{t}$. There exists a basis
$\\{{\boldsymbol{c}}_{1},{\boldsymbol{c}}_{2},\ldots,{\boldsymbol{c}}_{t}\\}$
for $\ker\mathbf{A}_{{\boldsymbol{\kappa}}}$ with
${\boldsymbol{c}}_{p}\in\mathbb{R}_{\geq}^{n}$ and

$\displaystyle\left\\{\begin{array}[]{cl}[{\boldsymbol{c}}_{p}]_{i}>0&\text{
if }{\boldsymbol{y}}_{i}\in V_{p},\\\\[5.0pt]
[{\boldsymbol{c}}_{p}]_{i}=0&\text{ otherwise.}\end{array}\right.$

According to the Matrix-Tree Theorem [CraciunDickensteinShiuSturmfels2009,
horn1972necessary], there is an explicit formula for the entries of
${\boldsymbol{c}}_{p}$. Each non-zero $[{\boldsymbol{c}}_{p}]_{i}$ is a
polynomial of $\kappa_{ij}$ with positive coefficients, given by the maximal
minors of $\mathbf{A}_{{\boldsymbol{\kappa}}}$
[CraciunDickensteinShiuSturmfels2009, StanleyFomin1999_MtxTree,
MirzaevGunawardena2013].

If $(V,E)$ is weakly reversible, then
$\delta=\dim(\ker\mathbf{Y}\cap\operatorname{im}\mathbf{A}_{{\boldsymbol{\kappa}}})$.
More generally,
$\dim(\ker\mathbf{Y}\cap\operatorname{im}\mathbf{A}_{{\boldsymbol{\kappa}}})=|V|-\ell-t$,
when $(V,E,{\boldsymbol{\kappa}})$ has $t$ terminal strongly connected
components [FeinbergHorn1977]. Therefore if $(V,E)$ is WR0, then
$\ker(\mathbf{Y}\mathbf{A}_{{\boldsymbol{\kappa}}})=\ker\mathbf{A}_{{\boldsymbol{\kappa}}}$,
and the matrix of net direction vectors $\mathbf{W}$ coincides with
$\mathbf{Y}\mathbf{A}_{{\boldsymbol{\kappa}}}$ (see Lemma 3.1).

For the purpose of this work, we assume that we are given $\mathbf{W}$ and the
matrix of source vertices $\mathbf{Y}_{s}$, but we do not know the
decomposition of $\mathbf{W}$ into the product
$\mathbf{Y}\mathbf{A}_{{\boldsymbol{\kappa}}}$, where the columns of
$\mathbf{Y}_{s}$ are also columns of $\mathbf{Y}$. Because
$\ker\mathbf{A}_{{\boldsymbol{\kappa}}}$ is well characterized
[FeinbergHorn1977, GunaNts, feinberg2019foundations], we make use of it in our
search for WR0 realizations.

$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\kappa_{12}$$\kappa_{21}$$\kappa_{34}$$\kappa_{43}$$\kappa_{56}$$\kappa_{67}$$\kappa_{75}$$\kappa_{87}$$\kappa_{84}$${\boldsymbol{y}}_{1}$${\boldsymbol{y}}_{2}$${\boldsymbol{y}}_{3}$${\boldsymbol{y}}_{4}$${\boldsymbol{y}}_{5}$${\boldsymbol{y}}_{6}$${\boldsymbol{y}}_{7}$${\boldsymbol{y}}_{8}$
Figure 2: A weighted E-graph with two connected components but three terminal
strongly connected components (boxed). Its Kirchoff matrix
$\mathbf{A}_{{\boldsymbol{\kappa}}}$ and a basis for
$\ker\mathbf{A}_{{\boldsymbol{\kappa}}}$ are given in Example 2.15.

######  Example 2.15.

Consider the weighted E-graph $(V,E,{\boldsymbol{\kappa}})$ in Figure 2. While
$(V,E)$ has two connected components, it has three terminal strongly connected
components (boxed in Figure 2). With the ordering of vertices as labelled in
the figure, the Kirchoff matrix of $(V,E,{\boldsymbol{\kappa}})$ is

$\displaystyle\mathbf{A}_{{\boldsymbol{\kappa}}}=\begin{pmatrix}-\kappa_{12}&\hphantom{-}\kappa_{21}\\\
\hphantom{-}\kappa_{12}&-\kappa_{21}\\\
&&-\kappa_{34}&\hphantom{-}\kappa_{43}\\\
&&\hphantom{-}\kappa_{34}&-\kappa_{43}&&&&\kappa_{84}\\\
&&&&-\kappa_{56}&\hphantom{-}0&\hphantom{-}\kappa_{75}\\\
&&&&\hphantom{-}\kappa_{56}&-\kappa_{67}&\hphantom{-}0\\\
&&&&\hphantom{-}0&\hphantom{-}\kappa_{67}&-\kappa_{75}&\kappa_{87}\\\
&&&&&&&-\kappa_{84}-\kappa_{87}\end{pmatrix}.$

A basis for its kernel is given by the vectors

$\displaystyle{\boldsymbol{c}}_{1}=\begin{pmatrix}\kappa_{21}\\\
\kappa_{12}\\\ 0\\\ 0\\\ 0\\\ 0\\\ 0\\\
0\end{pmatrix},\qquad{\boldsymbol{c}}_{2}=\begin{pmatrix}0\\\ 0\\\
\kappa_{43}\\\ \kappa_{34}\\\ 0\\\ 0\\\ 0\\\
0\end{pmatrix},\qquad{\boldsymbol{c}}_{3}=\begin{pmatrix}0\\\ 0\\\ 0\\\ 0\\\
\kappa_{67}\kappa_{75}\\\ \kappa_{56}\kappa_{75}\\\ \kappa_{56}\kappa_{67}\\\
0\end{pmatrix}.$

The supports of the basis vectors ${\boldsymbol{c}}_{p}$ are precisely the
terminal strongly connected components of $(V,E)$. If the graph is weakly
reversible, then the basis of $\ker\mathbf{A}_{{\boldsymbol{\kappa}}}$ given
in Theorem 2.14 provides a way to partition the set of vertices.

##  3 Main results

In this section, we present Algorithm 1 (see also Figure 3) that searches for
a weakly reversible and deficiency zero (WR0) realization of a given system of
polynomial differential equations

$\displaystyle\frac{d{\boldsymbol{x}}}{dt}=\sum_{i=1}^{m}{\boldsymbol{x}}^{{\boldsymbol{y}}_{i}}{\boldsymbol{w}}_{i},$
(5)

where
${\boldsymbol{y}}_{1},\ldots,{\boldsymbol{y}}_{m}\in\mathbb{Z}_{\geq}^{n}$ are
distinct, and ${\boldsymbol{w}}_{1}$,
${\boldsymbol{w}}_{2},\ldots,{\boldsymbol{w}}_{m}\in\mathbb{R}^{n}\setminus\\{{\boldsymbol{0}}\\}$.
Whenever (5) admits a WR0 realization, the system is complex-balanced and
enjoys all the properties listed in Theorem 2.12. Moreover, if no WR0
realization exists for (5), our algorithm would conclude as much. Whenever a
WR0 realization exists, the set of positive steady states has a log-linear
structure that allows us to easily find the steady states of (5), as outlined
in Theorem 3.12. Finally, our algorithm is valid even if the
${\boldsymbol{w}}_{i}$’s are only known up to a positive scalar multiple; we
prove this in Theorem 3.13.

###  3.1 Algorithm for WR0 realization

$\mathbf{Y}_{s}=\begin{pmatrix}{\boldsymbol{y}}_{1},\ldots,{\boldsymbol{y}}_{m}\end{pmatrix}\in\mathbb{Z}_{\geq}^{n\times
m}$$\mathbf{W}=\begin{pmatrix}{\boldsymbol{w}}_{1},\ldots,{\boldsymbol{w}}_{m}\end{pmatrix}\in\mathbb{R}^{n\times
m}$Find minimal set of
generators$\\{{\boldsymbol{c}}_{1},\ldots,{\boldsymbol{c}}_{\ell}\\}\subset\mathbb{R}_{\geq}^{m}$for
the cone $\ker\mathbf{W}\cap\mathbb{R}_{\geq}^{m}$Does
$\\{\mathrm{supp}({\boldsymbol{c}}_{p})\\}_{p=1}^{\ell}$partition
$[1,m]\subseteq\mathbb{N}$?No WR0 realizationnoDefine candidate
connectedcomponents by
$V_{p}\coloneqq\mathrm{supp}({\boldsymbol{c}}_{p})$yesIs
$\\{{\boldsymbol{y}}_{i}\colon i\in V_{p}\\}$affinely independent?for p =
1,2,...,l no$\forall\,i\in V_{p}$, is
${\boldsymbol{w}}_{i}\in\mathrm{Cone}\\{{\boldsymbol{y}}_{j}-{\boldsymbol{y}}_{i}\colon
j\in V_{p}\\}$?yesUniquely decompose each ${\boldsymbol{w}}_{i}$as
$\displaystyle{\boldsymbol{w}}_{i}=\sum_{\begin{subarray}{c}j\neq i\\\ j\in
V_{p}\end{subarray}}\kappa_{ij}({\boldsymbol{y}}_{j}-{\boldsymbol{y}}_{i})$yesnoWR0
realization foundendfor Figure 3: Algorithm 1 for finding WR0 realization of a
polynomial dynamical system
$\dot{{\boldsymbol{x}}}=\sum_{i=1}^{m}{\boldsymbol{x}}^{{\boldsymbol{y}}_{i}}{\boldsymbol{w}}_{i}$.

The inputs of Algorithm 1 are the source vertices and their net direction
vectors via $\mathbf{Y}_{s}$ and $\mathbf{W}$ respectively. To find a WR0
realization $(V,E,{\boldsymbol{\kappa}})$ is to find a matrix decomposition of
$\mathbf{W}=\mathbf{Y}_{s}\mathbf{A}_{{\boldsymbol{\kappa}}}$, where
$\mathbf{A}_{{\boldsymbol{\kappa}}}$ encodes the graph structure of $(V,E)$.
In the following lemma, we prove properties that can be expected should a WR0
realization exists.

Recall that a set $X$ is a _polyhedral cone_ if
$X=\\{{\boldsymbol{x}}\colon\mathbf{M}{\boldsymbol{x}}\leq{\boldsymbol{0}}\\}$
for some matrix $\mathbf{M}$. Such a cone is convex. It is _pointed_ , or
_strongly convex_ , if it does not contain a positive dimensional linear
subspace. Note that a cone contained in the positive orthant
$\mathbb{R}_{\geq}^{m}$ is always pointed. A pointed polyhedral cone admits a
unique (up to scalar multiple) minimal set of generators
[CoxLittleSchenck2011].

Algorithm 1 (WR0 algorithm)

0: Matrices
$\mathbf{Y}_{s}=\begin{pmatrix}{\boldsymbol{y}}_{1},\ldots,{\boldsymbol{y}}_{m}\end{pmatrix}$
and
$\mathbf{W}=\begin{pmatrix}{\boldsymbol{w}}_{1},\ldots,{\boldsymbol{w}}_{m}\end{pmatrix}$
that define
$\dot{{\boldsymbol{x}}}=\sum_{i=1}^{m}{\boldsymbol{x}}^{{\boldsymbol{y}}_{i}}{\boldsymbol{w}}_{i}$.

0: Either return the unique WR0 realization, or print that such a realization
does not exist.

1: Find the minimal set of generators
$\\{{\boldsymbol{c}}_{1},{\boldsymbol{c}}_{2},\ldots,{\boldsymbol{c}}_{\ell}\\}$
of the pointed convex cone $\ker\mathbf{W}\cap\mathbb{R}_{\geq}^{m}$.

2: if the sets
$\mathrm{supp}({\boldsymbol{c}}_{1}),\mathrm{supp}({\boldsymbol{c}}_{2}),\dots,\mathrm{supp}({\boldsymbol{c}}_{\ell})$
do not form a partition of $\\{1,2,\dots,m\\}$ then

3: Print: WR0 realization does not exist. Exit.

4: else

5: Define the candidate connected components to be
$V_{p}\coloneqq\mathrm{supp}({\boldsymbol{c}}_{p})$, for $p=1,2,\dots,\ell.$

6: for $p=1,2,\dots,\ell$ do

7: if the vectors $\\{{\boldsymbol{y}}_{i}\colon i\in V_{p}\\}$ are not
affinely independent then

8: Print: WR0 realization does not exist. Exit.

9: else

10: for each $i\in V_{p}$ do

11: if
${\boldsymbol{w}}_{i}\notin\mathrm{Cone}\\{{\boldsymbol{y}}_{j}-{\boldsymbol{y}}_{i}\colon
j\in V_{p}\\}$ then

12: Print: WR0 realization does not exist. Exit.

13: else

14: Uniquely decompose
${\boldsymbol{w}}_{i}=\sum_{j}\kappa_{ij}({\boldsymbol{y}}_{j}-{\boldsymbol{y}}_{i})$
with $\kappa_{ij}\geq 0$.

15: Add $\\{{\boldsymbol{y}}_{i}\to{\boldsymbol{y}}_{j}\colon\kappa_{ij}>0\\}$
to edge set $E$.

16: end if

17: end for

18: end if

19: end for

20: end if

21: Print: The WR0 realization does exist, and has $\ell$ connected
components.

22: Print: The connected components of this realization are given by
$\\{V_{p}\\}_{p=1}^{\ell}$.

23: Print: The edges are given by $E$, with weights $\kappa_{ij}$.

######  Lemma 3.1.

Suppose a polynomial dynamical system
$\dot{{\boldsymbol{x}}}={\boldsymbol{f}}({\boldsymbol{x}})$ admits a WR0
realization $(V,E,{\boldsymbol{\kappa}})$ with $\ell$ connected components.
Let $\mathbf{Y}_{s}$ be the matrix of source vertices, and $\mathbf{W}$ the
matrix of net direction vectors of the polynomial dynamical system
$\dot{{\boldsymbol{x}}}={\boldsymbol{f}}({\boldsymbol{x}})$. Let $S$ be the
associated linear space, and $\mathbf{A}_{{\boldsymbol{\kappa}}}$ the Kirchoff
matrix of the weighted E-graph $(V,E,{\boldsymbol{\kappa}})$. Then we have:

  1. (i)

$\mathbf{W}=\mathbf{Y}_{s}\mathbf{A}_{{\boldsymbol{\kappa}}}$ and
$\ker\mathbf{W}=\ker\mathbf{A}_{{\boldsymbol{\kappa}}}$,

  2. (ii)

$S=\operatorname{im}\mathbf{W}$, and the rank of $\mathbf{W}$ is $|V|-\ell$,

  3. (iii)

$\ker\mathbf{W}\cap\mathbb{R}_{\geq}^{m}$ is a pointed polyhedral cone, and

  4. (iv)

a minimal set of generators for $\ker\mathbf{W}\cap\mathbb{R}_{\geq}^{m}$ has
$\ell$ elements, whose supports correspond to the connected components of
$(V,E)$.

###### Proof.

Because $(V,E)$ is weakly reversible, all vertices in $V$ are sources.
Moreover, because the deficiency of $(V,E)$ is zero, by [CraciunJinYu_unique,
Proposition 3.5], the net direction vector from any ${\boldsymbol{y}}_{i}$ is
non-zero, and the set of source vertices corresponds exactly to the columns of
$\mathbf{Y}_{s}$.

  1. (i)

By definition,
${\boldsymbol{f}}({\boldsymbol{x}})=\mathbf{W}{\boldsymbol{x}}^{\mathbf{Y}_{s}}$,
and by dynamical equivalence,
${\boldsymbol{f}}({\boldsymbol{x}})=\mathbf{Y}_{s}\mathbf{A}_{{\boldsymbol{\kappa}}}{\boldsymbol{x}}^{\mathbf{Y}_{s}}$.
Since the coefficients of polynomial functions are uniquely determined,
$\mathbf{W}=\mathbf{Y}_{s}\mathbf{A}_{{\boldsymbol{\kappa}}}$. Because
$\dim(\ker\mathbf{Y}_{s}\cap\operatorname{im}\mathbf{A}_{{\boldsymbol{\kappa}}})=\delta=0$,
we have $\ker\mathbf{W}=\ker\mathbf{A}_{{\boldsymbol{\kappa}}}$.

  2. (ii)

Note that

$\displaystyle\operatorname{rank}\mathbf{W}=\operatorname{rank}\mathbf{A}_{{\boldsymbol{\kappa}}}=|V|-\ell=\dim
S,$

where the first and last equalities follow from
$0=\delta=\dim(\ker\mathbf{Y}_{s}\cap\operatorname{im}\mathbf{A}_{{\boldsymbol{\kappa}}})=|V|-\ell-\dim
S$, and the second equality follows from weak reversibility and Theorem 2.14.
Clearly $\operatorname{im}\mathbf{W}\subseteq S$, so
$\operatorname{im}\mathbf{W}=S$.

  3. (iii)

The set $\ker\mathbf{W}\cap\mathbb{R}_{\geq}^{m}$ is the solution to
$\mathbf{W}{\boldsymbol{\nu}}\geq{\boldsymbol{0}}$,
$-\mathbf{W}{\boldsymbol{\nu}}\geq{\boldsymbol{0}}$, and
$\mathbf{Id}\,{\boldsymbol{\nu}}\geq{\boldsymbol{0}}$; thus the set is a
polyhedral cone. That $\ker\mathbf{W}\cap\mathbb{R}_{\geq}^{m}$ is pointed
follows from it being a subset of $\mathbb{R}_{\geq}^{m}$.

  4. (iv)

Let $\mathcal{B}=\\{{\boldsymbol{c}}_{1}$,
${\boldsymbol{c}}_{2},\ldots,{\boldsymbol{c}}_{\ell}\\}$ be a basis of
$\ker\mathbf{A}_{{\boldsymbol{\kappa}}}$ as in Theorem 2.14, where
${\boldsymbol{c}}_{p}\geq{\boldsymbol{0}}$, and each
$V_{p}=\\{{\boldsymbol{y}}_{i}\colon
i\in\mathrm{supp}({\boldsymbol{c}}_{p})\\}$ is a connected component of
$(V,E)$. Clearly
$\mathcal{B}\subseteq\ker\mathbf{W}\cap\mathbb{R}_{\geq}^{m}$; we claim that
$\mathcal{B}$ is a _minimal_ set of _generators_ for the pointed cone.

Let ${\boldsymbol{\nu}}\in\ker\mathbf{W}\cap\mathbb{R}_{\geq}^{m}$ be
arbitrary. By (ii), $\mathcal{B}$ is a basis for $\ker\mathbf{W}$, so
decompose ${\boldsymbol{\nu}}$ accordingly:

$\displaystyle{\boldsymbol{\nu}}=\sum_{p=1}^{\ell}\lambda_{p}{\boldsymbol{c}}_{p},$

for some $\lambda_{p}\in\mathbb{R}$. By Theorem 2.14, each
${\boldsymbol{c}}_{p}$ is supported on the connected components of $(V,E)$,
which partition the set of vertices. In particular, for each $i=1,\ldots,m$,
there is exactly one $p(i)$ such that
$\nu_{i}=\lambda_{p(i)}[{\boldsymbol{c}}_{p(i)}]_{i}$. Since
${\boldsymbol{\nu}}$, ${\boldsymbol{c}}_{p(i)}\in\mathbb{R}_{\geq}^{m}$, it
must be the case that $\lambda_{p(i)}\geq 0$. In other words, $\mathcal{B}$
generates the cone $\ker\mathbf{W}\cap\mathbb{R}_{\geq}^{m}$. Because the
vectors in $\mathcal{B}$ have disjoint supports, $\mathcal{B}$ is minimal. ∎

######  Lemma 3.2.

If Algorithm 1 exits at lines 3, 8, or 12, then
$\dot{{\boldsymbol{x}}}=\sum_{i=1}^{m}{\boldsymbol{x}}^{{\boldsymbol{y}}_{i}}{\boldsymbol{w}}_{i}$
does not admit a WR0 realization.

###### Proof.

If the algorithm exits at line 3, then by the contrapositive of Lemma 3.1(iv)
no WR0 realization exists. Continuing with the algorithm, let
$\\{{\boldsymbol{c}}_{1},\ldots,{\boldsymbol{c}}_{\ell}\\}$ be a minimal set
of generators of $\ker\mathbf{W}\cap\mathbb{R}_{\geq}^{n}$, and partition the
vertices as $V_{p}\coloneqq\mathrm{supp}({\boldsymbol{c}}_{p})$. If instead
the algorithm exits at line 8, then again no WR0 realization exists because
WR0 realizations have affinely independent connected components
[CraciunJohnstonSzederkenyiTonelloTothYu2020, Theorem 9]. Finally, exiting at
line 12 means that some net direction vector ${\boldsymbol{w}}_{i}$ cannot be
decomposed as edges from ${\boldsymbol{y}}_{i}$ to other vertices in $V_{p}$,
which defines a connected component of a WR0 realization if it exists
according to Lemma 3.1(iv). ∎

######  Lemma 3.3.

Suppose Algorithm 1 reaches line 23. Then the connected components of $(V,E)$
are given by $V_{1}$, $V_{2},\ldots,V_{\ell}$.

###### Proof.

If the algorithm reaches line 23, a realization has been found with edges
among $V_{1},\ldots V_{\ell}$, i.e., the connected components are subsets of
$V_{p}$. We prove now that in fact, each $V_{p}$ is connected in $(V,E)$.

For any $p=1,\ldots,\ell$, let $U\coloneqq V_{p}$ be the support of
${\boldsymbol{c}}\coloneqq{\boldsymbol{c}}_{p}$. Suppose for a contradiction
that $V^{*}\subsetneq U$ is a connected component. Because the if statement in
line 7 is false, $U$ is affinely independent, so the linear subspaces

$\displaystyle
S(V^{*})=\\{{\boldsymbol{y}}_{i}-{\boldsymbol{y}}_{j}\colon{\boldsymbol{y}}_{i},{\boldsymbol{y}}_{j}\in
V^{*}\\}\quad\text{and}\quad S(U\setminus
V^{*})=\\{{\boldsymbol{y}}_{i}-{\boldsymbol{y}}_{j}\colon{\boldsymbol{y}}_{i},{\boldsymbol{y}}_{j}\in
U\setminus V^{*}\\}$

are linearly independent. Because
${\boldsymbol{c}}\in\ker\mathbf{W}\cap\mathbb{R}_{\geq}^{m}$ and
$\mathrm{supp}({\boldsymbol{c}})=U$, we have

$\displaystyle{\boldsymbol{0}}$ $\displaystyle=\sum_{i\in
V^{*}}c_{i}{\boldsymbol{w}}_{i}+\\!\\!\sum_{i\in U\setminus
V^{*}}\\!\\!c_{i}{\boldsymbol{w}}_{i},$ (6)

with $c_{i}>0$. Furthermore, the if statement in line 11 returning false
implies that each ${\boldsymbol{w}}_{i}$ in (6) can be further decomposed as
edges between vertices in $U$. For any ${\boldsymbol{y}}_{i}$ in $V^{*}$,
which is a connected component, the net direction vector
${\boldsymbol{w}}_{i}$ is a positive linear combination of edge vectors
between ${\boldsymbol{y}}_{i}$ and other vertices in $V^{*}$, so
${\boldsymbol{w}}_{i}\in S(V^{*})$. In particular,
${\boldsymbol{c}}^{*}\coloneqq\sum_{i\in V^{*}}c_{i}{\boldsymbol{w}}_{i}\in
S(V^{*})$. Similarly, any vertices in $U\setminus V^{*}$ are only connected to
other vertices in $U\setminus V^{*}$. Linear independence of $S(V^{*})$ and
$S(U\setminus V^{*})$ means that the vectors ${\boldsymbol{c}}^{*}$ and
${\boldsymbol{c}}-{\boldsymbol{c}}^{*}$ are linearly independent. Both
${\boldsymbol{c}}^{*}$ and ${\boldsymbol{c}}-{\boldsymbol{c}}^{*}$ lie in the
cone $\ker\mathbf{W}\cap\mathbb{R}_{\geq}^{m}$, so
$\\{{\boldsymbol{c}}_{1},\ldots,{\boldsymbol{c}}_{\ell}\\}$ is not a set of
generators, which is a contradiction. ∎

######  Lemma 3.4.

Suppose Algorithm 1 reaches line 23. Then the deficiency of $(V,E)$ is zero.

###### Proof.

The falsity of the if statement in line 7 and Lemma 3.3 imply that the
connected components $V_{1}$, $V_{2},\ldots,V_{\ell}$ are affinely
independent. To prove $\delta=0$, it remains to show that $S(V_{1})$,
$S(V_{2}),\ldots,S(V_{\ell})$ are linearly independent subspaces
[CraciunJohnstonSzederkenyiTonelloTothYu2020, Theorem 9].

We claim that the minimal set of generators
$\\{{\boldsymbol{c}}_{1},{\boldsymbol{c}}_{2},\ldots,{\boldsymbol{c}}_{\ell}\\}$
forms a basis for $\ker\mathbf{W}$. Let ${\boldsymbol{c}}\in\ker\mathbf{W}$ be
arbitrary. If ${\boldsymbol{c}}$ has non-negative components, then it is a
linear combination of ${\boldsymbol{c}}_{p}$’s. If
${\boldsymbol{c}}\not\in\mathbb{R}_{\geq}^{m}$, then there exist sufficiently
large constants $\mu_{p}>0$ so that

$\displaystyle{\boldsymbol{c}}+\sum_{p=1}^{\ell}\mu_{p}{\boldsymbol{c}}_{p}\in\mathbb{R}_{\geq}^{m}.$

This vector is a non-negative combination of ${\boldsymbol{c}}_{1}$,
${\boldsymbol{c}}_{2},\ldots,{\boldsymbol{c}}_{\ell}$; thus,
${\boldsymbol{c}}$ is a linear combination of the generating vectors. Since
$\mathbf{W}\in\mathbb{R}^{n\times|V|}$, we have
$\operatorname{rank}\mathbf{W}=|V|-\ell$.

Let $S$ be the associated linear space of $(V,E)$, i.e.,
$S=\mathrm{span}\\{{\boldsymbol{y}}_{j}-{\boldsymbol{y}}_{i}\colon{\boldsymbol{y}}_{i}\to{\boldsymbol{y}}_{j}\in
E\\}$. The falsity of the if statement in line 11 implies that
$\operatorname{im}\mathbf{W}\subseteq S$, so $\dim S\geq|V|-\ell$. Hence, the
deficiency of $(V,E)$ is $\delta=|V|-\ell-\dim S\leq 0$. Because $\delta$ is
always non-negative, we conclude that $\delta=0$. ∎

######  Lemma 3.5.

Suppose Algorithm 1 reaches line 23. Then $(V,E)$ is weakly reversible.

###### Proof.

Without loss of generality, we reorder the vertices according to the connected
components of the deficiency zero realization $(V,E,{\boldsymbol{\kappa}})$
found in line 23. Let ${\boldsymbol{y}}_{1}$,
${\boldsymbol{y}}_{2},\ldots,{\boldsymbol{y}}_{m_{1}}$ form the first
connected component;
${\boldsymbol{y}}_{m_{1}+1},\ldots,{\boldsymbol{y}}_{m_{2}}$ form the second
connected component, and so on. Accordingly, reorder the columns of the matrix
$\mathbf{W}$ of net direction vectors, the indices of the generators
${\boldsymbol{c}}_{p}$, and their supports $V_{p}$. Under the new ordering of
vertices,
$\mathrm{supp}({\boldsymbol{c}}_{p})=V_{p}=\\{{\boldsymbol{y}}_{m_{p-1}+1},\ldots,{\boldsymbol{y}}_{m_{p}}\\}$
for $p=1$, $2,\ldots,\ell$, where we take $m_{0}=0$.

With the new ordering of vertices, the Kirchoff matrix
$\mathbf{A}_{{\boldsymbol{\kappa}}}$ of $(V,E)$ is block-diagonal. Denote by
$\mathbf{A}_{{\boldsymbol{\kappa}}}^{(p)}$ the $p$th block, which encodes the
connectivity of the connected component defined by $V_{p}$. Indeed,
$\mathbf{A}_{{\boldsymbol{\kappa}}}^{(p)}$ is the Kirchoff matrix of the $p$th
connected component when viewed as a weighted E-graph in the subgraph sense.

We now prove that the first component $(V_{1},E_{1})$ is strongly connected.
Let $\mathbf{Y}^{(1)}$ be the first $m_{1}$ columns of $\mathbf{Y}_{s}$ and
$\mathbf{W}^{(1)}$ be the first $m_{1}$ columns of $\mathbf{W}$, so that
$\mathbf{Y}^{(1)}\mathbf{A}_{{\boldsymbol{\kappa}}}^{(1)}=\mathbf{W}^{(1)}$.
Let
${\boldsymbol{c}}_{1}={\boldsymbol{c}}_{1}^{\prime}\oplus{\boldsymbol{0}}\in\mathbb{R}_{\geq}^{n}$,
with ${\boldsymbol{c}}_{1}^{\prime}\in\mathbb{R}_{>}^{m_{1}}$ spanning the
one-dimensional subspace $\ker\mathbf{W}^{(1)}$. Finally, let
$S_{1}=\operatorname{span}\\{{\boldsymbol{y}}_{j}-{\boldsymbol{y}}_{i}\colon{\boldsymbol{y}}_{i}\in
V_{1}\\}$.

Because $(V_{1},E_{1})$ is connected and $V_{1}$ is affinely independent,
$\dim S_{1}=|V_{1}|-1$, so $\delta_{1}=0$. Moreover, if $t$ denotes the number
of terminal strongly connected components in $(V_{1},E_{1})$, then
$\dim(\ker\mathbf{Y}^{(1)}\cap\operatorname{im}\mathbf{A}_{{\boldsymbol{\kappa}}}^{(1)})=|V_{1}|-t-\dim
S_{1}$ [FeinbergHorn1977] [GunaNts]*Proposition 3.1. Since $|V_{1}|-t-\dim
S_{1}\leq\delta_{1}$, we conclude that
$\ker(\mathbf{Y}^{(1)}\mathbf{A}_{{\boldsymbol{\kappa}}}^{(1)})=\ker\mathbf{A}_{{\boldsymbol{\kappa}}}^{(1)}$,
and ${\boldsymbol{c}}_{1}^{\prime}$ also spans
$\ker\mathbf{A}_{{\boldsymbol{\kappa}}}^{(1)}$. By Theorem 2.14,
${\boldsymbol{c}}_{1}^{\prime}$ is supported on the terminal strongly
connected component, which in this case is all of $V_{1}$. Therefore,
$(V_{1},E_{1})$ is in fact strongly connected.

An analogous claim can be made about the other connected components.
Consequently $(V,E)$ is weakly reversible. ∎

######  Remark 3.6.

In the proof of Lemma 3.3 and Lemma 3.4, we have proved that the associated
linear subspace $S(V_{p})$ is in fact the span of the net direction vectors
belonging to said connected component. Thus, there are multiple ways of
generating $S(V_{p})$:

$\displaystyle\operatorname{span}\\{{\boldsymbol{y}}_{j}-{\boldsymbol{y}}_{i}\colon{\boldsymbol{y}}_{i},{\boldsymbol{y}}_{j}\in
V_{p}\\}=\operatorname{span}\\{{\boldsymbol{y}}_{j}-{\boldsymbol{y}}_{i}\colon{\boldsymbol{y}}_{i}\to{\boldsymbol{y}}_{j}\in
E_{p}\\}=\operatorname{im}\mathbf{W}^{(p)}.$

The lemmas above provide the technical parts that we need to prove the main
result of this paper.

######  Theorem 3.7.

Given a system of differential equations

$\displaystyle\frac{d{\boldsymbol{x}}}{dt}=\sum_{i=1}^{m}{\boldsymbol{x}}^{{\boldsymbol{y}}_{i}}{\boldsymbol{w}}_{i},$

with distinct ${\boldsymbol{y}}_{i}\in\mathbb{Z}_{\geq}^{n}$, and
${\boldsymbol{w}}_{i}\in\mathbb{R}^{n}\setminus\\{{\boldsymbol{0}}\\}$.
Algorithm 1 returns the unique WR0 realization of the dynamical system if it
exists, or concludes that no WR0 realization exists.

###### Proof.

There are two possible scenarios: either (1) the algorithm exits at lines 3,
8, or 12 by failing one of the if statements, or (2) the algorithm
successfully reaches line 23. In the first scenario, Lemma 3.2 implies that no
WR0 realization exists. In the second scenario, the realization has connected
components $V_{1}$, $V_{2},\ldots,V_{\ell}$ according to Lemma 3.3. The
realization is weakly reversible and deficiency zero by Lemmas 3.4 and 3.5
respectively. The uniqueness of the realization follows from
[CraciunJinYu_unique]. ∎

######  Remark 3.8.

The uniqueness of the WR0 realization is also a consequence of Algorithm 1.
This is due to the affine and linear independences, as well as the structure
of $\ker\mathbf{W}=\ker\mathbf{A}_{{\boldsymbol{\kappa}}}$.

If a WR0 realization exists, the polynomial dynamical system is complex-
balanced. Therefore, if a system passes Algorithm 1, it automatically inherits
all the algebraic and dynamical properties of complex-balanced system. Weak
reversibility implies that a positive steady state exists [Boros2019]. The
remaining statements in the theorem below are easy consequence of Theorems
2.12 and 3.7.

######  Theorem 3.9.

Suppose the system of differential equations

$\displaystyle\frac{d{\boldsymbol{x}}}{dt}=\sum_{i=1}^{m}{\boldsymbol{x}}^{{\boldsymbol{y}}_{i}}{\boldsymbol{w}}_{i},$
(7)

with distinct ${\boldsymbol{y}}_{i}\in\mathbb{Z}_{\geq}^{n}$ and
${\boldsymbol{w}}_{i}\in\mathbb{R}^{n}\setminus\\{{\boldsymbol{0}}\\}$, passes
Algorithm 1. Let $\mathbf{W}$ be the matrix of net direction vectors and
$S=\operatorname{im}\mathbf{W}$. Then the following holds.

  1. (i)

A positive steady state ${\boldsymbol{x}}^{*}$ exists.

  2. (ii)

There is exactly one steady state within every invariant polyhedron
$({\boldsymbol{x}}_{0}+S)\cap\mathbb{R}_{>}^{n}$ for any
${\boldsymbol{x}}_{0}\in\mathbb{R}_{>}^{n}$, and it is complex-balanced.

  3. (iii)

Any positive steady state ${\boldsymbol{x}}$ satisfies
$\ln{\boldsymbol{x}}-\ln{\boldsymbol{x}}^{*}\in S^{\perp}$.

  4. (iv)

The function

$\displaystyle L({\boldsymbol{x}})=\sum_{i=1}^{n}x_{i}(\ln x_{i}-\ln
x^{*}_{i}-1),$

defined on $\mathbb{R}_{>}^{n}$, is a strict Lyapunov function of (7) within
every invariant polyhedron $({\boldsymbol{x}}_{0}+S)\cap\mathbb{R}_{>}^{n}$,
with a global minimum at the corresponding complex-balanced steady state.

  5. (v)

Every positive steady state is locally asymptotically stable with respect to
its invariant polyhedron.

######  Example 3.10.

Consider the system of differential equations

$\begin{split}\frac{dx_{1}}{dt}&=-12x_{1}+x^{2}_{3},\\\
\frac{dx_{2}}{dt}&=14x_{1}-4x^{2}_{2}+8x^{2}_{3},\\\
\frac{dx_{3}}{dt}&=10x_{1}+4x^{2}_{2}-10x^{2}_{3}.\end{split}$ (8)

We have $n=3$ for the three state variables, and $m=3$ for the three distinct
monomials. The matrices of source vertices and net direction vectors are

$\displaystyle\mathbf{Y}_{s}$
$\displaystyle=\begin{pmatrix}{\boldsymbol{y}}_{1}&{\boldsymbol{y}}_{2}&{\boldsymbol{y}}_{3}\end{pmatrix}=\begin{pmatrix}1&0&0\\\
0&2&0\\\ 0&0&2\end{pmatrix},$ $\displaystyle\mathbf{W}$
$\displaystyle=\begin{pmatrix}{\boldsymbol{w}}_{1}&{\boldsymbol{w}}_{2}&{\boldsymbol{w}}_{3}\end{pmatrix}=\begin{pmatrix}-12&\hphantom{-}0&\hphantom{-}1\\\
\hphantom{-}14&-4&\hphantom{-}8\\\
\hphantom{-}10&\hphantom{-}4&-10\end{pmatrix}$

respectively, which are inputs to Algorithm 1. Let
$V=\\{{\boldsymbol{y}}_{1},{\boldsymbol{y}}_{2},{\boldsymbol{y}}_{3}\\}\subset\mathbb{Z}_{\geq}^{n}$.
A generator for the cone $\ker\mathbf{W}\cap\mathbb{R}_{\geq}^{m}$ is
${\boldsymbol{c}}=(48/1441,120/131,576/1441)^{\top}$. In the notation of
Algorithm 1, $V_{1}=[1,3]$. Clearly $V$ is affinely independent and the net
direction vectors admit the following unique decompositions:

$\displaystyle{\boldsymbol{w}}_{1}$
$\displaystyle=7({\boldsymbol{y}}_{2}-{\boldsymbol{y}}_{1})+5({\boldsymbol{y}}_{3}-{\boldsymbol{y}}_{1}),$
$\displaystyle{\boldsymbol{w}}_{2}$
$\displaystyle=2({\boldsymbol{y}}_{3}-{\boldsymbol{y}}_{2}),$
$\displaystyle{\boldsymbol{w}}_{3}$
$\displaystyle=({\boldsymbol{y}}_{1}-{\boldsymbol{y}}_{3})+4({\boldsymbol{y}}_{2}-{\boldsymbol{y}}_{3}).$

Therefore (8) admits a WR0 realization, whose weighted E-graph is shown in
Figure 4(b).

This implies that the system (8) has exactly one steady state within each
invariant triangle given by $2x_{1}+x_{2}+x_{3}=C$ for some $C>0$, and this
steady state is a global attractor within each such triangle. From Theorem
3.9, we know the steady state set admits a monomial parametrization of the
form $(a_{1}s^{2},a_{2}s,a_{3}s)$ for some constants $a_{i}>0$. In fact, the
set of steady states is given by

$\displaystyle(x_{1}^{*},x_{2}^{*},x_{3}^{*})=\left(3s^{2},\,\frac{\sqrt{330}}{2}s,\,6s\right),$

and an explanation for the coefficients above will be provided in Theorem
3.12.

${\boldsymbol{y}}_{1}$
${\boldsymbol{y}}_{2}$${\boldsymbol{y}}_{3}$$18$$4$$10$ (a)

${\boldsymbol{y}}_{1}$
${\boldsymbol{y}}_{2}$${\boldsymbol{y}}_{3}$$7$$5$$1$$4$$2$ (b)

${\boldsymbol{y}}_{1}$
${\boldsymbol{y}}_{2}$${\boldsymbol{y}}_{3}$$4$$4$$10$ (c)

Figure 4: (a) A weighted E-graph realizing (8) from Example 3.10, which admits
(b) a WR0 realization. (c) A weighted E-graph realizing (9) from Example 3.11
that does not admit a WR0 realization.

######  Example 3.11.

Consider the system of differential equations

$\begin{split}\frac{dx_{1}}{dt}&=-\frac{1}{2}x_{1}+x^{2}_{3},\\\
\frac{dx_{2}}{dt}&=-2x_{1}-4x^{2}_{2}+8x^{2}_{3},\\\
\frac{dx_{3}}{dt}&=3x_{1}+4x^{2}_{2}-10x^{2}_{3}.\end{split}$ (9)

Again, we have $n=3$ and $m=3$. The monomials are the same as those in the
previous example. The difference lies in the first column of the matrix of net
direction vectors

$\displaystyle\mathbf{W}$
$\displaystyle=\begin{pmatrix}{\boldsymbol{w}}_{1}&{\boldsymbol{w}}_{2}&{\boldsymbol{w}}_{3}\end{pmatrix}=\begin{pmatrix}-\frac{1}{2}&\hphantom{-}0&\hphantom{-}1\\\
-2&-4&\hphantom{-}8\\\ \hphantom{-}3&\hphantom{-}4&-10\end{pmatrix},$

whose kernel is spanned by ${\boldsymbol{c}}=(2,1,1)^{\top}$. As in the
previous example, the vertices ${\boldsymbol{y}}_{1}$, ${\boldsymbol{y}}_{2}$,
and ${\boldsymbol{y}}_{3}$ are affinely independent. However,
${\boldsymbol{w}}_{1}\not\in\mathrm{Cone}\\{{\boldsymbol{y}}_{j}-{\boldsymbol{y}}_{1}\colon
j=2,3\\}$, so no WR0 realization exists.

###  3.2 The set of positive steady states of a WR0 realization

Algorithm 1 determines whether a given polynomial dynamical system admits a
WR0 realization. If it does, its steady state set is in fact log-linear. In
this section, we write down a system of linear equations whose solution set is
in bijection with the set of positive steady states; this provides an explicit
parametrization of the set of positive steady states.

For any ${\boldsymbol{z}}\in\mathbb{R}^{n}$ and
${\boldsymbol{x}}\in\mathbb{R}_{>}^{n}$, define the component-wise operations
$\exp{\boldsymbol{z}}=(e^{z_{1}},e^{z_{2}},\ldots,e^{z_{n}})^{\top}$ and
$\log({\boldsymbol{x}})=(\log x_{1},\log x_{2},\ldots,\log x_{n})^{\top}$. We
extend these operations to sets. If $Z\subseteq\mathbb{R}^{n}$, then
$\exp(Z)=\\{\exp{\boldsymbol{z}}\colon{\boldsymbol{z}}\in Z\\}$, and if
$X\subseteq\mathbb{R}_{>}^{n}$, then
$\log(X)=\\{\log{\boldsymbol{x}}\colon{\boldsymbol{x}}\in X\\}$.

Assume that the polynomial dynamical system

$\displaystyle\frac{d{\boldsymbol{x}}}{dt}=\sum_{i=1}^{m}{\boldsymbol{x}}^{{\boldsymbol{y}}_{i}}{\boldsymbol{w}}_{i},$
(10)

with distinct ${\boldsymbol{y}}_{i}\in\mathbb{Z}_{\geq}^{n}$ and
${\boldsymbol{w}}_{i}\in\mathbb{R}^{n}\setminus\\{{\boldsymbol{0}}\\}$, passes
Algorithm 1, i.e., it admits a WR0 realization $(V,E,{\boldsymbol{\kappa}})$.
Without loss of generality, assume the vertices are ordered according to
connected components in $(V,E)$, i.e., the first $m_{1}$ vertices belong to
the connected component $(V_{1},E_{1})$, the next $m_{2}$ vertices belong to
the connected component $(V_{2},E_{2})$, and so forth. Let
$\\{{\boldsymbol{c}}_{1},{\boldsymbol{c}}_{2},\ldots,{\boldsymbol{c}}_{\ell}\\}$
be a minimal set of generators of $\ker\mathbf{W}\cap\mathbb{R}_{\geq}^{m}$,
ordered in an analogous way. From Algorithm 1, we know that the supports of
the vectors ${\boldsymbol{c}}_{1}$,
${\boldsymbol{c}}_{2},\ldots,{\boldsymbol{c}}_{\ell}$ correspond to the
connected components of $(V,E)$.

Let
${\boldsymbol{c}}_{1}=(\alpha_{1},\alpha_{2},\ldots,\alpha_{m_{1}},0,\ldots,0)^{\top}$.
Define matrix $\mathbf{D}_{1}\in\mathbb{R}^{(m_{1}-1)\times n}$ whose rows are
the affine vectors from ${\boldsymbol{y}}_{1}$ to the remaining vertices of
$V_{1}$, and define vector ${\boldsymbol{J}}_{1}\in\mathbb{R}^{m_{1}-1}$ using
the log-differences of the components of ${\boldsymbol{c}}_{1}$, i.e.,

$\displaystyle\mathbf{D}_{1}=\begin{pmatrix}{\boldsymbol{y}}_{2}-{\boldsymbol{y}}_{1}\\\
{\boldsymbol{y}}_{3}-{\boldsymbol{y}}_{1}\\\ \vdots\\\
{\boldsymbol{y}}_{m_{1}}-{\boldsymbol{y}}_{1}\end{pmatrix}\quad\text{and}\quad{\boldsymbol{J}}_{1}=\begin{pmatrix}\log(\alpha_{2}/\alpha_{1})\\\
\log(\alpha_{3}/\alpha_{1})\\\ \vdots\\\
\log(\alpha_{m_{1}}/\alpha_{1})\end{pmatrix}.$

For the connected component $(V_{p},E_{p})$, define $\mathbf{D}_{p}$ and
${\boldsymbol{J}}_{p}$ in a similar fashion. Define

$\displaystyle\mathbf{D}=\begin{pmatrix}\mathbf{D}_{1}\\\
\hline\cr\mathbf{D}_{2}\\\ \hline\cr\vdots\\\
\hline\cr\mathbf{D}_{\ell}\end{pmatrix}\in\mathbb{R}^{(m-\ell)\times
n}\quad\text{and}\quad{\boldsymbol{J}}=\begin{pmatrix}{\boldsymbol{J}}_{1}\\\
\hline\cr{\boldsymbol{J}}_{2}\\\ \hline\cr\vdots\\\
\hline\cr{\boldsymbol{J}}_{\ell}\end{pmatrix}\in\mathbb{R}^{m-\ell}.$ (11)

######  Theorem 3.12.

Suppose the system of differential equation (10) admits a WR0 realization
$(V,E,{\boldsymbol{\kappa}})$, and let
$\mathbf{D}\in\mathbb{R}^{(m-\ell)\times n}$ and
${\boldsymbol{J}}\in\mathbb{R}_{>}^{m-\ell}$ be defined as in (11). Then the
system $\mathbf{D}{\boldsymbol{z}}={\boldsymbol{J}}$ is solvable. Let
${\boldsymbol{z}}^{*}+\ker\mathbf{D}$ be its solution set. Then the set of
positive steady states of (10) is $\exp({\boldsymbol{z}}^{*}+\ker\mathbf{D})$.

###### Proof.

First we prove that the linear system
$\mathbf{D}{\boldsymbol{z}}={\boldsymbol{J}}$ is solvable. Consider
$\mathbf{D}_{1}$. The vertices ${\boldsymbol{y}}_{1}$,
${\boldsymbol{y}}_{2},\ldots,{\boldsymbol{y}}_{m_{1}}$ in the first connected
component are affinely independent, so the rows of $\mathbf{D}_{1}$ are
linearly independent. Moreover, as noted in Remark 3.6 the row-space of
$\mathbf{D}_{1}$ is the associated linear subspace $S(V_{1})$. Therefore
$\operatorname{rank}\mathbf{D}_{1}=m_{1}-1$, and the matrix $\mathbf{D}_{1}$
is surjective onto $\mathbb{R}^{m_{1}-1}$.

Similarly, for each $p=2,\ldots,\ell$, the row-space of the matrix
$\mathbf{D}_{p}$ is $S(V_{p})$, and the matrix $\mathbf{D}_{p}$ is surjective.
In addition, because the realization $(V,E,{\boldsymbol{\kappa}})$ has
deficiency zero, $S(V_{1})$, $S(V_{2}),\ldots,S(V_{\ell})$ are linearly
independent; in other words, the $m-\ell$ rows of the matrix $\mathbf{D}$ are
linearly independent. Consequently, $\mathbf{D}$ is surjective, and the system
$\mathbf{D}{\boldsymbol{z}}={\boldsymbol{J}}$ is solvable.

Let ${\boldsymbol{z}}^{*}+\ker\mathbf{D}$ be the set of solution to
$\mathbf{D}{\boldsymbol{z}}={\boldsymbol{J}}$. We next show that each solution
can be related to a positive steady state of (10), which by definition
satisfies

$\displaystyle{\boldsymbol{0}}=\sum_{i=1}^{m}{\boldsymbol{x}}^{{\boldsymbol{y}}_{i}}{\boldsymbol{w}}_{i}.$

In other words,
$({\boldsymbol{x}}^{{\boldsymbol{y}}_{1}},\ldots,{\boldsymbol{x}}^{{\boldsymbol{y}}_{m}})^{\top}$
lies in the steady state flux cone $\ker\mathbf{W}\cap\mathbb{R}_{>}^{m}$.
Decomposing this vector with respect to the generators of the cone allows us
to focus on one connected component at a time.

For simplicity of notation, consider the first connected component. At steady
state, for some constant $\lambda>0$, we have
${\boldsymbol{x}}^{{\boldsymbol{y}}_{j}}=\lambda\alpha_{j}$ for $j=1$,
$2,\ldots,m_{1}$. Thus

$\displaystyle{\boldsymbol{x}}^{{\boldsymbol{y}}_{j}-{\boldsymbol{y}}_{1}}=\frac{\alpha_{j}}{\alpha_{1}}$

for $j=2,3,\ldots,m_{1}$. Taking the logarithm of both sides, we obtain the
system $\mathbf{D}_{1}{\boldsymbol{z}}={\boldsymbol{J}}_{1}$ with
${\boldsymbol{z}}=\log{\boldsymbol{x}}$.

Repeating this computation for each connected component, we conclude that
${\boldsymbol{x}}$ is a positive steady state for (10) if and only if
${\boldsymbol{x}}$ solves $\mathbf{D}{\boldsymbol{z}}={\boldsymbol{J}}$ with
${\boldsymbol{z}}=\log{\boldsymbol{x}}$. This leads us to the characterization
of the set of positive steady states for (10) as
$\exp({\boldsymbol{z}}^{*}+\ker\mathbf{D})$, where
${\boldsymbol{z}}^{*}+\ker\mathbf{D}$ is the set of solutions to
$\mathbf{D}{\boldsymbol{z}}={\boldsymbol{J}}$. ∎

###  3.3 Extension to polynomial systems with unspecified coefficients

If instead of (10), we need to analyze

$\displaystyle\frac{d{\boldsymbol{x}}}{dt}=\sum_{i=1}^{m}a_{i}{\boldsymbol{x}}^{{\boldsymbol{y}}_{i}}{\boldsymbol{w}}_{i}$
(12)

for some unknown $a_{i}>0$, it turns out that the answer as to whether a WR0
realization exists is the same:

######  Theorem 3.13.

For any $a_{i}>0$, the system (12) admits a WR0 realization
$(V,E,{\boldsymbol{\kappa}})$ if and only if the system (10) admits a WR0
realization $(V,E,{\boldsymbol{\kappa}}^{*})$. Moreover,
$\kappa_{ij}=a_{i}\kappa^{*}_{ij}$.

###### Proof.

The forward implication is trivial. We focus our attention on the other
direction. For any $i$, $j$, let $\kappa_{ij}=a_{i}\kappa^{*}_{ij}$, so
$\kappa_{ij}>0$ if and only if $\kappa^{*}_{ij}>0$. In other words, the
weighted E-graph $(V,E,{\boldsymbol{\kappa}})$ shares the same set of edges as
$(V,E,{\boldsymbol{\kappa}}^{*})$. Because the deficiency is characterized by
affine and linear independence of the connected components, and the two graphs
share the same structure, $(V,E,{\boldsymbol{\kappa}})$ is weakly reversible
and deficiency zero if and only if $(V,E,{\boldsymbol{\kappa}}^{*})$ is.

Suppose $(V,E,{\boldsymbol{\kappa}}^{*})$ is a realization of (10). Then in
$(V,E,{\boldsymbol{\kappa}})$, the net direction vector from
${\boldsymbol{y}}_{i}$ can be expanded using the realization
$(V,E,{\boldsymbol{\kappa}}^{*})$, since

$\displaystyle a_{i}{\boldsymbol{w}}_{i}=a_{i}\sum_{(i,j)\in
E}\kappa^{*}_{ij}({\boldsymbol{y}}_{j}-{\boldsymbol{y}}_{i})=\sum_{(i,j)\in
E}\kappa_{ij}({\boldsymbol{y}}_{j}-{\boldsymbol{y}}_{i}).$

Therefore, $(V,E,{\boldsymbol{\kappa}})$ realizes (12). ∎

###  3.4 Deficiency zero realizations that are not weakly reversible

If a polynomial dynamical system admits a deficiency zero realization that is
not weakly reversible, then its dynamics is also greatly restricted: it can
have no positive steady states, no oscillations, and no chaotic dynamics
[horn1972necessary, Feinberg1972, feinberg2019foundations]. Actually, such
realizations are special examples of mass-action system that are not
consistent [AndersonBrunnerCraciunJohnston2020]. An E-graph $(V,E)$ is said to
be consistent if there exist real numbers $\alpha_{ij}>0$ such that

$\displaystyle\sum_{(i,j)\in
E}\alpha_{ij}({\boldsymbol{y}}_{j}-{\boldsymbol{y}}_{i})={\boldsymbol{0}}.$

It is easy to see that a polynomial dynamical system of the form (10) has a
realization that is not consistent if and only if

$\ker\mathbf{W}\cap\mathbb{R}^{m}_{>}=\emptyset.$ (13)

If a polynomial dynamical system has a realization that is not consistent,
then it cannot have realization that is weakly reversible, because weakly
reversible systems must have at least one positive steady state [Boros2019].
Therefore, if Algorithm 1 is accompanied by a preprocessing step that checks
condition (13), then that step will decide whether our given system (10) has a
realization that is not consistent; in particular, this step will also find
all cases where our given system has a deficiency zero realization that is not
weakly reversible.

## Acknowledgements

The authors were supported in part by the National Science Foundation under
grants DMS–1816238 and DMS–2051568. G.C. was also partially supported by the
Simons Foundation.

## References