‘Glocal’ Robustness Analysis and Model Discrimination for Circadian Oscillators

To characterize the behavior and robustness of cellular circuits with many unknown parameters is a major challenge for systems biology. Its difficulty rises exponentially with the number of circuit components. We here propose a novel analysis method to meet this challenge. Our method identifies the region of a high-dimensional parameter space where a circuit displays an experimentally observed behavior. It does so via a Monte Carlo approach guided by principal component analysis, in order to allow efficient sampling of this space. This ‘global’ analysis is then supplemented by a ‘local’ analysis, in which circuit robustness is determined for each of the thousands of parameter sets sampled in the global analysis. We apply this method to two prominent, recent models of the cyanobacterial circadian oscillator, an autocatalytic model, and a model centered on consecutive phosphorylation at two sites of the KaiC protein, a key circadian regulator. For these models, we find that the two-sites architecture is much more robust than the autocatalytic one, both globally and locally, based on five different quantifiers of robustness, including robustness to parameter perturbations and to molecular noise. Our ‘glocal’ combination of global and local analyses can also identify key causes of high or low robustness. In doing so, our approach helps to unravel the architectural origin of robust circuit behavior. Complementarily, identifying fragile aspects of system behavior can aid in designing perturbation experiments that may discriminate between competing mechanisms and different parameter sets.     

Efficient Characterization of Parametric Uncertainty of Complex (Bio)chemical Networks

Parametric uncertainty is a particularly challenging and relevant aspect of systems analysis in domains such as systems biology where, both for inference and for assessing prediction uncertainties, it is essential to characterize the system behavior globally in the parameter space. However, current methods based on local approximations or on Monte-Carlo sampling cope only insufficiently with high-dimensional parameter spaces associated with complex network models. Here, we propose an alternative deterministic methodology that relies on sparse polynomial approximations. We propose a deterministic computational interpolation scheme which identifies most significant expansion coefficients adaptively. We present its performance in kinetic model equations from computational systems biology with several hundred parameters and state variables, leading to numerical approximations of the parametric solution on the entire parameter space. The scheme is based on adaptive Smolyak interpolation of the parametric solution at judiciously and adaptively chosen points in parameter space. As Monte-Carlo sampling, it is “non-intrusive” and well-suited for massively parallel implementation, but affords higher convergence rates. This opens up new avenues for large-scale dynamic network analysis by enabling scaling for many applications, including parameter estimation, uncertainty quantification, and systems design.     

Homoclinic and heteroclinic orbits to a cycle in a tri-trophic food chain

 The asymptotic behavior of a tri-trophic food chain model is studied. The analysis is carried out numerically, by finding both local and global bifurcations of equilibria and limit cycles. The existence of transversal homoclinic orbits to a limit cycle is shown. The appearance of homoclinic orbits, by moving through a homoclinic bifurcation point, is associated with the sudden disappearance of a chaotic attractor. A homoclinic bifurcation curve, which bounds a region of extinction, is continued through a two-dimensional parameter space. Heteroclinic orbits from an equilibrium to a limit cycle are computed. The existence of these heteroclinic orbits has important consequences on the domains of attraction. Continuation of non-transversal heteroclinic orbits through parameter space shows the existence of two codimension-two bifurcations points, where the saddle cycle is non-hyperbolic. The results are summarized by dividing the parameter space in subregions with different asymptotic behavior. Received: 25 February 1998 / Revised version: 19 August 1998     

Estimating the Stochastic Bifurcation Structure of Cellular Networks

