The ability to sense the environment is heterogeneously distributed in cell populations

Channel capacity of signaling networks quantifies their fidelity in sensing extracellular inputs. Low estimates of channel capacities for several mammalian signaling networks suggest that cells can barely detect the presence/absence of environmental signals. However, given the extensive heterogeneity and temporal stability of cell state variables, we hypothesize that the sensing ability itself may depend on the state of the cells. In this work, we present an information-theoretic framework to quantify the distribution of sensing abilities from single-cell data. Using data on two mammalian pathways, we show that sensing abilities are widely distributed in the population and most cells achieve better resolution of inputs compared to an 'average cell'. We verify these predictions using live-cell imaging data on the IGFR/FoxO pathway. Importantly, we identify cell state variables that correlate with cells' sensing abilities. This information-theoretic framework will significantly improve our understanding of how cells sense in their environment.

Analytical refractory period distribution for a class of time-variant biochemical systems with second-order reactions

Refractory period (RP), the waiting time between signals, can induce complex signaling dynamics, such as acceleration, adaptation, and oscillation, within many cellular biochemical networks. However, its underlying molecular mechanisms are still unclear. Rigorously estimating the RP distribution may be essential to identify its causal regulatory mechanisms. Traditional methods of estimating the RP distribution depend on solving the underlying Chemical Master Equations (CMEs), the dominant modeling formalism of biochemical systems. However, exact solutions of the CME are only known for simple reaction systems with zero- and first-order reactions or specific systems with second-order reactions. General solutions still need to be derived for systems with bimolecular reactions. It is even more challenging if large state-space and nonconstant reaction rates are involved. Here, we developed a direct method to gain the analytical RP distribution for a class of second-order reaction systems with nonconstant reaction rates and large state space. Instead of using the CME, we used an equivalent path-wise representation, which is the solution to a transformed martingale problem of the CME. This allowed us to bypass solving a CME. We then applied the method to derive the analytical RP distribution of a real complex biochemical network with second-order reactions, the Drosophila phototransduction cascade. Our approach provides an alternative to the CMEs in deriving the analytical RP distributions of a class of second-order reaction systems. Since the bimolecular reactions are common in biological systems, our approach could enhance understanding real-world biochemical processes.

Effect of dedifferentiation on noise propagation in cellular hierarchy

Many fast renewing tissues have a hierarchical structure. Tissue-specific stem cells are at the root of this cellular hierarchy, which give rive to a whole range of specialized cells via cellular differentiation. However, increasing evidence shows that the hierarchical structure can be broken due to cellular dedifferentiation in which cells at differentiated stages can revert to the stem cell stage. Dedifferentiation has significant impacts on many aspects of hierarchical tissues. Here we investigate the effect of dedifferentiation on noise propagation by developing a stochastic model composed of different cell types. The moment equations are derived, via which we systematically investigate how the noise in the cell number is changed by dedifferentiation. Our results suggest that dedifferentiation have different effects on the noises in the numbers of stem cells and nonstem cells. Specifically, the noise in the number of stem cells is significantly reduced by increasing dedifferentiation probability. Due to the dual effect of dedifferentiation on nonstem cells, however, more complex changes could happen to the noise in the number of nonstem cells by increasing dedifferentiation probability. Furthermore, it is found that even though dedifferentiation could turn part of the noise propagation process into a noise-amplifying step, it is very unlikely to turn the entire process into a noise-amplifying cascade.

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.

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.

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.

Absolutely robust controllers for chemical reaction networks

In this work, we design a type of controller that consists of adding a specific set of reactions to an existing mass-action chemical reaction network in order to control a target species. This set of reactions is effective for both deterministic and stochastic networks, in the latter case controlling the mean as well as the variance of the target species. We employ a type of network property called absolute concentration robustness (ACR). We provide applications to the control of a multisite phosphorylation model as well as a receptor-ligand signalling system. For this framework, we use the so-called deficiency zero theorem from chemical reaction network theory as well as multiscaling model reduction methods. We show that the target species has approximately Poisson distribution with the desired mean. We further show that ACR controllers can bring robust perfect adaptation to a target species and are complementary to a recently introduced antithetic feedback controller used for stochastic chemical reactions.

Maximum Entropy Framework for Predictive Inference of Cell Population Heterogeneity and Responses in Signaling Networks

