Multivariate moment closure techniques for stochastic kinetic models

Holimap: an accurate and efficient method for solving stochastic gene network dynamics

Gene-gene interactions are crucial to the control of sub-cellular processes but our understanding of their stochastic dynamics is hindered by the lack of simulation methods that can accurately and efficiently predict how the distributions of gene product numbers vary across parameter space. To overcome these difficulties, here we present Holimap (high-order linear-mapping approximation), an approach that approximates the protein or mRNA number distributions of a complex gene regulatory network by the distributions of a much simpler reaction system. We demonstrate Holimap's computational advantages over conventional methods by applying it to predict the stochastic time-dependent dynamics of various gene networks, including transcriptional networks ranging from simple autoregulatory loops to complex randomly connected networks, post-transcriptional networks, and post-translational networks. Holimap is ideally suited to study how the intricate network of gene-gene interactions results in precise coordination and control of gene expression.

A Provably Convergent Control Closure Scheme for the Method of Moments of the Chemical Master Equation

In this article, we introduce a novel moment closure scheme based on concepts from model predictive control (MPC) to accurately describe the time evolution of the statistical moments of the solution of the chemical master equation (CME). The method of moments, a set of ordinary differential equations frequently used to calculate the first nm moments, is generally not closed since lower-order moments depend on higher-order moments. To overcome this limitation, we interpret the moment equations as a nonlinear dynamical system, where the first nm moments serve as states, and the closing moments serve as the control input. We demonstrate the efficacy of our approach using three example systems and show that it outperforms existing closure schemes. For polynomial systems, which encompass all mass-action systems, we provide probability bounds for the error between true and estimated moment trajectories. We achieve this by combining the convergence properties of a priori moment estimates from stochastic simulations with guarantees for nonlinear reference tracking MPC. Our proposed method offers an effective solution to accurately predict the time evolution of moments of the CME, which has wide-ranging implications for many fields, including biology, chemistry, and engineering.

Bye bye, linearity, bye: quantification of the mean for linear CRNs in a random environment

Molecular reactions within a cell are inherently stochastic, and cells often differ in morphological properties or interact with a heterogeneous environment. Consequently, cell populations exhibit heterogeneity both due to these intrinsic and extrinsic causes. Although state-of-the-art studies that focus on dissecting this heterogeneity use single-cell measurements, the bulk data that shows only the mean expression levels is still in routine use. The fingerprint of the heterogeneity is present also in bulk data, despite being hidden from direct measurement. In particular, this heterogeneity can affect the mean expression levels via bimolecular interactions with low-abundant environment species. We make this statement rigorous for the class of linear reaction systems that are embedded in a discrete state Markov environment. The analytic expression that we provide for the stationary mean depends on the reaction rate constants of the linear subsystem, as well as the generator and stationary distribution of the Markov environment. We demonstrate the effect of the environment on the stationary mean. Namely, we show how the heterogeneous case deviates from the quasi-steady state (Q.SS) case when the embedded system is fast compared to the environment.

Extensions of mean-field approximations for environmentally-transmitted pathogen networks

Many pathogens spread via environmental transmission, without requiring host-to-host direct contact. While models for environmental transmission exist, many are simply constructed intuitively with structures analogous to standard models for direct transmission. As model insights are generally sensitive to the underlying model assumptions, it is important that we are able understand the details and consequences of these assumptions. We construct a simple network model for an environmentally-transmitted pathogen and rigorously derive systems of ordinary differential equations (ODEs) based on different assumptions. We explore two key assumptions, namely homogeneity and independence, and demonstrate that relaxing these assumptions can lead to more accurate ODE approximations. We compare these ODE models to a stochastic implementation of the network model over a variety of parameters and network structures, demonstrating that with fewer restrictive assumptions we are able to achieve higher accuracy in our approximations and highlighting more precisely the errors produced by each assumption. We show that less restrictive assumptions lead to more complicated systems of ODEs and the potential for unstable solutions. Due to the rigour of our derivation, we are able to identify the reason behind these errors and propose potential resolutions.

NLoed: A Python Package for Nonlinear Optimal Experimental Design in Systems Biology