High throughput measurement of gene expression at single-cell resolution, combined with systematic perturbation of environmental or cellular variables, provides information that can be used to generate novel insight into the properties of gene regulatory networks by linking cellular responses to external parameters. In dynamical systems theory, this information is the subject of bifurcation analysis, which establishes how system-level behaviour changes as a function of parameter values within a given deterministic mathematical model. Since cellular networks are inherently noisy, we generalize the traditional bifurcation diagram of deterministic systems theory to stochastic dynamical systems. We demonstrate how statistical methods for density estimation, in particular, mixture density and conditional mixture density estimators, can be employed to establish empirical bifurcation diagrams describing the bistable genetic switch network controlling galactose utilization in yeast Saccharomyces cerevisiae. These approaches allow us to make novel qualitative and quantitative observations about the switching behavior of the galactose network, and provide a framework that might be useful to extract information needed for the development of quantitative network models.     

Robustness as a measure of plausibility in models of biochemical networks

Theory, experiment, and observation suggest that biochemical networks which are conserved across species are robust to variations in concentrations and kinetic parameters. Here, we exploit this expectation to propose an approach to model building and selection. We represent a model as a mapping from parameter space to behavior space, and utilize bifurcation analysis to study the robustness of each region of steady-state behavior to parameter variations. The hypothesis that potential errors in models will result in parameter sensitivities is tested by analysis of two models of the biochemical oscillator underlying the Xenopus cell cycle. Our analysis successfully identifies known weaknesses in the older model and suggests areas for further investigation in the more recent, more plausible model. It also correctly highlights why the more recent model is more plausible.     

Stability of discrete age-structured and aggregated delay-difference population models

Sufficiency conditions for local stability are derived for a class of density dependent Leslie matrix models. Four of the recruitment functions in common use in fisheries management are then considered. In two of these oscillating instability can never occur (Beverton and Holt and Cushing forms). In the other two (Deriso-Schnute and Shepherd forms) undamped oscillations are possible within the region of parameter space described here. An algorithm is developed for calculating necessary and sufficient local stability conditions for a simplified form of the general age-structured model. The complete spectrum of stability states (monotonic stability; monotonic instability; oscillating-stable; oscillating-unstable) and the bifurcation periods are given for selected examples of this model. The examples cover a large portion of the parameter space of interest in resource management. It is shown that in perfectly deterministic systems which are observed with error, oscillating instabilities may be missed, and such systems could be erroneously assumed to be stable.     

From global to local: exploring the relationship between parameters and behaviors in models of electrical excitability

Models of electrical activity in excitable cells involve nonlinear interactions between many ionic currents. Changing parameters in these models can produce a variety of activity patterns with sometimes unexpected effects. Further more, introducing new currents will have different effects depending on the initial parameter set. In this study we combined global sampling of parameter space and local analysis of representative parameter sets in a pituitary cell model to understand the effects of adding K ⁺ conductances, which mediate some effects of hormone action on these cells. Global sampling ensured that the effects of introducing K ⁺ conductances were captured across a wide variety of contexts of model parameters. For each type of K ⁺ conductance we determined the types of behavioral transition that it evoked. Some transitions were counterintuitive, and may have been missed without the use of global sampling. In general, the wide range of transitions that occurred when the same current was applied to the model cell at different locations in parameter space highlight the challenge of making accurate model predictions in light of cell-to-cell heterogeneity. Finally, we used bifurcation analysis and fast/slow analysis to investigate why specific transitions occur in representative individual models. This approach relies on the use of a graphics processing unit (GPU) to quickly map parameter space to model behavior and identify parameter sets for further analysis. Acceleration with modern low-cost GPUs is particularly well suited to exploring the moderate-sized (5-20) parameter spaces of excitable cell and signaling models.     

A generalized reaction diffusion model for spatial structure formed by motile cells

A non-linear stability analysis using a multi-scale perturbation procedure is carried out on a model of a generalized reaction diffusion mechanism which involves only a single equation but which nevertheless exhibits bifurcation to non-uniform states. The patterns generated by this model by variation in a parameter related to the scalar dimensions of domain of definition, indicate its capacity to represent certain key morphogenetic features of multicellular systems formed by motile cells.     

Bifurcation discovery tool