Predictive models of signaling networks are essential for understanding cell population heterogeneity and designing rational interventions in disease. However, using computational models to predict heterogeneity of signaling dynamics is often challenging because of the extensive variability of biochemical parameters across cell populations. Here, we describe a maximum entropy-based framework for inference of heterogeneity in dynamics of signaling networks (MERIDIAN). MERIDIAN estimates the joint probability distribution over signaling network parameters that is consistent with experimentally measured cell-to-cell variability of biochemical species. We apply the developed approach to investigate the response heterogeneity in the EGFR/Akt signaling network. Our analysis demonstrates that a significant fraction of cells exhibits high phosphorylated Akt (pAkt) levels hours after EGF stimulation. Our findings also suggest that cells with high EGFR levels predominantly contribute to the subpopulation of cells with high pAkt activity. We also discuss how MERIDIAN can be extended to accommodate various experimental measurements.

Discrete and continuous models of probability flux of switching dynamics: Uncovering stochastic oscillations in a toggle-switch system

The probability flux and velocity in stochastic reaction networks can help in characterizing dynamic changes in probability landscapes of these networks. Here, we study the behavior of three different models of probability flux, namely, the discrete flux model, the Fokker-Planck model, and a new continuum model of the Liouville flux. We compare these fluxes that are formulated based on, respectively, the chemical master equation, the stochastic differential equation, and the ordinary differential equation. We examine similarities and differences among these models at the nonequilibrium steady state for the toggle switch network under different binding and unbinding conditions. Our results show that at a strong stochastic condition of weak promoter binding, continuum models of Fokker-Planck and Liouville fluxes deviate significantly from the discrete flux model. Furthermore, we report the discovery of stochastic oscillation in the toggle-switch system occurring at weak binding conditions, a phenomenon captured only by the discrete flux model.

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.

OPTIMAL PEST REGULATION TACTICS FOR A STOCHASTIC PROCESS MODEL WITH IMPULSIVE CONTROLS USING REGRESSION ANALYSIS— TAKING COTTON APHIDS AS AN EXAMPLE

Aphids, the sap-sucking insects, often feeding in clusters on new plant growth, have resulted in large amounts of resources and efforts being spent attempting to control their activities. Taking cotton aphids as an example, this paper presents optimal control problems governed by stochastic models with impulsive interferences. Differing from the moment closure equation methods which are computationally intractable when the model contains excessive species, a new computational approach is employed to solve this problem. The key of the approach is to establish a functional relationship between the control variables involving the releasing rates of sterile insects and the spraying rates of pesticide and corresponding states on aphids and sterile aphids. Then the log-linear regression model is proposed to link the control variables with the moments (including the mean and variance) of states. Using training sample simulated from Gillespie algorithm, the regression coefficients for constraints and the objective function are estimated by least squares method. Simulation shows the error of the prediction of this model is relatively low and control results based on regression model are superior to the method based on the moment closure equations in terms of the control cost. Finally, the relative impacts of the prices and area of the field on optimal tactics are explored.

Chemical Kinetics Roots and Methods to Obtain the Probability Distribution Function Evolution of Reactants and Products in Chemical Networks Governed by a Master Equation

In this paper first, we review the physical root bases of chemical reaction networks as a Markov process in multidimensional vector space. Then we study the chemical reactions from a microscopic point of view, to obtain the expression for the propensities for the different reactions that can happen in the network. These chemical propensities, at a given time, depend on the system state at that time, and do not depend on the state at an earlier time indicating that we are dealing with Markov processes. Then the Chemical Master Equation (CME) is deduced for an arbitrary chemical network from a probability balance and it is expressed in terms of the reaction propensities. This CME governs the dynamics of the chemical system. Due to the difficulty to solve this equation two methods are studied, the first one is the probability generating function method or z-transform, which permits to obtain the evolution of the factorial moment of the system with time in an easiest way or after some manipulation the evolution of the polynomial moments. The second method studied is the expansion of the CME in terms of an order parameter (system volume). In this case we study first the expansion of the CME using the propensities obtained previously and splitting the molecular concentration into a deterministic part and a random part. An expression in terms of multinomial coefficients is obtained for the evolution of the probability of the random part. Then we study how to reconstruct the probability distribution from the moments using the maximum entropy principle. Finally, the previous methods are applied to simple chemical networks and the consistency of these methods is studied.

An algebraic method to calculate parameter regions for constrained steady-state distribution in stochastic reaction networks

Steady state is an essential concept in reaction networks. Its stability reflects fundamental characteristics of several biological phenomena such as cellular signal transduction and gene expression. Because biochemical reactions occur at the cellular level, they are affected by unavoidable fluctuations. Although several methods have been proposed to detect and analyze the stability of steady states for deterministic models, these methods cannot be applied to stochastic reaction networks. In this paper, we propose an algorithm based on algebraic computations to calculate parameter regions for constrained steady-state distribution of stochastic reaction networks, in which the means and variances satisfy some given inequality constraints. To evaluate our proposed method, we perform computer simulations for three typical chemical reactions and demonstrate that the results obtained with our method are consistent with the simulation results.

