Weakly reversible deficiency one realizations of polynomial dynamical systems: an algorithmic perspective

Gheorghe Craciun

Abhishek Deshpande

Jiaxin Jin

[1]Department of Mathematics and Department of Biomolecular Chemistry, University of Wisconsin-Madison
[2]Center for Computational Natural Sciences and Bioinformatics,

International Institute of Information Technology Hyderabad
[3]Department of Mathematics, The Ohio State University

Given a dynamical system with polynomial right-hand side, can it be generated by a reaction network that possesses certain properties? This question is important because some network properties may guarantee specific dynamical properties, such as existence or uniqueness of equilibria, persistence, permanence, or global stability. Here we focus on this problem in the context of weakly reversible deficiency one networks. In particular, we describe an algorithm for deciding if a polynomial dynamical system admits a weakly reversible deficiency one realization, and identifying one if it does exist. In addition, we show that weakly reversible deficiency one realizations can be partitioned into mutually exclusive Type I and Type II realizations, where Type I realizations guarantee existence and uniqueness of positive steady states, while Type II realizations are related to stoichiometric generators, and therefore to multistability.

§ INTRODUCTION

By a  polynomial dynamical system we mean a dynamical system of the form

\begin{equation}\label{eq:poly-intro}
\begin{split}
\frac{dx_1}{dt} &= p_1(x_1, \ldots, x_n), \\
\frac{dx_2}{dt} &= p_2(x_1, \ldots, x_n), \\
&\qquad \quad \vdots  \\
\frac{dx_n}{dt} &= p_n(x_1, \ldots, x_n), \\
\end{split}
\end{equation}

where each $p_i(x_1,\ldots, x_n)$ is a polynomial in the variables  $x_1,\ldots, x_n$.
Such systems can exhibit exotic behaviors like multistability, presence of oscillations, and chaos due to the underlying nonlinearities.
We are especially interested in the dynamics of these systems when restricted to the positive orthant, because such systems are very common models of biological interaction networks, population dynamics models, or models for the transmission of infectious diseases [1, 2, 3, 4].

Very often,  polynomial dynamical systems are generated by reaction networks. It is often convenient to study the graphical structure of these networks to make inference about their dynamics. It is also possible to study the inverse problem, i.e., for some given polynomial dynamical systems, ask what reaction networks can generate them. Due to the phenomenon of dynamical equivalence [5, 6], such a network may not be unique, i.e., there  exist multiple reaction networks that  generate the same dynamics.

A key quantity in the study of these networks is its deficiency. In particular, networks possessing low deficiency have been studied in reaction network theory using the Deficiency Zero and Deficiency One theorems [2, 6, 7, 8, 9, 10]. In particular, the Deficiency One Theorem [2, 10] guarantees  uniqueness of the steady state within each linear invariant subspace; this, together with the existence result in [11] completely characterizes the steady states of weakly reversible networks that satisfy the hypotheses of the  Deficiency One Theorem.  Further, weakly reversible networks (i.e., networks where each reaction is part of a cycle) are related to dynamical properties like persistence, permanence, and the existence of a globally attracting steady state [12, 13, 14, 15].

The problem of identifying weakly reversible deficiency zero realizations has been addressed in [16]. Here we analyze the realizability problem for weakly reversible reaction networks with deficiency one. In particular, given a polynomial dynamical system, we describe an algorithm to identify if there exists a weakly reversible deficiency one reaction network that generates this dynamical system.

Moreover, if  weakly reversible deficiency one realizations do exist, our algorithm uses the geometry of some convex cones generated using the net reaction vectors to construct one such realization explicitly.

Structure of this paper.
In Section <ref>, we recall some basic notions from reaction network theory.
Primarily, we introduce dynamical equivalence and the matrix of net reaction vectors.
In Section <ref>, we give a short primer on weakly reversible deficiency one reaction networks, and define Type I and Type II weakly reversible deficiency one realizations.
In Section <ref>, we analyze the pointed cone $\text{ker}(\bW)\cap\mathbb{R}^m_{\geq 0}$ and its minimal set of generators in the context of weakly reversible deficiency one networks.
Moreover, we prove there cannot exist a dynamical equivalence between such networks of two types in Theorem <ref>.
In Section <ref>, we state the main algorithm of our paper: Algorithm <ref>, which checks the existence of a weakly reversible deficiency one realization and returns a realization if it exists.
In Section <ref>, we summarize our findings in this paper and flesh out directions for future work.

We let $\mathbb{R}_{\geq 0}^n$ and $\mathbb{R}_{>0}^n$ denote the set of vectors in $\mathbb{R}^n$ with non-negative and positive entries respectively.
Given two vectors $\bx = (\bx_1, \ldots, \bx_n)^{\intercal}\in \RR^n_{>0}$ and $\by = (\by_1, \ldots, \by_n)^{\intercal} \in \RR^n$, we use the vector operation as follows:

\begin{equation} \notag
\bx^{\by} := \bx_1^{y_{1}} \ldots \bx_n^{y_{n}}.
\end{equation}

Given a positive integer $m$, we denote $[m] := \{1, \ldots, m \}$.

§ REACTION NETWORKS

§.§ Terminology

A reaction network $G = (V, E)$, also called the Euclidean embedded graph (E-graph), is a directed graph in $\RR^n$, where $V \subsetneq \mathbb{R}^n$ represents a finite set of vertices, and $E \subseteq V \times V$ represents a finite set of edges.
In this paper, there are neither self-loops nor isolated vertices in $G$.

* Let $V = \{ \by_1, \ldots, \by_m \}$, and denote the number of vertices in $G$ by $m$.

* We denote a directed edge by $(\by_i, \by_j) \in E$, which represents a reaction in the network. Here $\by_i$ and $\by_j$ are called the source vertex and target vertex respectively.
Moreover, we denote the reaction vector associated with the edge $\by_i \rightarrow \by_j$ by $\by_j - \by_i \in\mathbb{R}^n$.

Let $G=(V, E)$ be a Euclidean embedded graph.
The stoichiometric subspace of $G$ is the vector space spanned by the reaction vectors as follows:

\begin{equation} \notag
S = \text{span} \{\by' - \by\, |\, \by \rightarrow \by' \in E \}.
\end{equation}

Given a subset of vertices $V_0 \subseteq V$, the stoichiometric subspace defined by $V_0$ is

