# Related Articles

Competence is a transiently differentiated state that certain bacterial cells reach when faced with a stressful environment. Entrance into competence can be attributed to the excitability of the dynamics governing the genetic circuit that regulates this cellular behavior. Like many biological behaviors, entrance into competence is a stochastic event. In this case cellular noise is responsible for driving the cell from a vegetative state into competence and back. In this work we present a novel numerical method for the analysis of stochastic biochemical events and use it to study the excitable dynamics responsible for competence in Bacillus subtilis. Starting with a Finite State Projection (FSP) solution of the chemical master equation (CME), we develop efficient numerical tools for accurately computing competence probability. Additionally, we propose a new approach for the sensitivity analysis of stochastic events and utilize it to elucidate the robustness properties of the competence regulatory genetic circuit. We also propose and implement a numerical method to calculate the expected time it takes a cell to return from competence. Although this study is focused on an example of cell-differentiation in Bacillus subtilis, our approach can be applied to a wide range of stochastic phenomena in biological systems.

Author Summary

When exposed to stress, organisms react by taking actions that help them protect their DNA. ComK protein is a key regulator which activates hundreds of genes, including the genes encoding the DNA-uptake and recombination systems. In Bacillus subtilis, stress in the environment activates a sequence of chemical reactions that, driven by cellular noise, stochastically increases the level of ComK in some bacterial cells driving them from their original vegetative state into a competent state. Entrance into and exit from competence are stochastic switching events that the cell undergoes. In this work, we present a novel numerical method that allows the analysis of stochastic events in biological systems. We illustrate our method by computing the probability with which Bacillus subtilis enters in competence. We also present a method to analyze the sensitivity of stochastic events. We use this method to study the sensitivity of the probability of entrance in competence with respect to various gene expressions and degradation rates. We finally present a numerical method to calculate the expected time it takes a cell to return from competence. Although we studied the competence regulatory genetic circuit, our approach can be applied to a variety of stochastic events in biological systems.

doi:10.1371/journal.pcbi.1000985

PMCID: PMC2978674
PMID: 21085679

Background

The chemical master equation (CME) is a system of ordinary differential equations that describes the evolution of a network of chemical reactions as a stochastic process. Its solution yields the probability density vector of the system at each point in time. Solving the CME numerically is in many cases computationally expensive or even infeasible as the number of reachable states can be very large or infinite. We introduce the sliding window method, which computes an approximate solution of the CME by performing a sequence of local analysis steps. In each step, only a manageable subset of states is considered, representing a "window" into the state space. In subsequent steps, the window follows the direction in which the probability mass moves, until the time period of interest has elapsed. We construct the window based on a deterministic approximation of the future behavior of the system by estimating upper and lower bounds on the populations of the chemical species.

Results

In order to show the effectiveness of our approach, we apply it to several examples previously described in the literature. The experimental results show that the proposed method speeds up the analysis considerably, compared to a global analysis, while still providing high accuracy.

Conclusions

The sliding window method is a novel approach to address the performance problems of numerical algorithms for the solution of the chemical master equation. The method efficiently approximates the probability distributions at the time points of interest for a variety of chemically reacting systems, including systems for which no upper bound on the population sizes of the chemical species is known a priori.

doi:10.1186/1752-0509-4-42

PMCID: PMC2867774
PMID: 20377904

Intrinsic noise is a common phenomenon in biochemical reaction networks and may affect the occurence and amplitude of sustained oscillations in the states of the network. To evaluate properties of such oscillations in the time domain, it is usually required to conduct long-term stochastic simulations, using for example the Gillespie algorithm. In this paper, we present a new method to compute the amplitude distribution of the oscillations without the need for long-term stochastic simulations. By the derivation of the method, we also gain insight into the structural features underlying the stochastic oscillations. The method is applicable to a wide class of non-linear stochastic differential equations that exhibit stochastic oscillations. The application is exemplified for the MAPK cascade, a fundamental element of several biochemical signalling pathways. This example shows that the proposed method can accurately predict the amplitude distribution for the stochastic oscillations even when using further computational approximations.

PACS Codes: 87.10.Mn, 87.18.Tt, 87.18.Vf

MSC Codes: 92B05, 60G10, 65C30

doi:10.1186/1757-5036-2-10

PMCID: PMC2796987
PMID: 19919689

Cellular processes are “noisy”. In each cell, concentrations of molecules are subject to random fluctuations due to the small numbers of these molecules and to environmental perturbations. While noise varies with time, it is often measured at steady state, for example by flow cytometry. When interrogating aspects of a cellular network by such steady-state measurements of network components, a key need is to develop efficient methods to simulate and compute these distributions. We describe innovations in stochastic modeling coupled with approaches to this computational challenge: first, an approach to modeling intrinsic noise via solution of the chemical master equation, and second, a convolution technique to account for contributions of extrinsic noise. We show how these techniques can be combined in a streamlined procedure for evaluation of different sources of variability in a biochemical network. Evaluation and illustrations are given in analysis of two well-characterized synthetic gene circuits, as well as a signaling network underlying the mammalian cell cycle entry.

