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.     

An evaluation of plant-borne factors promoting the solubilization of alkaline rock phosphates

Summary Four experiments are described in which the uptake patterns of various crop species were shown to enhance the solubilization of alkaline rock phosphates. In most cases, such enhancements resulted from soil acidification brought about by the alkaline uptake patterns. In one case, the enhancement may have stemmed primarily from a high Ca-uptake rate, which might be responsible for a shift in mass-action equilibria favoring the solubilization of rock phosphates. High uptake of Ca is expected to serve two purposes: 1. it acts in creating alkaline uptake patterns causing a decline in soil pH, and 2. it promotes a shift in mass-action equilibria.     

Species Coexistence and Periodicity in Host-Host-Pathogen Models

Models for the transmission of an infectious disease in one and two host populations with and without self-regulation are analyzed. Many unusual behaviors such as multiple positive equilibria and periodic solutions occur in previous models that use the mass-action (density-dependent) incidence. In contrast, the models formulated using the frequency-dependent (standard) incidence have the behavior of a classic endemic model, since below the threshold, the disease dies out, and above the threshold, the disease persists and the infectious fractions approach an endemic equilibrium. The results given here reinforce previous examples in which there are major differences in behavior between models using mass-action and frequency-dependent incidences.     

Genetically engineered Escherichia coli FBR5: Part I. Comparison of high cell density bioreactors for enhanced ethanol production from xylose

Five reactor systems (free cell batch, free cell continuous, entrapped cell immobilized, adsorbed cell packed bed, and cell recycle membrane reactors) were compared for ethanol production from xylose using Escherichia coli FBR5. In the free cell batch and free cell continuous reactors (continuous stirred tank reactor‐CSTR) productivities of 0.84 gL^?1 h^?1 and 1.77 gL^?1 h^?1 were achieved, respectively. A cell recycle membrane reactor resulted in the highest productivity of 55.56 gL^?1 h^?1, which is an increase of 66‐fold (e.g., 6614%) over the batch reactor. Calcium alginate gel CSTR resulted in a productivity of 2.04 gL^?1 h^?1 whereas adsorbed cell packed bed reactor resulted in a productivity of 4.39 gL^?1 h^?1. In the five reactor systems, ethanol concentrations ranged from 18.9 to 40.30 gL^?1 with metabolic yields from 0.44 to 0.51. Published 2012 American Institute of Chemical Engineers Biotechnol. Prog., 2012     

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.     

Stochastic kinetics of reproduction of virions inside a cell

We analyze intracellular viral kinetics in the framework of the model incorporating viral genome replication, mRNA synthesis and degradation, protein synthesis and degradation, capsid assembly, and virion release from a cell. Due to the existence of the critical concentration of viral capsid proteins and other features of reproduction of virions inside a cell, the kinetics is demonstrated to exhibit three distinct initial stages. Specifically, (i) the exponential growth of the viral genome, mRNA and protein concentrations is followed by (ii) the transient stage to (iii) the steady-state regime. The formation of mature virions starts during the transient stage. Comparison of the kinetics, obtained by using the mass-action law and Monte Carlo (MC) technique, indicates that they are nearly identical during the initial exponential growth of the viral intermediates and also during the steady-state stage. The transition from the initial stage to the steady-state regime occurs however somewhat faster in the determenistic case even if the steady-state populations of virions and genomes are appreciable (e.g., about 250 and 500, respectively).     

Role of substrate inhibition kinetics in enzymatic chemical oscillations.

