The complex chemical Langevin equation

Noise facilitates entrainment of a population of uncoupled limit cycle oscillators

Many biological oscillators share two properties: they are subject to stochastic fluctuations (noise) and they must reliably adjust their period to changing environmental conditions (entrainment). While noise seems to distort the ability of single oscillators to entrain, in populations of uncoupled oscillators noise allows population-level entrainment for a wider range of input amplitudes and periods. Here, we investigate how this effect depends on the noise intensity and the number of oscillators in the population. We have found that, if a population consists of a sufficient number of oscillators, increasing noise intensity leads to faster entrainment after a phase change of the input signal (jet lag) and increases sensitivity to low-amplitude input signals.

The chemical Langevin equation for biochemical systems in dynamic environments

Modelling and simulation of complex biochemical reaction networks form cornerstones of modern biophysics. Many of the approaches developed so far capture temporal fluctuations due to the inherent stochasticity of the biophysical processes, referred to as intrinsic noise. Stochastic fluctuations, however, predominantly stem from the interplay of the network with many other - and mostly unknown - fluctuating processes, as well as with various random signals arising from the extracellular world; these sources contribute extrinsic noise. Here we provide a computational simulation method to probe the stochastic dynamics of biochemical systems subject to both intrinsic and extrinsic noise. We develop an extrinsic chemical Langevin equation-a physically motivated extension of the chemical Langevin equation- to model intrinsically noisy reaction networks embedded in a stochastically fluctuating environment. The extrinsic CLE is a continuous approximation to the Chemical Master Equation (CME) with time-varying propensities. In our approach, noise is incorporated at the level of the CME, and can account for the full dynamics of the exogenous noise process, irrespective of timescales and their mismatches. We show that our method accurately captures the first two moments of the stationary probability density when compared with exact stochastic simulation methods, while reducing the computational runtime by several orders of magnitude. Our approach provides a method that is practical, computationally efficient and physically accurate to study systems that are simultaneously subject to a variety of noise sources.

Steady-state joint distribution for first-order stochastic reaction kinetics

While the analytical solution for the marginal distribution of a stochastic chemical reaction network has been extensively studied, its joint distribution, i.e., the solution of a high-dimensional chemical master equation, has received much less attention. Here we develop an alternative method of computing the exact joint distributions of a wide class of first-order stochastic reaction systems in steady-state conditions. The effectiveness of our method is validated by applying it to four gene expression models of biological significance, including models with 2A peptides, nascent mRNA, gene regulation, translational bursting, and alternative splicing.

Three-dimensional stochastic simulation of chemoattractant-mediated excitability in cells

During the last decade, a consensus has emerged that the stochastic triggering of an excitable system drives pseudopod formation and subsequent migration of amoeboid cells. The presence of chemoattractant stimuli alters the threshold for triggering this activity and can bias the direction of migration. Though noise plays an important role in these behaviors, mathematical models have typically ignored its origin and merely introduced it as an external signal into a series of reaction-diffusion equations. Here we consider a more realistic description based on a reaction-diffusion master equation formalism to implement these networks. In this scheme, noise arises naturally from a stochastic description of the various reaction and diffusion terms. Working on a three-dimensional geometry in which separate compartments are divided into a tetrahedral mesh, we implement a modular description of the system, consisting of G-protein coupled receptor signaling (GPCR), a local excitation-global inhibition mechanism (LEGI), and signal transduction excitable network (STEN). Our models implement detailed biochemical descriptions whenever this information is available, such as in the GPCR and G-protein interactions. In contrast, where the biochemical entities are less certain, such as the LEGI mechanism, we consider various possible schemes and highlight the differences between them. Our simulations show that even when the LEGI mechanism displays perfect adaptation in terms of the mean level of proteins, the variance shows a dose-dependence. This differs between the various models considered, suggesting a possible means for determining experimentally among the various potential networks. Overall, our simulations recreate temporal and spatial patterns observed experimentally in both wild-type and perturbed cells, providing further evidence for the excitable system paradigm. Moreover, because of the overall importance and ubiquity of the modules we consider, including GPCR signaling and adaptation, our results will be of interest beyond the field of directed migration.