Author Summary

Variability from one cell to another is a pronounced and universal trend in living organisms; much of this variability is related to varying concentrations of proteins and other chemical species across the cells. Understanding this variability is necessary if we are to fully understand cellular functions, particularly the ways in which cells differ from each other and in which cells with the same origin behave in different ways (e.g. in human development and cancer). When using a chemical model for some aspect of cellular function, one needs to consider two sources of variability: intrinsic variability, which results from the reactions proceeding as in the model but naturally varying because of the finite number of molecules in the cells and their random behavior; and extrinsic variability, which results from other kinds of variation not accounted for in the specific, deterministic model. We present new methods to model and compute both kinds of variability, to facilitate the study of cellular variability as a whole. Our methods provide advantages in speed, accuracy, and scope of mechanisms modeled, and we apply them to experimental data, demonstrating the nature of intrinsic and extrinsic noise in those systems.

doi:10.1371/journal.pcbi.1002209

PMCID: PMC3192818
PMID: 22022252

The aim of this study is to demonstrate that in molecular dynamical systems with the underlying bi- or multistability, the type of noise determines the most strongly attracting steady state or stochastic attractor. As an example we consider a simple stochastic model of autoregulatory gene with a nonlinear positive feedback, which in the deterministic approximation has two stable steady state solutions. Three types of noise are considered: transcriptional and translational – due to the small number of gene product molecules and the gene switching noise – due to gene activation and inactivation transitions. We demonstrate that the type of noise in addition to the noise magnitude dictates the allocation of probability mass between the two stable steady states. In particular, we found that when the gene switching noise dominates over the transcriptional and translational noise (which is characteristic of eukaryotes), the gene preferentially activates, while in the opposite case, when the transcriptional noise dominates (which is characteristic of prokaryotes) the gene preferentially remains inactive. Moreover, even in the zero-noise limit, when the probability mass generically concentrates in the vicinity of one of two steady states, the choice of the most strongly attracting steady state is noise type-dependent. Although the epigenetic attractors are defined with the aid of the deterministic approximation of the stochastic regulatory process, their relative attractivity is controlled by the type of noise, in addition to noise magnitude. Since noise characteristics vary during the cell cycle and development, such mode of regulation can be potentially employed by cells to switch between alternative epigenetic attractors.

doi:10.1016/j.jtbi.2012.10.004

PMCID: PMC3677133
PMID: 23063780

Gene expression; Bistability; Stochastic processes; Epigenetic attractors

Background

In a gene regulatory network (GRN), gene expressions are affected by noise, and stochastic fluctuations exist in the interactions among genes. These stochastic interactions are context dependent, thus it becomes important to consider noise in a context-sensitive manner in a network model. As a logical model, context-sensitive probabilistic Boolean networks (CSPBNs) account for molecular and genetic noise in the temporal context of gene functions. In a CSPBN with n genes and k contexts, however, a computational complexity of O(nk222n
) (or O(nk2
n
)) is required for an accurate (or approximate) computation of the state transition matrix (STM) of the size (2
n
∙ k) × (2
n
∙ k) (or 2
n
× 2
n
). The evaluation of a steady state distribution (SSD) is more challenging. Recently, stochastic Boolean networks (SBNs) have been proposed as an efficient implementation of an instantaneous PBN.

Results

The notion of stochastic Boolean networks (SBNs) is extended for the general model of PBNs, i.e., CSPBNs. This yields a novel structure of context-sensitive SBNs (CSSBNs) for modeling the stochasticity in a GRN. A CSSBN enables an efficient simulation of a CSPBN with a complexity of O(nLk2
n
) for computing the state transition matrix, where L is a factor related to the required sequence length in CSSBN for achieving a desired accuracy. A time-frame expanded CSSBN can further efficiently simulate the stationary behavior of a CSPBN and allow for a tunable tradeoff between accuracy and efficiency. The CSSBN approach is more efficient than an analytical method and more accurate than an approximate analysis.

Conclusions

Context-sensitive stochastic Boolean networks (CSSBNs) are proposed as an efficient approach to modeling the effects of gene perturbation and intervention in gene regulatory networks. A CSSBN analysis provides biologically meaningful insights into the oscillatory dynamics of the p53-Mdm2 network in a context-switching environment. It is shown that random gene perturbation has a greater effect on the final distribution of the steady state of a network compared to context switching activities. The CSSBN approach can further predict the steady state distribution of a glioma network under gene intervention. Ultimately, this will help drug discovery and develop effective drug intervention strategies.

doi:10.1186/1752-0509-8-60

PMCID: PMC4062525
PMID: 24886608

Gene regulatory networks; Boolean networks; Stochastic Boolean networks; Context dependent; Gene perturbation; Intervention; Context switch; Steady state distribution; p53 network; glioma network