\begin{equation} \notag
S(V_0) = \text{span} \{ \by' - \by\, |\, \by \rightarrow \by' \in E \ \text{and} \ \by', \by \in V_0 \}.
\end{equation}

Furthermore, given a positive vector $\bx_0 \in\mathbb{R}_{>0}^n$, the polyhedron $(\bx_0 + S ) \cap \mathbb{R}^n_{>0 }$ is called the stoichiometric compatibility class of $\bx_0$.

Let $G=(V, E)$ be a Euclidean embedded graph.

* The set of vertices $V$ is partitioned by its connected components, also called linkage classes.
Every connected component is denoted by the set of vertices belonging to it.

* A connected component $L \subseteq V$ is strongly connected, if every edge is part of an oriented cycle.
Moreover, a strongly connected component $L \subseteq V$ is terminal strongly connected, if for every vertex $\by\in L$ and $\by\rightarrow\by'\in E$, we have $\by'\in L$.

* $G=(V,E)$ is weakly reversible, if all connected components are strongly connected.

For any weakly reversible reaction network $G=(V, E)$, every vertex $\by \in V$ is a source and a target vertex.
Moreover, every connected component of $G$ is strongly connected, and terminal strongly connected.

Let $G = (V, E)$ be a Euclidean embedded graph, which contains $m$ vertices and $\ell$ connected components. Denote the dimension of the stoichiometric subspace of $G$ by $s = \dim (S)$, the deficiency of $G$ is a non-negative integer as follows:

\begin{equation*}
\delta = m - \ell - s.
\end{equation*}

When considering a linkage class with $V_i \subseteq V$, we define the deficiency of a linkage class as

\begin{equation*}
\delta_i = |V_i| - 1 - \dim S(V_i).
\end{equation*}

One can check that

\begin{equation} \label{delta inequality}
\delta \geq \sum_{i=1}^\ell \delta_i,
\end{equation}

where the equality holds when the stoichiometric subspaces of all linkage classes are linearly independent.

(a) This reaction network consists of two linkage classes and contains three terminal strongly connected components (shown in circles).  It has a stoichiometric subspace of dimension 2 and a deficiency $\delta = n - \ell -s = 7- 2-2=3$. (b) This reaction network is weakly reversible and contains one terminal strongly connected component. It has a stoichiometric subspace of dimension 2 and a deficiency $\delta = n - \ell -s = 3- 1-2=0$.

Figure <ref> shows two reaction networks represented as Euclidean embedded graphs.
Given a reaction network, it can generate an extensive variety  of dynamical systems. Here, we focus on mass-action kinetics, which has been studied in [1, 3, 4, 17, 18, 19].

Let $G=(V,E)$ be a Euclidean embedded graph, we
let $\bk = (k_{\by_i \rightarrow \by_j})_{\by_i \rightarrow \by_j \in E} \in \mathbb{R}_{>0}^{E}$
denote the vector of reaction rate constants, where $k_{\by_i \rightarrow \by_j}$ or $k_{ij}$
is called the reaction rate constant of the edge $\by_i \rightarrow \by_j \in E$. Furthermore, the associated mass-action system generated by $(G, \bk)$
on $\RR^n_{>0}$ is

\begin{equation} \label{eq:mass_action}
\frac{\mathrm{d} \bx}{\mathrm{d} t}= \sum_{\by_i \rightarrow \by_j \in E}k_{\by_i \rightarrow \by_j} \bx^{\by_i}(\by_j - \by_i).
\end{equation}

A point $\bx^* \in \mathbb{R}_{>0}^n$ is called a positive steady state, if

\begin{equation} \label{eq:ss}
%\frac{\mathrm{d} \bx}{\mathrm{d} t}
\sum_{\by_i \rightarrow \by_j \in E}k_{\by_i \rightarrow \by_j} (\bx^*)^{\by_i}(\by_j - \by_i)
= \mathbf{0}.
\end{equation}

From [6], it is known that every mass-action system admits the following matrix decomposition:

\begin{equation} \label{eq:mass_action_horn_jackson}
\frac{d\bx}{dt} = \bY \bA_{\bk}{\bx}^{\bY},
\end{equation}

where $\bY$ is called the matrix of vertices, whose columns are the vertices, defined as

\begin{equation*}
\bY = (\by_1, \ \by_2, \ \ldots, \ \by_m),
\end{equation*}

${\bx}^{\bY}$ is the vector
of monomials given by

\begin{equation*}
{\bx}^{\bY} = ({\bx}^{\by_1},{\bx}^{\by_2}, \ldots, {\bx}^{\by_m})^{\intercal},
\end{equation*}

and $\bA_{\bk}$ is the negative transpose of the Laplacian of $(G, \bk)$, defined as

\begin{equation*}
[\bA_{\bk}]_{ji} =
\begin{cases}
k_{\by_i\rightarrow\by_j}, & \text{if } \by_i\rightarrow\by_j \in E, \\[5pt]
- \sum\limits_{\by_i\rightarrow\by_j \in E} k_{\by_i\rightarrow\by_j}, & \text{if } i=j, \\[5pt]
0,  & \text{otherwise}.
\end{cases}
\end{equation*}

Here, $\bA_{\bk}$ is called the Kirchoff matrix, whose column sums are zero from the definition.

The properties of the kernel of $\bA_{\bk}$ are well known in reaction network theory [2, 3, 20].
Below we collect some of the most important properties.

Let $(G, \bk)$ be a mass-action system, and $T_1, T_2, \ldots, T_t$ be the terminal strongly connected components of $G$. Then there exists a basis $\{\bc_1, \bc_2, \ldots, \bc_t\}$ for $\ker (\bA_{\bk})$, such that

\begin{equation*}
\bc_q =
\begin{cases}
\begin{array}{cl}
[\bc_q]_i  > 0, & \text{ if } \by_i \in T_q, \\[5pt]
[\bc_q]_i  = 0, & \text{ otherwise.}
\end{array}
\end{cases}
\end{equation*}

Consider a mass-action system $(G, \bk)$, then $G$ is weakly reversible if and only if the kernel of the Kirchoff matrix $\bA_{\bk}$ contains a positive vector.

§.§ Net reaction vectors and dynamical equivalence

Inspired by the matrix decomposition in (<ref>), we introduce the key concept: net reaction vector, and illustrate a new matrix decomposition in terms of net reaction vectors.

Let $(G, \bk)$ be a mass-action system, and $V_s = \{ \by_1, \by_2, \ldots, \by_{m_s}\} \subseteq V$ be the set of source vertices of $G$. For each source vertex $\by_i \in V_s$, we denote the net reaction vector $\bw_i$ corresponding to $\by_i$ by

\begin{equation}
\bw_i = \sum\limits_{\by_i\rightarrow\by_j\in E}k_{\by_i\rightarrow\by_j}(\by_j - \by_i).
\end{equation}

Further, we denote the matrix of net reaction vectors of $G$ as follows:

\begin{equation}
\bW = \left(\bw_1, \bw_2, \ldots, \bw_{m_s} \right).
\end{equation}

It is convenient to refer to $\bw_j$ even when $\by_j \not\in V_s$, in which case we consider $\bw_j$ represents an empty sum, i.e., $\bw_j = \mathbf{0}$.

From Definition <ref>, every net reaction vector $\bw_i$ corresponding to $\by_i$ can be expressed as

\begin{equation}
\bw_i = \sum\limits_{\by_i\rightarrow\by_j\in E}k_{\by_i\rightarrow\by_j} \by_j
- \bigg( \sum\limits_{\by_i\rightarrow\by_j\in E}k_{\by_i\rightarrow\by_j} \bigg) \by_i.
\end{equation}

Using a direct computation, we rewrite the mass-action system in (<ref>) as

\begin{equation} \label{eq:mass_action_net_reaction}
\frac{d\bx}{dt} = \bW {\bx}^{\bY_s},
\end{equation}

where $\bY_s$ is called the matrix of source vertices, whose columns are source vertices, defined as

\begin{equation*}
\bY_s = (\by_1, \by_2, \ldots, \by_{m_s}),
\end{equation*}

and ${\bx}^{\bY_s}$ is the vector
of monomials given by

\begin{equation*}
{\bx}^{\bY_s} = ({\bx}^{\by_1},{\bx}^{\by_2}, \ldots, {\bx}^{\by_{m_s}})^{\intercal}.
\end{equation*}

Note that for any weakly reversible mass-action system $(G, \bk)$, we have $V_s = V$ and $\bY_s = \bY$.
Moreover, we derive that $\bW = \bY \bA_{\bk}$, which follows from matrix decomposition in (<ref>).

Let $(G,\bk)$ and $(\bar{G},\bar{\bk})$ be two mass-action systems. Then $(G,\bk)$ and $(\bar{G},\bar{\bk})$ are called dynamically equivalent, if for any $\bx \in \RR^n_{>0}$,

\begin{equation} \label{eq:eqDE}
\sum\limits_{\by\rightarrow \by'\in E}k_{\by\rightarrow \by'}{\bx}^{\by}(\by' - \by) =  \sum\limits_{\bar{\by}\rightarrow \bar{\by}'\in \bar{E}}\bar{k}_{\bar{\by}\rightarrow \bar{\by}'}{\bx}^{\bar{\by}}(\bar{\by}' - \bar{\by}).
\end{equation}

From Equation (<ref>), we achieve a necessary and sufficient condition for dynamical equivalence
between $(G,\bk)$ and $(\bar{G},\bar{\bk})$:
for every vertex $\by_0 \in V_s \cup \bar{V}_s$,

\begin{equation} \label{eq:DE}
\sum_{\by_0 \rightarrow \by \in E} k_{\by_0 \rightarrow \by} (\by - \by_0) =
\sum_{\by_0 \rightarrow \by' \in \bar{E}} \bar{k}_{\by_0 \rightarrow \by'}  (\by' - \by_0),
\end{equation}

From Definition <ref>, this is equivalent to

\begin{equation}
\bw_0 =  \Bar{\bw}_0.
\end{equation}

Note that if either $\by_0 \not\in V_s$ or $\by_0 \not\in \bar{V}_s$, then one side of Equation (<ref>) gives an empty sum, i.e., $\bw_0 = \mathbf{0}$ or $\bar{\bw}_0 = \mathbf{0}$.

Two dynamically equivalent mass-action systems are presented in Figure <ref>. The mass-action systems (a) $(G, \bk)$ and (b) $(G', \bk')$ share the vertices

\begin{equation} \notag
\by_1 = \begin{pmatrix} 0 \\ 0 \end{pmatrix}, \quad
\by_2 = \begin{pmatrix} 0 \\ 2 \end{pmatrix}, \quad
\by_3 = \begin{pmatrix} 2 \\ 0 \end{pmatrix}.
\end{equation}

The reaction network $G'$ has an additional vertex

\begin{equation} \notag
\by_4 = \begin{pmatrix}  1 \\ 1 \end{pmatrix}.
\end{equation}

Given the rate constants in Figure <ref>, we note that $\by_1$ is the only source vertex in both $G$ and $G'$.
Thus, it suffices to check whether two systems satisfy Equation (<ref>) on the vertex $\by_1$.

For the system $(G, \bk)$, we have

\begin{equation}
\sum_{\by_1 \rightarrow \by \in E} k_{\by_1 \rightarrow \by} (\by - \by_1) =
k_{12} \begin{pmatrix} 0 \\ 2 \end{pmatrix}
+ k_{13} \begin{pmatrix} 2 \\ 0 \end{pmatrix}
= \begin{pmatrix} 4 \\ 4 \end{pmatrix}.
\end{equation}

For the system $(G', \bk')$, we have

\begin{equation}
\sum_{\by_1 \rightarrow \by' \in E'} k_{\by_1 \rightarrow \by'} (\by' - \by_1) =
k'_{12} \begin{pmatrix} 0 \\ 2 \end{pmatrix}
+ k'_{13} \begin{pmatrix} 2 \\ 0 \end{pmatrix}
+ k'_{14} \begin{pmatrix} 1 \\ 1 \end{pmatrix}
= \begin{pmatrix} 4 \\ 4 \end{pmatrix}.
\end{equation}

This shows two systems have the same net reaction vector corresponding to the source vertex $\by_1$, and are hence dynamically equivalent.

Examples of dynamically equivalent networks (a) and (b).

§ WEAKLY REVERSIBLE DEFICIENCY ONE NETWORKS

Deficiency analysis [9, 10, 7, 21] forms an integral component of reaction network theory. The dynamics generated by reaction networks with low deficiency has been studied extensively using the Deficiency zero and Deficiency one theorems [2, 8, 9, 10]. In particular, properties like the existence of a unique equilibrium within each stoichiometric compatibility class, local asymptotic stability of the equilibrium owing to the existence of a Lyapunov function have been established. In this paper, we focus on weakly reversible deficiency one reaction networks. Such networks are ubiquitous in applications, and some noteworthy examples are listed below.

[Edelstein network,

\[
X_1 \rightleftharpoons 2X_1,  \ X_1 + X_2 \rightleftharpoons X_3 \rightleftharpoons X_2
\]

This is a weakly reversible reaction network with deficiency $\delta = 5 - 2- 2 =1$.

[Symmetry breaking network, [2]]

\[
L+2R+P \rightleftharpoons R+Q, \ R+2L+P  \rightleftharpoons 3L+Q, \ P \rightleftharpoons 0  \rightleftharpoons Q
\]

This is a weakly reversible reaction network with the stoichiometric subspace given by:

\begin{equation} \notag
S = \text{span} \{(-1,-1,1,-1)^{\intercal}, \ (1,-1,1,-1)^{\intercal}, \ (0,-1,0,0)^{\intercal}, \ (0,0,-1,0)^{\intercal} \},
\end{equation}

It is a three-dimensional stoichiometric subspace. The network has deficiency $\delta = 7 - 3 - 3 =1$.

From inequality in Equation (<ref>), deficiency one networks can be classified into the following types [Without loss of generality, we always assume $\delta_{1} = \ldots = \delta_{\ell-1} =0, \ \delta_{\ell} = 1$ in Type I networks, and $\delta_{1} = \ldots = \delta_{\ell} =0$ in Type II networks in the rest of this paper.]:

* $\delta = 1 = \delta_1 + \delta_2 + \cdots + \delta_{\ell}$.

We call this a Type I network.

* $\delta = 1 > \delta_1 + \delta_2 + \cdots + \delta_{\ell}$.

We call this a Type II network.

Weakly reversible deficiency one networks for which $\delta = 1 = \delta_1 + \delta_2 + \cdots + \delta_{\ell}$ (Type I)
fall into the regime of the Deficiency one Theorem, which we state below.

Consider a reaction network $G$ consisting of $\ell$ linkage classes $L_1,L_2,\cdots, L_{\ell}$. Let us assume that $G$ satisfies the following conditions:

* $\delta_i\leq 1$.

* $\sum\limits_{i=1}^{\ell} \delta_i = \delta$.

* Each linkage class $L_i$ contains exactly one terminal strongly connected component.

If there exists a $\bk$ for which the mass-action system $(G,\bk)$ possesses a positive equilibrium, then every stoichiometric compatibility class has exactly one positive equilibrium. If $G$ is weakly reversible, then for all values of $\bk$ the mass-action system $(G,\bk)$ possesses a positive equilibrium.

We also state a theorem [11] that guarantees the existence of positive steady states for weakly reversible  systems.

For weakly reversible mass-action systems, there exists a positive steady state within each stoichiometric compatibility class.

Using Theorem <ref> in conjunction with the Deficiency one Theorem, we conclude that for any weakly reversible deficiency one network of Type I, there exists a unique equilibrium within each stoichiometric compatibility class for all values of the rate constants $\bk$.

For weakly reversible deficiency one networks of Type II, all linkage classes have deficiency zero and they possess the following geometric property:

Consider a reaction network $G$. Let $L_1$ be a linkage class of $G$. Then $L_1$ has deficiency zero if and only if its vertices are affinely independent.

Recall that a set $X$ is a polyhedral cone if $X = \{\bx: M \bx \leq \textbf{0} \text{ for some matrix } M \}$. Such a cone is convex. It is pointed, or strongly convex if it does not contain a positive dimensional linear subspace. A pointed polyhedral cone admits a unique (up to scalar multiple) minimal set of generators, and these generating vectors are called extreme vectors [24].

Consider a mass-action system $(G,\bk)$ with vertices $\{\by_i\}_{i=1}^m$. Let $\bW$ be the matrix of net reaction vectors of $G$, then we have:

* $\text{ker}(\bW) \cap \RR^m_{\geq 0}$ is a pointed polyhedral cone.

* There exists the minimal set of generators for $\text{ker}(\bW) \cap \RR^m_{\geq 0}$.

* It is clear that the set $\text{ker}(\bW) \cap \RR^m_{\geq 0}$ is the solution to $\bW \nu \geq \textbf{0}, -\bW \nu \geq \textbf{0}$, and $\bI_{m} \nu\geq\textbf{0}$, and the set is a polyhedral cone. From the definition, a cone contained in the positive orthant $\RR^m_{\geq 0}$ is always pointed.
Therefore, we deduce that $\text{ker}(\bW) \cap \RR^m_{\geq 0}$ is a pointed polyhedral cone.

* Since $\text{ker}(\bW)\cap\mathbb{R}^m_{\geq 0}$ is a pointed cone, by the Minkowski-Weyl theorem [24, 25], we have

\begin{equation}
\text{ker}(\bW)\cap\mathbb{R}^m_{\geq 0}
= \sum\limits_{j=1}^r \zeta_i \bd_i,
\end{equation}

where $\{\bd_j\}_{j=1}^r $ is the unique (up to scalar multiple) minimal set of generators of the cone $\text{ker}(\bW)\cap\mathbb{R}^m_{\geq 0}$.

Furthermore, using the rank-nullity theorem, we have

\begin{equation} \notag
\text{dim}(\text{ker}(\bW)) = \text{dim}(\text{ker}(\bA_{\bk})) + \text{dim}(\text{ker}(\bY)\cap\text{Im}(\bA_{\bk})).
\end{equation}

Thus the minimal set of generators can be divided into two groups. The next definition illustrates this point.

Consider a mass-action system $(G,\bk)$ and let $\bW$ be the matrix of net reaction vectors of $G$. An extreme vector $\bd_i$ of the cone $\text{ker}(\bW)\cap\mathbb{R}^m_{\geq 0}$ is called

* a cyclic generator, if $\bd_i\in\text{ker}(\bA_{\bk})$.

* a stoichiometric generator, if $\bA_{\bk} \bd_i \in \text{ker}(\bY) \backslash
\{\mathbf{0}\}$.

Here we give an example where a reaction network possesses both cyclic and stoichiometric generators.

Consider the network shown in Figure <ref>. This weakly reversible reaction network has two linkage classes, and the deficiency of the entire network is one (i.e. $\delta = 1$). Moreover, the net reaction vector matrix follows:

\begin{equation}
\bW = \begin{pmatrix}
1 & 1   & -2	 & 1 & -1\\
1  & -1 & 0   & 0  & 0
\end{pmatrix},
\end{equation}

\begin{equation}
\text{ker}(\bW)) = \text{span} \big\{ (1,1,1,0,0)^{\intercal}, \ (1,1,0,-2,0)^{\intercal}, \
(1,1,0,0,2)^{\intercal} \big\}.
% \left\{\begin{pmatrix}
% 1 \\
% 1 \\
% 1 \\
% 0 \\
% 0
% \end{pmatrix}, \ \
% \begin{pmatrix}
% -1 \\
% -1 \\
% 0 \\
% 2\\
% 0
% \end{pmatrix}, \ \
% \begin{pmatrix}
% 1 \\
% 1 \\
% 0 \\
% 0 \\
% 2
% \end{pmatrix}\right\}.
\end{equation}

Therefore, we can compute the minimal set of generators of $\text{ker}(\bW)\cap\mathbb{R}^5_{\geq 0}$:

(i) Cyclic generators:

\begin{equation}
\bd_1 = (1,1,1,0,0)^{\intercal}, \ \
\bd_2 = (0,0,0,1,1)^{\intercal}.
\end{equation}

(ii) Stoichiometric generators:

\begin{equation}
\bd_3 = (1,1,0,0,2)^{\intercal}, \ \
\bd_4 = (0,0,1,2,0)^{\intercal}.
\end{equation}

The mass-action system corresponds to Example <ref>, which has both cyclic and stoichiometric generators.

In general, both cyclic and stoichiometric generators can be studied by flux mode analysis.
Further, Conradi et al. [26] defined subnetworks generated by stoichiometric generators, and showed that under some  conditions if these subnetworks exhibit multistationarity, then so does the original network.

As remarked before, weakly reversible deficiency one realizations of Type I satisfy the conditions of the Deficiency One Theorem. This implies that there exists a unique equilibrium within each stoichiometric compatibility class for all values of the rate constants $\bk$ of these realizations. Weakly reversible deficiency one realizations of Type II are also important since the subnetworks generated by the stoichiometric generators can help answer questions about multistationarity. It is therefore important to identify and analyze weakly reversible deficiency one realizations.

§ THE POINTED CONE KERWRM>0

The goal of this section is to analyze the pointed cone $\text{ker}(\bW)\cap\mathbb{R}^m_{\geq 0}$ for weakly reversible deficiency one reaction networks. Specifically, we focus on the extreme vectors of $\text{ker}(\bW)\cap\mathbb{R}^m_{\geq 0}$.

Consider a weakly reversible mass-action system $(G,\bk)$ with vertices $\{\by_i\}_{i=1}^m$.

Let $\bW$ be the matrix of net reaction vectors of $G$, and $\{\bd_1, \ldots, \bd_r\}$ be the minimal set of generators of $\text{ker}(\bW) \cap \mathbb{R}^m_{\geq 0}$, then

\begin{equation} \label{lem:supp_whole}
\bigcup\limits_{i=1}^r \text{supp}(\bd_i) = [m].
\end{equation}

For contradiction, assume there exists $j \in [m]$, such that $j \notin \bigcup\limits_{i=1}^r \text{supp}(\bd_i)$.  Then for any $\bv = (v_1, \ldots, v_m)^{\intercal} \in\text{ker}(\bW)\cap\mathbb{R}^m_{\geq 0}$, we obtain that

\begin{equation} \label{vj = 0}
\bv_j=0.
\end{equation}

Since $(G,\bk)$ is a weakly reversible mass-action system, by Proposition <ref> there exists a positive vector in the kernel of the Kirchoff matrix $\bA_{\bk}$.
Note that weakly reversibility indicates $\bW = \bY\bA_{\bk}$. Thus we have

\begin{equation} \notag
\text{ker}(\bA_{\bk}) \subseteq \text{ker}(\bY\bA_{\bk}) = \text{ker}(\bW).
\end{equation}

This implies the existence of a positive vector in $\text{ker}(\bW)\cap\mathbb{R}^m_{\geq 0}$, contradicting Equation (<ref>).

Consider a weakly reversible mass-action system $(G,\bk)$ with vertices $\{\by_i\}_{i=1}^m$ and stoichiometric subspace $S$. Let

$\bW$ be the matrix of net reaction vectors of $G$, then

\begin{equation} \label{lem:image_stoich}
\text{Im}(\bW) = S.
\end{equation}

The following lemma concerns the dimension of  $\text{ker}(\bW)$

in various cases.

Consider a weakly reversible mass-action system $(G,\bk)$ with vertices $\{\by_i\}_{i=1}^m$ and stoichiometric subspace $S$. Let $\bW$ be the matrix of net reaction vectors of $G$.

*  If $G$ has deficiency $\delta$ and a single linkage class (i.e. $\ell = 1$), we have

\begin{equation}
\text{dim}(\text{ker}(\bW)) = \delta+ 1.
\end{equation}

Moreover, if $\delta = 0$, then for any $\bz\in\text{ker}(\bW) \backslash \{\mathbf{0}\}$, $\text{supp}(\bz)= [m]$.

*  If $G$ has deficiency one and $\ell \geq 1$ linkage classes, we have

\begin{equation}
\text{dim}(\text{ker}(\bW)) = \ell + 1.
\end{equation}

Since $G$ has deficiency $\delta$ and one linkage class, we have

\begin{equation} \notag
\dim (S) = s = m- (\delta + 1).
\end{equation}

By Lemma <ref>, $\text{rank}(\bW) = \dim(\text{Im}(\bW)) = s$. Using the rank-nullity theorem, we obtain

\begin{equation} \notag
\text{dim}(\text{ker}(\bW)) = \delta+ 1.
\end{equation}

Furthermore, if $\delta = 0$, we deduce that

\begin{equation} \label{dim ker_w delta 0}
\text{dim}(\text{ker}(\bW)) = 1.
\end{equation}

As a weakly reversible mass-action system, $(G,\bk)$ possesses a strictly positive steady state $\hat{\bx} \in \mathbb{R}^n_{>0}$ by Theorem <ref>.
Using Equation (<ref>), we get

\begin{eqnarray} \label{eq:ker_w delta 0}
\sum\limits_{j\in [m]}{\hat{\bx}}^{\by_j}\sum\limits_{\by_j\to \by_j' \in \mathcal{R}} k_{\by_j\to \by_j'}(\by_j'- \by_j) = \sum\limits_{j\in [m]}{\hat{\bx}}^{\by_j} \bw_j = 0.
\end{eqnarray}

Note that $({\hat{\bx}}^{\by_1}, {\hat{\bx}}^{\by_2},\ldots, {\hat{\bx}}^{\by_m}) \in \mathbb{R}^m_{>0}$, and it spans $\text{ker}(\bW)$ due to Equation (<ref>).

Since the deficiency of $G$ is one, we get

\begin{equation} \notag
\dim (S) = s = m - \ell - \delta = m - (\ell + 1).
\end{equation}

From Lemma <ref>, we conclude that

\begin{equation} \notag
\text{dim}(\text{ker}(\bW)) = m - \dim(\text{Im}(\bW)) = m - \dim (S) = \ell + 1.
\end{equation}

Here we start with the minimal set of generators of $\text{ker}(\bW) \cap \mathbb{R}^m_{\geq 0}$, when weakly reversible mass-action systems contain a single linkage class.

Consider a weakly reversible mass-action system $(G,\bk)$ that has deficiency $\delta$ and a single linkage class $L = \{ \by_1, \ldots, \by_m \}$.
Let $\bW$ be the matrix of net reaction vectors of $G$, and $\{\bd_1, \ldots, \bd_r\}$ be the minimal set of generators of $\text{ker}(\bW) \cap \mathbb{R}^m_{\geq 0}$, then

\begin{equation} \label{eq: generator on deficiency}
r \geq \delta + 1.
\end{equation}

Moreover, if $\delta = 1$, then $r = 2$.
Assume $\{\bd_1,\bd_2\}$ is the minimal set of generators, then

\begin{equation}
\label{support of deficiency_one_cases part1}
\text{supp}(\bd_1) \subsetneq [m], \ \
\text{supp}(\bd_2) \subsetneq [m], \ \
\text{supp}(\bd_1) \cup \text{supp}(\bd_2) = [m].
\end{equation}

Since $G$ has deficiency $\delta$ and one linkage class, from Lemma <ref>.<ref>
we obtain

\begin{equation} \notag
\text{dim}(\text{ker}(\bW)) = \delta +1.
\end{equation}

Using Equation (<ref>) in Lemma <ref>, we set $\bd = (\bx^{\by_1}, \bx^{\by_2},\ldots, \bx^{\by_m})$ where $\bx \in \mathbb{R}^n_{>0}$ is a steady state for the system, and obtain $\bd \in \text{ker}(\bW) \cap \mathbb{R}^m_{\geq 0}$.
Then there exists a basis of $\text{ker} (\bW)$ that contains $\bd$ as follows.

\begin{equation} \notag
B = \{ \bd, \be_1, \ldots, \be_{\delta} \}.
\end{equation}

Since $\bd \in \mathbb{R}^m_{>0}$, for any weights
$\lambda_1, \ldots, \lambda_{\delta}$, we can always find a sufficiently large $\lambda > 0$, such that

\begin{equation} \notag
\sum\limits^{\delta}_{i=1} \lambda_i \be_i + \lambda \bd \in  \mathbb{R}^m_{>0}.
\end{equation}

Thus, we conclude

\begin{equation} \notag
r \geq \text{dim}(\text{ker}(\bW))= \delta +1.
\end{equation}

Furthermore, if the system has deficiency one (i.e. $\delta = 1$), we derive that

\begin{equation} \notag
\text{dim}(\text{ker}(\bW)) = \delta +1 = 2,
\end{equation}

and thus $\text{ker} (\bW) \cap \mathbb{R}^m_{\geq 0}$ is a two-dimensional pointed cone.
Therefore, the cone $\text{ker} (\bW) \cap \mathbb{R}^m_{\geq 0}$ must have two generators, i.e., $r=2$.

Now assume $\{\bd_1,\bd_2\}$ is the minimal set of generators of $\text{ker} (\bW) \cap \mathbb{R}^m_{\geq 0}$ when the system has deficiency one. Using $\bd = (\bx^{\by_1}, \bx^{\by_2},\ldots, \bx^{\by_m}) \in \mathbb{R}^m_{>0}$, we derive that

\begin{equation} \notag
\text{supp}(\bd_1) \cup \text{supp}(\bd_2) = [m].
\end{equation}

Suppose $\text{supp}(\bd_1) = [m]$, thus $\bd_1 \in \mathbb{R}^m_{>0}$. Then we can find a sufficiently large $\lambda > 0$, such that

\begin{equation} \notag
\bv = \lambda \bd_1 - \bd_2 \in  \text{ker} (\bW) \cap \mathbb{R}^m_{>0}.
\end{equation}

Note that $\bd_1$ and $\bd_2$ are linearly independent, this contradicts with $\{\bd_1,\bd_2\}$ being the generating set of $\text{ker} (\bW) \cap \mathbb{R}^m_{\geq 0}$. Thus, we derive that $\text{supp}(\bd_1) \subsetneq [m]$. Similarly, we can show $\text{supp}(\bd_2) \subsetneq [m]$, and conclude (<ref>).

Consider a weakly reversible mass-action system $(G,\bk)$ with a single linkage class $L = \{ \by_1, \ldots, \by_m \}$, and let $\bW$ be the matrix of net reaction vectors of $G$. Then there exists a vector $\bd \in \mathbb{R}^m_{\geq 0}$ generating the cone $\text{ker}(\bW) \cap \mathbb{R}^m_{\geq 0}$ if and only if the system has deficiency zero.

First, suppose the system has deficiency zero.
From Equation (<ref>) in Lemma <ref>, we set $\bd = (\bx^{\by_1}, \bx^{\by_2},\ldots, \bx^{\by_m})$ where $\bx \in \mathbb{R}^n_{>0}$ is a steady state for the system $(G,\bk)$, and obtain

\begin{equation} \notag
\text{ker}(\bW) = \text{span}\{ \bd \}.
\end{equation}

One can check $\bd \in \mathbb{R}^m_{>0}$, and hence $\bd$ generates $\text{ker}(\bW) \cap \mathbb{R}^m_{\geq 0}$.

On the other hand, consider a vector $\bd$ that generates $\text{ker}(\bW) \cap \mathbb{R}^m_{\geq 0}$. Using (<ref>) and deficiency is non-negative, we get that

\begin{equation}
0 \leq \delta \leq 1 - 1 = 0.
\end{equation}

Thus, we conclude the deficiency of the system is zero.

The remark below follows from Lemmas <ref> and <ref>.

Consider a weakly reversible mass-action system $(G,\bk)$ with a single linkage class $L = \{\by_i\}_{i=1}^m$.
Let $\bW$ be the matrix of net reaction vectors of $G$.
Suppose two vectors $\bd_1, \bd_2$ form the minimal set of generators of $\text{ker}(\bW) \cap \mathbb{R}^m_{\geq 0}$, then $(G,\bk)$ has deficiency one.

Next, we work on the minimal set of generators of $\text{ker}(\bW)\cap\mathbb{R}^m_{\geq 0}$ when the weakly reversible deficiency one networks have multiple linkage classes.

Consider a weakly reversible deficiency one mass-action system $(G,\bk)$ of Type I that has $\ell > 1$ linkage classes, denoted by $L_1, \ldots, L_{\ell}$. Let $\bW$ be the matrix of net reaction vectors of $G$, and $\{\bW_p \}^{\ell}_{p=1}$ be the matrix of net vectors corresponding to linkage classes $\{L_p \}^{\ell}_{p=1}$, then

\begin{equation} \label{dim of deficiency_one_cases part2}
\text{dim}(\text{ker}(\bW_1)) + \cdots + \text{dim}(\text{ker}(\bW_{\ell}))  =  \text{dim}(\text{ker}(\bW)) = \ell + 1,
\end{equation}

\begin{equation} \notag
\text{dim}(\text{ker}(\bW_i)) =
\begin{cases}
1, & \text{for } 1 \leq i \leq \ell-1, \\[5pt]
2, & \text{for } i = \ell.
\end{cases}
\end{equation}

Moreover, for any $1 \leq i \leq \ell-1$ and $\bz \in \text{ker}(\bW_i) \backslash \{\mathbf{0}\}$, we have $\text{supp} (\bz) = L_i$.

* There exist $\ell + 1$ vectors $\bd_1, \ldots, \bd_{\ell + 1}$, which form the minimal set of generators of the cone $\text{ker} (\bW)\cap\mathbb{R}^m_{\geq 0}$, such that

\begin{equation}
\label{support of deficiency_one_cases part2 a}
\text{supp}(\bd_i) = L_i, \ \text{for } 1 \leq i \leq \ell-1,
\end{equation}

\begin{equation}
\label{support of deficiency_one_cases part2 b}
\text{supp}(\bd_{\ell}) \subsetneq L_{\ell}, \ \
\text{supp}(\bd_{\ell+1}) \subsetneq L_{\ell}, \ \
\text{supp}(\bd_{\ell})\cup \text{supp}(\bd_{\ell +1}) = L_{\ell}.
\end{equation}

$(a)$ From the assumption, the system $G$ is of Type I with $\delta_{1} = \cdots = \delta_{\ell-1} =0$ and $\delta_{\ell} = 1$.
Using Lemma <ref>, we get

\begin{equation} \notag
\text{dim}(\text{ker}(\bW_{\ell})) = \delta_{\ell} + 1 = 2,
\ \
\text{dim}(\text{ker}(\bW_i)) = \delta_{i} + 1 = 1, \ \text{for } 1 \leq i \leq \ell-1.
\end{equation}

Further, for any $1 \leq i \leq \ell-1$ and $\bz \in \text{ker}(\bW_i) \backslash \{\mathbf{0}\}$,

\begin{equation} \notag
\text{supp}(\bz) = L_i.
\end{equation}

Note that $G$ has deficiency one, thus
$\text{dim}(\text{ker}(\bW)) = \ell + 1$,
and we derive (<ref>).

Now we construct the minimal set of generators of $\text{ker}(\bW) \cap \mathbb{R}^m_{\geq 0}$, denoted by $\{\bd_1, \ldots,\bd_{r}\}$. It will follow from the construction that $r = \ell + 1=  \text{dim}(\text{ker}(\bW))$.

Since $(G,\bk)$ is a weakly reversible mass-action system, it possesses a strictly positive steady state $\hat{\bx}\in \mathbb{R}^n_{>0}$ by Theorem <ref>. Following Equation (<ref>) in Lemma <ref>, we can build $\ell-1$ vectors $\bd_1, \ldots, \bd_{\ell - 1}$. We define $\bd_1 = (\bd_{1,1},\ldots,\bd_{1,m}) \in \text{ker}(\bW) \cap \mathbb{R}^m_{\geq 0}$, such that

\begin{eqnarray} \label{eq:defn_generators_1 a}
\bd_{1,i} =
\begin{cases}
\hat{\bx}^{\by_i}, & \text{for } i \in L_1, \\
0, & \text{for } i \notin L_1.
\end{cases}
\end{eqnarray}

It is clear that $\text{supp}(\bd_1) = L_1$. Analogously, for $i=1, \ldots, \ell - 1$, we can define $\bd_i$ corresponding to the linkage classes $L_i$ with $\text{supp}(\bd_i) = L_i$.

Note that $G$ is of Type I and the linkage class $L_{\ell}$ has deficiency one.
Let $L_{\ell} = \big\{ \by_{\ell_i} \big\}^{m_{\ell}}_{i = 1}$ with $m_{\ell} = |L_{\ell}|$. From Lemma <ref>, the cone $\text{dim}(\text{ker}(\bW_{\ell}))\cap\mathbb{R}^{m_{\ell}}_{\geq 0}$ has two generators, denoted by $\{ \bg_1, \bg_2 \}$. Suppose $\bg_1 = \big( \bg_{1,i} \big)_{i \in L_{\ell}}$ and   $\bg_2 = \big( \bg_{2,i} \big)_{i \in L_{\ell}}$, then we define $\bd_{\ell} = (\bd_{\ell, 1}, \ldots,\bd_{\ell, m})$, $\bd_{\ell+1} = (\bd_{\ell+1,1}, \ldots, \bd_{\ell+1,m})$ as

\begin{equation}
\label{eq:defn_generators_1 b}
\begin{split}
\bd_{\ell,i}
= \begin{cases}
\bg_{1,i}, & \text{for } i \in L_{\ell}, \\
0, & \text{for } i \notin L_{\ell}.
\end{cases}
\ \ \text{and } \
\bd_{\ell + 1,i}
= \begin{cases}
\bg_{2,i}, & \text{for } i \in L_{\ell}, \\
0, & \text{for } i \notin L_{\ell}.
\end{cases}
\end{split}
\end{equation}

Note that both $\bd_{\ell}, \bd_{\ell+1} \in \text{ker}(\bW) \cap \mathbb{R}^m_{\geq 0}$, and satisfy Equation (<ref>).

We claim that the vectors $\bd_1, \ldots , \bd_{\ell}$ form a set of generators for $\text{ker}(\bW) \cap \mathbb{R}^m_{\geq 0}$.
From Equations (<ref>) and (<ref>), we deduce that the vectors $\{ \bd_i \}^{\ell + 1}_{i=1}$ are linearly independent. Together with $\text{dim}(\text{ker}(\bW))=\ell + 1$, we derive that the set $\{\bd_1, \ldots ,\bd_{\ell + 1}\}$ is a basis for $\text{ker}(\bW)$. Thus, any vector $\bv \in \text{ker}(\bW) \cap \mathbb{R}^m_{\geq 0}$ can be expressed as

\begin{equation} \label{linear combination 1}
\bv = a_1 \bd_1 + a_2 \bd_2 + \cdots + a_{\ell+1} \bd_{\ell+1} \in \mathbb{R}^m_{\geq 0},
\end{equation}

where $a_1, \ldots, a_{\ell+1} \in \mathbb{R}$.
So it suffices to prove all $\{ a_i \}^{\ell+1}_{i=1}$ are non-negative.
Recall $\{ L_i\}^{\ell}_{i=1}$ are linkage classes with $\text{supp}(\bd_i) = L_i$, and $\text{supp}(\bd_{\ell})$, $\text{supp}(\bd_{\ell+1}) \subseteq L_{\ell}$, then we obtain

\begin{equation} \notag
a_i \geq 0, \ \text{for } i=1, \ldots, \ell-1.
\end{equation}

Moreover, we set $\hat{\bv} = a_{\ell} \bg_1 + a_{\ell+1} \bg_2$. From $\sum\limits^{\ell}_{i=1} \text{dim}(\text{ker}(\bW_i)) = \text{dim}(\text{ker}(\bW))$, we derive

\begin{equation} \notag
\hat{\bv} \in \text{ker}(\bW_{\ell})\cap\mathbb{R}^{m_{\ell}}_{\geq 0}.
\end{equation}

Since $\bg_1, \bg_2$ form the generators of the cone $\text{dim}(\text{ker}(\bW_{\ell}))\cap\mathbb{R}^{m_{\ell}}_{\geq 0}$, we have

\begin{equation} \notag
a_{\ell} \geq 0, \ \
a_{\ell+1} \geq 0.
\end{equation}

Therefore, we prove the claim.

Finally, we show $\{\bd_1, \bd_2, \ldots, \bd_{\ell +1}\}$ is the minimal set of generators for $\text{ker}(\bW) \cap \mathbb{R}^m_{\geq 0}$.
Note from Equations (<ref>) and (<ref>), $\bd_{\ell + 1}$ cannot be generated by $\{\bd_i\}^{\ell}_{i=1}$, thus it suffices to show $\{\bd_1,\bd_2, \ldots ,\bd_{\ell}\}$ are all extreme vectors.

Suppose not, there exists $1 \leq j \leq \ell$, such that $\bd_j$ is not an extreme vector.
Then we can find two vectors $\gamma, \theta \in \text{ker}(\bW) \cap \mathbb{R}^m_{\geq 0}$ and $0 < \lambda < 1$, such that

\begin{equation} \label{d_decompose_1}
\lambda \gamma + (1 - \lambda) \theta = \bd_j,
\end{equation}

where $\gamma \neq \nu \theta$ for any constant $\nu$.
From Equation (<ref>), we can express $\gamma$ and $\theta$ as the conical combination of $\{\bd_1,\bd_2, \ldots ,\bd_{\ell + 1}\}$ as

\begin{equation} \notag
\gamma = \sum\limits^{\ell+1}_{i=1} \gamma_{i}\bd_i
\ \ \text{and } \
\theta = \sum\limits^{\ell+1}_{i=1} \theta_{i} \bd_i,
\end{equation}

where  $\gamma_{i}, \ \theta_{i} \geq 0$, for $i = 1, \ldots, \ell+1$.

If $j \neq \ell$, from $\text{supp}(\bd_j) = L_j$ and Equation (<ref>), we derive that $\gamma_{i} = \theta_{i} = 0$ for $1 \leq i \leq \ell+1, i \neq j$.
This implies $\gamma = \gamma_j d_j$ and $\theta = \theta_j d_j$, which contradicts with $\gamma \neq \nu \theta$.

If $j = \ell$, we deduce that $\gamma_{i} = \theta_{i} = 0$ for $1 \leq i \leq \ell+1$ such that $i \neq j$ in a similar way. This implies $\gamma = \gamma_{\ell}\bd_{\ell}$ and $\theta = \theta_{\ell}\bd_{\ell}$, which also contradicts with $\gamma \neq \nu \theta$.
Therefore we conclude that $\{\bd_1, \ldots, \bd_{\ell + 1}\}$ is the minimal set of generators of the cone $\text{ker}(\bW) \cap \mathbb{R}^m_{\geq 0}$.

Here we illustrate an example where Lemma <ref> can be verified.

Consider a weakly reversible deficiency one mass-action system

shown in Figure <ref>.
This reaction network has two linkage classes.  One linkage class has deficiency zero, and the other has deficiency one (i.e. $\delta_1 = 0, \ \delta_2 = 1$), and the deficiency of the entire network is one (i.e. $\delta = 1$). Therefore, we have

\begin{equation}
1 = \delta = \delta_1 + \delta_2.
\end{equation}

For all reactions $\by\rightarrow\by' \in E$, we assume
$k_{\by\rightarrow\by'} = 1$, and get

\begin{equation}
\bW_1 =
\begin{pmatrix}
1 & -1\\
0  & 0
\end{pmatrix}, \ \
\bW_2 =
\begin{pmatrix}
0 &  0 & 0\\
1  &  0  &-1
\end{pmatrix}, \ \
\bW =
\begin{pmatrix}
1 & -1 & 0 &  0 & 0 \\
0  & 0 & 1  &  0  & -1
\end{pmatrix}.
\end{equation}

So we can derive that

\begin{equation}
\text{ker}(\bW_1)) = \text{span}
\left\{
\begin{pmatrix}
1 \\
1 \\
\end{pmatrix}
\right\}, \
\text{ker}(\bW_2)) = \text{span}
\left\{ \begin{pmatrix}
0 \\
1 \\
0 \\
\end{pmatrix}, \ \
\begin{pmatrix}
1 \\
0 \\
1 \\
\end{pmatrix}\right\},
\end{equation}

\begin{equation}
\text{ker}(\bW)) = \text{span}
\left\{\begin{pmatrix}
1 \\
1 \\
0 \\
0 \\
0 \\
\end{pmatrix}, \ \
\begin{pmatrix}
0 \\
0 \\
0 \\
1 \\
0 \\
\end{pmatrix}, \ \
\begin{pmatrix}
0 \\
0 \\
1 \\
0 \\
1 \\
\end{pmatrix}\right\}.
\end{equation}

For any vector $\bz_1 \in \text{ker}(\bW_1) \backslash \{\mathbf{0}\}$, we have

\begin{equation}
\text{supp}(\bz_1)= L_1.
% \ \
% \text{dim}(\text{ker}(\bW_1)) + \text{dim}(\text{ker}(\bW_2)) = \text{dim}(\text{ker}(\bW)) = 3.
\end{equation}

Then, we compute the minimal set of generators of $\text{ker}(\bW)\cap\mathbb{R}^5_{\geq 0}$:

\begin{equation}
\bd_2 =
\begin{pmatrix}
1 \\
1 \\
0 \\
0 \\
0 \\
\end{pmatrix}, \ \
\bd_2 =
\begin{pmatrix}
0 \\
0 \\
0 \\
1 \\
0 \\
\end{pmatrix}, \ \
\bd_3 = \begin{pmatrix}
0 \\
0 \\
1 \\
0 \\
1 \\
\end{pmatrix}
\end{equation}

This shows that the number of extreme vectors: $r = 3$, and

\begin{equation}
r = \text{dim}(\text{ker}(\bW)) = \ell +1,
\end{equation}

where $\ell = 2$.
Moreover, for $q = 2, 3$,

\begin{equation}
\text{supp}(\bd_1) = L_1, \ \
\text{supp}(\bd_{q}) \cap L_2 \subsetneq L_2, \ \
\text{supp}(\bd_{2})\cup \text{supp}(\bd_{3}) = L_{2}.
\end{equation}

Therefore, we verify Lemma <ref>.

A weakly reversible deficiency one mass-action system of Type I from Example <ref>

Consider a weakly reversible deficiency one mass-action system $(G,\bk)$ of Type II that has $\ell > 1$ linkage classes denoted by $L_1, \ldots, L_{\ell}$. Let $\bW$ be the matrix of net reaction vectors of $G$, and $\{\bW_p \}^{\ell}_{p=1}$ be the matrix of net reaction vectors corresponding to linkage classes $\{L_p \}^{\ell}_{p=1}$, then

\begin{equation} \label{dim of deficiency_one_cases part3}
\text{dim}(\text{ker}(\bW_1)) + \cdots + \text{dim}(\text{ker}(\bW_{\ell}))  =  \text{dim}(\text{ker}(\bW)) - 1 = \ell,
\end{equation}

\begin{equation} \notag
\text{dim}(\text{ker}(\bW_i)) = 1, \ \text{for } 1 \leq i \leq \ell.
\end{equation}

Moreover, for any $1 \leq i \leq \ell$ and $\bz \in \text{ker}(\bW_i) \backslash \{\mathbf{0}\}$, we have $\text{supp} (\bz) = L_i$.

(b) There exist $\ell + 2$ vectors $\bd_1, \ldots, \bd_{\ell + 2}$, which form the minimal set of generators of the cone $\text{ker} (\bW)\cap\mathbb{R}^m_{\geq 0}$, such that for $i = 1, \ldots, \ell$,

\begin{equation}
\label{support of deficiency_one_cases part 3}
\text{supp}(\bd_i) = L_i, \ \
\emptyset \neq \text{supp}(\bd_{\ell+1}) \cap L_i \subsetneq L_{i}, \ \
\emptyset \neq \text{supp}(\bd_{\ell+2}) \cap L_i \subsetneq L_{i}.
\end{equation}

From the assumption, the system $G$ is of Type II with $\delta_{1} = \cdots = \delta_{\ell} =0$.
Using Lemma <ref>, we get for $i=1, \ldots, \ell$,

\begin{equation} \notag
\text{dim}(\text{ker}(\bW_i)) = \delta_{i} + 1 = 1.
\end{equation}

Further, for any $1 \leq i \leq \ell$ and $\bz \in \text{ker}(\bW_i) \backslash \{\mathbf{0}\}$,

\begin{equation} \notag
\text{supp}(\bz) = L_i.
\end{equation}

Note that $G$ has deficiency one, thus
$\text{dim}(\text{ker}(\bW)) = \ell + 1$,
and we derive (<ref>).

$(b)$ Now we construct the minimal set of generators of $\text{ker}(\bW) \cap \mathbb{R}^m_{\geq 0}$, denoted by $\{\bd_1, \ldots,\bd_{r}\}$. It will follow from the construction that $r = \ell + 2 = \text{dim}(\text{ker}(\bW)) + 1$.

Since $(G,\bk)$ is a weakly reversible mass-action system, it possesses a strictly positive steady state $\hat{\bx}\in \mathbb{R}^n_{>0}$ by Theorem <ref>. Following Equation (<ref>) in Lemma <ref>, we can build $\ell$ vectors $\bd_1, \ldots, \bd_{\ell}$. We define $\bd_1 = (\bd_{1,1},\ldots,\bd_{1,m}) \in \text{ker}(\bW) \cap \mathbb{R}^m_{\geq 0}$, such that

\begin{eqnarray} \label{eq:defn_generators_2 a}
\bd_{1,i} =
\begin{cases}
\bx^{\by_i}, & \text{for } i \in L_1, \\
0, & \text{for } i \notin L_1,
\end{cases}
\end{eqnarray}

It is clear that $\text{supp}(\bd_1) = L_1$. Analogously, for $i=1, \ldots, \ell$, we can define $\bd_i$ corresponding to the linkage classes $L_i$, with $\text{supp}(\bd_i) = L_i$.

Now we show that there exists a non-zero vector $\bd_{\ell + 1} \in \text{ker}(\bW) \cap \mathbb{R}^m_{\geq 0}$, such that

\begin{equation} \label{support l+1}
\text{supp}(\bd_{\ell + 1}) \cap L_i \subsetneq L_i, \ \text{for } i = 1, \ldots, \ell.
\end{equation}

From Equation (<ref>), there exists a vector $\tilde{\bd} \in \text{ker}(\bW) \backslash \{\mathbf{0}\}$, which is linearly independent from $\{\bd_i\}^{\ell}_{i=1}$.
Since $\bd_i \in \text{ker}(\bW) \cap \mathbb{R}^m_{\geq 0}$ with $\text{supp}(\bd_i) = L_i$, we set for $i = 1, \ldots, \ell$,

\begin{equation} \label{max weight l+1}
\alpha_i = \max_{k \in L_i} \Big\{ - \frac{\tilde{\bd}_k}{\bd_{i,k}} \Big\}.
\end{equation}

Then we define $\bd_{\ell + 1}$ as

\begin{equation} \label{eq:defn_generators_2 b}
\bd_{\ell + 1}
= \sum\limits^{\ell}_{i=1} \alpha_i \bd_i + \tilde{\bd}.
\end{equation}

For any $1 \leq j \leq  \ell$ and $\theta \in L_j$, we obtain that

\begin{equation} \notag
\bd_{\ell + 1, \theta} = \alpha_j \bd_{i, \theta} + \tilde{\bd}_{\theta}
\geq  - \frac{\tilde{\bd}_\theta}{\bd_{i,\theta}} \bd_{i, \theta} + \tilde{\bd}_{\theta}
= 0,
\end{equation}

and the inequality holds when $\theta = k \in L_j$ in Equation (<ref>).
Moreover, the linear independence between $\tilde{\bd}$ and $\{\bd_i\}^{\ell}_{i=1}$ implies that $\bd_{\ell + 1}$ is non-zero.
Thus, we show $\bd_{\ell+1} \in \text{ker}(\bW) \cap \mathbb{R}^m_{\geq 0}$, and it satisfies Equation (<ref>).

Furthermore, we claim that there exist at least two linkage classes: $L_i$, $L_j$ with $1 \leq i, j \leq \ell$ and $i \neq j$, such that

\begin{equation} \label{support l+1 on two classes}
\text{supp}(\bd_{\ell + 1}) \cap L_i \neq \emptyset, \ \text{supp}(\bd_{\ell + 1}) \cap L_j \neq \emptyset.
\end{equation}

Suppose not, we assume that only the linkage class $L_1$ satisfies $\text{supp}(\bd_{\ell + 1}) \cap L_1 \neq \emptyset$. This implies that

\begin{equation} \notag
\text{supp}(\bd_{\ell + 1}) \cap L_1 \subsetneq L_1.
\end{equation}

Using $\text{dim}(\text{ker}(\bW_1))=1$, we get that  $\bd_{\ell + 1}$ must be a scalar multiple of $\bd_{1}$, contradicting Equation (<ref>).

Next, we construct another non-zero vector $\bd_{\ell + 2} \in \text{ker}(\bW) \cap \mathbb{R}^m_{\geq 0}$, such that

\begin{equation} \label{support l+2}
\text{supp}(\bd_{\ell + 2}) \cap L_i \subsetneq L_i, \ \text{for} \ i = 1, \ldots, \ell.
\end{equation}

Given $\bd_{1}, \ldots, \bd_{\ell}, \bd_{\ell + 1} \in \mathbb{R}^m_{\geq 0}$, we set for $i = 1, \ldots, \ell$,

\begin{equation} \label{max weight l+2}
\beta_i = \max_{k \in L_i}\Big\{ \frac{\bd_{\ell+1,k}}{\bd_{i,k}} \Big\}.
\end{equation}

It is clear that $\beta_i\geq 0$ for $1\leq i\leq \ell$, then we define $\bd_{\ell + 2}$ as

\begin{equation} \label{eq:defn_generators_2 c}
\bd_{\ell + 2}= \sum\limits^{\ell}_{i=1} \beta_i \bd_i - \bd_{\ell + 1}.
\end{equation}

For any $1 \leq j \leq  \ell$ and $\theta \in L_j$, we get

\begin{equation} \notag
\bd_{\ell + 2, \theta} = \beta_j \bd_{i, \theta} - \bd_{\ell + 1, \theta}
\geq \frac{\bd_{\ell+1,\theta}}{\bd_{i,\theta}} \bd_{i, \theta} - \bd_{\ell + 1, \theta}
= 0
\end{equation}

The inequality holds when $\theta = k \in L_j$ in Equation (<ref>). Moreover, the linear independence between
$\bd_{\ell+1}$ and $\{\bd_i\}^{\ell}_{i=1}$ implies that $\bd_{\ell + 2}$ is non-zero. Thus, we show $\bd_{\ell+2} \in \text{ker}(\bW) \cap \mathbb{R}^m_{\geq 0}$, and it satisfies Equation (<ref>). Similarly as in (<ref>), there also exist at least two linkage classes: $L_i$, $L_j$ with $1 \leq i, j \leq \ell$ and $i \neq j$, such that

\begin{equation} \label{support l+2 on two classes}
\text{supp}(\bd_{\ell + 2}) \cap L_i \neq \emptyset, \ \
\text{supp}(\bd_{\ell + 2}) \cap L_j \neq \emptyset.
\end{equation}

We claim that the vectors $\bd_1, \ldots , \bd_{\ell+2}$ form a set of generators of $\text{ker}(\bW) \cap \mathbb{R}^m_{\geq 0}$. Using Equations (<ref>) and (<ref>), we deduce the vectors $\{ \bd_i \}^{\ell + 1}_{i=1}$ are linearly independent. Together with $\text{dim}(\text{ker}(\bW))=\ell + 1$, we get that the set $\{\bd_1, \ldots ,\bd_{\ell + 1}\}$ is a basis for $\text{ker}(\bW)$. Thus, any vector $\bv \in \text{ker}(\bW) \cap \mathbb{R}^m_{\geq 0}$ can be expressed as

\begin{equation} \label{linear combination 2a}
\bv = a_1 \bd_1 + a_2 \bd_2 + \cdots + a_{\ell+1} \bd_{\ell+1} \in \mathbb{R}^m_{\geq 0},
\end{equation}

where $a_1, \ldots, a_{\ell+1} \in \mathbb{R}$.
Recall $\{ L_i\}^{\ell}_{i=1}$ are linkage classes with $\text{supp}(\bd_i) = L_i$, for $i=1, \ldots, \ell$, and $\text{supp}(\bd_{\ell + 1})$ in Equation (<ref>), then we obtain

\begin{equation} \notag
a_i \geq 0, \ \text{for } i=1, \ldots, \ell.
\end{equation}

If $a_{\ell+1} \geq 0$, it is clear that $\bv$ can be expressed as a conical combination of $\{\bd_1,\bd_2, \ldots ,\bd_{\ell + 1}\}$ from Equation (<ref>).  Otherwise, if $a_{\ell+1} < 0$, we rewrite $\bv$ as

\begin{equation} \label{linear combination 2b}
\begin{split}
\bv
& = a_1 \bd_1 + \cdots + a_{\ell} \bd_{\ell} + a_{\ell+1} ( \sum\limits^{\ell}_{i=1} \beta_i \bd_i - \bd_{\ell +2} )
\\& = a_1\bd_1 + \cdots  + a_{\ell} \bd_{\ell} + a_{\ell+1} \sum\limits^{\ell}_{i=1} \beta_i\bd_i - a_{\ell+1} \bd_{\ell +2}
\\& = (a_1 + a_{\ell+1} \beta_1) \bd_1 + \cdots  + (a_{\ell} + a_{\ell+1} \beta_{\ell}) \bd_{\ell} - a_{\ell+1} \bd_{\ell +2}.
\end{split}
\end{equation}

Using $\bv \in \text{ker}(\bW) \cap \mathbb{R}^m_{\geq 0}$ and Equation (<ref>), we get that for $i = 1, \ldots, \ell$,

\begin{equation} \notag
a_i + a_{\ell+1} \beta_i \geq 0,
\end{equation}

which implies that $\bv$ can be generated by $\{\bd_1, \ldots, \bd_{\ell}, \bd_{\ell+2}\}$.

Finally, we show $\{\bd_1,\bd_2, \ldots ,\bd_{\ell + 2}\}$ is the minimal set of generators for $\text{ker}(\bW) \cap \mathbb{R}^m_{\geq 0}$.  Note that $\bd_1, \ldots ,\bd_{\ell + 1}$ form a basis for $\text{ker}(\bW)$ and $\bd_{\ell + 2} = \sum\limits^{\ell}_{i=1} \beta_i \bd_i - \bd_{\ell + 1}$, thus $\bd_{\ell + 2}$ cannot be generated by $\{\bd_i\}^{\ell+1}_{i=1}$.
So it suffices to show $\{\bd_1,\bd_2, \ldots ,\bd_{\ell + 1}\}$ are all extreme vectors.

Suppose not, there exists $1 \leq j \leq \ell+1$, such that $\bd_j$ is not an extreme vector.
Then we can find two vectors $\gamma, \theta \in \text{ker}(\bW) \cap \mathbb{R}^m_{\geq 0}$ and $0 < \lambda < 1$, such that

\begin{equation} \label{d_decompose_2}
\lambda \gamma + (1 - \lambda) \theta = \bd_j,
\end{equation}

where $\gamma \neq \nu \theta$ for any constant $\nu$. Then we write $\gamma$ and $\theta$ as the combination of
$\{\bd_i\}^{\ell + 1}_{i=1}$,

\begin{equation} \notag
\gamma = \sum\limits^{\ell+1}_{i=1} \gamma_{i}\bd_i, \ \
\theta = \sum\limits^{\ell+1}_{i=1} \theta_{i} \bd_i.
\end{equation}

Since $\gamma, \theta \in \mathbb{R}^m_{\geq 0}$,

we have for $i = 1, \ldots, \ell$,

\begin{equation} \notag
\gamma_{i} \geq 0, \ \theta_{i} \geq 0.
\end{equation}

If $j \neq \ell + 1$, from Equation (<ref>), we can derive that $\gamma_{i} = \theta_{i} = 0$ when $1 \leq i \leq \ell, i \neq j$. Since $\gamma, \theta \in \mathbb{R}^m_{\geq 0}$, and $\text{supp}(\bd_{\ell +1}) \cap L_i \subsetneq L_i$ for $i=1, \ldots, \ell$, we derive $\gamma_{\ell + 1} = \theta_{\ell + 1} = 0$.
This implies $\gamma = \gamma_j d_j$ and $\theta = \theta_j d_j$, and this contradicts with $\gamma \neq \nu \theta$.

If $j = \ell + 1$, in a similar way we can deduce that $\gamma_{i} = \theta_{i} = 0$ for $1 \leq i \leq \ell$. This implies $\gamma = \gamma_{\ell+1}\bd_{\ell+1}$ and $\theta = \theta_{\ell+1}\bd_{\ell+1}$, which also contradicts with $\gamma \neq \nu \theta$.
Therefore we conclude that $\{\bd_1, \ldots, \bd_{\ell + 2}\}$ is the minimal set of generators of the cone $\text{ker}(\bW) \cap \mathbb{R}^m_{\geq 0}$.

We also illustrate an example where Lemma <ref> can be verified.

Consider a weakly reversible deficiency one mass-action system

shown in Figure <ref>.
This reaction network has two deficiency zero linkage classes (i.e. $\delta_1 = \delta_2 = 0$), and the deficiency of the entire network is one (i.e. $\delta = 1$).
Therefore, we have

\begin{equation}
1 = \delta > \delta_1 + \delta_2 = 0 + 0.
\end{equation}

For all reactions $\by\rightarrow\by' \in E$, we assume
$k_{\by\rightarrow\by'} = 1$, and get

\begin{equation}
\bW_1 = \begin{pmatrix}
1 & -1\\
0  & 0
\end{pmatrix}, \ \
\bW_2 = \begin{pmatrix}
1 & -1\\
0  & 0
\end{pmatrix}, \ \
\bW = \begin{pmatrix}
1 & -1 & 1 &  -1 \\
0  & 0 & 0  &  0
\end{pmatrix}.
\end{equation}

So we can derive that

\begin{equation}
\text{ker}(\bW_1)) = \text{ker}(\bW_2)) = \text{span}
\left\{
\begin{pmatrix}
1 \\
1 \\
\end{pmatrix}
\right\},
\end{equation}

