Bistability and oscillations in chemical reaction networks

Bifurcation theory is one of the most widely used approaches for analysis of dynamical behaviour of chemical and biochemical reaction networks. Some of the interesting qualitative behaviour that are analyzed are oscillations and bistability (a situation where a system has at least two coexisting stable equilibria). Both phenomena have been identified as central features of many biological and biochemical systems. This paper, using the theory of stoichiometric network analysis (SNA) and notions from algebraic geometry, presents sufficient conditions for a reaction network to display bifurcations associated with these phenomena. The advantage of these conditions is that they impose fewer algebraic conditions on model parameters than conditions associated with standard bifurcation theorems. To derive the new conditions, a coordinate transformation will be made that will guarantee the existence of branches of positive equilibria in the system. This is particularly useful in mathematical biology, where only positive variable values are considered to be meaningful. The first part of the paper will be an extended introduction to SNA and algebraic geometry-related methods which are used in the coordinate transformation and set up of the theorems. In the second part of the paper we will focus on the derivation of bifurcation conditions using SNA and algebraic geometry. Conditions will be derived for three bifurcations: the saddle-node bifurcation, a simple branching point, both linked to bistability, and a simple Hopf bifurcation. The latter is linked to oscillatory behaviour. The conditions derived are sufficient and they extend earlier results from stoichiometric network analysis as can be found in (Aguda and Clarke in J Chem Phys 87:3461–3470, 1987; Clarke and Jiang in J Chem Phys 99:4464–4476, 1993; Gatermann et al. in J Symb Comput 40:1361–1382, 2005). In these papers some necessary conditions for two of these bifurcations were given. A set of examples will illustrate that algebraic conditions arising from given sufficient bifurcation conditions are not more difficult to interpret nor harder to calculate than those arising from necessary bifurcation conditions. Hence an increasing amount of information is gained at no extra computational cost. The theory can also be used in a second step for a systematic bifurcation analysis of larger reaction networks. We have added a dedication of the paper to K. Gatermann.     

Quasi-steady-state approximation for chemical reaction networks

 The parameter embedding leading to the quasi-steady-state approximation of Heinrich [9] is investigated within the theory of invariant manifolds of Fenichel [4] in order to clarify the essential assumptions needed for this reduction to a low dimensional system. In particular, the concept of pool-variables can be avoided in this generalized approach. Moreover, the dominating influence of the slow subnetwork over the complementary fast subnetwork is interpreted geometrically and in chemical terms and this can be seen as an “enslaving” of the fast subsystem by the slow subsystem. Finally, the results are applied to a system of slime mould communication [6, 7, 13] and to a maltose transport system [2, 3]. Received: 16 June 1997     

The flow topology of chemical reaction networks

Analysis of the flow topology of dynamic systems modelling chemical reaction networks leads to the following conclusions: (1) a connected chemical system with a negative feedback circuit containing all except one of the internal species (phase determining intermediates) can exhibit chaotic oscillations; (2) a chemical system in which the only negative feedback circuit contains all of the internal species can exhibit periodic oscillations; (3) a chemical system containing no negative feedback circuits cannot exhibit oscillatory behavior but only multiple steady states. In this analysis a negative feedback circuit refers to a circuit containing an odd number of inhibitory relationships between internal species.     

Graph-theoretic methods for the analysis of chemical and biochemical networks. I. Multistability and oscillations in ordinary differential equation models

A chemical mechanism is a model of a chemical reaction network consisting of a set of elementary reactions that express how molecules react with each other. In classical mass-action kinetics, a mechanism implies a set of ordinary differential equations (ODEs) which govern the time evolution of the concentrations. In this article, ODE models of chemical kinetics that have the potential for multiple positive equilibria or oscillations are studied. We begin by considering some methods of stability analysis based on the digraph of the Jacobian matrix. We then prove two theorems originally given by A. N. Ivanova which correlate the bifurcation structure of a mass-action model to the properties of a bipartite graph with nodes representing chemical species and reactions. We provide several examples of the application of these theorems.     

Concordant chemical reaction networks and the Species-Reaction Graph