As stochastic simulations become increasingly common in biological research, tools for analysis of such systems are in demand. The deterministic analogue to stochastic models, a set of probability moment equations equivalent to the Chemical Master Equation (CME), offers the possibility of a priori analysis of systems without the need for computationally costly Monte Carlo simulations. Despite the drawbacks of the method, in particular non-linearity in even the simplest of cases, the use of moment equations combined with moment-closure techniques has been used effectively in many fields. The techniques currently available to generate moment equations rely upon analytical expressions that are not efficient upon scaling. Additionally, the resulting moment-dependent matrix is lower diagonal and demands massive memory allocation in extreme cases. Here it is demonstrated that by utilizing factorial moments and the probability generating function (the Z-transform of the probability distribution) a recursive algorithm is produced. The resulting method is scalable and particularly efficient when high-order moments are required. The matrix produced is banded and often demands substantially less memory resources.

doi:10.1016/j.ces.2012.08.031

PMCID: PMC3501206
PMID: 23175571

Stochastic Simulation; Moments and Probability; Computational Chemistry; Mathematical Modelling; Kinetics; Numerical Analysis

The accepted stochastic descriptions of biochemical dynamics under well-mixed conditions are given by the Chemical Master Equation and the Stochastic Simulation Algorithm, which are equivalent. The latter is a Monte-Carlo method, which, despite enjoying broad availability in a large number of existing software packages, is computationally expensive due to the huge amounts of ensemble averaging required for obtaining accurate statistical information. The former is a set of coupled differential-difference equations for the probability of the system being in any one of the possible mesoscopic states; these equations are typically computationally intractable because of the inherently large state space. Here we introduce the software package intrinsic Noise Analyzer (iNA), which allows for systematic analysis of stochastic biochemical kinetics by means of van Kampen’s system size expansion of the Chemical Master Equation. iNA is platform independent and supports the popular SBML format natively. The present implementation is the first to adopt a complementary approach that combines state-of-the-art analysis tools using the computer algebra system Ginac with traditional methods of stochastic simulation. iNA integrates two approximation methods based on the system size expansion, the Linear Noise Approximation and effective mesoscopic rate equations, which to-date have not been available to non-expert users, into an easy-to-use graphical user interface. In particular, the present methods allow for quick approximate analysis of time-dependent mean concentrations, variances, covariances and correlations coefficients, which typically outperforms stochastic simulations. These analytical tools are complemented by automated multi-core stochastic simulations with direct statistical evaluation and visualization. We showcase iNA’s performance by using it to explore the stochastic properties of cooperative and non-cooperative enzyme kinetics and a gene network associated with circadian rhythms. The software iNA is freely available as executable binaries for Linux, MacOSX and Microsoft Windows, as well as the full source code under an open source license.

doi:10.1371/journal.pone.0038518

PMCID: PMC3373587
PMID: 22723865

Background

Stochastic biochemical reaction networks are commonly modelled by the chemical master equation, and can be simulated as first order linear differential equations through a finite state projection. Due to the very high state space dimension of these equations, numerical simulations are computationally expensive. This is a particular problem for analysis tasks requiring repeated simulations for different parameter values. Such tasks are computationally expensive to the point of infeasibility with the chemical master equation.

Results

In this article, we apply parametric model order reduction techniques in order to construct accurate low-dimensional parametric models of the chemical master equation. These surrogate models can be used in various parametric analysis task such as identifiability analysis, parameter estimation, or sensitivity analysis. As biological examples, we consider two models for gene regulation networks, a bistable switch and a network displaying stochastic oscillations.

Conclusions

The results show that the parametric model reduction yields efficient models of stochastic biochemical reaction networks, and that these models can be useful for systems biology applications involving parametric analysis problems such as parameter exploration, optimization, estimation or sensitivity analysis.

doi:10.1186/1752-0509-6-81

PMCID: PMC3532330
PMID: 22748204

Stochastic biochemical network; Model reduction; Reduced basis; Genetic regulatory network; Computational efficiency; Parameter estimation

Biological systems possess the ability to adapt quickly andadequately to both environmental and internal changes. This vital ability cannot be explained in terms ofconventional stochastic processes because such processes arecharacterized by atrade-off between flexibility and accuracy, that is, they either show shorttransition times (large Kramers escape rates) to broad steady-statedistributions or long transition times to sharply peaked distributions. To develop a stochastic theory for systemsexhibiting both flexibility and accuracy, we study systems under the impact of white noise multiplied with anaccordant statistical measure, here the probability density. Thisresults in negative feedback and circular causality: the more probable a stable state the lessit will be affected by noise and, conversely, the less a stable state is affected by noisethe more probable it is. Using nonlinear Fokker-Planckequations, steady states are computed via transformations ofsolutions of the corresponding linear Fokker-Planck equations. Transients reveal rapidly evolving and sharply peaked probability densities and thus mimic systems characterized by both flexibility and accuracy.

doi:10.1023/A:1016256613673

PMCID: PMC3456821
PMID: 23345756

Circular causality; flexibility-accuracy trade-offs; nonlinear Fokker-Planck equations; statistical feedback