\begin{equation}
\text{ker}(\bW)) = \text{span}
\left\{\begin{pmatrix}
1 \\
1 \\
0 \\
0 \\
\end{pmatrix}, \ \
\begin{pmatrix}
-1 \\
0 \\
1 \\
0 \\
\end{pmatrix}, \ \
\begin{pmatrix}
1 \\
0 \\
0 \\
1 \\
\end{pmatrix}\right\}.
\end{equation}

For any vectors $\bz_1 \in \text{ker}(\bW_1) \backslash \{\mathbf{0}\}$ and $\bz_2 \in \text{ker}(\bW_2) \backslash \{\mathbf{0}\}$, we have

\begin{equation}
\text{supp}(\bz_1)= L_1, \ \
\text{supp}(\bz_2)= L_2.
% \ \text{dim}(\text{ker}(\bW_1)) + \text{dim}(\text{ker}(\bW_2)) = \text{dim}(\text{ker}(\bW)) -1 = 2.
\end{equation}

Then, we compute the minimal set of generators of $\text{ker}(\bW)\cap\mathbb{R}^4_{\geq 0}$:

\begin{equation}
\bd_1 =
\begin{pmatrix}
1 \\
1 \\
0 \\
0 \\
\end{pmatrix}, \ \
\bd_2 =
\begin{pmatrix}
0 \\
0 \\
1 \\
\end{pmatrix}, \ \
\bd_3 =
\begin{pmatrix}
0 \\
1 \\
1 \\
0 \\
\end{pmatrix}, \ \
\bd_4 =
\begin{pmatrix}
1 \\
0 \\
0 \\
1 \\
\end{pmatrix}.
\end{equation}

