Chemical reaction networks can only express ODEs driven by nonnegative polynomials, but the scalar ODE $y' = p(y)$ one wants to simulate allows $y$ to be signed and allows $p$ to contain negative monomials. Dual-railing is the standard fix: split $y = u - v$ with $u, v \ge 0$, each tracked by its own species. The subtle question is how to prevent $u$ and $v$ from drifting off to infinity together. A natural open problem asks whether a single constant $k$ can always suppress that drift, uniformly over initial conditions. This post closes the first nontrivial case.

The Constant-$k$ Boundedness Question

A recent paper by Haisler, Huang, Migunov, Mohammed, and Provence introduced selective dual-railing: given a scalar polynomial $p(y) \in \mathbb{Z}[y]$ with the property that $y' = p(y)$ is bounded on some interval $I \subseteq \mathbb{R}$ containing $0$, one builds a pair of ODEs in $(u, v)$ such that $u - v$ solves the original equation, $u$ and $v$ stay nonnegative, and a cross-cancellation term $-k \cdot u v$ keeps the pair itself bounded.

Their construction works, but it leaves a natural question open:

Open problem. Is there a constant $k > 0$ — independent of $y(0) \in I$ — such that the dual-railed system stays bounded on $[0, \infty)$?

The word “constant” is the whole point. If $k$ were allowed to depend on $y(0)$, the problem would be trivial: pick $k$ large enough for the initial data you happen to have. What the question demands is a single $k = k(p)$, determined by the polynomial alone, that handles every admissible starting value.

The simplest nontrivial instance is the scalar cubic

The scalar ODE $y' = 1 - y^3$, with $y(0) \in [0, 1]$, is globally bounded and converges to the fixed point $y = 1$. Its dual-railed version is

where the positive and negative monomials of $1 - (u - v)^3 = (1 + 3u^2 v + v^3) - (u^3 + 3 u v^2)$ have been routed into $u$’s and $v$’s rates respectively. Subtracting the two equations gives $u' - v' = 1 - (u - v)^3$, recovering $y' = 1 - y^3$ along $y = u - v$. The $-k u v$ term is the annihilator that has to earn its keep.

Zero-Initialization Non-Collapse

This post focuses on the boundary case $u(0) = v(0) = 0$. The corresponding scalar initial condition is $y(0) = 0$, at the edge of the admissible interval $[0, 1]$. Why single out zero initialization? Because the open problem wants a single $k$ that works across all admissible $y(0)$, and there is no canonical way to split a positive $y(0)$ between $u(0)$ and $v(0)$ — so the most symmetric case, $u(0) = v(0) = 0$, must at least not blow up on its own. Call this requirement zero-initialization non-collapse. It is the first necessary condition for the constant-$k$ statement.

Theorem. For $p(y) = 1 - y^3$ and every $k > 6$, the dual-railed system with $u(0) = v(0) = 0$ has a unique global solution on $[0, \infty)$, and $\max(u(t), v(t)) \le k/3$ for all $t \ge 0$.

The threshold $k = 6$ is not arbitrary. It is the positive real root of the cubic

which factors as

The double factor $(k + 3)^2$ at $k = -3$ is a separate saddle-node bifurcation in the unphysical branch; the one that matters here sits at $k = 6$.

The Proof in Three Moves

1. Symmetrize via $\sigma = u + v$

Let $\sigma := u + v$ and $y := u - v$. Two algebraic identities do most of the work:

• $u^3 + v^3 + 3 u v (u + v) = (u + v)^3 = \sigma^3$,
• $4 u v = \sigma^2 - y^2$.

Adding the two dual-railed ODEs collapses all the cubic monomials into $\sigma^3$, and the annihilator $-2 k \cdot u v$ becomes $-\tfrac{k}{2}(\sigma^2 - y^2)$. Thus

Meanwhile $y = u - v$ satisfies the original scalar ODE $y' = 1 - y^3$, $y(0) = 0$, and a one-dimensional barrier argument shows $y(t) \in [0, 1]$ for all $t$. So while analyzing $\sigma$ we may treat $y^2 \in [0, 1]$ as a bounded exogenous parameter.

2. The Saddle-Node at $k = 6$

For fixed $y \in [0, 1]$, regard $\sigma$ as the only dynamical variable:

The equilibria are the zeros of $Q_k(\,\cdot\,; y)$. The worst case is $y = 1$, where forcing is strongest:

Its discriminant, up to positive factors, is exactly $(k - 6)(k + 3)^2$.

• For $k < 6$, the discriminant is negative: no pair of positive equilibria exists, so $\sigma$ has no wall to keep it from running off to infinity.
• For $k = 6$, the discriminant vanishes with a double root at $\sigma = 2$ — a textbook saddle-node.
• For $k > 6$, two positive equilibria $0 < \sigma_- < \sigma_+ \le k/2$ appear. The smaller one, $\sigma_-(y)$, bounds a forward-invariant interval $[0, \sigma_-(y)]$ for $\sigma$.

3. A Strict Explicit Barrier at $\sigma = k/3$

The Lean formalization does not use $\sigma_-(1)$ directly (it is defined implicitly by a cubic and interacts poorly with constructive real analysis). Instead it uses the cruder explicit overestimate $\sigma = k/3$. The key inequality is

Plug in the sharpest allowed $y^2 \le 1$:

The same factorization resurfaces. For $k > 6$ the right-hand side is strictly negative, so at $\sigma = k/3$ the drift is strictly inward. Combined with $\sigma(0) = 0$ and the outward drift $Q_k(0; y) = 1 + \tfrac{k}{2} y^2 > 0$ at the lower boundary, this traps $\sigma(t) \in [0, k/3]$ for all $t \ge 0$ via a first-crossing argument: if $t_1$ were the first time $\sigma(t_1) = k/3$, then $\sigma'(t_1) < 0$ contradicts $t_1$ being a first crossing.

The rest is Picard–Lindelöf. Since $(u, v)$ lies in the compact box $\mathcal{B} := \{(u, v) \ge 0 : u + v \le k/3\}$ and the dual-railed vector field is polynomial (hence Lipschitz on compacts), global existence and uniqueness on $[0, \infty)$ follow.

Lean 4 Formalization

The entire proof has been mechanized in Lean 4 with Mathlib:

• File: Ripple/DualRail/ScalarCubic.lean, 2045 lines.
• Top-level theorem scalar_cubic_bounded at line 2013.
• Six sublemmas: scalar_cubic_nonneg (816), scalar_cubic_dual_rail_identity (838), scalar_cubic_original_bounded (1146), scalar_cubic_sigma_drift (1483), scalar_cubic_sigma_bound (1834), scalar_cubic_picard (1850).
• sorry count: $0$. Axioms beyond Mathlib: none.
• Commit: b88b18b.

The mechanization closely follows the paper proof. The trickiest step was not the algebraic identity — ring disposes of that — but the strict-barrier argument that $\sigma(t) \le k/3$ holds for all $t \ge 0$, which uses neither a Lyapunov function nor Grönwall but the first-crossing argument above plus ODE uniqueness. A $d$-dimensional Picard wrapper in Ripple.Core.ODEGlobal takes a forward-invariant box as input and produces the global solution directly, side-stepping the need to encode invariants as sublevel sets of a Lyapunov function.

Nonzero Initial Data: the Universal Ratio $K^\ast \approx 3.158$

The full statement of the open problem still needs $y(0) = U_0 \in (0, 1]$. Numerical experiments reveal a clean scaling law: starting from $u(0) = U_0$, $v(0) = 0$, the critical ratio for boundedness satisfies

Rescaling by $U_0$, the $+1$ forcing term drops out (it is three orders below the leading monomial), leaving a universal two-dimensional autonomous system

Since $\tilde{y}$ integrates explicitly, $\tilde{y}(\tau) = (1 + 2\tau)^{-1/2}$, one reduces $\tilde{\sigma}$ to a non-autonomous scalar cubic with explicit rational forcing,

$K^\ast$ is the value of $K$ for which the initial point $(1, 1)$ lies on the one-dimensional center manifold of the non-hyperbolic fixed point $(K/2, 0)$. Phrased as a shooting problem, $K^\ast$ is the eigenvalue of a boundary-value problem.

This constant appears to have no closed form. It differs from $\pi$ by $0.016$ and from $\sqrt{10}$ by $0.004$. An integer-relation search (PSLQ) up to degree $10$ with coefficients bounded by $10^{15}$ found no algebraic relation. If a closed form exists, it likely involves transcendental special-function values — possibly from the Painlevé hierarchy — but that remains speculation at present.

Three thresholds, two limits. It is worth separating three critical values of $k$ that arise for the scalar cubic:

1. Zero-initialization saddle-node, $k = 6$. This is where the cubic $Q_k(\sigma; 1)$ just barely develops a forward-invariant equilibrium; it is the threshold proved in this note for $y(0) = 0$.
2. Rescaled static saddle-node, $K = 3\sqrt{3} \approx 5.196$. After rescaling by $U_0$, if one treats $\tilde{y} \equiv 1$ as a frozen parameter, the corresponding cubic has a saddle-node exactly at $3\sqrt{3}$.
3. Rescaled dynamic eigenvalue, $K^\ast \approx 3.158$. This is the actual critical $K$ once one accounts for the fact that $\tilde{y}(\tau)$ decays like $(1 + 2\tau)^{-1/2}$. Because $\tilde{y}$ ultimately drops to zero, the problem is strictly easier than the static version, which is why $K^\ast < 3\sqrt{3}$.

The same stratification appears for the scalar quintic $p(y) = 1 - y^5$: the zero-init saddle-node (degree-$5$ critical polynomial, no clean factorization) sits at $k \approx 13.28$; the rescaled static saddle-node is $K = \sqrt{3125/27} \approx 10.76$; and the rescaled dynamic eigenvalue is numerically $K^\ast_{(5)} \approx 5.41$. The ratios $K_{\text{static}} / K^\ast_{\text{dynamic}}$ are $1.65$ for the cubic and $1.99$ for the quintic, suggesting a degree-dependent factor that might have a clean asymptotic description.

What Comes Next

with no factorization analogous to the cubic’s $(k-6)(k+3)^2$. A saddle-node threshold still exists, but its value is no longer rational. The longer-range goal is to write the annihilation rate $k$ as an explicit function of the polynomial $p$ and of $\sup_I |y|$. The present proof suggests $k = \Theta(\lVert p \rVert \cdot \sup|y|^{(\deg p - 1)/\deg p})$, though this is only a scaling guess at the moment.

Closing the full scalar-cubic case of the open problem will also require handling nonzero initial data, which in turn requires turning the empirical scaling law $K^\ast \approx 3.158$ into a rigorous asymptotic. Whether $K^\ast$ is algebraic is a separate research question in its own right.

The full research notes, with theorem statements, proofs, numerical details, and cross-references to the Lean file, are available as PDF.

Truth emerges more readily from error than from confusion.

— Francis Bacon, Novum Organum