PaperlayerView exact revision on arXiv

Source-linked AI summary

Approximation and inference methods for stochastic biochemical kinetics - a tutorial review

David Schnoerr, Guido Sanguinetti, Ramon Grima

arXiv:1608.06582v2q-bio.QMcond-mat.stat-mechphysics.bio-phq-bio.MNstat.ML

TL;DR

Stochastic biochemical systems are difficult to analyze because the Chemical Master Equation usually lacks analytic solutions and exact simulations are computationally expensive. The review synthesizes modelling, approximation, and Bayesian inference methods, compares approximations numerically, and discusses their limitations and applicability.

  • Problem

    The Chemical Master Equation generally lacks analytic solutions, while exact stochastic simulation is computationally expensive, complicating analysis and inference for stochastic biochemical systems.

  • Method

    The review introduces deterministic and stochastic chemical-kinetics models, surveys CME approximations and Bayesian inference methods, and compares selected approximations in a numerical case study.

  • Results

    The case study found the complex-valued CLE most accurate for moments, while SSE-1 and 2MA outperformed the LNA; CLE accuracy depended strongly on implementation and enzyme number.

  • Takeaways & Limitations

    The review provides a self-contained introduction and overview of approximation and inference methods for stochastic chemical kinetics.

  • Takeaways & Limitations

    Forward-backward inference requires transition probabilities that generally require solving the CME, and open systems may need state-space truncation that introduces difficult-to-quantify bias.

Abstract

from arXiv · show

Stochastic fluctuations of molecule numbers are ubiquitous in biological systems. Important examples include gene expression and enzymatic processes in living cells. Such systems are typically modelled as chemical reaction networks whose dynamics are governed by the Chemical Master Equation. Despite its simple structure, no analytic solutions to the Chemical Master Equation are known for most systems. Moreover, stochastic simulations are computationally expensive, making systematic analysis and statistical inference a challenging task. Consequently, significant effort has been spent in recent decades on the development of efficient approximation and inference methods. This article gives an introduction to basic modelling concepts as well as an overview of state of the art methods. First, we motivate and introduce deterministic and stochastic methods for modelling chemical networks, and give an overview of simulation and exact solution methods. Next, we discuss several approximation methods, including the chemical Langevin equation, the system size expansion, moment closure approximations, time-scale separation approximations and hybrid methods. We discuss their various properties and review recent advances and remaining challenges for these methods. We present a comparison of several of these methods by means of a numerical case study and highlight some of their respective advantages and disadvantages. Finally, we discuss the problem of inference from experimental data in the Bayesian framework and review recent methods developed the literature. In summary, this review gives a self-contained introduction to modelling, approximations and inference methods for stochastic chemical kinetics.

1 Introduction

Stochastic chemical kinetics is difficult to analyze because biological systems fluctuate, the CME usually lacks analytic solutions, and exact simulation is expensive. The review addresses these gaps with an accessible, self-contained introduction and comparison of approximation and inference methods.

  • Motivation: Biological molecule-number fluctuations make stochastic systems harder to analyze than deterministic counterparts.These fluctuations can produce different behaviors in genetically identical cells.
  • Motivation: The CME generally has no analytic solution, while exact stochastic simulation becomes computationally expensive for larger systems.The stochastic simulation algorithm simulates every chemical reaction event.
  • Review gap: Existing reviews mainly emphasize simulation or require substantial prior mathematical knowledge, leaving an accessible overview and comparison of approximation and inference methods unavailable.The article is intended for non-experts new to stochastic modelling.
  • Scope: The review introduces deterministic and stochastic modelling, CME solution and simulation methods, and several approximation approaches.Covered approximations include the chemical Langevin equation, system size expansion, moment closure, time-scale separation, and hybrid methods.
  • Scope: A numerical case study compares the approximation methods, while a Bayesian section reviews inference methods for experimental data.The review presents both modelling approximations and statistical inference in one framework.

2 Stochasticity in biological systems