MOTIVATION: Biochemical networks often yield interesting behavior such as switching, oscillation and chaotic dynamics. This article describes a tool that is capable of searching for bifurcation points in arbitrary ODE-based reaction networks by directing the user to regions in the parameter space, where such interesting dynamical behavior can be observed. RESULTS: We have implemented a genetic algorithm that searches for Hopf bifurcations, turning points and bistable switches. The software is implemented as a Systems Biology Workbench (SBW) enabled module and accepts the standard SBML model format. The interface permits a user to choose the parameters to be searched, admissible parameter ranges, and the nature of the bifurcation to be sought. The tool will return the parameter values for the model for which the particular behavior is observed. AVAILABILITY: The software, tutorial manual and test models are available for download at the following website: http:/www.sys-bio.org/ under the bifurcation link. The software is an open source and licensed under BSD.     

Inverse bifurcation analysis: application to simple gene systems

Background  

Bifurcation analysis has proven to be a powerful method for understanding the qualitative behavior of gene regulatory networks. In addition to the more traditional forward problem of determining the mapping from parameter space to the space of model behavior, the inverse problem of determining model parameters to result in certain desired properties of the bifurcation diagram provides an attractive methodology for addressing important biological problems. These include understanding how the robustness of qualitative behavior arises from system design as well as providing a way to engineer biological networks with qualitative properties.     

Memory switches in chemical reaction space

Just as complex electronic circuits are built from simple Boolean gates, diverse biological functions, including signal transduction, differentiation, and stress response, frequently use biochemical switches as a functional module. A relatively small number of such switches have been described in the literature, and these exhibit considerable diversity in chemical topology. We asked if biochemical switches are indeed rare and if there are common chemical motifs and family relationships among such switches. We performed a systematic exploration of chemical reaction space by generating all possible stoichiometrically valid chemical configurations up to 3 molecules and 6 reactions and up to 4 molecules and 3 reactions. We used Monte Carlo sampling of parameter space for each such configuration to generate specific models and checked each model for switching properties. We found nearly 4,500 reaction topologies, or about 10% of our tested configurations, that demonstrate switching behavior. Commonly accepted topological features such as feedback were poor predictors of bistability, and we identified new reaction motifs that were likely to be found in switches. Furthermore, the discovered switches were related in that most of the larger configurations were derived from smaller ones by addition of one or more reactions. To explore even larger configurations, we developed two tools: the “bistabilizer,” which converts almost-bistable systems into bistable ones, and frequent motif mining, which helps rank untested configurations. Both of these tools increased the coverage of our library of bistable systems. Thus, our systematic exploration of chemical reaction space has produced a valuable resource for investigating the key signaling motif of bistability.     

The effects of the functional response on the bifurcation behavior of a mite predator–prey interaction model

 A predator–prey interaction model based on a system of differential equations with temperature-dependent parameters chosen appropriately for a mite interaction on apple trees is analyzed to determine how the type of functional response influences bifurcation and stability behavior. Instances of type I, II, III, and IV functional responses are considered, the last of which incorporates prey interference with predation. It is shown that the model systems with the type I, II, and III functional responses exhibit qualitatively similar bifurcation and stability behavior over the interval of definition of the temperature parameter. Similar behavior is found in the system with the type IV functional response at low levels of prey interference. Higher levels of interference are destabilizing, as illustrated by the prevalence of bistability and by the presence of three attractors for some values of the model parameters. All four systems are capable of modeling population oscillations and outbreaks. Received 12 March 1996; received in revised form 25 October 1996     

Comprehensive,Population-Based Sensitivity Analysis of a Two-Mass Vocal Fold Model

