Chemical Reaction Systems with Toric Steady States

Mass-action chemical reaction systems are frequently used in computational biology. The corresponding polynomial dynamical systems are often large (consisting of tens or even hundreds of ordinary differential equations) and poorly parameterized (due to noisy measurement data and a small number of data points and repetitions). Therefore, it is often difficult to establish the existence of (positive) steady states or to determine whether more complicated phenomena such as multistationarity exist. If, however, the steady state ideal of the system is a binomial ideal, then we show that these questions can be answered easily. The focus of this work is on systems with this property, and we say that such systems have toric steady states. Our main result gives sufficient conditions for a chemical reaction system to have toric steady states. Furthermore, we analyze the capacity of such a system to exhibit positive steady states and multistationarity. Examples of systems with toric steady states include weakly-reversible zero-deficiency chemical reaction systems. An important application of our work concerns the networks that describe the multisite phosphorylation of a protein by a kinase/phosphatase pair in a sequential and distributive mechanism.     

THE THEORY OF DIFFUSION IN CELL MODELS : II. SOLUTION OF THE STEADY STATE FOR THREE DIFFUSING SUBSTANCES

The differential equations describing diffusion in cell models have been extended to include the simultaneous penetration of water and two salts. These equations have been solved for the steady state. Values for the concentrations in the steady state which may be computed from the equations compare favorably with the experimental values obtained by Osterhout, Kamerling, and Stanley. Moreover, it has been shown elsewhere that the solution for the steady state is essential to a discussion of the volume change or "growth" of phase C in the models and, by analogy, in living cells.     

Logical identification of all steady states: The concept of feedback loop characteristic states

Biological regulatory systems can be described in terms of non-linear differential equations or in logical terms (using an “infinitely non-linear” approximation). Until recently, only part of the steady states of a system could be identified on logical grounds. The reason was that steady states frequently have one or more variable located on a threshold (see below); those steady states were not detected because so far no logical status was assigned to threshold values. This is why we introduced logical scales with values 0,¹θ, 1²θ, 2, ..., in which¹θ,²θ, ... are the logical values assigned to the successive thresholds of the scale. We thus have, in addition to the regular logical states,singular states in which one or more variables is located on a threshold. This permits identifyingall the steady states on logical grounds. It was noticed that each feedback loop (or reunion of disjointed loops) can be characterized by a logical state located at the thresholds at which the variables of the loop operate. This led to the concept ofloop-characteristic state, which, as we will see, enormously simplifies the analysis.The core of this paper is a formal demonstration that among the singular states of a system, only loop-characteristic states can be steady. Reciprocally, given a loop-characteristic state, there are parameter values for which this state is steady; in this case, the loop is effective (i.e. it generates multistationarity if it is a positive loop, homeostasis if it is a negative loop). This not only results in the above-mentioned radical simplification of the identification of the steady states, but in an entirely new view of the relation between feedback loops and steady states.     

Convergence of genotypic frequencies for differential selfing and positive assortative mating at a biallelic locus

For a biallelic model of differential self-fertilization and differential positive assortative mating based on genotype, it is shown that the genotypic frequencies converge for all sets of mating system parameters. Overdominance and underdominance with respect to the parameters are necessary but not sufficient conditions for global convergence to a polymorphic equilibrium and local attractiveness of both the fixation states, respectively. There are cases of overdominance and underdominance for which one fixation state is globally attractive. The relationship of the result to those known from the classical viability selection model are briefly discussed. For the multiallelic version, it is shown that after the first generation all of the homozygote frequencies are always in excess of the corresponding Hardy-Weinberg proportions if at least one homozygote rate of self-fertilization or assortment probability is positive.     

Membrane transport models with fast and slow reactions: General analytical solution for a single relaxation