Modeling in systems and synthetic biology relies on accurate parameter estimates and predictions. Accurate model calibration relies, in turn, on data and on how well suited the available data are to a particular modeling task. Optimal experimental design (OED) techniques can be used to identify experiments and data collection procedures that will most efficiently contribute to a given modeling objective. However, implementation of OED is limited by currently available software tools that are not well suited for the diversity of nonlinear models and non-normal data commonly encountered in biological research. Moreover, existing OED tools do not make use of the state-of-the-art numerical tools, resulting in inefficient computation. Here, we present the NLoed software package and demonstrate its use with in vivo data from an optogenetic system in Escherichia coli. NLoed is an open-source Python library providing convenient access to OED methods, with particular emphasis on experimental design for systems biology research. NLoed supports a wide variety of nonlinear, multi-input/output, and dynamic models and facilitates modeling and design of experiments over a wide variety of data types. To support OED investigations, the NLoed package implements maximum likelihood fitting and diagnostic tools, providing a comprehensive modeling workflow. NLoed offers an accessible, modular, and flexible OED tool set suited to the wide variety of experimental scenarios encountered in systems biology research. We demonstrate NLoed's capabilities by applying it to experimental design for characterization of a bacterial optogenetic system.

Quasi-Entropy Closure: a fast and reliable approach to close the moment equations of the Chemical Master Equation

MOTIVATION

The Chemical Master Equation is a stochastic approach to describe the evolution of a (bio)chemical reaction system. Its solution is a time-dependent probability distribution on all possible configurations of the system. As this number is typically large, the Master Equation is often practically unsolvable. The Method of Moments reduces the system to the evolution of a few moments, which are described by ordinary differential equations. Those equations are not closed, since lower order moments generally depend on higher order moments. Various closure schemes have been suggested to solve this problem. Two major problems with these approaches are first that they are open loop systems, which can diverge from the true solution, and second, some of them are computationally expensive.

RESULTS

Here we introduce Quasi-Entropy Closure, a moment-closure scheme for the Method of Moments. It estimates higher order moments by reconstructing the distribution that minimizes the distance to a uniform distribution subject to lower order moment constraints. Quasi-Entropy Closure can be regarded as an advancement of Zero-Information Closure, which similarly maximizes the information entropy. Results show that both approaches outperform truncation schemes. Quasi-Entropy Closure is computationally much faster than Zero-Information Closure, although both methods consider solutions on the space of configurations and hence do not completely overcome the curse of dimensionality. In addition, our scheme includes a plausibility check for the existence of a distribution satisfying a given set of moments on the feasible set of configurations. All results are evaluated on different benchmark problems.

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

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.

Revisiting moment-closure methods with heterogeneous multiscale population models

Stochastic chemical kinetics at the single-cell level give rise to heterogeneous populations of cells even when all individuals are genetically identical. This heterogeneity can lead to nonuniform behaviour within populations, including different growth characteristics, cell-fate dynamics, and response to stimuli. Ultimately, these diverse behaviours lead to intricate population dynamics that are inherently multiscale: the population composition evolves based on population-level processes that interact with stochastically distributed single-cell states. Therefore, descriptions that account for this heterogeneity are essential to accurately model and control chemical processes. However, for real-world systems such models are computationally expensive to simulate, which can make optimisation problems, such as optimal control or parameter inference, prohibitively challenging. Here, we consider a class of multiscale population models that incorporate population-level mechanisms while remaining faithful to the underlying stochasticity at the single-cell level and the interplay between these two scales. To address the complexity, we study an order-reduction approximations based on the distribution moments. Since previous moment-closure work has focused on the single-cell kinetics, extending these techniques to populations models prompts us to revisit old observations as well as tackle new challenges. In this extended multiscale context, we encounter the previously established observation that the simplest closure techniques can lead to non-physical system trajectories. Despite their poor performance in some systems, we provide an example where these simple closures outperform more sophisticated closure methods in accurately, efficiently, and robustly solving the problem of optimal control of bioproduction in a microbial consortium model.

Quantifying biochemical reaction rates from static population variability within incompletely observed complex networks

Quantifying biochemical reaction rates within complex cellular processes remains a key challenge of systems biology even as high-throughput single-cell data have become available to characterize snapshots of population variability. That is because complex systems with stochastic and non-linear interactions are difficult to analyze when not all components can be observed simultaneously and systems cannot be followed over time. Instead of using descriptive statistical models, we show that incompletely specified mechanistic models can be used to translate qualitative knowledge of interactions into reaction rate functions from covariability data between pairs of components. This promises to turn a globally intractable problem into a sequence of solvable inference problems to quantify complex interaction networks from incomplete snapshots of their stochastic fluctuations.

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.