Stochasticity arises throughout gene expression and produces intrinsic fluctuations within cells as well as extrinsic differences between cells. Experiments show that these fluctuations can substantially affect cellular behavior, motivating mathematical models that connect measured noise to underlying processes.

  • Sources of stochasticity: Gene expression involves stochastic transcription, translation, and degradation because molecular encounters and reactions occur randomly.RNA polymerase binding and movement, ribosome–mRNA encounters, and degradation are described as stochastic processes.
  • Gene expression: In eukaryotic cells, transcription occurs in the nucleus and mature mRNA must undergo processing and diffusion before cytosolic translation.This adds steps compared with the corresponding prokaryotic gene-expression mechanism.
  • Experimental measurements: Measurements show strong temporal protein fluctuations and substantial contributions from both intrinsic and extrinsic noise.The observations come from fluorescent time series and dual-reporter measurements in Escherichia coli.
  • Sources of stochasticity: Intrinsic noise denotes temporal molecule-number fluctuations caused by stochastic chemical processes, whereas extrinsic noise reflects physiological or environmental differences between cells.Examples of extrinsic variation include differing RNA polymerase or ribosome numbers and nutrient concentrations.
  • Experimental measurements: The dual reporter technique separates intrinsic and extrinsic noise by measuring two identically regulated fluorescent genes exposed to shared external effects.The reporters encode distinguishable fluorescent proteins under identical promoters.
  • Functional consequences: Stochastic fluctuations can influence cell function, including stochastic cell-fate decisions that produce different states among genetically identical cells.Such diversity is described as potentially beneficial in fluctuating environments.

3 Stochastic chemical kinetics

Chemical reaction networks can be described deterministically with rate equations or stochastically with discrete molecule counts. The section develops the reaction-network representation, rate functions, stoichiometry, and a negative-feedback gene example.

  • Reaction-network representation: Chemical reaction networks represent complicated biological mechanisms using effective chemical reaction events, such as modelling transcription or translation as single reactions.This provides a simplified mathematical representation of multistep biological processes.
  • Reaction-network representation: Reaction order is determined by the number of reactant molecules; unimolecular and bimolecular reactions have orders one and two, respectively.A system is linear when all reactions have order at most one, and open when molecules can be generated.
  • Deterministic modelling: The Law of Mass Action makes reaction rates proportional to products of reactant concentrations in the macroscopic description.The stoichiometric matrix combines reaction-specific net changes with macroscopic rates to determine concentration dynamics.
  • Negative-feedback example: The example gene network uses negative autoregulation: protein production occurs in the unbound state, while promoter binding switches the gene off.Protein degradation and promoter binding are represented as additional reactions.
  • Deterministic modelling: Deterministic rate equations describe mean concentrations through ordinary differential equations and ignore molecule-number fluctuations.They are relatively straightforward to analyze and are typically accurate for systems with large molecule numbers.
  • Negative-feedback example: A conservation law for total gene number reduces the example from three species variables to a lower-dimensional system used throughout the article.The conserved quantity is the sum of genes in the on and off states.

3.2 Stochastic methods

Stochastic chemical kinetics retains discrete molecule numbers and models reactions as a continuous-time Markov jump process under well-mixed, dilute conditions. The CME is usually analytically intractable, motivating simulation, moment equations, and approximations.

  • Stochastic description: Under well-mixed and dilute conditions, the stochastic state is the discrete molecule-count vector and spatial locations need not be modelled.Reaction probabilities over an infinitesimal interval are determined by propensity functions.
  • Chemical Master Equation: Mass-action propensity functions depend on available reactant combinations and provide the microscopic counterpart of macroscopic mass-action rate equations.Zeroth-, first-, and second-order reactions have propensity forms based on volume, molecule counts, and reactant combinations.
  • Chemical Master Equation: The CME is an infinite coupled system of linear ordinary differential equations for state probabilities and generally has no analytic solution.Even the relatively simple gene example lacks a known time-dependent CME solution.
  • Moment equations: Moment equations are derived from the CME by multiplying by products of molecule counts and summing over all states.The resulting equations describe time evolution of moments such as means and variances.
  • Moment equations: For linear systems, moment equations close exactly because moments of order m depend only on moments of order m or lower.Finite sets of linear ordinary differential equations then yield exact moments up to a chosen order.
  • Moment equations: For nonlinear systems, higher-order reactions couple each moment equation to higher moments, producing an infinite hierarchy that motivates moment-closure approximations.In the gene example, higher-moment terms are proportional to the bimolecular reaction rate constant k2.

