Polynomial Dynamical Systems, Reaction Networks, and Toric Differential Inclusions

TL;DR

This paper links polynomial dynamical systems to reaction network models using Euclidean embedded graphs, introduces toric differential inclusions, and explores their implications for longstanding conjectures in dynamical systems.

Contribution

It establishes a geometric framework connecting polynomial systems with reaction networks and introduces toric differential inclusions as a tool for analyzing global stability conjectures.

Findings

Polynomial dynamical systems can be represented by Euclidean embedded graphs.

Toric differential inclusions have a distinctive geometric structure.

Embedding polynomial systems into toric inclusions offers new approaches to global stability conjectures.

Abstract

Some of the most common mathematical models in biology, chemistry, physics, and engineering, are polynomial dynamical systems, i.e., systems of differential equations with polynomial right-hand sides. Inspired by notions and results that have been developed for the analysis of reaction networks in biochemistry and chemical engineering, we show that any polynomial dynamical system on the positive orthant $\mathbb R^n_{> 0}$ can be regarded as being generated by an oriented graph embedded in $\mathbb R^n$, called $\mathit{Euclidean \ embedded \ graph}$. This allows us to recast key conjectures about reaction network models (such as the Global Attractor Conjecture, or the Persistence Conjecture) into more general versions about some important classes of polynomial dynamical systems. Then, we introduce $\mathit{toric \ differential \ inclusions}$, which are piecewise constant autonomous dynamical systems with a remarkable geometric structure. We show that if a Euclidean embedded graph $G$ has some reversibility properties, then any polynomial dynamical system generated by $G$ can be embedded into a toric differential inclusion. We discuss how this embedding suggests an approach for the proof of the Global Attractor Conjecture and Persistence Conjecture.

Equations

$$\frac{dx}{dt}=\sum _{i=1}^m{x}^{s_i}{v}_i$$

$$\frac{dx}{dt}=\sum _{i=1}^m{k}_i{x}^{s_i}{v}_i$$

$$\frac{dx}{dt}=\sum _{i=1}^m{k}_i(t)x^{s_i}{v}_i$$

$$\frac{dx}{dt}=\sum _{e\in E}{k}_e{x}^{s(e)}v(e)$$

$$\frac{dx}{dt}=\sum _{e\in E}{k}_e(t)x^{s(e)}v(e)$$

$$V=\{s_i|i=1,...,m\}\cup \{s_i+v_i|i=1,...,m\},$$

$$E=\{(s_i,s_i+v_i)|i=1,...,m\}.$$

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

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

$$\frac{dx}{dt}=k_1(t)x_1^2\left(\genfrac{}{}{0ex}{}{-2}{1}\right)+k_2(t)x_2\left(\genfrac{}{}{0ex}{}{2}{-1}\right),$$

$$\frac{dx}{dt}=k_1(t)x^{s_1}\left(\genfrac{}{}{0ex}{}{-2}{1}\right)+k_2(t)x^{s_2}\left(\genfrac{}{}{0ex}{}{2}{-1}\right),$$

$$x_i(t)>{\epsilon }_0\mathrm{for}\mathrm{all}i\in \{1,...,n\}\mathrm{and}\mathrm{for}\mathrm{all}t\in I\cap [0,\infty ).$$

$$(x_0+S)\cap {\mathbb{R}}_{>0}^n,\mathrm{for}\mathrm{some}{x}_0\in {\mathbb{R}}_{>0}^n.$$

$$\underset{t\to \infty }{\lim }x(t)={\bar{x}}_0.$$

$$\sum _{e=(s,s^{\prime })\in E}{k}_e{\bar{x}}^s=\sum _{e=(s^{\prime },s)\in E}{k}_e{\bar{x}}^{s^{\prime }}.$$

$$\frac{dx}{dt}=\sum _{s\rightleftharpoons s^{\prime }\in E}\left(k_{s\to s^{\prime }}(t)x^s-k_{s^{\prime }\to s}(t)x^{s^{\prime }}\right)(s^{\prime }-s),$$

$$\frac{dx}{dt}=\left(k_{s\to s^{\prime }}(t)x^s-k_{s^{\prime }\to s}(t)x^{s^{\prime }}\right)(s^{\prime }-s).$$