Two chemical kinetic models are investigated using standard nonlinear dynamics techniques to determine the conditions under which substrate inhibition kinetics can lead to oscillations. The first model is a classical substrate inhibition scheme based on Michaelis-Menten kinetics and involves a single substrate. Only when this reaction takes place in a flow reactor (i.e., both substrate and product are taken to follow reversible flow terms) are oscillations observed; however, the range of parameter values over which such oscillations occur is so narrow it is experimentally unobservable. A second model based on a general mechanism applied to the kinetics of many pH-dependent enzymes is also studied. This second model includes both substrate inhibition kinetics as well as autocatalysis through the activation of the enzyme by hydrogen ion. We find that it is the autocatalysis that is always responsible for oscillatory behavior in this scheme. The substrate inhibition terms affect the steady-state behavior but do not lead to oscillations unless product inhibition or multiple substrates are present; this is a general conclusion we can draw from our studies of both the classical substrate inhibition scheme and the pH-dependent enzyme mechanism. Finally, an analysis of the nullclines for these two models allows us to prove that the nullcline slopes must have a negative value for oscillatory behavior to exist; this proof can explain our results. From our analysis, we conclude with a brief discussion of other enzymes that might be expected to produce oscillatory behavior based on a pH-dependent substrate inhibition mechanism.     

Antagonism and bistability in protein interaction networks

A protein interaction network (PIN) is a set of proteins that modulate one another's activities by regulated synthesis and degradation, by reversible binding to form complexes, and by catalytic reactions (e.g., phosphorylation and dephosphorylation). Most PINs are so complex that their dynamical characteristics cannot be deduced accurately by intuitive reasoning alone. To predict the properties of such networks, many research groups have turned to mathematical models (differential equations based on standard biochemical rate laws, e.g., mass-action, Michaelis-Menten, Hill). When using Michaelis-Menten rate expressions to model PINs, care must be exercised to avoid making inconsistent assumptions about enzyme-substrate complexes. We show that an appealingly simple model of a PIN that functions as a bistable switch is compromised by neglecting enzyme-substrate intermediates. When the neglected intermediates are put back into the model, bistability of the switch is lost. The theory of chemical reaction networks predicts that bistability can be recovered by adding specific reaction channels to the molecular mechanism. We explore two very different routes to recover bistability. In both cases, we show how to convert the original 'phenomenological' model into a consistent set of mass-action rate laws that retains the desired bistability properties. Once an equivalent model is formulated in terms of elementary chemical reactions, it can be simulated accurately either by deterministic differential equations or by Gillespie's stochastic simulation algorithm.     

Ideal free distributions, evolutionary games, and population dynamics in multiple-species environments

In this article, we develop population game theory, a theory that combines the dynamics of animal behavior with population dynamics. In particular, we study interaction and distribution of two species in a two-patch environment assuming that individuals behave adaptively (i.e., they maximize Darwinian fitness). Either the two species are competing for resources or they are in a predator-prey relationship. Using some recent advances in evolutionary game theory, we extend the classical ideal free distribution (IFD) concept for single species to two interacting species. We study population dynamical consequences of two-species IFD by comparing two systems: one where individuals cannot migrate between habitats and one where migration is possible. For single species, predator-prey interactions, and competing species, we show that these two types of behavior lead to the same population equilibria and corresponding species spatial distributions, provided interspecific competition is patch independent. However, if differences between patches are such that competition is patch dependent, then our predictions strongly depend on whether animals can migrate or not. In particular, we show that when species are settled at their equilibrium population densities in both habitats in the environment where migration between habitats is blocked, then the corresponding species spatial distribution need not be an IFD. Thus, when species are given the opportunity to migrate, they will redistribute to reach an IFD (e.g., under which the two species can completely segregate), and this redistribution will also influence species population equilibrial densities. Alternatively, we also show that when two species are distributed according to the IFD, the corresponding population equilibrium can be unstable.     

Optimizing the complexity of phytoplankton functional group modeling: An allometric approach

