Species X_1, ..., X_d have nonnegative concentrations x_i, and x = (x_1, ..., x_d). A complex is a vector y = (y_1, ..., y_d) in Z_{>=0}^d. The network has a finite set C of complexes and a finite set R of reactions: directed edges y -> z with y, z in C and y != z. Each reaction y -> z carries a positive constant rate kappa_{y->z} > 0. The mass-action system is

    dx_i/dt = f_i(x) = sum_{y->z in R} kappa_{y->z} x^y (z_i - y_i),   (1)

for i = 1, ..., d, with monomials x^y = prod_{i=1}^d x_i^{y_i}. Write f = (f_1, ..., f_d). The vector field f is polynomial on R^d, hence smooth, and each trajectory extends to a unique maximal forward solution x : [0, T_max) -> R^d with T_max in (0, +infinity]. If T_max < infinity, then ||x(t)|| -> infinity as t -> T_max.

The positive orthant is forward invariant. On the hyperplane {x_i = 0}, the i-th component reduces to a sum over reactions with y_i = 0, since reactions with y_i >= 1 contribute a monomial divisible by x_i; each surviving term is nonnegative. Consequently x(0) in R_{>0}^d implies x_i(t) > 0 for all t in [0, T_max).

The stoichiometric subspace S is the real span of the vectors z - y, y -> z in R. The vector field takes values in S, so x(t) stays in the affine class x(0) + S.

The network is weakly reversible if every y -> z in R admits a directed path y_0 = z, y_1, ..., y_k = y with y_{j-1} -> y_j a reaction of R. Equivalently: every arrow is undone by a path of arrows. This is a property of the reaction graph alone and involves neither the rate constants nor the initial condition.

Theorem. Let C be finite, let R be weakly reversible, let kappa_{y->z} > 0 be arbitrary positive constants, and let x(0) in R_{>0}^d. Then the solution of (1) exists for all t >= 0, and there is eps in (0, 1), depending only on C, R, kappa, and x(0), such that

    eps <= x_i(t) <= 1/eps    for all t >= 0 and i = 1, ..., d.   (2)

pdf: pdfrefs/anderson-2011.pdf, pdfrefs/cnp-2013.pdf

The upper bound is the boundedness conjecture of Anderson (2011). The conjunction with the lower bound is the persistence conjecture stated by Feinberg (1987, Remark 6.1.E) and separated into an independent question by Anderson (2011). Neither bound is assumed in the statement; both follow from weak reversibility.

No assertion is made about convergence. The trajectory is trapped in a fixed compact box but need not approach a single point; convergence to an equilibrium requires complex balance or a further structural hypothesis.

Earlier work covers only restricted subclasses. Anderson (2011) proved boundedness for single-linkage systems. Gopalkrishnan, Miller, and Shiu (2013, 2014) proved permanence for the larger class of strongly endotactic systems, and Boros and Hofbauer (2020) gave an alternative proof of the single-linkage case. Craciun, Nazarov, and Pantea (2013) covered all two-species systems. Pantea (2012) proved persistence of bounded weakly reversible trajectories when dim S = 2, for arbitrary species and linkage counts. Craciun (2019, 2026) developed a toric differential inclusion whose weakly reversible case gives persistence conditional on a preexisting upper bound. Boros (2019) proved the existence of a positive equilibrium in every positive compatibility class, reproduced below as a consequence of the convexity of the invariant region. August and Barahona (2010, Theorem 6) stated the general conclusion but their argument fails at the point noted below.

pdf: pdfrefs/gms-2013.pdf, pdfrefs/gms-2014.pdf, pdfrefs/cnp-2013.pdf, pdfrefs/pantea-2012.pdf, pdfrefs/craciun-2019.pdf, pdfrefs/craciun-2026.pdf, pdfrefs/boros-2019.pdf, pdfrefs/anderson-2011.pdf

The obstruction. A single reaction need not point inward across any proposed face of the orthant. A single linear functional r . x - c therefore gives a derivative r . f(x) whose sign can fail; individual summands of r . f(x) can be negative. Weak reversibility supplies return paths whose summands cancel the negative ones once monomials are ordered correctly. The whole strategy is to convert "return paths exist" into a geometric inequality r_lambda . f(x) >= 0 valid at every active label lambda and every x in a compact box.

Why a naive monomial argument fails. August and Barahona (2010) proposed to infer monomial dominance from a maximal sum of exponents together with a "sufficiently large product of concentrations". In 2A <-> B with both rates one and a = [A], b = [B], the state (a, b) = (T, T^3) for T > 1 has ab = T^4 -> infinity but a^2 = T^2 < T^3 = b, and

    d/dt (a + b) = b - a^2 = T^3 - T^2 > 0.

So the sum grows at (T, T^3) even though the smaller monomial a^2 is at the level of b. The remedy is not to compare arbitrary monomials, but to compare them with a definite slack determined by an ACTIVE slope: for the correct active slope r_lambda, a strict ordering r_lambda . y < r_lambda . z forces the sharp inequality x^z <= h^H x^y with a fixed H > 0, indexed on a small parameter h in (0, 1).