Though the term noise usually carries negative connotations, it can also contribute positively to the characteristic dynamics of a system. In biological systems, where noise arises from the stochastic interactions between molecules, its study is usually confined to genetic regulatory systems in which copy numbers are small and fluctuations large. However, noise can have important roles when the number of signaling molecules is large. The extension of pseudopods and the subsequent motion of amoeboid cells arises from the noise-induced trigger of an excitable system. Chemoattractant signals bias this triggering thereby directing cell motion. To date, this paradigm has not been tested by mathematical models that account accurately for the noise that arises in the corresponding reactions. In this study, we employ a reaction-diffusion master equation approach to investigate the effects of noise. Using a modular approach and a three-dimensional cell model with specific subdomains attributed to the cell membrane and cortex, we explore the spatiotemporal dynamics of the system. Our simulations recreate many experimentally-observed cell behaviors thereby supporting the biased-excitable network hypothesis.

Effects of noise and time delay on E2F's expression level in a bistable Rb-E2F gene's regulatory network

The bistable Rb-E2F gene regulatory network plays a central role in regulating cellular proliferation-quiescence transition. Based on Gillespie's chemical Langevin method, the stochastic bistable Rb-E2F gene's regulatory network with time delays is proposed. It is found that under the moderate intensity of internal noise, delay in the Cyclin E synthesis rate can greatly increase the average concentration value of E2F. When the delay is considered in both E2F-related positive feedback loops, within a specific range of delay (3-13) hr , the average expression of E2F is significantly increased. Also, this range is in the scope with that experimentally given by Dong et al. [65]. By analysing the quasi-potential curves at different delay times, simulation results show that delay regulates the dynamic behaviour of the system in the following way: small delay stabilises the bistable system; the medium delay is conducive to a high steady-state, making the system fluctuate near the high steady-state; large delay induces approximately periodic transitions between high and low steady-state. Therefore, by regulating noise and time delay, the cell itself can control the expression level of E2F to respond to different situations. These findings may provide an explanation of some experimental result intricacies related to the cell cycle.

Identifiability analysis for stochastic differential equation models in systems biology

Mathematical models are routinely calibrated to experimental data, with goals ranging from building predictive models to quantifying parameters that cannot be measured. Whether or not reliable parameter estimates are obtainable from the available data can easily be overlooked. Such issues of parameter identifiability have important ramifications for both the predictive power of a model, and the mechanistic insight that can be obtained. Identifiability analysis is well-established for deterministic, ordinary differential equation (ODE) models, but there are no commonly adopted methods for analysing identifiability in stochastic models. We provide an accessible introduction to identifiability analysis and demonstrate how existing ideas for analysis of ODE models can be applied to stochastic differential equation (SDE) models through four practical case studies. To assess structural identifiability, we study ODEs that describe the statistical moments of the stochastic process using open-source software tools. Using practically motivated synthetic data and Markov chain Monte Carlo methods, we assess parameter identifiability in the context of available data. Our analysis shows that SDE models can often extract more information about parameters than deterministic descriptions. All code used to perform the analysis is available on Github.

Stochastic time-dependent enzyme kinetics: Closed-form solution and transient bimodality

We derive an approximate closed-form solution to the chemical master equation describing the Michaelis-Menten reaction mechanism of enzyme action. In particular, assuming that the probability of a complex dissociating into an enzyme and substrate is significantly larger than the probability of a product formation event, we obtain expressions for the time-dependent marginal probability distributions of the number of substrate and enzyme molecules. For delta function initial conditions, we show that the substrate distribution is either unimodal at all times or else becomes bimodal at intermediate times. This transient bimodality, which has no deterministic counterpart, manifests when the initial number of substrate molecules is much larger than the total number of enzyme molecules and if the frequency of enzyme-substrate binding events is large enough. Furthermore, we show that our closed-form solution is different from the solution of the chemical master equation reduced by means of the widely used discrete stochastic Michaelis-Menten approximation, where the propensity for substrate decay has a hyperbolic dependence on the number of substrate molecules. The differences arise because the latter does not take into account enzyme number fluctuations, while our approach includes them. We confirm by means of a stochastic simulation of all the elementary reaction steps in the Michaelis-Menten mechanism that our closed-form solution is accurate over a larger region of parameter space than that obtained using the discrete stochastic Michaelis-Menten approximation.

Chemical Langevin equation: A path-integral view of Gillespie's derivation

In 2000, Gillespie rehabilitated the chemical Langevin equation (CLE) by describing two conditions that must be satisfied for it to yield a valid approximation of the chemical master equation (CME). In this work, we construct an original path-integral description of the CME and show how applying Gillespie's two conditions to it directly leads to a path-integral equivalent to the CLE. We compare this approach to the path-integral equivalent of a large system size derivation and show that they are qualitatively different. In particular, both approaches involve converting many sums into many integrals, and the difference between the two methods is essentially the difference between using the Euler-Maclaurin formula and using Riemann sums. Our results shed light on how path integrals can be used to conceptualize coarse-graining biochemical systems and are readily generalizable.