Elucidating patterns and mechanisms that shape phytoplankton assemblages is a popular area of research for empirical and theoretical ecologists. Despite the daunting complexity of phytoplankton dynamics, much of our current understanding has been based on simple models describing food-web interactions with few differential equations. Skeptical views in the literature raise concerns about the increasing model complexity and advice to seek parsimony rather than simplicity. To address this controversy (simple versus complex models), we propose the introduction of an extra layer of causality into plankton models by connecting algal processes (maximum growth rates, nutrient kinetics, settling velocities, metabolic rates) with species-specific morphological features (cell volume, surface-to-volume ratio, shape). In this study, we demonstrate the capacity of a size-based plankton model to reproduce observed water quality patterns (phosphate, total phosphorus, nitrate, total ammonia, total nitrogen, chlorophyll a, and total zooplankton biomass) in the Hamilton Harbour, Ontario. Consistent with empirical evidence, our modeling analysis showed that small algal species have a distinct competitive advantage in summer epilimnetic environments across the range of cell volume and nutrient loading conditions examined; especially, when they are characterized by higher optimal temperature for growth. Strong top-down pressure mediated by high zooplankton abundance effectively controls the standing biomass of phytoplankton species that can otherwise realize high growth rates under the conditions typically prevailing in the end-of-summer epilimnetic environments (e.g., higher temperature optima, higher tolerance in low water clarity). Under high zooplankton control, the secondary variations of phytoplankton are modulated by the ambient phosphorus levels and the size-based strategies for resources procurement, such as the regulation of nutrient transport kinetics. By contrast, when the summer algal assemblage is released by the zooplankton grazing, the exceedance of critical phytoplankton biomass levels and the likelihood of harmful algal blooms are determined by the multitude of factors that shape inter-specific competition patterns (e.g., relative abundance of competing species, nutrient uptake kinetics). Our study evaluates the strengths and weaknesses of this approach and identifies future directions that would provide operational models founded upon concepts of allometry.     

Transfer Functions for Protein Signal Transduction: Application to a Model of Striatal Neural Plasticity

We present a novel formulation for biochemical reaction networks in the context of protein signal transduction. The model consists of input-output transfer functions, which are derived from differential equations, using stable equilibria. We select a set of “source” species, which are interpreted as input signals. Signals are transmitted to all other species in the system (the “target” species) with a specific delay and with a specific transmission strength. The delay is computed as the maximal reaction time until a stable equilibrium for the target species is reached, in the context of all other reactions in the system. The transmission strength is the concentration change of the target species. The computed input-output transfer functions can be stored in a matrix, fitted with parameters, and even recalled to build dynamical models on the basis of state changes. By separating the temporal and the magnitudinal domain we can greatly simplify the computational model, circumventing typical problems of complex dynamical systems. The transfer function transformation of biochemical reaction systems can be applied to mass-action kinetic models of signal transduction. The paper shows that this approach yields significant novel insights while remaining a fully testable and executable dynamical model for signal transduction. In particular we can deconstruct the complex system into local transfer functions between individual species. As an example, we examine modularity and signal integration using a published model of striatal neural plasticity. The modularizations that emerge correspond to a known biological distinction between calcium-dependent and cAMP-dependent pathways. Remarkably, we found that overall interconnectedness depends on the magnitude of inputs, with higher connectivity at low input concentrations and significant modularization at moderate to high input concentrations. This general result, which directly follows from the properties of individual transfer functions, contradicts notions of ubiquitous complexity by showing input-dependent signal transmission inactivation.     

Validation of fractal-like kinetic models by time-resolved binding kinetics of dansylamide and carbonic anhydrase in crowded media

Kinetic studies of biochemical reactions are typically carried out in a dilute solution that rarely contains anything more than reactants, products, and buffers. In such studies, mass-action-based kinetic models are used to analyze the progress curves. However, intracellular compartments are crowded by macromolecules. Therefore, we investigated the adequacy of the proposed generalizations of the mass-action model, which are meant to describe reactions in crowded media. To validate these models, we measured time-resolved kinetics for dansylamide binding to carbonic anhydrase in solutions crowded with polyethylene glycol and Ficoll. The measured progress curves clearly show the effects of crowding. The fractal-like model proposed by Savageau was used to fit these curves. In this model, the association rate coefficient k_a allometrically depends on concentrations of reactants. We also considered the fractal kinetic model proposed by Schnell and Turner, in which k_a depends on time according to a Zipf-Mandelbrot distribution, and some generalizations of these models. We found that the generalization of the mass-action model, in which association and dissociation rate coefficients are concentration-dependent, represents the preferred model. Other models based on time-dependent rate coefficients were inadequate or not preferred by model selection criteria.     