3.3 Stochastic simulations

The stochastic simulation algorithm (SSA) generates exact sample paths of chemical reaction networks by explicitly simulating reaction events, but its computational cost limits applicability. Standard SSA assumptions also fail for time-varying propensities unless approximations or modified integration strategies are used.

  • The SSA is a Monte Carlo method that simulates exact sample paths of the stochastic process described by the CME.
  • The direct SSA updates reaction time and system state by sampling the next event and then selecting its reaction type.
  • SSA computational cost rises because every reaction event is simulated explicitly, especially with large molecule-number fluctuations or many reactions per unit time.
  • Standard SSA relies on constant propensities between events, producing exponentially distributed inter-reaction times.
  • With stochastic time-varying propensities, inter-reaction times are not exponentially distributed, so the standard algorithms are no longer valid.
  • Assuming propensities remain constant is an approximation that can become inaccurate when extrinsic fluctuations are strong between reaction events.

3.4 Exact results

Exact CME solutions exist for restricted reaction systems and selected special cases, but most practically relevant systems lack analytic solutions. Finite-state matrix methods are principled yet often computationally expensive, reinforcing the need for approximation methods.

  • For finite state spaces, the CME can be written as a finite matrix ODE and solved through matrix exponentiation in principle.
  • Even finite-state CME solutions can be impractical because matrix dimension is often large and exponentiation computationally expensive.
  • Moment equations close at finite order for linear reaction systems, whereas nonlinear systems couple moments to higher orders.
  • Analytic time-dependent solutions are available for some linear systems with at most one product molecule per reaction, but this excludes biologically common reactions such as A → A + B.
  • Non-linear systems: Detailed balance requires each reaction's flow in every state to equal the flow of its corresponding reverse reaction; under mass-action kinetics, this is equivalent to deterministic detailed balance.
  • Non-linear systems: Weakly reversible, complex-balanced mass-action systems possess steady-state CME solutions given by product-Poisson distributions adjusted for conservation laws.
  • Most practical reaction systems satisfy neither linearity nor detailed or complex balance, and known exact examples generally involve few species or steady states.

4 Approximation methods

Because analytic CME solutions and exact stochastic simulations are often unavailable or computationally infeasible, approximation methods are developed to make stochastic chemical-kinetics analysis tractable.

  • Approximation methods address the lack of analytic CME solutions and the computational infeasibility of exact stochastic simulation for many practical systems.
  • The review introduces CFPE/CLE approximations, system size expansion, and further approximation families for stochastic chemical kinetics.

4.1 The chemical Langevin equation

The chemical Langevin equation (CLE) is a diffusion approximation to the Chemical Master Equation (CME), equivalent to a chemical Fokker–Planck equation and often more efficient to simulate. Its accuracy improves with molecule abundance, but real-valued formulations face boundary problems that motivate complex-valued and hybrid alternatives.

  • Definition and derivation: The CLE and chemical Fokker–Planck equation form equivalent diffusion approximations of the CME, with the CLE generating process realisations.The CLE is an Itô stochastic differential equation, while the chemical Fokker–Planck equation describes the corresponding process distribution.
  • Definition and derivation: Different factorizations of the diffusion matrix produce different CLE representations, although the most commonly used representation is only one such choice.The diffusion matrix remains unchanged under corresponding sign changes in the noise-factor columns.
  • Stochastic simulations: CLE simulations scale with the number of species rather than reaction events, often making them more efficient than CME simulations when molecule counts are not too small.Euler–Maruyama simulation introduces a time-step trade-off: smaller dt improves approximation accuracy but increases computational cost.
  • Properties and recent developments: The CLE is typically accurate at moderate or high molecule numbers and becomes exact in the thermodynamic limit.For linear reaction systems, its moments up to order two agree exactly with those of the CME; for general systems, higher-order coupling can produce differences.
  • Properties and recent developments: Real-valued CLE formulations can become ill-defined when square roots receive negative arguments, and ad hoc boundary modifications may be inaccurate for nonlinear systems.The problem is independent of the factorization of the noise matrix and is expected for the majority of reaction systems.
  • Properties and recent developments: The complex-valued CLE restores second-order moment exactness for linear systems and can be highly accurate for some nonlinear systems, but its results must be projected to real space.In adaptive hybrid methods, the complex CLE reduces automatically to the real-valued CLE when negative-concentration probabilities are avoided.
  • Properties and recent developments: The CLE can capture multimodality that persists at large volumes and may also reproduce noise-induced multimodality, though not for every reaction system.Tensor-based methods have also been proposed for direct numerical solution of the chemical Fokker–Planck equation, supporting sensitivity and bifurcation analysis.