This indicates the number of extreme vectors: $r = 4$ and

\begin{equation}
r = \text{dim}(\text{ker}(\bW)) +1 = \ell +2,
\end{equation}

where $\ell = 2$.
Moreover, for $p=1, 2$,

\begin{equation}
\text{supp}(\bd_p) = L_p, \ \
\text{supp}(\bd_{\ell +1}) \cap L_p \subsetneq L_p, \ \
\text{supp}(\bd_{\ell +2}) \cap L_p \subsetneq L_p.
\end{equation}

Therefore, we verify Lemma <ref>.

A weakly reversible deficiency one mass-action system of Type II from Example <ref>

To conclude this section, we show that given a mass-action system that admits weakly reversible deficiency one realizations, then these realizations must be of the same type.

First, we recall a special result from [28]:

Consider two weakly reversible mass-action systems $(G, \bk)$ and $(G', \bk')$ having a deficiency of one and the same number of linkage classes. Let $G$ be of Type I and $G'$ be of Type II. Then $(G, \bk)$ and $(G', \bk')$ cannot be dynamically equivalent.

After Theorem <ref>, we are ready to prove the  more general result as follows:

Given two weakly reversible deficiency one mass-action systems: $(G, \bk)$ of Type I and $(G', \bk')$ of Type II, then $(G, \bk)$ and $(G', \bk')$ cannot be dynamically equivalent.

