When Enzymes Hit Their Limits — A Rigorous Stochastic Proof of the Total Quasi-Steady State Approximation
Two mathematicians have delivered the first fully rigorous probabilistic derivation of the total QSSA for the Michaelis-Menten enzyme model, settling a long-standing gap in the mathematical foundations of systems biology.
Deep inside living cells, enzymes grab substrate molecules and turn them into products millions of times per second. The mathematical story of how that process unfolds at scale dates back to Leonor Michaelis and Maud Menten in 1913. But a rigorous stochastic foundation for the most accurate version of their simplifying approximation has remained unproven for over a century. A new paper by Arnab Ganguly of Louisiana State University and Wasiur R. KhudaBukhsh of the University of Nottingham closes that gap entirely, deriving the total quasi-steady state approximation from first principles as both a Functional Law of Large Numbers and a Functional Central Limit Theorem.
Why Enzyme Reactions Need Approximations at All
The Michaelis-Menten reaction system describes one of the most fundamental processes in biochemistry. A substrate molecule S binds with an enzyme E to form a complex C. That complex either breaks apart to release S and E again, or it converts into a product P while releasing the enzyme back into circulation. In symbolic form this looks deceptively simple. The substrate and enzyme bind at rate \(\kappa_1\), the complex dissociates at rate \(\kappa_{-1}\), and the product formation step happens at rate \(\kappa_2\).
In reality, a single cell can contain thousands of enzyme molecules and millions of substrate molecules all reacting simultaneously. Tracking every individual molecule exactly is computationally impossible. So biochemists and mathematicians use approximations. The idea is that the complex C reaches a kind of instantaneous steady state much faster than the substrate changes, allowing the system to be reduced to a single equation for the substrate alone.
That shortcut is called a quasi-steady state approximation, or QSSA. The standard version, known as sQSSA, works well when the enzyme concentration is low compared to the substrate. But many real biological situations violate that assumption. Signaling cascades, gene regulatory networks, and drug-target interactions all involve enzymes and substrates at comparable concentrations. For those cases, a different approximation called the total QSSA, or tQSSA, has been known to perform better. It has been argued to be valid over a much wider range of parameter values than the sQSSA.
The sQSSA has had a rigorous stochastic derivation since at least 2006. The tQSSA has been used extensively in practice and justified heuristically, but a mathematically rigorous derivation from the stochastic model had never been carried out. Ganguly and KhudaBukhsh provide exactly that, along with a precise characterization of the fluctuations around the approximation.
What Makes the tQSSA Different
The key move in the total QSSA is a change of variables. Instead of tracking the substrate S directly, you introduce the combined quantity \(T = S + C\), which is the total amount of substrate in any form. This is a more natural variable because it is conserved at a higher level: whatever substrate is not free is locked up in the complex, and vice versa.
Setting the time derivative of C approximately equal to zero in the equations for T and C leads to a quadratic equation for the steady-state value of C in terms of T. The physically meaningful root of that quadratic gives an expression for C as a function of T and the enzyme total \(e_0\). Plugging this back in yields the tQSSA equation for how T evolves over time.
This expression looks more complicated than the sQSSA equation, but it works across a far broader range of conditions. The question that Ganguly and KhudaBukhsh set out to answer is whether this heuristic reduction is actually justified when the underlying process is stochastic rather than deterministic. Their answer is yes, and they prove it with full mathematical rigor.
The Stochastic Model — Poisson Random Measures and Scaling
The starting point of the rigorous analysis is a stochastic model of the MM system. Rather than writing ordinary differential equations, Ganguly and KhudaBukhsh model the species copy numbers as a pure jump Markov process and express the dynamics as stochastic differential equations driven by Poisson Random Measures. This is a natural and well-established framework for modeling chemical reactions at the individual molecule level.
The key technical device is a scaling regime. You introduce a large parameter n and rescale both the species abundances and the reaction rates. The scaling used for the tQSSA sets all species at the same abundance level (all of order n), makes the binding and product formation reactions slow compared to the unbinding reaction, and does not alter time. In this regime the unbinding reaction is fast at order \(n^2\), while binding and product formation operate at order n and order 1 respectively.
Under this scaling the process \(Z_V = Z_S + Z_C\) plays the role of the slow variable, analogous to the total substrate T in the deterministic picture. The complex \(Z_C\) is the fast variable, making rapid jumps at rate of order \(n^2\). The goal is to show that as n grows large, the slow process converges to the deterministic tQSSA trajectory, with the fluctuations around that trajectory being Gaussian and satisfying a stochastic differential equation.
The Technical Heart — Why This Problem Is Genuinely Hard
One might expect that deriving a stochastic averaging result for such a simple-looking model would be routine. It is not. The paper’s introduction is careful to explain exactly what makes the standard tools fail here, and the explanation is illuminating.
First, the fast process Z_C does not have a unique invariant distribution. The generator of the limiting fast ODE has two equilibrium points for any fixed value of the slow variable. One of those equilibria is stable and one is unstable, but this means there are infinitely many invariant probability measures: any convex combination of Dirac masses at the two equilibrium points is invariant. Standard stochastic averaging theorems typically require a unique invariant distribution. Here that assumption fails fundamentally.
Second, the fast and slow components are not cleanly separated in the usual sense. In standard multiscale reaction network frameworks, the fast dynamics are typically modeled by a finite-state Markov chain. Here the fast process takes values in a continuous state space and the coupling between fast and slow components is intricate.
The resolution of the first difficulty is one of the most elegant pieces of the paper. Even though there are infinitely many invariant distributions for the fast generator, the only one supported on the physically accessible region of state space — the interval from zero to the current value of the slow variable — is the Dirac mass at the stable equilibrium. The physically meaningful constraint on the state space is what selects the right invariant distribution, not uniqueness in the abstract sense.
“A rigorous mathematical derivation of the tQSSA from the stochastic perspective is notably absent in the literature despite its practical importance.” — Ganguly & KhudaBukhsh · J. Math. Anal. Appl. 561 (2026)
Theorem One — The Functional Law of Large Numbers
The first main result establishes convergence of the slow process to the deterministic tQSSA trajectory. To state it precisely, you need the notion of relative compactness for a sequence of stochastic processes indexed by n. The proof proceeds in three steps.
The first step establishes that the sequence of rescaled processes and their associated occupation measures is relatively compact. The occupation measure \(\Gamma^{(n)}\) captures how much time the pair \((Z_C, Z_V)\) spends in each region of state space up to time t. Even though the fast process itself is not relatively compact because of its rapid jumps, the sequence of occupation measures is well-behaved and forms a compact family.
The second step identifies what any limit point of the occupation measures must look like. For any test function g and any limit measure \(\Gamma\), the integral of the fast generator against \(\Gamma\) must equal zero. This is a form of the invariance equation. Combined with the support constraint mentioned above, Lemma 4.2 shows that the conditional distribution of \(Z_C\) given the slow component must be a Dirac mass at the stable equilibrium of the fast ODE.
The third and final step uses this identification of the limiting occupation measure to show that the slow process itself must solve the tQSSA differential equation. The key convergence argument uses the martingale decomposition of the slow process combined with the Lenglart-Rebolledo inequality to show that the martingale part vanishes in the limit.
Assume the initial data converge appropriately. Then as \(n \to \infty\), the rescaled slow process \(Z_V^{(n)}\) converges weakly to a limit \(Z_V\) that solves the random ODE
$$\frac{d}{dt} Z_V = -\kappa_2 \, z^-_{C,\star}(K_2, Z_V),$$
where \(z^-_{C,\star}(K_2, Z_V)\) is the stable equilibrium of the fast ODE at the current slow state. If the initial conditions are deterministic, the limit is deterministic and the convergence holds in probability in the uniform metric on \([0,T]\).
The formula for \(z^-_{C,\star}(K_2, Z_V)\) is exactly the quadratic root expression that appears in the deterministic tQSSA. This confirms that the tQSSA trajectory is not just a heuristic but the genuine large-population limit of the stochastic process.
Theorem Two — The Functional Central Limit Theorem
Knowing that the process converges to the tQSSA trajectory is valuable, but it leaves open how large the deviations from that trajectory typically are. The second main result answers this question by characterizing the fluctuations as the solution of an explicit stochastic differential equation.
Define the fluctuation process \(U_V^{(n)} = \sqrt{n}(Z_V^{(n)} – Z_V)\). This is the deviation from the tQSSA trajectory scaled by \(\sqrt{n}\) so that it remains of order one. Under appropriate assumptions on the initial data and the scaling of the enzyme total, the paper proves that \(U_V^{(n)}\) converges weakly to a limiting process \(U_V\) satisfying an Ito stochastic differential equation.
Under suitable moment conditions, the fluctuation process \(U_V^{(n)}\) converges weakly in \(D([0,T], \mathbb{R})\) to a limiting process \(U_V\) satisfying the SDE
$$dU_V(t) = \kappa_2 \, \partial z^-_{C,\star}(K_2, Z_V(t)) \, U_V(t) \, dt + \sqrt{\kappa_2 z^-_{C,\star}(K_2, Z_V(t))} \, dW(t) + \text{(drift from fluctuating enzyme total)},$$
where \(W\) is a standard Brownian motion. The diffusion coefficient at each time t equals the square root of the stable equilibrium complex concentration times the product formation rate.
This SDE has a clear structure. The first term is a linear restoring force that pulls fluctuations back toward zero at a rate determined by how sensitively the stable complex concentration depends on the slow variable. The second term is a Brownian noise term whose amplitude depends on the current stable equilibrium. The third term accounts for fluctuations in the enzyme total and contributes a deterministic drift to the limiting equation.
The Functional CLT says that at large but finite n, deviations of the stochastic process from the tQSSA trajectory are approximately Gaussian with a variance that evolves according to a known SDE. This has direct implications for statistical inference: it provides the theoretical foundation for deriving consistent estimators of reaction rate constants from experimental data on the reduced model.
The Poisson Equation — The Engine Behind the Fluctuation Analysis
The proof of the FCLT relies on a key auxiliary result about the solution to a Poisson equation for the fast generator. For each fixed value of the slow variable \(z_V\), you need to find a function F such that applying the fast generator to F gives the negative of the deviation of the fast process from its steady state.
The Poisson equation for the fast generator \(\mathcal{B}_{K_2, z_V}\) takes the form
Because the fast generator is the infinitesimal generator of a one-dimensional ODE with two roots, this equation can be solved exactly. The solution involves a logarithm of the distance from the fast state to the unstable equilibrium, normalized by the gap between the two equilibria. The explicit formula allows careful estimation of the derivatives of F with respect to both the fast and slow variables, which are essential inputs to the FCLT proof.
The role of the Poisson equation in fluctuation theory is analogous to the role of the averaging lemma in the law of large numbers proof. Just as the occupation measure argument pins down the leading-order behavior of the process, the Poisson equation argument pins down the fluctuations around that leading order. The two results together provide a complete asymptotic picture of the stochastic tQSSA.
Layer 1 — The Limit
The FLLN (Theorem 4.1) establishes that the stochastic process converges to the deterministic tQSSA ODE trajectory as molecule numbers grow. This is the zeroth-order picture.
Layer 2 — The Fluctuations
The FCLT (Theorem 5.1) characterizes the Gaussian fluctuations around the tQSSA trajectory at the next order in \(1/\sqrt{n}\). This gives the first-order correction.
Layer 3 — The Applications
Together, the two theorems justify accelerated simulation algorithms and statistical inference for biochemical networks that use the tQSSA as a model reduction tool.
Why This Matters for Systems Biology and Drug Development
The practical stakes of this work extend well beyond pure mathematics. The Michaelis-Menten framework is the foundational model of enzyme kinetics used in pharmacology, systems biology, synthetic biology, and biochemical engineering. Virtually every model of drug action at the molecular level involves some variant of the MM equations. The rate constants \(\kappa_1\), \(\kappa_{-1}\), and \(\kappa_2\) are quantities that experimentalists measure and that determine how quickly a drug is metabolized or how efficiently an enzyme functions in a metabolic pathway.
When scientists fit these rate constants to experimental data, they typically use reduced models like the tQSSA rather than the full system. The validity of statistical inference from such a procedure depends on understanding the relationship between the reduced model and the true stochastic process. The FCLT established in this paper is precisely the tool needed to make that relationship rigorous.
The paper notes, for example, that establishing the consistency of estimators of rate constants obtained from the reduced model requires a mathematically rigorous derivation of the tQSSA. Without Theorem 5.1, you cannot formally prove that your estimated rate constants will converge to the true values as you collect more data. With it, you can.
The same principle applies to simulation algorithms. Modern methods for simulating chemical reaction networks often split the network into fast and slow subsystems and apply different simulation strategies to each. The tQSSA provides the slow-subsystem equation that replaces the fast dynamics. Knowing that this replacement is mathematically valid, and knowing the size of the errors it introduces, allows developers of simulation software to provide rigorous error bounds on their algorithms.
The Challenge of Non-Unique Invariant Distributions
One of the most technically interesting contributions of this paper is the way it handles the non-uniqueness of the invariant distribution for the fast generator. This is not a minor technicality. Standard references on stochastic averaging, including the influential work of Kang and Kurtz, explicitly assume that the fast process has a unique invariant distribution. When this assumption fails, the general theory does not apply and a new argument must be constructed from scratch.
The key insight is geometric and physical rather than purely analytical. The fast ODE for the complex concentration \(z_C\) has two equilibrium points given by the roots of a quadratic. The larger root, the unstable equilibrium, always lies strictly above the maximum possible value of \(z_C\) given the conservation laws. The complex concentration can never reach the unstable equilibrium under physical initial conditions. Therefore, the part of the state space that the process actually visits contains only one equilibrium point of the fast ODE.
This observation, formalized as Lemma 4.2, says that the only invariant distribution of the fast generator supported on the physically relevant interval is the Dirac mass at the stable equilibrium. Once you restrict attention to the correct part of state space, the uniqueness that the general theory requires is recovered. The paper’s proof of the FLLN hinges on this argument in its third step.
Let \(\mu\) be a probability measure supported on \([0, z_V]\) that is invariant for the fast generator \(\mathcal{B}_{K,z_V}\). Then \(\mu\) must equal the Dirac mass at the stable equilibrium \(z^-_{C,\star}(K, z_V)\). The key fact is that the unstable equilibrium lies strictly above \(z_V\), so it is outside the physical region, and the only root of the fast ODE in \([0, z_V]\) is the stable one.
Relation to Earlier Work on Stochastic Enzyme Kinetics
The stochastic derivation of the sQSSA was carried out by Ball, Kurtz, Popovic, and Rempala in their landmark 2006 paper in the Annals of Applied Probability. That work established a rigorous framework for multiscale approximations using a generator-based approach and martingale problems. The present paper takes a more direct probabilistic route, working explicitly with the SDEs driven by Poisson Random Measures rather than passing through abstract generator theory.
This choice is deliberate and has several advantages. The PRM-based approach makes the connection between the stochastic model and the limiting ODE more transparent. The fluctuation analysis, in particular, benefits from having explicit formulas for the martingale terms rather than abstract existence results. The Burkholder-Davis-Gundy inequality and the Lenglart-Rebolledo inequality are used repeatedly to convert moment estimates on martingales into probability estimates on their suprema.
The paper also clarifies why earlier works on stochastic averaging for reaction networks cannot be directly applied to this model. The Kang-Kurtz framework in particular assumes that the reaction propensities for the fast reactions are of the same order as the abundance of the abundant species, a condition that does not hold in the tQSSA scaling regime. The present paper builds a custom averaging argument tailored to the specific structure of the MM system.
Numerical Confirmation and the Doob-Gillespie Algorithm
The paper includes a numerical comparison that illustrates the quality of the tQSSA approximation in practice. Using the Doob-Gillespie algorithm — the standard exact simulation method for discrete stochastic chemical reaction networks — the authors generate exact sample paths of the stochastic MM system at moderate values of n. These are compared visually to the tQSSA ODE trajectory.
With parameters \(\kappa_1 = 1\), \(\kappa_{-1} = 1\), \(\kappa_2 = 0.75\), \(K_2 = 0.1\), \(K_1 = 1.0\), and \(n = 1000\), the tQSSA trajectory tracks the stochastic simulations closely across the entire time horizon. The figure in the paper shows the two curves essentially overlapping, providing visual confirmation that the approximation is accurate at realistic parameter values long before the large-n limit is formally achieved.
This is not a substitute for the mathematical proof, but it is a useful sanity check. It confirms that the scaling regime chosen for the analysis is not so extreme as to be biologically unrealistic, and that the theorems are likely to be informative at the population sizes encountered in actual biological systems.
What Comes Next
The results of this paper open several directions for future research. One natural extension is to more complex enzyme reaction networks, where a single enzyme can bind multiple substrates or where several enzymes compete for the same substrate. The MM system is the simplest possible case, and understanding whether the same averaging phenomena persist in more complex networks is an important open question.
A second direction concerns statistical inference. Now that the FCLT is available, it becomes possible to derive explicit formulas for the asymptotic distributions of maximum likelihood estimators for the rate constants obtained from tQSSA-reduced data. This would provide practitioners with confidence intervals and hypothesis testing procedures that are grounded in the actual stochastic dynamics of the system.
A third direction involves large deviations. The FLLN and FCLT together describe the typical behavior and the typical fluctuations. What happens in the tail, when the process deviates far from the tQSSA trajectory, is governed by large deviation principles. Establishing a large deviation principle for the stochastic MM system in the tQSSA scaling regime would complete the asymptotic picture and is a natural sequel to the present work.
Read the Full Paper
Published in the Journal of Mathematical Analysis and Applications, Volume 561 (2026). Full proofs, lemmas, and the numerical simulation are available via the journal website.
A. Ganguly, W.R. KhudaBukhsh. Asymptotic analysis of the total quasi-steady state approximation for the Michaelis-Menten enzyme kinetic reactions. Journal of Mathematical Analysis and Applications, 561 (2026) 130551. https://doi.org/10.1016/j.jmaa.2026.130551
This article is an independent editorial analysis of a peer-reviewed open-access paper published under the CC BY 4.0 license. Mathematical statements paraphrase the original results; for complete proofs and precise formulations, consult the published paper. The paper was received 26 May 2025 and made available online 27 February 2026. Research by A. Ganguly was supported in part by NSF DMS-2246815 and the Simons Foundation. W.R. KhudaBukhsh was supported by the University of Nottingham International Research Collaboration Fund, the London Mathematical Society Scheme 4 grant 42360, and EPSRC grant EP/Y027795/1.
Selected References
- [1] K. Ball, T.G. Kurtz, L. Popovic, G.A. Rempala. Asymptotic analysis of multiscale approximations to reaction networks. Ann. Appl. Probab. 16 (2006) 1925–1961.
- [2] J. Eilertsen, S. Schnell, S. Walcher. The unreasonable effectiveness of the total quasi-steady state approximation, and its limitations. J. Theor. Biol. 583 (2024) 111770.
- [3] A. Ganguly, W.R. KhudaBukhsh. Functional limit theorems and parameter inference for multiscale stochastic models of enzyme kinetics. arXiv:2409.06565, 2024.
- [4] H.W. Kang, W.R. KhudaBukhsh, H. Koeppl, G.A. Rempala. Quasi-steady-state approximations derived from the stochastic model of enzyme kinetics. Bull. Math. Biol. 81 (2019) 1303–1336.
- [5] H.W. Kang, T.G. Kurtz. Separation of time-scales and model reduction for stochastic reaction networks. Ann. Appl. Probab. 23 (2013) 529–583.
- [6] H.W. Kang, T.G. Kurtz, L. Popovic. Central limit theorems and diffusion approximations for multiscale Markov chain models. Ann. Appl. Probab. 24 (2014) 721–759.
- [7] J.K. Kim, J.J. Tyson. Misuse of the Michaelis-Menten rate law for protein interaction networks and its remedy. PLoS Comput. Biol. 16 (2020) e1008258.
- [8] A.R. Tzafriri, E.R. Edelman. The total quasi-steady-state approximation is valid for reversible enzyme kinetics. J. Theor. Biol. 226 (2004) 303–313.
- [9] W. Whitt. Proofs of the martingale FCLT. Probab. Surv. 4 (2007) 268–302.
- [10] D.J. Wilkinson. Stochastic Modelling for Systems Biology. Chapman and Hall/CRC, 2018.
