Linear systems approach to analysis of complex dynamic behaviours in biochemical networks

Central functions in the cell are often linked to complex dynamic behaviours, such as sustained oscillations and multistability, in a biochemical reaction network. Determination of the specific mechanisms underlying such behaviours is important, e.g. to determine sensitivity, robustness, and modelling requirements of given cell functions. In this work we adopt a systems approach to the analysis of complex behaviours in intracellular reaction networks, described by ordinary differential equations with known kinetic parameters. We propose to decompose the overall system into a number of low complexity subsystems, and consider the importance of interactions between these in generating specific behaviours. Rather than analysing the network in a state corresponding to the complex non-linear behaviour, we move the system to the underlying unstable steady state, and focus on the mechanisms causing destabilisation of this steady state. This is motivated by the fact that all complex behaviours in unforced systems can be traced to destabilisation (bifurcation) of some steady state, and hence enables us to use tools from linear system theory to qualitatively analyse the sources of given network behaviours. One important objective of the present study is to see how far one can come with a relatively simple approach to the analysis of highly complex biochemical networks. The proposed method is demonstrated by application to a model of mitotic control in Xenopus frog eggs, and to a model of circadian oscillations in Drosophila. In both examples we are able to identify the subsystems, and the related interactions, which are instrumental in generating the observed complex non-linear behaviours.     

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.     

Shape, Size, and Robustness: Feasible Regions in the Parameter Space of Biochemical Networks

The concept of robustness of regulatory networks has received much attention in the last decade. One measure of robustness has been associated with the volume of the feasible region, namely, the region in the parameter space in which the system is functional. In this paper, we show that, in addition to volume, the geometry of this region has important consequences for the robustness and the fragility of a network. We develop an approximation within which we could algebraically specify the feasible region. We analyze the segment polarity gene network to illustrate our approach. The study of random walks in the parameter space and how they exit the feasible region provide us with a rich perspective on the different modes of failure of this network model. In particular, we found that, between two alternative ways of activating Wingless, one is more robust than the other. Our method provides a more complete measure of robustness to parameter variation. As a general modeling strategy, our approach is an interesting alternative to Boolean representation of biochemical networks.     

Robustness and topology of the yeast cell cycle Boolean network

Yeast cell cycle Boolean network was used as a case study of robustness to protein noise. Robustness was interpreted as involving stability of G1 steady state and sequence of gene expression from cell cycle START to stationary G1. A robustness measure to evaluate robustness strength of a network was proposed. Robust putative networks corresponding to the same steady state and sequence of gene expression of wild-type network were sampled. Architecture of wild-type yeast cell cycle network can be revealed by average topology profile of sampled robust putative networks.     

Kinetic constraints for formation of steady states in biochemical networks

The constraint-based analysis has emerged as a useful tool for analysis of biochemical networks. This work introduces the concept of kinetic constraints. It is shown that maximal reaction rates are appropriate constraints only for isolated enzymatic reactions. For biochemical networks, it is revealed that constraints for formation of a steady state require specific relationships between maximal reaction rates of all enzymes. The constraints for a branched network are significantly different from those for a cyclic network. Moreover, the constraints do not require Michaelis-Menten constants for most enzymes, and they only require the constants for the enzymes at the branching or cyclic point. Reversibility of reactions at system boundary or branching point may significantly impact on kinetic constraints. When enzymes are regulated, regulations may impose severe kinetic constraints for the formation of steady states. As the complexity of a network increases, kinetic constraints become more severe. In addition, it is demonstrated that kinetic constraints for networks with co-regulation can be analyzed using the approach. In general, co-regulation enhances the constraints and therefore larger fluctuations in fluxes can be accommodated in the networks with co-regulation. As a first example of the application, we derive the kinetic constraints for an actual network that describes sucrose accumulation in the sugar cane culm, and confirm their validity using numerical simulations.     

Robust filtering circuit design for stochastic gene networks under intrinsic and extrinsic molecular noises

How to design a robust gene network to tolerate more intrinsic kinetic parameter variations and to attenuate more extrinsic environmental noises to achieve a desired filtering level will be an important topic for systems biology and synthetic biology. At present, there is no good systematic design method to achieve robust gene network design. In this study, a gene network suffering from intrinsic kinetic parameter fluctuations and extrinsic environmental noises is modeled as a Langevin equation with state-dependent stochastic noises. Based on the nonlinear stochastic filtering theory, a systematic gene circuit design method is proposed to make gene networks improve their robustness to tolerate more intrinsic noises and to attenuate extrinsic noises to a prescribed filtering level. The robust gene network design principles have not only yielded a comprehensive design theory of robust gene networks, but also gained valuable insights into the molecular noise filtering of gene networks from the systematic perspective.     