In simulations of chemical systems, the main task is to find an exact or approximate solution of the chemical master equation (CME) that satisfies certain constraints with respect to computation time and accuracy. While Brownian motion simulations of single molecules are often too time consuming to represent the mesoscopic level, the classical Gillespie algorithm is a stochastically exact algorithm that provides satisfying results in the representation of calcium microdomains. Gillespie's algorithm can be approximated via
the tau-leap method and the chemical Langevin equation (CLE). Both methods lead to a substantial acceleration in computation time and a relatively small decrease in accuracy. Elimination of the noise terms leads to the classical, deterministic reaction rate equations (RRE). For complex multiscale systems, hybrid simulations are increasingly
proposed to combine the advantages of stochastic and deterministic algorithms. An often used exemplary cell type in this context are striated muscle cells (e.g., cardiac and skeletal muscle cells).
The properties of these cells are well described and they express many common calcium-dependent signaling
pathways. The purpose of the present paper is to provide an overview of the aforementioned simulation approaches and their mutual relationships in the spectrum ranging from stochastic to deterministic algorithms.

doi:10.1155/2011/572492

PMCID: PMC3216318
PMID: 22131814

We develop the stochastic, chemical master equation as a unifying approach to the dynamics of biochemical reaction systems in a mesoscopic volume under a living environment. A living environment provides a continuous chemical energy input that sustains the reaction system in a nonequilibrium steady state with concentration fluctuations. We discuss the linear, unimolecular single-molecule enzyme kinetics, phosphorylation-dephosphorylation cycle (PdPC) with bistability, and network exhibiting oscillations. Emphasis is paid to the comparison between the stochastic dynamics and the prediction based on the traditional approach based on the Law of Mass Action. We introduce the difference between nonlinear bistability and stochastic bistability, the latter has no deterministic counterpart. For systems with nonlinear bistability, there are three different time scales: (a) individual biochemical reactions, (b) nonlinear network dynamics approaching to attractors, and (c) cellular evolution. For mesoscopic systems with size of a living cell, dynamics in (a) and (c) are stochastic while that with (b) is dominantly deterministic. Both (b) and (c) are emergent properties of a dynamic biochemical network; We suggest that the (c) is most relevant to major cellular biochemical processes such as epi-genetic regulation, apoptosis, and cancer immunoediting. The cellular evolution proceeds with transitions among the attractors of (b) in a “punctuated equilibrium” manner.

doi:10.3390/ijms11093472

PMCID: PMC2956107
PMID: 20957107

chemical kinetics; chemical master equation; stochastic dynamics; biochemical reaction; systems biology

Background

Stochastic fluctuations in molecular numbers have been in many cases shown to be crucial for the understanding of biochemical systems. However, the systematic study of these fluctuations is severely hindered by the high computational demand of stochastic simulation algorithms. This is particularly problematic when, as is often the case, some or many model parameters are not well known. Here, we propose a solution to this problem, namely a combination of the linear noise approximation with optimisation methods. The linear noise approximation is used to efficiently estimate the covariances of particle numbers in the system. Combining it with optimisation methods in a closed-loop to find extrema of covariances within a possibly high-dimensional parameter space allows us to answer various questions. Examples are, what is the lowest amplitude of stochastic fluctuations possible within given parameter ranges? Or, which specific changes of parameter values lead to the increase of the correlation between certain chemical species? Unlike stochastic simulation methods, this has no requirement for small numbers of molecules and thus can be applied to cases where stochastic simulation is prohibitive.

Results

We implemented our strategy in the software COPASI and show its applicability on two different models of mitogen-activated kinases (MAPK) signalling -- one generic model of extracellular signal-regulated kinases (ERK) and one model of signalling via p38 MAPK. Using our method we were able to quickly find local maxima of covariances between particle numbers in the ERK model depending on the activities of phospho-MKKK and its corresponding phosphatase. With the p38 MAPK model our method was able to efficiently find conditions under which the coefficient of variation of the output of the signalling system, namely the particle number of Hsp27, could be minimised. We also investigated correlations between the two parallel signalling branches (MKK3 and MKK6) in this model.

Conclusions

Our strategy is a practical method for the efficient investigation of fluctuations in biochemical models even when some or many of the model parameters have not yet been fully characterised.

doi:10.1186/1752-0509-6-86

PMCID: PMC3814289
PMID: 22805626

Linear noise approximation; Optimisation; Mitogen-activated kinases signalling; COPASI; Intrinsic noise; Stochastic biochemical models; Systems biology

