Weakly reversible single linkage class realizations of polynomial dynamical systems: an algorithmic perspective

TL;DR

This paper presents an algorithm to identify weakly reversible single linkage class realizations of polynomial dynamical systems, which guarantees persistent long-term behavior of the modeled species across various parameters.

Contribution

It introduces a novel algorithmic approach to find $W ext{R}^1$ realizations, aiding the analysis of stability and persistence in polynomial dynamical systems.

Findings

Algorithm successfully finds $W ext{R}^1$ realizations when they exist.

Ensures robustness of long-term dynamics and species persistence.

Applicable to systems in ecology, biochemistry, and epidemiology.

Abstract

Systems of differential equations with polynomial right-hand sides are very common in applications. In particular, when restricted to the positive orthant, they appear naturally (according to the law of mass-action kinetics) in ecology, population dynamics, as models of biochemical interaction networks, and models of the spread of infectious diseases. On the other hand, their mathematical analysis is very challenging in general; in particular, it is very difficult to answer questions about the long-term dynamics of the variables (species) in the model, such as questions about persistence and extinction. Even if we restrict our attention to mass-action systems, these questions still remain challenging. On the other hand, if a polynomial dynamical system has a weakly reversible single linkage class ($W\!R^1$) realization, then its long-term dynamics is known to be remarkably robust: all the variables are persistent (i.e., no species goes extinct), irrespective of the values of the parameters in the model. Here we describe an algorithm for finding $W\!R^1$ realizations of polynomial dynamical systems, whenever such realizations exist.

Tables

Table 1. Table 1:

Reactions in $G$	Reactions in $\tilde{G}$
$X_1+X_2\stackrel{k_1}{\to }2X_1+X_2$	$X_1+2X_2\stackrel{k_1}{\to }2X_1+X_2$ $X_1+X_2+X_3\stackrel{k_1}{\to }2X_1+X_2$
$X_2+X_3\stackrel{k_2}{\to }2X_2+X_3$	$X_2+2X_3\stackrel{k_2}{\to }2X_2+X_3$ $X_1+X_2+X_3\stackrel{k_2}{\to }2X_2+X_3$
$X_3+X_1\stackrel{k_3}{\to }2X_3+X_1$	$2X_1+X_3\stackrel{k_3}{\to }2X_3+X_1$ $X_1+X_2+X_3\stackrel{k_3}{\to }2X_3+X_1$
$X_1+X_2\stackrel{k_4}{\to }{X}_1+X_2+X_3$	$X_1+2X_2\stackrel{k_4}{\to }{X}_1+X_2+X_3$ $2X_1+X_2\stackrel{k_4}{\to }{X}_1+X_2+X_3$
$X_2+X_3\stackrel{k_5}{\to }{X}_1+X_2+X_3$	$X_2+2X_3\stackrel{k_5}{\to }{X}_1+X_2+X_3$ $2X_2+X_3\stackrel{k_5}{\to }{X}_1+X_2+X_3$
$X_1+X_3\stackrel{k_6}{\to }{X}_1+X_2+X_3$	$X_1+2X_3\stackrel{k_6}{\to }{X}_1+X_3+X_3$ $2X_1+X_3\stackrel{k_6}{\to }{X}_1+X_2+X_3$

Equations

$$\begin{array}{rl}\frac{d{x}_1}{dt}& =p_1(x_1,\dots ,x_n),\\ \frac{d{x}_2}{dt}& =p_2(x_1,\dots ,x_n),\\ & ⋮\\ \frac{d{x}_n}{dt}& =p_n(x_1,\dots ,x_n),\end{array}$$

$${\mathbf{x}}^{\mathbf{y}}=x_1^{y_1}\dots x_n^{y_n},$$

$$S=\mathrm{span}\{{\mathbf{y}}^{\prime }-\mathbf{y}|\mathbf{y}\to {\mathbf{y}}^{\prime }\in \mathrm{E}\}.$$

$$\delta =m-\ell -s.$$

$$\mathbf{v}\cdot (\hat{\mathbf{y}}-\tilde{\mathbf{y}})>0,\text{and }\mathbf{v}\cdot \tilde{\mathbf{y}}<\mathbf{v}\cdot \mathbf{y}.$$

$$\mathbf{v}\cdot (\hat{\mathbf{y}}-\tilde{\mathbf{y}})>0,\mathbf{v}\cdot \tilde{\mathbf{y}}<\mathbf{v}\cdot \mathbf{y},\text{and }\mathbf{v}\cdot \tilde{\mathbf{y}}\le \mathbf{v}\cdot {\mathbf{y}}^{\ast }.$$

$$\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t}=\sum _{{\mathbf{y}}_i\to {\mathbf{y}}_j\in E}{k}_{{\mathbf{y}}_i\to {\mathbf{y}}_j}{\mathbf{x}}^{{\mathbf{y}}_i}({\mathbf{y}}_j-{\mathbf{y}}_i).$$

$$\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t}=\sum _{{\mathbf{y}}_i\to {\mathbf{y}}_j\in E}{k}_{{\mathbf{y}}_i\to {\mathbf{y}}_j}({\mathbf{x}}^{\ast }{)}^{{\mathbf{y}}_i}({\mathbf{y}}_j-{\mathbf{y}}_i)=0.$$

$$\frac{d\mathbf{x}}{dt}=\mathbf{Y}{\mathbf{A}}_{\mathbf{k}}{\mathbf{x}}^{\mathbf{Y}},$$

$$\mathbf{Y}=({\mathbf{y}}_1,{\mathbf{y}}_2,\dots ,{\mathbf{y}}_m),$$