Strategy. Fix h in (0, 1) small. Points in the positive box B_h = [h, h^{-1}]^d are parametrized by

    X_h(p) = ( h^{p_1}, ..., h^{p_d} ),    p in [-1, 1]^d,

so that x^y = h^{p . y}. The construction defines

    F(x) = min over labels lambda in Lambda of ( r_lambda . x + c_lambda(h) ),

where Lambda is a finite family of labels, each with a fixed slope r_lambda in [-1, 1]^d and an offset c_lambda(h) that depends on h, chosen by the affine approximation lemma (Lemma 1). Active slopes of F at X_h(p) are forced to be E-close to p, where the tolerance E can be prescribed independently at each slope and may be discontinuous. This freedom is what allows the active slope to see monomial comparisons with the correct sign: for finitely many test vectors v in a set D containing {e_1, ..., e_d} and every difference z - y with y, z in C distinct, one has sgn(r_lambda . v) = sgn(p . v) and |p . v| >= H.

Nondecrease of F along the trajectory uses a finite cut-flux argument (Lemma 2): the sum over reactions of kappa_{y->z} x^y (r_lambda . z - r_lambda . y) is nonnegative whenever the monomials satisfy the sharp inequality above and a rate-separation condition (h^H sum kappa < min kappa). Then a right-derivative formula for finite affine minima (Lemma 3) turns the active-label inequality into nondecrease of F o x(t) inside B_h, via the classical active-branch directional-derivative principle (Bertsekas 2009).

pdf: https://web.mit.edu/dimitrib/www/Convex_Theory_Entire_Book.pdf

The maximum set K = { x : F(x) = M } with M = F(x(0)) is a closed convex polyhedron. Once every label active at x(0) has slope zero (arranged by choosing h small enough so that log(x_i(0)) / log(h) is small compared to H), the value M = c_{lambda_0}(h) is attained on a constant branch, K avoids the boundary of B_h, and convexity places K in the interior of B_h. Positivity, boundedness, and forward invariance of K follow. The last step, and the reason the invariance is nontrivial at ties, is that at an active label switching between branches, the right derivative of F o x is the minimum over all active branches and is therefore controlled by the cut-flux inequality for every active branch.

The affine approximation lemma. For p in R^m and 0 < h < 1, write X_h(p) = (h^{p_1}, ..., h^{p_m}) in R^m_{>0}. A label is a pair (r, c) where r in R^m is a fixed slope and c = c(h) is a real-valued offset defined for all sufficiently small h > 0. For a finite family Lambda of labels, a label lambda = (r_lambda, c_lambda) is active at x in R^m for a given h if

    r_lambda . x + c_lambda(h) = min over eta in Lambda of ( r_eta . x + c_eta(h) ).

Every minimizer is active: ties count.

Lemma 1. Let m >= 1, let a < b be real numbers, and let E : [a, b]^m -> (0, infinity) be any positive function. Then there is a finite nonempty family Lambda of labels (r_lambda, c_lambda) with r_lambda in [a, b]^m and c_lambda(h) = o(h^a) as h -> 0+, such that for every pair of sequences h_nu -> 0+ and p^nu in R^m converging to p in [a, b]^m, if a fixed lambda in Lambda is active at X_{h_nu}(p^nu) for every nu, then

    || p - r_lambda ||_infinity < E(r_lambda).   (3)

No regularity of E is assumed, and the approximating vectors p^nu need not belong to [a, b]^m.

Two aspects of the statement serve separate purposes. The slopes r_lambda form a fixed finite set, not a family that changes with h. The conclusion concerns every label that is active along a sequence, so tied labels are counted. The freedom for p^nu to leave the cube is needed inside the induction: a projected sequence can approach a face of a smaller cube from either side.


Proof of Lemma 1. Induction on m. Fix m, a < b, E, and write mu(u) = min_i u_i. Fix t in [a, b] and construct a finite family Lambda_t of labels with the following three properties:

  L1. every slope r_lambda in Lambda_t has mu(r_lambda) = t;

  L2. the label (t*1, 0) is in Lambda_t, where 1 = (1, ..., 1) in R^m;

  L3. if p in [a, b]^m has mu(p) = t, if h_nu -> 0+ and p^nu -> p, and if a
      fixed lambda in Lambda_t is active at X_{h_nu}(p^nu) for every nu,
      then || p - r_lambda ||_infinity < E(r_lambda).

Base case. If m = 1 or t = b, take Lambda_t = {(t*1, 0)}. Then mu(p) = t with p in [a, b]^m forces p = t*1, so || p - r_lambda ||_infinity = 0 < E(r_lambda). Assume henceforth m > 1 and t < b.

Types. A label in Lambda_t has type k (0 <= k <= m-1) if exactly k of its slope entries are strictly greater than t and the remaining m - k entries equal t. Types are constructed in increasing k.

Parameters. Choose t < beta_{m-1} < beta_{m-2} < ... < beta_1 < b so that

    beta_k - t < E(r)   for every slope r of type k - 1.   (4)

