Scaling methods for accelerating kinetic Monte Carlo simulations of chemical reaction networks.

Various kinetic Monte Carlo algorithms become inefficient when some of the population sizes in a system are large, which gives rise to a large number of reaction events per unit time. Here, we present a new acceleration algorithm based on adaptive and heterogeneous scaling of reaction rates and stoichiometric coefficients. The algorithm is conceptually related to the commonly used idea of accelerating a stochastic simulation by considering a subvolume λΩ (0 < λ < 1) within a system of interest, which reduces the number of reaction events per unit time occurring in a simulation by a factor 1/λ at the cost of greater error in unbiased estimates of first moments and biased overestimates of second moments. Our new approach offers two unique benefits. First, scaling is adaptive and heterogeneous, which eliminates the pitfall of overaggressive scaling. Second, there is no need for an a priori classification of populations as discrete or continuous (as in a hybrid method), which is problematic when discreteness of a chemical species changes during a simulation. The method requires specification of only a single algorithmic parameter, Nc, a global critical population size above which populations are effectively scaled down to increase simulation efficiency. The method, which we term partial scaling, is implemented in the open-source BioNetGen software package. We demonstrate that partial scaling can significantly accelerate simulations without significant loss of accuracy for several published models of biological systems. These models characterize activation of the mitogen-activated protein kinase ERK, prion protein aggregation, and T-cell receptor signaling.

Related papers:

title: Generalizing Gillespie’s Direct Method to Enable Network-Free Simulations
abstract: Gillespie’s direct method for stochastic simulation of chemical kinetics is a staple of computational systems biology research. However, the algorithm requires explicit enumeration of all reactions and all chemical species that may arise in the system. In many cases, this is not feasible due to the combinatorial explosion of reactions and species in biological networks. Rule-based modeling frameworks provide a way to exactly represent networks containing such combinatorial complexity, and generalizations of Gillespie’s direct method have been developed as simulation engines for rule-based modeling languages. Here, we provide both a high-level description of the algorithms underlying the simulation engines, termed network-free simulation algorithms, and how they have been applied in systems biology research. We also define a generic rule-based modeling framework and describe a number of technical details required for adapting Gillespie’s direct method for network-free simulation. Finally, we briefly discuss potential avenues for advancing network-free simulation and the role they continue to play in modeling dynamical systems in biology.

title: Rule‐based modeling: a computational approach for studying biomolecular site dynamics in cell signaling systems
abstract: Rule‐based modeling was developed to address the limitations of traditional approaches for modeling chemical kinetics in cell signaling systems. These systems consist of multiple interacting biomolecules (e.g., proteins), which themselves consist of multiple parts (e.g., domains, linear motifs, and sites of phosphorylation). Consequently, biomolecules that mediate information processing generally have the potential to interact in multiple ways, with the number of possible complexes and posttranslational modification states tending to grow exponentially with the number of binary interactions considered. As a result, only large reaction networks capture all possible consequences of the molecular interactions that occur in a cell signaling system, which is problematic because traditional modeling approaches for chemical kinetics (e.g., ordinary differential equations) require explicit network specification. This problem is circumvented through representation of interactions in terms of local rules. With this approach, network specification is implicit and model specification is concise. Concise representation results in a coarse graining of chemical kinetics, which is introduced because all reactions implied by a rule inherit the rate law associated with that rule. Coarse graining can be appropriate if interactions are modular, and the coarseness of a model can be adjusted as needed. Rules can be specified using specialized model‐specification languages, and recently developed tools designed for specification of rule‐based models allow one to leverage powerful software engineering capabilities. A rule‐based model comprises a set of rules, which can be processed by general‐purpose simulation and analysis tools to achieve different objectives (e.g., to perform either a deterministic or stochastic simulation).  2014, 6:13–36. doi: 10.1002/wsbm.1245

title: Kinetic Monte Carlo method for rule-based modeling of biochemical networks.
abstract: We present a kinetic Monte Carlo method for simulating chemical transformations specified by reaction rules, which can be viewed as generators of chemical reactions, or equivalently, definitions of reaction classes. A rule identifies the molecular components involved in a transformation, how these components change, conditions that affect whether a transformation occurs, and a rate law. The computational cost of the method, unlike conventional simulation approaches, is independent of the number of possible reactions, which need not be specified in advance or explicitly generated in a simulation. To demonstrate the method, we apply it to study the kinetics of multivalent ligand-receptor interactions. We expect the method will be useful for studying cellular signaling systems and other physical systems involving aggregation phenomena.