Stochastic dynamics of predator-prey interactions

The interaction between a consumer (such as, a predator or a parasitoid) and a resource (such as, a prey or a host) forms an integral motif in ecological food webs, and has been modeled since the early 20th century starting from the seminal work of Lotka and Volterra. While the Lotka-Volterra predator-prey model predicts a neutrally stable equilibrium with oscillating population densities, a density-dependent predator attack rate is known to stabilize the equilibrium. Here, we consider a stochastic formulation of the Lotka-Volterra model where the prey's reproduction rate is a random process, and the predator's attack rate depends on both the prey and predator population densities. Analysis shows that increasing the sensitivity of the attack rate to the prey density attenuates the magnitude of stochastic fluctuations in the population densities. In contrast, these fluctuations vary non-monotonically with the sensitivity of the attack rate to the predator density with an optimal level of sensitivity minimizing the magnitude of fluctuations. Interestingly, our systematic study of the predator-prey correlations reveals distinct signatures depending on the form of the density-dependent attack rate. In summary, stochastic dynamics of nonlinear Lotka-Volterra models can be harnessed to infer density-dependent mechanisms regulating predator-prey interactions. Moreover, these mechanisms can have contrasting consequences on population density fluctuations, with predator-dependent attack rates amplifying stochasticity, while prey-dependent attack rates countering to buffer fluctuations.

Uncertainty propagation for deterministic models of biochemical networks using moment equations and the extended Kalman filter

Differential equation models of biochemical networks are frequently associated with a large degree of uncertainty in parameters and/or initial conditions. However, estimating the impact of this uncertainty on model predictions via Monte Carlo simulation is computationally demanding. A more efficient approach could be to track a system of low-order statistical moments of the state. Unfortunately, when the underlying model is nonlinear, the system of moment equations is infinite-dimensional and cannot be solved without a moment closure approximation which may introduce bias in the moment dynamics. Here, we present a new method to study the time evolution of the desired moments for nonlinear systems with polynomial rate laws. Our approach is based on solving a system of low-order moment equations by substituting the higher-order moments with Monte Carlo-based estimates from a small number of simulations, and using an extended Kalman filter to counteract Monte Carlo noise. Our algorithm provides more accurate and robust results compared to traditional Monte Carlo and moment closure techniques, and we expect that it will be widely useful for the quantification of uncertainty in biochemical model predictions.

Improved estimations of stochastic chemical kinetics by finite-state expansion

Stochastic reaction networks are a fundamental model to describe interactions between species where random fluctuations are relevant. The master equation provides the evolution of the probability distribution across the discrete state space consisting of vectors of population counts for each species. However, since its exact solution is often elusive, several analytical approximations have been proposed. The deterministic rate equation (DRE) gives a macroscopic approximation as a compact system of differential equations that estimate the average populations for each species, but it may be inaccurate in the case of nonlinear interaction dynamics. Here we propose finite-state expansion (FSE), an analytical method mediating between the microscopic and the macroscopic interpretations of a stochastic reaction network by coupling the master equation dynamics of a chosen subset of the discrete state space with the mean population dynamics of the DRE. An algorithm translates a network into an expanded one where each discrete state is represented as a further distinct species. This translation exactly preserves the stochastic dynamics, but the DRE of the expanded network can be interpreted as a correction to the original one. The effectiveness of FSE is demonstrated in models that challenge state-of-the-art techniques due to intrinsic noise, multi-scale populations and multi-stability.

MomentClosure.jl: automated moment closure approximations in Julia

SUMMARY

MomentClosure.jl is a Julia package providing automated derivation of the time-evolution equations of the moments of molecule numbers for virtually any chemical reaction network using a wide range of moment closure approximations. It extends the capabilities of modelling stochastic biochemical systems in Julia and can be particularly useful when exact analytic solutions of the chemical master equation are unavailable and when Monte Carlo simulations are computationally expensive.

AVAILABILITY AND IMPLEMENTATION

MomentClosure.jl is freely accessible under the MIT licence. Source code and documentation are available at https://github.com/augustinas1/MomentClosure.jl.

Compartor: A toolbox for the automatic generation of moment equations for dynamic compartment populations

Summary

