Graph-based, dynamics-preserving reduction of (bio)chemical systems.

Complex dynamical systems are often governed by equations containing many unknown parameters whose precise values may or may not be important for the system's dynamics. In particular, for chemical and biochemical systems, there may be some reactions or subsystems that are inessential to understanding the bifurcation structure and consequent behavior of a model, such as oscillations, multistationarity and patterning. Due to the size, complexity and parametric uncertainties of many (bio)chemical models, a dynamics-preserving reduction scheme that is able to isolate the necessary contributors to particular dynamical behaviors would be useful. In this contribution, we describe model reduction methods for mass-action (bio)chemical models based on the preservation of instability-generating subnetworks known as critical fragments. These methods focus on structural conditions for instabilities and so are parameter-independent. We apply these results to an existing model for the control of the synthesis of the NO-detoxifying enzyme Hmp in Escherichia coli that displays bistability.

Parametric Sensitivity Analysis of Oscillatory Delay Systems with an Application to Gene Regulation.

A parametric sensitivity analysis for periodic solutions of delay-differential equations is developed. Because phase shifts cause the sensitivity coefficients of a periodic orbit to diverge, we focus on sensitivities of the extrema, from which amplitude sensitivities are computed, and of the period. Delay-differential equations are often used to model gene expression networks. In these models, the parametric sensitivities of a particular genotype define the local geometry of the evolutionary landscape. Thus, sensitivities can be used to investigate directions of gradual evolutionary change. An oscillatory protein synthesis model whose properties are modulated by RNA interference is used as an example. This model consists of a set of coupled delay-differential equations involving three delays. Sensitivity analyses are carried out at several operating points. Comments on the evolutionary implications of the results are offered.

Analytic delay distributions for a family of gene transcription models.

Models intended to describe the time evolution of a gene network must somehow include transcription, the DNA-templated synthesis of RNA, and translation, the RNA-templated synthesis of proteins. In eukaryotes, the DNA template for transcription can be very long, often consisting of tens of thousands of nucleotides, and lengthy pauses may punctuate this process. Accordingly, transcription can last for many minutes, in some cases hours. There is a long history of introducing delays in gene expression models to take the transcription and translation times into account. Here we study a family of detailed transcription models that includes initiation, elongation, and termination reactions. We establish a framework for computing the distribution of transcription times, and work out these distributions for some typical cases. For elongation, a fixed delay is a good model provided elongation is fast compared to initiation and termination, and there are no sites where long pauses occur. The initiation and termination phases of the model then generate a nontrivial delay distribution, and elongation shifts this distribution by an amount corresponding to the elongation delay. When initiation and termination are relatively fast, the distribution of elongation times can be approximated by a Gaussian. A convolution of this Gaussian with the initiation and termination time distributions gives another analytic approximation to the transcription time distribution. If there are long pauses during elongation, because of the modularity of the family of models considered, the elongation phase can be partitioned into reactions generating a simple delay (elongation through regions where there are no long pauses), and reactions whose distribution of waiting times must be considered explicitly (initiation, termination, and motion through regions where long pauses are likely). In these cases, the distribution of transcription times again involves a nontrivial part and a shift due to fast elongation processes.

Probabilistic models of uORF-mediated ATF4 translation control.

ATF4 is a key transcription factor that activates transcription of genes needed to respond to cellular stress. Although the mRNA encoding ATF4 is present at constant levels in the cell during the initial response, translation of ATF4 increases under conditions of cellular stress while the global translation rate decreases. We study two models for the control system that regulates the translation of ATF4, both based on the Vattem-Wek hypothesis. This hypothesis is based on a race to reload, following the translation of a small upstream open reading frame (uORF), the ternary complex that brings the initiator tRNA to the ribosome as the 40S subunit scans along the mRNA, encountering first a start codon for an inhibitory uORF whose reading frame overlaps the start of the ATF4 coding sequence. We develop a pair of simple, analytic, probabilistic models, one of which assumes all nucleotide triplets have identical kinetic properties, while the other recognizes the existence of triplets at which the ternary complex loads more efficiently. We also consider two different functions representing the dependence of the rate of initiation at uORF1 on the ternary complex concentration. In keeping with the theme of this Special Issue, we studied the properties of these models in a Maple document, which can easily be modified to consider different parameters, translation rate initiation functions, and so on.

On the quasi-steady-state approximation in an open Michaelis-Menten reaction mechanism.

The conditions for the validity of the standard quasi-steady-state approximation in the Michaelis-Menten mechanism in a closed reaction vessel have been well studied, but much less so the conditions for the validity of this approximation for the system with substrate inflow. We analyze quasi-steady-state scenarios for the open system attributable to singular perturbations, as well as less restrictive conditions. For both settings we obtain distinguished invariant manifolds and time scale estimates, and we highlight the special role of singular perturbation parameters in higher order approximations of slow manifolds. We close the paper with a discussion of distinguished invariant manifolds in the global phase portrait.

