Polynomial Dynamical Systems, Reaction Networks, and Toric Differential Inclusions

# Polynomial Dynamical Systems, Reaction Networks, and Toric Differential Inclusions

Gheorghe Craciun
Department of Mathematics and
Department of Biomolecular Chemistry
e-mail: craciun@math.wisc.edu

July 20, 2019
###### 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 can be regarded as being generated by an oriented graph embedded in , called 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 toric differential inclusions, which are piecewise constant autonomous dynamical systems with a remarkable geometric structure. We show that if a Euclidean embedded graph has some reversibility properties, then any polynomial dynamical system generated by 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.

## 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 can be represented as

 dxdt=m∑i=1xsivi (1)

where , are some vectors in called exponent vectors, denotes the monomial , and are vectors in . A solution of (1) is a function that satisfies (1), where is an interval in .

Note that, since the coordinates are positive, the monomials are well-defined even if the coordinates of the exponent vectors are arbitrary real numbers (i.e., are not necessarily in ). 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

 dxdt=m∑i=1kixsivi (2)

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

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

 dxdt=m∑i=1ki(t)xsivi (3)

where the functions are positive and uniformly bounded, i.e., there exists some such that for all . We will refer to models of the form (3) as variable- 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 (i.e., no variable “goes extinct”), irrespective of the choices of uniformly bounded external factors ?

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-) 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 whose vertices are labeled by distinct elements of for some . With an abuse of notation, we identify the set with the set of vertex labels, i.e., we assume that . Moreover, we associate to each edge its edge vector . Also, we define its source vertex to be , and its target vertex to be .

Given an Euclidean embedded graph , the polynomial dynamical systems generated by are the dynamical systems on given by

 dxdt=∑e∈Ekexs(e)v(e) (4)

for some positive constants . Note that if then (4) is just mass-action kinetics for a chemical reaction network represented by , i.e, one where there is a reaction of the form for each edge . (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 is proportional to ; see [9, 16, 20, 29] for more details.)

More generally, the variable-k polynomial dynamical systems generated by are the (nonautonomous) dynamical systems on given by

 dxdt=∑e∈Eke(t)xs(e)v(e) (5)

such that there exists some for which we have for all and for all .

Let us note that for any variable- polynomial dynamical system (3) we can construct an E-graph that generates it, and is not unique. Assuming that the ordered pairs are distinct, the simplest way to construct such a is to choose the set of vertices

 V={si | i=1,...,m}∪{si+vi | i=1,...,m},

and the set of edges

 E={(si,si+vi) | i=1,...,m}.

If we want to obtain a different E-graph that generates (3), we can, for example, write one of the vectors as a positive linear combination of two different nonzero vectors, and use these new vectors to obtain a graph with 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 is reversible if for any edge in the reverse edge also belongs to . Also we say that is weakly reversible if any edge is part of an oriented cycle in , or equivalently, any connected component of 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 is reversible, and the graphs in Fig. 1 and 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- reversible and weakly reversible polynomial dynamical systems.

Example 1. Consider the dynamical system given by

 dx1dt = −2k1(t)x21+2k2(t)x2 (6) dx2dt = k1(t)x21−k2(t)x2

for some functions with for all .

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

 dxdt=k1(t)x21(−21)+k2(t)x2(2−1), (7)

where . In turn, this can be written in the form (3), as follows:

 dxdt=k1(t)xs1(−21)+k2(t)xs2(2−1), (8)

where and . Then, the simplest E-graph that generates the dynamical system (6) has two edges, one edge going from to , and the other edge going from to . But, note that we happen to have and , so the graph actually has only two vertices, and is reversible. The graph is shown in Fig. 2 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- polynomial dynamical system in is persistent if, for any solution defined on an interval that contains , there exists some such that

 xi(t)>ε0  for all i∈{1,...,n} and for all t∈I∩[0,∞).

In other words, the system is called persistent if, for any solution with positive initial condition, there exists a positive lower bound for all the variables , and for all the future times at which the solution is defined (we cannot say that there exists a positive lower bound for all 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-) polynomial dynamical system is given by (4) or (5), then its edge space is the linear span of the set of edge vectors . Then we define its affine invariant sets to be the sets of the form

 (x0+S)∩Rn>0,  for some x0∈Rn>0.

These are indeed invariant spaces for solutions of (4) or (5) on the domain , because all the vectors that appear on the right-hand side of these equations are linear combinations of , so are contained in .