Summary Membrane transport models are usually expressed on the basis of chemical kinetics. The states of a transporter are related by rate constants, and the time-dependent changes of these states are given by linear differential equations of first order. To calculate the time-dependent transport equation, it is necessary to solve a system of differential equations which does not have a general analytical solution if there are more than five states. Since transport measurements in a complex system rarely provide all the time constants because some of them are too rapid, it is more appropriate to obtain approximate analytical solutions, assuming that there are fast and slow reaction steps. The states of the fast steps are related by equilibrium constants, thus permitting the elimination of their differential equations and leaving only those for the slow steps. With a system having only two slow steps, a single differential equation is obtained and the state equations have a single relaxation. Initial conditions for the slow reactions are determined after the perturbation which redistribute the states related by fast reactions. Current and zero-trans uptake equations are calculated. Curve fitting programs can be used to implement the general procedure and obtain the model parameters.     

Mathematical properties of pump-leak models of cell volume control and electrolyte balance

Homeostatic control of cell volume and intracellular electrolyte content is a fundamental problem in physiology and is central to the functioning of epithelial systems. These physiological processes are modeled using pump-leak models, a system of differential algebraic equations that describes the balance of ions and water flowing across the cell membrane. Despite their widespread use, very little is known about their mathematical properties. Here, we establish analytical results on the existence and stability of steady states for a general class of pump-leak models. We treat two cases. When the ion channel currents have a linear current-voltage relationship, we show that there is at most one steady state, and that the steady state is globally asymptotically stable. If there are no steady states, the cell volume tends to infinity with time. When minimal assumptions are placed on the properties of ion channel currents, we show that there is an asymptotically stable steady state so long as the pump current is not too large. The key analytical tool is a free energy relation satisfied by a general class of pump-leak models, which can be used as a Lyapunov function to study stability.     

Persistence despite perturbations for interacting populations

Two definitions of persistence despite perturbations in deterministic models are presented. The first definition, persistence despite frequent small perturbations, is shown to be equivalent to the existence of a positive attractor i.e. an attractor bounded away from extinction. The second definition, persistence despite rare large perturbations, is shown to be equivalent to permanence i.e. a positive attractor whose basin of attraction includes all positive states. Both definitions set up a natural dichotomy for classifying models of interacting populations. Namely, a model is either persistent despite perturbations or not. When it is not persistent, it follows that all initial conditions are prone to extinction due to perturbations of the appropriate type. For frequent small perturbations, this method of classification is shown to be generically robust: there is a dense set of models for which persistent (respectively, extinction prone) models lies within an open set of persistent (resp. extinction prone) models. For rare large perturbations, this method of classification is shown not to be generically robust. Namely, work of Josef Hofbauer and the author have shown there are open sets of ecological models containing a dense sets of permanent models and a dense set of extinction prone models. The merits and drawbacks of these different definitions are discussed.     

Erratum to: Mathematical properties of pump-leak models of cell volume control and electrolyte balance

Homeostatic control of cell volume and intracellular electrolyte content is a fundamental problem in physiology and is central to the functioning of epithelial systems. These physiological processes are modeled using pump-leak models, a system of differential algebraic equations that describes the balance of ions and water flowing across the cell membrane. Despite their widespread use, very little is known about their mathematical properties. Here, we establish analytical results on the existence and stability of steady states for a general class of pump-leak models. We treat two cases. When the ion channel currents have a linear current-voltage relationship, we show that there is at most one steady state, and that the steady state is globally asymptotically stable. If there are no steady states, the cell volume tends to infinity with time. When minimal assumptions are placed on the properties of ion channel currents, we show that there is an asymptotically stable steady state so long as the pump current is not too large. The key analytical tool is a free energy relation satisfied by a general class of pump-leak models, which can be used as a Lyapunov function to study stability.     

Threshold behaviour of a SIR epidemic model with age structure and immigration

We consider a SIR age-structured model with immigration of infectives in all epidemiological compartments; the population is assumed to be in demographic equilibrium between below-replacement fertility and immigration; the spread of the infection occurs through a general age-dependent kernel. We analyse the equations for steady states; because of immigration of infectives a steady state with a positive density of infectives always exists; however, a quasi-threshold theorem is proved, in the sense that, below the threshold, the density of infectives is close to 0, while it is away from 0, above the threshold; furthermore, conditions that guarantee uniqueness of steady states are obtained. Finally, we present some numerical examples, inspired by the Italian demographic situation, that illustrate the threshold-like behaviour, and other features of the stationary solutions and of the transient. Supported in part by FIRB project RBAU01K7M2 “Metodi dell’Analisi Matematica in Biologia, Medicina e Ambiente” of the Italian Ministero Istruzione Università e Ricerca     