$$[{\mathbf{A}}_{\mathbf{k}}{]}_{ji}=\{\begin{array}{ll}{k}_{{\mathbf{y}}_i\to {\mathbf{y}}_j},& \text{if }{\mathbf{y}}_i\to {\mathbf{y}}_j\in E,\\ -\sum _{{\mathbf{y}}_i\to {\mathbf{y}}_j\in E}{k}_{{\mathbf{y}}_i\to {\mathbf{y}}_j},& \text{if }i=j,\\ 0& \mathrm{otherwise},\end{array}$$

$${\mathbf{x}}^{\mathbf{Y}}=({\mathbf{x}}^{{\mathbf{y}}_1},{\mathbf{x}}^{{\mathbf{y}}_2},\dots ,{\mathbf{x}}^{{\mathbf{y}}_m}{)}^{⊺}.$$

$$\boldsymbol{e}_{p}=\begin{cases}\begin{array}[]{cl}[\boldsymbol{e}_{p}]_{i}>0,&\text{ if }\boldsymbol{y}_{i}\in T_{p},\\[5.0pt] [\boldsymbol{e}_{p}]_{i}=0,&\text{ otherwise.}\end{array}\end{cases}$$

$$\begin{array}{rl}& X\equiv {\mathbf{y}}_1=(1,0{)}^{⊺},2X\equiv {\mathbf{y}}_2=(2,0{)}^{⊺},3X\equiv {\mathbf{y}}_3=(3,0{)}^{⊺},\\ & 4X+Y\equiv {\mathbf{y}}_4=(4,1{)}^{⊺},5X+Y\equiv {\mathbf{y}}_5=(5,1{)}^{⊺},5X+2Y\equiv {\mathbf{y}}_6=(5,2{)}^{⊺}.\end{array}$$

$${\mathbf{A}}_{\mathbf{k}}=\left[\begin{array}{cccccc}-k_{12}& k_{21}& 0& 0& 0& 0\\ k_{12}& -k_{21}& 0& 0& 0& 0\\ 0& 0& -k_{34}& k_{43}& 0& 0\\ 0& 0& k_{34}& -k_{43}-k_{45}& 0& 0\\ 0& 0& 0& k_{45}& -k_{56}& k_{65}\\ 0& 0& 0& 0& k_{56}& -k_{65}\end{array}\right].$$

$${\mathbf{e}}_1=(k_{21},k_{12},0,0,0,0{)}^{⊺},{\mathbf{e}}_2=(0,0,0,0,k_{65},k_{56}{)}^{⊺}.$$

$$\text{supp}({\mathbf{e}}_1)=\{1,2\},\text{and }\text{supp}({\mathbf{e}}_2)=\{5,6\}.$$

$${\mathbf{w}}_i=\sum _{{\mathbf{y}}_i\to {\mathbf{y}}_j\in E}{k}_{{\mathbf{y}}_i\to {\mathbf{y}}_j}({\mathbf{y}}_j-{\mathbf{y}}_i),$$

$$\mathbf{W}=\left({\mathbf{w}}_1,{\mathbf{w}}_2,\dots ,{\mathbf{w}}_{m_s}\right).$$

$${\mathbf{w}}_i=\sum _{{\mathbf{y}}_i\to {\mathbf{y}}_j\in E}{k}_{{\mathbf{y}}_i\to {\mathbf{y}}_j}{\mathbf{y}}_j-\left(\sum _{{\mathbf{y}}_i\to {\mathbf{y}}_j\in E}{k}_{{\mathbf{y}}_i\to {\mathbf{y}}_j}\right){\mathbf{y}}_i.$$

$$\frac{d\mathbf{x}}{dt}=\mathbf{W}{\mathbf{x}}^{{\mathbf{Y}}_s},$$

$${\mathbf{x}}^{{\mathbf{Y}}_s}=({\mathbf{x}}^{{\mathbf{y}}_1},{\mathbf{x}}^{{\mathbf{y}}_2},\dots ,{\mathbf{x}}^{{\mathbf{y}}_{m_s}}{)}^{⊺}.$$

$$\mathrm{Im}(\mathbf{W})=S.$$

$$\mathbf{v}\in S,\text{and }\mathbf{v}\perp \mathbf{W}.$$

$$V_{\max }=\{{\mathbf{y}}_i\in V:\mathbf{v}\cdot {\mathbf{y}}_i=\underset{j}{\max }(\mathbf{v}\cdot {\mathbf{y}}_j)\}.$$

$$\mathbf{v}\cdot {\mathbf{w}}_1=\sum _{{\mathbf{y}}_j\in V}{k}_{{\mathbf{y}}_1\to {\mathbf{y}}_j}\mathbf{v}\cdot ({\mathbf{y}}_j-{\mathbf{y}}_1)\le k_{{\mathbf{y}}_1\to {\mathbf{y}}_2}\mathbf{v}\cdot ({\mathbf{y}}_2-{\mathbf{y}}_1)<0.$$

$$\underset{t\to \infty }{\mathrm{liminf}}{\mathbf{x}}_i(t)>0,\text{for }i=1,2,\dots ,n.$$

$$\mathbf{x}(t)\in D.$$

$$r_{i,j}=\{\begin{array}{ll}1,& \text{for }j\in S_i,j\ne i,\\ 0,& \text{for }j\notin S_i,\\ -\sum _{l\ne i}{r}_{i,l},& \text{for }j=i.\end{array}$$

$${\mathbf{B}}_i\mathbf{v}={\mathbf{w}}_i.$$

Full text

1 Introduction

By a system of differential equations with polynomial right-hand sides (or simply a polynomial dynamical system), we mean a dynamical system of the form

[TABLE]

where each $p_{i}(x_{1},\ldots,x_{n})$ is a polynomial in the variables $x_{1},\ldots,x_{n}$. In general, such systems are very difficult to analyze due to nonlinearities and feedbacks 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 [1]. Even the simplest object associated to (1), its steady state set, can give rise to highly nontrivial questions in real algebraic geometry.

Polynomial dynamical systems show up very often as standard models (based on mass-action kinetics) in biology, chemistry, population dynamics, infectious disease models, and many other areas of applications. In such models the variables $x_{i}$ represent populations, concentrations, or other quantities that cannot become negative, so the domain of (1) is restricted to the positive orthant. For example, in a biochemical network we may have the reaction $X_{1}+X_{2}\to X_{3}$, which consumes $X_{1}$ and $X_{2}$ and produces $X_{3}$; according to mass-action kinetics, this reaction contributes a negative monomial term of the form “$-kx_{1}x_{2}$” on the right-hand side of $\frac{dx_{1}}{dt}$ and $\frac{dx_{2}}{dt}$, and a positive monomial term “$kx_{1}x_{2}$” on the right-hand side of $\frac{dx_{3}}{dt}$, where $x_{1},x_{2},x_{3}$ denote the concentrations of the chemical species $X_{1},X_{2},X_{3}$. The parameter $k$ is called reaction rate constant. A reaction network consists of a set of such reactions, and if we add all these terms for all the reactions in the network (each one with its own reaction rate constant) we obtain standard dynamical system models for the network. In general, one cannot just rely on numerical simulations to deduce the dynamical properties of these models, because the values of the reaction rate constants cannot usually be estimated accurately. Therefore, it becomes very important to relate the structural properties of the reaction network with dynamical properties that can be generated by it [2, 3, 4, 5, 6, 7, 8, 9].

Alternatively, one may start with a system of the form (1) obtained from fitting some experimental data, with little or no information on the generating reaction network. In general, if a polynomial dynamical system is generated by some reaction network, then there are actually infinitely many other networks that also generate it [4]. This lack of unique identifiability of an underlying network can actually be leveraged to analyze the dynamics of a system of the form (1): if a network with certain properties can be found to generate it, then we may be able to immediately infer its dynamical behavior.

Some of the most important questions for polynomial systems (1) are related to the long-term dynamics of its solutions, which is usually analyzed in terms of the mathematical properties of persistence and permanence. The property of persistence means that no species can “go extinct”, i.e., for any solution $\boldsymbol{x}(t)$ of the system, we have $\displaystyle{\lim\inf}_{t\to\infty}\boldsymbol{x}_{i}(t)>0$ for all species $i$. The (stronger) property of permanence means that the system has a globally attracting compact set.

A class of networks whose long-term dynamics is best understood is the family of weakly reversible single linkage class networks [2]; here we call them simply “$W\!R^{1}$* networks*”, and we will refer to polynomial systems (1) that have $W\!R^{1}$ realizations as “$W\!R^{1}$* systems*”. Specifically, $W\!R^{1}$ systems have been shown to be persistent and permanent in a very robust way, which even allows for the explicit construction of globally attracting invariant sets [6, 10]. Moreover, complex balanced $W\!R^{1}$ systems have been shown to be globally stable, i.e., they have a globally attracting point within each linear invariant subspace [11].

Not only are the long-term dynamical properties of $W\!R^{1}$ systems well understood, but also their persistence and permanence properties hold for any choices of parameter values, 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 this paper, we describe an efficient algorithm for determining whether a given polynomial dynamical system (1) admits a $W\!R^{1}$ realization, and for finding such a realization whenever it exists.

Structure of the paper. In Section 2, we introduce some basic terminology of reaction networks. Primarily, we present the notion of net reaction vectors, which play a key role in the main algorithm. In Section 3, we propose Algorithm 1 to find if there exists a weakly reversible reaction network consisting of a single connected component that generates a given dynamical system. In Section 4, we discuss some special cases of weakly reversible realizations with a single linkage class and go through the steps in Algorithm 1 using several examples. Moreover, we illustrate how to implement this algorithm in practice. In Section 5, we summarize our findings in this paper and outline directions for future work.

Notation. We denote by $\mathbb{R}_{\geq 0}^{n}$ and $\mathbb{R}_{>0}^{n}$ the set of vectors in $\mathbb{R}^{n}$ with non-negative and positive entries respectively. Given two vectors $\boldsymbol{x}\in\mathbb{R}_{>0}^{n}$ and $\boldsymbol{y}\in\mathbb{R}^{n}$, we use the following notation for a monomial with exponents given by $\boldsymbol{y}$:

[TABLE]

where $\boldsymbol{x}=(x_{1},\ldots,x_{n})^{\intercal}$ and $\boldsymbol{y}=(y_{1},\ldots,y_{n})^{\intercal}$.

2 Reaction networks

Definition 2.1**.**

A reaction network, also called a Euclidean embedded graph (E-graph), is a directed graph $G=(V,E)$ in $\mathbb{R}^{n}$, where $V\subset\mathbb{R}^{n}$ is a finite set of vertices, $E\subseteq V\times V$ represents the set of edges, and such that there are neither self-loops nor isolated vertices in $G$. We denote the number of vertices by $m$, and let $V=\{\boldsymbol{y}_{1},\ldots,\boldsymbol{y}_{m}\}$. A directed edge $(\boldsymbol{y}_{i},\boldsymbol{y}_{j})\in E$ represents a reaction in the network, and is also denoted by $\boldsymbol{y}_{i}\rightarrow\boldsymbol{y}_{j}$. Moreover, we define the reaction vector associated with the edge $\boldsymbol{y}_{i}\rightarrow\boldsymbol{y}_{j}$ as $\boldsymbol{y}_{j}-\boldsymbol{y}_{i}\in\mathbb{R}^{n}$. Here $\boldsymbol{y}_{i}$ and $\boldsymbol{y}_{j}$ denote the source vertex and target vertex respectively.

Definition 2.2**.**

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:

[TABLE]

Moreover, for any positive vector $\boldsymbol{x}_{0}\in\mathbb{R}_{>0}^{n}$, the affine polyhedron $(\boldsymbol{x}_{0}+S)\cap\mathbb{R}^{n}_{>0}$ is called the stoichiometric compatibility class of $\boldsymbol{x}_{0}$.

Definition 2.3**.**

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

1. (a)

The set of vertices $V$ is partitioned by its connected components (also called linkage classes), which correspond to the subset of vertices belonging to that connected component. 2. (b)

A connected component $L\subseteq V$ is strongly connected, if every edge is part of an oriented cycle. Further, a strongly connected component $L\subseteq V$ is called a terminal strongly connected component if for every vertex $\boldsymbol{y}\in L$ and $\boldsymbol{y}\rightarrow\boldsymbol{y}^{\prime}\in E$, we have $\boldsymbol{y}^{\prime}\in L$. 3. (c)

$G=(V,E)$ is said to be weakly reversible, if every connected component is strongly connected, i.e., every edge is part of an oriented cycle.

Remark 2.4**.**

For a weakly reversible reaction network $G=(V,E)$, every vertex $\boldsymbol{y}\in V$ is a source and a target vertex. Furthermore, every linkage class is a strong linkage class, as well as a terminal strong linkage class.

Definition 2.5**.**

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

[TABLE]

Definition 2.6**.**

Consider a Euclidean embedded graph $G=(V,E)$, and denote by $V_{S}\subseteq V$ the set of source vertices in $G$. Then $G$ is said to be endotactic [12, 13, 14], if for every $\boldsymbol{v}\in\mathbb{R}^{n}$ and $\boldsymbol{y}\to\boldsymbol{y}^{\prime}\in E$ satisfying $\boldsymbol{v}\cdot(\boldsymbol{y}^{\prime}-\boldsymbol{y})<0$, there exists $\tilde{\boldsymbol{y}}\to\hat{\boldsymbol{y}}\in E$, such that

[TABLE]

Moreover, $G$ is said to be strongly endotactic [6, 13], if for every $\boldsymbol{v}\in\mathbb{R}^{n}$ and $\boldsymbol{y}\to\boldsymbol{y}^{\prime}\in E$ satisfying $\boldsymbol{v}\cdot(\boldsymbol{y}^{\prime}-\boldsymbol{y})<0$, there exists $\tilde{\boldsymbol{y}}\to\hat{\boldsymbol{y}}\in E$, such that for every ${\boldsymbol{y}}^{*}\in V_{S}$,

[TABLE]

Remark 2.7** ([13]).**

It can be shown that weakly reversible reaction networks are endotactic. Furthermore, if a network is weakly reversible and consists of a single linkage class, then it is strongly endotactic.

Figure 1 shows two examples of reaction networks. A reaction network can generate a wide range of dynamical systems. We are interested in mass-action kinetics, which has been extensively studied in  [2, 8, 15, 16, 17, 18].

Definition 2.8**.**

Let $G=(V,E)$ be a Euclidean embedded graph, we denote the vector of reaction rate constants by $\boldsymbol{k}=(k_{\boldsymbol{y}_{i}\rightarrow\boldsymbol{y}_{j}})_{\boldsymbol{y}_{i}\rightarrow\boldsymbol{y}_{j}\in E}\in\mathbb{R}_{>0}^{E}$, and $k_{\boldsymbol{y}_{i}\rightarrow\boldsymbol{y}_{j}}$ or $k_{ij}$ is called the reaction rate constant on the edge $\boldsymbol{y}_{i}\rightarrow\boldsymbol{y}_{j}$. The associated mass-action system generated by $(G,\boldsymbol{k})$ on $\mathbb{R}^{n}_{>0}$ is given by:

[TABLE]

Definition 2.9**.**

Consider the associated mass-action system generated by $(G,\boldsymbol{k})$ in (2). A point $\boldsymbol{x}^{*}\in\mathbb{R}_{>0}^{n}$ is called a positive steady state if

[TABLE]

It is well known that every mass-action system admits a matrix decomposition [19]. Hence, we can illustrate the mass-action system (2) in the following vectorial representation:

[TABLE]

where $\boldsymbol{Y}$ is a matrix whose columns are the vertices, defined as

[TABLE]

and $\boldsymbol{A}_{\boldsymbol{k}}$ is the negative transpose of the graph Laplacian of $(G,\boldsymbol{k})$, defined as

[TABLE]

and ${\boldsymbol{x}}^{\boldsymbol{Y}}$ is the vector of monomials given by

[TABLE]

In general, $\boldsymbol{Y}$ is called the matrix of vertices, and $\boldsymbol{A}_{\boldsymbol{k}}$ is called the Kirchoff matrix.

Here, we list one of the most important properties of the Kirchoff matrix $\boldsymbol{A}_{\boldsymbol{k}}$.

Theorem 2.10** ([20]).**

Let $(G,\boldsymbol{k})$ 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 $\{\boldsymbol{e}_{1},\boldsymbol{e}_{2},\ldots,\boldsymbol{e}_{t}\}$ for $\ker(\boldsymbol{A}_{\boldsymbol{k}})$, such that

[TABLE]

Example 2.11**.**

We will revisit the network shown in Figure 1(a) to verify Theorem 2.10. First, we set all vertices in the network as follows:

[TABLE]

The Kirchoff matrix of the network is given by:

[TABLE]

Using a direct computation, the following vectors form a basis for $\rm{ker}(\boldsymbol{A}_{\boldsymbol{k}})$:

[TABLE]

Thus, we have

[TABLE]

It is clear that the supports of two basis vectors relate to two terminal strongly connected components $\{X,2X\}$ and $\{5X+Y,5X+2Y\}$.

Motivated by the matrix decomposition in (4), we introduce a crucial concept: net reaction vector, and another matrix decomposition in terms of net reaction vectors, which play an important role in finding a realization.

Definition 2.12**.**

Consider a mass-action system $(G,\boldsymbol{k})$, and let $V_{S}=\{\boldsymbol{y}_{1},\boldsymbol{y}_{2},\ldots,\boldsymbol{y}_{m_{s}}\}\subseteq V$ be the set of source vertices of $G$. For each source vertex $\boldsymbol{y}_{i}\in V_{S}$, the net reaction vector $\boldsymbol{w}_{i}$ corresponding to $\boldsymbol{y}_{i}$ is given by:

[TABLE]

Moreover, we denote the matrix of net reaction vectors as follows:

[TABLE]

Following Definition 2.12, for each source vertex $\boldsymbol{y}_{i}\in V_{S}$, we can rewrite the corresponding net reaction vector $\boldsymbol{w}_{i}$ as

[TABLE]

Using a direct computation, we can rewrite the matrix decomposition in (4) as

[TABLE]

where ${\boldsymbol{x}}^{\boldsymbol{Y}_{s}}$ is the vector of monomials given by

[TABLE]

Further, we let $\boldsymbol{Y}_{s}=(\boldsymbol{y}_{1},\ \boldsymbol{y}_{2},\ \ldots,\ \boldsymbol{y}_{m_{s}})$ denote the matrix of source vertices, whose columns are the source vertices.

The following Lemma concerns the matrix of net reaction vectors when the mass-action system is weakly reversible.

Lemma 2.13**.**

Consider a weakly reversible mass-action system $(G,\boldsymbol{k})$ with vertices $\{\boldsymbol{y}_{i}\}_{i=1}^{m}$ and stoichiometric subspace $S$. Let $\{\boldsymbol{w}_{i}\}_{i=1}^{m}$ be the net reaction vectors of $G$, and $\boldsymbol{W}=(\boldsymbol{w}_{1},\boldsymbol{w}_{2},\ldots,\boldsymbol{w}_{m})$ be the matrix of net reaction vectors. Then we have

[TABLE]

Proof.

It is clear that $\rm{Im}(\boldsymbol{W})\subseteq S$ from Definition 2.12. Suppose that $\rm{Im}(\boldsymbol{W})\subset S$, then there exists a non-zero vector $\boldsymbol{v}$, such that

[TABLE]

Since $\mathbf{0}\neq\boldsymbol{v}\in S$, there exists a reaction $\boldsymbol{y}_{i}\to\boldsymbol{y}_{j}\in E$ such that $\boldsymbol{v}\cdot(\boldsymbol{y}_{j}-\boldsymbol{y}_{i})\neq 0$. This implies that the set $\{\boldsymbol{v}\cdot\boldsymbol{y}_{i}\}_{i=1}^{m}$ has at least two different numbers. Now let $V_{\max}$ be the subset of vertices which maximizes the dot product as follows.

[TABLE]

Since $G$ is weakly reversible, there exists an edge from a vertex in $V_{\max}$ to a vertex not belonging to it. Without loss of generality, let $\boldsymbol{y}_{1}\in V_{\max}$, and $\boldsymbol{y}_{1}\to\boldsymbol{y}_{2}\notin V_{\max}$ be this edge. Note that for all $i=1$, $2,\ldots,m$, we have $\boldsymbol{v}\cdot(\boldsymbol{y}_{i}-\boldsymbol{y}_{1})\leq 0$. Thus, we obtain

[TABLE]

This contradicts with $\boldsymbol{v}\perp\boldsymbol{W}$, and the result follows. ∎

At the end of this section, we introduce some important dynamical properties.

Definition 2.14**.**

Let $(G,\boldsymbol{k})$ be a mass-action system. Then $(G,\boldsymbol{k})$ is called persistent, if every solution $\boldsymbol{x}(t)$ with initial condition $\boldsymbol{x}(0)\in\mathbb{R}^{n}_{>0}$ satisfies the following:

[TABLE]

Definition 2.15**.**

Let $(G,\boldsymbol{k})$ be a mass-action system. Then $(G,\boldsymbol{k})$ is called permanent, if given any stoichiometric compatibility class $A$ and any solution $\boldsymbol{x}(t)$ with initial condition $\boldsymbol{x}(0)\in A$, there exists a time $T$ and a compact set $D\subset A$, such that for all $t>T$,

[TABLE]

Definition 2.16**.**

Let $(G,\boldsymbol{k})$ be a mass-action system. A point $\tilde{\boldsymbol{x}}\in\mathbb{R}^{n}_{>0}$ is said to be a global attractor within its stoichiometric compatibility class, if $\displaystyle\lim_{t\to\infty}\boldsymbol{x}(t)=\tilde{\boldsymbol{x}}$.

3 Main result

The goal of this section is to present the main algorithm of this paper: Algorithm 1, which searches for the existence of a weakly reversible realization consisting of a single linkage class. In particular, this algorithm outputs a maximal realization, whenever it exists. The input of Algorithm 1 is the matrix of source vertices $\boldsymbol{Y}_{s}=(\boldsymbol{y}_{1},\boldsymbol{y}_{2},\ldots,\boldsymbol{y}_{m})$, and the matrix of net reaction vectors $\boldsymbol{W}=(\boldsymbol{w}_{1},\boldsymbol{w}_{2},\ldots,\boldsymbol{w}_{m})$.

3.1 Algorithm for weakly reversible realization with a single linkage class

Here, we sketch the key idea behind Algorithm 1: given a weakly reversible realization $G=(V,E)$ with a single linkage class, adding new reactions among the vertices in $V$ on the realization $G$ preserves the properties of weak reversibility and single linkage class. We present this algorithm below and give proof of its correctness.

Now we show the correctness of Algorithm 1 via the following two Lemmas.

Lemma 3.1**.**

Suppose Algorithm 1 reaches line 23 and satisfies the conditions on line 23, then there exists a weakly reversible realization consisting of a single linkage class that generates the dynamical system $\dot{\boldsymbol{x}}=\sum\limits_{i=1}^{m}\boldsymbol{x}^{\boldsymbol{y}_{i}}\boldsymbol{w}_{i}$.

Proof.

Since the algorithm passes the condition on line 5  for $i=1,\ldots,m$, all net reaction vectors $\{\boldsymbol{w}_{i}\}^{m}_{i=1}$ can be realized by conical combinations of the vectors $B_{i,k}=(\boldsymbol{y}_{k}-\boldsymbol{y}_{i})$, with $1\leq k\leq m$.

For $i=1,\ldots,m$, we denote $S_{i}$ by the union of supports on some vectors $\boldsymbol{v}\in\mathbb{R}^{m}_{\geq 0}$, such that

[TABLE]

Suppose there exist $a_{i}$ distinct vectors $\boldsymbol{v}^{1},\ldots,\boldsymbol{v}^{a_{i}}\in\mathbb{R}^{m}_{\geq 0}$, with $S_{i}=\bigcup\limits^{a_{i}}_{q=1}\text{supp}(\boldsymbol{v}^{q})$. Then we consider the following vector:

[TABLE]

and it satisfies

[TABLE]

Recall that each $v^{i}$ represents one realization corresponding to the net reaction vector $\boldsymbol{w}_{i}$. Here we choose the vector $\tilde{\boldsymbol{v}}$ in (12), where we have weighted all vectors $\{\boldsymbol{v}^{q}\}^{a_{i}}_{q=1}$ equally in the graph. Hence, it is clear that the reaction $\boldsymbol{y}_{i}\to\boldsymbol{y}_{j}$ represented by $j\in S_{i}$ is included in the realization. Further, we note that scaling the reaction rates neither affects weak reversibility nor the number of linkage classes. Hence, the Kirchoff matrix can be designed from line 22.

Using Theorem 2.10, we know that the kernel of the Kirchoff matrix has a basis consisting of non-negative vectors whose supports are the terminal strongly connected components. Recall that since the algorithm satisfies the condition on line 23, we have

[TABLE]

This implies that all vertices corresponding to the Kirchoff matrix $Q$ are in the same terminal strongly connected component. Therefore, this realization is weakly reversible and consists of a single linkage class. ∎

Lemma 3.2**.**

Suppose there exists a weakly reversible realization consisting of a single linkage class on $\dot{\boldsymbol{x}}=\sum\limits_{i=1}^{m}\boldsymbol{x}^{\boldsymbol{y}_{i}}\boldsymbol{w}_{i}$, then Algorithm 1 must satisfy the conditions on lines 5 and 23.

Proof.

From the existence of a realization, all net reaction vectors $\{\boldsymbol{w}_{i}\}^{m}_{i=1}$ can be realized by conical combinations of the vectors $(\boldsymbol{y}_{k}-\boldsymbol{y}_{i})^{m}_{k=1}$. Thus, the algorithm must satisfy the condition on line 5, for $i=1,\ldots,m$. Now it suffices for us to show

[TABLE]

Here we claim that the realization produced by Algorithm 1 consists of the maximum number of reactions. To realize the system $\dot{\boldsymbol{x}}=\sum_{i=1}^{m}\boldsymbol{x}^{\boldsymbol{y}_{i}}\boldsymbol{w}_{i}$, we need to find a vector $\boldsymbol{v}$ for each vertex $\boldsymbol{y}_{i}$, such that

[TABLE]

First, for $i=1,\ldots,m$, we get a vector $\boldsymbol{v}^{*}\in\mathbb{R}^{m}_{\geq 0}$ from line 5, which solves Equation (15). After setting $\boldsymbol{v}^{*}_{i}=1$, we define the initial support set $S_{i}$ as follows:

[TABLE]

Next, we build an inner loop on $j=1,\ldots,m$. If $j\in S_{i}$, this implies that we already incorporated the reaction $\boldsymbol{y}_{i}\to\boldsymbol{y}_{j}$ in the realization. Otherwise, for each $j\notin S_{i}$, we further check whether there exists a vector $\boldsymbol{v}\in\mathbb{R}^{m}_{\geq 0}$, such that

[TABLE]

Once we find such vector $\boldsymbol{v}$, we update the set $S_{i}$ as

[TABLE]

This implies that whenever $j\in\text{supp}(\boldsymbol{v})$, we have $j\in S_{i}$.

After going through the whole inner loop, we obtain the complete version of set $S_{i}$. Then we follow the construction in (12), and it is clear that the reaction $\boldsymbol{y}_{i}\to\boldsymbol{y}_{j}$ represented by $j\in S_{i}$ is included in the realization.

Now suppose there is a vector $\hat{\boldsymbol{v}}$, solving Equation (15) and $\text{supp}(\hat{\boldsymbol{v}})\not\subseteq S_{i}$. This implies that there exists $j\in\text{supp}(\hat{\boldsymbol{v}})$ with $j\notin S_{i}$, a contradiction. Thus, the set $S_{i}$ contains the maximal number of positive entries.

From the claim above and line 22, we deduce that for any realization of the system, all reactions between $\{\boldsymbol{y}_{i}\}^{m}_{i=1}$ are included in the realization given by the Kirchoff matrix $Q$. Note that adding more reactions among the current vertices of a weakly reversible single linkage class network will preserve the properties of weak reversibility and the single linkage class condition. This implies that if there exists a weakly reversible realization consisting of a single linkage class, then the realization generated by $Q$ will also be weakly reversible and consist of a single linkage class. By Theorem 2.10, we conclude (14). ∎

The following remark is a direct consequence of Lemma 3.2.

Remark 3.3**.**

If Algorithm 1 fails at lines 5 or 23, then $\dot{\boldsymbol{x}}=\sum\limits_{i=1}^{m}\boldsymbol{x}^{\boldsymbol{y}_{i}}\boldsymbol{w}_{i}$ does not admit a weakly reversible realization with a single linkage class.

4 Special cases and the implementation of Algorithm 1

After showing Algorithm 1, we focus on some special cases and the implementation of the algorithm. We will discuss various properties of weakly reversible realizations consisting of a single linkage class but having different deficiencies, and the corresponding implementation of the algorithm.

The following Lemma allows us to compute the deficiency of the realization obtained from Algorithm 1.

Lemma 4.1**.**

Suppose that the dynamical system $\dot{\boldsymbol{x}}=\sum\limits_{i=1}^{m}\boldsymbol{x}^{\boldsymbol{y}_{i}}\boldsymbol{w}_{i}$ with $m$ vertices, and the matrix of net reaction vectors $\boldsymbol{W}=(\boldsymbol{w}_{1},\ldots,\boldsymbol{w}_{m})$ passes Algorithm 1 and outputs a weakly reversible realization consisting of a single linkage class. Then the deficiency of this realization is $m-1-\rm{Im}(\boldsymbol{W})$.

Proof.

From Lemma 2.13, we have $\rm{Im}(W)=S$. Therefore, the deficiency of realization obtained from Algorithm 1 is

[TABLE]

∎

4.1 Weakly reversible deficiency zero realizations consisting of a single linkage class

We first consider the case when a dynamical system admits a weakly reversible deficiency zero realization consisting of a single linkage class.

It is well known that weakly reversible deficiency zero networks are complex-balanced for any choice of positive rate constants [19]. In addition, for complex-balanced dynamical systems consisting of a single linkage class, there exists a globally attracting positive steady state within each stoichiometric compatibility class [11]. This leads to the following Lemma.

Lemma 4.2**.**

For a weakly reversible deficiency zero reaction network consisting of a single linkage class, every stoichiometric compatibility class admits a globally attracting positive steady state.

This is our primary motivation for finding weakly reversible deficiency zero realizations consisting of a single linkage class. Now we state the upcoming Lemma that relates these realizations to the existence of a vector in line 5 of Algorithm 1.

Lemma 4.3**.**

Consider a weakly reversible deficiency zero reaction network $G$ consisting of a single linkage class $L=\{\boldsymbol{y}_{1},\ldots,\boldsymbol{y}_{m}\}$. Let $\{\boldsymbol{w}_{1},\ldots,\boldsymbol{w}_{m}\}$ denote the net reaction vectors corresponding to these vertices. Define the matrix $\boldsymbol{B}_{i}\in\mathbb{R}_{n\times m}$, with $k^{\rm th}$ column $B_{i,k}:=(\boldsymbol{y}_{k}-\boldsymbol{y}_{i})$ for $1\leq k\leq m$. For each vertex $\boldsymbol{y}_{i}\in L$, there exists a unique vector $\boldsymbol{v}\in\mathbb{R}^{m}_{\geq 0}$, such that

[TABLE]

Proof.

For each vertex $\boldsymbol{y}_{i}\in L$, we write the stoichiometric subspace as

[TABLE]

Since $G$ has a single linkage class and deficiency zero, we obtain

[TABLE]

This shows that

[TABLE]

From weak reversibility and Lemma 2.13, we deduce that

[TABLE]

Since $\boldsymbol{B}_{i,i}=\mathbf{0}$, it is easy to see that $\mathbf{e_{i}}\in\ker(\boldsymbol{B}_{i})$ where $\mathbf{e_{i}}$ represents the unit vector in the $i$-th coordinate. Then, we have

[TABLE]

Since the net reaction vectors $\{\boldsymbol{w}_{i}\}^{m}_{i=1}$ come from the dynamics generated by network $G$, all of them can be realized. Applying (20), we conclude that Equation (19) has a unique solution for each vertex $y_{i}\in L$. ∎

Remark 4.4**.**

It is worth mentioning that if the dynamical system $\dot{\boldsymbol{x}}=\sum_{i=1}^{m}\boldsymbol{x}^{\boldsymbol{y}_{i}}\boldsymbol{w}_{i}$ admits a weakly reversible realization consisting of a single linkage class $L$, such that for each vertex $y_{i}\in L$, the Equation 5 in Algorithm 1 has a unique solution, such realization still can have a positive deficiency. For example, the network in Example 4.9 has a unique solution to Equation 5 for each vertex $y_{i}\in L$, but it has deficiency one.

Example 4.5**.**

Consider the matrices corresponding to the source vertices and net reaction vectors given by

[TABLE]

respectively, which are inputs to Algorithm 1. These inputs generate the following system of differential equations

[TABLE]

We have $n=2$ for two state variables $x,y$, and $m=3$ for two distinct monomials.

Next, applying line 5 in algorithm on $\boldsymbol{Y}_{s}=(\boldsymbol{y}_{1},\boldsymbol{y}_{2},\boldsymbol{y}_{3})$, we obtain

[TABLE]

and

[TABLE]

where $\boldsymbol{v}^{*}_{i}\in\mathbb{R}^{3}_{\geq 0}$ and $\boldsymbol{B}_{i}\boldsymbol{v}^{*}_{i}=\boldsymbol{w}_{i}$, for $i=1,2,3$.

Then, we can compute that for $i=1,2,3$,

[TABLE]

and derive

[TABLE]

Note that $\rm{dim(ker}(B_{i}))=1$. Together with Equation (23), we deduce that $\boldsymbol{v}^{*}_{i}$ is the unique solution to the equations $\boldsymbol{B}_{i}\boldsymbol{v}=\boldsymbol{w}_{i}$ and $\boldsymbol{v}_{i}=1$ for $i=1,2,3$. Therefore, we do not need to execute the inner loop given by lines 11-19 in Algorithm 1.

Following line 22, we construct the Kirchoff matrix:

[TABLE]

It is easy to check that $\ker(Q)=\text{span}\{(1,1,1)^{\intercal}\}$, and we deduce that

[TABLE]

Therefore, we conclude that the system given by (22) admits a weakly reversible realization with a single linkage class, whose E-graph is shown in Figure 3.

Example 4.6**.**

Consider the matrices of source vertices and net reaction vectors given by

[TABLE]

respectively, which are inputs to Algorithm 1. These inputs generate the following system of differential equations

[TABLE]

We have $n=2$ for two state variables $x,y$, and $m=2$ for two distinct monomials.

Next, following line 5 in algorithm on $\boldsymbol{Y}_{s}=(\boldsymbol{y}_{1},\boldsymbol{y}_{2})$, we obtain

[TABLE]

However, there does not exist a positive vector $\boldsymbol{v}^{*}$, which solves $\boldsymbol{B}_{1}\boldsymbol{v}^{*}=\boldsymbol{w}_{1}$. Therefore, there exists no weakly reversible realization consisting of a single linkage class that generates the dynamical system given by Equation 25.

4.2 Weakly reversible deficiency one realizations consisting of a single linkage class

In this section, we analyze the case when a dynamical system admits a weakly reversible deficiency one realization consisting of a single linkage class.

If a reaction network satisfies the conditions of the Deficiency One Theorem, then every stoichiometric compatibility class contains a unique positive steady state (if it exists) [21, 9]. On the other hand, for any weakly reversible network, there always exists a positive steady state within every stoichiometric compatibility class [22]. It is easy to check that every weakly reversible deficiency one network with a single linkage class must satisfy all conditions in the Deficiency One Theorem. Therefore, we get the following Lemma.

Lemma 4.7**.**

For a weakly reversible deficiency one network consisting of a single linkage class, there exists a unique positive steady state within every stoichiometric compatibility class.

This explains the importance of discovering weakly reversible deficiency one realizations with a single linkage class. Moreover, we introduce the next Lemma showing the existence of a vector in line 5 of Algorithm 1.

Lemma 4.8**.**

Consider a weakly reversible and deficiency one reaction network $G$ consisting of a single linkage class given by $L=\{\boldsymbol{y}_{1},\ldots,\boldsymbol{y}_{m}\}$. Let $\{\boldsymbol{w}_{1},\ldots,\boldsymbol{w}_{m}\}$ denote the net reaction vectors corresponding to these vertices. Define the matrix $\boldsymbol{B}_{i}\in\mathbb{R}_{n\times m}$, with $k^{\rm th}$ column $B_{i,k}:=(\boldsymbol{y}_{k}-\boldsymbol{y}_{i})$ for $1\leq k\leq m$. For each vertex $\boldsymbol{y}_{i}\in L$, the following system

[TABLE]

has at most two linearly independent solutions.

Proof.

For each vertex $y_{i}\in L$, we denote the stoichiometric subspace by $S$, such that

[TABLE]

Since $G$ has a single linkage class and deficiency one, we obtain

[TABLE]

This shows that

[TABLE]

From the Rank-Nullity Theorem, $\dim(\ker(\boldsymbol{B}_{i}))+\dim(\rm{Im}(\boldsymbol{B}_{i}))=m$. Since $G$ is weakly reversible, using Lemma 2.13, we get

[TABLE]

thus we obtain that $\dim(\ker(\boldsymbol{B}_{i}))=m-\dim(S)=2$.

Since $\boldsymbol{B}_{i,i}=\mathbf{0}$, we deduce $\ker(\boldsymbol{B}_{i})$ has one vector $\boldsymbol{u}\in\mathbb{R}^{m}$ such that $\boldsymbol{u}_{i}\neq 0$. Then, we have

[TABLE]

Since the net reaction vectors $\{\boldsymbol{w}_{i}\}^{m}_{i=1}$ come from the dynamics generated by the network $G$, all of them can be realized. Together with (32), the conclusion follows. ∎

Example 4.9**.**

Consider the matrices corresponding to the source vertices and net reaction vectors given by

[TABLE]

respectively, which are inputs to Algorithm 1. These inputs generate the following system of differential equations

[TABLE]

We have $n=2$ for two state variables $x,y$, and $m=4$ for four distinct monomials.

Next, applying line 5 in algorithm on $\boldsymbol{Y}_{s}=(\boldsymbol{y}_{1},\boldsymbol{y}_{2},\boldsymbol{y}_{3},\boldsymbol{y}_{4})$, we obtain

[TABLE]

and

[TABLE]

where $\boldsymbol{v}^{*}_{i}\in\mathbb{R}^{4}_{\geq 0}$ and $\boldsymbol{B}_{i}\boldsymbol{v}^{*}_{i}=\boldsymbol{w}_{i}$, for $1\leq i\leq 4$.

Then, we get the initial $S_{i}$ for $1\leq i\leq 4$,

[TABLE]

After executing the inner loop in lines 11-19, we do not have any update on $S_{i}$.

Now we follow line 22, and construct the Kirchoff matrix:

[TABLE]

It is easy to check that $\ker(Q)=\text{span}\{(1,1,1,1)^{\intercal}\}$, which shows

[TABLE]

Therefore, we conclude that (34) admits a weakly reversible realization with a single linkage class, whose E-graph is shown in Figure 4.

Example 4.10**.**

Consider the matrices corresponding to the source vertices and net reaction vectors given by

[TABLE]

respectively, which are inputs to Algorithm 1. These inputs generate the following differential equation

[TABLE]

We have $n=1$ for the state variables $x$, and $m=3$ for three distinct monomials.

Next, applying line 5 in algorithm on $\boldsymbol{Y}_{s}=(\boldsymbol{y}_{1},\boldsymbol{y}_{2},\boldsymbol{y}_{3})$, we obtain

[TABLE]

and

[TABLE]

where $\boldsymbol{v}^{*}_{i}\in\mathbb{R}^{3}_{\geq 0}$ and $\boldsymbol{B}_{i}\boldsymbol{v}^{*}_{i}=\boldsymbol{w}_{i}$, for $i=1,2,3$.

Then, we get the initial $S_{i}$ for $i=1,2,3$,

[TABLE]

Following the inner loop in lines 11-19, we can compute that

[TABLE]

After updating $S_{i}$ with $\boldsymbol{v}^{i}$ for $i=1,2,3$, we derive

[TABLE]

Now we follow line 22, and construct the Kirchoff matrix:

[TABLE]

It is easy to check that $\ker(Q)=\text{span}\{(1,1,1)^{\intercal}\}$, and we deduce that

[TABLE]

Therefore, we conclude (36) admits a weakly reversible realization with a single linkage class, whose E-graph is shown in Figure 5.

4.3 Weakly reversible realizations with a single linkage class have high deficiency

Now we list some properties of weakly reversible realizations of arbitrary positive deficiency consisting of a single linkage class.

Our motivation comes from autocatalytic networks, which are often associated with the context of the origin of life models [23, 24, 25, 26]. Owing to their autocatalytic nature, the concentrations of species in these networks can go unbounded. The crucial component in their analysis is the dynamics corresponding to the relative concentration of species. Given species $X_{1},X_{2},\ldots,X_{n}$ with concentrations $x_{1},x_{2},\ldots,x_{n}$, the relative concentration corresponding to species $X_{i}$ is given by $(\sum_{i=1}^{n}x_{i})^{-1}x_{i}$. It can be shown that for certain autocatalytic networks, the dynamics corresponding to the relative concentration of species can be generated by a reaction network [27].

In particular, we present an example of an autocatalytic network such that the network corresponding to the relative concentration of species is weakly reversible and consists of a single linkage class. Table 1 illustrates this fact. The left column of the table describes the reactions in $G$, which is an autocatalytic network. The right column of the table describes the reactions in $\tilde{G}$, which is the network corresponding to the relative concentration of species in $G$. The reactions in $\tilde{G}$ are obtained in the following way: for every reaction in $G$, there exists a corresponding pair of reactions in $\tilde{G}$ that is generated using  [27, Theorem 3.5]. In particular, the reactions in $\tilde{G}$ are generated by adding all possible species to the reactants of the corresponding reaction in $G$.

The network $\tilde{G}$ is depicted in Figure 6.(a). Note that the deficiency of $\tilde{G}$ is given by $\delta=7-1-3=3$. Using some modifications, we can construct a network shown in Figure 6.(b) which generates the dynamics as Figure 6.(a). Figure 6.(b) is a weakly reversible network consisting of a single linkage class. By [10], the dynamics generated by it is permanent. This implies that the dynamics generated by $\tilde{G}$ is also permanent.

4.4 Implementation of Algorithm 1

In this section, we discuss the implementation aspects of Algorithm 1. The algorithm is designed to find a weakly reversible realization consisting of a single linkage class for $\dot{\boldsymbol{x}}=\sum\limits_{i=1}^{m}\boldsymbol{x}^{\boldsymbol{y}_{i}}\boldsymbol{w}_{i}$, and it has three key steps:

Check for the existence of a vector $\boldsymbol{v}^{*}\in\mathbb{R}^{m}_{\geq 0}$, such that for $i=1,\ldots,m$,

[TABLE] 2. 2.

Check for the existence of a vector $\boldsymbol{v}\in\mathbb{R}^{m}_{\geq 0}$, such that for $i,j=1,\ldots,m$,

[TABLE] 3. 3.

Check $\dim(\ker(Q))=1$, and $\text{supp}(\ker(Q))=\{1,\ldots,m\}$.

In step 1, we compute the positive vector solving $\boldsymbol{B}_{i}\boldsymbol{v}^{*}=\boldsymbol{w}_{i}$ and consider the implementation as a sequence of linear programming problems. For $i=1,\ldots,m$, set the matrix $\boldsymbol{B}_{i}\in\mathbb{R}_{n\times m}$ as in line 4,

[TABLE]

From Lemma 3.2, if there exists some number $1\leq i\leq m$, such that there is no solution for (37), then the implementation fails. Therefore, no weakly reversible realization with a single linkage class exists.

In step 2, we compute the positive vector solving $\boldsymbol{B}_{i}\boldsymbol{v}=\boldsymbol{w}_{i}$ and $\boldsymbol{v}_{j}>0$. For each $j=1,\ldots,m$, if $j\notin S_{i}$, we do the following:

[TABLE]

It is clear that the solution to (38) is the desired vector if its $j$-th component is positive. Meanwhile, if the $j$-th component of the solution is zero, then implementation fails. Furthermore, we restrict $\boldsymbol{x}_{j}\leq 1$ to avoid the risk that $\boldsymbol{x}_{j}$ can be arbitrarily large.

Here we explain why adding the restriction on $j$-th component in (38) does not change the solvability of the problem. From $j\notin S_{i}$, there must exist a vector $\boldsymbol{x}^{*}\in\mathbb{R}^{m}_{\geq 0}$, such that

[TABLE]

Suppose there is a vector $\boldsymbol{x}\in\mathbb{R}^{m}_{\geq 0}$, which solves

[TABLE]

Then we can always find a sufficient small constant $\epsilon$ with $\boldsymbol{x}^{\epsilon}:=(1-\epsilon)\boldsymbol{x}^{*}+\epsilon\boldsymbol{x}$, such that

[TABLE]

This implies that if (39) admits a solution, there must exist another solution for (38).

In step 3, the implementation needs a rank-revealing factorization; we need to find a basis of $\ker(Q)$, and then we can check the number of vectors in this basis and their support. This is again solving a linear programming problem.

5 Discussion

Weakly reversible networks consisting of a single linkage class form an important class of networks, owing to the robust properties of the dynamical systems they generate. In particular, the dynamics produced by these networks (according to mass-action kinetics) is known to be persistent and permanent for all choices of reaction rate parameters [6, 10].

We describe an algorithm that determines if there exists a weakly reversible realization consisting of a single linkage class that generates a given polynomial dynamical system. Our input consists of two matrices: a matrix of source monomials and a matrix containing the corresponding net reaction vectors. The algorithm outputs a maximal weakly reversible realization consisting of a single linkage class (if one exists), which generates the dynamical system formed by the inputs. We also describe approaches for efficient implementations of this algorithm; in particular, we show that all the key steps in our algorithm reduce to solving simple linear programming problems.

Other approaches for finding weakly reversible realizations of polynomial dynamical systems are based primarily on mixed integer programming methods [28, 29, 30, 31]. The algorithm we describe here uses a simpler greedy approach, which works specifically because we are looking for realizations consisting of a single linkage class.

At the same time, our algorithm lays down the foundation for some future work. In particular, extending our algorithm to check the existence of more general realizations (e.g., weakly reversible realizations with multiple linkage classes that satisfy other desirable properties) is a potential avenue worthy of exploration. More specifically, the problem of finding weakly reversible realizations that satisfy the conditions of the Deficiency One Theorem [9] is a possibility that will explore in a follow-up paper [32].

Acknowledgements

This work was supported in part by the National Science Foundation grant DMS-2051568.

Data availability

Data sharing not applicable to this article as no datasets were generated or analyzed during the current study.