We say that a (variable-) polynomial dynamical system on is permanent if, for each affine invariant set there exists a compact set such that any solution with can be extended for all , and there exists some such that for all . It follows that, for a permanent system, all solutions that start in can be extended for all , and, for large enough , 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 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 if there exists a point such that any solution with can be extended for all and we have

 limt→∞x(t)=¯x0.

Also, we say that a polynomial dynamical system (2) is vertex balanced if there exists an E-graph that generates our system as in (4), and there exists a point such that for any vertex we have

 ∑e=(s,s′)∈Eke¯xs=∑e=(s′,s)∈Eke¯xs′. (9)

In other words, if we think of the positive number as the rate of a flow along the edge (i.e., a flow from the vertex to the vertex ), then condition (9) says that, if , then at each vertex of the graph , 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- weakly reversible polynomial dynamical system is persistent.

Extended Permanence Conjecture. Any variable- 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 has a single connected component [2], and this result has been generalized by Gopalkrishnan, Miller, and Shiu for the case where the graph 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 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- 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 ), but are also piecewise constant in .

## 3 Toric differential inclusions

Given an E-graph , let us write if is an edge of ; also let us write if both and are edges of .

Then, if is reversible, the dynamical system (5) can be written as

 dxdt=∑s⇌s′∈E(ks→s′(t)xs−ks′→s(t)xs′)(s′−s), (10)

by grouping together pairs of terms given by an edge and its reverse . In particular, if consists of a single reversible edge , then we obtain

 dxdt=(ks→s′(t)xs−ks′→s(t)xs′)(s′−s). (11)

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 of the system along the same line (parallel to the vector ), but in opposite directions. Recall that for all . Then, the domain can be partitioned into three regions: the region where the inequality holds (which implies ), the region where the inequality holds (which implies ), and an uncertainty region where neither one of these two inequalities are satisfied, and either one of the two terms and may win the tug-of-war, or there can be a tie, due to the fact that and may take any values between and .

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

Indeed, let us define by the line through the origin in generated by the vector , and by the ray starting from the origin in the direction , and by the ray starting from the origin in the direction . Also, let us denote by the hyperplane through the origin orthogonal to . For some define the set-valued function at , as follows:

 FH,δ(X)=⎧⎪⎨⎪⎩l+s′−s if \ \ dist(X,H)>δ and X⋅(s′−s)>0l−s′−s if \ \ dist(X,H)>δ and X⋅(s′−s)<0ls′−s if \ \ dist(X,H)≤δ

Recall that for all . We have:

###### Lemma 3.1.

The dynamical system (11) (which is given by an E-graph that consists of a single reversible edge ) is embedded in the differential inclusion

 dxdt∈FH,δ(logx), (12)

where .

###### Proof.

If there is nothing to be proved, because we already know that the right-hand side of (11), i.e., the vector

 (ks→s′(t)xs−ks′→s(t)xs′)(s′−s)

belongs to .

If and , we need to show that the right-hand side of (11) belongs to , i.e., we need to show that

 ks→s′(t)xs

For this, it is sufficient to show that , which is equivalent to , and, by taking logarithm on both sides of this inequality, it can also be written as

 (logx)⋅(s′−s)>−2logε. (13)

On the other hand, is just the dot product between the vector and the unit vector that is orthogonal to and on the same side of as . But, since , this unit vector is just . Then the inequality implies

 (logx)⋅s′−s||s′−s||>δ. (14)

Finally, given that and , we can see that the inequalities (13) and (14) are equivalent.

The case where and is analogous, so this concludes the proof. ∎