Mathematical identification of critical reactions in the interlocked feedback model

Dynamic simulations are necessary for understanding the mechanism of how biochemical networks generate robust properties to environmental stresses or genetic changes. Sensitivity analysis allows the linking of robustness to network structure. However, it yields only local properties regarding a particular choice of plausible parameter values, because it is hard to know the exact parameter values in vivo. Global and firm results are needed that do not depend on particular parameter values. We propose mathematical analysis for robustness (MAR) that consists of the novel evolutionary search that explores all possible solution vectors of kinetic parameters satisfying the target dynamics and robustness analysis. New criteria, parameter spectrum width and the variability of solution vectors for parameters, are introduced to determine whether the search is exhaustive. In robustness analysis, in addition to single parameter sensitivity analysis, robustness to multiple parameter perturbation is defined. Combining the sensitivity analysis and the robustness analysis to multiple parameter perturbation enables identifying critical reactions. Use of MAR clearly identified the critical reactions responsible for determining the circadian cycle in the Drosophila interlocked circadian clock model. In highly robust models, while the parameter vectors are greatly varied, the critical reactions with a high sensitivity are uniquely determined. Interestingly, not only the per-tim loop but also the dclk-cyc loop strongly affect the period of PER, although the dclk-cyc loop hardly changes its amplitude and it is not potentially influential. In conclusion, MAR is a powerful method to explore wide parameter space without human-biases and to link a robust property to network architectures without knowing the exact parameter values. MAR identifies the reactions critically responsible for determining the period and amplitude in the interlocked feedback model and suggests that the circadian clock intensively evolves or designs the kinetic parameters so that it creates a highly robust cycle.     

Robustness of MetaNet graph models: predicting control of urea production in humans

Urea production in human liver was described by a MetaNet graph, a flowchart-like representation of metabolic pathways that includes parameters for the kinetic constants of the constituent enzymes. Formal operations on the graph facilitate the identification of ligand-binding equilibria that participate in feedback regulation in the network of biochemical reactions. The state of the biochemical network is specified by the concentrations of the intermediates. At any particular time, the influence of an identified locus of regulation is proportional to the respective fractional saturation of the corresponding binding site. Enzymes that make or consume the feedback chemicals share in the control of the strength of the feedback signal in proportion to their fractional saturation. This model predicts control of urea production by the processes that deliver amino groups to the urea cycle enzymes more than by the cycle enzymes themselves. Mitochondrial membrane transport processes are important for transmission of information through the network, but irreversible enzymes and processes far from equilibrium control the strength of the feedback signal. Systematic variation of the parameter values by amounts comparable to the expected variability of their measured values indicated a high probability of invariance in the identities of the predicted control points. The properties of the model are consistent with those of error-tolerant scale-free networks. These results demonstrate the robustness of a MetaNet model's predictions with respect to uncertainties in the values of its parameters.     

Ensemble Modeling for Robustness Analysis in engineering non-native metabolic pathways

Metabolic pathways in cells must be sufficiently robust to tolerate fluctuations in expression levels and changes in environmental conditions. Perturbations in expression levels may lead to system failure due to the disappearance of a stable steady state. Increasing evidence has suggested that biological networks have evolved such that they are intrinsically robust in their network structure. In this article, we presented Ensemble Modeling for Robustness Analysis (EMRA), which combines a continuation method with the Ensemble Modeling approach, for investigating the robustness issue of non-native pathways. EMRA investigates a large ensemble of reference models with different parameters, and determines the effects of parameter drifting until a bifurcation point, beyond which a stable steady state disappears and system failure occurs. A pathway is considered to have high bifurcational robustness if the probability of system failure is low in the ensemble. To demonstrate the utility of EMRA, we investigate the bifurcational robustness of two synthetic central metabolic pathways that achieve carbon conservation: non-oxidative glycolysis and reverse glyoxylate cycle. With EMRA, we determined the probability of system failure of each design and demonstrated that alternative designs of these pathways indeed display varying degrees of bifurcational robustness. Furthermore, we demonstrated that target selection for flux improvement should consider the trade-offs between robustness and performance.     