The inference of reaction rate parameters in biochemical network models from time series concentration data is a central task in computational systems biology. Under the assumption of well mixed conditions the network dynamics are typically described by the chemical master equation, the Fokker Planck equation, the linear noise approximation or the macroscopic rate equation. The inverse problem of estimating the parameters of the underlying network model can be approached in deterministic and stochastic ways, and available methods often compare individual or mean concentration traces obtained from experiments with theoretical model predictions when maximizing likelihoods, minimizing regularized least squares functionals, approximating posterior distributions or sequentially processing the data. In this article we assume that the biological reaction network can be observed at least partially and repeatedly over time such that sample moments of species molecule numbers for various time points can be calculated from the data. Based on the chemical master equation we furthermore derive closed systems of parameter dependent nonlinear ordinary differential equations that predict the time evolution of the statistical moments. For inferring the reaction rate parameters we suggest to not only compare the sample mean with the theoretical mean prediction but also to take the residual of higher order moments explicitly into account. Cost functions that involve residuals of higher order moments may form landscapes in the parameter space that have more pronounced curvatures at the minimizer and hence may weaken or even overcome parameter sloppiness and uncertainty. As a consequence both deterministic and stochastic parameter inference algorithms may be improved with respect to accuracy and efficiency. We demonstrate the potential of moment fitting for parameter inference by means of illustrative stochastic biological models from the literature and address topics for future research.

doi:10.1371/journal.pone.0043001

PMCID: PMC3416831
PMID: 22900079

Control in the natural environment is difficult in part because of uncertainty in the effect of actions. Uncertainty can be due to added motor or sensory noise, unmodeled dynamics, or quantization of sensory feedback. Biological systems are faced with further difficulties, since control must be performed by networks of cooperating neurons and neural subsystems. Here, we propose a new mathematical framework for modeling and simulation of distributed control systems operating in an uncertain environment. Stochastic Differential Operators can be derived from the stochastic differential equation describing a system, and they map the current state density into the differential of the state density. Unlike discrete-time Markov update operators, stochastic differential operators combine linearly for a large class of linear and nonlinear systems, and therefore the combined effects of multiple controllable and uncontrollable subsystems can be predicted. Design using these operators yields systems whose statistical behavior can be specified throughout state space. The relationship to Bayesian estimation and discrete-time Markov processes is described.

doi:10.1162/NECO_a_00151

PMCID: PMC3086813
PMID: 21521040

Stochastic control; Bayes rule; Markov operators

Background

Fluorescent and luminescent gene reporters allow us to dynamically quantify changes in molecular species concentration over time on the single cell level. The mathematical modeling of their interaction through multivariate dynamical models requires the deveopment of effective statistical methods to calibrate such models against available data. Given the prevalence of stochasticity and noise in biochemical systems inference for stochastic models is of special interest. In this paper we present a simple and computationally efficient algorithm for the estimation of biochemical kinetic parameters from gene reporter data.

Results

We use the linear noise approximation to model biochemical reactions through a stochastic dynamic model which essentially approximates a diffusion model by an ordinary differential equation model with an appropriately defined noise process. An explicit formula for the likelihood function can be derived allowing for computationally efficient parameter estimation. The proposed algorithm is embedded in a Bayesian framework and inference is performed using Markov chain Monte Carlo.

Conclusion

The major advantage of the method is that in contrast to the more established diffusion approximation based methods the computationally costly methods of data augmentation are not necessary. Our approach also allows for unobserved variables and measurement error. The application of the method to both simulated and experimental data shows that the proposed methodology provides a useful alternative to diffusion approximation based methods.

doi:10.1186/1471-2105-10-343

PMCID: PMC2774326
PMID: 19840370

Phenotypic differences of genetically identical cells under the same environmental conditions have been attributed to the inherent stochasticity of biochemical processes. Various mechanisms have been suggested, including the existence of alternative steady states in regulatory networks that are reached by means of stochastic fluctuations, long transient excursions from a stable state to an unstable excited state, and the switching on and off of a reaction network according to the availability of a constituent chemical species. Here we analyse a detailed stochastic kinetic model of two-component system signalling in bacteria, and show that alternative phenotypes emerge in the absence of these features. We perform a bifurcation analysis of deterministic reaction rate equations derived from the model, and find that they cannot reproduce the whole range of qualitative responses to external signals demonstrated by direct stochastic simulations. In particular, the mixed mode, where stochastic switching and a graded response are seen simultaneously, is absent. However, probabilistic and equation-free analyses of the stochastic model that calculate stationary states for the mean of an ensemble of stochastic trajectories reveal that slow transcription of either response regulator or histidine kinase leads to the coexistence of an approximate basal solution and a graded response that combine to produce the mixed mode, thus establishing its essential stochastic nature. The same techniques also show that stochasticity results in the observation of an all-or-none bistable response over a much wider range of external signals than would be expected on deterministic grounds. Thus we demonstrate the application of numerical equation-free methods to a detailed biochemical reaction network model, and show that it can provide new insight into the role of stochasticity in the emergence of phenotypic diversity.

Author Summary