title: Efficient modeling, simulation and coarse-graining of biological complexity with NFsim
abstract: Managing the overwhelming numbers of molecular states and interactions is a fundamental obstacle to building predictive models of biological systems. Here we introduce the Network-Free Stochastic Simulator (NFsim), a general-purpose modeling platform that overcomes the combinatorial nature of molecular interactions. Unlike standard simulators that represent molecular species as variables in equations, NFsim uses a biologically intuitive representation: objects with binding and modification sites acted on by reaction rules. During simulations, rules operate directly on molecular objects to produce exact stochastic results with performance that scales independently of the reaction network size. Reaction rates can be defined as arbitrary functions of molecular states to provide powerful coarse-graining capabilities, for example to merge Boolean and kinetic representations of biological networks. NFsim enables researchers to simulate many biological systems that were previously inaccessible to general-purpose software, as we illustrate with models of immune system signaling, microbial signaling, cytoskeletal assembly and oscillating gene expression.

title: Specification, annotation, visualization and simulation of a large rule-based model for ERBB receptor signaling
abstract: BackgroundMathematical/computational models are needed to understand cell signaling networks, which are complex. Signaling proteins contain multiple functional components and multiple sites of post-translational modification. The multiplicity of components and sites of modification ensures that interactions among signaling proteins have the potential to generate myriad protein complexes and post-translational modification states. As a result, the number of chemical species that can be populated in a cell signaling network, and hence the number of equations in an ordinary differential equation model required to capture the dynamics of these species, is prohibitively large. To overcome this problem, the rule-based modeling approach has been developed for representing interactions within signaling networks efficiently and compactly through coarse-graining of the chemical kinetics of molecular interactions.ResultsHere, we provide a demonstration that the rule-based modeling approach can be used to specify and simulate a large model for ERBB receptor signaling that accounts for site-specific details of protein-protein interactions. The model is considered large because it corresponds to a reaction network containing more reactions than can be practically enumerated. The model encompasses activation of ERK and Akt, and it can be simulated using a network-free simulator, such as NFsim, to generate time courses of phosphorylation for 55 individual serine, threonine, and tyrosine residues. The model is annotated and visualized in the form of an extended contact map.ConclusionsWith the development of software that implements novel computational methods for calculating the dynamics of large-scale rule-based representations of cellular signaling networks, it is now possible to build and analyze models that include a significant fraction of the protein interactions that comprise a signaling network, with incorporation of the site-specific details of the interactions. Modeling at this level of detail is important for understanding cellular signaling.

title: Approximate accelerated stochastic simulation of chemically reacting systems
abstract: The stochastic simulation algorithm (SSA) is an essentially exact procedure for numerically simulating the time evolution of a well-stirred chemically reacting system. Despite recent major improvements in the efficiency of the SSA, its drawback remains the great amount of computer time that is often required to simulate a desired amount of system time. Presented here is the “τ-leap” method, an approximate procedure that in some circumstances can produce significant gains in simulation speed with acceptable losses in accuracy. Some primitive strategies for control parameter selection and error mitigation for the τ-leap method are described, and simulation results for two simple model systems are exhibited. With further refinement, the τ-leap method should provide a viable way of segueing from the exact SSA to the approximate chemical Langevin equation, and thence to the conventional deterministic reaction rate equation, as the system size becomes larger.

title: Stochastic simulation of chemical kinetics.
abstract: Stochastic chemical kinetics describes the time evolution of a well-stirred chemically reacting system in a way that takes into account the fact that molecules come in whole numbers and exhibit some degree of randomness in their dynamical behavior. Researchers are increasingly using this approach to chemical kinetics in the analysis of cellular systems in biology, where the small molecular populations of only a few reactant species can lead to deviations from the predictions of the deterministic differential equations of classical chemical kinetics. After reviewing the supporting theory of stochastic chemical kinetics, I discuss some recent advances in methods for using that theory to make numerical simulations. These include improvements to the exact stochastic simulation algorithm (SSA) and the approximate explicit tau-leaping procedure, as well as the development of two approximate strategies for simulating systems that are dynamically stiff: implicit tau-leaping and the slow-scale SSA.