A Detailed Modular Analysis of Heat-Shock Protein Dynamics under Acute and Chronic Stress and Its Implication in Anxiety Disorders

Physiological and psychological stresses cause anxiety disorders such as depression and post-traumatic stress disorder (PTSD) and induce drastic changes at a molecular level in the brain. To counteract this stress, the heat-shock protein (HSP) network plays a vital role in restoring the homeostasis of the system. To study the stress-induced dynamics of heat-shock network, we analyzed three modules of the HSP90 network-namely trimerization reactions, phosphorylation-dephosphorylation reactions, and the conversion of HSP90 from an open to a closed conformation-and constructed a corresponding nonlinear differential equation model based on mass action kinetics laws. The kinetic parameters of the model were obtained through global optimization, and sensitivity analyses revealed that the most sensitive parameters are the kinase and phosphatase that drive the phosphorylation-dephosphorylation reactions. Bifurcation analysis carried out with the estimated kinetic parameters of the model with stress as bifurcation parameter revealed the occurrence of "mushroom", a type of complex dynamics in which S-shaped and Z-shaped hysteretic bistable forms are present together. We mapped the molecular events responsible for generating the mushroom dynamics under stress and interpreted the occurrence of the S-shaped hysteresis to a normal level of stress, and the Z-shaped hysteresis to the HSP90 variations under acute and chronic stress in the fear conditioned system, and further, we hypothesized that this can be extended to stress-related disorders such as depression and PTSD in humans. Finally, we studied the effect of parameter variations on the mushroom dynamics to get insight about the role of phosphorylation-dephosphorylation parameters in HSP90 network in bringing about complex dynamics such as isolas, where the stable steady states in a bistable system are isolated and separated from each other and not connected by an unstable steady state.     

On the estimation of robustness and filtering ability of dynamic biochemical networks under process delays, internal parametric perturbations and external disturbances

Inherently, biochemical regulatory networks suffer from process delays, internal parametrical perturbations as well as external disturbances. Robustness is the property to maintain the functions of intracellular biochemical regulatory networks despite these perturbations. In this study, system and signal processing theories are employed for measurement of robust stability and filtering ability of linear and nonlinear time-delay biochemical regulatory networks. First, based on Lyapunov stability theory, the robust stability of biochemical network is measured for the tolerance of additional process delays and additive internal parameter fluctuations. Then the filtering ability of attenuating additive external disturbances is estimated for time-delay biochemical regulatory networks. In order to overcome the difficulty of solving the Hamilton Jacobi inequality (HJI), the global linearization technique is employed to simplify the measurement procedure by a simple linear matrix inequality (LMI) method. Finally, an example is given in silico to illustrate how to measure the robust stability and filtering ability of a nonlinear time-delay perturbative biochemical network. This robust stability and filtering ability measurement for biochemical network has potential application to synthetic biology, gene therapy and drug design.     

Ground State Robustness as an Evolutionary Design Principle in Signaling Networks

The ability of an organism to survive depends on its capability to adapt to external conditions. In addition to metabolic versatility and efficient replication, reliable signal transduction is essential. As signaling systems are under permanent evolutionary pressure one may assume that their structure reflects certain functional properties. However, despite promising theoretical studies in recent years, the selective forces which shape signaling network topologies in general remain unclear. Here, we propose prevention of autoactivation as one possible evolutionary design principle. A generic framework for continuous kinetic models is used to derive topological implications of demanding a dynamically stable ground state in signaling systems. To this end graph theoretical methods are applied. The index of the underlying digraph is shown to be a key topological property which determines the so-called kinetic ground state (or off-state) robustness. The kinetic robustness depends solely on the composition of the subdigraph with the strongly connected components, which comprise all positive feedbacks in the network. The component with the highest index in the feedback family is shown to dominate the kinetic robustness of the whole network, whereas relative size and girth of these motifs are emphasized as important determinants of the component index. Moreover, depending on topological features, the maintenance of robustness differs when networks are faced with structural perturbations. This structural off-state robustness, defined as the average kinetic robustness of a network''s neighborhood, turns out to be useful since some structural features are neutral towards kinetic robustness, but show up to be supporting against structural perturbations. Among these are a low connectivity, a high divergence and a low path sum. All results are tested against real signaling networks obtained from databases. The analysis suggests that ground state robustness may serve as a rationale for some structural peculiarities found in intracellular signaling networks.     

