A multi-time-scale analysis of chemical reaction networks: I. Deterministic systems

Revisiting kinetic Monte Carlo algorithms for time-dependent processes: From open-loop control to feedback control

Simulating stochastic systems with feedback control is challenging due to the complex interplay between the system's dynamics and the feedback-dependent control protocols. We present a single-step-trajectory probability analysis to time-dependent stochastic systems. Based on this analysis, we revisit several time-dependent kinetic Monte Carlo (KMC) algorithms designed for systems under open-loop-control protocols. Our analysis provides a unified alternative proof to these algorithms, summarized into a pedagogical tutorial. Moreover, with the trajectory probability analysis, we present a novel feedback-controlled KMC algorithm that accurately captures the dynamics systems controlled by an external signal based on the measurements of the system's state. Our method correctly captures the system dynamics and avoids the artificial Zeno effect that arises from incorrectly applying the direct Gillespie algorithm to feedback-controlled systems. This work provides a unified perspective on existing open-loop-control KMC algorithms and also offers a powerful and accurate tool for simulating stochastic systems with feedback control.

Universal slow dynamics of chemical reaction networks

Understanding the emergent behavior of chemical reaction networks (CRNs) is a fundamental aspect of biology and its origin from inanimate matter. A closed CRN monotonically tends to thermal equilibrium, but when it is opened to external reservoirs, a range of behaviors is possible, including transition to a new equilibrium state, a nonequilibrium state, or indefinite growth. This study shows that slowly driven CRNs are governed by the conserved quantities of the closed system, which are generally far fewer in number than the species. Considering both deterministic and stochastic dynamics, a universal slow-dynamics equation is derived with singular perturbation methods and is shown to be thermodynamically consistent. The slow dynamics is highly robust against microscopic details of the network, which may be unknown in practical situations. In particular, nonequilibrium states of realistic large CRNs can be sought without knowledge of bulk reaction rates. The framework is successfully tested against a suite of networks of increasing complexity and argued to be relevant in the treatment of open CRNs as chemical machines.

Multi-Time-Scale Simulation of Complex Reactive Mixtures: How Do Polyoxometalates Form?

Understanding the dynamics of reactive mixtures still challenges both experiments and theory. A relevant example can be found in the chemistry of molecular metal-oxide nanoclusters, also known as polyoxometalates. The high number of species potentially involved, the interconnectivity of the reaction network, and the precise control of the pH and concentrations needed in the synthesis of such species make the theoretical/computational treatment of such processes cumbersome. This work addresses this issue relying on a unique combination of recently developed computational methods that tackle the construction, kinetic simulation, and analysis of complex chemical reaction networks. By using the Bell-Evans-Polanyi approximation for estimating activation energies, and an accurate and robust linear scaling for correcting the computed pKa values, we report herein multi-time-scale kinetic simulations for the self-assembly processes of polyoxotungstates that comprise 22 orders of magnitude, from tens of femtoseconds to months of reaction time. This very large time span was required to reproduce very fast processes such as the acid/base equilibria (at 10-12 s), relatively slow reactions such as the formation of key clusters such as the metatungstate (at 103 s), and the very slow assembly of the decatungstate (at 106 s). Analysis of the kinetic data and of the reaction network topology shed light onto the details of the main reaction mechanisms, which explains the origin of kinetic and thermodynamic control followed by the reaction. Simulations at alkaline pH fully reproduce experimental evidence since clusters do not form under those conditions.

Algorithmic criteria for the validity of quasi-steady state and partial equilibrium models: the Michaelis-Menten reaction mechanism

We present "on the fly" algorithmic criteria for the accuracy and stability (non-stiffness) of reduced models constructed with the quasi-steady state and partial equilibrium approximations. The criteria comprise those introduced in Goussis (Combust Theor Model 16:869-926, 2012) that addressed the case where each fast time scale is due to one reaction and a new one that addresses the case where a fast time scale is due to more than one reactions. The development of these criteria is based on the ability to approximate accurately the fast and slow subspaces of the tangent space. Their validity is assessed on the basis of the Michaelis-Menten reaction mechanism, for which extensive literature is available regarding the validity of the existing various reduced models. The criteria predict correctly the regions in both the parameter and phase spaces where each of these models is valid. The findings are supported by numerical computations at indicative points in the parameter space. Due to their algorithmic character, these criteria can be readily employed for the reduction of large and complex mathematical models.

On rapid oscillations driving biological processes at disparate timescales

