Boundedness and persistence of weakly reversible mass-action systems

1 Introduction

Weak reversibility is a condition on the directed graph of a chemical reaction network: each reaction has a directed return path. Two central questions ask whether this condition alone prevents species extinction under mass-action kinetics and whether every positive trajectory remains bounded. We prove the joint boundedness and persistence conclusion for arbitrary positive constant reaction rates, any number of connected components in the reaction graph, and every species dimension.

1.1 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

Since the exponents are nonnegative integers, f is a polynomial vector field. A forward trajectory in \R_{>0}^d is called persistent if

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

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.1 applies whether or not P is bounded. It therefore includes the persistence question for bounded entire compatibility classes. The estimate (1.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 (1.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 3.3: each x^0>0 belongs to a compact convex polytope K\subset\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 [5]: every positive stoichiometric class contains a positive equilibrium (Corollary 3.4).

1.2 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 [16], Horn’s complex-balancing criteria [15], and Feinberg’s deficiency-zero analysis [10] 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.1 under the additional complex-balance assumption; Remark 3.5 gives the argument.

Feinberg [12, 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 [1, Section 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.1 establishes both conclusions for arbitrary positive constant reaction rates.

August and Barahona [3, 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

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 [1] proved boundedness for weakly reversible systems with a single linkage class. Gopalkrishnan, Miller, and Shiu [14] proved permanence for the larger class of strongly endotactic systems, which includes every weakly reversible single-linkage system; Boros and Hofbauer [6] gave a later proof of the single-linkage result. Craciun, Nazarov, and Pantea [9] proved boundedness, persistence, and permanence for two-species weakly reversible systems as part of their more general endotactic result. Pantea [18] 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 [8, Section 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 [9], and projected polyhedra in Pantea [18]. Anderson’s monomial tiers and return-path domination [2], the projection and multiscale arguments of Gopalkrishnan, Miller, and Shiu [13, 14], and Craciun’s ordered-complex embedding into toric differential inclusions [7] 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.

1.3 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

These points fill the box B_h=[h,h^{-1}]^d. We construct a concave piecewise-affine function

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.1 and Corollary 2.2 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 \norm{p-r_\lambda}_\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,

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

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.

2 Affine labels with a prescribed tolerance

We construct a finite minimum of affine functions whose minimizing slopes approximate the exponents of a positive point. The permitted error is prescribed separately at each slope. Allowing this prescription to be discontinuous will let us impose sign conditions on finitely many linear forms in the application.

For an integer m\geq 1, u\in\R^m, and 0<h<1, write

A label is a pair (r,c) consisting of a vector r\in\R^m, independent of h, and a real-valued offset c(h) defined for all sufficiently small positive h. For a finite nonempty family \Lambda of labels, write its members as \lambda=(r_\lambda,c_\lambda). For fixed sufficiently small h>0 and x\in\R^m, a label \lambda\in\Lambda is active at x if

Thus every minimizing label is active when a tie occurs. Whenever finitely many families are combined, we restrict their offsets to a common interval 0<h<h_* on which they are all defined.

Two aspects of the statement will be used separately. The slopes form a fixed finite set, rather than a set that changes with h. Also, the conclusion concerns every label that stays active along a sequence, including labels tied with other minimizers. The freedom for p^\nu to leave the cube is needed by the induction: a projected sequence may approach a face of a smaller cube from either side.

The sequential conclusion gives a uniform estimate on the entire cube, with the same fixed finite family. This is the form used for reaction networks.

In Section 3, we choose the tolerance so that every nonzero dot product of an active slope with a prescribed vector forces the exponent vector to have a dot product of the same sign. The freedom to choose a discontinuous E permits this requirement near the corresponding hyperplanes.

3 From affine labels to persistent trajectories

We now apply Lemma 2.1 to the mass-action vector field. Its uniform form in Corollary 2.2 lets us choose the tolerance so that every strict ordering made by an active slope gives a monomial comparison with a definite gap. Weak reversibility converts those comparisons into nonnegative derivatives, and an initially active constant branch supplies the compact set that traps the trajectory.

3.1 Uniform approximation and sign constraints

Fix henceforth a network, positive rate constants, and x^0\in\R_{>0}^d as in Theorem 1.1. Write e_1,\ldots,e_d for the standard basis of \R^d, and set

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

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

Apply that Lemma with m=d, a=-1, b=1, and the tolerance (3.2). Denote its fixed family by \{(r_\lambda,c_\lambda):\lambda\in\Lambda\} and put

Choose h_0\in(0,1) small enough that all offsets are defined and Corollary 2.2 applies for 0<h<h_0. Before choosing a particular h, define

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

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

Consequently,

In particular, \abs{p\cdot v}<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

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

and, when \mathcal R\ne\varnothing, also satisfying

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

Since H\le1, Equation (3.6) places x^0 in the interior of B_h.

Every label active at x^0 has slope zero. In fact, if its ith slope coordinate were nonzero, Equation (3.4) applied to v=e_i would give \abs{p_i^0}\ge H, contradicting (3.6). The finite nonempty family has at least one active label; fix such a label \lambda_0 and put M=F(x^0). Then

Thus a constant branch is forced by the sign constraints, even though it was not prescribed in the application of Lemma 2.1.

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 \norm{p}_\infty=1. Equality F(x)=M would make \lambda_0 active, so Corollary 2.2 would instead give \norm{p}_\infty<E(0)=1/2. We have proved

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.

3.2 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 [2, Lemma 4.7], and quantitative inward reaction estimates underlie the invariant polygons of Craciun, Nazarov, and Pantea [9]. Ordered-cut regrouping also appears in the proof of Craciun’s toric-inclusion embedding theorem [7, preprint, Theorem 4.1]. Lemma 3.1 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 3.1 with

The return paths are supplied by weak reversibility, Equation (3.5) supplies (3.10), and Equation (3.7) supplies the rate inequality when there are reactions. Thus

This includes f=0 for an empty reaction set. The comparison in (3.10) 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 (3.1).

3.3 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 [4, Proposition A.3.2(a)] for the corresponding maximum formula. The short proof includes all ties and the monotonicity implication needed here.

3.4 The invariant polytope and trajectory bounds

All active branches have nonnegative derivative in B_h, while the constant branch in (3.8) 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.1.

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

3.5 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 [5, preprint, Theorem 1].