Detecting early-warning signals in periodically forced systems with noise

Early-warning signals for imminent regime shifts in multi-stable systems are highly desirable because it is often difficult to revert a system to the previous state once a transition has occurred. In this paper, two indicators, the phase lag and amplitude difference of the system's response, are extended to detect early-warning signals of a periodically driven, bistable complex system with noise. Our results show that both indicators can announce a regime shift of a complex system with small noise, namely, the critical point of the regime shift near a bifurcation point of the corresponding deterministic system. However, they fail to early indicate the regime shift in the case of large noise where the shift is far from the original bifurcation point. Based on the moment-expanding scheme, we reduce a large noise to a small one, and then both indicators work well again. We illustrate this approach via a parameterized lake eutrophication model verified by data. The regime shift to eutrophication could be detected in advance by studying the phase lag and amplitude difference of phosphorus concentrations. A basic statistical test is performed for the robustness of the proposed indicators. This approach provides a theoretical basis to prevent ecological environment deteriorations.

Estimation methods for heterogeneous cell population models in systems biology

Heterogeneity among individual cells is a characteristic and relevant feature of living systems. A range of experimental techniques to investigate this heterogeneity is available, and multiple modelling frameworks have been developed to describe and simulate the dynamics of heterogeneous populations. Measurement data are used to adjust computational models, which results in parameter and state estimation problems. Methods to solve these estimation problems need to take the specific properties of data and models into account. The aim of this review is to give an overview on the state of the art in estimation methods for heterogeneous cell population data and models. The focus is on models based on the population balance equation, but stochastic and individual-based models are also discussed. It starts with a brief discussion of common experimental approaches and types of measurement data that can be obtained in this context. The second part describes computational modelling frameworks for heterogeneous populations and the types of estimation problems occurring for these models. The third part starts with a discussion of observability and identifiability properties, after which the computational methods to solve the various estimation problems are described.

Quantifying the phosphorylation timescales of receptor-ligand complexes: a Markovian matrix-analytic approach

Cells interact with the extracellular environment by means of receptor molecules on their surface. Receptors can bind different ligands, leading to the formation of receptor-ligand complexes. For a subset of receptors, called receptor tyrosine kinases, binding to ligand enables sequential phosphorylation of intra-cellular residues, which initiates a signalling cascade that regulates cellular function and fate. Most mathematical modelling approaches employed to analyse receptor signalling are deterministic, especially when studying scenarios of high ligand concentration or large receptor numbers. There exist, however, biological scenarios where low copy numbers of ligands and/or receptors need to be considered, or where signalling by a few bound receptor-ligand complexes is enough to initiate a cellular response. Under these conditions stochastic approaches are appropriate, and in fact, different attempts have been made in the literature to measure the timescales of receptor signalling initiation in receptor-ligand systems. However, these approaches have made use of numerical simulations or approximations, such as moment-closure techniques. In this paper, we study, from an analytical perspective, the stochastic times to reach a given signalling threshold for two receptor-ligand models. We identify this time as an extinction time for a conveniently defined auxiliary absorbing continuous time Markov process, since receptor-ligand association/dissociation events can be analysed in terms of quasi-birth-and-death processes. We implement algorithmic techniques to compute the different order moments of this time, as well as the steady-state probability distribution of the system. A novel feature of the approach introduced here is that it allows one to quantify the role played by each kinetic rate in the timescales of signal initiation, and in the steady-state probability distribution of the system. Finally, we illustrate our approach by carrying out numerical studies for the vascular endothelial growth factor and one of its receptors, the vascular endothelial growth factor receptor of human endothelial cells.

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.

Identifiability and Reconstruction of Biochemical Reaction Networks from Population Snapshot Data

Inference of biochemical network models from experimental data is a crucial problem in systems and synthetic biology that includes parameter calibration but also identification of unknown interactions. Stochastic modelling from single-cell data is known to improve identifiability of reaction network parameters for specific systems. However, general results are lacking, and the advantage over deterministic, population-average approaches has not been explored for network reconstruction. In this work, we study identifiability and propose new reconstruction methods for biochemical interaction networks. Focusing on population-snapshot data and networks with reaction rates affine in the state, for parameter estimation, we derive general methods to test structural identifiability and demonstrate them in connection with practical identifiability for a reporter gene in silico case study. In the same framework, we next develop a two-step approach to the reconstruction of unknown networks of interactions. We apply it to compare the achievable network reconstruction performance in a deterministic and a stochastic setting, showing the advantage of the latter, and demonstrate it on population-snapshot data from a simulated example.