In a recent paper it was shown that, for chemical reaction networks possessing a subtle structural property called concordance, dynamical behavior of a very circumscribed (and largely stable) kind is enforced, so long as the kinetics lies within the very broad and natural weakly monotonic class. In particular, multiple equilibria are precluded, as are degenerate positive equilibria. Moreover, under certain circumstances, also related to concordance, all real eigenvalues associated with a positive equilibrium are negative. Although concordance of a reaction network can be decided by readily available computational means, we show here that, when a nondegenerate network’s Species-Reaction Graph satisfies certain mild conditions, concordance and its dynamical consequences are ensured. These conditions are weaker than earlier ones invoked to establish kinetic system injectivity, which, in turn, is just one ramification of network concordance. Because the Species-Reaction Graph resembles pathway depictions often drawn by biochemists, results here expand the possibility of inferring significant dynamical information directly from standard biochemical reaction diagrams.     

Atomic motion in enzymatic reaction coordinates

Atomic excursions of reactants in enzymatic catalytic sites can be estimated from high-resolution crystal structures of enzyme complexes with substrates, transition state analog inhibitors and products. Transition state structures, defined from kinetic isotope effect studies, are compared to crystallographic structures to validate the properties of the transition state analog. Atomic excursions in enzymatic catalytic sites can differ from those in solution and define the role of the enzymatic catalyst in directing atomic motion.     

A Petri net approach to the study of persistence in chemical reaction networks

Persistence is the property, for differential equations in R(n), that solutions starting in the positive orthant do not approach the boundary of the orthant. For chemical reactions and population models, this translates into the non-extinction property: provided that every species is present at the start of the reaction, no species will tend to be eliminated in the course of the reaction. This paper provides checkable conditions for persistence of chemical species in reaction networks, using concepts and tools from Petri net theory, and verifies these conditions on various systems which arise in the modeling of cell signaling pathways.     

Computing weakly reversible linearly conjugate chemical reaction networks with minimal deficiency

Mass-action kinetics is frequently used in systems biology to model the behavior of interacting chemical species. Many important dynamical properties are known to hold for such systems if their underlying networks are weakly reversible and have a low deficiency. In particular, the Deficiency Zero and Deficiency One Theorems guarantee strong regularity with regards to the number and stability of positive equilibrium states. It is also known that chemical reaction networks with distinct reaction structure can admit mass-action systems with the same qualitative dynamics. The theory of linear conjugacy encapsulates the cases where this relationship is captured by a linear transformation. In this paper, we propose a mixed-integer linear programming algorithm capable of determining the minimal deficiency weakly reversible reaction network which admits a mass-action system which is linearly conjugate to a given reaction network.     

Multiple equilibria in complex chemical reaction networks: extensions to entrapped species models

In two earlier papers, means were provided to decide the capacity of complex chemical reaction networks, taken with mass-action kinetics, to admit multiple equilibria in the context of the isothermal homogeneous continuous flow stirred tank reactor (CFSTR). In such a reactor, all species are deemed to be in the outflow, a fact which has an important bearing on the nature of the governing equations. In contrast, one can imagine CFSTR-like models of the cell in which certain large molecules (e.g., enzymes) remain entrapped within the cell, whereas smaller ones (e.g., metabolites) are free to diffuse through the cell boundary. Although such models bear a strong physical resemblance to the classical CFSTR picture, there are substantive differences in the corresponding mathematics. Without a presumption of mass-action kinetics, this research is intended to indicate a general way in which results about uniqueness of equilibria in the classical CFSTR context extend to entrapped species models.     

Recent developments in methods for identifying reaction coordinates

In the study of rare events in complex systems with many degrees of freedom, a key element is to identify the reaction coordinates of a given process. Over recent years, a number of methods and protocols have been developed to extract the reaction coordinates based on limited information from molecular dynamics simulations. In this review, we provide a brief survey over a number of major methods developed in the past decade, some of which are discussed in greater detail, to provide an overview of the problems that are partially solved and challenges that still remain. A particular emphasis has been placed on methods for identifying reaction coordinates that are related to the committor.     

Validation of archived chemical shifts through atomic coordinates