Many biochemical processes in living organisms take place inside compartments that can interact with each other and remodel over time. In a recent work, we have shown how the stochastic dynamics of a compartmentalized biochemical system can be effectively studied using moment equations. With this technique, the time evolution of a compartment population is summarized using a finite number of ordinary differential equations, which can be analyzed very efficiently. However, the derivation of moment equations by hand can become time-consuming for systems comprising multiple reactants and interactions. Here we present Compartor, a toolbox that automatically generates the moment equations associated with a user-defined compartmentalized system. Through the moment equation method, Compartor renders the analysis of stochastic population models accessible to a broader scientific community.

Availability and implementation

Compartor is provided as a Python package and is available at https://pypi.org/project/compartor/. Source code and usage tutorials for Compartor are available at https://github.com/zechnerlab/Compartor.

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.

Moment dynamics for gene regulation with rational rate laws

This aim of this paper is mainly to investigate the performance of two typical moment closure schemes in gene regulatory master equations of rational rate laws. When the reaction rate is polynomial, the error bounds between the authentic and approximate moments obtained by schemes of Gaussian moment closure and log-normal moment closure are explicitly given. When the reaction rate is not polynomial, it is shown that the two schemes both behave well in the absence of active-inactive state switch, but in the presence of active-inactive state switch the log-normal closure scheme is far superior to the Gaussian closure scheme in capturing the asymptotic ensemble statistics. Moreover, the accuracy of the log-normal closure method is further confirmed by steady-state analytic results and the conditional Gaussian closure method. It is also disclosed that optimal negative feedback exists in suppressing protein noise in the presence of the on-off switch control.

Stochastic reaction networks in dynamic compartment populations

Many biochemical processes in living systems take place in compartmentalized environments, where individual compartments can interact with each other and undergo dynamic remodeling. Studying such processes through mathematical models poses formidable challenges because the underlying dynamics involve a large number of states, which evolve stochastically with time. Here we propose a mathematical framework to study stochastic biochemical networks in compartmentalized environments. We develop a generic population model, which tracks individual compartments and their molecular composition. We then show how the time evolution of this system can be studied effectively through a small number of differential equations, which track the statistics of the population. Our approach is versatile and renders an important class of biological systems computationally accessible.

Compartmentalization of biochemical processes underlies all biological systems, from the organelle to the tissue scale. Theoretical models to study the interplay between noisy reaction dynamics and compartmentalization are sparse, and typically very challenging to analyze computationally. Recent studies have made progress toward addressing this problem in the context of specific biological systems, but a general and sufficiently effective approach remains lacking. In this work, we propose a mathematical framework based on counting processes that allows us to study dynamic compartment populations with arbitrary interactions and internal biochemistry. We derive an efficient description of the dynamics in terms of differential equations which capture the statistics of the population. We demonstrate the relevance of our approach by analyzing models inspired by different biological processes, including subcellular compartmentalization and tissue homeostasis.

Stochastic master equation for early protein aggregation in the transthyretin amyloid disease

It is significant to understand the earliest molecular events occurring in the nucleation of the amyloid aggregation cascade for the prevention of amyloid related diseases such as transthyretin amyloid disease. We develop chemical master equation for the aggregation of monomers into oligomers using reaction rate law in chemical kinetics. For this stochastic model, lognormal moment closure method is applied to track the evolution of relevant statistical moments and its high accuracy is confirmed by the results obtained from Gillespie's stochastic simulation algorithm. Our results show that the formation of oligomers is highly dependent on the number of monomers. Furthermore, the misfolding rate also has an important impact on the process of oligomers formation. The quantitative investigation should be helpful for shedding more light on the mechanism of amyloid fibril nucleation.

Enhancement of gene expression noise from transcription factor binding to genomic decoy sites

The genome contains several high-affinity non-functional binding sites for transcription factors (TFs) creating a hidden and unexplored layer of gene regulation. We investigate the role of such “decoy sites” in controlling noise (random fluctuations) in the level of a TF that is synthesized in stochastic bursts. Prior studies have assumed that decoy-bound TFs are protected from degradation, and in this case decoys function to buffer noise. Relaxing this assumption to consider arbitrary degradation rates for both bound/unbound TF states, we find rich noise behaviors. For low-affinity decoys, noise in the level of unbound TF always monotonically decreases to the Poisson limit with increasing decoy numbers. In contrast, for high-affinity decoys, noise levels first increase with increasing decoy numbers, before decreasing back to the Poisson limit. Interestingly, while protection of bound TFs from degradation slows the time-scale of fluctuations in the unbound TF levels, the decay of bound TFs leads to faster fluctuations and smaller noise propagation to downstream target proteins. In summary, our analysis reveals stochastic dynamics emerging from nonspecific binding of TFs and highlights the dual role of decoys as attenuators or amplifiers of gene expression noise depending on their binding affinity and stability of the bound TF.