Abstracting the dynamics of biological pathways using information theory: a case study of apoptosis pathway

Motivation

Quantitative models are increasingly used in systems biology. Usually, these quantitative models involve many molecular species and their associated reactions. When simulating a tissue with thousands of cells, using these large models becomes computationally and time limiting.

Results

In this paper, we propose to construct abstractions using information theory notions. Entropy is used to discretize the state space and mutual information is used to select a subset of all original variables and their mutual dependencies. We apply our method to an hybrid model of TRAIL-induced apoptosis in HeLa cell. Our abstraction, represented as a Dynamic Bayesian Network (DBN), reduces the number of variables from 92 to 10, and accelerates numerical simulation by an order of magnitude, yet preserving essential features of cell death time distributions.

Availability and Implementation

This approach is implemented in the tool DBNizer, freely available at http://perso.crans.org/genest/DBNizer .

Contact

gregory.batt@inria.fr or bgenest@irisa.fr.

Supplementary information

Supplementary data are available at Bioinformatics online.

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.

Maximum Entropy Prediction of Non-Equilibrium Stationary Distributions for Stochastic Reaction Networks with Oscillatory Dynamics

Many chemical reaction networks in biological systems present complex oscillatory dynamics. In systems such as regulatory gene networks, cell cycle, and enzymatic processes, the number of molecules involved is often far from the thermodynamic limit. Although stochastic models based on the probabilistic approach of the Chemical Master Equation (CME) have been proposed, studies in the literature have been limited by the challenges of solving the CME and the lack of computational power to perform large-scale stochastic simulations. In this paper, we show that the infinite set of stationary moment equations describing the stochastic Brusselator and Schnakenberg oscillatory reactions networks can be truncated and solved using maximization of the entropy of the distributions. The results from our numerical experiments compare with the distributions obtained from well-established kinetic Monte Carlo methods and suggest that the accuracy of the prediction increases exponentially with the closure order chosen for the system. We conclude that maximum entropy models can be used as an efficient closure scheme alternative for moment equations to predict the non-equilibrium stationary distributions of stochastic chemical reactions with oscillatory dynamics. This prediction is accomplished without any prior knowledge of the system dynamics and without imposing any biased assumptions on the mathematical relations among species involved.

Diagnostics for assessing the linear noise and moment closure approximations

Solving the chemical master equation exactly is typically not possible, so instead we must rely on simulation based methods. Unfortunately, drawing exact realisations, results in simulating every reaction that occurs. This will preclude the use of exact simulators for models of any realistic size and so approximate algorithms become important. In this paper we describe a general framework for assessing the accuracy of the linear noise and two moment approximations. By constructing an efficient space filling design over the parameter region of interest, we present a number of useful diagnostic tools that aids modellers in assessing whether the approximation is suitable. In particular, we leverage the normality assumption of the linear noise and moment closure approximations.

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.

Splitting nodes and linking channels: A method for assembling biocircuits from stochastic elementary units

Akin to electric circuits, we construct biocircuits that are manipulated by cutting and assembling channels through which stochastic information flows. This diagrammatic manipulation allows us to create a method which constructs networks by joining building blocks selected so that (a) they cover only basic processes; (b) it is scalable to large networks; (c) the mean and variance-covariance from the Pauli master equation form a closed system; and (d) given the initial probability distribution, no special boundary conditions are necessary to solve the master equation. The method aims to help with both designing new synthetic signaling pathways and quantifying naturally existing regulatory networks.

Sparse Regression Based Structure Learning of Stochastic Reaction Networks from Single Cell Snapshot Time Series

Stochastic chemical reaction networks constitute a model class to quantitatively describe dynamics and cell-to-cell variability in biological systems. The topology of these networks typically is only partially characterized due to experimental limitations. Current approaches for refining network topology are based on the explicit enumeration of alternative topologies and are therefore restricted to small problem instances with almost complete knowledge. We propose the reactionet lasso, a computational procedure that derives a stepwise sparse regression approach on the basis of the Chemical Master Equation, enabling large-scale structure learning for reaction networks by implicitly accounting for billions of topology variants. We have assessed the structure learning capabilities of the reactionet lasso on synthetic data for the complete TRAIL induced apoptosis signaling cascade comprising 70 reactions. We find that the reactionet lasso is able to efficiently recover the structure of these reaction systems, ab initio, with high sensitivity and specificity. With only < 1% false discoveries, the reactionet lasso is able to recover 45% of all true reactions ab initio among > 6000 possible reactions and over 10²⁰⁰⁰ network topologies. In conjunction with information rich single cell technologies such as single cell RNA sequencing or mass cytometry, the reactionet lasso will enable large-scale structure learning, particularly in areas with partial network structure knowledge, such as cancer biology, and thereby enable the detection of pathological alterations of reaction networks. We provide software to allow for wide applicability of the reactionet lasso.