4.2 The system size expansion

The system size expansion approximates stochastic chemical kinetics by separating deterministic dynamics from fluctuations and expanding in inverse system size. Its lowest-order form yields the linear noise approximation, while higher orders provide moment corrections efficiently but have important scope and validity limitations.

  • 4.2.1 Derivation: The system size expansion separates molecule numbers into a deterministic rate-equation component and fluctuations about the deterministic mean.It introduces fluctuation variables scaled by the inverse square root of system volume.
  • 4.2.2 The linear noise approximation: Truncating at zeroth order gives the linear noise approximation, whose Fokker–Planck equation has linear drift, constant diffusion, and a multivariate normal solution under suitable initial conditions.The resulting molecule-number mean follows the rate equations, while covariance dynamics are obtained from ordinary differential equations.
  • 4.2.3 Higher order corrections: Higher-order truncations produce ordinary differential equations for moments, with order Ω^-1 corrections to the mean and order Ω^-2 corrections to the covariance.These corrections are termed effective mesoscopic rate equations and Inverse Omega Square, respectively.
  • 4.2.4 Properties and recent developments: For moments rather than higher-order distributions, the system size expansion is generally significantly more efficient than stochastic simulation of the CME or CLE.It reduces the calculation to finite sets of ordinary differential equations without ensemble averaging.
  • 4.2.4 Properties and recent developments: The expansion cannot generally handle deterministically multistable systems because it is centered on the deterministic mean, except for short-time analyses.Higher-order PDE truncations can also lack probabilistic interpretation because their solutions may not be positive-definite.

4.3 Moment closure approximations

Moment closure approximations replace higher-order moments with functions of lower-order moments, converting the CME’s infinite moment hierarchy into finite ODE systems. They are efficient and easy to implement, but accuracy and physical validity depend strongly on system size, dynamics, and closure choice.

  • 4.3.1 Derivation: Nonlinear reaction systems generate an infinite hierarchy because moment equations of a given order depend on higher-order moments.Moment closures truncate this hierarchy by expressing moments above order M using lower-order moments.
  • 4.3.1 Derivation: Moment closure methods produce finite coupled ODE systems by assuming relationships among higher- and lower-order moments.Normal closure sets cumulants above order M to zero; 2MA and 3MA denote second- and third-order normal closures.
  • 4.3.1 Derivation: Common alternatives include Poisson, log-normal, central-moment-neglect, and derivative-matching closures, each imposing a different higher-moment approximation.Poisson closure uses mean-valued diagonal cumulants and zero mixed cumulants above the closure order, while derivative matching targets initial-time moment derivatives.
  • 4.3.1 Derivation: Moment equations become closed and integrable when the required higher-order moments are expressed in terms of lower-order moments.In the gene-system example, normal closure removes dependence on third- and higher-order moments.
  • 4.3.2 Properties and recent developments: Moment closures are computationally efficient because they solve finite ODE sets without ensemble averaging, but increasing closure order does not generally guarantee higher accuracy.For monostable systems, higher orders may become more accurate at sufficiently large system sizes; convergence is not generally expected at small sizes.
  • 4.3.2 Properties and recent developments: Moment closures can yield unphysical behavior, including negative means or variances, negative higher central moments, and diverging trajectories.Numerical studies report such behavior below critical volumes and non-physical oscillations or multistability in some larger-volume dynamical systems.

4.4 Construction of distributions from moments