Mathematical models for hollow-fiber enzyme reactors

Two analytically solved mathematical models are presented for a reactor ystem employing immobilized whole cells as a biocatalyst. The whole cells are entrapped or pumped through the shell side of the dialyzer reactor unit. The reactant mixture is circulated through the cialyzer tube side. Nutrient diffuses across the hollow fiber membrane from the tube side to the shell side, where it reacts to form product, which then back diffuses into the reactant mixture stream. The use of a high recirculation ratio of nutrient through the dialyzer tubes to nutrient feed rate to the entire system, allows the system to be modeled as a continuous-flow stirred-tank reactor. The first analysis details the development of an effectiveness-factor correlation for first- and zero-order kinetics. The second analysis presents the solution to an unsteady-state-system mass balance with Michaelis–menten kinetics.     

Quantum-like interference effect in gene expression: glucose-lactose destructive interference

In this note we illustrate on a few examples of cells and proteins behavior that microscopic biological systems can exhibit a complex probabilistic behavior which cannot be described by classical probabilistic dynamics. These examples support authors conjecture that behavior of microscopic biological systems can be described by quantum-like models, i.e., models inspired by quantum-mechanics. At the same time we do not couple quantum-like behavior with quantum physical processes in bio-systems. We present arguments that such a behavior can be induced by information complexity of even smallest bio-systems, their adaptivity to context changes. Although our examples of the quantum-like behavior are rather simple (lactose-glucose interference in E. coli growth, interference effect for differentiation of tooth stem cell induced by the presence of mesenchymal cell, interference in behavior of PrP(C) and PrP(Sc) prions), these examples may stimulate the interest in systems biology to quantum-like models of adaptive dynamics and lead to more complex examples of nonclassical probabilistic behavior in molecular biology.     

Phagocytosis in Dictyostelium: Nibbling, Eating and Cannibalism

ABSTRACT. Phagocytosis is a highly conserved biological process that serves numerous functions in a wide variety of organisms. Over the past few decades Dictyostelium has proven to be an excellent organism for investigations in cell biology and this is certainly no less the case for a study of phagocytosis. This review examines three distinct phagocytic activities which have been characterized in Dictyostelium. The first, "vegetative phagocytosis," represents the classical eukaryotic microbial uptake of food particles (bacteria). The second, a predatory form of phagocytosis, arises when one species such as Dictyostelium caveatum attacks another species of slime mold, engulfing small pieces of the target prey. This has been termed "cell nibbling." The third phagocytic process is "sexual cannibalistic phagocytosis." In this situation a zygote giant cell, having arisen from the fusion of gametic amoebae, attracts unfused nonzygotic amoebae of the same species and engulfs them as a food source. While cell nibbling has not been actively studied, vegetative and sexual cannibalistic phagocytosis have received varying amounts of attention leading to the idea that some of the elements (e.g., glycoprotein receptors and a G_αs subunit) involved in certain of these phagocytic events may be the same. On the other hand, some unique events (e.g., filopodial induction in prey by D. caveatum ) are also worthy of further investigation. Among other things, the presence of self-nonself recognition, the existence of opsonin-like substances and the presence of signal transduction elements (e.g., an A2-like receptor that negatively modulates sexual phagocytosis) once considered to be extant only in higher organisms suggest that much can be learned about phagocytosis in general by further studies in the classic, eukaryotic microbe Dictyostelium discoideum and related species.     

Models for Bounded Systems with Continuous Dynamics