Extrinsic Noise and Heavy-Tailed Laws in Gene Expression

Noise in gene expression is one of the hallmarks of life at the molecular scale. Here we derive analytical solutions to a set of models describing the molecular mechanisms underlying transcription of DNA into RNA. Our ansatz allows us to incorporate the effects of extrinsic noise-encompassing factors external to the transcription of the individual gene-and discuss the ramifications for heterogeneity in gene product abundance that has been widely observed in single cell data. Crucially, we are able to show that heavy-tailed distributions of RNA copy numbers cannot result from the intrinsic stochasticity in gene expression alone, but must instead reflect extrinsic sources of variability.

Bounding the stationary distributions of the chemical master equation via mathematical programming

The stochastic dynamics of biochemical networks are usually modeled with the chemical master equation (CME). The stationary distributions of CMEs are seldom solvable analytically, and numerical methods typically produce estimates with uncontrolled errors. Here, we introduce mathematical programming approaches that yield approximations of these distributions with computable error bounds which enable the verification of their accuracy. First, we use semidefinite programming to compute increasingly tighter upper and lower bounds on the moments of the stationary distributions for networks with rational propensities. Second, we use these moment bounds to formulate linear programs that yield convergent upper and lower bounds on the stationary distributions themselves, their marginals, and stationary averages. The bounds obtained also provide a computational test for the uniqueness of the distribution. In the unique case, the bounds form an approximation of the stationary distribution with a computable bound on its error. In the nonunique case, our approach yields converging approximations of the ergodic distributions. We illustrate our methodology through several biochemical examples taken from the literature: Schlögl's model for a chemical bifurcation, a two-dimensional toggle switch, a model for bursty gene expression, and a dimerization model with multiple stationary distributions.

Comprehensive Review of Models and Methods for Inferences in Bio-Chemical Reaction Networks

The key processes in biological and chemical systems are described by networks of chemical reactions. From molecular biology to biotechnology applications, computational models of reaction networks are used extensively to elucidate their non-linear dynamics. The model dynamics are crucially dependent on the parameter values which are often estimated from observations. Over the past decade, the interest in parameter and state estimation in models of (bio-) chemical reaction networks (BRNs) grew considerably. The related inference problems are also encountered in many other tasks including model calibration, discrimination, identifiability, and checking, and optimum experiment design, sensitivity analysis, and bifurcation analysis. The aim of this review paper is to examine the developments in literature to understand what BRN models are commonly used, and for what inference tasks and inference methods. The initial collection of about 700 documents concerning estimation problems in BRNs excluding books and textbooks in computational biology and chemistry were screened to select over 270 research papers and 20 graduate research theses. The paper selection was facilitated by text mining scripts to automate the search for relevant keywords and terms. The outcomes are presented in tables revealing the levels of interest in different inference tasks and methods for given models in the literature as well as the research trends are uncovered. Our findings indicate that many combinations of models, tasks and methods are still relatively unexplored, and there are many new research opportunities to explore combinations that have not been considered-perhaps for good reasons. The most common models of BRNs in literature involve differential equations, Markov processes, mass action kinetics, and state space representations whereas the most common tasks are the parameter inference and model identification. The most common methods in literature are Bayesian analysis, Monte Carlo sampling strategies, and model fitting to data using evolutionary algorithms. The new research problems which cannot be directly deduced from the text mining data are also discussed.

Solving Stochastic Reaction Networks with Maximum Entropy Lagrange Multipliers

The time evolution of stochastic reaction networks can be modeled with the chemical master equation of the probability distribution. Alternatively, the numerical problem can be reformulated in terms of probability moment equations. Herein we present a new alternative method for numerically solving the time evolution of stochastic reaction networks. Based on the assumption that the entropy of the reaction network is maximum, Lagrange multipliers are introduced. The proposed method derives equations that model the time derivatives of these Lagrange multipliers. We present detailed steps to transform moment equations to Lagrange multiplier equations. In order to demonstrate the method, we present examples of non-linear stochastic reaction networks of varying degrees of complexity, including multistable and oscillatory systems. We find that the new approach is as accurate and significantly more efficient than Gillespie’s original exact algorithm for systems with small number of interacting species. This work is a step towards solving stochastic reaction networks accurately and efficiently.