We consider a generic biological process described by a dynamical system, subject to an input signal with a high-frequency periodic component. The rapid oscillations of the input signal induce inherently multiscale dynamics, motivating order-reduction techniques. It is intuitive that the system behaviour is well approximated by its response to the averaged input signal. However, changes to the high-frequency component that preserve the average signal are beyond the reach of such intuitive reasoning. In this study, we explore system response under the influence of such an input signal by exploiting the timescale separation between high-frequency input variations and system response time. Employing the asymptotic method of multiple scales, we establish that, in some circumstances, the intuitive approach is simply the leading-order asymptotic contribution. We focus on higher-order corrections that capture the response to the details of the high-frequency component beyond its average. This approach achieves a reduction in system complexity while providing valuable insight into the structure of the response to the oscillations. We develop the general theory for nonlinear systems, while highlighting the important case of systems affine in the state and the input signal, presenting examples of both discrete and continuum state spaces. Importantly, this class of systems encompasses biochemical reaction networks described by the chemical master equation and its continuum approximations. Finally, we apply the framework to a nonlinear system describing mRNA translation and protein expression previously studied in the literature. The analysis shines new light on several aspects of the system quantification and both extends and simplifies results previously obtained.

A new hierarchical approach to multi-level model abstraction for simplifying ODE models of biological networks and a case study: The G1/S Checkpoint/DNA damage signalling pathways of mammalian cell cycle

Model reduction is an important topic in studies of biological systems. By reducing the complexity of large models through multi-level models while keeping the essence (biological meaning) of the model, model reduction can help answer many important questions about these systems. In this paper, we present a new reduction method based on hierarchical representation and a lumping approach. We used G1/S checkpoint pathway represented in Ordinary Differential Equations (ODE) in Iwamoto et al. (2011) as a case study to present this reduction method. The approach consists of two parts; the first part represents the biological network as a hierarchy (multiple levels) based on protein binding relations, which allowed us to model the biological network at different levels of abstraction. The second part applies different levels (level 1, 2 and 3) of lumping the species together depending on the level of the hierarchy, resulting in a reduced and transformed model for each level. The model at each level is a representation of the whole system and can address questions pertinent to that level. We develop and simulate reduced models for levels-1, 2 and 3 of lumping for the G1/S checkpoint pathway and evaluate the biological plausibility of the proposed method by comparing the results with the original ODE model of Iwamoto et al. (2011). The results for continuous dynamics of the G1/S checkpoint pathway with or without DNA-damage for reduced models of level- 1, 2 and 3 of lumping are in very good agreement and consistent with the original model results and with biological findings. Therefore, the reduced models (level-1, 2 and 3) can be used to study cell cycle progression in G1 and the effects of DNA damage on it. It is suitable for reducing complex ODE biological network models while retaining accurate continuous dynamics of the system. The 3 levels of the reduced models respectively achieved 20%, 26% and 31% reduction of the base model. Moreover, the reduced model is more efficient to run (30%, 44% and 52% time reduction for the three levels) and generate solutions than the original ODE model. Simplification of complex mathematical models is possible and the proposed reduction method has the potential to make an impact across many fields of biomedical research.

Bridging Mechanistic and Phenomenological Models of Complex Biological Systems

The inherent complexity of biological systems gives rise to complicated mechanistic models with a large number of parameters. On the other hand, the collective behavior of these systems can often be characterized by a relatively small number of phenomenological parameters. We use the Manifold Boundary Approximation Method (MBAM) as a tool for deriving simple phenomenological models from complicated mechanistic models. The resulting models are not black boxes, but remain expressed in terms of the microscopic parameters. In this way, we explicitly connect the macroscopic and microscopic descriptions, characterize the equivalence class of distinct systems exhibiting the same range of collective behavior, and identify the combinations of components that function as tunable control knobs for the behavior. We demonstrate the procedure for adaptation behavior exhibited by the EGFR pathway. From a 48 parameter mechanistic model, the system can be effectively described by a single adaptation parameter τ characterizing the ratio of time scales for the initial response and recovery time of the system which can in turn be expressed as a combination of microscopic reaction rates, Michaelis-Menten constants, and biochemical concentrations. The situation is not unlike modeling in physics in which microscopically complex processes can often be renormalized into simple phenomenological models with only a few effective parameters. The proposed method additionally provides a mechanistic explanation for non-universal features of the behavior.

A multi-time-scale analysis of chemical reaction networks: II. Stochastic systems