Virtually all biological processes are driven by biochemical reactions. However, their quantitative description in terms of stochastic chemical reaction networks is often precluded by the computational difficulty of structure learning, i.e. the identification of biologically active reaction networks among the combinatorially many possible topologies. This work describes the reactionet lasso, a structure learning approach that takes advantage of novel, information-rich single cell data and a tractable problem formulation to achieve structure learning for problem instances hundreds of orders of magnitude larger than previously reported. This approach opens the prospect of obtaining quantitative and predictive reaction models in many areas of biology and medicine, and in particular areas such as cancer biology, which are characterized by significant system alterations and many unknown reactions.

Effects of cell-cycle-dependent expression on random fluctuations in protein levels

Expression of many genes varies as a cell transitions through different cell-cycle stages. How coupling between stochastic expression and cell cycle impacts cell-to-cell variability (noise) in the level of protein is not well understood. We analyse a model where a stable protein is synthesized in random bursts, and the frequency with which bursts occur varies within the cell cycle. Formulae quantifying the extent of fluctuations in the protein copy number are derived and decomposed into components arising from the cell cycle and stochastic processes. The latter stochastic component represents contributions from bursty expression and errors incurred during partitioning of molecules between daughter cells. These formulae reveal an interesting trade-off: cell-cycle dependencies that amplify the noise contribution from bursty expression also attenuate the contribution from partitioning errors. We investigate the existence of optimum strategies for coupling expression to the cell cycle that minimize the stochastic component. Intriguingly, results show that a zero production rate throughout the cell cycle, with expression only occurring just before cell division, minimizes noise from bursty expression for a fixed mean protein level. By contrast, the optimal strategy in the case of partitioning errors is to make the protein just after cell division. We provide examples of regulatory proteins that are expressed only towards the end of the cell cycle, and argue that such strategies enhance robustness of cell-cycle decisions to the intrinsic stochasticity of gene expression.

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.

Experimental Design for Stochastic Models of Nonlinear Signaling Pathways Using an Interval-Wise Linear Noise Approximation and State Estimation

BACKGROUND

Computational modeling is a key technique for analyzing models in systems biology. There are well established methods for the estimation of the kinetic parameters in models of ordinary differential equations (ODE). Experimental design techniques aim at devising experiments that maximize the information encoded in the data. For ODE models there are well established approaches for experimental design and even software tools. However, data from single cell experiments on signaling pathways in systems biology often shows intrinsic stochastic effects prompting the development of specialized methods. While simulation methods have been developed for decades and parameter estimation has been targeted for the last years, only very few articles focus on experimental design for stochastic models.

METHODS

The Fisher information matrix is the central measure for experimental design as it evaluates the information an experiment provides for parameter estimation. This article suggest an approach to calculate a Fisher information matrix for models containing intrinsic stochasticity and high nonlinearity. The approach makes use of a recently suggested multiple shooting for stochastic systems (MSS) objective function. The Fisher information matrix is calculated by evaluating pseudo data with the MSS technique.

RESULTS

The performance of the approach is evaluated with simulation studies on an Immigration-Death, a Lotka-Volterra, and a Calcium oscillation model. The Calcium oscillation model is a particularly appropriate case study as it contains the challenges inherent to signaling pathways: high nonlinearity, intrinsic stochasticity, a qualitatively different behavior from an ODE solution, and partial observability. The computational speed of the MSS approach for the Fisher information matrix allows for an application in realistic size models.

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.

Intercellular Variability in Protein Levels from Stochastic Expression and Noisy Cell Cycle Processes