It is a surprising fact that genetically identical bacteria, living in identical conditions, can develop in completely different ways: for example, one subpopulation might grow very fast and another very slowly. These different phenotypes are thought to be one reason why bacteria that cause disease can survive antibiotic treatment or become persistent. This diversity of behaviour is usually attributed to the existence of multiple stable phenotypic states, or to the coexistence of one stable state with another unstable excited state, or finally to the possibility of the whole biochemical system that controls the phenotype being switched on and off. In this paper we describe a different scenario that leads to phenotypic diversity in two-component system signalling, a very common mechanism that bacteria use to sense external signals and control their response to changes in their environment. We use probability theory and equation-free computational analysis to calculate the average number of molecules of each chemical species present in the two-component system and hence show that sporadic production of either of two key chemical components required for signalling can delay the response to the external signal in some bacterial cells and so lead to the emergence of two distinct cell populations.

doi:10.1371/journal.pcbi.1002396

PMCID: PMC3386199
PMID: 22761552

The Chemical Master Equation (CME) is a cornerstone of stochastic analysis and simulation of models of biochemical reaction networks. Yet direct solutions of the CME have remained elusive. Although several approaches overcome the infinite dimensional nature of the CME through projections or other means, a common feature of proposed approaches is their susceptibility to the curse of dimensionality, i.e. the exponential growth in memory and computational requirements in the number of problem dimensions. We present a novel approach that has the potential to “lift” this curse of dimensionality. The approach is based on the use of the recently proposed Quantized Tensor Train (QTT) formatted numerical linear algebra for the low parametric, numerical representation of tensors. The QTT decomposition admits both, algorithms for basic tensor arithmetics with complexity scaling linearly in the dimension (number of species) and sub-linearly in the mode size (maximum copy number), and a numerical tensor rounding procedure which is stable and quasi-optimal. We show how the CME can be represented in QTT format, then use the exponentially-converging -discontinuous Galerkin discretization in time to reduce the CME evolution problem to a set of QTT-structured linear equations to be solved at each time step using an algorithm based on Density Matrix Renormalization Group (DMRG) methods from quantum chemistry. Our method automatically adapts the “basis” of the solution at every time step guaranteeing that it is large enough to capture the dynamics of interest but no larger than necessary, as this would increase the computational complexity. Our approach is demonstrated by applying it to three different examples from systems biology: independent birth-death process, an example of enzymatic futile cycle, and a stochastic switch model. The numerical results on these examples demonstrate that the proposed QTT method achieves dramatic speedups and several orders of magnitude storage savings over direct approaches.

Author Summary

Stochastic models of chemical networks are necessary to quantitatively describe random fluctuations and other probabilistic phenomena within living cells. The Chemical Master Equation (CME) describes the time evolution of molecular abundance probabilities in these models, and is a basis for many stochastic simulation and analysis methods. Yet the CME is difficult to solve directly except for very simple structures. Indeed current approaches are susceptible to the curse of dimensionality, that is, the exponential growth of memory and computational requirements in the number of problem dimensions. In this paper, we propose a novel approach that has the potential to overcome the curse of dimensionality. It is based on the use of the recently proposed Quantized Tensor Train (QTT) formatted numerical linear algebra for numerical representation of tensors, using algorithms for basic tensor arithmetics with complexity scaling linearly in the number of reacting species considered, and sub-linearly in the maximum allowed copy number per species. We present this approach and demonstrate its effectiveness by applying it to three problems from systems biology. Numerical experiments are reported which show that several orders of magnitude memory savings is typically afforded by the new approach presented here.

doi:10.1371/journal.pcbi.1003359

PMCID: PMC3953644
PMID: 24626049

True steady states are a rare occurrence in living organisms, yet their knowledge is essential for quasi-steady state approximations, multistability analysis, and other important tools in the investigation of chemical reaction networks (CRN) used to describe molecular processes on the cellular level. Here we present an approach that can provide closed form steady-state solutions to complex systems, resulting from CRN with binary reactions and mass-action rate laws. We map the nonlinear algebraic problem of finding steady states onto a linear problem in a higher dimensional space. We show that the linearized version of the steady state equations obeys the linear conservation laws of the original CRN. We identify two classes of problems for which complete, minimally parameterized solutions may be obtained using only the machinery of linear systems and a judicious choice of the variables used as free parameters. We exemplify our method, providing explicit formulae, on CRN describing signal initiation of two important types of RTK receptor-ligand systems, VEGF and EGF-ErbB1.

doi:10.1109/TCBB.2013.41

PMCID: PMC4090023
PMID: 24334389

Chemical reaction networks; cell signaling; VEGF; EGF; linear conservation laws; analytical solution; bilinear systems; minimally parameterized solutions