We consider stochastic descriptions of chemical reaction networks in which there are both fast and slow reactions, and for which the time scales are widely separated. We develop a computational algorithm that produces the generator of the full chemical master equation for arbitrary systems, and show how to obtain a reduced equation that governs the evolution on the slow time scale. This is done by applying a state space decomposition to the full equation that leads to the reduced dynamics in terms of certain projections and the invariant distributions of the fast system. The rates or propensities of the reduced system are shown to be the rates of the slow reactions conditioned on the expectations of fast steps. We also show that the generator of the reduced system is a Markov generator, and we present an efficient stochastic simulation algorithm for the slow time scale dynamics. We illustrate the numerical accuracy of the approximation by simulating several examples. Graph-theoretic techniques are used throughout to describe the structure of the reaction network and the state-space transitions accessible under the dynamics.

The relationship between stochastic and deterministic quasi-steady state approximations

Background

The quasi steady-state approximation (QSSA) is frequently used to reduce deterministic models of biochemical networks. The resulting equations provide a simplified description of the network in terms of non-elementary reaction functions (e.g. Hill functions). Such deterministic reductions are frequently a basis for heuristic stochastic models in which non-elementary reaction functions are used to define reaction propensities. Despite their popularity, it remains unclear when such stochastic reductions are valid. It is frequently assumed that the stochastic reduction can be trusted whenever its deterministic counterpart is accurate. However, a number of recent examples show that this is not necessarily the case.

Results

Here we explain the origin of these discrepancies, and demonstrate a clear relationship between the accuracy of the deterministic and the stochastic QSSA for examples widely used in biological systems. With an analysis of a two-state promoter model, and numerical simulations for a variety of other models, we find that the stochastic QSSA is accurate whenever its deterministic counterpart provides an accurate approximation over a range of initial conditions which cover the likely fluctuations from the quasi steady-state (QSS). We conjecture that this relationship provides a simple and computationally inexpensive way to test the accuracy of reduced stochastic models using deterministic simulations.

Conclusions

The stochastic QSSA is one of the most popular multi-scale stochastic simulation methods. While the use of QSSA, and the resulting non-elementary functions has been justified in the deterministic case, it is not clear when their stochastic counterparts are accurate. In this study, we show how the accuracy of the stochastic QSSA can be tested using their deterministic counterparts providing a concrete method to test when non-elementary rate functions can be used in stochastic simulations.

Electronic supplementary material

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

Designing Conservation Relations in Layered Synthetic Biomolecular Networks

In Synthetic Biology, biomolecular networks are designed and constructed to perform specified tasks. Design strategies for these networks tend to center on tuning the parameters of mathematical models to achieve a specified behavior, and implementing these parameters experimentally. This design strategy often assumes a fixed network structure that defines the possible behaviors, which may be too restrictive for our purposes. This paper investigates the extent to which the state space of a synthetic network can also be designed and shaped by parametric tuning. We exploit timescale separation to implement new, nonlinear, tunable conservation relations that hold for all times beyond a fast transient. We demonstrate an application of this design strategy by flexibly constraining the possible behaviors of a gene regulatory network through the design of fast protein interactions.

Stochastic representations of ion channel kinetics and exact stochastic simulation of neuronal dynamics

In this paper we provide two representations for stochastic ion channel kinetics, and compare the performance of exact simulation with a commonly used numerical approximation strategy. The first representation we present is a random time change representation, popularized by Thomas Kurtz, with the second being analogous to a "Gillespie" representation. Exact stochastic algorithms are provided for the different representations, which are preferable to either (a) fixed time step or (b) piecewise constant propensity algorithms, which still appear in the literature. As examples, we provide versions of the exact algorithms for the Morris-Lecar conductance based model, and detail the error induced, both in a weak and a strong sense, by the use of approximate algorithms on this model. We include ready-to-use implementations of the random time change algorithm in both XPP and Matlab. Finally, through the consideration of parametric sensitivity analysis, we show how the representations presented here are useful in the development of further computational methods. The general representations and simulation strategies provided here are known in other parts of the sciences, but less so in the present setting.

Model reduction by manifold boundaries

Understanding the collective behavior of complex systems from their basic components is a difficult yet fundamental problem in science. Existing model reduction techniques are either applicable under limited circumstances or produce "black boxes" disconnected from the microscopic physics. We propose a new approach by translating the model reduction problem for an arbitrary statistical model into a geometric problem of constructing a low-dimensional, submanifold approximation to a high-dimensional manifold. When models are overly complex, we use the observation that the model manifold is bounded with a hierarchy of widths and propose using the boundaries as submanifold approximations. We refer to this approach as the manifold boundary approximation method. We apply this method to several models, including a sum of exponentials, a dynamical systems model of protein signaling, and a generalized Ising model. By focusing on parameters rather than physical degrees of freedom, the approach unifies many other model reduction techniques, such as singular limits, equilibrium approximations, and the renormalization group, while expanding the domain of tractable models. The method produces a series of approximations that decrease the complexity of the model and reveal how microscopic parameters are systematically "compressed" into a few macroscopic degrees of freedom, effectively building a bridge between the microscopic and the macroscopic descriptions.