Inside individual cells, expression of genes is inherently stochastic and manifests as cell-to-cell variability or noise in protein copy numbers. Since proteins half-lives can be comparable to the cell-cycle length, randomness in cell-division times generates additional intercellular variability in protein levels. Moreover, as many mRNA/protein species are expressed at low-copy numbers, errors incurred in partitioning of molecules between two daughter cells are significant. We derive analytical formulas for the total noise in protein levels when the cell-cycle duration follows a general class of probability distributions. Using a novel hybrid approach the total noise is decomposed into components arising from i) stochastic expression; ii) partitioning errors at the time of cell division and iii) random cell-division events. These formulas reveal that random cell-division times not only generate additional extrinsic noise, but also critically affect the mean protein copy numbers and intrinsic noise components. Counter intuitively, in some parameter regimes, noise in protein levels can decrease as cell-division times become more stochastic. Computations are extended to consider genome duplication, where transcription rate is increased at a random point in the cell cycle. We systematically investigate how the timing of genome duplication influences different protein noise components. Intriguingly, results show that noise contribution from stochastic expression is minimized at an optimal genome-duplication time. Our theoretical results motivate new experimental methods for decomposing protein noise levels from synchronized and asynchronized single-cell expression data. Characterizing the contributions of individual noise mechanisms will lead to precise estimates of gene expression parameters and techniques for altering stochasticity to change phenotype of individual cells.

Inside individual cells, gene products often occur at low molecular counts and are subject to considerable stochastic fluctuations (noise) in copy numbers over time. An important consequence of noisy expression is that the level of a protein can vary considerably even among genetically identical cells exposed to the same environment. Such non-genetic phenotypic heterogeneity is physiologically relevant and critically influences diverse cellular processes. In addition to noise sources inherent in gene product synthesis, recent experimental studies have uncovered additional noise mechanisms that critically effect expression. For example, the time within the cell cycle when a gene duplicates, and the time taken to complete cell cycle are governed by random processes. The key contribution of this work is development of novel mathematical results quantifying how cell cycle-related noise sources combine with stochastic expression to drive intercellular variability in protein molecular counts. Derived formulas lead to many counterintuitive results, such as increasing randomness in the timing of cell division can lower noise in the level of a protein. Finally, these results inform experimental strategies to systematically dissect the contributions of different noise sources in the expression of a gene of interest.

Effects of distribution of infection rate on epidemic models

A goal of many epidemic models is to compute the outcome of the epidemics from the observed infected early dynamics. However, often, the total number of infected individuals at the end of the epidemics is much lower than predicted from the early dynamics. This discrepancy is argued to result from human intervention or nonlinear dynamics not incorporated in standard models. We show that when variability in infection rates is included in standard susciptible-infected-susceptible (SIS) and susceptible-infected-recovered (SIR) models the total number of infected individuals in the late dynamics can be orders lower than predicted from the early dynamics. This discrepancy holds for SIS and SIR models, where the assumption that all individuals have the same sensitivity is eliminated. In contrast with network models, fixed partnerships are not assumed. We derive a moment closure scheme capturing the distribution of sensitivities. We find that the shape of the sensitivity distribution does not affect R_{0} or the number of infected individuals in the early phases of the epidemics. However, a wide distribution of sensitivities reduces the total number of removed individuals in the SIR model and the steady-state infected fraction in the SIS model. The difference between the early and late dynamics implies that in order to extrapolate the expected effect of the epidemics from the initial phase of the epidemics, the rate of change in the average infectivity should be computed. These results are supported by a comparison of the theoretical model to the Ebola epidemics and by numerical simulation.

Inference for Stochastic Chemical Kinetics Using Moment Equations and System Size Expansion

Quantitative mechanistic models are valuable tools for disentangling biochemical pathways and for achieving a comprehensive understanding of biological systems. However, to be quantitative the parameters of these models have to be estimated from experimental data. In the presence of significant stochastic fluctuations this is a challenging task as stochastic simulations are usually too time-consuming and a macroscopic description using reaction rate equations (RREs) is no longer accurate. In this manuscript, we therefore consider moment-closure approximation (MA) and the system size expansion (SSE), which approximate the statistical moments of stochastic processes and tend to be more precise than macroscopic descriptions. We introduce gradient-based parameter optimization methods and uncertainty analysis methods for MA and SSE. Efficiency and reliability of the methods are assessed using simulation examples as well as by an application to data for Epo-induced JAK/STAT signaling. The application revealed that even if merely population-average data are available, MA and SSE improve parameter identifiability in comparison to RRE. Furthermore, the simulation examples revealed that the resulting estimates are more reliable for an intermediate volume regime. In this regime the estimation error is reduced and we propose methods to determine the regime boundaries. These results illustrate that inference using MA and SSE is feasible and possesses a high sensitivity.