Moment approximations can be converted into full distributions using maximum entropy, which matches selected moments while choosing the highest-entropy distribution. Combined with moment closure, this approach captures skewed steady-state protein distributions more accurately at higher closure order.

  • Maximum-entropy construction: Maximum entropy constructs a distribution whose first K moments match approximate moments obtained from methods such as the system size expansion.The construction is formulated as a nonlinear constrained optimization problem over distributions.
  • Maximum-entropy construction: The optimization uses Lagrange multipliers and reduces to an unconstrained problem with a normalization constant that can be solved numerically.The resulting optimization can be handled with standard numerical methods.
  • Numerical illustration: For bursty protein production, closure orders K = 3 and K = 5 combined with maximum entropy accurately capture the skewed marginal steady-state distribution.The exact reference is computed using the SSA.
  • Numerical illustration: Higher-order moment closure gives increasing accuracy for the reconstructed protein distribution.The comparison uses the central-moment-neglect closure at orders K = 3 and K = 5.
  • Related approach: Maximum entropy has also been used to close moment equations directly, but this is computationally expensive for time-dependent approximations and efficient and accurate for steady-state approximations.Time-dependent use requires iterative moment solving and multivariate optimization at small time steps.

4.5 Software

The review identifies freely available software for exact stochastic simulation and for several approximation methods, including system size expansion and moment closure.

  • Software: Dizzy, COPASI, StochKit, and StochPy are available packages for exact stochastic simulations.These include Java-based, stand-alone, and Python implementations.
  • Software: MOCA supports system size expansion beyond the linear-noise approximation and extends applicability to non-polynomial and time-dependent propensity functions.It does not require programming skills.
  • Software: MEANS extends moment closure methods to non-mass-action propensity functions, while CERENA implements exact simulations, system size expansion, and moment closure methods.Both packages are described as Python or Matlab tools, respectively.

4.6 Other approximations

Other approximations improve efficiency by exploiting favorable regimes, including approximate time stepping, time-scale separation, and hybrid discrete-continuous modeling. Their accuracy depends on assumptions, parameter regimes, species populations, and implementation choices.

  • Approximation methods: CLE, system size expansion, and moment closure methods are popular because they are relatively easy to implement, efficient, and often accurate, but can fail in some scenarios.The system size expansion additionally requires monostability.
  • Finite state projection: Finite state projection truncates an infinite state space and applies matrix exponentiation to approximate the distribution on the truncated subspace.Its computational feasibility depends on how large the truncated space must be for reasonable accuracy.
  • Tau-leaping: Tau-leaping advances the system in steps τ, sampling multiple reaction events together instead of simulating every event individually.The step must be small enough that propensity functions remain approximately constant.
  • Tau-leaping: Tau-leaping is more efficient when many reactions occur per step, but increasing τ lowers accuracy and low molecule numbers may require very small steps.For systems with very low particle numbers, tau-leaping can become less efficient than exact simulation.
  • Time-scale separation: Time-scale separation can reduce biochemical models when fast reactions or species can be eliminated while retaining their effects on slower dynamics.Deterministic reductions use quasi-equilibrium or quasi-steady-state assumptions, which are not generally equivalent.
  • Time-scale separation: In Michaelis-Menten kinetics, the quasi-equilibrium assumption yields the Michaelis-Menten product-production equation under the stated fast-binding and unbinding regime.The derivation uses enzyme conservation and the definition KM = k2/k1.
  • Time-scale separation: Stochastic quasi-steady-state reductions require stronger system conditions than deterministic reductions.The stochastic construction separates fast and slow reactions and defines species according to their participation in those reactions.
  • Hybrid methods: Hybrid methods must decide which species are discrete or continuous, and heuristic designs make performance assessment and convergence proofs difficult.Some hybrid methods nevertheless have derived error bounds or convergence results.

5 Comparison of approximation methods