Modern molecular biology has always been a great source of inspiration for computational science. Half a century ago, the challenge from understanding macromolecular dynamics has led the way for computations to be part of the tool set to study molecular biology. Twenty-five years ago, the demand from genome science has inspired an entire generation of computer scientists with an interest in discrete mathematics to join the field that is now called bioinformatics. In this paper, we shall lay out a new mathematical theory for dynamics of biochemical reaction systems in a small volume (i.e., mesoscopic) in terms of a stochastic, discrete-state continuous-time formulation, called the chemical master equation (CME). Similar to the wavefunction in quantum mechanics, the dynamically changing probability landscape associated with the state space provides a fundamental characterization of the biochemical reaction system. The stochastic trajectories of the dynamics are best known through the simulations using the Gillespie algorithm. In contrast to the Metropolis algorithm, this Monte Carlo sampling technique does not follow a process with detailed balance. We shall show several examples how CMEs are used to model cellular biochemical systems. We shall also illustrate the computational challenges involved: multiscale phenomena, the interplay between stochasticity and nonlinearity, and how macroscopic determinism arises from mesoscopic dynamics. We point out recent advances in computing solutions to the CME, including exact solution of the steady state landscape and stochastic differential equations that offer alternatives to the Gilespie algorithm. We argue that the CME is an ideal system from which one can learn to understand “complex behavior” and complexity theory, and from which important biological insight can be gained.

doi:10.1007/s11390-010-9312-6

PMCID: PMC4079062

biochemical networks; cellular signaling; epigenetics; master equation; nonlinear reactions; stochastic modeling

Background

Computing the long term behavior of regulatory and signaling networks is critical in understanding how biological functions take place in organisms. Steady states of these networks determine the activity levels of individual entities in the long run. Identifying all the steady states of these networks is difficult due to the state space explosion problem.

Methodology

In this paper, we propose a method for identifying all the steady states of Boolean regulatory and signaling networks accurately and efficiently. We build a mathematical model that allows pruning a large portion of the state space quickly without causing any false dismissals. For the remaining state space, which is typically very small compared to the whole state space, we develop a randomized traversal method that extracts the steady states. We estimate the number of steady states, and the expected behavior of individual genes and gene pairs in steady states in an online fashion. Also, we formulate a stopping criterion that terminates the traversal as soon as user supplied percentage of the results are returned with high confidence.

Conclusions

This method identifies the observed steady states of boolean biological networks computationally. Our algorithm successfully reported the G1 phases of both budding and fission yeast cell cycles. Besides, the experiments suggest that this method is useful in identifying co-expressed genes as well. By analyzing the steady state profile of Hedgehog network, we were able to find the highly co-expressed gene pair GL1-SMO together with other such pairs.

Availability

Source code of this work is available at http://bioinformatics.cise.ufl.edu/palSteady.html twocolumnfalse]

doi:10.1371/journal.pone.0007992

PMCID: PMC2779454
PMID: 19956604

A theory for an non-equilibrium phase transition in a driven biochemical network is presented. The theory is based on the chemical master equation (CME) formulation of mesoscopic biochemical reactions and the mathematical method of large deviations. The large deviations theory provides an analytical tool connecting the macroscopic multi-stability of an open chemical system with the multi-scale dynamics of its mesoscopic counterpart. It shows a corresponding non-equilibrium phase transition among multiple stochastic attractors. As an example, in the canonical phosphorylation–dephosphorylation system with feedback that exhibits bistability, we show that the non-equilibrium steady-state (NESS) phase transition has all the characteristics of classic equilibrium phase transition: Maxwell construction, a discontinuous first-derivative of the ‘free energy function’, Lee–Yang's zero for a generating function and a critical point that matches the cusp in nonlinear bifurcation theory. To the biochemical system, the mathematical analysis suggests three distinct timescales and needed levels of description. They are (i) molecular signalling, (ii) biochemical network nonlinear dynamics, and (iii) cellular evolution. For finite mesoscopic systems such as a cell, motions associated with (i) and (iii) are stochastic while that with (ii) is deterministic. Both (ii) and (iii) are emergent properties of a dynamic biochemical network.

doi:10.1098/rsif.2010.0202

PMCID: PMC3024822
PMID: 20466813

non-equilibrium phase transition; bistability; chemical master equation; large deviation; Maxwell construction

A defining characteristic of living cells is the ability to respond dynamically to external stimuli while maintaining homeostasis under resting conditions. Capturing both of these features in a single kinetic model is difficult because the model must be able to reproduce both behaviors using the same set of molecular components. Here, we show how combining small, well-defined steady-state networks provides an efficient means of constructing large-scale kinetic models that exhibit realistic resting and dynamic behaviors. By requiring each kinetic module to be homeostatic (at steady state under resting conditions), the method proceeds by (i) computing steady-state solutions to a system of ordinary differential equations for each module, (ii) applying principal component analysis to each set of solutions to capture the steady-state solution space of each module network, and (iii) combining optimal search directions from all modules to form a global steady-state space that is searched for accurate simulation of the time-dependent behavior of the whole system upon perturbation. Importantly, this stepwise approach retains the nonlinear rate expressions that govern each reaction in the system and enforces constraints on the range of allowable concentration states for the full-scale model. These constraints not only reduce the computational cost of fitting experimental time-series data but can also provide insight into limitations on system concentrations and architecture. To demonstrate application of the method, we show how small kinetic perturbations in a modular model of platelet P2Y1 signaling can cause widespread compensatory effects on cellular resting states.

Author Summary