The public archives containing protein information in the form of NMR chemical shift data at the BioMagResBank (BMRB) and of 3D structure coordinates at the Protein Data Bank are continuously expanding. The quality of the data contained in these archives, however, varies. The main issue for chemical shift values is that they are determined relative to a reference frequency. When this reference frequency is set incorrectly, all related chemical shift values are systematically offset. Such wrongly referenced chemical shift values, as well as other problems such as chemical shift values that are assigned to the wrong atom, are not easily distinguished from correct values and effectively reduce the usefulness of the archive. We describe a new method to correct and validate protein chemical shift values in relation to their 3D structure coordinates. This method classifies atoms using two parameters: the per‐atom solvent accessible surface area (as calculated from the coordinates) and the secondary structure of the parent amino acid. Through the use of Gaussian statistics based on a large database of 3220 BMRB entries, we obtain per‐entry chemical shift corrections as well as Z scores for the individual chemical shift values. In addition, information on the error of the correction value itself is available, and the method can retain only dependable correction values. We provide an online resource with chemical shift, atom exposure, and secondary structure information for all relevant BMRB entries ( http://www.ebi.ac.uk/pdbe/nmr/vasco ) and hope this data will aid the development of new chemical shift‐based methods in NMR. Proteins 2010. © 2010 Wiley‐Liss, Inc.     

Distance-type reaction coordinates for modelling activated processes

A suitably defined distance is the simplest parameter for measuring the difference between two positions, orientations, and/or conformations of a molecular system. Distances also provide a first guess for the reaction coordinates of activated processes. It is shown here that mass-weighted distances possess remarkable mechanical and statistical mechanical properties. They allow us to restrict motions to internal coordinates of a molecule in a simple way where this demand makes sense. Moreover, the computation of free energy changes and rates is facilitated by simple explicit formulae. The numerical treatment of a rate process in a peptide, the ring flip of a phenylalanine, demonstrates the practical application of our results. It also indicates the role of internal friction in macromolecules and the need to consider transmission coefficients.     

Mutation in autocatalytic reaction networks

A class of kinetic equations describing catalysed and template induced replication, and mutation is introduced. This ODE in its most general form is split into two vector fields, a replication and a mutation field. The mutation field is considered as a perturbation of the replicator equation. The perturbation expansion is a Taylor series in a mutation parameter . First, second and higher order contributions are computed by means of the conventional Rayleigh-Schrödinger approach. Qualitative shifts in the positions of rest points and limit cycles on the boundary of the physically meaningful part of concentration space are predicted from flow topologies. The results of the topological analysis are summarized in two theorems which turned out to be useful in applications: the rest point migration theorem (RPM) and the limit cycle migration theorem (LCM). Quantitative expressions for the shifts of rest points are computed directly from the perturbation expansion. The concept is applied to a collection of selected examples from biophysical chemistry and biology.     

Analyzing molecular reaction networks

Pathways are typically the central concept in the analysis of biochemical reaction networks. A pathway can be interpreted as a chain of enzymatical reactions performing a specific biological function. A common way to study metabolic networks are minimal pathways that can operate at steady state called elementary modes. The theory of chemical organizations has recently been used to decompose biochemical networks into algebraically closed and self-maintaining subnetworks termed organizations. The aim of this paper is to elucidate the relation between these two concepts. Whereas elementary modes represent the boundaries of the potential behavior of the network, organizations define metabolite compositions that are likely to be present in biological feasible situations. Hence, steady state organizations consist of combinations of elementary modes. On the other hand, it is possible to assign a unique (and possibly empty) set of organizations to each elementary mode, indicating the metabolites accompanying the active pathway in a feasible steady state.     

Noise-induced modulation of the relaxation kinetics around a non-equilibrium steady state of non-linear chemical reaction networks

Stochastic effects from correlated noise non-trivially modulate the kinetics of non-linear chemical reaction networks. This is especially important in systems where reactions are confined to small volumes and reactants are delivered in bursts. We characterise how the two noise sources confinement and burst modulate the relaxation kinetics of a non-linear reaction network around a non-equilibrium steady state. We find that the lifetimes of species change with burst input and confinement. Confinement increases the lifetimes of all species that are involved in any non-linear reaction as a reactant. Burst monotonically increases or decreases lifetimes. Competition between burst-induced and confinement-induced modulation may hence lead to a non-monotonic modulation. We quantify lifetime as the integral of the time autocorrelation function (ACF) of concentration fluctuations around a non-equilibrium steady state of the reaction network. Furthermore, we look at the first and second derivatives of the ACF, each of which is affected in opposite ways by burst and confinement. This allows discriminating between these two noise sources. We analytically derive the ACF from the linear Fokker-Planck approximation of the chemical master equation in order to establish a baseline for the burst-induced modulation at low confinement. Effects of higher confinement are then studied using a partial-propensity stochastic simulation algorithm. The results presented here may help understand the mechanisms that deviate stochastic kinetics from its deterministic counterpart. In addition, they may be instrumental when using fluorescence-lifetime imaging microscopy (FLIM) or fluorescence-correlation spectroscopy (FCS) to measure confinement and burst in systems with known reaction rates, or, alternatively, to correct for the effects of confinement and burst when experimentally measuring reaction rates.