Previous vocal fold modeling studies have generally focused on generating detailed data regarding a narrow subset of possible model configurations. These studies can be interpreted to be the investigation of a single subject under one or more vocal conditions. In this study, a broad population-based sensitivity analysis is employed to examine the behavior of a virtual population of subjects and to identify trends between virtual individuals as opposed to investigating a single subject or model instance. Four different sensitivity analysis techniques were used in accomplishing this task. Influential relationships between model input parameters and model outputs were identified, and an exploration of the model’s parameter space was conducted. Results indicate that the behavior of the selected two-mass model is largely dominated by complex interactions, and that few input-output pairs have a consistent effect on the model. Results from the analysis can be used to increase the efficiency of optimization routines of reduced-order models used to investigate voice abnormalities. Results also demonstrate the types of challenges and difficulties to be expected when applying sensitivity analyses to more complex vocal fold models. Such challenges are discussed and recommendations are made for future studies.     

A non-linear analysis for spatial structure in a reaction-diffusion model

A non-linear stability analysis using a multi-scale perturbation procedure is carried out on the practical Thomas reaction-diffusion mechanism which exhibits bifurcation to non-uniform states. The analytical results compare favourably with the numerical solutions. The sequential patterns generated by this model by variations in a parameter related to the reaction-diffusion domain indicate its capacity to represent certain key morphogenetic features required in a recent model by Kauffman for pattern formation in theDrosophila embryo.     

Spontaneous secondary spiking in excitable cells

Kepler & Marder (1993, Biol. Cybern.68, 209-214) proposed a model describing the electrical activity of a crab neuron in which a train of directly induced action potentials is sometimes followed by one or more spontaneous action potentials, referred to as spontaneous secondary spikes. We reduce their five-dimensional model to three dimensions in two different ways in order to gain insight into the mechanism underlying the spontaneous spikes. We then treat a slowly varying current as a parameter in order to give a qualitative explanation of the phenomenon using phase-plane and bifurcation analysis. We demonstrate that a three-dimensional model, consisting of a two-dimensional excitable system plus a slow inward current, is sufficient to produce the behaviour observed in the original model. The exact dynamics of the excitable system are not important, but the relative time constant and amplitude of the slow inward current are crucial. Using the numerical bifurcation analysis package AUTO (Doedel & Kernevez, 1986, AUTO: Software for Continuation and Bifurcation Problems in Ordinary Differential Equations. California Institute of Technology), we compute bifurcation diagrams using the maximum amplitude of the slow inward current as the bifurcation parameter. The full and reduced models have a stable resting potential for all values of the bifurcation parameter. At a critical value of the bifurcation parameter, a stable tonic firing mode arises via a saddle-node of periodics bifurcation. Whether or not the models can exhibit transient or continuous spontaneous spiking depends on their position in parameter space relative to this saddle-node of periodics.     

Instability helps virtual flies to mate

In this paper we perform a bifurcation analysis for a discrete time dynamical system, describing the behavior of a virtual fly, developed by Böddeker and Egelhaaf (2003). Like real blowflies, the virtual counterparts exhibit a dichotomous behavior: they catch small targets but follow big objects at a constant distance. We consider this model for targets on linear and on circular trajectories. Then we transform the system into a ‘‘frozen’’ form, such that the position of the target is fixed. It turns out that the loss of stability of a fixed point in the frozen system due to a Neimark-Sacker bifurcation, explains the dichotomous behavior of the virtual fly.     

A High-Throughput Screening Approach to Discovering Good Forms of Biologically Inspired Visual Representation

While many models of biological object recognition share a common set of “broad-stroke” properties, the performance of any one model depends strongly on the choice of parameters in a particular instantiation of that model—e.g., the number of units per layer, the size of pooling kernels, exponents in normalization operations, etc. Since the number of such parameters (explicit or implicit) is typically large and the computational cost of evaluating one particular parameter set is high, the space of possible model instantiations goes largely unexplored. Thus, when a model fails to approach the abilities of biological visual systems, we are left uncertain whether this failure is because we are missing a fundamental idea or because the correct “parts” have not been tuned correctly, assembled at sufficient scale, or provided with enough training. Here, we present a high-throughput approach to the exploration of such parameter sets, leveraging recent advances in stream processing hardware (high-end NVIDIA graphic cards and the PlayStation 3''s IBM Cell Processor). In analogy to high-throughput screening approaches in molecular biology and genetics, we explored thousands of potential network architectures and parameter instantiations, screening those that show promising object recognition performance for further analysis. We show that this approach can yield significant, reproducible gains in performance across an array of basic object recognition tasks, consistently outperforming a variety of state-of-the-art purpose-built vision systems from the literature. As the scale of available computational power continues to expand, we argue that this approach has the potential to greatly accelerate progress in both artificial vision and our understanding of the computational underpinning of biological vision.     