In this manuscript, we introduce efficient methods for parameter estimation for stochastic processes. The stochasticity of chemical reactions can influence the average behavior of the considered system. For some biological systems, a microscopic, stochastic description is computationally intractable but a macroscopic, deterministic description too inaccurate. This inaccuracy manifests itself in an error in parameter estimates, which impede the predictive power of the proposed model. Until now, no rigorous analysis on the magnitude of the estimation error exists. We show by means of two simulation examples that using mesoscopic descriptions based on the system size expansions and moment-closure approximations can reduce this estimation error compared to inference using a macroscopic description. This reduction is most pronounced in an intermediate volume regime where the influence of stochasticity on the average behavior is moderately strong. For the JAK/STAT pathway where experimental data is available, we show that one parameter that was not structurally identifiable when using a macroscopic description becomes structurally identifiable when using a mesoscopic description for parameter estimation.

Data-driven optimizations for model checking of multi-valued regulatory networks

The model checking method has been long since established as an important tool for modelling and reverse engineering of biological systems. However, due to a high complexity of both the method and the biological systems, this approach often requires a vast amount of computational resources. In this article we show that by reducing the expressivity of the method we can gain performance while still being able to use all biologically relevant data. We utilize this approach to conduct a study of mutations in the EGFR signalling, motivated by a paper from Klinger et al. (2013). Here we aim at constructing approximated models of multiple cell-lines from sizeable sets of experimental data. Due to cancerous mutations in each cell line, there is a high degree of parameter uncertainty and the study would not be practically tractable without the performance optimizations described here.

CERENA: ChEmical REaction Network Analyzer--A Toolbox for the Simulation and Analysis of Stochastic Chemical Kinetics

Gene expression, signal transduction and many other cellular processes are subject to stochastic fluctuations. The analysis of these stochastic chemical kinetics is important for understanding cell-to-cell variability and its functional implications, but it is also challenging. A multitude of exact and approximate descriptions of stochastic chemical kinetics have been developed, however, tools to automatically generate the descriptions and compare their accuracy and computational efficiency are missing. In this manuscript we introduced CERENA, a toolbox for the analysis of stochastic chemical kinetics using Approximations of the Chemical Master Equation solution statistics. CERENA implements stochastic simulation algorithms and the finite state projection for microscopic descriptions of processes, the system size expansion and moment equations for meso- and macroscopic descriptions, as well as the novel conditional moment equations for a hybrid description. This unique collection of descriptions in a single toolbox facilitates the selection of appropriate modeling approaches. Unlike other software packages, the implementation of CERENA is completely general and allows, e.g., for time-dependent propensities and non-mass action kinetics. By providing SBML import, symbolic model generation and simulation using MEX-files, CERENA is user-friendly and computationally efficient. The availability of forward and adjoint sensitivity analyses allows for further studies such as parameter estimation and uncertainty analysis. The MATLAB code implementing CERENA is freely available from http://cerenadevelopers.github.io/CERENA/.

Identifying early-warning signals of critical transitions with strong noise by dynamical network markers

Identifying early-warning signals of a critical transition for a complex system is difficult, especially when the target system is constantly perturbed by big noise, which makes the traditional methods fail due to the strong fluctuations of the observed data. In this work, we show that the critical transition is not traditional state-transition but probability distribution-transition when the noise is not sufficiently small, which, however, is a ubiquitous case in real systems. We present a model-free computational method to detect the warning signals before such transitions. The key idea behind is a strategy: “making big noise smaller” by a distribution-embedding scheme, which transforms the data from the observed state-variables with big noise to their distribution-variables with small noise, and thus makes the traditional criteria effective because of the significantly reduced fluctuations. Specifically, increasing the dimension of the observed data by moment expansion that changes the system from state-dynamics to probability distribution-dynamics, we derive new data in a higher-dimensional space but with much smaller noise. Then, we develop a criterion based on the dynamical network marker (DNM) to signal the impending critical transition using the transformed higher-dimensional data. We also demonstrate the effectiveness of our method in biological, ecological and financial systems.

A termination criterion for parameter estimation in stochastic models in systems biology

Parameter estimation procedures are a central aspect of modeling approaches in systems biology. They are often computationally expensive, especially when the models take stochasticity into account. Typically parameter estimation involves the iterative optimization of an objective function that describes how well the model fits some measured data with a certain set of parameter values. In order to limit the computational expenses it is therefore important to apply an adequate stopping criterion for the optimization process, so that the optimization continues at least until a reasonable fit is obtained, but not much longer. In the case of stochastic modeling, at least some parameter estimation schemes involve an objective function that is itself a random variable. This means that plain convergence tests are not a priori suitable as stopping criteria. This article suggests a termination criterion suited to optimization problems in parameter estimation arising from stochastic models in systems biology. The termination criterion is developed for optimization algorithms that involve populations of parameter sets, such as particle swarm or evolutionary algorithms. It is based on comparing the variance of the objective function over the whole population of parameter sets with the variance of repeated evaluations of the objective function at the best parameter set. The performance is demonstrated for several different algorithms. To test the termination criterion we choose polynomial test functions as well as systems biology models such as an Immigration-Death model and a bistable genetic toggle switch. The genetic toggle switch is an especially challenging test case as it shows a stochastic switching between two steady states which is qualitatively different from the model behavior in a deterministic model.