title: Breakdown of fast-slow analysis in an excitable system with channel noise.
abstract: We consider a stochastic version of an excitable system based on the Morris-Lecar model of a neuron, in which the noise originates from stochastic sodium and potassium ion channels opening and closing. One can analyze neural excitability in the deterministic model by using a separation of time scales involving a fast voltage variable and a slow recovery variable, which represents the fraction of open potassium channels. In the stochastic setting, spontaneous excitation is initiated by ion channel noise. If the recovery variable is constant during initiation, the spontaneous activity rate can be calculated using Kramer's rate theory. The validity of this assumption in the stochastic model is examined using a systematic perturbation analysis. We find that, in most physically relevant cases, this assumption breaks down, requiring an alternative to Kramer's theory for excitable systems with one deterministic fixed point. We also show that an exit time problem can be formulated in an excitable system by considering maximum likelihood trajectories of the stochastic process.

title: Transcriptional Bursting Diversifies the Behaviour of a Toggle Switch: Hybrid Simulation of Stochastic Gene Expression
abstract: Hybrid models for gene expression combine stochastic and deterministic representations of the underlying biophysical mechanisms. According to one of the simplest hybrid formalisms, protein molecules are produced in randomly occurring bursts of a randomly distributed size while they are degraded deterministically. Here, we use this particular formalism to study two key regulatory motifs—the autoregulation loop and the toggle switch. The distribution of burst times is determined and used as a basis for the development of exact simulation algorithms for gene expression dynamics. For the autoregulation loop, the simulations are compared to an analytic solution of a master equation. Simulations of the toggle switch reveal a number of qualitatively distinct scenarios with implications for the modelling of cell-fate selection.

title: Stochastic hybrid model of spontaneous dendritic NMDA spikes.
abstract: Following recent advances in imaging techniques and methods of dendritic stimulation, active voltage spikes have been observed in thin dendritic branches of excitatory pyramidal neurons, where the majority of synapses occur. The generation of these dendritic spikes involves both Na(+) ion channels and M-methyl-D-aspartate receptor (NMDAR) channels. During strong stimulation of a thin dendrite, the resulting high levels of glutamate, the main excitatory neurotransmitter in the central nervous system and an NMDA agonist, modify the current-voltage (I-V) characteristics of an NMDAR so that it behaves like a voltage-gated Na(+) channel. Hence, the NMDARs can fire a regenerative dendritic spike, just as Na(+) channels support the initiation of an action potential following membrane depolarization. However, the duration of the dendritic spike is of the order 100 ms rather than 1 ms, since it involves slow unbinding of glutamate from NMDARs rather than activation of hyperpolarizing K(+) channels. It has been suggested that dendritic NMDA spikes may play an important role in dendritic computations and provide a cellular substrate for short-term memory. In this paper, we consider a stochastic, conductance-based model of dendritic NMDA spikes, in which the noise originates from the stochastic opening and closing of a finite number of Na(+) and NMDA receptor ion channels. The resulting model takes the form of a stochastic hybrid system, in which membrane voltage evolves according to a piecewise deterministic dynamics that is coupled to a jump Markov process describing the opening and closing of the ion channels. We formulate the noise-induced initiation and termination of a dendritic spike in terms of a first-passage time problem, under the assumption that glutamate unbinding is negligible, which we then solve using a combination of WKB methods and singular perturbation theory. Using a stochastic phase-plane analysis we then extend our analysis to take proper account of the combined effects of glutamate unbinding and noise on the termination of a spike.

title: Bursting noise in gene expression dynamics: linking microscopic and mesoscopic models
abstract: The dynamics of short-lived mRNA results in bursts of protein production in gene regulatory networks. We investigate the propagation of bursting noise between different levels of mathematical modelling and demonstrate that conventional approaches based on diffusion approximations can fail to capture bursting noise. An alternative coarse-grained model, the so-called piecewise deterministic Markov process (PDMP), is seen to outperform the diffusion approximation in biologically relevant parameter regimes. We provide a systematic embedding of the PDMP model into the landscape of existing approaches, and we present analytical methods to calculate its stationary distribution and switching frequencies.