Novel recurrent neural network for modelling biological networks: Oscillatory p53 interaction dynamics

Understanding the control of cellular networks consisting of gene and protein interactions and their emergent properties is a central activity of Systems Biology research. For this, continuous, discrete, hybrid, and stochastic methods have been proposed. Currently, the most common approach to modelling accurate temporal dynamics of networks is ordinary differential equations (ODE). However, critical limitations of ODE models are difficulty in kinetic parameter estimation and numerical solution of a large number of equations, making them more suited to smaller systems. In this article, we introduce a novel recurrent artificial neural network (RNN) that addresses above limitations and produces a continuous model that easily estimates parameters from data, can handle a large number of molecular interactions and quantifies temporal dynamics and emergent systems properties. This RNN is based on a system of ODEs representing molecular interactions in a signalling network. Each neuron represents concentration change of one molecule represented by an ODE. Weights of the RNN correspond to kinetic parameters in the system and can be adjusted incrementally during network training. The method is applied to the p53-Mdm2 oscillation system – a crucial component of the DNA damage response pathways activated by a damage signal. Simulation results indicate that the proposed RNN can successfully represent the behaviour of the p53-Mdm2 oscillation system and solve the parameter estimation problem with high accuracy. Furthermore, we presented a modified form of the RNN that estimates parameters and captures systems dynamics from sparse data collected over relatively large time steps. We also investigate the robustness of the p53-Mdm2 system using the trained RNN under various levels of parameter perturbation to gain a greater understanding of the control of the p53-Mdm2 system. Its outcomes on robustness are consistent with the current biological knowledge of this system. As more quantitative data become available on individual proteins, the RNN would be able to refine parameter estimation and mapping of temporal dynamics of individual signalling molecules as well as signalling networks as a system. Moreover, RNN can be used to modularise large signalling networks.     

CADLIVE optimizer: web-based parameter estimation for dynamic models

ABSTRACT: Computer simulation has been an important technique to capture the dynamics of biochemical networks. In most networks, however, few kinetic parameters have been measured in vivo because of experimental complexity. We develop a kinetic parameter estimation system, named the CADLIVE Optimizer, which comprises genetic algorithms-based solvers with a graphical user interface. This optimizer is integrated into the CADLIVE Dynamic Simulator to attain efficient simulation for dynamic models.     

Characterizing multistationarity regimes in biochemical reaction networks

Switch like responses appear as common strategies in the regulation of cellular systems. Here we present a method to characterize bistable regimes in biochemical reaction networks that can be of use to both direct and reverse engineering of biological switches. In the design of a synthetic biological switch, it is important to study the capability for bistability of the underlying biochemical network structure. Chemical Reaction Network Theory (CRNT) may help at this level to decide whether a given network has the capacity for multiple positive equilibria, based on their structural properties. However, in order to build a working switch, we also need to ensure that the bistability property is robust, by studying the conditions leading to the existence of two different steady states. In the reverse engineering of biological switches, knowledge collected about the bistable regimes of the underlying potential model structures can contribute at the model identification stage to a drastic reduction of the feasible region in the parameter space of search. In this work, we make use and extend previous results of the CRNT, aiming not only to discriminate whether a biochemical reaction network can exhibit multiple steady states, but also to determine the regions within the whole space of parameters capable of producing multistationarity. To that purpose we present and justify a condition on the parameters of biochemical networks for the appearance of multistationarity, and propose an efficient and reliable computational method to check its satisfaction through the parameter space.     

RankAggreg, an R package for weighted rank aggregation

Background

Determining the parameters of a mathematical model from quantitative measurements is the main bottleneck of modelling biological systems. Parameter values can be estimated from steady-state data or from dynamic data. The nature of suitable data for these two types of estimation is rather different. For instance, estimations of parameter values in pathway models, such as kinetic orders, rate constants, flux control coefficients or elasticities, from steady-state data are generally based on experiments that measure how a biochemical system responds to small perturbations around the steady state. In contrast, parameter estimation from dynamic data requires time series measurements for all dependent variables. Almost no literature has so far discussed the combined use of both steady-state and transient data for estimating parameter values of biochemical systems.

Results

In this study we introduce a constrained optimization method for estimating parameter values of biochemical pathway models using steady-state information and transient measurements. The constraints are derived from the flux connectivity relationships of the system at the steady state. Two case studies demonstrate the estimation results with and without flux connectivity constraints. The unconstrained optimal estimates from dynamic data may fit the experiments well, but they do not necessarily maintain the connectivity relationships. As a consequence, individual fluxes may be misrepresented, which may cause problems in later extrapolations. By contrast, the constrained estimation accounting for flux connectivity information reduces this misrepresentation and thereby yields improved model parameters.