Summary .  Models for natural nonlinear processes, such as population dynamics, have been given much attention in applied mathematics. For example, species competition has been extensively modeled by differential equations. Often, the scientist has preferred to model the underlying dynamical processes (i.e., theoretical mechanisms) in continuous time. It is of both scientific and mathematical interest to implement such models in a statistical framework to quantify uncertainty associated with the models in the presence of observations. That is, given discrete observations arising from the underlying continuous process, the unobserved process can be formally described while accounting for multiple sources of uncertainty (e.g., measurement error, model choice, and inherent stochasticity of process parameters). In addition to continuity, natural processes are often bounded; specifically, they tend to have nonnegative support. Various techniques have been implemented to accommodate nonnegative processes, but such techniques are often limited or overly compromising. This article offers an alternative to common differential modeling practices by using a bias-corrected truncated normal distribution to model the observations and latent process, both having bounded support. Parameters of an underlying continuous process are characterized in a Bayesian hierarchical context, utilizing a fourth-order Runge–Kutta approximation.     

Trichloroethene and cis‐1,2‐dichloroethene concentration‐dependent toxicity model simulates anaerobic dechlorination at high concentrations. II: Continuous flow and attached growth reactors

A model that was used to describe toxicity from high concentrations of chlorinated aliphatic hydrocarbons (CAHs) on reductively dechlorinating cultures in batch reactors (Sabalowsky and Semprini (in press)) was extended here to simulate observations in continuous flow suspended and attached growth reactors. The reductively dechlorinating anaerobic Evanite subculture (EV‐cDCE) was fed trichloroethene (TCE) and excess electron donor to accumulate cis‐1,2‐dichloroethene (cDCE) in a continuous flow stirred tank reactor (CFSTR); and an attached growth recirculating packed column (RPC). A concentration‐dependent toxicity model used to simulate the results of batch reactors in part I (Sabalowsky and Semprini (in press) Biotechnol Bioeng) also simulated well the observations for the CFSTR and RPC growth modes. The toxicity model incorporates cDCE and TCE toxicity coefficients that directly increase the cell decay coefficient in proportion with cDCE and TCE concentrations. Simulated estimates of the cDCE and TCE toxicity coefficients indicate reductively dechlorinating cells are most sensitive to high concentrations of cDCE and TCE in batch‐fed growth, followed by CFSTR, with attached growth being least sensitive. The greater toxicity of TCE than cDCE, and ratio of the modeled toxicity coefficients, agrees with previously proposed models relating toxicity to partitioning in the cell wall (K_M/B), proportional to octanol‐water partitioning (K_OW) coefficients. Biotechnol. Bioeng. 2010;107: 540–549. © 2010 Wiley Periodicals, Inc.     

Dynamical properties of Discrete Reaction Networks

Reaction networks are commonly used to model the dynamics of populations subject to transformations that follow an imposed stoichiometry. This paper focuses on the efficient characterisation of dynamical properties of Discrete Reaction Networks (DRNs). DRNs can be seen as modeling the underlying discrete nondeterministic transitions of stochastic models of reaction networks. In that sense, a proof of non-reachability in a given DRN has immediate implications for any concrete stochastic model based on that DRN, independent of the choice of kinetic laws and constants. Moreover, if we assume that stochastic kinetic rates are given by the mass-action law (or any other kinetic law that gives non-vanishing probability to each reaction if the required number of interacting substrates is present), then reachability properties are equivalent in the two settings. The analysis of two types of global dynamical properties of DRNs is addressed: irreducibility, i.e., the ability to reach any discrete state from any other state; and recurrence, i.e., the ability to return to any initial state. Our results consider both the verification of such properties when species are present in a large copy number, and in the general case. The necessary and sufficient conditions obtained involve algebraic conditions on the network reactions which in most cases can be verified using linear programming. Finally, the relationship of DRN irreducibility and recurrence with dynamical properties of stochastic and continuous models of reaction networks is discussed.