title: Stochastic switching in biology: From genotype to phenotype
abstract: There has been a resurgence of interest in non-equilibrium stochastic processes in recent years, driven in part by the observation that the number of molecules (genes, mRNA, proteins) involved in gene expression are often of order 1–1000. This means that deterministic mass-action kinetics tends to break down, and one needs to take into account the discrete, stochastic nature of biochemical reactions. One of the major consequences of molecular noise is the occurrence of stochastic biological switching at both the genotypic and phenotypic levels. For example, individual gene regulatory networks can switch between graded and binary responses, exhibit translational/ transcriptional bursting, and support metastability (noise-induced switching between states that are stable in the deterministic limit). If random switching persists at the phenotypic level then this can confer certain advantages to cell populations growing in a changing environment, as exemplified by bacterial persistence in response to antibiotics. Gene expression at the single-cell level can also be regulated by changes in cell density at the population level, a process known as quorum sensing. In contrast to noise-driven phenotypic switching, the switching mechanism in quorum sensing is stimulus-driven and thus noise tends to have a detrimental effect. A common approach to modeling stochastic gene expression is to assume a large but finite system and to approximate the discrete processes by continuous processes using a systemsize expansion. However, there is a growing need to have some familiarity with the theory of stochastic processes that goes beyond the standard topics of chemical master equations, the system-size expansion, Langevin equations  and the Fokker–Planck equation. Examples include stochastic hybrid systems (piecewise deterministic Markov processes), large deviations and the Wentzel–Kramers–Brillouin (WKB) method, adiabatic reductions, and queuing/renewal theory. The major aim of this review is to provide a selfcontained survey of these mathematical methods, mainly within the context of biological switching processes at both the genotypic and phenotypic levels. P C Bressloff Stochastic switching in biology: from genotype to phenotype Printed in the UK 133001

title: Feynman-Kac formula for stochastic hybrid systems.
abstract: We derive a Feynman-Kac formula for functionals of a stochastic hybrid system evolving according to a piecewise deterministic Markov process. We first derive a stochastic Liouville equation for the moment generator of the stochastic functional, given a particular realization of the underlying discrete Markov process; the latter generates transitions between different dynamical equations for the continuous process. We then analyze the stochastic Liouville equation using methods recently developed for diffusion processes in randomly switching environments. In particular, we obtain dynamical equations for the moment generating function, averaged with respect to realizations of the discrete Markov process. The resulting Feynman-Kac formula takes the form of a differential Chapman-Kolmogorov equation. We illustrate the theory by calculating the occupation time for a one-dimensional velocity jump process on the infinite or semi-infinite real line. Finally, we present an alternative derivation of the Feynman-Kac formula based on a recent path-integral formulation of stochastic hybrid systems.

title: Efficient analysis of stochastic gene dynamics in the non-adiabatic regime using piecewise deterministic Markov processes
abstract: Single-cell experiments show that gene expression is stochastic and bursty, a feature that can emerge from slow switching between promoter states with different activities. In addition to slow chromatin and/or DNA looping dynamics, one source of long-lived promoter states is the slow binding and unbinding kinetics of transcription factors to promoters, i.e. the non-adiabatic binding regime. Here, we introduce a simple analytical framework, known as a piecewise deterministic Markov process (PDMP), that accurately describes the stochastic dynamics of gene expression in the non-adiabatic regime. We illustrate the utility of the PDMP on a non-trivial dynamical system by analysing the properties of a titration-based oscillator in the non-adiabatic limit. We first show how to transform the underlying chemical master equation into a PDMP where the slow transitions between promoter states are stochastic, but whose rates depend upon the faster deterministic dynamics of the transcription factors regulated by these promoters. We show that the PDMP accurately describes the observed periods of stochastic cycles in activator and repressor-based titration oscillators. We then generalize our PDMP analysis to more complicated versions of titration-based oscillators to explain how multiple binding sites lengthen the period and improve coherence. Last, we show how noise-induced oscillation previously observed in a titration-based oscillator arises from non-adiabatic and discrete binding events at the promoter site.

title: Stochastic Differential Equations With Markovian Switching
abstract: This textbook provides the first systematic presentation of the theory of stochastic differential equations with Markovian switching. It presents the basic principles at an introductory level but emphasizes current advanced level research trends. The material takes into account all the features of Ito equations, Markovian switching, interval systems and time-lag. The theory developed is applicable in different and complicated situations in many branches of science and industry.