At the step choosing beta_k, types 0, ..., k-1 are finite and every E(r) is positive, so (4) imposes only finitely many strict conditions and such a choice is possible.

Type 0. The single label (t*1, 0).

Type k (1 <= k <= m-1). For each k-element subset U of {1, ..., m} and each s in [beta_k, b]^U, define r(s) in R^m by r(s)_i = s_i on U and r(s)_i = t off U. Apply the induction hypothesis in dimension k, to the cube [beta_k, b]^U and tolerance s -> E(r(s)), yielding a finite family of k-dimensional labels (s_j, d_j) with s_j in [beta_k, b]^U and d_j(h) = o(h^{beta_k}). Put in Lambda_t the m-dimensional label

    r = r(s_j),    c_lambda(h) = -h^{beta_k} + d_j(h).   (5)

All k-element subsets U and all resulting local labels are combined; the family is finite.

L1 holds because every slope entry of a type-k label is either s_i >= beta_k > t on U, or t off U, so min_i r_i = t. L2 is the type-0 label. Every offset is o(h^{beta_k}) for its type k and hence o(h^{beta_{k-1}}), ..., o(h^{beta_1}); in particular, every offset is o(h^t), and c_lambda = 0 trivially for the type-0 label.

Verify L3. Let p in [a, b]^m with mu(p) = t, let h_nu -> 0+ and p^nu -> p, and let lambda = (r, c) be a type-k label of Lambda_t active at X_{h_nu}(p^nu) for every nu. Let U be its raised set: the k indices with r_i > t, so r_i >= beta_k on U. Write h for h_nu.

Claim 1. If k >= 1 and i in U, then p_i >= beta_k.

Compare lambda to the center (t*1, 0). Activity gives

    sum_{i in U} (r_i - t) h^{p_i^nu} + c_lambda(h) <= 0.

Substituting (5),

    sum_{i in U} (r_i - t) h^{p_i^nu} <= h^{beta_k} - d_j(h).   (6)

Every r_i - t >= beta_k - t > 0. If some p_i < beta_k with i in U, fix delta > 0 with p_i^nu <= beta_k - delta for all large nu; then the i-th term of the left side of (6), after division by h^{beta_k}, is at least (beta_k - t) h^{-delta} -> +infinity, while the right side divided by h^{beta_k} is 1 + o(1). Contradiction. So p_i >= beta_k.

Claim 2. If k <= m - 2 and j not in U, then p_j <= beta_{k+1}.

Suppose p_j > beta_{k+1} for some j not in U. Construct a type-(k+1) competitor with raised set U U {j}, slope equal to r on U, an extra raised coordinate j of value s in [beta_{k+1}, b], and offset c'(h) = -h^{beta_{k+1}} + d'(h) as in (5). Activity of lambda forces competitor minus lambda >= 0. On U the slopes cancel; coordinate j contributes (s - t) h^{p_j^nu}; the offsets contribute

    c'(h) - c_lambda(h) = h^{beta_k} - h^{beta_{k+1}} + d'(h) - d_j(h)
                        = -h^{beta_{k+1}} (1 + o(1)),   (7)

since h^{beta_k} = o(h^{beta_{k+1}}), d'(h) = o(h^{beta_{k+1}}), and d_j(h) = o(h^{beta_k}) = o(h^{beta_{k+1}}). (For k = 0 the same computation applies with c_lambda = 0.) Also p_j > beta_{k+1} gives h^{p_j^nu} = o(h^{beta_{k+1}}). Hence

    competitor - lambda = o(h^{beta_{k+1}}) - h^{beta_{k+1}} (1 + o(1)) < 0

for all small h. Contradiction, so p_j <= beta_{k+1}.

Claim 3. If k = m - 1, the unique j not in U satisfies p_j = t. By Claim 1, every i in U has p_i >= beta_{m-1} > t; with mu(p) = t, the minimum is attained only at j, giving p_j = t.