A practical guide to pseudo-marginal methods for computational inference in systems biology

For many stochastic models of interest in systems biology, such as those describing biochemical reaction networks, exact quantification of parameter uncertainty through statistical inference is intractable. Likelihood-free computational inference techniques enable parameter inference when the likelihood function for the model is intractable but the generation of many sample paths is feasible through stochastic simulation of the forward problem. The most common likelihood-free method in systems biology is approximate Bayesian computation that accepts parameters that result in low discrepancy between stochastic simulations and measured data. However, it can be difficult to assess how the accuracy of the resulting inferences are affected by the choice of acceptance threshold and discrepancy function. The pseudo-marginal approach is an alternative likelihood-free inference method that utilises a Monte Carlo estimate of the likelihood function. This approach has several advantages, particularly in the context of noisy, partially observed, time-course data typical in biochemical reaction network studies. Specifically, the pseudo-marginal approach facilitates exact inference and uncertainty quantification, and may be efficiently combined with particle filters for low variance, high-accuracy likelihood estimation. In this review, we provide a practical introduction to the pseudo-marginal approach using inference for biochemical reaction networks as a series of case studies. Implementations of key algorithms and examples are provided using the Julia programming language; a high performance, open source programming language for scientific computing (https://github.com/davidwarne/Warne2019_GuideToPseudoMarginal).

Chemical Langevin Equation for Complex Reactions

The stochastic dynamical behaviors of an elementary reaction system can be investigated by the chemical Langevin equation (CLE). However, most of the reactions in engineering belong to complex reactions. It is not appropriate to describe the random evolution process of a complex chemical reaction system by directly using CLE because the foundation of the deviation of CLE is the equilibrium equation of the number of molecules in elementary reaction systems. In the study, the chemical Langevin equation for complex reactions (CLE-CR) is proposed based on the random process theory by introducing the extent of reactions to express the reaction rates of complex reactions. The reaction rates of complex reactions are regarded as some random variables following Poisson distribution. To illustrate the essential consistency of CLE-CR and CLE, the physical meaning of the propensity function in CLE is comprehensively discussed. A numerical example from chemical engineering is employed to demonstrate the effectiveness of CLE-CR and the solving procedure. The results show that CLE-CR can be conveniently applied into engineering to investigate the stochastic dynamical behaviors of complex reaction systems, giving the probabilistic information of the concentration evolution of chemical constituents.

Neural field models for latent state inference: Application to large-scale neuronal recordings

Large-scale neural recording methods now allow us to observe large populations of identified single neurons simultaneously, opening a window into neural population dynamics in living organisms. However, distilling such large-scale recordings to build theories of emergent collective dynamics remains a fundamental statistical challenge. The neural field models of Wilson, Cowan, and colleagues remain the mainstay of mathematical population modeling owing to their interpretable, mechanistic parameters and amenability to mathematical analysis. Inspired by recent advances in biochemical modeling, we develop a method based on moment closure to interpret neural field models as latent state-space point-process models, making them amenable to statistical inference. With this approach we can infer the intrinsic states of neurons, such as active and refractory, solely from spiking activity in large populations. After validating this approach with synthetic data, we apply it to high-density recordings of spiking activity in the developing mouse retina. This confirms the essential role of a long lasting refractory state in shaping spatiotemporal properties of neonatal retinal waves. This conceptual and methodological advance opens up new theoretical connections between mathematical theory and point-process state-space models in neural data analysis.

Developing statistical tools to connect single-neuron activity to emergent collective dynamics is vital for building interpretable models of neural activity. Neural field models relate single-neuron activity to emergent collective dynamics in neural populations, but integrating them with data remains challenging. Recently, latent state-space models have emerged as a powerful tool for constructing phenomenological models of neural population activity. The advent of high-density multi-electrode array recordings now enables us to examine large-scale collective neural activity. We show that classical neural field approaches can yield latent state-space equations and demonstrate that this enables inference of the intrinsic states of neurons from recorded spike trains in large populations.

Spatial Stochastic Intracellular Kinetics: A Review of Modelling Approaches

Models of chemical kinetics that incorporate both stochasticity and diffusion are an increasingly common tool for studying biology. The variety of competing models is vast, but two stand out by virtue of their popularity: the reaction-diffusion master equation and Brownian dynamics. In this review, we critically address a number of open questions surrounding these models: How can they be justified physically? How do they relate to each other? How do they fit into the wider landscape of chemical models, ranging from the rate equations to molecular dynamics? This review assumes no prior knowledge of modelling chemical kinetics and should be accessible to a wide range of readers.

Insights into noise modulation in oligomerization systems of increasing complexity

Understanding under which conditions the increase of systems complexity is evolutionarily advantageous, and how this trend is related to the modulation of the intrinsic noise, are fascinating issues of utmost importance for synthetic and systems biology. To get insights into these matters, we analyzed a series of chemical reaction networks with different topologies and complexity, described by mass-action kinetics. We showed, analytically and numerically, that the global level of fluctuations at the steady state, measured by the sum over all species of the Fano factors of the number of molecules, is directly related to the network's deficiency. For zero-deficiency systems, this sum is constant and equal to the rank of the network. For higher deficiencies, additional terms appear in the Fano factor sum, which are proportional to the net reaction fluxes between the molecular complexes. We showed that the system's global intrinsic noise is reduced when all fluxes flow from lower to higher degree oligomers, or equivalently, towards the species of higher complexity, whereas it is amplified when the fluxes are directed towards lower complexity species.

Reduction of multiscale stochastic biochemical reaction networks using exact moment derivation

Biochemical reaction networks (BRNs) in a cell frequently consist of reactions with disparate timescales. The stochastic simulations of such multiscale BRNs are prohibitively slow due to high computational cost for the simulations of fast reactions. One way to resolve this problem uses the fact that fast species regulated by fast reactions quickly equilibrate to their stationary distribution while slow species are unlikely to be changed. Thus, on a slow timescale, fast species can be replaced by their quasi-steady state (QSS): their stationary conditional expectation values for given slow species. As the QSS are determined solely by the state of slow species, such replacement leads to a reduced model, where fast species are eliminated. However, it is challenging to derive the QSS in the presence of nonlinear reactions. While various approximation schemes for the QSS have been developed, they often lead to considerable errors. Here, we propose two classes of multiscale BRNs which can be reduced by deriving an exact QSS rather than approximations. Specifically, if fast species constitute either a feedforward network or a complex balanced network, the reduced model based on the exact QSS can be derived. Such BRNs are frequently observed in a cell as the feedforward network is one of fundamental motifs of gene or protein regulatory networks. Furthermore, complex balanced networks also include various types of fast reversible bindings such as bindings between transcriptional factors and gene regulatory sites. The reduced models based on exact QSS, which can be calculated by the computational packages provided in this work, accurately approximate the slow scale dynamics of the original full model with much lower computational cost.

Modelling biochemical reaction systems by stochastic differential equations with reflection

In this paper, we gave a new framework for modelling and simulating biochemical reaction systems by stochastic differential equations with reflection not in a heuristic way but in a mathematical way. The model is computationally efficient compared with the discrete-state Markov chain approach, and it ensures that both analytic and numerical solutions remain in a biologically plausible region. Specifically, our model mathematically ensures that species numbers lie in the domain D, which is a physical constraint for biochemical reactions, in contrast to the previous models. The domain D is actually obtained according to the structure of the corresponding chemical Langevin equations, i.e., the boundary is inherent in the biochemical reaction system. A variant of projection method was employed to solve the reflected stochastic differential equation model, and it includes three simple steps, i.e., Euler-Maruyama method was applied to the equations first, and then check whether or not the point lies within the domain D, and if not perform an orthogonal projection. It is found that the projection onto the closure D¯ is the solution to a convex quadratic programming problem. Thus, existing methods for the convex quadratic programming problem can be employed for the orthogonal projection map. Numerical tests on several important problems in biological systems confirmed the efficiency and accuracy of this approach.

Co-operative intermolecular kinetics of 2-oxoglutarate dependent dioxygenases may be essential for system-level regulation of plant cell physiology

Can the stimulus-driven synergistic association of 2-oxoglutarate dependent dioxygenases be influenced by the kinetic parameters of binding and catalysis?In this manuscript, I posit that these indices are necessary and specific for a particular stimulus, and are key determinants of a dynamic clustering that may function to mitigate the effects of this trigger. The protein(s)/sequence(s) that comprise this group are representative of all major kingdoms of life, and catalyze a generic hydroxylation, which is, in most cases accompanied by a specialized conversion of the substrate molecule. Iron is an essential co-factor for this transformation and the response to waning levels is systemic, and mandates the simultaneous participation of molecular sensors, transporters, and signal transducers. Here, I present a proof-of-concept model, that an evolving molecular network of 2OG-dependent enzymes can maintain iron homeostasis in the cytosol of root hair cells of members of the family Gramineae by actuating a non-reductive compensatory chelation by the phytosiderophores. Regression models of empirically available kinetic data (iron and alpha-ketoglutarate) were formulated, analyzed, and compared. The results, when viewed in context of the superfamily responding as a unit, suggest that members can indeed, work together to accomplish system-level function. This is achieved by the establishment of transient metabolic conduits, wherein the flux is dictated by kinetic compatibility of the participating enzymes. The approach adopted, i.e., predictive mathematical modeling, is integral to the hypothesis-driven acquisition of experimental data points and, in association with suitable visualization aids may be utilized for exploring complex plant biochemical systems.

Approximate probability distributions of the master equation

Master equations are common descriptions of mesoscopic systems. Analytical solutions to these equations can rarely be obtained. We here derive an analytical approximation of the time-dependent probability distribution of the master equation using orthogonal polynomials. The solution is given in two alternative formulations: a series with continuous and a series with discrete support, both of which can be systematically truncated. While both approximations satisfy the system size expansion of the master equation, the continuous distribution approximations become increasingly negative and tend to oscillations with increasing truncation order. In contrast, the discrete approximations rapidly converge to the underlying non-Gaussian distributions. The theory is shown to lead to particularly simple analytical expressions for the probability distributions of molecule numbers in metabolic reactions and gene expression systems.

General transient solution of the one-step master equation in one dimension

Exact analytical solutions of the master equation are limited to special cases and exact numerical methods are inefficient. Even the generic one-dimensional, one-step master equation has evaded exact solution, aside from the steady-state case. This type of master equation describes the dynamics of a continuous-time Markov process whose range consists of positive integers and whose transitions are allowed only between adjacent sites. The solution of any master equation can be written as the exponential of a (typically huge) matrix, which requires the calculation of the eigenvalues and eigenvectors of the matrix. Here we propose a linear algebraic method for simplifying this exponential for the general one-dimensional, one-step process. In particular, we prove that the calculation of the eigenvectors is actually not necessary for the computation of exponential, thereby we dramatically cut the time of this calculation. We apply our new methodology to examples from birth-death processes and biochemical networks. We show that the computational time is significantly reduced compared to existing methods.

Noise-induced multistability in chemical systems: Discrete versus continuum modeling

The noisy dynamics of chemical systems is commonly studied using either the chemical master equation (CME) or the chemical Fokker-Planck equation (CFPE). The latter is a continuum approximation of the discrete CME approach. It has recently been shown that for a particular system, the CFPE captures noise-induced multistability predicted by the CME. This phenomenon involves the CME's marginal probability distribution changing from unimodal to multimodal as the system size decreases below a critical value. We here show that the CFPE does not always capture noise-induced multistability. In particular we find simple chemical systems for which the CME predicts noise-induced multistability, whereas the CFPE predicts monostability for all system sizes.

The Spatial Chemical Langevin Equation and Reaction Diffusion Master Equations: moments and qualitative solutions

BACKGROUND

It has been established that stochastic effects play an important role in spatio-temporal biochemical networks. A popular method of representing such stochastic systems is the Reaction Diffusion Master Equation (RDME). However, simulating sample paths from the RDME can be computationally expensive, particularly at large populations. Here we investigate an uncommon, but much faster alternative: the Spatial Chemical Langevin Equation (SCLE).

METHODS

We investigate moment equations and correlation functions analytically, then we compare sample paths and moments of the SCLE to the RDME and associated deterministic solutions. Sample paths are generated computationally by the Next Subvolume method (RDME) and the Euler-Maruyama method (SCLE), while a deterministic solution is obtained with an Euler method. We consider the Gray-Scott model, a well-known pattern generating system, and a predator-prey system with spatially inhomogeneous parameters as sample applications.

RESULTS

For linear reaction networks, it is well known that the first order moments of all three approaches match, that the RDME and SCLE match to the second moment, and that all approaches diverge at third order moments. For non-linear reaction networks, differential equations governing moments do not form a closed system, but a general moment equation can be compared term wise. All approaches match at the leading order, and the RDME and SCLE match at the second leading order. As expected, the SCLE captures many dynamics of the RDME where deterministic methods fail to represent them. However, areas of the parameter space in the Gray-Scott model exist where either the SCLE and RDME give qualitatively different predictions, or the RDME predicts patterns, while the SCLE does not.

CONCLUSIONS

The SCLE provides a fast alternative to existing methods for simulation of spatial stochastic biochemical networks, capturing many aspects of dynamics represented by the RDME. This becomes very useful in search of quantitative parameters yielding desired qualitative solutions. However, there exist parameter sets where both the qualitative and quantitative behaviour of the SCLE can differ when compared to the RDME, so care should be taken in its use for applications demanding greater accuracy.