On local bifurcations in neural field models with transmission delays

Neural field models with transmission delays may be cast as abstract delay differential equations (DDE). The theory of dual semigroups (also called sun-star calculus) provides a natural framework for the analysis of a broad class of delay equations, among which DDE. In particular, it may be used advantageously for the investigation of stability and bifurcation of steady states. After introducing the neural field model in its basic functional analytic setting and discussing its spectral properties, we elaborate extensively an example and derive a characteristic equation. Under certain conditions the associated equilibrium may destabilise in a Hopf bifurcation. Furthermore, two Hopf curves may intersect in a double Hopf point in a two-dimensional parameter space. We provide general formulas for the corresponding critical normal form coefficients, evaluate these numerically and interpret the results.     

Steady-State Kinetic Modeling Constrains Cellular Resting States and Dynamic Behavior

A defining characteristic of living cells is the ability to respond dynamically to external stimuli while maintaining homeostasis under resting conditions. Capturing both of these features in a single kinetic model is difficult because the model must be able to reproduce both behaviors using the same set of molecular components. Here, we show how combining small, well-defined steady-state networks provides an efficient means of constructing large-scale kinetic models that exhibit realistic resting and dynamic behaviors. By requiring each kinetic module to be homeostatic (at steady state under resting conditions), the method proceeds by (i) computing steady-state solutions to a system of ordinary differential equations for each module, (ii) applying principal component analysis to each set of solutions to capture the steady-state solution space of each module network, and (iii) combining optimal search directions from all modules to form a global steady-state space that is searched for accurate simulation of the time-dependent behavior of the whole system upon perturbation. Importantly, this stepwise approach retains the nonlinear rate expressions that govern each reaction in the system and enforces constraints on the range of allowable concentration states for the full-scale model. These constraints not only reduce the computational cost of fitting experimental time-series data but can also provide insight into limitations on system concentrations and architecture. To demonstrate application of the method, we show how small kinetic perturbations in a modular model of platelet P2Y₁ signaling can cause widespread compensatory effects on cellular resting states.     

Seed dispersal in a patchy environment with global age dependence

A model of seed population dynamics proposed by S. A. Levin, A. Hastings, and D. Cohen is presented and analyzed. With the environment considered as a mosaic of patches, patch age is used along with time as an independent variable. Local dynamics depend not only on the local state, but also on the global environment via dispersal modelled by an integral over all patch ages. Basic technical properties of the time varying solutions are examined; necessary and sufficient conditions for nontrivial steady states are given; and general sufficient conditions for global asymptotic stability of these steady states are established. Primary tools of analysis include a hybrid Picard iteration, fixed point methods, monotonicity of solution structure, and upper and lower solutions for differential equations.This work was supported in part by National Science Foundation Grants MCS-7903497 and MCS-790349701     

Maximum sustainable yield of populations with continuous age-structure.

We address the problem of finding the harvesting policy that will maximize the yield and maintain a population in a steady state. The population is characterized by continuous age classes and therefore follows differential equations. Here, we assume that the equations are linear (no density dependence). Two possible constraints are considered: either recruitment or total population are fixed to a constant. Under these conditions, the optimal policy is to harvest the fraction theta of a younger age class ? and to harvest totally an older age class b. The optimal solution (theta, ?, b) can be calculated explicitly if the fecundity and mortality schedules are given. The solution is compared to the simpler strategy of harvesting all individuals beyond a single age class a. It is shown that the latter strategy can be much less profitable than harvesting two age classes because it cannot take account of the different values of individuals according to their age.     

Pheromone trapping models for pest control: Effects of mating patterns and immigration

Several models are presented which examine pest population behaviour with the release of female sex pheromones for the attraction and annihilation of males. These models include male polygamy and female monogamy, various mating frequencies, delayed mating of females, immigration of one or all individual types, and differential survivorship of males and females. In all the models there are two steady states, a stable s.s. at the origin and an unstable s.s. in the positive domain for a given value of pheromone release rate. In all the models, control relies on the reduced ability of males to fertilize virgin females following trapping and male annihilation. As such, control is very sensitive to mating frequency, being very difficult when males mate frequently. Control is also very difficult with the immigration of even a moderate number of fertilized females. Control is much easier when mating is delayed, especially if survivorship is low, or with density dependent population regulation.     

