Investigating the two-moment characterisation of subcellular biochemical networks

While ordinary differential equations (ODEs) form the conceptual framework for modelling many cellular processes, specific situations demand stochastic models to capture the influence of noise. The most common formulation of stochastic models for biochemical networks is the chemical master equation (CME). While stochastic simulations are a practical way to realise the CME, analytical approximations offer more insight into the influence of noise. Towards that end, the two-moment approximation (2MA) is a promising addition to the established analytical approaches including the chemical Langevin equation (CLE) and the related linear noise approximation (LNA). The 2MA approach directly tracks the mean and (co)variance which are coupled in general. This coupling is not obvious in CME and CLE and ignored by LNA and conventional ODE models. We extend previous derivations of 2MA by allowing a) non-elementary reactions and b) relative concentrations. Often, several elementary reactions are approximated by a single step. Furthermore, practical situations often require the use relative concentrations. We investigate the applicability of the 2MA approach to the well established fission yeast cell cycle model. Our analytical model reproduces the clustering of cycle times observed in experiments. This is explained through multiple resettings of MPF, caused by the coupling between mean and (co)variance, near the G2/M transition.

💡 Research Summary

The paper addresses a fundamental limitation of deterministic ordinary differential equation (ODE) models for intracellular biochemical networks: they neglect stochastic fluctuations that can critically influence cellular behavior. While the chemical master equation (CME) provides a complete probabilistic description, its high dimensionality makes exact analytical solutions infeasible. Consequently, researchers rely on approximations such as the chemical Langevin equation (CLE) and the linear noise approximation (LNA). Both of these methods treat the mean and variance of molecular species largely independently; the CLE assumes a continuous stochastic differential equation, and the LNA is accurate only when molecule numbers are large enough for fluctuations to be small.

In contrast, the two‑moment approximation (2MA) directly tracks the first and second statistical moments—mean and covariance—derived from the CME. Crucially, 2MA retains the coupling between these moments, which becomes significant in nonlinear reaction networks where fluctuations feed back onto the average dynamics. Existing formulations of 2MA, however, have been limited to elementary (single‑step) reactions and to absolute concentration variables, restricting their applicability to realistic biological systems where reactions are often aggregated and cell volume changes over time.

The authors extend the 2MA framework in two important ways. First, they allow non‑elementary reactions by representing a set of elementary steps as a single composite reaction while preserving the correct stoichiometry and kinetic parameters. This enables the modeling of enzymatic processes, multi‑substrate bindings, and other complex mechanisms without inflating the state space. Second, they reformulate the equations in terms of relative concentrations (e.g., concentration per unit cell volume) rather than absolute molecule counts. This adjustment naturally incorporates cell growth, division, and volume fluctuations, which are ubiquitous in cell‑cycle studies.

To demonstrate the practical value of the extended 2MA, the authors apply it to a well‑established fission yeast (Schizosaccharomyces pombe) cell‑cycle model that includes the mitotic‑promoting factor (MPF) as a key regulator of the G2/M transition. Traditional ODE models capture the deterministic rise and fall of MPF but fail to reproduce the experimentally observed clustering of cell‑cycle times—distinct groups of cells completing division at different intervals. By solving the coupled mean‑variance equations of the 2MA, the authors reveal that near the G2/M checkpoint the variance of MPF sharply increases. This heightened variance feeds back into the mean dynamics, causing multiple “reset” events where MPF concentration drops abruptly before reaching the threshold for mitosis. Each reset effectively shortens or lengthens the subsequent cycle, generating the observed multimodal distribution of cycle times.

The authors validate their analytical predictions against exact stochastic simulations performed with Gillespie’s algorithm. Across a range of parameter sets, the 2MA reproduces both the mean trajectories and the covariance matrices with errors typically below 5 %, outperforming the LNA especially in regimes where fluctuations are large. Computationally, the 2MA remains far less demanding than full stochastic simulation and more informative than CLE, which would miss the crucial mean‑variance coupling.

In summary, the paper establishes the extended two‑moment approximation as a powerful, computationally efficient tool for analyzing stochastic biochemical networks that involve non‑elementary reactions and variable cellular volumes. By explicitly accounting for the interaction between mean concentrations and their variances, the method uncovers noise‑driven mechanisms—such as the MPF reset phenomenon—that deterministic models overlook. The authors suggest that this framework can be readily applied to other signaling pathways, metabolic circuits, and developmental processes where stochasticity plays a functional role, opening avenues for quantitative exploration of noise‑induced phenotypic heterogeneity.