Measuring edge importance: a quantitative analysis of the stochastic shielding approximation for random processes on graphs

Mathematical models of cellular physiological mechanisms often involve random walks on graphs representing transitions within networks of functional states. Schmandt and Galán recently introduced a novel stochastic shielding approximation as a fast, accurate method for generating approximate sample paths from a finite state Markov process in which only a subset of states are observable. For example, in ion-channel models, such as the Hodgkin-Huxley or other conductance-based neural models, a nerve cell has a population of ion channels whose states comprise the nodes of a graph, only some of which allow a transmembrane current to pass. The stochastic shielding approximation consists of neglecting fluctuations in the dynamics associated with edges in the graph not directly affecting the observable states. We consider the problem of finding the optimal complexity reducing mapping from a stochastic process on a graph to an approximate process on a smaller sample space, as determined by the choice of a particular linear measurement functional on the graph. The partitioning of ion-channel states into conducting versus nonconducting states provides a case in point. In addition to establishing that Schmandt and Galán's approximation is in fact optimal in a specific sense, we use recent results from random matrix theory to provide heuristic error estimates for the accuracy of the stochastic shielding approximation for an ensemble of random graphs. Moreover, we provide a novel quantitative measure of the contribution of individual transitions within the reaction graph to the accuracy of the approximate process.

Time-scale separation--Michaelis and Menten's old idea, still bearing fruit

Michaelis and Menten introduced to biochemistry the idea of time-scale separation, in which part of a system is assumed to be operating sufficiently fast compared to the rest so that it may be taken to have reached a steady state. This allows, in principle, the fast components to be eliminated, resulting in a simplified description of the system's behaviour. Similar ideas have been widely used in different areas of biology, including enzyme kinetics, protein allostery, receptor pharmacology, gene regulation and post-translational modification. However, the methods used have been independent and ad hoc. In the present study, we review the use of time-scale separation as a means to simplify the description of molecular complexity and discuss recent work setting out a single framework that unifies these separate calculations. The framework offers new capabilities for mathematical analysis and helps to do justice to Michaelis and Menten's insights about individual enzymes in the context of multi-enzyme biological systems.

Rigorous elimination of fast stochastic variables from the linear noise approximation using projection operators

The linear noise approximation (LNA) offers a simple means by which one can study intrinsic noise in monostable biochemical networks. Using simple physical arguments, we have recently introduced the slow-scale LNA (ssLNA), which is a reduced version of the LNA under conditions of timescale separation. In this paper we present the first rigorous derivation of the ssLNA using the projection operator technique and show that the ssLNA follows uniquely from the standard LNA under the same conditions of timescale separation as those required for the deterministic quasi-steady-state approximation. We also show that the large molecule number limit of several common stochastic model reduction techniques under timescale separation conditions constitutes a special case of the ssLNA.

The effect of the signalling scheme on the robustness of pattern formation in development

Pattern formation in development is a complex process which involves spatially distributed signals called morphogens that influence gene expression and thus the phenotypic identity of cells. Usually different cell types are spatially segregated, and the boundary between them may be determined by a threshold value of some state variable. The question arises as to how sensitive the location of such a boundary is to variations in properties, such as parameter values, that characterize the system. Here, we analyse both deterministic and stochastic reaction-diffusion models of pattern formation with a view towards understanding how the signalling scheme used for patterning affects the variability of boundary determination between cell types in a developing tissue.

A "trimer of dimers"-based model for the chemotactic signal transduction network in bacterial chemotaxis

The network that controls chemotaxis in Escherichia coli is one of the most completely characterized signal transduction systems to date. Receptor clustering accounts for characteristics such as high sensitivity, precise adaptation over a wide dynamic range of ligand concentrations, and robustness to variations in the amounts of intracellular proteins. To gain insights into the structure-function relationship of receptor clusters and understand the mechanism behind the high-performance signaling, we develop and analyze a model for a single trimer of dimers. This new model extends an earlier model (Spiro et al. in Proc. Natl. Acad. Sci. 94:7263-7268, 1997) to incorporate the recent experimental findings that the core structure of receptor clusters is the trimer of receptor dimers. We show that the model can reproduce most of the experimentally-observed behaviors, including excitation, adaptation, high sensitivity, and robustness to parameter variations. In addition, the model makes a number of new predictions as to how the adaptation time varies with the expression level of various proteins involved in signal transduction. Our results provide a more mechanistically-based description of the structure-function relationship for the signaling system, and show the key role of the interaction among dimer members of the trimer in the chemotactic response of cells.