The case study compares approximation accuracy across enzyme systems, saturation levels, molecule numbers, and burst sizes. Complex-valued CLE is generally most accurate, while LNA and real-valued CLE can fail under strong fluctuations, despite higher computational cost for CLE sampling.

  • Comparison setup: Few studies have directly compared the accuracy of CLE, system size expansion, and moment closure methods, motivating the numerical case study.The study examines a Michaelis–Menten-type enzyme system and an extension with transcription and translation.
  • Enzymatic protein degradation: For E0 = 2500, CLE-R, CLE-C, SSE-1, and 2MA approximate the mean well, whereas LNA shows significant deviations.For variance, all methods deviate more; CLE implementations perform best, followed by 2MA and SSE-1, with LNA least accurate.
  • Enzymatic protein degradation: Approximation errors increase for larger saturation parameter α, as the system approaches instability at α = 1 and substrate fluctuations become larger and more skewed.For α → 0, most enzymes are free, reducing the nonlinear effect of the bimolecular reaction S+E →C.
  • Enzymatic protein degradation: At E0 = 60, CLE-C reproduces the SSA mean and highly skewed protein distribution, while CLE-R and LNA give inaccurate means and distributions.The result indicates that rejecting steps to maintain positivity can be a poor treatment of the CLE boundary problem.
  • Discussion: CLE-C remains highly accurate for skewed distributions, whereas LNA cannot capture them because it predicts a Gaussian distribution.SSE-1 and 2MA perform similarly overall, with 2MA more accurate for the variance in Figure 9.
  • Discussion: Method choice depends on the target: CLE approximates full processes and distributions, while moment closure and system size expansion are more limited but can be substantially cheaper.System size expansion is systematic and accurate for large volumes, but is not applicable to deterministically multi-stable systems; moment closure is more flexible but ad hoc.

6 Inference

Inference for stochastic reaction networks targets joint uncertainty over trajectories and parameters, but exact forward propagation is usually unavailable because it requires solving the CME. The review therefore presents forward–backward methods and approximation-based alternatives, including approaches for open and structured networks.

  • Inference formulation: The inference problem arises because reaction parameters are often only approximately known, while distributions can change qualitatively with parameter values.The Bayesian setup treats the stochastic trajectory and parameters jointly.
  • Inference formulation: Bayesian inference computes the joint posterior p(x0:T, θ|y) over system trajectories and parameters, providing information about both state paths and parameter uncertainty.The posterior normalization is usually analytically intractable because it involves large sums or high-dimensional integrals.
  • Forward–backward inference: The Forward–Backward algorithm factorizes single-time posterior marginals into filtering distributions for past data and likelihoods of future data.The forward and backward factors can be computed recursively under the Markov assumption.
  • Forward–backward inference: Filtering alternates forward propagation between measurements with Bayesian measurement updates, yielding p(xt|yi≤t).Forward filtering combined with backward sampling can generate posterior trajectories for MCMC-based joint state and parameter inference.
  • Parameter inference: For parameter-only inference, either forward or backward recursion computes the marginal likelihood p(y), which can be optimized or combined with a prior to quantify parameter uncertainty.Optimizing p(y|θ) gives maximum-likelihood estimates; combining it with p(θ) yields a parameter posterior.
  • Methods for general networks: Forward–backward inference is difficult for chemical reaction networks because transition probabilities require solving the CME, generally impossible analytically.Numerical integration can work for closed systems with low molecule numbers, but open systems often require state-space truncation that introduces difficult-to-quantify bias.
  • Approximate and structured inference: Mesoscopic approximations such as CLE, LNA, and second-order normal moment closure reduce computational demands, but moment closure generally provides only moments rather than a full marginal distribution.Hybrid inference methods exploit prior structural knowledge, including promoter-state change points driving linear protein-dynamics SDEs.
  • Open challenges: The field lacks standard software and systematic comparisons of inference methods, limiting broader diffusion of the available approaches.

7 Conclusions

The review presents a self-contained overview of modelling, approximation, simulation, and Bayesian inference methods for stochastic chemical kinetics based on the CME. Its case study compares approximations with exact stochastic simulations, while the conclusions note that spatial systems violate common CME assumptions and remain harder to analyse.

  • The review introduces key approximation and inference methods for stochastic chemical kinetics and surveys recent developments.
  • In the case study, CLE-C was most accurate for moments and accurately approximated highly skewed steady-state distributions.CLE-R was less accurate for moments, while CLE-R and LNA did not accurately capture steady-state distributions.
  • SSE-1 and 2MA performed similarly and significantly better than the LNA for the evaluated moments.
  • The CME is valid for well-mixed, sufficiently dilute systems in which diffusion is the fastest time scale.
  • Many biological systems do not meet these assumptions, requiring spatial models whose analysis and approximation are substantially less developed and computationally expensive.SRDPs are mainly analysed algorithmically because their governing equations are difficult to solve or approximate.
  • The review aims to help scientists from other disciplines enter the field and stimulate further research in the presented areas.
Loading 1608.06582v2…