For contradiction, assume that the two weakly reversible deficiency one mass-action systems $(G=(V, E), \bk)$ of Type I, and $(G'=(V', E'), \bk')$ of Type II are dynamically equivalent. By Remark <ref>, they have the same set of non-zero net reaction vectors. Using $\text{Im}(\bW) = S$ from Lemma <ref>, we get that $(G, \bk)$ and $(G', \bk')$ share the same stoichiometric subspace.

Now we claim that $(G, \bk)$ and $(G', \bk')$ have the same number of vertices. For contradiction, suppose there exists a vertex $\by \in V'$ such that $\by \notin V$. Let $\bw_{\by}$ and $\bw'_{\by}$ represent the net reaction vectors corresponding to the vertex $\by$ in $G$ and $G'$. From Remark <ref>, we deduce that

\begin{equation} \label{bw'y = 0}
\bw'_{\by} = \sum\limits_{\by \rightarrow\by_j\in E'} k'_{\by \rightarrow\by_j}(\by_j - \by) = \mathbf{0}.
\end{equation}

Since the network $G'$ is of Type II, each linkage of $G'$ has deficiency zero. By Proposition <ref>, we get that its vertices are affinely independent within each linkage class. This implies that the reaction vectors $\{ \by_i - \by \}_{\by \rightarrow\by_i \in E'}$ are linearly independent,  contradicting Equation (<ref>).

Assume that there exists a vertex $\by \in L \subseteq V$ such that $\by \notin V'$  (where $L$ is some linkage class in $G$).  Following the steps in the first part, we have that for any $\by' \in V'$, $\bw'_{\by} \neq \mathbf{0}$. This shows that $V' \subsetneq V$.
Moreover, from $\by \notin V'$, we get that

\begin{equation} \notag
\bw_{\by} = \mathbf{0}.
\end{equation}

This implies that the vertices in the linkage class $L$ are not affinely independent. Therefore the deficiency of linkage class $L$ is one. Since $(G, \bk)$ and $(G', \bk')$ have the same stoichiometric subspace and deficiency, we deduce that $G$ has at least one more linkage class than $G'$. From the Pigeonhole Principle, there exists at least one linkage class in $G'$ that is split into different linkage classes in $G$. Let us call this linkage class as $L'_1 \subsetneq V'$. Using Lemma <ref>, we have

\begin{equation} \notag
\text{dim}(\text{ker}(\bW'_1)) \geq 1.
\end{equation}

where $\bW'_1$ is the matrix of net reaction vectors on $L'_1$. This implies that the stoichiometric subspaces corresponding to the linkage classes in $G$ are not linearly independent, contradicting the fact that $G$ is of Type I.

Since $(G, \bk)$ and $(G', \bk')$ have the same stoichiometric subspace, number of vertices, and deficiency we obtain that $(G, \bk)$ and $(G', \bk')$ possess the same number of linkage classes. Finally, applying Theorem <ref>, we get that $(G, \bk)$ and $(G', \bk')$ cannot be dynamically equivalent, which leads to a contradiction.

The following remark is a direct consequence of Theorem <ref>.

For any mass-action system $(G, \bk)$, it at most has one type of weakly reversible deficiency one realization, i.e. either Type I or Type II.

§ MAIN RESULTS

This section aims to present the main algorithm of this paper, which checks the existence of a weakly reversible deficiency one realization and outputs one if it exists.
In this algorithm, the inputs are the matrices of source vertices and net reaction vectors via

\[
\bY_s = (\by_1,\by_2,\ldots,\by_m)
\ \ \text{and } \
\bW = (\bw_1,\bw_2,\ldots,\bw_m)
\]

For the sake of simplicity, we temporarily  let $d_{\bW}=\{\bd_1,\bd_2,\ldots,\bd_r\}$ denote the minimal set of generators of $\text{ker}(\bW)\cap\mathbb{R}^m_{\geq 0}$ in this section.

To build the main algorithm, we need an algorithm to search for a weakly reversible realization with a single linkage class.
We use the algorithm in [27] and summarize its main idea as follows.

First, the algorithm in [27] checks whether there exists a reaction network realization that generates the given dynamical system such that all the target vertices are among the source vertices, without imposing the restrictions that $(i)$ the network should be weakly reversible, and $(ii)$ there should be only one linkage class. Next, if such a realization exists, the algorithm greedily searches for a maximal realization (a realization containing the maximum number of reactions) that generates the same dynamical system, while still imposing the restriction that all target vertices are among the source vertices.
The algorithm uses the fact that if the initial realization was weakly reversible and consisted of a single linkage class, then the maximal realization found using this procedure preserves weakly reversibility and a single linkage class.
Finally, based on this maximal realization, the algorithm constructs a Kirchoff matrix $Q$ and checks whether $\dim(\text{ker}(Q))=1$ and $\text{supp}(\text{ker}(Q))=\{1,2,3,\ldots,m\}$.
If both conditions are satisfied, then the maximal realization is weakly reversible and consists of a single linkage class. Otherwise,  there is no such realization that generates the given polynomial dynamical system.

For more details on this algorithm and its implementation and complexity, please see [27]. In what follows, we will refer to the algorithm in [27] as Alg-WR$^{\ell=1}$.

§.§ Algorithm for weakly reversible and deficiency one realization

Now we state the main algorithm. The key idea is to find a proper decomposition on $\bW = \bY \bA_{\bk}$, which allows a weakly reversible and deficiency one realization. We apply

Alg-WR$^{\ell=1}$ to ensure weakly reversibility and  the single linkage class condition, and use results in Section <ref> to guarantee that the deficiency of the network is one.

(Check the existence of a weakly reversible deficiency one realization)

Input: The matrices of source vertices $\bY_s = (\by_1, \ldots, \by_m)$, and net reaction vectors $\bW = (\bw_1, \ldots, \bw_m)$ that generate the dynamical system $\dot{\bx} = \sum\limits_{i=1}^m \bx^{\by_i} \bw_i$.

Output: A weakly reversible deficiency one realization if exists or output that it does not exist.

[1]

Set flag = 0 and $\text{dim}(\text{ker}(\bW)) = \bw^*$;

Find the minimal set of generators $d_{\bW}=\{\bd_1,\bd_2,\ldots,\bd_r\}$ of the pointed cone $\text{ker}(\bW)\cap\mathbb{R}^m_{\geq 0}$.

$r < 2$ or
$\displaystyle\bigcup_{i=1}^r \text{supp} (\bd_i) \neq [m]$:
Exit the main program;

r = 2

Pass $\bY_s, \bW$ through

Alg-WR$^{\ell=1}$ outputs that a weakly reversible realization consisting of a single linkage class does not exist

Exit the main program;

Exit the main program;

$i = 1, 2, \ldots, r-1$
$j = i+1, i+2, \ldots, r$
$S_1 = \{\bd_i, \bd_j\}$.

$S_2 = \bd_{\bW} \setminus S_1$ and [For the simplicity of notations, we use a new symbol $\hbd$ to represent the vectors in $S_2$]set $S_2 := \{ \hbd_p \}^{r-2}_{p=1}$.

$r = \bw^*$ and the support of $S_1$ and every member of $S_2$ are disjoint

Define linkage classes to be $\{L_p\}^{r-1}_{p=1}$, where $L_p:= \{\text{supp}(\hbd_p), \hbd_p \in S_2\}$ for $1 \leq p \leq r-2$, and $L_{r-1} := \{\text{supp}(\bd_i) \cup \text{supp}(\bd_j) \}$.

Let $\bY_p$ denotes the vertices in linkage class $L_p$, and $\bW_p$ denotes the matrix of net reaction vectors corresponding to $\bY_p$.

$p=1$ to $r-1$

Pass $\bY_p,\bW_p$ through

Alg-WR$^{\ell=1}$ outputs that a weakly reversible realization consisting of a single linkage class does not exist

Go to line <ref>;

Exit the main program;

$r = \bw^* + 1$ and the support of the members of $S_2$ partition $[m]$

Define linkage classes to be $\{L_p\}^{r-2}_{p=1}$, where $L_p:= \{\text{supp}(\hbd_p), \hbd_p \in S_2\}$ for $1 \leq p \leq r-2$.

Let $\bY_p$ denotes the vertices in linkage class $L_p$, and $\bW_p$ denotes the matrix of net reaction vectors corresponding to $\bY_p$.

$\text{dim}(\text{ker}(\bW_1)) = \text{dim}(\text{ker}(\bW_2))  = \cdots = \text{dim}(\text{ker}(\bW_{r-2}))  =  1$

$p=1$ to $r-2$

Pass $\bY_p,\bW_p$ through

Alg-WR$^{\ell=1}$ outputs that a weakly reversible realization consisting of a single linkage class does not exist

Go to line <ref>;

Exit the main program;

End of main program

flag = 0

Print: No weakly reversible and deficiency one realization exists.

flag = 1

Print: Weakly reversible and deficiency one realization consisting of a single linkage class exists.

flag = 2

Print: Weakly reversible and deficiency one realization of Type I exists.

flag = 3

Print: Weakly reversible and deficiency one realization of Type II exists.

Node/.style=rectangle, rounded corners, draw=black, thick, fill=blue!10, fill opacity = 1, minimum width=4.5cm, minimum height=1cm, outer sep=0pt,
Edge/.style=very thick, style=very thick, arrows=-Stealth[length=7.5pt,width=7.5pt],

Algorithm <ref> for finding a weakly reversible deficiency one realization that generates a given polynomial dynamical system $\dot{\bx} = \sum\limits_{i=1}^m \bx^{\by_i} \bw_i$.

Now we show the correctness of Algorithm <ref> via the following two lemmas.

Suppose Algorithm <ref> exits with a positive flag value, then there exists a weakly reversible deficiency one realization of the dynamical system $\dot{\bx} = \sum\limits_{i=1}^m \bx^{\by_i} \bw_i$. Moreover, we have

(a) If flag = 1, the system admits a weakly reversible deficiency one realization consisting of a single linkage class.

(b) If flag = 2, the system admits a weakly reversible deficiency one realization of Type I.

(c) If flag = 3, the system admits a weakly reversible deficiency one realization of Type II.

From flag $= 1$, we obtain that $r = 2$ with $\text{supp}(\bd_1) \cup \text{supp}(\bd_2) = [m]$. Moreover, the input matrices $\bY_s$ and $\bW$ pass through Alg-WR$^{\ell=1}$.

Then there exists a weakly reversible realization with a single linkage class that generates the dynamical system $\dot{\bx} = \sum\limits_{i=1}^m \bx^{\by_i} \bw_i$. Using Remark <ref>, we conclude its deficiency is one.

From flag $= 2$, we get $r = \text{dim} (\text{ker}(\bW)) > 2$, and $r-1$ linkage classes $\{L_1, \ldots, L_{r-1}\}$ as follows.
There exists some $1 \leq i < j \leq r$,

\begin{equation} \notag
S_1 = \{\bd_i,\bd_j\} \ \ \text{and } \
S_2 = d_{\bW} \setminus S_1.
\end{equation}

For the simplicity of notations, we rename $r-2$ vectors in $S_2$ as
$S_2 := \{ \hbd_p \}^{r-2}_{p=1}$ and set

\begin{equation} \notag
\begin{split}
& L_p = \{\text{supp}(\hbd_p): \ \hbd_p \in S_2 \}, \ \text{for} \ 1 \leq p \leq r-2,
\\& L_{r-1} =\{\text{supp}(\bd_i) \cup \text{supp}(\bd_j) \},
\end{split}
\end{equation}

with $\{L_1, L_2, \ldots, L_{r-1} \}$ partition $[m]$.

Moreover, for any $1 \leq q \leq r-1$, the
matrices of source vertices and net reaction vectors $\bY_q, \bW_q$ related to the linkage class $L_q$ pass through

Thus each linkage class $L_q$ admits a weakly reversible realization.
Together with $\{L_q\}^{r-1}_{q=1}$ have disjoint supports in $[m]$, we have

\begin{equation}
\begin{split} \label{flag=2 ker estimate}
& \text{ker}(\bW_p) = \{\text{span}(\hbd_p) \}, \ \text{for} \ 1 \leq p \leq r-2,
\\& \text{ker}(\bW_{r-1}) = \text{span} \{ \bd_i, \bd_j \}.
\end{split}
\end{equation}

Using Lemma <ref> and Remark <ref> on the realization under Alg-WR$^{\ell=1}$,

we get

\begin{equation} \label{flag=2 deficiency}
\delta_1 = \cdots = \delta_{r-1} = 0
\ \ \text{and } \
\delta_{r-1} = 1,
\end{equation}

where $\delta_q$ represents the deficiency of linkage class $L_q$.

Now we compute the deficiency of the whole realization $\delta$. From (<ref>), we obtain

\begin{equation} \label{flag=2 ker sum}
\text{dim}(\text{ker}(\bW_1) + \text{dim}(\text{ker}(\bW_2) + \cdots + \text{dim}(\text{ker}(\bW_{r-1})
= r = \text{dim}(\text{ker}(\bW).
\end{equation}

Applying Lemma <ref> and Lemma <ref> on Equation (<ref>), we deduce for $p = 1, \ldots, r-1$,

\begin{equation} \notag
\text{dim}(\text{ker}(\bW_q)) = 1 + \delta_{q} = |L_q| - s_q \ \ \text{and } \
\sum\limits_{q=1}^{r-1} s_{q} = s,
\end{equation}

where $s_q$ and $s$ represent the stoichiometric subspace for linkage class $L_q$ and whole network respectively.
Then we do the summation from $q = 1$ to $r-1$, and get

\begin{equation} \notag
\sum\limits_{q=1}^{r-1} (|L_q| - s_q) = m - s = (r-1) + \sum\limits_{q=1}^{r-1} \delta_{q}.
\end{equation}

From Equation (<ref>), we conclude that

\begin{equation} \notag
\delta = m - s - (r-1) = \sum\limits_{q=1}^{r-1} \delta_{q} = 1,
\end{equation}

and the system admits a weakly reversible deficiency one realization of Type I.

From flag $= 3$, we get $r = \text{dim} (\text{ker}(\bW)) + 1 > 2$, and $r-2$ linkage classes $\{L_1, \ldots, L_{r-2}\}$ as follows.
There exists some $1 \leq i < j \leq r$,

\begin{equation} \notag
S_1 = \{\bd_i,\bd_j\}
\ \ \text{and } \
S_2 = d_{\bW} \setminus S_1,
\end{equation}

Similarly, we rename $r-2$ vectors in $S_2$ as $S_2 := \{ \hbd_p \}^{r-2}_{p=1}$ and set

\begin{equation} \notag
L_p = \{\text{supp}(\hbd_p): \ \hbd_p \in S_2 \}, \ \text{for} \ 1 \leq p \leq r-2,
\end{equation}

with $\{L_1, L_2, \ldots, L_{r-2} \}$ partition $[m]$.

Moreover, for any $1 \leq p \leq r-2$,  the
matrices of source vertices and net reaction vectors $\bY_p, \bW_p$
related to the linkage class $L_p$ pass through Alg-WR$^{\ell=1}$.

Thus each linkage class $L_q$ admits a weakly reversible realization with $\text{dim}(\text{ker}(\bW_p)) = 1$.
Applying that $\{L_p\}^{r-2}_{p=1}$ partition $[m]$, we have

\begin{equation}
\label{flag=3 ker estimate}
\text{ker}(\bW_p) = \{\text{span}(\hbd_p) : \hbd_p \in S_2 \}, \ \text{for} \ 1 \leq p \leq r-2.
\end{equation}

Using Lemma <ref> on the realization under

we get

\begin{equation} \label{flag=3 deficiency}
\delta_1 = \cdots = \delta_{r-2} = 0,
\end{equation}

where $\delta_p$ represents the deficiency of linkage class $L_p$.

Now we compute the deficiency of the whole realization $\delta$. From (<ref>), we obtain

\begin{equation} \label{flag=3 ker sum}
\text{dim}(\text{ker}(\bW_1) + \text{dim}(\text{ker}(\bW_2) + \cdots + \text{dim}(\text{ker}(\bW_{r-2})
= r-2 = \text{dim}(\text{ker}(\bW) - 1.
\end{equation}

Applying Lemma <ref> and Lemma <ref> on Equation (<ref>), we deduce for $p = 1, \ldots, r-2$,

\begin{equation} \notag
\text{dim}(\text{ker}(\bW_p)) = 1 + \delta_{p} = |L_p| - s_p
\ \ \text{and } \
\sum\limits_{p=1}^{r-2} s_{p} = s + 1,
\end{equation}

where $s_p$ and $s$ represent the stoichiometric subspace for linkage class $L_p$ and whole network respectively.
Summing from $p = 1$ to $r-2$, we get

\begin{equation} \notag
\sum\limits_{p=1}^{r-1} (|L_p| - s_p) = m - (s + 1) = (r-2) + \sum\limits_{q=1}^{r-1} \delta_{q}.
\end{equation}

From Equation (<ref>), we conclude that

\begin{equation} \notag
\delta = m - s - (r-2) = \sum\limits_{q=1}^{r-1} \delta_{q} + 1 = 1,
\end{equation}

and the system admits a weakly reversible deficiency one realization of Type II.

Suppose the dynamical system $\dot{\bx} =\displaystyle \sum_{i=1}^m \bx^{\by_i} \bw_i$ admits a weakly reversible deficiency one realization, then Algorithm <ref> must set the flag value to be either $1,2$ or $3$.

Note that every weakly reversible deficiency one network belongs to the following:

* Weakly reversible deficiency one realization consisting of a single linkage class.

* Weakly reversible deficiency one realization of Type I, with two or more  linkage classes.

* Weakly reversible deficiency one realization of Type II.

Therefore, we split our proof into the above three cases.

Case 1:
Suppose the system admits a weakly reversible deficiency one realization consisting of a single linkage class.
From Lemma <ref>, we obtain that $r = 2$ and the input $\bY$ and $\bW$ pass through

from the weakly reversibility. Therefore, Algorithm <ref> will exit with flag $=1$.

Case 2:
Suppose the system admits a weakly reversible deficiency one realization of Type I with $\ell > 1$ linkage classes, denoted by $L_1, L_2, \ldots, L_{\ell}$.
From Lemma <ref>, we have

\begin{equation} \notag
\text{dim}(\text{ker}(\bW_{\ell}))= 2
\ \ \text{and } \
\text{dim}(\text{ker}(\bW_p)) = 1, \ \text{for} \ 1 \leq p \leq \ell-1,
\end{equation}

\begin{equation} \notag
\text{dim}(\text{ker}(\bW)) = \text{dim}(\text{ker}(\bW_1)) + \cdots + \text{dim}(\text{ker}(\bW_{\ell})) = \ell + 1 = r.
\end{equation}

Moreover, there exist $\ell + 1$ vectors  $\bd_1, \ldots, \bd_{\ell + 1}$ forming the minimal set of generators of $\text{ker}(\bW)\cap\mathbb{R}^m_{\geq 0}$, such that for
$p=1, \ldots, \ell-1$,

\begin{equation} \notag
\text{supp}(\bd_p) = L_p,
\end{equation}

\begin{equation} \notag
\text{supp}(\bd_{\ell}) \subsetneq L_{\ell}, \ \
\text{supp}(\bd_{\ell+1}) \subsetneq L_{\ell}, \ \
\text{supp}(\bd_{\ell})\cup \text{supp}(\bd_{\ell +1}) = L_{\ell}.
\end{equation}

Thus, when $i = \ell$ and $j = \ell +1$,
(i.e. $\bd_i = \bd_{\ell}$ and $\bd_j = \bd_{\ell+1}$), Algorithm <ref> will exit with flag $=2$.

Case 3:
Suppose the system admits a weakly reversible deficiency one realization of Type II with $\ell > 1$ linkage classes, denoted by $L_1, L_2, \ldots, L_{\ell}$.
From Lemma <ref>, we have

\begin{equation} \notag
\text{dim}(\text{ker}(\bW_p)) = 1, \ \text{for} \ 1 \leq p \leq \ell,
\end{equation}

\begin{equation} \notag
\text{dim}(\text{ker}(\bW)) - 1 = \text{dim}(\text{ker}(\bW_1)) + \cdots + \text{dim}(\text{ker}(\bW_{\ell})) = \ell = r - 2.
\end{equation}

Moreover, there exist $\ell + 2$ vectors  $\bd_1, \ldots, \bd_{\ell + 2}$ forming the minimal set of generators of $\text{ker}(\bW)\cap\mathbb{R}^m_{\geq 0}$, such that for
$p=1, \ldots, \ell$,

\begin{equation} \notag
\text{supp}(\bd_p) = L_p, \ \
\text{supp}(\bd_{\ell +1}) \cap L_p \subsetneq L_p, \ \
\text{supp}(\bd_{\ell +2}) \cap L_p \subsetneq L_p.
\end{equation}

Again when we pick $i = \ell$ and $j = \ell +1$, Algorithm <ref> will exit with flag $=3$.

Lastly, we show every mass-action system admitting a weakly reversible deficiency one realization has a unique flag value after applying Algorithm <ref>.
Following Remark <ref>, we deduce that
if flag $=3$ after passing the same mass-action system through the algorithm, the flag value cannot equal $1$ or $2$.
From Lemma <ref> and Lemma <ref>, we have $r = 2$ if flag $=1$, and $r = \ell + 1 > 2$ if flag $=2$.
Thus, it is also impossible that the flag equals both $1$ and $2$ on the same mass-action system.
Therefore, we show the uniqueness and prove this lemma.

The following remark is a direct consequence of Lemma <ref>.

If Algorithm <ref> sets the value of flag to 0, then $\dot{\bx} = \sum\limits_{i=1}^m \bx^{\by_i} \bw_i$ does not admit a weakly reversible deficiency one realization.

Consider the system of differential equations

\begin{equation} \label{eq: alg2 pass}
\begin{split}
\dot{x} & = - x +x^2, \\
\dot{y} & =  x^3 - x^3y^2.
\end{split}
\end{equation}

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

\begin{equation}
\bY_s =  \begin{pmatrix}
1  & 2  & 3 & 3 &3\\
0 	& 0 & 0 & 1  & 2
\end{pmatrix}, \ \text{and } \
\bW =  \begin{pmatrix}
1  & -1  & 0 & 0 & 0\\
0 	& 0 & 1 & 0 & -1
\end{pmatrix}.
\end{equation}

respectively, which are inputs to Algorithm <ref>.

Then, we can compute that $\text{dim}(\text{ker}(\bW)) = 3$, and extreme vectors of $\text{ker}(\bW)\cap\mathbb{R}^m_{\geq 0}$ is given by

\begin{equation} \notag
\bd_1 =
\begin{pmatrix}
1  \\
1 \\
0  \\
0 \\
0 \\
\end{pmatrix}, \ \
\bd_2=\begin{pmatrix}
0  \\
0 \\
0  \\
1 \\
0 \\
\end{pmatrix}, \ \
\bd_3=\begin{pmatrix}
0  \\
0 \\
1  \\
0 \\
1 \\
\end{pmatrix}.
\end{equation}

This shows that $r =3$, and the algorithm enters line <ref>.

Next, when we pick $i = 2$, $S_1= \{\bd_2,\bd_3\}$ and $S_2= \{\bd_1\}$. Note that $r = \text{dim}(\text{ker}(\bW)) = 3$, and the support of $S_1$ and every member of $S_2$ are disjoint, the algorithm defines candidate linkage classes are follows:

\begin{equation} \notag
L_1 = \{\text{supp}(\bd_1)\} = \{1, 2\}, \ \
L_2 = \{\text{supp}(\bd_2) \cup \text{supp}(\bd_3) \} = \{3, 4, 5\}.
\end{equation}

Following the candidate linkage classes $L_1, L_2$, we derive the corresponding matrices of source vertices and net direction vectors:

\begin{equation} \notag
\begin{split}
& \bY_1 =  \begin{pmatrix}
1  & 2 \\
0  & 0
\end{pmatrix}, \ \
\bW_1 =  \begin{pmatrix}
1  & -1  \\
0  &  0
\end{pmatrix}, \ \text{and } \
\bY_2 =  \begin{pmatrix}
3 & 3 & 3 \\
0 & 1 & 2
\end{pmatrix}, \ \
\bW_2 =  \begin{pmatrix}
0  &  0 & 0 \\
1  &  0 & -1
\end{pmatrix}.
\end{split}
\end{equation}

After that, we pass two pairs $(\bY_1, \bW_1)$ and $(\bY_2, \bW_2)$ through

Both pairs pass successfully through

i.e., a weakly reversible single linkage class exists for both arrangements. Finally, the algorithm sets flag $=2$ on line <ref>, and exits.
Therefore, (<ref>) admits a weakly reversible deficiency one realization of Type I, whose E-graph is shown in Figure <ref>.

A weakly reversible deficiency one mass-action system from Example <ref>

Consider the system of differential equations

\begin{equation} \label{ex: alg2 fail}
\begin{split}
\dot{x} & = - x +x^2, \\
\dot{y} & =  0.
\end{split}
\end{equation}

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

\begin{equation}
\bY_s =  \begin{pmatrix}
1  & 2  \\
0 	& 0
\end{pmatrix}, \ \text{and } \
\bW = \begin{pmatrix}
-1  & 1 \\
0 	& 0
\end{pmatrix}.
\end{equation}

respectively, which are inputs to Algorithm <ref>.

Then, we compute that $\text{dim}(\text{ker}(\bW)) = 1$, and the extreme vector of $\text{ker}(\bW)\cap\mathbb{R}^2_{\geq 0}$ is

\begin{equation} \notag
\bd_1 =
\begin{pmatrix}
1  \\
\end{pmatrix}.
\end{equation}

This shows that $r=1$, then the algorithm satisfies the condition on line <ref> and exits the program with initial flag $= 0$.
Therefore, there doesn't exist any weakly reversible deficiency one realization for this system.

§.§ Implementation of Algorithm <ref>

In this section, we discuss how to implement Algorithm <ref>.
The algorithm is designed to find a weakly reversible deficiency one realization that generates the dynamical system $\dot{\bx} = \sum\limits_{i=1}^m \bx^{\by_i} \bw_i$, and it has three key steps:

* Compute $\dim (\ker(\bW))$ and $\dim (\ker(\bW_i))$.

* Find the extreme vectors of the cone $\text{ker}(\bW)\cap\mathbb{R}^m_{\geq 0}$.

* Pass pairs of the matrices
$\bY_s, \bW$ or $\bY_i, \bW_i$ through Alg-WR$^{\ell=1}$.

In Step 1, the implementation needs a rank-revealing factorization; we need to find a basis of $\bW$ or $\bW_{i}$, and then we can check the number of vectors in this basis.
This is equivalent to solving a linear programming problem.

In Step 2, we note that by the Minkowski-Weyl theorem [24, 25], there exists two representations of a polyhedral cone $C$ given by:

(a) H-representation:

There exists a matrix $A$, such that the cone $C$ can be written as

\begin{equation*}
C = \{A x\leq 0\}.
\end{equation*}

(b) V-representation:

The cone $C$ has the minimal set of generators $\{d_i\}$, such that

\begin{equation*}
C = \displaystyle\sum_{i=1}^r \lambda_i d_i,
\end{equation*}

where $\lambda_i \geq 0$.

To find the extreme vectors of the cone $\text{ker}
(\bW)\cap\mathbb{R}^m_{\geq 0}$,
we need a way to convert from the H-representation to the V-representation. There are two popular ways of performing this conversion:

(a) Double description method: This is an example of an incremental method, where the conversion from H-representation to V-representation is performed assuming that the solution to a smaller problem is already known [29]. In particular, let  $C(A):=\{A x\leq 0\}$. Let $J$ be a subset of the row indices of $A$. We will denote by $A_J$ the submatrix of $A$ obtained by selecting the $J$ rows of $A$. Let us assume that we have found the minimal set of generators for the cone $C(A_J)$. We will denote by $E$ the generating matrix whose columns are the extreme vectors of $C(A_J)$.  The double description algorithm selects an index $h$ that is not present in $J$ and constructs the generating matrix $E'$ that corresponds to the $A_{J+h}$. This is repeated for several iterations until the generating matrix for $C(A)$ is found. This algorithm is useful in cases where the inputs are degenerate and the dimension of the cone is small.

(b) Pivoting methods: In this method, the extreme vectors of the cone are found by the reverse search technique, where the simplex algorithm (that uses pivots iteratively) is run in reverse for the linear programming problem $Ax\leq 0$. The reverse search method determines the extreme vectors of the cone by building a tree in a depth-first-search fashion. This method was developed by Avis and Fukuda [30]. It is particularly useful for non-degenerate inputs where it runs in time polynomial of the input size.

In Step 3, we apply Alg-WR$^{\ell=1}$, and this step can be done by solving a sequence of linear programming problems. More details can be found in section 4.4 in [27].

§ DISCUSSION

Weakly reversible deficiency one networks are ubiquitous in biochemistry, and are known to have the capacity to exhibit sophisticated dynamics. Some notable examples include the Edelstein network, as in Example <ref>.
To better understand their dynamics, we divide them into two categories:  (i) Type I networks, where all linkage classes have deficiency zero except one linkage class having deficiency one, and (ii)  Type II networks, where all linkage classes have deficiency zero. The crucial quantity in the analysis of such networks is the pointed cone $\text{ker}(\bW)\cap\mathbb{R}^m_{\geq 0}$, where $\bW$ is the matrix formed by the net reaction vectors. In particular, the extreme vectors of this cone can be divided into two classes: cyclic generators and stoichiometric generators. Networks of Type I possess only cyclic generators and satisfy the conditions of the Deficiency One Theorem. Consequently, for Type I networks, there exists a unique steady state within every stoichiometric compatibility class. For Type II networks, the set of stoichiometric generators is not empty. The stoichiometric generators define subnetworks, such that if these subnetworks possess multiple steady states, then the original network also allows multiple steady states [26].

In addition, we show that networks of different types cannot be dynamically equivalent. Theorem <ref> establishes this fact, and this implies that any mass-action system, at most, has one type of weakly reversible deficiency one realization, either Type I or Type II. In Section <ref> we analyze in depth the extreme vectors of the cone $\text{ker}(\bW)\cap\mathbb{R}^m_{\geq 0}$ for weakly reversible deficiency one networks. In particular, we show that for Type I networks with $\ell$ linkage classes, there exist $\ell +1$ generators, while for Type II networks with $\ell$ linkage classes, there exist $\ell +2$ generators. Lemmas <ref> and <ref> establish these facts.

In Section <ref> we describe our main result: the construction and the proof of correctness of Algorithm <ref>.

This algorithm takes as input a matrix of source vertices and the corresponding matrix of net reaction vectors.

Algorithm <ref> uses

as a subroutine and determines whether or not there exists a weakly reversible deficiency one realization for this input.

It is interesting to put this algorithm in the context of existing algorithms in the literature. There has been seminal work in this direction [31, 32, 33, 34, 35, 36, 37]
based mostly on optimization methods that rely  on mixed integer linear programming to determine the existence of realizations of a certain type.

The algorithm in this paper uses a novel and straightforward geometric approach by focusing on the extreme vectors of the cone $\text{ker}(\bW)\cap\mathbb{R}^m_{\geq 0}$, instead of posing it as a constrained optimization problem.
Algorithm <ref> uses

and the properties of the extreme vectors of the cone $\text{ker}(\bW)\cap\mathbb{R}^m_{\geq 0}$ to determine the existence of weakly reversible deficiency one realizations. This geometric approach in both algorithms allows for a fully self-contained mathematical analysis of the correctness of these algorithms.

This work opens up interesting new avenues for future research. In particular, the relationship between the minimal set of generators of the cone $\text{ker}(\bW)\cap\mathbb{R}^m_{\geq 0}$ and the deficiency of the network can be explored in greater depth.

One could also explore the existence of mutually exclusive types of weakly reversible realizations for networks of higher deficiency.
Another possible direction would be to explore the geometry of this minimal set of generators for weakly reversible networks of higher deficiency.

[1]
E. Voit, H. Martens, and S. Omholt.
150 years of the mass action law.
PLOS Comput. Biol., 11(1):e1004012, 2015.

[2]
M. Feinberg.
Foundations of chemical reaction network theory.
Springer, 2019.

[3]
J. Gunawardena.
Chemical reaction network theory for in-silico biologists.
Notes available for download at http://vcp. med. harvard.
edu/papers/crnt. pdf, 2003.

[4]
P. Yu and G. Craciun.
Mathematical Analysis of Chemical Reaction Systems.
Isr. J. Chem., 58(6-7):733–741, 2018.

[5]
G. Craciun and C. Pantea.
Identifiability of chemical reaction networks.
J. Math. Chem., 44(1):244–259, 2008.

[6]
F. Horn and R. Jackson.
General mass action kinetics.
Arch. Ration. Mech. Anal., 47(2):81–116, 1972.

[7]
F. Horn.
Necessary and sufficient conditions for complex balancing in chemical
Arch. Ration. Mech. Anal., 49(3):172–186, 1972.

[8]
M. Feinberg.
Chemical oscillations, multiple equilibria, and reaction network
In Dynamics and modelling of reactive systems, pages 59–130.
Elsevier, 1980.

[9]
M. Feinberg.
Chemical reaction network structure and the stability of complex
isothermal reactors - I. the deficiency zero and deficiency one theorems.
Chem. Eng. Sci., 42(10):2229–2268, 1987.

[10]
M. Feinberg.
The existence and uniqueness of steady states for a class of chemical
reaction networks.
Arch. Ration. Mech. Anal., 132(4):311–370, 1995.

[11]
B. Boros.
Existence of positive steady states for weakly reversible mass-action
SIAM J. Math. Anal., 51(1):435–449, 2019.

[12]
M. Gopalkrishnan, E. Miller, and A. Shiu.
A geometric approach to the global attractor conjecture.
SIAM J. Appl. Dyn. Syst., 13(2):758–797, 2014.

[13]
D Anderson.
A proof of the global attractor conjecture in the single linkage
class case.
SIAM J. Appl. Math., 71(4):1487–1508, 2011.

[14]
G. Craciun, F. Nazarov, and C. Pantea.
Persistence and permanence of mass-action and power-law dynamical
SIAM J. Appl. Math., 73(1):305–329, 2013.

[15]
B. Boros and J. Hofbauer.
Permanence of weakly reversible mass-action systems with a single
linkage class.
SIAM J. Appl. Dyn. Syst., 19(1):352–365, 2020.

[16]
G. Craciun, J. Jin, and P. Yu.
An algorithm for weakly reversible deficiency zero realizations of
polynomial dynamical systems.
arXiv preprint arXiv:2205.14267, 2022.

[17]
C. Guldberg and P. Waage.
Studies Concerning Affinity.
CM Forhandlinger: Videnskabs-Selskabet I Christiana,
35(1864):1864, 1864.

[18]
M. Feinberg.
Lectures on chemical reaction networks.
Notes of lectures given at the Mathematics Research Center,
University of Wisconsin, page 49, 1979.

[19]
L. Adleman, M. Gopalkrishnan, M. Huang, P. Moisset, and D. Reishus.
On the mathematics of the law of mass action.
In A Systems Theoretic Approach to Systems and Synthetic Biology
I: Models and System Characterizations, pages 3–46. Springer, 2014.

[20]
M. Feinberg and F. Horn.
Chemical mechanism structure and the coincidence of the
stoichiometric and kinetic subspaces.
Arch. Ration. Mech. Anal., 66(1):83–97, 1977.

[21]
M. Feinberg.
Multiple steady states for chemical reaction networks of deficiency
Arch. Ration. Mech. Anal., 132(4):371–406, 1995.

[22]
A. Ferragut, C. Valls, and C. Wiuf.
On the liouville integrability of edelstein’s reaction system in
Chaos, Solitons & Fractals, 108:129–135, 2018.

[23]
G. Craciun, M. Johnston, G. Szederkényi, E. Tonello, J. Tóth, and
P. Yu.
Realizations of kinetic differential equations.
arXiv preprint arXiv:1907.07266, 2019.

[24]
D. Cox, J. Little, and H. Schenck.
Toric varieties.
Graduate Studies in Mathematics, Providence, R.I, 2011.

[25]
R. Rockafellar.
Convex analysis, volume 18.
Princeton university press, 1970.

[26]
C. Conradi, D. Flockerzi, J. Raisch, and J. Stelling.
Subnetwork analysis reveals dynamic features of complex (bio)
chemical networks.
Proceedings of the National Academy of Sciences,
104(49):19175–19180, 2007.

[27]
G. Craciun, A. Deshpande, and J. Jin.
Weakly reversible single linkage class realizations of polynomial
dynamical systems: an algorithmic perspective.
arXiv preprint arXiv:2302.13119, 2023.

[28]
A. Deshpande.
Source-only realizations, weakly reversible deficiency one networks
and dynamical equivalence.
arXiv preprint arXiv:2205.00801, 2022.

[29]
T. Motzkin, H. Raiffa, G. Thompson, and R. Thrall.
The double description method.
Contributions to the Theory of Games, 2(28):51–73, 1953.

[30]
D. Avis and K. Fukuda.
A pivoting algorithm for convex hulls and vertex enumeration of
arrangements and polyhedra.
In Proceedings of the seventh annual symposium on Computational
geometry, pages 98–104, 1991.

[31]
M. Johnston, D. Siegel, and G. Szederkényi.
Computing weakly reversible linearly conjugate chemical reaction
networks with minimal deficiency.
Mathematical biosciences, 241(1):88–98, 2013.

[32]
G. Lipták, G. Szederkényi, and K. Hangos.
Kinetic feedback design for polynomial systems.
J. Process Control, 41:56–66, 2016.

[33]
G. Szederkényi, G. Lipták, J. Rudan, and K. Hangos.
Optimization-based design of kinetic feedbacks for nonnegative
polynomial systems.
In 2013 IEEE 9th International Conference on Computational
Cybernetics (ICCC), pages 67–72. IEEE, 2013.

[34]
J. Rudan, G. Szederkényi, and K. Hangos.
Efficiently computing alternative structures of large biochemical
reaction networks using linear programming.
MATCH Commun. Math. Comput. Chem, 71:71–92, 2014.

[35]
J. Rudan, G. Szederkényi, K. Hangos, and T. Péni.
Polynomial time algorithms to determine weakly reversible
realizations of chemical reaction networks.
J. Math. Chem., 52(5):1386–1404, 2014.

[36]
G. Szederkényi, K. Hangos, and Z. Tuza.
Finding weakly reversible realizations of chemical reaction networks
using optimization.
arXiv preprint arXiv:1103.4741, 2011.

[37]
G. Lipták, G. Szederkényi, and K. Hangos.
Computing zero deficiency realizations of kinetic systems.
Systems & Control Letters, 81:24–30, 2015.