A linear framework for time-scale separation in nonlinear biochemical systems

Cellular physiology is implemented by formidably complex biochemical systems with highly nonlinear dynamics, presenting a challenge for both experiment and theory. Time-scale separation has been one of the few theoretical methods for distilling general principles from such complexity. It has provided essential insights in areas such as enzyme kinetics, allosteric enzymes, G-protein coupled receptors, ion channels, gene regulation and post-translational modification. In each case, internal molecular complexity has been eliminated, leading to rational algebraic expressions among the remaining components. This has yielded familiar formulas such as those of Michaelis-Menten in enzyme kinetics, Monod-Wyman-Changeux in allostery and Ackers-Johnson-Shea in gene regulation. Here we show that these calculations are all instances of a single graph-theoretic framework. Despite the biochemical nonlinearity to which it is applied, this framework is entirely linear, yet requires no approximation. We show that elimination of internal complexity is feasible when the relevant graph is strongly connected. The framework provides a new methodology with the potential to subdue combinatorial explosion at the molecular level.

Guaranteed error bounds for structured complexity reduction of biochemical networks

Biological systems are typically modelled by nonlinear differential equations. In an effort to produce high fidelity representations of the underlying phenomena, these models are usually of high dimension and involve multiple temporal and spatial scales. However, this complexity and associated stiffness makes numerical simulation difficult and mathematical analysis impossible. In order to understand the functionality of these systems, these models are usually approximated by lower dimensional descriptions. These can be analysed and simulated more easily, and the reduced description also simplifies the parameter space of the model. This model reduction inevitably introduces error: the accuracy of the conclusions one makes about the system, based on reduced models, depends heavily on the error introduced in the reduction process. In this paper we propose a method to calculate the error associated with a model reduction algorithm, using ideas from dynamical systems. We first define an error system, whose output is the error between observables of the original and reduced systems. We then use convex optimisation techniques in order to find approximations to the error as a function of the initial conditions. In particular, we use the Sum of Squares decomposition of polynomials in order to compute an upper bound on the worst-case error between the original and reduced systems. We give biological examples to illustrate the theory, which leads us to a discussion about how these techniques can be used to model-reduce large, structured models typical of systems biology.

A new method for choosing the computational cell in stochastic reaction-diffusion systems

How to choose the computational compartment or cell size for the stochastic simulation of a reaction-diffusion system is still an open problem, and a number of criteria have been suggested. A generalized measure of the noise for finite-dimensional systems based on the largest eigenvalue of the covariance matrix of the number of molecules of all species has been suggested as a measure of the overall fluctuations in a multivariate system, and we apply it here to a discretized reaction-diffusion system. We show that for a broad class of first-order reaction networks this measure converges to the square root of the reciprocal of the smallest mean species number in a compartment at the steady state. We show that a suitably re-normalized measure stabilizes as the volume of a cell approaches zero, which leads to a criterion for the maximum volume of the compartments in a computational grid. We then derive a new criterion based on the sensitivity of the entire network, not just of the fastest step, that predicts a grid size that assures that the concentrations of all species converge to a spatially-uniform solution. This criterion applies for all orders of reactions and for reaction rate functions derived from singular perturbation or other reduction methods, and encompasses both diffusing and non-diffusing species. We show that this predicts the maximal allowable volume found in a linear problem, and we illustrate our results with an example motivated by anterior-posterior pattern formation in Drosophila, and with several other examples.

Reduced models of networks of coupled enzymatic reactions

The Michaelis-Menten equation has played a central role in our understanding of biochemical processes. It has long been understood how this equation approximates the dynamics of irreversible enzymatic reactions. However, a similar approximation in the case of networks, where the product of one reaction can act as an enzyme in another, has not been fully developed. Here we rigorously derive such an approximation in a class of coupled enzymatic networks where the individual interactions are of Michaelis-Menten type. We show that the sufficient conditions for the validity of the total quasi-steady state assumption (tQSSA), obtained in a single protein case by Borghans, de Boer and Segel can be extended to sufficient conditions for the validity of the tQSSA in a large class of enzymatic networks. Secondly, we derive reduced equations that approximate the network's dynamics and involve only protein concentrations. This significantly reduces the number of equations necessary to model such systems. We prove the validity of this approximation using geometric singular perturbation theory and results about matrix differentiation. The ideas used in deriving the approximating equations are quite general, and can be used to systematize other model reductions.