Conclusion of L3. The induction hypothesis applied to the k-labeled family in dimension k gives, on the U-coordinates, | p_i - r_i | < E(r) for i in U: by Claim 1, p_U^nu lies eventually in [beta_k - delta', b]^U for any delta' > 0, and the activity of (s_j, d_j) within its own family is preserved, since all labels of Lambda_t at the same level receive the same shift -h^{beta_k}. Off U, Claims 2 and 3 give | p_i - r_i | = | p_i - t | <= beta_{k+1} - t < E(r) for k <= m - 2, and | p_i - r_i | = 0 < E(r) for k = m - 1. Combining the two coordinate sets proves || p - r_lambda ||_infinity < E(r_lambda). L3 holds.

Step 3. Neighborhood of t. For lambda in Lambda_t let A_lambda be the set of all limits p in [a, b]^m of sequences p^nu in R^m with h_nu -> 0+ for which lambda is active at X_{h_nu}(p^nu) within Lambda_t. The set A_lambda is closed: if z_n in A_lambda and z_n -> z, choose from a realizing sequence for each z_n a pair (h_n, u_n) with h_n < 1/n, || u_n - z_n ||_infinity < 1/n, and lambda active at X_{h_n}(u_n); then u_n -> z, so z is also a limit of active points. The argument uses actual active points, so no continuity of c_lambda is required.

Define

    B_t = union over lambda = (r, c) in Lambda_t of
          { p in A_lambda : || p - r ||_infinity >= E(r) }.

B_t is compact: finite union of closed subsets of the compact [a, b]^m. The fixed-minimum verification shows B_t contains no point with mu(p) = t. Since mu is continuous, there is rho_t > 0 such that B_t contains no point with

    mu(p) in I_t := (t - rho_t, t + rho_t).   (8)

If B_t is nonempty, any rho_t smaller than the positive minimum of | mu(p) - t | on B_t works; if B_t is empty, any positive radius initially works. Thus the limiting-activity conclusion holds whenever p in [a, b]^m and mu(p) in I_t.

Shrink rho_t if necessary so that every offset of Lambda_t also satisfies

    c_lambda(h) = o(h^{t + rho_t}).   (9)

For a larger family this is achieved by requiring t + rho_t < min_{1 <= k <= m-1} beta_k; for a singleton family the offset is zero. The stronger decay uses the constructed form c_lambda = -h^{beta_k} + o(h^{beta_k}) with beta_k > t; it does not follow from c_lambda = o(h^t) alone. Shrinking the radius preserves (8).

Step 4. Gluing. The intervals I_t as t ranges over [a, b] cover [a, b], because t lies in I_t. A finite subcover is extracted. Intervals contained in another are removed, keeping one copy if two coincide. The surviving centers are ordered as a <= t_1 < t_2 < ... < t_N <= b, and their intervals are written I_{t_j} = (L_j, R_j). Noncontainment forces L_1 < ... < L_N and R_1 < ... < R_N (the two orders are tied together by R_j - L_j = 2 rho_{t_j} and by the ordering of centers). The equality R_j = L_{j+1} would leave the point L_{j+1} uncovered, so R_j > L_{j+1}: successive intervals overlap.

For 1 <= j < N, choose a cut

    q_j in ( max(t_j, L_{j+1}), min(t_{j+1}, R_j) ).   (10)

The open interval is nonempty because I_{t_j} contains its own center and overlaps I_{t_{j+1}}, which contains its own center. Set q_0 = a and q_N = b. Then

    t_j < q_j < t_{j+1}   (1 <= j < N),                (11)

    [q_{j-1}, q_j] subset I_{t_j}   (1 <= j <= N).    (12)

If N = 1, no cuts are required and [a, b] = I_{t_1}.

Gluing. Combine the selected families by taking their disjoint union and adding to each label of level t_j the constant shift:

    c~_lambda(h) = c_lambda(h) - sum_{l=1}^{j-1} h^{q_l},     (13)

with the empty sum interpreted as zero. Slopes are unchanged. Every new offset is o(h^a): by (9), c_lambda = o(h^{t_j + rho_{t_j}}) with t_j + rho_{t_j} > a, and every shift exponent satisfies q_l > t_l >= a, so h^{q_l} = o(h^a). A finite sum of o(h^a) terms remains o(h^a).

Fix a label from level t_j that is active in the glued family at X_{h_nu}(p^nu) along h_nu -> 0+ and p^nu -> p in [a, b]^m. The claim is

    q_{j-1} <= mu(p) <= q_j.   (14)

Upper bound. Suppose j < N and mu(p) > q_j. Then h^{p_i^nu} = o(h^{q_j}) in every coordinate. By (9) and the membership q_j in I_{t_j}, the original offset satisfies c_lambda = o(h^{q_j}). Compare the shifted center (t_{j+1}*1, 0) with its shift -sum_{l=1}^{j} h^{q_l} to the glued label lambda with its shift -sum_{l=1}^{j-1} h^{q_l}. The difference (center minus label) is

    (t_{j+1}*1 - r) . X_h(p^nu) - c_lambda(h) - h^{q_j}.   (15)

Every coefficient t_{j+1} - r_i is bounded by |b - a|, and every h^{p_i^nu} = o(h^{q_j}), so the first term is o(h^{q_j}). The second term is o(h^{q_j}). Hence the difference equals -h^{q_j} (1 + o(1)) < 0 for small h, contradicting activity of lambda.

Lower bound. Suppose j > 1 and mu(p) < q_{j-1}. Then some p_i < q_{j-1}. Compare the shifted center (t_{j-1}*1, 0) with its shift -sum_{l=1}^{j-2} h^{q_l} to the glued label. The difference (center minus label) is

    (t_{j-1}*1 - r) . X_h(p^nu) - c_lambda(h) + h^{q_{j-1}}.   (16)

Every r_i >= t_j (L1) and so t_{j-1} - r_i <= -(t_j - t_{j-1}) < 0. After division by h^{q_{j-1}}, the i-th contribution of the scalar product is -sum_i (r_i - t_{j-1}) h^{p_i^nu} / h^{q_{j-1}} <= -(t_j - t_{j-1}) h^{p_i^nu - q_{j-1}} for the specific i with p_i < q_{j-1}, hence -> -infinity. The second term divided by h^{q_{j-1}} vanishes by (9) and q_{j-1} < t_j + rho_{t_j}. The third term is +1. So the difference is negative for small h, contradicting activity.

Combining (14) with (12) gives mu(p) in I_{t_j}. The label is active within its original family, because all labels of level t_j received the same shift. The neighborhood guarantee at I_{t_j} then gives || p - r_lambda ||_infinity < E(r_lambda). Every contradiction constructed above produced a strictly smaller competing value, so the argument covers every active label even at ties. This completes the induction and proves Lemma 1.

Corollary. Uniform form. With E : [a, b]^m -> (0, infinity) any positive function and Lambda as in Lemma 1, there is h_0 > 0 such that for every 0 < h < h_0 and every p in [a, b]^m, every label lambda in Lambda active at X_h(p) satisfies

    || p - r_lambda ||_infinity < E(r_lambda).   (17)

Proof. If not, there are h_nu -> 0+ (with h_nu < any candidate h_0), points p^nu in [a, b]^m, and labels lambda_nu in Lambda active at X_{h_nu}(p^nu), with || p^nu - r_{lambda_nu} ||_infinity >= E(r_{lambda_nu}). By compactness of [a, b]^m and finiteness of Lambda, a subsequence has lambda_nu = lambda fixed and p^nu -> p in [a, b]^m. The limit contradicts the strict inequality of Lemma 1 at the fixed positive value E(r_lambda). The value E(r_lambda) is fixed on this subsequence, so no regularity of E is used.

pdf: pdfrefs/anderson-2011.pdf

Lemma 2 (cut-flux). The statement isolates the contribution of weak reversibility and is independent of mass-action structure. Let V be a finite set, E a finite set of directed edges with source and target maps s, t : E -> V, and a : V -> R a level. Assign weights w_v > 0 for v in V and rates kappa_e > 0 for e in E. Suppose every edge has a directed path from its target t(e) back to its source s(e). Suppose that delta > 0 satisfies

    a(u) < a(z)   implies   w_z <= delta w_u     (u, z in V),    (18)

and, when E is nonempty,

    delta * sum_{e in E} kappa_e  <  min_{e in E} kappa_e.   (19)

Then

    sum_{e in E} kappa_e w_{s(e)} ( a(t(e)) - a(s(e)) )  >=  0.    (20)

pdf: pdfrefs/anderson-2011.pdf

Proof. If E is empty, both sides vanish. Otherwise list the distinct values of a on V as alpha_1 < alpha_2 < ... < alpha_q. For 1 <= j < q, set

    V_j^- = { v in V : a(v) <= alpha_j },
    V_j^+ = { v in V : a(v) > alpha_j },

let E_j^+ be the edges from V_j^- to V_j^+ and E_j^- the edges from V_j^+ to V_j^-, and define the net upward flux

    Q_j = sum_{e in E_j^+} kappa_e w_{s(e)} - sum_{e in E_j^-} kappa_e w_{s(e)}.   (21)

If E_j^- is empty then Q_j >= 0 as a sum of positive contributions. If E_j^- is nonempty, take any e_in in E_j^-; its return path starts at t(e_in) in V_j^- and ends at s(e_in) in V_j^+, so it contains an edge e_0 in E_j^+. Put u = s(e_0) in V_j^-. Every source z of an edge in E_j^- satisfies a(z) > alpha_j >= a(u), so a(u) < a(z) and w_z <= delta w_u. Hence

    sum_{e in E_j^-} kappa_e w_{s(e)}
        <= delta w_u sum_{e in E_j^-} kappa_e
        <= delta w_u sum_{e in E} kappa_e
        < w_u min_{e in E} kappa_e
        <= kappa_{e_0} w_u.   (22)

The last term is a single summand of sum_{e in E_j^+} kappa_e w_{s(e)}, so Q_j > 0. In both cases, Q_j >= 0.

The selected edge e_0 need not share a linkage class with every edge in E_j^-: condition (18) applies to all ordered pairs, allowing one upward summand to dominate the entire downward sum.

Every edge e satisfies the telescoping identity

    a(t(e)) - a(s(e)) = sum_{j: e in E_j^+} (alpha_{j+1} - alpha_j) - sum_{j: e in E_j^-} (alpha_{j+1} - alpha_j).   (23)

An upward edge from level j_s to j_t > j_s crosses the intermediate cuts upward, and alpha_{j_t} - alpha_{j_s} equals the corresponding sum of gaps; a downward edge crosses them downward; an edge with j_t = j_s contributes zero on both sides. Weight (23) by kappa_e w_{s(e)} and sum over e: the sum of the upward-cut terms is sum_j (alpha_{j+1} - alpha_j) sum_{e in E_j^+} kappa_e w_{s(e)}, the sum of the downward-cut terms is analogous over E_j^-, and the difference is sum_j (alpha_{j+1} - alpha_j) Q_j >= 0 by (21) and (22). This is (20).

Setup. Let I be a finite nonempty index set with vectors b_i in R^d and reals d_i. Define Phi : R^d -> R by

    Phi(x) = min_{i in I} ( b_i . x + d_i ).   (24)

Phi is concave, continuous, and piecewise affine. For x in R^d, write A(x) = { i in I : b_i . x + d_i = Phi(x) } for the active index set.

Lemma 3 (derivative). Let gamma : [u, v] -> R^d be continuously differentiable. At every tau in [u, v), the right derivative of Phi o gamma exists and equals

    (d/dtau)^+ Phi(gamma(tau)) = min_{i in A(gamma(tau))} b_i . gamma'(tau).   (25)

If b_i . gamma'(tau) >= 0 for every tau in [u, v) and every i in A(gamma(tau)), then Phi o gamma is nondecreasing on [u, v].

pdf: https://web.mit.edu/dimitrib/www/Convex_Theory_Entire_Book.pdf

Proof. Fix tau in [u, v) and put A = A(gamma(tau)). Every inactive i has b_i . gamma(tau) + d_i > Phi(gamma(tau)) + c_i for some c_i > 0. Continuity of gamma then gives, for tau' > tau sufficiently close to tau, b_i . gamma(tau') + d_i >= Phi(gamma(tau')) + c_i / 2. So only active indices can determine Phi(gamma(tau')) for tau' near tau.