Small binding-site clearance delays are not negligible in gene expression modeling.

During the templated biopolymerization processes of transcription and translation, a macromolecular machine, either an RNA polymerase or a ribosome, binds to a specific site on the template. Due to the sizes of these enzymes, there is a waiting time before one clears the binding site and another can bind. These clearance delays are relatively short, and one might think that they could be neglected. However, in the case of transcription, these clearance delays are associated with conservation laws, resulting in surprisingly large effects on the bifurcation diagrams in models of gene expression networks. We study an example of this phenomenon in a model of a gene regulated by a non-coding RNA displaying bistability. Neglecting the binding-site clearance delays in this model can only be compensated for by making ad hoc, unphysical adjustments to the model's kinetic constants.

Heineken, Tsuchiya and Aris on the mathematical status of the pseudo-steady state hypothesis: A classic from volume 1 of Mathematical Biosciences.

Volume 1, Issue 1 of Mathematical Biosciences was the venue for a now-classic paper on the application of singular perturbation theory in enzyme kinetics, "On the mathematical status of the pseudo-steady state hypothesis of biochemical kinetics" by F. G. Heineken, H. M. Tsuchiya and R. Aris. More than 50 years have passed, and yet this paper continues to be studied and mined for insights. This perspective discusses both the strengths and weaknesses of the work presented in this paper. For many, the justification of the pseudo-steady-state approximation using singular perturbation theory is the main achievement of this paper. However, there is so much more material here, which laid the foundation for a great deal of research in mathematical biochemistry in the intervening decades. The parameterization of the equations, construction of the first-order uniform singular-perturbation solution, and an attempt to apply similar principles to the pseudo-equilibrium approximation are discussed in particular detail.

A mathematical model of the biochemical network underlying left-right asymmetry establishment in mammals.

The expression of the TGF-ß protein Nodal on the left side of vertebrate embryos is a determining event in the development of internal-organ asymmetry. We present a mathematical model for the control of the expression of Nodal and its antagonist Lefty consisting entirely of realistic elementary reactions. We analyze the model in the absence of Lefty and find a wide range of parameters over which bistability (two stable steady states) is observed, with one stable steady state a low-Nodal state corresponding to the right-hand developmental fate, and the other a high-Nodal state corresponding to the left. We find that bistability requires a transcription factor containing two molecules of phosphorylated Smad2. A numerical survey of the full model, including Lefty, shows the effects of Lefty on the potential for bistability, and on the conditions that lead to the system reaching one or the other steady state.

Temporal metabolic partitioning of the yeast and protist cellular networks: the cell is a global scale-invariant (fractal or self-similar) multioscillator.

Britton Chance, electronics expert when a teenager, became an enthusiastic student of biological oscillations, passing on this enthusiasm to many students and colleagues, including one of us (DL). This historical essay traces BC's influence through the accumulated work of DL to DL's many collaborators. The overall temporal organization of mass-energy, information, and signaling networks in yeast in self-synchronized continuous cultures represents, until now, the most characterized example of in vivo elucidation of time structure. Continuous online monitoring of dissolved gases by direct measurement (membrane-inlet mass spectrometry, together with NAD(P)H and flavin fluorescence) gives strain-specific dynamic information from timescales of minutes to hours as does two-photon imaging. The predominantly oscillatory behavior of network components becomes evident, with spontaneously synchronized cellular respiration cycles between discrete periods of increased oxygen consumption (oxidative phase) and decreased oxygen consumption (reductive phase). This temperature-compensated ultradian clock provides coordination, linking temporally partitioned functions by direct feedback loops between the energetic and redox state of the cell and its growing ultrastructure. Multioscillatory outputs in dissolved gases with 13 h, 40 min, and 4 min periods gave statistical self-similarity in power spectral and relative dispersional analyses: i.e., complex nonlinear (chaotic) behavior and a functional scale-free (fractal) network operating simultaneously over several timescales.

Observation of a chaotic multioscillatory metabolic attractor by real-time monitoring of a yeast continuous culture.

We monitored a continuous culture of the yeast Saccharomyces cerevisiae by membrane-inlet mass spectrometry. This technique allows very rapid simultaneous measurements (one point every 12 s) of several dissolved gases. During our experiment, the culture exhibited a multioscillatory mode in which the dissolved oxygen and carbon dioxide records displayed periodicities of 13 h, 36 min and 4 min. The 36- and 4-min modes were not visible at all times, but returned at regular intervals during the 13-h cycle. The 4-min mode, which has not previously been described in continuous culture, can also be seen when the culture displays simpler oscillatory behavior. The data can be used to visualize a metabolic attractor of this system, i.e. the set of dissolved gas concentrations which are consistent with the multioscillatory state. Computation of the leading Lyapunov exponent reveals the dynamics on this attractor to be chaotic.