title: Stochasticity in transcriptional regulation: origins, consequences, and mathematical representations.
abstract: Transcriptional regulation is an inherently noisy process. The origins of this stochastic behavior can be traced to the random transitions among the discrete chemical states of operators that control the transcription rate and to finite number fluctuations in the biochemical reactions for the synthesis and degradation of transcripts. We develop stochastic models to which these random reactions are intrinsic and a series of simpler models derived explicitly from the first as approximations in different parameter regimes. This innate stochasticity can have both a quantitative and qualitative impact on the behavior of gene-regulatory networks. We introduce a natural generalization of deterministic bifurcations for classification of stochastic systems and show that simple noisy genetic switches have rich bifurcation structures; among them, bifurcations driven solely by changing the rate of operator fluctuations even as the underlying deterministic system remains unchanged. We find stochastic bistability where the deterministic equations predict monostability and vice-versa. We derive and solve equations for the mean waiting times for spontaneous transitions between quasistable states in these switches.

title: Rule-based modeling of biochemical systems with BioNetGen.
abstract: Rule-based modeling involves the representation of molecules as structured objects and molecular interactions as rules for transforming the attributes of these objects. The approach is notable in that it allows one to systematically incorporate site-specific details about protein-protein interactions into a model for the dynamics of a signal-transduction system, but the method has other applications as well, such as following the fates of individual carbon atoms in metabolic reactions. The consequences of protein-protein interactions are difficult to specify and track with a conventional modeling approach because of the large number of protein phosphoforms and protein complexes that these interactions potentially generate. Here, we focus on how a rule-based model is specified in the BioNetGen language (BNGL) and how a model specification is analyzed using the BioNetGen software tool. We also discuss new developments in rule-based modeling that should enable the construction and analyses of comprehensive models for signal transduction pathways and similarly large-scale models for other biochemical systems.

title: Relaxation oscillations and hierarchy of feedbacks in MAPK signaling
abstract: We formulated a computational model for a MAPK signaling cascade downstream of the EGF receptor to investigate how interlinked positive and negative feedback loops process EGF signals into ERK pulses of constant amplitude but dose-dependent duration and frequency. A positive feedback loop involving RAS and SOS, which leads to bistability and allows for switch-like responses to inputs, is nested within a negative feedback loop that encompasses RAS and RAF, MEK, and ERK that inhibits SOS via phosphorylation. This negative feedback, operating on a longer time scale, changes switch-like behavior into oscillations having a period of 1 hour or longer. Two auxiliary negative feedback loops, from ERK to MEK and RAF, placed downstream of the positive feedback, shape the temporal ERK activity profile but are dispensable for oscillations. Thus, the positive feedback introduces a hierarchy among negative feedback loops, such that the effect of a negative feedback depends on its position with respect to the positive feedback loop. Furthermore, a combination of the fast positive feedback involving slow-diffusing membrane components with slower negative feedbacks involving faster diffusing cytoplasmic components leads to local excitation/global inhibition dynamics, which allows the MAPK cascade to transmit paracrine EGF signals into spatially non-uniform ERK activity pulses.

title: Dynamics of the nucleated polymerization model of prion replication.
abstract: The disease process for transmissible spongiform encephalopathies (TSEs), in one way or another, involves the conversion of a predominantly alpha-helical normal host-coded prion protein (PrP(C)) to an abnormally folded (predominantly beta sheet) protease resistant isoform (PrP(Sc)). Several alternative mechanisms have been proposed for this auto-catalytic process. Here the dynamical behavior of one of these models, the nucleated polymerization model, is studied by Monte Carlo discrete-event simulation of the explicit conversion reactions. These simulations demonstrate the characteristic dynamical behavior of this model for prion replication. Using estimates for the reaction rates and concentrations, time courses are estimated for concentration of PrP(Sc), PrP(Sc) aggregates, and PrP(C) as well as size distributions for the aggregates. The implications of these dynamics on protein misfolding cyclic amplification (PMCA) is discussed.

title: Stochastic effects and bistability in T cell receptor signaling.
abstract: The stochastic dynamics of T cell receptor (TCR) signaling are studied using a mathematical model intended to capture kinetic proofreading (sensitivity to ligand-receptor binding kinetics) and negative and positive feedback regulation mediated, respectively, by the phosphatase SHP1 and the MAP kinase ERK. The model incorporates protein-protein interactions involved in initiating TCR-mediated cellular responses and reproduces several experimental observations about the behavior of TCR signaling, including robust responses to as few as a handful of ligands (agonist peptide-MHC complexes on an antigen-presenting cell), distinct responses to ligands that bind TCR with different lifetimes, and antagonism. Analysis of the model indicates that TCR signaling dynamics are marked by significant stochastic fluctuations and bistability, which is caused by the competition between the positive and negative feedbacks. Stochastic fluctuations are such that single-cell trajectories differ qualitatively from the trajectory predicted in the deterministic approximation of the dynamics. Because of bistability, the average of single-cell trajectories differs markedly from the deterministic trajectory. Bistability combined with stochastic fluctuations allows for switch-like responses to signals, which may aid T cells in making committed cell-fate decisions.