Smooth bistable S-systems

S-systems have been used as models of biochemical systems for over 30 years. One of their hallmarks is that, although they are highly non-linear, their steady states are characterised by linear equations. This allows streamlined analyses of stability, sensitivities and gains as well as objective, mathematically controlled comparisons of similar model designs. Regular S-systems have a unique steady state at which none of the system variables is zero. This makes it difficult to represent switching phenomena, as they occur, for instance, in the expression of genes, cell cycle phenomena and signal transduction. Previously, two strategies were proposed to account for switches. One was based on a technique called recasting, which permits the modelling of any differentiable non-linearities, including bistability, but typically does not allow steady-state analyses based on linear equations. The second strategy formulated the switching system in a piece-wise fashion, where each piece consisted of a regular S-system. A representation gleaned from a simplified form of recasting is proposed and it is possible to divide the characterisation of the steady states into two phases, the first of which is linear, whereas the other is non-linear, but easy to execute. The article discusses a representative pathway with two stable states and one unstable state. The pathway model exhibits strong separation between the stable states as well as hysteresis.     

Variable elimination in post-translational modification reaction networks with mass-action kinetics

We define a subclass of chemical reaction networks called post-translational modification systems. Important biological examples of such systems include MAPK cascades and two-component systems which are well-studied experimentally as well as theoretically. The steady states of such a system are solutions to a system of polynomial equations. Even for small systems the task of finding the solutions is daunting. We develop a mathematical framework based on the notion of a cut (a particular subset of species in the system), which provides a linear elimination procedure to reduce the number of variables in the system to a set of core variables. The steady states are parameterized algebraically by the core variables, and graphical conditions for when steady states with positive core variables imply positivity of all variables are given. Further, minimal cuts are the connected components of the species graph and provide conservation laws. A criterion for when a (maximal) set of independent conservation laws can be derived from cuts is given.     

Mapping parameter spaces of biological switches

Since the seminal 1961 paper of Monod and Jacob, mathematical models of biomolecular circuits have guided our understanding of cell regulation. Model-based exploration of the functional capabilities of any given circuit requires systematic mapping of multidimensional spaces of model parameters. Despite significant advances in computational dynamical systems approaches, this analysis remains a nontrivial task. Here, we use a nonlinear system of ordinary differential equations to model oocyte selection in Drosophila, a robust symmetry-breaking event that relies on autoregulatory localization of oocyte-specification factors. By applying an algorithmic approach that implements symbolic computation and topological methods, we enumerate all phase portraits of stable steady states in the limit when nonlinear regulatory interactions become discrete switches. Leveraging this initial exact partitioning and further using numerical exploration, we locate parameter regions that are dense in purely asymmetric steady states when the nonlinearities are not infinitely sharp, enabling systematic identification of parameter regions that correspond to robust oocyte selection. This framework can be generalized to map the full parameter spaces in a broad class of models involving biological switches.     

Oscillations in Biochemical Reaction Networks Arising from Pairs of Subnetworks

Biochemical reaction models show a variety of dynamical behaviors, such as stable steady states, multistability, and oscillations. Biochemical reaction networks with generalized mass action kinetics are represented as directed bipartite graphs with nodes for species and reactions. The bipartite graph of a biochemical reaction network usually contains at least one cycle, i.e., a sequence of nodes and directed edges which starts and ends at the same species node. Cycles can be positive or negative, and it has been shown that oscillations can arise as a result of either a positive cycle or a negative cycle. In earlier work it was shown that oscillations associated with a positive cycle can arise from subnetworks with an odd number of positive cycles. In this article we formulate a similar graph-theoretic condition, which generalizes the negative cycle condition for oscillations. This new graph-theoretic condition for oscillations involves pairs of subnetworks with an even number of positive cycles. An example of a calcium reaction network with generalized mass action kinetics is discussed in detail.