Oscillations of the internal CO(2) concentration in tobacco leaves transferred to low CO(2).

Measurement of the internal CO(2) concentration (Ci) in tobacco leaves using a fast-response CO(2) exchange system showed that in the light, switching from 350 microLL(-1) to a low CO(2) concentration of 36.5 microLL(-1) (promoting high photorespiration) resulted in the Ci oscillating near the value of CO(2) compensation point (Gamma*). The oscillations are highly irregular, the range of Ci varying by 2-4 microLL(-1) in substomatal cavities with a period of a few seconds. The statistical properties of the time series became stationary after a transient of approximately 100s following transfer to low CO(2). Attractor reconstruction shows that the observed oscillations are not chaotic but exhibit stochastic behavior. The period of oscillations is consistent with the duration of photorespiratory post-illumination burst (PIB). We suggest that the observed oscillations may be due to a similar mechanism to that which leads to PIB, and may play a role in switching mitochondrial operation between oxidation of the photorespiratory glycine and of the tricarboxylic acid cycle substrates.

Optimal observability of sustained stochastic competitive inhibition oscillations at organellar volumes.

When molecules are present in small numbers, such as is frequently the case in cells, the usual assumptions leading to differential rate equations are invalid and it is necessary to use a stochastic description which takes into account the randomness of reactive encounters in solution. We display a very simple biochemical model, ordinary competitive inhibition with substrate inflow, which is only capable of damped oscillations in the deterministic mass-action rate equation limit, but which displays sustained oscillations in stochastic simulations. We define an observability parameter, which is essentially just the ratio of the amplitude of the oscillations to the mean value of the concentration. A maximum in the observability is seen as the volume is varied, a phenomenon we name system-size observability resonance by analogy with other types of stochastic resonance. For the parameters of this study, the maximum in the observability occurs at volumes similar to those of bacterial cells or of eukaryotic organelles.

Validation of an algorithm for delay stochastic simulation of transcription and translation in prokaryotic gene expression.

The quantitative modeling of gene transcription and translation requires a treatment of two key features: stochastic fluctuations due to the limited copy numbers of key molecules (genes, RNA polymerases, ribosomes), and delayed output due to the time required for biopolymer synthesis. Recently proposed algorithms allow for efficient simulations of such systems. However, it is critical to know whether the results of delay stochastic simulations agree with those from more detailed models of the transcription and translation processes. We present a generalization of previous delay stochastic simulation algorithms which allows both for multiple delays and for distributions of delay times. We show that delay stochastic simulations closely approximate simulations of a detailed transcription model except when two-body effects (e.g. collisions between polymerases on a template strand) are important. Finally, we study a delay stochastic model of prokaryotic transcription and translation which reproduces observations from a recent experimental study in which a single gene was expressed under the control of a repressed lac promoter in E. coli cells. This demonstrates our ability to quantitatively model gene expression using these new methods.

The mechanism distinguishability problem in biochemical kinetics: the single-enzyme, single-substrate reaction as a case study.

A theoretical analysis of the distinguishability problem of two rival models of the single enzyme-single substrate reaction, the Michaelis-Menten and Henri mechanisms, is presented. We also outline a general approach for analysing the structural indistinguishability between two mechanisms. The approach involves constructing, if possible, a smooth mapping between the two candidate models. Evans et al. [N.D. Evans, M.J. Chappell, M.J. Chapman, K.R. Godfrey, Structural indistinguishability between uncontrolled (autonomous) nonlinear analytic systems, Automatica 40 (2004) 1947-1953] have shown that if, in addition, either of the mechanisms satisfies a particular criterion then such a transformation always exists when the models are indistinguishable from their experimentally observable outputs. The approach is applied to the single enzyme-single substrate reaction mechanism. In principle, mechanisms can be distinguished using this analysis, but we show that our ability to distinguish mechanistic models depends both on the precise measurements made, and on our knowledge of the system prior to performing the kinetics experiments.

Reaction-diffusion models of development with state-dependent chemical diffusion coefficients.

Reaction-diffusion models are widely used to model developmental processes. The great majority of current models invoke constant diffusion coefficients. However, the diffusion of metabolites or signals through tissues is frequently such that this assumption may reasonably be questioned. We consider several different physical mechanisms leading to effective diffusion coefficients in biological tissues which vary with the local conditions, including models in which juxtacrine signaling results in the diffusion of a signal in the absence of material transport. We develop a mathematical formalism for transforming local transport laws into diffusive terms. This procedure is appropriate when the typical length scale over which the concentrations change significantly is much greater than the dimensions of a cell. We review previous developmental models which considered the possibility of state-dependent diffusion coefficients. We also provide a few new motivating examples.