Explosive route to chaos through a fractal torus in a generalized Lotka-Volterra model

The behavior of a model that generalizes the Lotka-Volterra problem into three dimensions is presented. The results show the analytic derivation of stability diagrams that describe the system's qualitative features. In particular, we show that for a certain value of the bifurcation parameter the system instantly jumps out of a steady state solution into a chaotic solution that portrays a fractal torus in the three-dimensional phase space. This scenario, is referred to as the explosive route to chaos and is attributed to the non-transversal saddle connection type bifurcation. The stability diagrams also present a region in which the Hopf type bifurcation leads to periodic and chaotic solutions. In addition, the bifurcation diagrams reveal a qualitative similarity to the data obtained in the Texas and Bordeaux experiments on the Belousov-Zhabotinskii chemical reaction. The paper is concluded by showing that the model can be useful for representing dynamics associated with biological and chemical phenomena.     

Steady-state solutions in mathematical models of atrial cell electrophysiology and their stability

The steady states of the Fenton-Karma, the Courtemanche and the Nygren cell models were studied by determining the fixed points of the dynamical system describing their cell kinetics. The linear stability of the fixed points was investigated, as well as their response to external stimuli. Symbolic calculations were carried out as far as possible in order to prove the existence of these fixed points. In the Fenton-Karma model, a unique stable fixed point was found, namely the resting state. In contrast, the Courtemanche model had an infinite number of fixed points. A bifurcation diagram was constructed by classifying these fixed points according to a conservation law. Initial conditions were identified, for which the dynamical behavior of the cell was auto-oscillatory. In its original formulation, the Nygren model had no fixed point. After having restored charge conservation, the system was found to have an infinite number of fixed points, resulting in a bifurcation diagram similar to that of the Courtemanche model. The approach proposed in this paper assists in the exploration of the high-dimensional parameter space of the cell models and the identification of the conditions leading to spontaneous pacemaker activity.     

Multistationary and Oscillatory Modes of Free Radicals Generation by the Mitochondrial Respiratory Chain Revealed by a Bifurcation Analysis

The mitochondrial electron transport chain transforms energy satisfying cellular demand and generates reactive oxygen species (ROS) that act as metabolic signals or destructive factors. Therefore, knowledge of the possible modes and bifurcations of electron transport that affect ROS signaling provides insight into the interrelationship of mitochondrial respiration with cellular metabolism. Here, a bifurcation analysis of a sequence of the electron transport chain models of increasing complexity was used to analyze the contribution of individual components to the modes of respiratory chain behavior. Our algorithm constructed models as large systems of ordinary differential equations describing the time evolution of the distribution of redox states of the respiratory complexes. The most complete model of the respiratory chain and linked metabolic reactions predicted that condensed mitochondria produce more ROS at low succinate concentration and less ROS at high succinate levels than swelled mitochondria. This prediction was validated by measuring ROS production under various swelling conditions. A numerical bifurcation analysis revealed qualitatively different types of multistationary behavior and sustained oscillations in the parameter space near a region that was previously found to describe the behavior of isolated mitochondria. The oscillations in transmembrane potential and ROS generation, observed in living cells were reproduced in the model that includes interaction of respiratory complexes with the reactions of TCA cycle. Whereas multistationarity is an internal characteristic of the respiratory chain, the functional link of respiration with central metabolism creates oscillations, which can be understood as a means of auto-regulation of cell metabolism.