Cells respond to extracellular signals through a complex coordination of interacting molecular components. Computational models can serve as powerful tools for prediction and analysis of signaling systems, but constructing large models typically requires extensive experimental datasets and computation. To facilitate the construction of complex signaling models, we present a strategy in which the models are built in a stepwise fashion, beginning with small “resting” networks that are combined to form larger models with complex time-dependent behaviors. Interestingly, we found that only a minor fraction of potential model configurations were compatible with resting behavior in an example signaling system. These reduced sets of configurations were used to limit the search for more complicated solutions that also captured the dynamic behavior of the system. Using an example model constructed by this approach, we show how a cell's resting behavior adjusts to changes in the kinetic rate processes of the system. This strategy offers a general and biologically intuitive framework for building large-scale kinetic models of steady-state cellular systems and their dynamics.

doi:10.1371/journal.pcbi.1000298

PMCID: PMC2637974
PMID: 19266013

Background

It is well known that the deterministic dynamics of biochemical reaction networks can be more easily studied if timescale separation conditions are invoked (the quasi-steady-state assumption). In this case the deterministic dynamics of a large network of elementary reactions are well described by the dynamics of a smaller network of effective reactions. Each of the latter represents a group of elementary reactions in the large network and has associated with it an effective macroscopic rate law. A popular method to achieve model reduction in the presence of intrinsic noise consists of using the effective macroscopic rate laws to heuristically deduce effective probabilities for the effective reactions which then enables simulation via the stochastic simulation algorithm (SSA). The validity of this heuristic SSA method is a priori doubtful because the reaction probabilities for the SSA have only been rigorously derived from microscopic physics arguments for elementary reactions.

Results

We here obtain, by rigorous means and in closed-form, a reduced linear Langevin equation description of the stochastic dynamics of monostable biochemical networks in conditions characterized by small intrinsic noise and timescale separation. The slow-scale linear noise approximation (ssLNA), as the new method is called, is used to calculate the intrinsic noise statistics of enzyme and gene networks. The results agree very well with SSA simulations of the non-reduced network of elementary reactions. In contrast the conventional heuristic SSA is shown to overestimate the size of noise for Michaelis-Menten kinetics, considerably under-estimate the size of noise for Hill-type kinetics and in some cases even miss the prediction of noise-induced oscillations.

Conclusions

A new general method, the ssLNA, is derived and shown to correctly describe the statistics of intrinsic noise about the macroscopic concentrations under timescale separation conditions. The ssLNA provides a simple and accurate means of performing stochastic model reduction and hence it is expected to be of widespread utility in studying the dynamics of large noisy reaction networks, as is common in computational and systems biology.

doi:10.1186/1752-0509-6-39

PMCID: PMC3532178
PMID: 22583770

Reaction networks are systems in which the populations of a finite number of species evolve through predefined interactions. Such networks are found as modeling tools in many biological disciplines such as biochemistry, ecology, epidemiology, immunology, systems biology and synthetic biology. It is now well-established that, for small population sizes, stochastic models for biochemical reaction networks are necessary to capture randomness in the interactions. The tools for analyzing such models, however, still lag far behind their deterministic counterparts. In this paper, we bridge this gap by developing a constructive framework for examining the long-term behavior and stability properties of the reaction dynamics in a stochastic setting. In particular, we address the problems of determining ergodicity of the reaction dynamics, which is analogous to having a globally attracting fixed point for deterministic dynamics. We also examine when the statistical moments of the underlying process remain bounded with time and when they converge to their steady state values. The framework we develop relies on a blend of ideas from probability theory, linear algebra and optimization theory. We demonstrate that the stability properties of a wide class of biological networks can be assessed from our sufficient theoretical conditions that can be recast as efficient and scalable linear programs, well-known for their tractability. It is notably shown that the computational complexity is often linear in the number of species. We illustrate the validity, the efficiency and the wide applicability of our results on several reaction networks arising in biochemistry, systems biology, epidemiology and ecology. The biological implications of the results as well as an example of a non-ergodic biological network are also discussed.

Author Summary

In many biological disciplines, computational modeling of interaction networks is the key for understanding biological phenomena. Such networks are traditionally studied using deterministic models. However, it has been recently recognized that when the populations are small in size, the inherent random effects become significant and to incorporate them, a stochastic modeling paradigm is necessary. Hence, stochastic models of reaction networks have been broadly adopted and extensively used. Such models, for instance, form a cornerstone for studying heterogeneity in clonal cell populations. In biological applications, one is often interested in knowing the long-term behavior and stability properties of reaction networks even with incomplete knowledge of the model parameters. However for stochastic models, no analytical tools are known for this purpose, forcing many researchers to use a simulation-based approach, which is highly unsatisfactory. To address this issue, we develop a theoretical and computational framework for determining the long-term behavior and stability properties for stochastic reaction networks. Our approach is based on a mixture of ideas from probability theory, linear algebra and optimization theory. We illustrate the broad applicability of our results by considering examples from various biological areas. The biological implications of our results are discussed as well.

doi:10.1371/journal.pcbi.1003669

PMCID: PMC4072526
PMID: 24968191