Conclusion

The method combines transient metabolic profiles and steady-state information and leads to the formulation of an inverse parameter estimation task as a constrained optimization problem. Parameter estimation and model selection are simultaneously carried out on the constrained optimization problem and yield realistic model parameters that are more likely to hold up in extrapolations with the model.     

Stability of Ensemble Models Predicts Productivity of Enzymatic Systems

Stability in a metabolic system may not be obtained if incorrect amounts of enzymes are used. Without stability, some metabolites may accumulate or deplete leading to the irreversible loss of the desired operating point. Even if initial enzyme amounts achieve a stable steady state, changes in enzyme amount due to stochastic variations or environmental changes may move the system to the unstable region and lose the steady-state or quasi-steady-state flux. This situation is distinct from the phenomenon characterized by typical sensitivity analysis, which focuses on the smooth change before loss of stability. Here we show that metabolic networks differ significantly in their intrinsic ability to attain stability due to the network structure and kinetic forms, and that after achieving stability, some enzymes are prone to cause instability upon changes in enzyme amounts. We use Ensemble Modelling for Robustness Analysis (EMRA) to analyze stability in four cell-free enzymatic systems when enzyme amounts are changed. Loss of stability in continuous systems can lead to lower production even when the system is tested experimentally in batch experiments. The predictions of instability by EMRA are supported by the lower productivity in batch experimental tests. The EMRA method incorporates properties of network structure, including stoichiometry and kinetic form, but does not require specific parameter values of the enzymes.     

S-system-based analysis of the robust properties common to many biochemical network models

Robustness is a key feature to characterize the adaptation of organisms to changes in their internal and external environments. A broad range of kinetic or dynamic models of biochemical systems have been developed. Robustness analyses are attractive for exploring some common properties of many biochemical models. To reveal such features, we transform different types of mathematical equations into a standard or intelligible formula and use the multiple parameter sensitivity (MPS) to identify some factors critically responsible for the total robustness to many perturbations. The MPS would be determined by the top quarter of the highly sensitive parameters rather than the single parameter with the maximum sensitivity. The MPS did not show any correlation to the network size. The MPS is closely related to the standard deviation of the sensitivity profile. A decrease in the standard deviation enhanced the total robustness, which shows the hallmark of distributed robustness that many factors (pathways) involve the total robustness.     

Akt versus p53 in a network of oncogenes and tumor suppressor genes regulating cell survival and death

The tumor suppressor protein, p53, and the oncoprotein, Akt, are involved in a cross talk that could be at the core of a cell's control machinery for switching between survival and death. This cross talk is a combination of reciprocally antagonistic pathways emanating from p53 and Akt, and also involves another tumor suppressor gene, PTEN, and another oncogene, Mdm2; such a connected network of cancer-relevant genes must be significant and demands a critical study. The p53-Akt network is shown in this report to possess the potential to exhibit bistability, a phenomenon in which two stable steady states of the system coexist for a fixed set of control parameter values. A hierarchy of qualitative networks and abstract kinetic models are analyzed and simulated on a computer to demonstrate the robustness of the bistable behavior, which, as argued in this study, is a likely candidate mechanism for a cellular survival-death switch. The analysis applies to cells that are neither p53-null nor Akt-null. The models presented here offer experimental predictions on the identity of control parameters of apoptotic thresholds and on network perturbations (including DNA damage and Akt inhibition) that are sufficient to generate switching between pro-survival and pro-death cellular states.     

Computational Multiple Steady States of the Penicillin G Acylase‐Catalyzed Hydrolysis in an Isothermal CFSTR

A family of an enzymatically catalyzed reaction network was studied, which involves the hydrolysis of penicillin G by penicillin G acylase in an isothermal continuous flow stirred tank reactor (CFSTR). This system consisted of 10 coupled non‐linear equations and was found to be capable of exhibiting computational multiple steady states. A set of kinetic parameters determined from the existing experimental data were used to compute a set of rate constants and two corresponding steady states. This suggested that multiple steady states may occur in the system studied. The phenomena of bistability, hysteresis and bifurcation were discussed. Moreover, the capacity of steady state multiplicity was extended to its family of reaction networks.