title: BioNetGen 2.2: advances in rule-based modeling
abstract: : BioNetGen is an open-source software package for rule-based modeling of complex biochemical systems. Version 2.2 of the software introduces numerous new features for both model specification and simulation. Here, we report on these additions, discussing how they facilitate the construction, simulation and analysis of larger and more complex models than previously possible.



Stable BioNetGen releases (Linux, Mac OS/X and Windows), with documentation, are available at http://bionetgen.org Source code is available at http://github.com/RuleWorld/bionetgen CONTACT: bionetgen.help@gmail.comSupplementary information: Supplementary data are available at Bioinformatics online.

title: Exact Stochastic Simulation of Coupled Chemical Reactions
abstract: There are two formalisms for mathematically describing the time behavior of a spatially homogeneous chemical system: The deterministic approach regards the time evolution as a continuous, wholly predictable process which is governed by a set of coupled, ordinary differential equations (the “reaction-rate equations”); the stochastic approach regards the time evolution as a kind of random-walk process which is governed by a single differential-difference equation (the “master equation”). Fairly simple kinetic theory arguments show that the stochastic formulation of chemical kinetics has a firmer physical basis than the deterministic formulation, but unfortunately the stochastic master equation is often mathematically intractable. There is, however, a way to make exact numerical calculations within the framework of the stochastic formulation without having to deal with the master equation directly. It is a relatively simple digital computer algorithm which uses a rigorously derived Monte Carlo procedure to numerically simulate the time evolution of the given chemical system. Like the master equation, this “stochastic simulation algorithm” correctly accounts for the inherent fluctuations and correlations that are necessarily ignored in the deterministic formulation. In addition, unlike most procedures for numerically solving the deterministic reaction-rate equations, this algorithm never approximates infinitesimal time increments df by finite time steps At. The feasibility and utility of the simulation algorithm are demonstrated by applying it to several well-known model chemical systems, including the Lotka model, the Brusselator, and the Oregonator.

title: a Power Series Expansion of the Master Equation
abstract: In order to solve the master equation by a systematic approximation method, an expansion in powers of some parameter is needed. The appropriate parameter is the reciprocal size of the system, defin...

title: Biological Modeling and Simulation: A Survey of Practical Models, Algorithms, and Numerical Methods
abstract: There are many excellent computational biology resources now available for learning about methods that have been developed to address specific biological systems, but comparatively little attention has been paid to training aspiring computational biologists to handle new and unanticipated problems. This text is intended to fill that gap by teaching students how to reason about developing formal mathematical models of biological systems that are amenable to computational analysis. It collects in one place a selection of broadly useful models, algorithms, and theoretical analysis tools normally found scattered among many other disciplines. It thereby gives the aspiring student a bag of tricks that will serve him or her well in modeling problems drawn from numerous subfields of biology. These techniques are taught from the perspective of what the practitioner needs to know to use them effectively, supplemented with references for further reading on more advanced use of each method covered. The text, which grew out of a class taught at Carnegie Mellon University, covers models for optimization, simulation and sampling, and parameter tuning. These topics provide a general framework for learning how to formulate mathematical models of biological systems, what techniques are available to work with these models, and how to fit the models to particular systems. Their application is illustrated by many examples drawn from a variety of biological disciplines and several extended case studies that show how the methods described have been applied to real problems in biology.

title: Efficient Exact Stochastic Simulation of Chemical Systems with Many Species and Many Channels
abstract: There are two fundamental ways to view coupled systems of chemical equations:  as continuous, represented by differential equations whose variables are concentrations, or as discrete, represented by stochastic processes whose variables are numbers of molecules. Although the former is by far more common, systems with very small numbers of molecules are important in some applications (e.g., in small biological cells or in surface processes). In both views, most complicated systems with multiple reaction channels and multiple chemical species cannot be solved analytically. There are exact numerical simulation methods to simulate trajectories of discrete, stochastic systems, (methods that are rigorously equivalent to the Master Equation approach) but these do not scale well to systems with many reaction pathways. This paper presents the Next Reaction Method, an exact algorithm to simulate coupled chemical reactions that is also efficient:  it (a) uses only a single random number per simulation event, and (b) ...