Reconstructing the hidden states in time course data of stochastic models

Parameter estimation is central for analyzing models in Systems Biology. The relevance of stochastic modeling in the field is increasing. Therefore, the need for tailored parameter estimation techniques is increasing as well. Challenges for parameter estimation are partial observability, measurement noise, and the computational complexity arising from the dimension of the parameter space. This article extends the multiple shooting for stochastic systems' method, developed for inference in intrinsic stochastic systems. The treatment of extrinsic noise and the estimation of the unobserved states is improved, by taking into account the correlation between unobserved and observed species. This article demonstrates the power of the method on different scenarios of a Lotka-Volterra model, including cases in which the prey population dies out or explodes, and a Calcium oscillation system. Besides showing how the new extension improves the accuracy of the parameter estimates, this article analyzes the accuracy of the state estimates. In contrast to previous approaches, the new approach is well able to estimate states and parameters for all the scenarios. As it does not need stochastic simulations, it is of the same order of speed as conventional least squares parameter estimation methods with respect to computational time.

Chemical master equation closure for computer-aided synthetic biology

With inexpensive DNA synthesis technologies, we can now construct biological systems by quickly piecing together DNA sequences. Synthetic biology is the promising discipline that focuses on the construction of these new biological systems. Synthetic biology is an engineering discipline, and as such, it can benefit from mathematical modeling. This chapter focuses on mathematical models of biological systems. These models take the form of chemical reaction networks. The importance of stochasticity is discussed and methods to simulate stochastic reaction networks are reviewed. A closure scheme solution is also presented for the master equation of chemical reaction networks. The master equation is a complete model of randomly evolving molecular populations. Because of its ambitious character, the master equation remained unsolved for all but the simplest of molecular interaction networks for over 70 years. With the first complete solution of chemical master equations, a wide range of experimental observations of biomolecular interactions may be mathematically conceptualized. We anticipate that models based on the closure scheme described herein may assist in rationally designing synthetic biological systems.

Conditional Moment Closure Schemes for Studying Stochastic Dynamics of Genetic Circuits

Inside individual cells, stochastic expression drives random fluctuations in gene product copy numbers, which corrupts functioning of both natural and synthetic genetic circuits. Dynamic models of genetic circuits are formulated stochastically using the chemical master equation framework. Since obtaining probability distributions can be computationally expensive in these models, noise is typically investigated through lower-order statistical moments (mean, variance, correlation, skewness, etc.) of mRNA/proteins levels. However, due to the nonlinearities in genetic circuits, this moment dynamics is typically not closed, in the sense that the time derivative of the lower-order statistical moments depends on high-order moments. Moment equations are closed by expressing higher-order moments as nonlinear functions of lower-order moments, a technique commonly referred to as moment closure. We provide a new moment closure scheme for studying stochastic dynamics of genetic circuits, where genes randomly toggle between transcriptionally active and inactive states. The method is based on conditioning protein levels on active states of genes and then expressing higher-order moments as functions of lower-order conditional moments. The conditional closure scheme is illustrated on different circuit motifs and found to outperform existing closure techniques. Rapid computation of stochasticity through closure methods will enable improved characterization and design of synthetic circuits that exhibit robust performance in spite of noisy expression of underlying genes.

Delay chemical master equation: direct and closed-form solutions

The stochastic simulation algorithm (SSA) describes the time evolution of a discrete nonlinear Markov process. This stochastic process has a probability density function that is the solution of a differential equation, commonly known as the chemical master equation (CME) or forward-Kolmogorov equation. In the same way that the CME gives rise to the SSA, and trajectories of the latter are exact with respect to the former, trajectories obtained from a delay SSA are exact representations of the underlying delay CME (DCME). However, in contrast to the CME, no closed-form solutions have so far been derived for any kind of DCME. In this paper, we describe for the first time direct and closed solutions of the DCME for simple reaction schemes, such as a single-delayed unimolecular reaction as well as chemical reactions for transcription and translation with delayed mRNA maturation. We also discuss the conditions that have to be met such that such solutions can be derived.