Let h = tau' - tau > 0. Write gamma(tau') = gamma(tau) + h gamma'(tau) + o(h). For an active i,

    b_i . gamma(tau') + d_i = Phi(gamma(tau)) + h b_i . gamma'(tau) + o(h).   (26)

Taking the minimum over i in A and using that only A can be active,

    Phi(gamma(tau')) - Phi(gamma(tau))
        = h min_{i in A} b_i . gamma'(tau) + o(h).   (27)

Divide by h and take h -> 0+. This proves (25). The monotonicity statement follows: the right derivative at every tau is nonnegative, so Phi o gamma is nondecreasing (right derivative nonnegativity on an interval implies monotonicity by the standard mean-value argument applied on compact subintervals).

The short proof is the active-branch directional-derivative principle specialized to a min of affine functions. The maximum of affine functions yields the analogous "max of active derivatives" formula, as in Bertsekas (2009, Prop. A.3.2(a)); the extension to finite minima is by taking Phi = - max (-b_i . x - d_i).

Application. Fix the network (C, R) with positive rates kappa_{y->z}, and fix x(0) in R_{>0}^d. Prepare the numerical data as follows.

Directions. Let e_1, ..., e_d be the standard basis of R^d and put

    D = { e_1, ..., e_d } union { z - y : y, z in C, y != z }.   (28)

D is finite and contains no zero vector.

Tolerance. For r in [-1, 1]^d define

    E(r) = min( {1/2} union { |r . v| / (2 ||v||_1) : v in D, r . v != 0 } ),   (29)

where ||v||_1 = sum_i |v_i|. Because D is finite and every |r . v| in the second set is positive, E(r) is positive for every r; E(0) = 1/2. E is discontinuous on the hyperplanes r . v = 0, which the affine lemma permits because it assumes no regularity of E.

Labels. Apply Lemma 1 with m = d, a = -1, b = 1, tolerance E as in (29). This yields a finite nonempty label family Lambda with slopes r_lambda in [-1, 1]^d and offsets c_lambda(h) = o(h^{-1}).

Threshold. Let h_0 > 0 be the threshold of the uniform corollary: for every 0 < h < h_0 and every p in [-1, 1]^d, every active label lambda of Lambda at X_h(p) satisfies

    || p - r_lambda ||_infinity < E(r_lambda).   (30)

Separation. Define

    H = min( {1} union { |r_lambda . v| / 2 : lambda in Lambda, v in D, r_lambda . v != 0 } ).   (31)

Finiteness gives H > 0; also H <= 1.

Order of choices: C, R -> D -> E -> Lambda -> (h_0, H) -> final h. In particular, H is fixed before any h.

Sign constraint. Fix h in (0, h_0) and p in [-1, 1]^d, let lambda in Lambda be active at x = X_h(p), and let v in D satisfy r_lambda . v != 0. Then

    sgn(p . v) = sgn(r_lambda . v)    and    |p . v| >= H.   (32)

Proof. By (30), || p - r_lambda ||_infinity < E(r_lambda). Compute

    |(p - r_lambda) . v| <= || p - r_lambda ||_infinity ||v||_1
                          < E(r_lambda) ||v||_1
                          <= |r_lambda . v| / 2   (by (29)).
So |p . v - r_lambda . v| < |r_lambda . v| / 2, which gives sgn(p . v) = sgn(r_lambda . v) (the same side of zero) and |p . v| >= |r_lambda . v| / 2. By (31), H <= |r_lambda . v| / 2, so |p . v| >= H. Done.

Monomial order. Let y, z in C be distinct and suppose r_lambda . y < r_lambda . z. Then for x = X_h(p) with 0 < h < h_0 and lambda active at x,

    x^z <= h^H x^y.   (33)

Proof. Because z - y is in D (28) and r_lambda . (z - y) > 0, (32) with v = z - y gives p . (z - y) >= H. Hence

    x^z / x^y = h^{p . (z - y)} <= h^H.   (34)

Done.

Construction of the invariant polytope. Fix h satisfying three conditions:

  (a) p^0 in R^d defined by p_i^0 = log x_i(0) / log h satisfies || p^0 ||_infinity < H;

  (b) if R is nonempty, h^H sum_{e in R} kappa_e < min_{e in R} kappa_e;

  (c) 0 < h < h_0.

All three hold for all sufficiently small h: X_h(p^0) = x(0), p^0 -> 0 as h -> 0+ (because log h -> -infinity), and h^H -> 0 since H > 0.

Set

    F(x) = min_{lambda in Lambda} ( r_lambda . x + c_lambda(h) ),   (35)

     B_h = [h, h^{-1}]^d.   (36)

F is the finite min of finitely many affine functions, hence concave, continuous, piecewise affine.

Step A. Suppose lambda in Lambda is active at x(0) = X_h(p^0). Then r_lambda = 0.

Proof. If for some i the entry r_lambda,i is nonzero, so is r_lambda . e_i. Apply (32) with v = e_i:

    sgn(p_i^0) = sgn(r_lambda,i)    and    |p_i^0| >= H.

But (a) gives |p_i^0| <= || p^0 ||_infinity < H. Contradiction, so r_lambda = 0.

Step B. Fix one active label lambda_0 at x(0). By Step A, r_lambda_0 = 0, so its value at x(0) is c_lambda_0(h). Because it is active, that value equals F(x(0)). Write

    M := F(x(0)) = c_lambda_0(h).   (37)

Since F is a minimum, F(x) <= M everywhere.

Step C. The boundary of B_h carries no maximum. Suppose x in dB_h with F(x) = M. Write x = X_h(p) with p in [-1, 1]^d and || p ||_infinity = 1. By (37), r_lambda_0 . x + c_lambda_0(h) = 0 + M = M = F(x); so lambda_0 is active at x. Corollary (17) gives || p - r_lambda_0 ||_infinity < E(r_lambda_0) = E(0) = 1/2. Since r_lambda_0 = 0, || p ||_infinity < 1/2, contradicting || p ||_infinity = 1.

Step D. Define

    K = { x in R^d : F(x) = M }.   (38)

By Step B, F(x) <= M everywhere, so K = { x : F(x) >= M } = intersection over lambda of closed half-spaces { x : r_lambda . x + c_lambda(h) >= M }. K is a closed convex polyhedron. Properties:

    (i) x(0) in K because F(x(0)) = M.

    (ii) K meets no boundary dB_h by Step C.

    (iii) K is contained in int B_h = (h, h^{-1})^d. If z in K lay outside int B_h, the segment [x(0), z] would lie in K by convexity, starting in int B_h and ending outside int B_h, so it would cross dB_h at some w in K, contradicting (ii).

    (iv) K is compact (closed in B_h, bounded), convex, and contained in R_{>0}^d by (iii).

Step E. K is forward invariant. Precisely, for every z in K, the solution of (1) starting at z exists for all t >= 0 and stays in K.

Proof. Three enclosures.

E1 (nonnegativity of active derivatives). For every x in B_h and every active lambda in Lambda at x,

    r_lambda . f(x) >= 0.   (39)

Apply the cut-flux lemma with V = C, E = R, a(y) = r_lambda . y, w_y = x^y, kappa_e = rate of reaction e, and delta = h^H. The return-path hypothesis is weak reversibility of R. Condition (18) is (33). Condition (19) is (b) of the h-choice. Hence

    sum_{e: y->z} kappa_e x^y ( r_lambda . z - r_lambda . y ) >= 0.   (40)

But this sum equals r_lambda . f(x). Done.

E2 (nondecrease of F along the trajectory inside B_h). Let x(t) solve (1) from z in K, on [0, T_max). Any interval [a, b] with x([a, b]) subset B_h satisfies E1, so the derivative lemma applied to Phi = F and gamma = x gives (d/dtau)^+ F(x(t)) >= 0 on [a, b]. Hence F o x is nondecreasing on [a, b].

E3 (value at M). F(x(t)) <= M by (37), and F(x(0)) = F(z) = M because z in K. Nondecrease forces F(x(t)) = M for all t with x(t) in B_h.

E4 (first-exit argument). z in K and K subset int B_h by (iii), so x(0) = z is in int B_h. Define

    T = sup { t in [0, T_max) : x([0, t]) subset K }.   (41)

Since x(0) in K and K is closed, T > 0 and x([0, T]) subset K. Thus x([0, T]) subset int B_h, and by E3, F(x(t)) = M there. If T < T_max, then x(T) in K subset int B_h, so continuity of x gives delta > 0 with x([T, T + delta]) subset int B_h. E2 and E3 on [T, T + delta] then give F(x(t)) = M, so x([T, T + delta]) subset K, contradicting the definition of T. Hence T = T_max.

E5 (extension past finite times). Suppose T_max < infinity. The trajectory lies in K subset B_h, so f, a polynomial, is bounded on K: sup_K ||f|| < infinity. The trajectory is Lipschitz, hence has a limit z* in K as t -> T_max. By the standard extension theorem for ODEs, the solution extends past T_max, contradiction. So T_max = infinity.

Done. Step E gives global existence and invariance.

Completing the theorem. Fix x(0) in R_{>0}^d and construct K as above. K is a compact convex subset of R_{>0}^d containing x(0), contained in B_h = [h, h^{-1}]^d, and forward invariant with global solutions. Hence x(t) exists for all t >= 0 and stays in K. Taking eps = h in (0, 1),

    eps <= x_i(t) <= 1/eps    (t >= 0, i = 1, ..., d).   (42)

Both bounds hold. This proves (2) and both parts of the theorem.

Positive equilibrium in every positive class. The invariant polytope is convex, compact, contained in R_{>0}^d, and forward invariant. Inside the affine class x(0) + S it therefore encloses a stationary point, and this supplies a direct proof of the known existence theorem of Boros (2019).

Corollary. Let (C, R) be a weakly reversible mass-action network with positive constant rates, and let x(0) in R_{>0}^d. Then there is an equilibrium x* in (x(0) + S) intersect R_{>0}^d.

pdf: pdfrefs/boros-2019.pdf

Proof. Let K be the invariant polytope of the previous paragraph and put

    Q = K intersect (x(0) + S).   (43)

The set Q is nonempty (it contains x(0)), compact, convex, and forward invariant: if y is in Q then phi_t(y) is in K by Step E and in y + S = x(0) + S by the definition of S.

For each integer n >= 1, the map phi_{1/n} : Q -> Q is continuous and Q is compact and convex, so Brouwer's fixed point theorem gives a fixed point z_n in Q:

    phi_{1/n}(z_n) = z_n.   (44)

The integral identity phi_s(z) - z = integral_0^s f(phi_r(z)) dr gives, at s = 1/n and z = z_n,

    0 = phi_{1/n}(z_n) - z_n = integral_0^{1/n} f(phi_r(z_n)) dr.   (45)

Write L = sup_Q ||f|| < infinity. The flow satisfies ||phi_r(z_n) - z_n|| <= L r <= L/n for 0 <= r <= 1/n. Uniform continuity of f on Q therefore gives

    ||f(z_n)|| = || n integral_0^{1/n} ( f(phi_r(z_n)) - f(z_n)) dr ||
              <= sup_{0 <= r <= 1/n} || f(phi_r(z_n)) - f(z_n) ||  ->  0.   (46)

A subsequence of z_n converges to some x* in Q (compact), and continuity of f gives f(x*) = 0. Since x* is in Q, it is in K and therefore in R_{>0}^d. Done.

Complex-balanced systems converge. A mass-action system is complex-balanced if there is a positive state x* at which incoming and outgoing fluxes agree at every complex:

    sum_{z: z->y in R} kappa_{z->y} (x*)^z = sum_{w: y->w in R} kappa_{y->w} (x*)^y    for every y in C.   (47)

This is stronger than f(x*) = 0 and implies weak reversibility of the reaction graph. Horn and Jackson (1972) proved that with complex-balanced positive rates and any positive stoichiometric class, there is exactly one complex-balanced state x*, and it is a positive equilibrium.

pdf: https://doi.org/10.1007/BF00251225

The Horn-Jackson entropy V : R_{>0}^d -> R, defined by

    V(x) = sum_{i=1}^d ( x_i log(x_i / x*_i) - x_i + x*_i ),   (48)

is smooth and nonnegative on R_{>0}^d, vanishes only at x*, and satisfies

    dV(x(t)) / dt <= 0   (49)

along any forward trajectory of (1), with equality only at complex-balanced states. These are the standard properties established by Horn and Jackson (1972) and Horn (1972), used by Feinberg (1987) in the deficiency-zero analysis.

Corollary. Under complex balance, for every x(0) in R_{>0}^d the trajectory x(t) converges to the unique positive equilibrium x* in (x(0) + S) intersect R_{>0}^d.

Proof. Complex balance implies weak reversibility, so Theorem (2) applies and x(t) stays in the compact positive set K. The trajectory therefore lies in Q = K intersect (x(0) + S). On Q, V is continuous, bounded below, and nonincreasing by (49). LaSalle's invariance principle identifies the omega-limit set of x(t) with a subset of the largest invariant set inside {x in Q : dV/dt = 0}. By (49) that set is {x*}, so x(t) -> x*.

Every weakly reversible deficiency-zero network with positive rates is complex-balanced, so its positive trajectories converge to a positive equilibrium. The deficiency is |C| minus the number of linkage classes minus the real dimension of S; the deficiency-zero convergence theorem is Feinberg (1987).

pdf: https://doi.org/10.1016/0009-2509(87)80096-5