title: Self-consistent proteomic field theory of stochastic gene switches.
abstract: We present a self-consistent field approximation approach to the problem of the genetic switch composed of two mutually repressing/activating genes. The protein and DNA state dynamics are treated stochastically and on an equal footing. In this approach the mean influence of the proteomic cloud created by one gene on the action of another is self-consistently computed. Within this approximation a broad range of stochastic genetic switches may be solved exactly in terms of finding the probability distribution and its moments. A much larger class of problems, such as genetic networks and cascades, also remain exactly solvable with this approximation. We discuss, in depth, certain specific types of basic switches used by biological systems and compare their behavior to the expectation for a deterministic switch.

title: Optimal path to epigenetic switching.
abstract: We use large deviation methods to calculate rates of noise-induced transitions between states in multistable genetic networks. We analyze a synthetic biochemical circuit, the toggle switch, and compare the results to those obtained from a numerical solution of the master equation.

title: Determining the stability of genetic switches: explicitly accounting for mRNA noise.
abstract: Cells use genetic switches to shift between alternate gene-expression states, e.g., to adapt to new environments or to follow a developmental pathway. Here, we study the dynamics of switching in a generic-feedback on-off switch. Unlike protein-only models, we explicitly account for stochastic fluctuations of mRNA, which have a dramatic impact on switch dynamics. Employing the WKB theory to treat the underlying chemical master equations, we obtain accurate results for the quasistationary distributions of mRNA and protein copy numbers and for the mean switching time, starting from either state. Our analytical results agree well with Monte Carlo simulations. Importantly, one can use the approach to study the effect of varying biological parameters on switch stability.

title: Stability and multiattractor dynamics of a toggle switch based on a two-stage model of stochastic gene expression
abstract: A toggle switch consists of two genes that mutually repress each other. This regulatory motif is active during cell differentiation and is thought to act as a memory device, being able to choose and maintain cell fate decisions. Commonly, this switch has been modeled in a deterministic framework where transcription and translation are lumped together. In this description, bistability occurs for transcription factor cooperativity, whereas autoactivation leads to a tristable system with an additional undecided state. In this contribution, we study the stability and dynamics of a two-stage gene expression switch within a probabilistic framework inspired by the properties of the Pu/Gata toggle switch in myeloid progenitor cells. We focus on low mRNA numbers, high protein abundance, and monomeric transcription-factor binding. Contrary to the expectation from a deterministic description, this switch shows complex multiattractor dynamics without autoactivation and cooperativity. Most importantly, the four attractors of the system, which only emerge in a probabilistic two-stage description, can be identified with committed and primed states in cell differentiation. To begin, we study the dynamics of the system and infer the mechanisms that move the system between attractors using both the quasipotential and the probability flux of the system. Next, we show that the residence times of the system in one of the committed attractors are geometrically distributed. We derive an analytical expression for the parameter of the geometric distribution, therefore completely describing the statistics of the switching process and elucidate the influence of the system parameters on the residence time. Moreover, we find that the mean residence time increases linearly with the mean protein level. This scaling also holds for a one-stage scenario and for autoactivation. Finally, we study the implications of this distribution for the stability of a switch and discuss the influence of the stability on a specific cell differentiation mechanism. Our model explains lineage priming and proposes the need of either high protein numbers or long-term modifications such as chromatin remodeling to achieve stable cell fate decisions. Notably, we present a system with high protein abundance that nevertheless requires a probabilistic description to exhibit multistability, complex switching dynamics, and lineage priming.

title: Tristability in cancer-associated microRNA-TF chimera toggle switch.
abstract: Cell fate decisions during embryonic development and tumorigenesis pose a major research challenge in modern developmental and cancer biology. Binary cell fate decisions are usually regulated by gene circuits incorporating either classical toggle switches with two mutually inhibiting transcription factor (TF) genes or chimera toggle switches with a mutually inhibiting pair of microRNA (miRNA) and TF gene. These circuits can explain binary cell fate decisions. Importantly, intermediate cell types can exist during the differentiation of both stem cells and cancer cells. It has been shown that TF-TF self-activating toggle switches (SATS) can have coexistence of three metastable states (tristability), yet the role of chimera toggle switches in opening these additional states remains elusive. Here we present a generalized framework for both the TF-TF SATS and miRNA-TF chimera SATS, starting from the TF-promoter and miRNA-mRNA binding/unbinding dynamics. We show that the chimera SATSs can also have tristability. We demonstrate that the dynamics of miRNA-TF SATS is qualitatively different from that of the TF-TF SATS because the nonlinear effects of translational silencing by miRNA are distinct from those of transcriptional repression. We discuss the possible relevance of these findings to fate decisions by cancer cells.