Invariant manifold methods for metabolic model reduction.

After the decay of transients, the behavior of a set of differential equations modeling a chemical or biochemical system generally rests on a low-dimensional surface which is an invariant manifold of the flow. If an equation for such a manifold can be obtained, the model has effectively been reduced to a smaller system of differential equations. Using perturbation methods, we show that the distinction between rapidly decaying and long-lived (slow) modes has a rigorous basis. We show how equations for attracting invariant (slow) manifolds can be constructed by a geometric approach based on functional equations derived directly from the differential equations. We apply these methods to two simple metabolic models. (c) 2001 American Institute of Physics.

Turing-Hopf instability in biochemical reaction networks arising from pairs of subnetworks.

Network conditions for Turing instability in biochemical systems with two biochemical species are well known and involve autocatalysis or self-activation. On the other hand general network conditions for potential Turing instabilities in large biochemical reaction networks are not well developed. A biochemical reaction network with any number of species where only one species moves is represented by a simple digraph and is modeled by a reaction-diffusion system with non-mass action kinetics. A graph-theoretic condition for potential Turing-Hopf instability that arises when a spatially homogeneous equilibrium loses its stability via a single pair of complex eigenvalues is obtained. This novel graph-theoretic condition is closely related to the negative cycle condition for oscillations in ordinary differential equation models and its generalizations, and requires the existence of a pair of subnetworks, each containing an even number of positive cycles. The technique is illustrated with a double-cycle Goodwin type model.

A biochemically semi-detailed model of auxin-mediated vein formation in plant leaves.

We present here a model intended to capture the biochemistry of vein formation in plant leaves. The model consists of three modules. Two of these modules, those describing auxin signaling and transport in plant cells, are biochemically detailed. We couple these modules to a simple model for PIN (auxin efflux carrier) protein localization based on an extracellular auxin sensor. We study the single-cell responses of this combined model in order to verify proper functioning of the modeled biochemical network. We then assemble a multicellular model from the single-cell building blocks. We find that the model can, under some conditions, generate files of polarized cells, but not true veins.

Feedforward non-Michaelis-Menten mechanism for CO(2) uptake by Rubisco: contribution of carbonic anhydrases and photorespiration to optimization of photosynthetic carbon assimilation.

Rubisco, the most abundant protein serving as the primary engine generating organic biomass on Earth, is characterized by a low catalytic constant (in higher plants approx. 3s(-1)) and low specificity for CO(2) leading to photorespiration. We analyze here why this enzyme evolved as the main carbon fixation engine. The high concentration of Rubisco exceeding the concentration of its substrate CO(2) by 2-3 orders of magnitude makes application of Michaelis-Menten kinetics invalid and requires alternative kinetic approaches to describe photosynthetic CO(2) assimilation. Efficient operation of Rubisco is supported by a strong flux of CO(2) to the chloroplast stroma provided by fast equilibration of bicarbonate and CO(2) and forwarding the latter to Rubisco reaction centers. The main part of this feedforward mechanism is a thylakoidal carbonic anhydrase associated with photosystem II and pumping CO(2) from the thylakoid lumen in coordination with the rate of electron transport, water splitting and proton gradient across the thylakoid membrane. This steady flux of CO(2) limits photosynthesis at saturating CO(2) concentrations. At low ambient CO(2) and correspondingly limited capacity of the bicarbonate pool in the stroma, its depletion at the sites of Rubisco is relieved by utilizing O(2) instead of CO(2), i.e. by photorespiration, a process which supplies CO(2) back to Rubisco and buffers the redox state and energy level in the chloroplast. Thus, the regulation of Rubisco function aims to keep steady non-equilibrium levels of CO(2), NADPH/NADP and ATP/ADP in the chloroplast stroma and to optimize the condition of homeostatic photosynthetic flux of matter and energy.

Dynamics and mechanisms of oscillatory photosynthesis.

We classify mathematical models that can be used to describe photosynthetic oscillations using ideas from nonlinear dynamics, and discuss potential mechanisms for photosynthetic oscillations in the context of this classification. We then turn our attention to recent experiments with leaves transferred to a low CO2 atmosphere which revealed stochastic oscillations with a period of a few seconds. Rubisco is the enzyme that takes both CO2 and O2 as substrates correspondingly for photosynthetic assimilation and for photorespiration. Photosynthesis depletes CO2 and produces O2 while respiration and photorespiration work in the opposite direction, so the product of one process becomes the reactant of the other coupled process. We examine the possibility of oscillations of CO2 and O2 in the leaf in relation to photorespiration. We suggest that in the cell, oscillations with a period of a few seconds, corresponding to the time between photosynthetic CO2 fixation and photorespiratory CO2 release, underlie the dynamics of metabolism in C3 plants.