$$F_{H,\delta }(X)=\{\begin{array}{ll}\text{hfill}{l}_{s^{\prime }-s}^{+}\text{hfill}& \text{ if dist(X,H)>δ and X⋅(s′−s)>0}\\ \text{hfill}{l}_{s^{\prime }-s}^{-}\text{hfill}& \text{ if dist(X,H)>δ and X⋅(s′−s)<0}\\ \text{hfill}{l}_{s^{\prime }-s}\text{hfill}& \text{ if dist(X,H)≤δ}\end{array}$$

$$\frac{dx}{dt}\in F_{H,\delta }(\log x),$$

$$\left(k_{s\to s^{\prime }}(t)x^s-k_{s^{\prime }\to s}(t)x^{s^{\prime }}\right)(s^{\prime }-s)$$

$$k_{s\to s^{\prime }}(t)x^s<k_{s^{\prime }\to s}(t)x^{s^{\prime }}.$$

$$(\log x)\cdot (s^{\prime }-s)>-2\log \epsilon .$$

$$(\log x)\cdot \frac{s^{\prime }-s}{||s^{\prime }-s||}>\delta .$$

$$C^o=\{y\in {\mathbb{R}}^n|x\cdot y\le 0\text{ for all }x\in C\}.$$

$$F_{H,\delta }(X)=\{\begin{array}{ll}\text{hfill}{H}_{+}^o\text{hfill}& \text{ if dist(X,H)>δ and X∈H+}\\ \text{hfill}{H}_{-}^o\text{hfill}& \text{ if dist(X,H)>δ and X∈H−}\\ \text{hfill}{H}_{+}^o+H_{-}^o\text{hfill}& \text{ if dist(X,H)≤δ},\end{array}$$

$$F_{H,\delta }(X)=\{\begin{array}{ll}\text{hfill}{H}_{+}^o\text{hfill}& \text{ if dist(X,H+)≤δ and dist(X,H−)>δ}\\ \text{hfill}{H}_{-}^o\text{hfill}& \text{ if dist(X,H−)≤δ and dist(X,H+)>δ}\\ \text{hfill}{H}_{+}^o+H_{-}^o\text{hfill}& \text{ if dist(X,H+)≤δ and dist(X,H−)≤δ}\end{array}$$

$$F_{H,\delta }(X)=\sum _{\begin{array}{c}C\in \{H_{+},H_{-},H\}\\ dist(X,C)\le \delta \end{array}}{C}^o$$

$$C_1^o+C_2^o=(C_1\cap C_2{)}^o$$

$$F_{H,\delta }(X)={\left(\bigcap _{\begin{array}{c}C\in \{H_{+},H_{-},H\}\\ dist(X,C)\le \delta \end{array}}C\right)}^o$$

$$\frac{dx}{dt}\in F_{\mathcal{F},\delta }(\log x),$$

$$F_{\mathcal{F},\delta }(X)=\sum _{\begin{array}{c}C\in \mathcal{F}\\ dist(X,C)\le \delta \end{array}}{C}^o.$$

Full text

1 Introduction

Many mathematical models in biology, chemistry, physics, and engineering are given by polynomial dynamical systems, or more generally, power-law dynamical systems [9, 10, 11, 12, 14, 15, 16, 17, 20, 21, 29, 30]. Almost always, these can be interpreted as population dynamics models, where the variables of interest are positive.

Any autonomous polynomial dynamical system (i.e., system of differential equations with polynomial right-hand side) on the strictly positive orthant $\mathbb{R}^{n}_{>0}$ can be represented as

[TABLE]

where $x=(x_{1},...,x_{n})\in\mathbb{R}^{n}_{>0}$, $s_{1},...,s_{m}$ are some vectors in $\mathbb{Z}^{n}_{\geq 0}$ called exponent vectors, $x^{s_{i}}$ denotes the monomial $x_{1}^{s_{i1}}x_{2}^{s_{i2}}...x_{n}^{s_{in}}$, and $v_{1},...,v_{m}$ are vectors in $\mathbb{R}^{n}$. A solution of (1) is a function $x:I\to\mathbb{R}^{n}_{>0}$ that satisfies (1), where $I$ is an interval in $\mathbb{R}$.

Note that, since the coordinates $x_{1},...,x_{n}$ are positive, the monomials $x^{s_{i}}$ are well-defined even if the coordinates of the exponent vectors $s_{i}$ are arbitrary real numbers (i.e., $s_{1},...,s_{m}$ are not necessarily in $\mathbb{Z}_{\geq 0}^{n}$). In that case we say that (1) is a power-law dynamical system. The approaches and results discussed in this paper apply not only to polynomial dynamical systems, but also to power-law dynamical systems. In this paper, whenever we say “polynomial dynamical system”, we mean “polynomial or power law dynamical system”.

In many applications there are also some positive parameter values in these systems, which may be difficult to estimate accurately (such as reaction rate constants in biochemistry and chemical engineering, or interaction rates in epidemiology and ecology). Then the dynamical system of interest may have the form

[TABLE]

where $k_{1},...,k_{m}$ are some positive constants. In this case, we may want to know if some properties of the solutions of the system (2) may hold for all choices of positive parameters $k_{i}$.

In other cases, the interaction network we need to model is part of a larger network that contains variables or “external factors” that influence our system, but are not contained in our system. In that case we cannot use an autonomous dynamical system as a model, but we may be able to use a nonautonomous dynamical system of the form

[TABLE]

where the functions $k_{i}$ are positive and uniformly bounded, i.e., there exists some $\varepsilon>0$ such that $\varepsilon\leq k_{i}(t)\leq\frac{1}{\varepsilon}$ for all $t$. We will refer to models of the form (3) as variable-$k$ polynomial dynamical systems.

In applications, there is great interest in understanding the global stability and persistence properties of dynamical systems of the form (2) and (3). For example, a natural question is the following: for what systems (3) is it true that all solutions have a positive lower bound for all $t>0$ (i.e., no variable “goes extinct”), irrespective of the choices of uniformly bounded external factors $k_{i}(t)$?

In this paper we describe an approach for analyzing such problems, even in the presence of unknown parameters (as in (2)) or external factors (as in (3)). In Section 2 we show that any polynomial dynamical system can be regarded as being generated by some “Euclidean embedded graph” (also called “E-graph”). In Section 3 we introduce the notion of “toric differential inclusion”, and we show that if an E-graph is reversible, then any (variable-$k$) polynomial dynamical system generated by it can be embedded into a toric differential inclusion. Then, in Section 4 we show that such an embedding still exist even if the reversibility restriction is relaxed significantly. In Section 5 we discuss how these embeddings may greatly simplify the analysis of some properties of polynomial dynamical systems.

2 Euclidean embedded graphs

A Euclidean embedded graph (or E-graph) is a finite oriented graph $\mathcal{G}=(V,E)$ whose vertices are labeled by distinct elements of $\mathbb{R}^{n}$ for some $n\geq 1$. With an abuse of notation, we identify the set $V$ with the set of vertex labels, i.e., we assume that $V\subset\mathbb{R}^{n}$. Moreover, we associate to each edge $e=(s,t)\in E$ its edge vector $v(e)=t-s$. Also, we define its source vertex to be $s(e)=s$, and its target vertex to be $t(e)=t$.

Given an Euclidean embedded graph $\mathcal{G}=(V,E)$, the polynomial dynamical systems generated by $\mathcal{G}$ are the dynamical systems on $\mathbb{R}^{n}_{>0}$ given by

[TABLE]

for some positive constants $k_{e}$. Note that if $V\subset\mathbb{Z}^{n}_{\geq 0}$ then (4) is just mass-action kinetics for a chemical reaction network represented by $\mathcal{G}$, i.e, one where there is a reaction of the form $s(e)\to t(e)$ for each edge $e\in E$. (Informally speaking, for mass-action systems the rate of each reaction is proportional to the product of the concentrations of all its reactants, i.e., the rate of the reaction $s(e)\to t(e)$ is proportional to $x^{s(e)}$; see [9, 16, 20, 29] for more details.)

More generally, the variable-k polynomial dynamical systems generated by $\mathcal{G}$ are the (nonautonomous) dynamical systems on $\mathbb{R}^{n}_{>0}$ given by

[TABLE]

such that there exists some $\varepsilon>0$ for which we have $\varepsilon\leq k_{e}(t)\leq\frac{1}{\varepsilon}$ for all $e\in E$ and for all $t$.

Let us note that for any variable-$k$ polynomial dynamical system (3) we can construct an E-graph $\mathcal{G}$ that generates it, and $\mathcal{G}$ is not unique. Assuming that the ordered pairs $(s_{1},v_{1}),(s_{2},v_{2}),...,(s_{m},v_{m})$ are distinct, the simplest way to construct such a $\mathcal{G}$ is to choose the set of vertices

[TABLE]

and the set of edges

[TABLE]

If we want to obtain a different E-graph that generates (3), we can, for example, write one of the vectors $v_{i}$ as a positive linear combination of two different nonzero vectors, and use these new vectors to obtain a graph with $m+1$ edges that also generates (3).

We will use E-graphs in order to try to identify the polynomial dynamical systems that are known to have (or are conjectured to have) important dynamical properties, such as persistence, permanence, and global stability. For this purpose, we first define some special kinds of E-graphs (namely, reversible and weakly reversible E-graphs), and then we focus our attention on polynomial dynamical systems that are generated by these special kinds of graphs.

We say that an E-graph $\mathcal{G}=(V,E)$ is reversible if for any edge $(s,t)$ in $E$ the reverse edge $(t,s)$ also belongs to $E$. Also we say that $\mathcal{G}$ is weakly reversible if any edge $(s,t)$ is part of an oriented cycle in $\mathcal{G}$, or equivalently, any connected component of $\mathcal{G}$ is strongly connected (a strongly connected directed graph is one where there exists a directed path between every pair of vertices). For example, graph in Fig. 1$(e)$ is reversible, and the graphs in Fig. 1$(c)$ and $(d)$ are weakly reversible. Of course, every reversible E-graph is also weakly reversible.

We say that a polynomial dynamical system is reversible if there exists some reversible E-graph that generates it, and, we say that a polynomial dynamical system is weakly reversible if there exists some weakly reversible E-graph that generates it. Analogously, we define variable-$k$ reversible and weakly reversible polynomial dynamical systems.

Example 1. Consider the dynamical system given by

[TABLE]

for some functions $k_{i}(t)$ with $\varepsilon<k_{i}(t)<\frac{1}{\varepsilon}$ for all $t$.

This system can be written in vector form, as follows:

[TABLE]

where $x=\binom{x_{1}}{x_{2}}$. In turn, this can be written in the form (3), as follows:

[TABLE]

where $s_{1}=\binom{2}{0}$ and $s_{2}=\binom{0}{1}$. Then, the simplest E-graph $\mathcal{G}$ that generates the dynamical system (6) has two edges, one edge going from $s_{1}$ to $s_{1}^{\prime}:=s_{1}+\binom{-2}{1}$, and the other edge going from $s_{2}$ to $s_{2}^{\prime}:=s_{2}+\binom{2}{-1}$. But, note that we happen to have $s_{1}^{\prime}=s_{2}$ and $s_{2}^{\prime}=s_{1}$, so the graph $\mathcal{G}$ actually has only two vertices, and is reversible. The graph $\mathcal{G}$ is shown in Fig. 2$(a)$ in Section 3. We will return to this kind of example in section 3, when we will see that the dynamics of this system can be understood by embedding it into a special kind of differential inclusion.

2.1 Persistent, permanenent, and globally stable polynomial dynamical systems

We say that a variable-$k$ polynomial dynamical system in $\mathbb{R}^{n}_{>0}$ is persistent if, for any solution $x(t)$ defined on an interval $I$ that contains $t=0$, there exists some $\varepsilon_{0}>0$ such that

[TABLE]

In other words, the system is called persistent if, for any solution $x(t)$ with positive initial condition, there exists a positive lower bound for all the variables $x_{i}(t)$, and for all the future times at which the solution is defined (we cannot say that there exists a positive lower bound for all $t>0$ for a technical reason: some solutions may blow up in finite time). Informally, persistence means that “no variable goes extinct”.

To define the permanence property we first need to point out that polynomial dynamical systems have some special invariant spaces. If a (variable-$k$) polynomial dynamical system is given by (4) or (5), then its edge space $S$ is the linear span of the set of edge vectors $\{v(e)\,|\,e\in E\}$. Then we define its affine invariant sets to be the sets of the form

[TABLE]

These are indeed invariant spaces for solutions of (4) or (5) on the domain $\mathbb{R}^{n}_{>0}$, because all the vectors that appear on the right-hand side of these equations are linear combinations of $v_{1},...,v_{m}$, so are contained in $S$.

We say that a (variable-$k$) polynomial dynamical system on $\mathbb{R}^{n}_{>0}$ is permanent if, for each affine invariant set $S_{x_{0}}=(x_{0}+S)\cap\mathbb{R}^{n}_{>0}$ there exists a compact set $K_{x_{0}}\subset S_{x_{0}}$ such that any solution $x(t)$ with $x(0)\in S_{x_{0}}$ can be extended for all $t>0$, and there exists some $t_{0}>0$ such that $x(t)\in K_{x_{0}}$ for all $t>t_{0}$. It follows that, for a permanent system, all solutions that start in $S_{x_{0}}$ can be extended for all $t>0$, and, for large enough $t$, they are bounded above and below by some positive constants (and these positive upper and lower bounds do not depend on the initial condition, while for persistent systems the lower bound $\varepsilon$ may depend on the initial condition). In particular, permanence implies persistence.

Also, we say that a polynomial dynamical system (2) has a globally attracting point within the affine invariant set $S_{x_{0}}$ if there exists a point $\bar{x}_{0}\in S_{x_{0}}$ such that any solution $x(t)$ with $x(0)\in S_{x_{0}}$ can be extended for all $t>0$ and we have

[TABLE]

Also, we say that a polynomial dynamical system (2) is vertex balanced if there exists an E-graph $\mathcal{G}=(V,E)$ that generates our system as in (4), and there exists a point $\bar{x}\in\mathbb{R}^{n}_{>0}$ such that for any vertex $s\in V$ we have

[TABLE]

In other words, if we think of the positive number $k_{e}x^{s(e)}$ as the rate of a flow along the edge $e$ (i.e., a flow from the vertex $s(e)$ to the vertex $t(e)$), then condition (9) says that, if $x=\bar{x}$, then at each vertex of the graph $G$, the sum of all the incoming flows equals the sum of all outgoing flows.

The notion of “vertex balanced polynomial dynamical system” (or “vertex balanced power-law dynamical system”) is a natural generalization of the notion of “toric dynamical system”, which in turn was a reformulation of the notion of “complex balanced mass-action system”, which was introduced by Fritz Horn and Roy Jackson in their seminal work on models of reaction networks with mass action kinetics  [21]. For more details see [8, 10, 16, 20]. This notion ultimately originates in the work of Boltzmann [5, 6]. For some recent connections between polynomial dynamical systems, reaction networks, and the Bolzmann equation, see [12].

2.2 Open problems

We can now formulate the following conjectures, inspired by analogous conjectures that have been formulated for mass-action systems [9, 10, 21, 22], and are widely regarded as the key open problems in this field.

Global Attractor Conjecture. Any vertex balanced polynomial dynamical system has a globally attracting point within any affine invariant set.

Extended Persistence Conjecture. Any variable-$k$ weakly reversible polynomial dynamical system is persistent.

Extended Permanence Conjecture. Any variable-$k$ weakly reversible polynomial dynamical system is permanent.

The global attractor conjecture is the oldest and best known of these conjectures, and has resisted efforts for a proof for over four decades, but proofs of many special cases have been obtained during this time, for example [1, 2, 9, 10, 24, 26, 27, 28]. The conjecture originates from the 1972 breakthrough work by Horn and Jackson [21], and was formulated by Horn in 1974 [22].

Recently, Craciun, Nazarov and Pantea [9] have proved the three-dimensional case of this conjecture, and Pantea has generalized this result for the case where the dimension of the linear invariant subspaces is at most three [24]. Using a different approach, Anderson has proved the conjecture under the additional hypothesis that the graph $G$ has a single connected component [2], and this result has been generalized by Gopalkrishnan, Miller, and Shiu for the case where the graph $G$ is strongly endotactic [19]. A proof of the global attractor conjecture in full generality (using as a main tool the embedding of weakly reversible polynomial dynamical systems into toric differential inclusions, which is the main topic of this paper) has been proposed in [8].

Note that all three conjectures above relate to weakly reversible polynomial dynamical systems. Indeed, it is known that if the vertex balance condition (9) is satisfied, then it follows that the E-graph $\mathcal{G}$ must be weakly reversible [10, 16, 21]. Moreover, all these conjectures are strongly related to some version of the persistence property; in particular, it is known that a proof of the Global Attractor Conjecture would follow if we could show that vertex balanced polynomial dynamical systems are persistent [8, 9, 27, 28].

In the next section we introduce toric differential inclusions, in order to facilitate the analysis of persistence properties of variable-$k$ weakly reversible polynomial dynamical systems. Indeed, we will see that the analysis of some properties of these nonautonomous systems can be reduced to the analysis of toric differential inclusions, which are not only autonomous (i.e., their right-hand sides are constant in $t$), but are also piecewise constant in $x$.

3 Toric differential inclusions

Given an E-graph $\mathcal{G}=(V,E)$, let us write $s\to s^{\prime}\in E$ if $(s,s^{\prime})$ is an edge of $\mathcal{G}$; also let us write $s\rightleftharpoons s^{\prime}\in E$ if both $(s,s^{\prime})$ and $(s^{\prime},s)$ are edges of $\mathcal{G}$.

Then, if $\mathcal{G}$ is reversible, the dynamical system (5) can be written as

[TABLE]

by grouping together pairs of terms given by an edge $s\to s^{\prime}$ and its reverse $s^{\prime}\to s$. In particular, if $\mathcal{G}$ consists of a single reversible edge $s\rightleftharpoons s^{\prime}$, then we obtain

[TABLE]

Note that we can understand the dynamics of the system (11), if we think of it as a “tug-of-war” between the forward and reverse terms, i.e., the positive and the negative monomials in (11). Indeed, both the forward and the reverse terms are trying to “pull” the state $x(t)$ of the system along the same line (parallel to the vector $s^{\prime}-s$), but in opposite directions. Recall that $\varepsilon<k_{e}(t)<\frac{1}{\varepsilon}$ for all $t$. Then, the domain $\mathbb{R}^{n}_{>0}$ can be partitioned into three regions: the region where the inequality $\varepsilon x^{s}>\frac{1}{\varepsilon}x^{s^{\prime}}$ holds (which implies $k_{s\to s^{\prime}}(t)x^{s}>k_{s^{\prime}\to s}(t)x^{s^{\prime}}$), the region where the inequality $\frac{1}{\varepsilon}x^{s}>\varepsilon x^{s^{\prime}}$ holds (which implies $k_{s\to s^{\prime}}(t)x^{s}>k_{s^{\prime}\to s}(t)x^{s^{\prime}}$), and an uncertainty region where neither one of these two inequalities are satisfied, and either one of the two terms $k_{s\to s^{\prime}}(t)x^{s}$ and $k_{s^{\prime}\to s}(t)x^{s^{\prime}}$ may win the tug-of-war, or there can be a tie, due to the fact that $k_{s\to s^{\prime}}(t)$ and $k_{s^{\prime}\to s}(t)$ may take any values between $\varepsilon$ and $\frac{1}{\varepsilon}$.

We will now show that the system (11) can be embedded into a piecewise constant differential inclusion defined using a partition of $\mathbb{R}^{n}$ into three corresponding (but simpler) regions, related to the ones above via a logarithmic transformation.

Indeed, let us define by $l_{s^{\prime}-s}$ the line through the origin in $\mathbb{R}^{n}$ generated by the vector ${s^{\prime}-s}$, and by $l_{s^{\prime}-s}^{-}\subset l_{s^{\prime}-s}$ the ray starting from the origin in the direction ${s^{\prime}-s}$, and by $l_{s^{\prime}-s}^{+}\subset l_{s^{\prime}-s}$ the ray starting from the origin in the direction ${s-s^{\prime}}$. Also, let us denote by $H$ the hyperplane through the origin orthogonal to ${s^{\prime}-s}$. For some $\delta>0$ define the set-valued function $F_{H,\delta}$ at $X\in\mathbb{R}^{n}$, as follows:

[TABLE]

Recall that $k_{s\to s^{\prime}}(t),k_{s^{\prime}\to s}(t)\in[\varepsilon,1/\varepsilon]$ for all $t$. We have:

Lemma 3.1**.**

The dynamical system (11) (which is given by an E-graph that consists of a single reversible edge $s\rightleftharpoons s^{\prime}$) is embedded in the differential inclusion

[TABLE]

where $\delta=\frac{2|\log\varepsilon|}{||s^{\prime}-s||}$.

Proof.

If $dist(\log x,H)\leq\delta$ there is nothing to be proved, because we already know that the right-hand side of (11), i.e., the vector

[TABLE]

belongs to $l_{s^{\prime}-s}$.

If $dist(\log x,H)>\delta$ and $(\log x)\cdot(s^{\prime}-s)>0$, we need to show that the right-hand side of (11) belongs to $l_{s^{\prime}-s}^{+}$, i.e., we need to show that

[TABLE]

For this, it is sufficient to show that $\frac{1}{\varepsilon}x^{s}<\varepsilon x^{s^{\prime}}$, which is equivalent to $x^{s^{\prime}-s}>\varepsilon^{-2}$, and, by taking logarithm on both sides of this inequality, it can also be written as

[TABLE]

On the other hand, $dist(\log x,H)$ is just the dot product between the vector $\log x$ and the unit vector that is orthogonal to $H$ and on the same side of $H$ as $\log x$. But, since $(\log x)\cdot(s^{\prime}-s)>0$, this unit vector is just $\frac{s^{\prime}-s}{||s^{\prime}-s||}$. Then the inequality $dist(\log x,H)>\delta$ implies

[TABLE]

Finally, given that $\delta=\frac{2|\log\varepsilon|}{||s^{\prime}-s||}$ and $\varepsilon\in(0,1)$, we can see that the inequalities (13) and (14) are equivalent.

The case where $dist(\log x,H)>\delta$ and $(\log x)\cdot(s^{\prime}-s)<0$ is analogous, so this concludes the proof. ∎

Therefore, if $\mathcal{G}$ consists of a single reversible edge and if $dist(\log x,H)>\delta$, then the right-hand side of the system (11) is a vector that is orthogonal to the hyperplane $H$ and points in the direction that goes from the point $\log x$ towards $H$ (see Fig. 2$(a)$-$(f)$ for some examples in $\mathbb{R}^{2}$). Next, we will use this observation in order to construct a generalization of Lemma 3.1 for all reversible E-graphs, by using the notion of polar cone.

Recall that a polyhedral cone $C\subset\mathbb{R}^{n}$ is the set of nonnegative linear combinations of a finite set of vectors in $\mathbb{R}^{n}$, or, equivalently, is a finite intersection of half-spaces in $\mathbb{R}^{n}$ [25]. For simplicity, we will often say cone instead of polyhedral cone, because the only cones we consider here are polyhedral cones.

Definition 3.1**.**

Consider a cone $C\subset\mathbb{R}^{n}$. The polar cone of $C$ is denoted $C^{o}$ and is given by

[TABLE]

The polar cone is just the negative of the better known dual cone. Also, note that if the cone $C$ is full-dimensional (i.e., its linear span is $\mathbb{R}^{n}$), then its polar cone $C^{o}$ is generated by the outer normal vectors of the codimension-1 faces of $C$. For more information about polyhedral cones and their polar (or dual) cones see [18, 25, 31].

Let us now return to the graph $\mathcal{G}$ that consists of a single reversible edge $s\rightleftharpoons s^{\prime}$. Recall that $H$ denotes the hyperplane through the origin and orthogonal to ${s^{\prime}-s}$. Denote by $H_{+}$ the closed half-space of $\mathbb{R}^{n}$ that is bounded by $H$ and contains the vector $s^{\prime}-s$, and denote by $H_{-}$ the closed half-space of $\mathbb{R}^{n}$ that is bounded by $H$ and contains the vector $s-s^{\prime}$. Note that the polar cone $H_{+}^{o}$ is equal to the ray $l_{s^{\prime}-s}^{+}$, and the polar cone $H_{-}^{o}$ is equal to the ray $l_{s^{\prime}-s}^{-}$. Then the set-valued function $F_{H,\delta}$ can be rewritten as

[TABLE]

where we define $A+B=\{a+b\,|\,a\in A\textrm{ and }b\in B\}$.

Then we can write $F_{H,\delta}(X)$ only in terms of the distance between $X$ and the half-spaces $H_{+}$ and $H_{-}$, as follows:

[TABLE]

In order to be able to construct a generalization of Lemma 3.1 for the case where the E-graph $\mathcal{G}$ contains several reversible edges, we want to regard the set of cones $\{H_{+},H_{-},H\}$ as a special case of a polyhedral fan, because if we have several reversible edges then we have to consider a cover of $\mathbb{R}^{n}$ given by several hyperplanes, the half-space they generate, and their intersections (see Fig. 2$(g)$-$(i)$). Recall the definition of a polyhedral fan [7, 18]:

Definition 3.2**.**

A finite set $\mathcal{F}$ of polyhedral cones in $\mathbb{R}^{n}$ is a polyhedral fan if the following two conditions are satisfied:

$(i)$* any face of a cone in $\mathcal{F}$ is also in $\mathcal{F}$,*

$(ii)$* the intersection of two cones in $\mathcal{F}$ is a face of both cones.*

We say that a polyhedral fan $\mathcal{F}$ is complete if  $\displaystyle\bigcup_{C\in\mathcal{F}}C=\mathbb{R}^{n}$.

For example, given a finite set $\mathcal{H}$ of hyperplanes that contain the origin, consider the set $\mathcal{F}_{\mathcal{H}}$ of all the possible intersections of half-spaces given by hyperplanes in $\mathcal{H}$. Then $\mathcal{F}_{\mathcal{H}}$ is a complete polyhedral fan. (Since we are only interested in complete polyhedral fans, from now on we will refer to them simply as fans.) This particular case of “hyperplane-generated fan” will be especially relevant for motivating our definition of toric differential inclusions.

Indeed, we can now write $F_{H,\delta}(X)$ as

[TABLE]

i.e., $F_{H,\delta}(X)$ consists of all possible sums of elements of polar cones $C^{o}$ such that the cone $C$ belongs to the fan $\{H_{+},H_{-},H\}$, and the distance between $X$ and $C$ is at most $\delta$. In general (see [25]) for any cones $C_{1},C_{2}$ we have

[TABLE]

so it follows that we can rewrite (16) as

[TABLE]

The characterizations (16) and (17) have the advantage that they can be easily carried over to the more general case where $\mathcal{G}$ consists of several reversible edges. In that case we have several tug-of-wars going on at the same time, but for each one of them we can specify the winning direction (if any) at $x$ by calculating the distance between $X=\log x$ and some hyperplane in $\mathbb{R}^{n}$. Depending on whether $X$ falls outside an uncertainty region or not, each reversible edge $s\rightleftharpoons s^{\prime}$ of $G$ contributes one or two vectors to $F_{H,\delta}(X)$ (if one, then it is either $s^{\prime}-s$ or $s-s^{\prime}$, and if two, then they are $\pm(s^{\prime}-s)$).

It follows that any system (10) can be embedded into a differential inclusion on $\mathbb{R}^{n}_{+}$ given by a set $\mathcal{H}$ of hyperplanes in $\mathbb{R}^{n}$ and a number $\delta>0$, as follows. For each $x\in\mathbb{R}^{n}_{>0}$ we define $F_{\mathcal{H},\delta}(\log x)$ to be the convex cone generated by vectors orthogonal to the hyperplanes of $\mathcal{H}$, in the direction that goes from the point $X=\log x$ towards each hyperplane, and also the opposite direction if $X$ is at distance $<\delta$ from some hyperplane. If $X$ does not belong to any uncertainty region, then $F_{\mathcal{H},\delta}(\log x)$ is defined to be exactly the polar cone $C^{o}$ of the (unique) cone $C\in\mathcal{F}_{\mathcal{H}}$ that contains $X$. If $X$ does belong to some uncertainty regions, then we can still describe $F_{\mathcal{H},\delta}(\log x)$ in terms of polar cones, by including not just the polar of the cone of $\mathcal{F}_{\mathcal{H}}$ that contains $X$, but also the polar of each cone of $\mathcal{F}_{\mathcal{H}}$ that is at distance $\leq\delta$ from $X$.

Of course, not every fan is generated by a set of hyperplanes as above. Nevertheless, we can generalize the construction described above to define a differential inclusion given by a general fan $\mathcal{F}$ in $\mathbb{R}^{n}$, as follows.

Definition 3.3**.**

Consider a polyhedral fan $\mathcal{F}$ in $\mathbb{R}^{n}$, and a number $\delta>0$. The toric differential inclusion generated by $\mathcal{F}$ and $\delta$ is the differential inclusion on $\mathbb{R}^{n}_{>0}$ given by

[TABLE]

where $F_{\mathcal{F},\delta}$ is a set-valued function defined as

[TABLE]

In other words, the toric differential inclusion (18) is a piecewise constant differential inclusion, and its right-hand side $F_{\mathcal{F},\delta}(\log x)$ is the cone generated by the sum of all the polar cones $C^{o}$ such that $C\in\mathcal{F}$ and $dist(\log x,C)\leq\delta$. Note that for every $x$ there is at least one such $C$, because the fan $\mathcal{F}$ is complete. If $C\in\mathcal{F}$ is a cone of dimension $<n$ and $dist(\log x,C)\leq\delta$, then we say that $x$ belongs to the uncertainty region of $C$.

Like before, we can also rewrite $F_{\mathcal{F},\delta}$ as

[TABLE]

From Lemma 3.1 and the discussion above, it follows that any reversible polynomial dynamical system can be embedded into a toric differential inclusion:

Proposition 3.2**.**

Consider a variable-$k$ reversible polynomial dynamical system (10). Then this system can be embedded into a toric differential inclusion.

Given a reversible polynomial dynamical system generated by the E-graph $\mathcal{G}$, the most natural such embedding is obtained if we choose $\mathcal{F}$ to be the fan generated by the set $\mathcal{H}$ of hyperplanes that are orthogonal to the edge vectors of $\mathcal{G}$, i.e.,

[TABLE]

and we choose $\delta$ as suggested by Lemma 3.1, i.e., $\displaystyle\delta=\max_{s\rightleftharpoons s^{\prime}\in\mathcal{G}}\frac{2|\log\varepsilon|}{||s-s^{\prime}||}$.

On the other hand, we are most interested in weakly reversible polynomial dynamical systems, since the conjectures we described in Section 2 refer to this larger class of dynamical systems. In the next section we address this problem, and we prove that variable-$k$ weakly reversible polynomial dynamical systems can be embedded into toric differential inclusions. This will imply that toric dynamical systems [10] can be embedded into toric differential inclusions, and is some of the motivation for calling these differential inclusions “toric”.

4 Embedding of variable-$k$ weakly reversible polynomial dynamical systems into toric differential inclusions

As we discussed in the previous section, the simplest examples of toric differential inclusions are generated by polyhedral fans $\mathcal{F}_{\mathcal{H}}$ that are determined by a finite set $\mathcal{H}$ of hyperplanes. We will refer to this class of toric differential inclusions as hyperplane-generated toric differential inclusions.

We have also seen in the previous section that any variable-$k$ reversible polynomial dynamical system in $\mathbb{R}^{n}_{>0}$ can be embedded into a (hyperplane-generated) toric differential inclusion. Here we show that the same is true for all variable-$k$ weakly reversible polynomial dynamical systems.

Theorem 4.1**.**

Consider a variable-$k$ weakly reversible polynomial dynamical system (5). Then this system can be embedded into a toric differential inclusion.

Proof.

Denote by $\mathcal{G}$ a weakly reversible E-graph that generates our system. Consider first the case where $\mathcal{G}$ consists of a single oriented cycle. Then the graph $\mathcal{G}$ is given by $s_{1}\to s_{2}\to...\to s_{r}\to s_{1}$, and the variable-$k$ weakly reversible dynamical system it generates has the form

[TABLE]

where $s_{r+1}=s_{1}$ and $\varepsilon\leq k_{i}(t)\leq\frac{1}{\varepsilon}$ for some ${\varepsilon}>0$.

Consider the set $\mathcal{L}$ of lines through the origin in the direction of vectors $s_{i}-s_{j}$ for all $i\neq j$, and denote by $\mathcal{H}$ the set of all hyperplanes that are orthogonal to a line in $\mathcal{L}$, i.e.,

[TABLE]

Denote by $\mathcal{F}_{\mathcal{H}}$ the polyhedral fan generated by the set of hyperplanes $\mathcal{H}$, and, for $\delta>0$, denote by $\mathcal{T}_{\mathcal{H},\delta}$ the corresponding hyperplane-generated toric differential inclusion. We will show that there exists $\delta_{0}>0$ such that the single-cycle variable-$k$ weakly reversible dynamical system (21) is embedded in the toric differential inclusion $\mathcal{T}_{\mathcal{H},\delta_{0}}$. (Note that when refer below to reversible edges $s_{i}\rightleftharpoons s_{j}$, we do not assume that these edges belong to $\mathcal{G}$; we only mention them as a tool in our construction of $\mathcal{T}_{\mathcal{H},\delta_{0}}$.)

Choose $\delta_{0}>0$ large enough such that the uncertainty regions given by the reversible edges $s_{i}\rightleftharpoons s_{j}$ and $\varepsilon$ are contained within the uncertainty regions of the toric differential inclusion $\mathcal{T}_{\mathcal{H},\delta_{0}}$. For example, according to Lemma 3.1 (see also Proposition 3.2), we can choose

[TABLE]

Let us first consider the case of a point $x\in\mathbb{R}_{>0}^{n}$ that does not belong to any uncertainty region of $\mathcal{T}_{\mathcal{H},\delta_{0}}$. Then the point $X=\log x$ is not contained in any hyperplane in $\mathcal{H}$, so there must exist a cone $C$ in $\mathcal{F}_{\mathcal{H}}$ such that $C$ has dimension $n$ and contains the point $X$ in its interior. We will show that the right-hand side of (21) is contained in the polar cone $C^{o}$.

Consider a vector $w$ in the interior of $C$, and consider the orthogonal projections of the vectors $s_{1},s_{2},...,s_{r}$ on the line $\ell_{w}$ that passes through the origin in the direction given by $w$. Then no two such projections are the same, because $w$ does not belong to any of the hyperplanes in $\mathcal{H}$. In other words, we have that $s_{i}\cdot w\neq s_{j}\cdot w$ whenever $s_{i}\neq s_{j}$. We now use these projections to give a second set of names to the vectors $s_{1},s_{2},...,s_{r}$, say $v_{1},v_{2},...,v_{r}$, to record the order in which these projections appear along the line $\ell_{w}$. More precisely, we choose the names $v_{1},v_{2},...,v_{r}$ in the order (from largest to smallest) of the values of $v_{l}\cdot w$, so for example, $v_{1}$ equals the $s_{i}$ that has the largest value of $s_{i}\cdot w$, $v_{2}$ equals the $s_{i}$ that has the second-largest value of $s_{i}\cdot w$, and so on. Note also that, since the signs of the dot products $(v_{l+1}-v_{l})\cdot w$ cannot change as $w$ is allowed to vary in the interior of $C$, it follows that the new names $v_{1},v_{2},...,v_{r}$ do not depend on the particular choice of vector $w$ in the interior of $C$. In other words, the dot products $(v_{l+1}-v_{l})\cdot w$ are all negative numbers, for all $w$ in the interior of $C$. Therefore, the vectors $v_{2}-v_{1},v_{3}-v_{2},...,v_{r}-v_{r-1}$ belong to $C^{o}$.

So, in order to show that the right-hand side of (21) is included in $C^{o}$, it is enough to show that it can be written as a positive linear combination of the vectors $v_{2}-v_{1},v_{3}-v_{2},...,v_{r}-v_{r-1}$.

If $s_{i_{1}}\!=v_{1},\ s_{i_{2}}\!=v_{2},\ ...,\ s_{i_{r}}\!=v_{r}$, then note that $(i_{1},i_{2},...,i_{r})$ is a permutation of $(1,2,...,r)$. If we denote the inverse permutation by $(j_{1},j_{2},...,j_{r})$, it follows that $s_{1}=v_{j_{1}}$, $s_{2}=v_{j_{2}}$, and so on.

Then we have $s_{2}-s_{1}=v_{j_{2}}-v_{j_{1}}$. If $j_{2}>j_{1}$ we write

[TABLE]

and if $j_{2}<j_{1}$ we write

[TABLE]

We do the same for $s_{3}-s_{2},\ s_{4}-s_{3}$, and so on. This way, we write each difference $s_{i+1}-s_{i}$ from the right-hand side of (21) in terms of the vectors $\pm(v_{l+1}-v_{l})$, with $l=1,2,...,r-1$. Therefore we can re-group terms to obtain

[TABLE]

where $\Phi_{l}$ is a sum of several terms of the form $k_{i}\,x^{v_{i}}$, with various signs.

Note now that the positive terms inside $\Phi_{l}$ correspond to edges of the form $v_{m}\to v_{n}$ with $m\leq l<n$, and negative terms inside $\Phi_{l}$ correspond to edges of the form $v_{m}\to v_{n}$ with $n\leq l<m$. This means that the positive terms inside $\Phi_{l}$ are a sum of terms of the form $k_{i}(t)\,x^{v_{i}}$ with $i\leq l$, and the negative terms inside $\Phi_{l}$ are a sum of terms of the form $k_{i}(t)\,x^{v_{i}}$ with $i>l$. In particular, since $\mathcal{G}$ is a cycle, it follows that $\Phi_{l}$ contains at least one positive term, and at least one negative term. Recall that $X=\log x$ is in the interior of $C$. Then the dot products of the form $(v_{l+1}-v_{l})\cdot X$ must be negative for all $l$, which implies that $x^{v_{l}}>x^{v_{l+1}}$ for all $l$. Moreover, since $x$ does not belong to any uncertainty region, and due to our choice of $\delta_{0}$ (see (22) and Lemma 3.1), this inequality remains the same even if we include the terms $k_{l}(t)$ and $k_{l+1}(t)$, and we obtain $k_{l}(t)x^{v_{l}}>k_{l+1}(t)x^{v_{l+1}}$. Therefore, we have

[TABLE]

Also, note that the number of positive terms inside $\Phi_{l}$ is the same as the number of negative terms inside $\Phi_{l}$, because the graph $\mathcal{G}$ is a cycle. Therefore, the sum of the positive terms inside $\Phi_{l}$ dominates the sum of the negative terms inside $\Phi_{l}$, for each $l$. In conclusion, the right-hand side of (23) (and therefore the right-hand side of (21)) is a positive linear combination of the vectors $v_{l+1}-v_{l}$, for $1\leq l\leq r-1$, so it belongs to $C^{o}$.

Consider now the case where $x$ does belong to one or more uncertainty regions of the toric differential inclusion $\mathcal{T}_{\mathcal{H},\delta_{0}}$. Recall that the cone of $\mathcal{T}_{\mathcal{H},\delta_{0}}$ at $x$ is $F_{\mathcal{H},\delta_{0}}(X)$ (see Definition 3.3). Then, by using (19) and the calculations in the proof of Lemma 3.1, we conclude that the set $F_{\mathcal{H},\delta_{0}}(X)$ is a cone generated by two types of vectors: vectors of the “first type”, which are of the form $\pm(s_{i}-s_{j})$ such that $s_{i}\rightleftharpoons s_{j}$ is a reversible edge and $x$ is in the uncertainty region of the hyperplane orthogonal to $s_{i}-s_{j}$, and vectors of the “second type”, which are of the form $s_{l}-s_{m}$ such that $s_{l}\rightleftharpoons s_{m}$ is a reversible edge and $x$ is not in the uncertainty region of the hyperplane orthogonal to $s_{l}-s_{m}$, and moreover $(s_{l}-s_{m})\cdot X<0$. Note that, without loss of generality we can assume that not all vectors are of the first type; otherwise we immediately obtain that the right-hand side of (21) belongs to $F_{\mathcal{H},\delta_{0}}(X)$.

Consider now a vector $w$ in the interior of the polar cone of the cone $F_{\mathcal{H},\delta_{0}}(X)$; then $w$ satisfies $(s_{i}-s_{j})\cdot w=0$ for vectors of first type, and $(s_{l}-s_{m})\cdot w<0$ for vectors of second type. We also consider the orthogonal projections of the vectors $s_{1},s_{2},...,s_{r}$ on the line $\ell_{w}$ that passes through the origin in the direction given by $w$. Unlike the previous case, in this case some projections will coincide; more precisely, if $s_{i}$ and $s_{j}$ are like in the first type above, then their projections will coincide because $s_{i}\cdot w=s_{j}\cdot w$.

Nevertheless, we can still give a second set of names to the vectors $s_{1},s_{2},...,s_{r}$, say $v_{1},v_{2},...,v_{r}$, to record the order in which these projections appear along the line $\ell_{w}$, with the caveat that we will have one or more cases where the projections coincide. Our ordering is chosen such that

[TABLE]

Note that if we allow $w$ to vary within the interior of the polar cone of the cone $F_{\mathcal{H},\delta_{0}}(X)$, the inequality (24) will still hold. This implies that whenever the projections on $\ell_{w}$ of $v_{l}$ and $v_{l+1}$ are distinct, the vector $v_{l+1}-v_{l}$ belongs to the cone $F_{\mathcal{H},\delta_{0}}(X)$; moreover, if their projections coincide, then both vectors $\pm(v_{l+1}-v_{l})$ belong to $F_{\mathcal{H},\delta_{0}}(X)$.

Now we proceed exactly like in the previous case. We can still re-group terms like in formula (23), and, in order to conclude that the right-hand side of (23) belongs to $F_{\mathcal{H},\delta_{0}}(X)$, we only have to check that $\Phi_{l}>0$ for values of $l$ where projections of $v_{l}$ and $v_{l+1}$ are distinct. But, in the same way as before, such $\Phi_{l}$ have an equal number of positive and negative terms of the form $k_{i}x^{v_{i}}$, and, also as before, the positive terms are larger than the negative terms.

Finally, if the weakly reversible graph $\mathcal{G}$ is not a single oriented cycle, then we write it as a union of cyclic graphs, $\displaystyle\mathcal{G}=\bigcup_{i=1}^{g}\mathcal{G}_{i}$, and we can argue as above for each such $\mathcal{G}_{i}$. We obtain that the variable-$k$ weakly reversible dynamical systems given by the cycle $\mathcal{G}_{i}$ are embedded in toric differential inclusions generated by some set of hyperplanes $\mathcal{H}_{i}$. Note now that the right-hand side of a variable-$k$ weakly reversible dynamical system given by $\mathcal{G}$ can be decomposed into a sum of terms, such that each term is of the form given by the right-hand side of a variable-$k$ weakly reversible dynamical system determined by $\mathcal{G}_{i}$. (We may have to use smaller $\varepsilon_{i}$ values for the terms in the decomposition, because the same edge of $\mathcal{G}$ may belong to several graphs $\mathcal{G}_{i}$.) Note also that if we define the set of hyperplanes

[TABLE]

then, for any $i$ and $\delta$, the toric differential inclusion given by the fan $\mathcal{F}_{\mathcal{H}_{i}}$ and $\delta$ is embedded in the toric differential inclusion given by the fan $\mathcal{F}_{\mathcal{H}}$ and $\delta$, because every cone $C_{i}\in\mathcal{F}_{\mathcal{H}_{i}}$ can be written as a union of cones from $\mathcal{F}_{\mathcal{H}}$, and whenever $C\subset\tilde{C}$ it follows that $C^{o}\supset\tilde{C}^{o}$. Then we conclude that any variable-$k$ mass-action system given by $\mathcal{G}$ can be embedded into a toric differential inclusion generated by the set of hyperplanes $\displaystyle\mathcal{H}$. ∎

5 Applications and conclusions

In this paper we have introduced toric differential inclusions, and we have shown that any polynomial dynamical system on the positive orthant is generated by an E-graph (which is not unique). Moreover, if this E-graph can be chosen to be weakly reversible, then the polynomial dynamical system can be embedded into a toric differential inclusion. Most importantly, toric differential inclusions have a rich geometric structure that can be used in the construction of invariant regions needed for the proof of important conjectures in this field [8, 9].

The idea of thinking about polynomial dynamical systems as being generated by E-graphs was inspired by the way reaction networks generate polynomial dynamical systems under the assumption of mass-action kinetics. Indeed, the set of polynomial dynamical systems generated by E-graphs $\mathcal{G}=(V,E)$ that satisfy $V\subset\mathbb{Z}_{\geq 0}^{n}$ is exactly the same as the set of all mass-action dynamical systems [3, 4].

Note also that while the formulations of the conjectures in Section 2.2 are more general than the usual formulations for mass-action systems (which restrict the exponent vectors to be nonnegative), there is a simple way to show that these two versions are actually equivalent, by “time-rescaling” via multiplication by a scalar field of the form $x^{(M\cdot{\mathbf{1})}}$, where $\mathbf{1}$ is a vector with all 1 components. This multiplication can be chosen such that it shifts the E-graph into the non-negative orthant, while preserving all trajectory curves [8].

On the other hand, our results on embedding weakly reversible polynomial dynamical systems into toric differential inclusions suggests other kinds of generalizations of the Persistence Conjecture and of the Permanence Conjecture, as follows. Note that, in the proof of Lemma 3.1, in order to be able to obtain that embedding into a differential inclusion, it is not really necessary to know that the values of $k_{s\to s^{\prime}}(t),\ k_{s^{\prime}\to s}(t)$ are contained in an interval of the form $[\varepsilon,\frac{1}{\varepsilon}]$; exactly the same calculations will work if we just know that

[TABLE]

for some $\varepsilon_{0}>0$. Similarly, in Proposition 3.2, we don’t need to assume that all the time-dependent parameter values $k_{s\to s^{\prime}}(t)$ are bounded away from zero and infinity; we just need to assume that for any reversible reaction the ratio of the corresponding parameter values is bounded away from zero and infinity, as in (25).

Even in Theorem 4.1, in the case of a weakly reversible E-graph $\mathcal{G}$, it is enough to assume that we have

[TABLE]

for any edges $s\to s$ and $\tilde{s}\to\tilde{s}^{\prime}$ that are in the same connected component of $\mathcal{G}$.

These observations suggest that, for example, if we could take advantage of the embedding into toric differential inclusions and prove this stronger version of the Persistence Conjecture, we would obtain the following interesting conclusion for the dynamics of chemical and biochemical reaction networks: if several (weakly reversible) reaction networks or pathways are coupled together, then the resulting dynamics will still be persistent even if the external factors that influence each pathway are widely different in size. For example, if we have two weakly reversible biochemical pathways, and the reaction rate parameters in one pathway are modulated by temperature, while in the other pathway they are modulated by a signaling protein whose concentration is unrelated to temperature, then by coupling together these two pathways we should still maintain the persistence property.

Also, generally speaking, the embedding of weakly reversible mass-action systems into toric differential inclusions provides a more rigorous interpretation for the standard “hand-waving” intuition behind the Persistence and Permanence Conjectures. Namely, for mass-action systems it is quite reasonable to think that, if every reaction is part of a cycle (i.e., if the reaction network is weakly reversible) then the chemical reactions should not be able to drive any concentration to zero, because if a reaction consumes a species, then soon a chain of reactions that follow it will work to replenish that species. This intuition is a bit too vague to be turned into a proof, but, as we can see in Fig. 2$(h)$, we now have a more concrete way to think about it: if we focus on some region where the right-hand side of the toric differential inclusion is constant, then this constant cone of directions seems to point “toward the middle” of the positive quadrant, i.e., away from the boundary, and also away from infinity.

In particular, given an embedding of a two-variable polynomial dynamical system into a toric differential inclusion, we can immediately obtain families of invariant regions for this dynamical system in $\mathbb{R}^{2}_{>0}$, as illustrated in Fig. 3. Such an embedding allows us to construct “zero-separating curves” (as shown in Fig. 3$(c)$) which prevent positive trajectories from approaching the origin. As shown in [9], when constructed for variable-$k$ polynomial dynamical systems, these curves are the key tool for a proof of the persistence of vertex-balanced dynamical systems in $\mathbb{R}^{3}_{>0}$. Similarly, given any positive initial point, we can use the embedding to construct a compact invariant region that contains that point (as shown in Fig. 3$(d)$).

Therefore, the embeddings of some polynomial dynamical systems into toric differential inclusions allow us to give very short proofs of the main results in [9] (i.e., a proof of the Persistence Conjecture in $\mathbb{R}^{2}_{>0}$, and of the Global Attractor Conjecture in $\mathbb{R}^{3}_{>0}$) and also to generalize the Persistence Conjecture result in $\mathbb{R}^{2}_{>0}$ as we described above; see [8] for more details. More importantly, similar geometric constructions of invariant regions based on toric differential inclusions can also be done in higher dimensions [8].

Outside the setting considered here, global convergence results for mass-action systems have been recently used to study reaction-diffusion equations via the methods of lines [23], have been adapted for analyzing some versions of discrete Boltzmann equations [12], and an entropy method has been used to study a large class of reaction-diffusion systems that arise from vertex-balanced networks [13]. Interestingly, reversibility and weak reversibility play a role in these works as well.

6 Acknowledgments

This work benefited from feedback from the organizers and participants of the Workshop on the Global Attractor Conjecture held at San Jose State University in March 2016. Exceptionally useful and detailed reviewer comments helped improve the presentation of these results. We also acknowledge support from the National Science Foundation under grants DMS-1412643 and DMS-1816238.