title: Automatic generation of cellular reaction networks with Moleculizer 1.0
abstract: Accurate simulation of intracellular biochemical networks is essential to furthering our understanding of biological system behavior. The number of protein complexes and of chemical interactions among them has traditionally posed significant problems for simulation algorithms. Here we describe an approach to the exact stochastic simulation of biochemical networks that emphasizes the contribution of protein complexes to these systems. This simulation approach starts from a description of monomeric proteins and specifications for binding, unbinding and other reactions. This manageable specification is reasonably intuitive for biologists. Rather than requiring the inclusion of all possible complexes and reactions from the outset, our approach incorporates new complexes and reactions only when needed as the simulation proceeds. As a result, the simulation generates much smaller reaction networks, which can be exported to other simulators for further analysis. We apply this approach to the automatic generation of reaction systems for the study of signal transduction networks.

title: Rule-based modeling of biochemical networks
abstract: A method for the automatic generation of mathematical/computational models that account comprehensively and precisely for the full spectrum of chemical species implied by user-specified activities, potential modifications and interactions of the molecular components of biomolecules is described. A computer-implemented system that includes software was used to generate models. The software has a user interface that allows a user to generate new models and modify existing models.

title: Phosphorylation Site Dynamics of Early T-cell Receptor Signaling
abstract: In adaptive immune responses, T-cell receptor (TCR) signaling impacts multiple cellular processes and results in T-cell differentiation, proliferation, and cytokine production. Although individual protein–protein interactions and phosphorylation events have been studied extensively, we lack a systems-level understanding of how these components cooperate to control signaling dynamics, especially during the crucial first seconds of stimulation. Here, we used quantitative proteomics to characterize reshaping of the T-cell phosphoproteome in response to TCR/CD28 co-stimulation, and found that diverse dynamic patterns emerge within seconds. We detected phosphorylation dynamics as early as 5 s and observed widespread regulation of key TCR signaling proteins by 30 s. Development of a computational model pointed to the presence of novel regulatory mechanisms controlling phosphorylation of sites with central roles in TCR signaling. The model was used to generate predictions suggesting unexpected roles for the phosphatase PTPN6 (SHP-1) and shortcut recruitment of the actin regulator WAS. Predictions were validated experimentally. This integration of proteomics and modeling illustrates a novel, generalizable framework for solidifying quantitative understanding of a signaling network and for elucidating missing links.

title: BioNetFit: a fitting tool compatible with BioNetGen, NFsim and distributed computing environments
abstract: Rule-based models are analyzed with specialized simulators, such as those provided by the BioNetGen and NFsim open-source software packages. Here, we present BioNetFit, a general-purpose fitting tool that is compatible with BioNetGen and NFsim. BioNetFit is designed to take advantage of distributed computing resources. This feature facilitates fitting (i.e. optimization of parameter values for consistency with data) when simulations are computationally expensive.



BioNetFit can be used on stand-alone Mac, Windows/Cygwin, and Linux platforms and on Linux-based clusters running SLURM, Torque/PBS, or SGE. The BioNetFit source code (Perl) is freely available (http://bionetfit.nau.edu).



Supplementary data are available at Bioinformatics online.



bionetgen.help@gmail.com.

title: Multistate modeling and simulation for regulatory networks
abstract: Many protein regulatory models contain chemical species best represented as having multiple states. Such models stem from the potential for multiple levels of phosphorylation or from the formation of multiprotein complexes. We seek to support such models by augmenting an existing modeling and simulation system. Interactions between multistate species can lead to a combinatorial explosion in the potential state space. This creates a challenge when using Gillespie's stochastic simulation algorithm (SSA). Both the network-free algorithm (NFA) and various rules-based methods have been proposed to more efficiently simulate such models. We show how to further improve NFA to integrate population-based and particle-based features. We then present a population-based scheme for the stochastic simulation of rule-based models. A complexity analysis is presented comparing the proposed simulation methods. We present numerical experiments for two sample models that demonstrate the power of the proposed simulation methods.