The problem and main result

A finite reaction network with \(d\ge1\) species consists of a finite set of complexes \(\mathcal C\subset\mathbb Z_{\ge0}^d\) and a set \(\mathcal R\subset\{(y,y'):y,y'\in\mathcal C,\ y\ne y'\}\) of directed reactions. We write a reaction as \(y\to y'\). The network is weakly reversible if, for each \(y\to y'\in\mathcal R\), there is a directed path from \(y'\) to \(y\). A connected component of the underlying undirected reaction graph is called a linkage class. There is no restriction on their number. Empty reaction sets and isolated complexes are allowed.

Assign a positive constant \(\kappa_{y\to y'}\) to each reaction. The mass-action equation is \[ \dot x=f(x):= \sum_{y\to y'\in\mathcal R} \kappa_{y\to y'}x^y(y'-y), \qquad x^y:=\prod_{i=1}^d x_i^{y_i}. \tag{1}\] Since the exponents are nonnegative integers, \(f\) is a polynomial vector field. A forward trajectory in \(\mathbb R_{>0}^d\) is called persistent if \[\liminf_{t\to\infty}x_i(t)>0\qquad(1\le i\le d).\] The bounds in this definition may depend on the trajectory.

In particular, every positive trajectory is bounded and persistent. To relate the statement to conservation laws, put \[S=\operatorname{span}_{\mathbb R}\{y'-y:y\to y'\in\mathcal R\}, \qquad P=(x^0+S)\cap\mathbb R_{\ge0}^d.\] The subspace \(S\) is the stoichiometric subspace, and \(P\) is the closed stoichiometric compatibility class of \(x^0\). Since \(f(x)\) takes values in \(S\), integration of the equation shows that a positive trajectory starting at \(x^0\) stays in \(P\). Theorem 1 applies whether or not \(P\) is bounded. It therefore includes the persistence question for bounded entire compatibility classes. The estimate (2) also shows directly that the omega-limit set is nonempty and disjoint from the orthant boundary: every sequence of times tending to infinity has a subsequence along which the trajectory converges inside the compact positive box in (2). We make no assertion in this generality that the trajectory converges to an equilibrium, or that the same \(\varepsilon\) works for all positive initial points. In particular, the theorem does not assert uniform eventual bounds over all positive initial points in a compatibility class.

The proof gives a geometric strengthening, stated precisely in Proposition 6: each \(x^0>0\) belongs to a compact convex polytope \(K\subset\mathbb R_{>0}^d\) such that every solution starting in \(K\) exists globally and remains in \(K\). The set \(K\) need not be a whole compatibility class. It is constructed in the ambient species space and may meet several such classes. Its convexity also yields a short proof of the known existence theorem of Boros (Boros 2019): every positive stoichiometric class contains a positive equilibrium (Corollary 7).

History and the role of the construction

The classical structural theory begins with complex balance. A system is complex-balanced if it admits a positive state at which the total incoming and outgoing reaction fluxes agree at every complex. Horn and Jackson’s equilibrium and stability theory (Horn and Jackson 1972), Horn’s complex-balancing criteria (Horn 1972), and Feinberg’s deficiency-zero analysis (Feinberg 1972) provide the foundations. Weak reversibility is a property of the graph, whereas complex balance also depends on the rates. The global-attractor question for complex-balanced systems asks whether every positive trajectory converges to the unique positive equilibrium in its stoichiometric class. This convergence follows from Theorem 1 under the additional complex-balance assumption; Remark 8 gives the argument.

Feinberg (Feinberg 1987, Remark 6.1.E) formulated the weaker expectation that a positive trajectory of a weakly reversible mass-action system cannot converge to a point where a species concentration is zero. Persistence excludes approach to the boundary along any sequence of times. Anderson (Anderson 2011b, sec. 1.1) separated the boundedness conjecture, that positive trajectories of weakly reversible mass-action systems are bounded, from the persistence conjecture for bounded trajectories. Together these imply persistence without an assumed upper bound. Theorem 1 establishes both conclusions for arbitrary positive constant reaction rates.

August and Barahona (August and Barahona 2010, Theorem 6) previously stated the general boundedness and persistence conclusion. Their proof, Part III on page 45 before Equation (11), infers monomial dominance from a maximal sum of exponents and a sufficiently large product of concentrations. That implication does not hold throughout the stated region. In \(2A\rightleftarrows B\) with both rates one, the concentrations \((a,b)=(T,T^3)\) satisfy \(ab\to\infty\) but \(a^2<b\), and \[\dot a+\dot b=b-a^2=T^3-T^2>0\qquad(T>1).\] This contradicts the asserted implication. It does not contradict their theorem: \(a+2b\) is conserved in this example. Our proof controls the relevant exponent pairings separately and establishes the conclusion independently.

Several earlier results give classwise permanence or permit time-dependent rates. Here permanence means that, for each fixed compatibility class, there is a compact subset of the positive orthant that eventually contains every positive trajectory in that class, with entry time depending on the initial state. Anderson (Anderson 2011b) proved boundedness for weakly reversible systems with a single linkage class. Gopalkrishnan, Miller, and Shiu (Gopalkrishnan et al. 2014) proved permanence for the larger class of strongly endotactic systems, which includes every weakly reversible single-linkage system; Boros and Hofbauer (Boros and Hofbauer 2020) gave a later proof of the single-linkage result. Craciun, Nazarov, and Pantea (Craciun et al. 2013) proved boundedness, persistence, and permanence for two-species weakly reversible systems as part of their more general endotactic result. Pantea (Pantea 2012) proved persistence of bounded weakly reversible trajectories when the stoichiometric subspace has dimension two, with arbitrary species and linkage counts. Two species and two-dimensional stoichiometry are different restrictions. These results also permit rates that vary in time while remaining bounded above and away from zero. Our theorem treats constant rates in all dimensions and with any number of linkage classes. Its bounds may depend on the initial point, so it does not assert the classwise uniformity of permanence.

The revised toric-differential-inclusion approach of Craciun (Craciun 2026, sec. 8, final remark, p. 81), in version 3 dated September 23, 2026, gives persistence for bounded positive trajectories of weakly reversible systems with uniformly positive and bounded variable rates. Its complex-balanced application combines the boundary construction with the Horn–Jackson Lyapunov bound. The general weakly reversible statement assumes a trajectory upper bound. Our construction supplies that bound together with persistence, through a compact convex forward-invariant polytope containing each prescribed positive initial point.

The geometric and analytical methods also have substantial antecedents. Invariant polygons and quantitative monomial comparisons appear in Craciun, Nazarov, and Pantea (Craciun et al. 2013), and projected polyhedra in Pantea (Pantea 2012). Anderson’s monomial tiers and return-path domination (Anderson 2011a), the projection and multiscale arguments of Gopalkrishnan, Miller, and Shiu (Gopalkrishnan et al. 2013, 2014), and Craciun’s ordered-complex embedding into toric differential inclusions (Craciun 2019) organize related comparisons. Our finite affine-label construction and directed-cut estimate are proved directly; none of these research results is an input to the proof.

Proof overview

The obstacle is to find a single finite family of linear inequalities that controls both small and large concentrations. An individual reaction need not point inward across a proposed face. Weak reversibility supplies return paths, so the inward estimate will concern the combined flux across cuts of the reaction graph. The normal of each active face must make strict comparisons that agree with the monomial sizes. Near a tie, an arbitrary approximation to a normal can reverse a comparison; the geometric construction must respect these ties.

For \(0<h<1\) and \(p\in[-1,1]^d\), write \[X_h(p)=(h^{p_1},\ldots,h^{p_d}).\] These points fill the box \(B_h=[h,h^{-1}]^d\). We construct a concave piecewise-affine function \[F_h(x)=\min_{\lambda\in\Lambda} \bigl(r_\lambda\cdot x+c_\lambda(h)\bigr),\] where \(\Lambda\) is a finite set indexing the affine functions, and the slopes \(r_\lambda\in[-1,1]^d\) are independent of \(h\). A label \(\lambda\) is active at \(x\) if its affine function attains this minimum; every minimizer is active at a tie.

Lemma 2 and Corollary 3 supply the family for any positive tolerance \(E:[-1,1]^d\to(0,\infty)\): for all sufficiently small \(h\), every label active at \(X_h(p)\) satisfies \(\lVert p-r_\lambda\rVert_\infty<E(r_\lambda)\), uniformly over the cube. No continuity of \(E\) is required. We can therefore make the tolerance small near the hyperplanes where a coordinate or a complex comparison vanishes, while keeping it positive on those hyperplanes. This forces each nonzero slope comparison to have the same sign as the corresponding exponent comparison, with a fixed positive gap. In particular, for some \(H>0\) independent of \(h\), \[r_\lambda\cdot y<r_\lambda\cdot z \quad\Longrightarrow\quad x^z\le h^H x^y \qquad(x=X_h(p),\ \lambda\text{ active}).\] The parameter order is essential: first the tolerance, then the fixed finite slopes, then their positive separation \(H\), and only then the final value of \(h\). A finite covering by nearby slopes alone would not ensure that the slopes attaining the minimum have the required accuracy.

Order the complexes by an active slope. Across a cut between two successive distinct levels, weak reversibility supplies an upward reaction whenever there is a downward reaction. The monomial comparison makes one upward rate dominate the sum of all downward rates once \(h\) is small enough. Summing over the cuts gives \(r_\lambda\cdot f(x)\ge0\) for every active label. Consequently \(F_h\) is nondecreasing along any trajectory segment in \(B_h\). At a tie, its right derivative is the minimum of the derivatives of all active branches, so controlling a preferred active label would not suffice. Equal slope levels need not have equal monomials; edges within one level contribute zero to the normal derivative.

The coordinate comparisons force every label active at the prescribed \(x^0\) to have slope zero when \(h\) is sufficiently small. Fixing such a label gives a constant branch with value \(M=F_h(x^0)\) and hence the global bound \(F_h\le M\). The maximum set \[K=\{x\in\mathbb R^d:F_h(x)=M\}\] is an intersection of finitely many closed affine half-spaces and contains \(x^0\). A zero slope cannot be active on \(\partial B_h\), so \(K\) misses that boundary. Convexity now places all of \(K\) inside \(B_h\): a segment from \(x^0\) to any point outside the box would meet its boundary. Thus \(K\) is a compact positive polytope. Along a solution starting in \(K\), nondecrease of \(F_h\) and the ceiling \(M\) force the solution to remain in the maximum set. A first-exit argument and compact continuation prove invariance and global existence, giving the trajectory bounds. No active-label estimate outside \(B_h\) is needed.

Section 2 proves the approximation lemma by induction on dimension and finite gluing over the minimum exponent, then obtains its uniform form. Section 3 imposes the sign constraints, proves the cut and derivative estimates, constructs the invariant polytope, and deduces the trajectory bounds, positive-equilibrium existence and complex-balanced convergence.

Uniform approximation and sign constraints

Fix henceforth a network, positive rate constants, and \(x^0\in\mathbb R_{>0}^d\) as in Theorem 1. Write \(e_1,\ldots,e_d\) for the standard basis of \(\mathbb R^d\), and set \[ D=\{e_1,\ldots,e_d\} \cup\{z-y:y,z\in\mathcal C,\ y\ne z\}. \tag{6}\] Every vector in \(D\) is nonzero. We include differences of all pairs of complexes, not only reaction vectors. The flux comparison below compares the source of one upward reaction with the sources of every downward reaction. Such sources need not be adjacent in the reaction graph and may even belong to different linkage classes.

For \(r\in[-1,1]^d\), define \[ E(r)=\min\left( \left\{\frac12\right\} \cup \left\{ \frac{\lvert r\cdot v\rvert}{2\lVert v\rVert_1}: v\in D,\ r\cdot v\ne0 \right\}\right), \qquad \lVert v\rVert_1=\sum_{i=1}^d\lvert v_i\rvert. \tag{7}\] This is a positive number for every \(r\), since the displayed minimum is over a finite nonempty set of positive numbers; in particular, \(E(0)=1/2\). The definition is allowed to be discontinuous at the hyperplanes \(r\cdot v=0\). This is the reason for allowing arbitrary positive tolerances in Lemma 2.

Apply that Lemma with \(m=d\), \(a=-1\), \(b=1\), and the tolerance (7). Denote its fixed family by \(\{(r_\lambda,c_\lambda):\lambda\in\Lambda\}\) and put \[F_h(x)=\min_{\lambda\in\Lambda} \bigl(r_\lambda\cdot x+c_\lambda(h)\bigr) \qquad(x\in\mathbb R^d).\] Choose \(h_0\in(0,1)\) small enough that all offsets are defined and Corollary 3 applies for \(0<h<h_0\). Before choosing a particular \(h\), define \[ H=\min\left( \{1\}\cup \left\{ \frac12\lvert r_\lambda\cdot v\rvert: \lambda\in\Lambda,\ v\in D,\ r_\lambda\cdot v\ne0 \right\}\right)>0. \tag{8}\] The positivity of \(H\) again follows from finiteness. Both the slopes and \(H\) are independent of the eventual choice of \(h\). Thus the order of choices is \[D\ \longrightarrow\ E\ \longrightarrow\ \Lambda \ \longrightarrow\ (h_0,H)\ \longrightarrow\ h.\] In particular, \(H\) is not a separation parameter that is allowed to tend to zero together with \(h\).

Suppose that \(0<h<h_0\), \(p\in[-1,1]^d\), and \(\lambda\) is active at \(X_h(p)\). For every \(v\in D\) with \(r_\lambda\cdot v\ne0\), we have \[\begin{align*} \lvert(p-r_\lambda)\cdot v\rvert &\le \lVert p-r_\lambda\rVert_\infty\lVert v\rVert_1\\ &<E(r_\lambda)\lVert v\rVert_1 \le\frac12\lvert r_\lambda\cdot v\rvert. \end{align*}\] Consequently, \[ \operatorname{sgn}(p\cdot v) =\operatorname{sgn}(r_\lambda\cdot v), \qquad \lvert p\cdot v\rvert\ge H. \tag{9}\] In particular, \(\lvert p\cdot v\rvert<H\) forces \(r_\lambda\cdot v=0\). Thus a comparison too close to a tie in exponent space is made an exact tie by the active slope. The converse is not required: equal slope levels may contain complexes with different monomial values. For \(x=X_h(p)\) and any \(y,z\in\mathcal C\), this implies \[ r_\lambda\cdot y<r_\lambda\cdot z \quad\Longrightarrow\quad x^z\le h^H x^y. \tag{10}\] Indeed, \(z-y\in D\) and \(p\cdot(z-y)\ge H\), while \(x^z/x^y=h^{p\cdot(z-y)}\) and \(0<h<1\).

We can now fix a single \(h\in(0,h_0)\) satisfying \[ p_i^0=\frac{\log x_i^0}{\log h}, \qquad \lVert p^0\rVert_\infty<H, \tag{11}\] and, when \(\mathcal R\ne\varnothing\), also satisfying \[ h^H\sum_{e\in\mathcal R}\kappa_e <\min_{e\in\mathcal R}\kappa_e. \tag{12}\] This choice is possible because \(p^0\to0\) and \(h^H\to0\) as \(h\to0+\), with \(H>0\) already fixed. No minimum of reaction rates is needed when \(\mathcal R\) is empty. Write \[B_h=[h,h^{-1}]^d, \qquad F=F_h.\] Since \(H\le1\), Equation (11) places \(x^0\) in the interior of \(B_h\).

Every label active at \(x^0\) has slope zero. In fact, if its \(i\)th slope coordinate were nonzero, Equation (9) applied to \(v=e_i\) would give \(\lvert p_i^0\rvert\ge H\), contradicting (11). The finite nonempty family has at least one active label; fix such a label \(\lambda_0\) and put \(M=F(x^0)\). Then \[ r_{\lambda_0}=0, \qquad M=F(x^0)=c_{\lambda_0}(h), \qquad F(x)\le M\quad(x\in\mathbb R^d). \tag{13}\] Thus a constant branch is forced by the sign constraints, even though it was not prescribed in the application of Lemma 2.

There is no point \(x\in\partial B_h\) at which \(F(x)=M\). To see this, write \(x=X_h(p)\) with \(p\in[-1,1]^d\). The boundary condition is \(\lVert p\rVert_\infty=1\). Equality \(F(x)=M\) would make \(\lambda_0\) active, so Corollary 3 would instead give \(\lVert p\rVert_\infty<E(0)=1/2\). We have proved \[ \{x\in B_h:F(x)=M\}\cap\partial B_h=\varnothing. \tag{14}\] Thus the maximum set \(\{F=M\}\) contains \(x^0\) and misses \(\partial B_h\). To prove its invariance, we next show that a trajectory segment in \(B_h\) cannot decrease \(F\). The constant branch alone does not prove this: another branch tied with it might otherwise become smaller. We must control the derivatives of all active branches. The next two arguments establish that control from weak reversibility.

Flux across the ordered complex levels

The following finite graph statement isolates the role of weak reversibility. Its weights will be the monomials \(x^y\), and its levels will be the scalar products \(r_\lambda\cdot y\).

Return-path monomial domination appears in Anderson’s tier argument (Anderson 2011a, Lemma 4.7), and quantitative inward reaction estimates underlie the invariant polygons of Craciun, Nazarov, and Pantea (Craciun et al. 2013). Ordered-cut regrouping also appears in the proof of Craciun’s toric-inclusion embedding theorem (Craciun 2019, preprint, Theorem 4.1). Lemma 4 uses a direct finite estimate: one upward edge dominates the full downward flux. No tier extraction, cycle decomposition, or differential inclusion is required in the proof above.

For any \(x\in B_h\) and any label \(\lambda\) active at \(x\), apply Lemma 4 with \[V=\mathcal C,\quad \mathcal E=\mathcal R,\quad a(y)=r_\lambda\cdot y,\quad w_y=x^y,\quad \delta=h^H.\] The return paths are supplied by weak reversibility, Equation (10) supplies (15), and Equation (12) supplies the rate inequality when there are reactions. Thus \[ r_\lambda\cdot f(x) =\sum_{y\to y'\in\mathcal R} \kappa_{y\to y'}x^y\,r_\lambda\cdot(y'-y) \ge0 \qquad(x\in B_h,\ \lambda\text{ active at }x). \tag{17}\] This includes \(f=0\) for an empty reaction set. The comparison in (15) is required for all ordered vertex pairs, including nonadjacent pairs in one component and pairs in different components; this explains the full set of complex differences in (6).

Derivatives of finite affine minima

Although \(F\) need not be differentiable where its active labels change, its right derivative along a differentiable curve is explicit. We record the finite-affine minimum case of the classical active-branch directional-derivative principle; see, for example, Bertsekas (Bertsekas 2009, Proposition A.3.2(a)) for the corresponding maximum formula. The short proof includes all ties and the monotonicity implication needed here.

The invariant polytope and trajectory bounds

All active branches have nonnegative derivative in \(B_h\), while the constant branch in (13) gives \(F\le M\) everywhere. We now show that the maximum set of \(F\) is the compact invariant polytope announced in the introduction. Its invariance will give Theorem 1.

Nonnegative derivatives on all active inequalities express the familiar inward-tangency condition for invariance, whose classical general form goes back to Nagumo (Nagumo 1942). Here we prove invariance directly by a first-exit argument.

Equilibria and complex-balanced convergence

The convexity of the invariant set also recovers the known existence of positive equilibria in every positive stoichiometric class, proved by Boros (Boros 2019, preprint, Theorem 1).