Therefore, if consists of a single reversible edge and if , then the right-hand side of the system (11) is a vector that is orthogonal to the hyperplane and points in the direction that goes from the point towards (see Fig. 2- for some examples in ). 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 is the set of nonnegative linear combinations of a finite set of vectors in , or, equivalently, is a finite intersection of half-spaces in [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 . The polar cone of is denoted and is given by

 Co={y∈Rn| x⋅y≤0 for all x∈C}. (15)

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

Let us now return to the graph that consists of a single reversible edge . Recall that denotes the hyperplane through the origin and orthogonal to . Denote by the closed half-space of that is bounded by and contains the vector , and denote by the closed half-space of that is bounded by and contains the vector . Note that the polar cone is equal to the ray , and the polar cone is equal to the ray . Then the set-valued function can be rewritten as

 FH,δ(X)=⎧⎪⎨⎪⎩Ho+ if \ \ dist(X,H)>δ and X∈H+Ho− if \ \ dist(X,H)>δ and X∈H−Ho++Ho− if \ \ dist(X,H)≤δ,

where we define .

Then we can write only in terms of the distance between and the half-spaces and , as follows:

 FH,δ(X)=⎧⎪⎨⎪⎩Ho+ if \ \ dist(X,H+)≤δ and dist(X,H−)>δHo− if \ \ dist(X,H−)≤δ and dist(X,H+)>δHo++Ho− if \ \ dist(X,H+)≤δ and dist(X,H−)≤δ

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

###### Definition 3.2.

A finite set of polyhedral cones in is a polyhedral fan if the following two conditions are satisfied:

any face of a cone in is also in ,

the intersection of two cones in is a face of both cones.

We say that a polyhedral fan is complete if  .

For example, given a finite set of hyperplanes that contain the origin, consider the set of all the possible intersections of half-spaces given by hyperplanes in . Then 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 as

 FH,δ(X)=∑C∈{H+,H−,H}dist(X,C)≤δCo (16)

i.e., consists of all possible sums of elements of polar cones such that the cone belongs to the fan , and the distance between and is at most . In general (see [25]) for any cones we have

 Co1+Co2=(C1∩C2)o

so it follows that we can rewrite (16) as

 FH,δ(X)=⎛⎜ ⎜⎝⋂C∈{H+,H−,H}dist(X,C)≤δC⎞⎟ ⎟⎠o (17)

The characterizations (16) and (17) have the advantage that they can be easily carried over to the more general case where 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 by calculating the distance between and some hyperplane in . Depending on whether falls outside an uncertainty region or not, each reversible edge of contributes one or two vectors to (if one, then it is either or , and if two, then they are ).

It follows that any system (10) can be embedded into a differential inclusion on given by a set of hyperplanes in and a number , as follows. For each we define to be the convex cone generated by vectors orthogonal to the hyperplanes of , in the direction that goes from the point towards each hyperplane, and also the opposite direction if is at distance from some hyperplane. If does not belong to any uncertainty region, then is defined to be exactly the polar cone of the (unique) cone that contains . If does belong to some uncertainty regions, then we can still describe in terms of polar cones, by including not just the polar of the cone of that contains , but also the polar of each cone of that is at distance from .

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 in , as follows.

###### Definition 3.3.

Consider a polyhedral fan in , and a number . The toric differential inclusion generated by and is the differential inclusion on given by

 dxdt∈FF,δ(logx), (18)

where is a set-valued function defined as

 FF,δ(X)=∑C∈Fdist(X,C)≤δCo. (19)

In other words, the toric differential inclusion (18) is a piecewise constant differential inclusion, and its right-hand side is the cone generated by the sum of all the polar cones such that and . Note that for every there is at least one such , because the fan is complete. If is a cone of dimension and , then we say that belongs to the uncertainty region of .

Like before, we can also rewrite as

 FF,δ(X)=⎛⎜ ⎜⎝⋂C∈Fdist(X,C)≤δC⎞⎟ ⎟⎠o. (20)

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- 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 , the most natural such embedding is obtained if we choose to be the fan generated by the set of hyperplanes that are orthogonal to the edge vectors of , i.e.,

 H={(s′−s)⊥| s⇌s′∈G},

and we choose as suggested by Lemma 3.1, i.e., .

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- 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 that are determined by a finite set 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- reversible polynomial dynamical system in can be embedded into a (hyperplane-generated) toric differential inclusion. Here we show that the same is true for all variable- weakly reversible polynomial dynamical systems.

###### Theorem 4.1.

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

###### Proof.

Denote by a weakly reversible E-graph that generates our system. Consider first the case where consists of a single oriented cycle. Then the graph is given by , and the variable- weakly reversible dynamical system it generates has the form

 dxdt=r∑i=1ki(t)xsi(si+1−si), (21)

where and for some .

Consider the set of lines through the origin in the direction of vectors for all , and denote by the set of all hyperplanes that are orthogonal to a line in , i.e.,

 H={ℓ⊥|ℓ∈L}.

Denote by the polyhedral fan generated by the set of hyperplanes , and, for , denote by the corresponding hyperplane-generated toric differential inclusion. We will show that there exists such that the single-cycle variable- weakly reversible dynamical system (21) is embedded in the toric differential inclusion . (Note that when refer below to reversible edges , we do not assume that these edges belong to ; we only mention them as a tool in our construction of .)

Choose large enough such that the uncertainty regions given by the reversible edges and are contained within the uncertainty regions of the toric differential inclusion . For example, according to Lemma 3.1 (see also Proposition 3.2), we can choose

 δ0=maxi≠j2|logε|||si−sj||. (22)

Let us first consider the case of a point that does not belong to any uncertainty region of . Then the point is not contained in any hyperplane in , so there must exist a cone in such that has dimension and contains the point in its interior. We will show that the right-hand side of (21) is contained in the polar cone .

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

So, in order to show that the right-hand side of (21) is included in , it is enough to show that it can be written as a positive linear combination of the vectors .

If , then note that is a permutation of . If we denote the inverse permutation by , it follows that , , and so on.

Then we have . If we write

 s2−s1=j2−1∑l=j1(vl+1−vl),

and if we write

 s2−s1=−j1−1∑l=j2(vl+1−vl),

We do the same for , and so on. This way, we write each difference from the right-hand side of (21) in terms of the vectors , with . Therefore we can re-group terms to obtain

 dxdt=r−1∑l=1Φl(vl+1−vl), (23)

where is a sum of several terms of the form , with various signs.

Note now that the positive terms inside correspond to edges of the form with , and negative terms inside correspond to edges of the form with . This means that the positive terms inside are a sum of terms of the form with , and the negative terms inside are a sum of terms of the form with . In particular, since is a cycle, it follows that contains at least one positive term, and at least one negative term. Recall that is in the interior of . Then the dot products of the form must be negative for all , which implies that for all . Moreover, since does not belong to any uncertainty region, and due to our choice of  (see (22) and Lemma 3.1), this inequality remains the same even if we include the terms and , and we obtain . Therefore, we have

 k1(t)xv1>k2(t)xv2>...>kr(t)xvr.

Also, note that the number of positive terms inside is the same as the number of negative terms inside , because the graph is a cycle. Therefore, the sum of the positive terms inside dominates the sum of the negative terms inside , for each . In conclusion, the right-hand side of (23) (and therefore the right-hand side of (21)) is a positive linear combination of the vectors , for , so it belongs to .

Consider now the case where does belong to one or more uncertainty regions of the toric differential inclusion . Recall that the cone of at is (see Definition 3.3). Then, by using (19) and the calculations in the proof of Lemma 3.1, we conclude that the set is a cone generated by two types of vectors: vectors of the “first type”, which are of the form such that is a reversible edge and is in the uncertainty region of the hyperplane orthogonal to , and vectors of the “second type”, which are of the form such that is a reversible edge and is not in the uncertainty region of the hyperplane orthogonal to , and moreover . 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 .

Consider now a vector in the interior of the polar cone of the cone ; then satisfies for vectors of first type, and for vectors of second type. We also consider the orthogonal projections of the vectors on the line that passes through the origin in the direction given by . Unlike the previous case, in this case some projections will coincide; more precisely, if and are like in the first type above, then their projections will coincide because .

Nevertheless, we can still give a second set of names to the vectors , say , to record the order in which these projections appear along the line , with the caveat that we will have one or more cases where the projections coincide. Our ordering is chosen such that

 (vl+1−vl)⋅w≤0. (24)

Note that if we allow to vary within the interior of the polar cone of the cone , the inequality (24) will still hold. This implies that whenever the projections on of and are distinct, the vector belongs to the cone ; moreover, if their projections coincide, then both vectors belong to .

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 , we only have to check that for values of where projections of and are distinct. But, in the same way as before, such have an equal number of positive and negative terms of the form , and, also as before, the positive terms are larger than the negative terms.

Finally, if the weakly reversible graph is not a single oriented cycle, then we write it as a union of cyclic graphs, , and we can argue as above for each such . We obtain that the variable- weakly reversible dynamical systems given by the cycle are embedded in toric differential inclusions generated by some set of hyperplanes . Note now that the right-hand side of a variable- weakly reversible dynamical system given by 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- weakly reversible dynamical system determined by . (We may have to use smaller values for the terms in the decomposition, because the same edge of may belong to several graphs .) Note also that if we define the set of hyperplanes

 H=g⋃i=1Hi

then, for any and , the toric differential inclusion given by the fan and is embedded in the toric differential inclusion given by the fan and , because every cone can be written as a union of cones from , and whenever it follows that . Then we conclude that any variable- mass-action system given by can be embedded into a toric differential inclusion generated by the set of hyperplanes . ∎

## 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 that satisfy 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 , where 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 are contained in an interval of the form ; exactly the same calculations will work if we just know that

 ε0≤ks→s′(t)ks′→s(t)≤1ε0 (25)

for some . Similarly, in Proposition 3.2, we don’t need to assume that all the time-dependent parameter values 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 , it is enough to assume that we have

 ε