Linear mapping approximation of gene regulatory networks with stochastic dynamics

The presence of protein-DNA binding reactions often leads to analytically intractable models of stochastic gene expression. Here we present the linear-mapping approximation that maps systems with protein-promoter interactions onto approximately equivalent systems with no binding reactions. This is achieved by the marriage of conditional mean-field approximation and the Magnus expansion, leading to analytic or semi-analytic expressions for the approximate time-dependent and steady-state protein number distributions. Stochastic simulations verify the method's accuracy in capturing the changes in the protein number distributions with time for a wide variety of networks displaying auto- and mutual-regulation of gene expression and independently of the ratios of the timescales governing the dynamics. The method is also used to study the first-passage time distribution of promoter switching, the sensitivity of the size of protein number fluctuations to parameter perturbation and the stochastic bifurcation diagram characterizing the onset of multimodality in protein number distributions.

Efficient simulation of intrinsic, extrinsic and external noise in biochemical systems

Motivation

Biological cells operate in a noisy regime influenced by intrinsic, extrinsic and external noise, which leads to large differences of individual cell states. Stochastic effects must be taken into account to characterize biochemical kinetics accurately. Since the exact solution of the chemical master equation, which governs the underlying stochastic process, cannot be derived for most biochemical systems, approximate methods are used to obtain a solution.

Results

In this study, a method to efficiently simulate the various sources of noise simultaneously is proposed and benchmarked on several examples. The method relies on the combination of the sigma point approach to describe extrinsic and external variability and the τ-leaping algorithm to account for the stochasticity due to probabilistic reactions. The comparison of our method to extensive Monte Carlo calculations demonstrates an immense computational advantage while losing an acceptable amount of accuracy. Additionally, the application to parameter optimization problems in stochastic biochemical reaction networks is shown, which is rarely applied due to its huge computational burden. To give further insight, a MATLAB script is provided including the proposed method applied to a simple toy example of gene expression.

Availability and implementation

MATLAB code is available at Bioinformatics online.

Supplementary information

Supplementary data are available at Bioinformatics online.

Optimization-based synthesis of stochastic biocircuits with statistical specifications

Model-guided design has become a standard approach to engineering biomolecular circuits in synthetic biology. However, the stochastic nature of biomolecular reactions is often overlooked in the design process. As a result, cell-cell heterogeneity causes unexpected deviation of biocircuit behaviours from model predictions and requires additional iterations of design-build-test cycles. To enhance the design process of stochastic biocircuits, this paper presents a computational framework to systematically specify the level of intrinsic noise using well-defined metrics of statistics and design highly heterogeneous biocircuits based on the specifications. Specifically, we use descriptive statistics of population distributions as an intuitive specification language of stochastic biocircuits and develop an optimization-based computational tool that explores parameter configurations satisfying design requirements. Sensitivity analysis methods are also performed to ensure the robustness of a biocircuit design against extrinsic perturbations. These design tools are formulated with convex optimization programs to enable rigorous and efficient quantification of the statistics. We demonstrate these features by designing a stochastic negative feedback biocircuit that satisfies multiple statistical constraints and perform an in-depth study of noise propagation and regulation in negative feedback pathways.

An efficient moments-based inference method for within-host bacterial infection dynamics

Over the last ten years, isogenic tagging (IT) has revolutionised the study of bacterial infection dynamics in laboratory animal models. However, quantitative analysis of IT data has been hindered by the piecemeal development of relevant statistical models. The most promising approach relies on stochastic Markovian models of bacterial population dynamics within and among organs. Here we present an efficient numerical method to fit such stochastic dynamic models to in vivo experimental IT data. A common approach to statistical inference with stochastic dynamic models relies on producing large numbers of simulations, but this remains a slow and inefficient method for all but simple problems, especially when tracking bacteria in multiple locations simultaneously. Instead, we derive and solve the systems of ordinary differential equations for the two lower-order moments of the stochastic variables (mean, variance and covariance). For any given model structure, and assuming linear dynamic rates, we demonstrate how the model parameters can be efficiently and accurately estimated by divergence minimisation. We then apply our method to an experimental dataset and compare the estimates and goodness-of-fit to those obtained by maximum likelihood estimation. While both sets of parameter estimates had overlapping confidence regions, the new method produced lower values for the division and death rates of bacteria: these improved the goodness-of-fit at the second time point at the expense of that of the first time point. This flexible framework can easily be applied to a range of experimental systems. Its computational efficiency paves the way for model comparison and optimal experimental design.

Control mechanisms for stochastic biochemical systems via computation of reachable sets

Controlling the behaviour of cells by rationally guiding molecular processes is an overarching aim of much of synthetic biology. Molecular processes, however, are notoriously noisy and frequently nonlinear. We present an approach to studying the impact of control measures on motifs of molecular interactions that addresses the problems faced in many biological systems: stochasticity, parameter uncertainty and nonlinearity. We show that our reachability analysis formalism can describe the potential behaviour of biological (naturally evolved as well as engineered) systems, and provides a set of bounds on their dynamics at the level of population statistics: for example, we can obtain the possible ranges of means and variances of mRNA and protein expression levels, even in the presence of uncertainty about model parameters.

Statistical physics approaches to subnetwork dynamics in biochemical systems

We apply a Gaussian variational approximation to model reduction in large biochemical networks of unary and binary reactions. We focus on a small subset of variables (subnetwork) of interest, e.g. because they are accessible experimentally, embedded in a larger network (bulk). The key goal is to write dynamical equations reduced to the subnetwork but still retaining the effects of the bulk. As a result, the subnetwork-reduced dynamics contains a memory term and an extrinsic noise term with non-trivial temporal correlations. We first derive expressions for this memory and noise in the linearized (Gaussian) dynamics and then use a perturbative power expansion to obtain first order nonlinear corrections. For the case of vanishing intrinsic noise, our description is explicitly shown to be equivalent to projection methods up to quadratic terms, but it is applicable also in the presence of stochastic fluctuations in the original dynamics. An example from the epidermal growth factor receptor signalling pathway is provided to probe the increased prediction accuracy and computational efficiency of our method.

Exact lower and upper bounds on stationary moments in stochastic biochemical systems

In the stochastic description of biochemical reaction systems, the time evolution of statistical moments for species population counts is described by a linear dynamical system. However, except for some ideal cases (such as zero- and first-order reaction kinetics), the moment dynamics is underdetermined as lower-order moments depend upon higher-order moments. Here, we propose a novel method to find exact lower and upper bounds on stationary moments for a given arbitrary system of biochemical reactions. The method exploits the fact that statistical moments of any positive-valued random variable must satisfy some constraints that are compactly represented through the positive semidefiniteness of moment matrices. Our analysis shows that solving moment equations at steady state in conjunction with constraints on moment matrices provides exact lower and upper bounds on the moments. These results are illustrated by three different examples-the commonly used logistic growth model, stochastic gene expression with auto-regulation and an activator-repressor gene network motif. Interestingly, in all cases the accuracy of the bounds is shown to improve as moment equations are expanded to include higher-order moments. Our results provide avenues for development of approximation methods that provide explicit bounds on moments for nonlinear stochastic systems that are otherwise analytically intractable.

Master equations and the theory of stochastic path integrals

This review provides a pedagogic and self-contained introduction to master equations and to their representation by path integrals. Since the 1930s, master equations have served as a fundamental tool to understand the role of fluctuations in complex biological, chemical, and physical systems. Despite their simple appearance, analyses of master equations most often rely on low-noise approximations such as the Kramers-Moyal or the system size expansion, or require ad-hoc closure schemes for the derivation of low-order moment equations. We focus on numerical and analytical methods going beyond the low-noise limit and provide a unified framework for the study of master equations. After deriving the forward and backward master equations from the Chapman-Kolmogorov equation, we show how the two master equations can be cast into either of four linear partial differential equations (PDEs). Three of these PDEs are discussed in detail. The first PDE governs the time evolution of a generalized probability generating function whose basis depends on the stochastic process under consideration. Spectral methods, WKB approximations, and a variational approach have been proposed for the analysis of the PDE. The second PDE is novel and is obeyed by a distribution that is marginalized over an initial state. It proves useful for the computation of mean extinction times. The third PDE describes the time evolution of a 'generating functional', which generalizes the so-called Poisson representation. Subsequently, the solutions of the PDEs are expressed in terms of two path integrals: a 'forward' and a 'backward' path integral. Combined with inverse transformations, one obtains two distinct path integral representations of the conditional probability distribution solving the master equations. We exemplify both path integrals in analysing elementary chemical reactions. Moreover, we show how a well-known path integral representation of averaged observables can be recovered from them. Upon expanding the forward and the backward path integrals around stationary paths, we then discuss and extend a recent method for the computation of rare event probabilities. Besides, we also derive path integral representations for processes with continuous state spaces whose forward and backward master equations admit Kramers-Moyal expansions. A truncation of the backward expansion at the level of a diffusion approximation recovers a classic path integral representation of the (backward) Fokker-Planck equation. One can rewrite this path integral in terms of an Onsager-Machlup function and, for purely diffusive Brownian motion, it simplifies to the path integral of Wiener. To make this review accessible to a broad community, we have used the language of probability theory rather than quantum (field) theory and do not assume any knowledge of the latter. The probabilistic structures underpinning various technical concepts, such as coherent states, the Doi-shift, and normal-ordered observables, are thereby made explicit.

MEANS: python package for Moment Expansion Approximation, iNference and Simulation

MOTIVATION

Many biochemical systems require stochastic descriptions. Unfortunately these can only be solved for the simplest cases and their direct simulation can become prohibitively expensive, precluding thorough analysis. As an alternative, moment closure approximation methods generate equations for the time-evolution of the system's moments and apply a closure ansatz to obtain a closed set of differential equations; that can become the basis for the deterministic analysis of the moments of the outputs of stochastic systems.

RESULTS

We present a free, user-friendly tool implementing an efficient moment expansion approximation with parametric closures that integrates well with the IPython interactive environment. Our package enables the analysis of complex stochastic systems without any constraints on the number of species and moments studied and the type of rate laws in the system. In addition to the approximation method our package provides numerous tools to help non-expert users in stochastic analysis.

AVAILABILITY AND IMPLEMENTATION

https://github.com/theosysbio/means

CONTACTS

m.stumpf@imperial.ac.uk or e.lakatos13@imperial.ac.uk

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

Reconstructing dynamic molecular states from single-cell time series

The notion of state for a system is prevalent in the quantitative sciences and refers to the minimal system summary sufficient to describe the time evolution of the system in a self-consistent manner. This is a prerequisite for a principled understanding of the inner workings of a system. Owing to the complexity of intracellular processes, experimental techniques that can retrieve a sufficient summary are beyond our reach. For the case of stochastic biomolecular reaction networks, we show how to convert the partial state information accessible by experimental techniques into a full system state using mathematical analysis together with a computational model. This is intimately related to the notion of conditional Markov processes and we introduce the posterior master equation and derive novel approximations to the corresponding infinite-dimensional posterior moment dynamics. We exemplify this state reconstruction approach using both in silico data and single-cell data from two gene expression systems in Saccharomyces cerevisiae, where we reconstruct the dynamic promoter and mRNA states from noisy protein abundance measurements.

Inferring extrinsic noise from single-cell gene expression data using approximate Bayesian computation

Background

Gene expression is known to be an intrinsically stochastic process which can involve single-digit numbers of mRNA molecules in a cell at any given time. The modelling of such processes calls for the use of exact stochastic simulation methods, most notably the Gillespie algorithm. However, this stochasticity, also termed “intrinsic noise”, does not account for all the variability between genetically identical cells growing in a homogeneous environment.

Despite substantial experimental efforts, determining appropriate model parameters continues to be a challenge. Methods based on approximate Bayesian computation can be used to obtain posterior parameter distributions given the observed data. However, such inference procedures require large numbers of simulations of the model and exact stochastic simulation is computationally costly.

In this work we focus on the specific case of trying to infer model parameters describing reaction rates and extrinsic noise on the basis of measurements of molecule numbers in individual cells at a given time point.

Results

To make the problem computationally tractable we develop an exact, model-specific, stochastic simulation algorithm for the commonly used two-state model of gene expression. This algorithm relies on certain assumptions and favourable properties of the model to forgo the simulation of the whole temporal trajectory of protein numbers in the system, instead returning only the number of protein and mRNA molecules present in the system at a specified time point. The computational gain is proportional to the number of protein molecules created in the system and becomes significant for systems involving hundreds or thousands of protein molecules.

Conclusions

We employ this simulation algorithm with approximate Bayesian computation to jointly infer the model’s rate and noise parameters from published gene expression data. Our analysis indicates that for most genes the extrinsic contributions to noise will be small to moderate but certainly are non-negligible.

Electronic supplementary material

The online version of this article (doi:10.1186/s12918-016-0324-x) contains supplementary material, which is available to authorized users.