- Methodology article
- Open Access
Validation of a model of the GAL regulatory system via robustness analysis of its bistability characteristics
- Luca Salerno†^{1},
- Carlo Cosentino†^{1}Email author,
- Alessio Merola^{1},
- Declan G Bates^{2} and
- Francesco Amato^{1}
https://doi.org/10.1186/1752-0509-7-39
© Salerno et al.; licensee BioMed Central Ltd. 2013
- Received: 27 September 2012
- Accepted: 26 April 2013
- Published: 17 May 2013
Abstract
Background
In Saccharomyces cerevisiæ, structural bistability generates a bimodal expression of the galactose uptake genes (GAL) when exposed to low and high glucose concentrations. This indicates that yeast cells can decide between using either the limited amount of glucose or growing on galactose under changing environmental conditions. A crucial requirement for any plausible mechanistic model of this system is that it reproduces the robustness of the bistable response observed in vivo against inter-individual parametric variability and fluctuating environmental conditions.
Results
We show how a control-theoretic analysis of the robustness of a model of the GAL regulatory network may be used to establish the model’s plausibility in characterizing the persistent memory of different carbon sources, without the need for extensive simulations. Chemical Reaction Network Theory is used to establish that the proposed network model is compatible with structural bistability. The robustness of each of the two operative conditions against fluctuations of the species concentrations is demonstrated by studying the Domains of Attraction of the corresponding equilibrium points. Finally, we use a global robustness analysis method based on Semi-Definite Programming to evaluate the modification of the bistable steady states induced by multiple parametric variations throughout bounded regions of the parameter space.
Conclusions
Our analysis provides convincing evidence for the robustness, and hence plausibility, of the GAL regulatory network model. The proposed workflow also demonstrates the power of analytical methods from control theory to provide a direct quantitative characterization of the dynamics of multistable biomolecular regulatory systems without recourse to extensive computer simulations.
Keywords
- Galactose network
- Bistability
- Robustness
- Domain of attraction
- Bifurcation
- Local sensitivity
- Global sensitivity
Background
Although yeasts, in common with most cellular organisms, can derive energy from a variety of different molecules, glucose is well-known to be their preferred source, because it provides more energy than any other saccharide. Therefore, yeasts have evolved a complex genetic network to make sure they can consume as much glucose as possible when it is available [1]. In [2], the authors experimentally investigated the regulation of galactose metabolism in S. Cerevisiæ, which is mediated by several positive and negative feedback loops acting at the transcriptional level. To probe the system for multistability, two identical cell populations were grown on different media, with and without galactose, respectively. In engineering terms, this amounts to initializing the system at two different operating conditions. Starting from these conditions, the two populations were then exposed to identical galactose concentrations for a period long enough to guarantee the attainment of steady-state conditions. For intermediate levels of the input (galactose concentrations), the two populations settled on different steady states, thus confirming the multistable nature of the system. These and other experimental results have revealed that the GAL system exhibits bistable dynamics and that such bistability generates a persistent memory of the type of carbon source consumed by the cell in the past.
In previous work, classical mathematical tools, such as bifurcation analysis, have been used to examine the dynamics of the GAL regulatory network, see e.g. [3, 4]. Recently, we showed, using a control-theoretic analysis, that the GAL network simplified mass-action models proposed in [5, 6], do not reproduce the bistable behavior exhibited by the experimental studies of Acar and coworkers; this finding motivated us to propose a new model of the GAL system, [7]. In this paper we extend this model and provide a thorough characterization of its dynamical properties, with the aim of validating it as a plausible mechanistic explanation of the persistent memory property. Our approach starts with the analysis of the model’s bistable dynamics as a structural property, arising from the topology of the reaction network. Afterwards, we focus on the study of the robustness of bistability both against fluctuations of the concentrations of the molecular species, caused by endogenous stochastic noise or by exogenous perturbations, and in the face of parametric uncertainties. The principle underpinning these analyses is that the quality of a model cannot be solely evaluated by its capability to reproduce a particular set of experimental measurements. Indeed, a common problem in modeling biological networks is that alternative, structurally different models can fit experimental data equally well [8]. In order to represent a plausible mechanistic description of a biological phenomenon, the model must also replicate an essential feature of biological systems, that is robustness against inter-individual parametric variability and in vivo fluctuating environmental conditions [9–12].
Our characterization of the robustness properties of the model starts with an analysis of the Domains of Attractions (DA’s) of the bistable system. Roughly speaking, the DA of an equilibrium point x_{ e } is a region $\mathcal{D}$ in the state space, such that ${x}_{e}\in \mathcal{D}$ and every state trajectory crossing $\mathcal{D}$ converges asymptotically to x_{ e }. DA analysis is crucial for establishing whether the proposed model provides a plausible explanation of the phenomenon under investigation, since the system is actually able to operate around a given equilibrium point with some degree of robustness in the face of both intrinsic stochastic noise and exogenous perturbations only if that equilibrium point possesses a nontrivial DA. Note that the estimation of the DA is, in general, a difficult problem for systems of nontrivial dimension. In our approach we show how, for any mass-action model, it is possible to apply a convex optimization-based method, devised in a purely theoretical context by our group in [13, 14], to test whether an assigned polytopic subset of the state space belongs to the DA of an equilibrium point.
We next consider the robustness of the model’s bistable dynamics in the face of uncertain parameter values. Many examples can be found in the literature of studies applying local sensitivity and bifurcations analysis as tools for characterizing the parametric robustness of biological systems, e.g. [15–17]; however, these tools suffer from a significant limitation due to their inability to take into account more than one or two simultaneous parameters variation at the same time. Multi-parametric sensitivity analysis of biomodels is typically performed by resorting to extensive sampling of the admissible parameter space, [18, 19], which requires a large large computational effort and can only provide probabilistic conclusions. To overcome these limitations, besides applying local sensitivity and bifurcations analysis, we employ a global sensitivity analysis method proposed in [20, 21]. This method is aimed at computing an outer approximation of the region of the state space that contains all the equilibrium points of a given biosystem for all admissible values of the parameters. In our analysis, we devise a straightforward way to adapt this method to provide robustness certificates for bistability in the face of parametric uncertainty.
Thus, beyond our primary goal of validating a new model of the bistable GAL regulatory network, we also present what should be a widely applicable and effective procedure for investigating the plausibility of dynamical models of multistable biomolecular circuits, without recourse to large-scale numerical simulations.
Results
A new model of the GAL regulatory network in S. Cerevisiæ
The mathematical model of this regulatory network considered here is based on mass-action kinetics and represents an extended version of the model proposed in [7]. In the new version of the model, the reaction of reversible dissociation of the complex Gal80p:Gal3p* is explicitly included, since this was found to be essential to ensure robustness of the bistable dynamics, according to the analysis procedure that will be illustrated in the next sections. Moreover, in the new model also the Gal1p protein dynamics are taken into account, since the concentration of this protein is taken as a measure of Gal4p activity in the experiments reported in [2]. model consists of the following set of nonlinear ordinary differential equations (ODE’s),
Steady states of the mass-action model (1), with the parameters values given in Table2
x _{ e1 } | Species | Description | x _{ e2 } |
---|---|---|---|
172.8212 | G _{3} | Gal3p protein | 2711.1839 |
172.8208 | G _{int} | internalized galactose | 2711.2003 |
1.0 | G _{3a} | active Gal3p protein | 318.5443 |
1.0 | G _{4} | Gal4p protein | 19.9479 |
1.0 | G _{80} | Gal80p protein | 0.3061 |
21.0604 | G _{4,80} | Gal4p:Gal80p complex | 2.1126 |
7.5945 | G _{80,3a} | Gal80p:Gal3p active complex | 589.1342 |
1.0 | G _{2} | Gal2p protein | 19.9479 |
1.0 | G _{1} | Gal1p protein | 19.9479 |
1000.0 | G _{ e x } | external galactose | 1000.0 |
A set of parameters values that renders system (1) bistable
Parameter | Value | Parameter | Value |
---|---|---|---|
k _{1} | 0.1814 [ μ M·h] | k _{11} | 85.8185 [1/h] |
k _{2} | 8.4586E-4 [1/ μ M·h] | k _{r 11} | 4.7482E-2 [1/ μ M·h] |
k _{r 2} | 16.6691 [1/h] | μ _{12} | 1.0 [1/h] |
μ _{3} | 1.0 [1/h] | μ _{13} | 1.0 [1/h] |
k _{4} | 3.0749 [1/ μ M·h] | μ _{14} | 1.0 [1/h] |
k _{r 4} | 0.1317 [1/h] | μ _{15} | 1.0 [1/h] |
k _{5} | 22.0604 [ μ M/h] | μ _{16} | 1.0 [1/h] |
μ _{6} | 1.0 [1/h] | μ _{17} | 1.0 [1/h] |
k _{7} | 29.6549 [1/h] | μ _{18} | 1.0 [1/h] |
k _{8} | 181.4157 [1/h] | k _{19} | 36.6342 [1/ μ M·h] |
k _{9} | 1.0 [1/h] | k _{r 19} | 222.0536 [1/h] |
k _{10} | 1.0 [1/h] |
where the species concentrations, namely the state variables listed in Table 1, are denoted by x:=(G_{3}G_{int}G_{3a}G_{4}G_{80}G_{4,80}G_{80,3a}G_{2}G_{1})^{ T }, the parameters by $p:={({k}_{1}{k}_{2}\cdots {k}_{r19})}^{T}\in {\mathbb{R}}^{23}$ (the full list of parameters is reported in Table 2), $N\in {\mathbb{R}}^{9\times 23}$ is the stoichiometric matrix and $v(x,p)\in {\mathbb{R}}^{23}$ is the vector of reaction rates. Since the parameters are inherently positive, positive values of x will result in positive values of v(x,p), i.e. $x\in {\mathbb{R}}_{+}^{9}\Rightarrow v(x,p)\in {\mathbb{R}}_{+}^{23}$. In the next section, we present the results of our analysis of the bistable dynamical properties of this model.
Structural analysis of the proposed model’s network topology confirms bistable dynamics
The first step of our procedure consists in the analysis of the topology of the regulatory network, to determine whether its structure can admit a bistable behavior. Subsequently, we will determine a possible realization (i.e. a set of parameter values) of model (1) that exhibits bistability.
The persistence of cellular memory exhibited by the galactose regulatory network is a system-level property which results from the interactions of several species in multiple nested feedback loops. Two coupled positive feedback loops involve the galactose permease Gal2p and the signaling protein Gal3p, while a negative feedback loop involves the inhibitor Gal80p. Recall that, according to [22], the existence of a positive-feedback loop, or a mutually inhibitory, double-negative-feedback loop, is a necessary condition for the occurrence of multistability.
The existence of a solution to Eqs. (1b). Such scaling yields an equilibrium concentration of 1mM for G_{ e x }; note also that it does not affect other equations, since the kinetic constant k_{1} does not appear elsewhere in the model.
Characterization of the domains of attraction confirms robustness of the bistable equilibria
Subsequently to the determination of the asymptotically stable equilibrium points, a primary goal in the characterization of the behavior of a system is that of estimating the DA’s of such points. Accurate estimates of the DA’s provide valuable information about the ability of a system to reject perturbations driving the system away from its steady state condition. At the same time, the boundaries of the DA’s constitute the concentration thresholds for the activation of the switching mechanism between different operative conditions.
The methodology proposed in [13], which allows to check whether an assigned box in the state space belongs to the DA of an equilibrium, has been employed in our study. It is worth noticing that the main result of [13] leads to a Linear Matrix Inequality (LMI) feasibility problem, which can be solved efficiently via off-the-shelf numerical algorithms.
In order to find the largest possible estimates of the DA’s of x_{e 1} and x_{e 2}, namely ${\stackrel{~}{\mathcal{D}}}_{1}$ and ${\stackrel{~}{\mathcal{D}}}_{2}$, our procedure takes two small initial polytopic regions, surrounding the equilibrium points, and then iteratively stretches them along the different dimensions of the state space until the feasibility conditions are no longer verified, thus obtaining two inner approximations of the DA’s.
Local and global analysis confirm robustness of the bistable equilibria to parametric uncertainty
In this section we provide further support for the plausibility of the proposed model of the GAL regulatory system by characterizing its robustness with respect to parametric uncertainties.
The underlying principle is that, in view of the large inter-individual variability of biochemical parameters, for a model to be considered plausible it is not sufficient to reproduce the qualitative behavior of the biological system for just one set of parameter values; instead, this behavior must be exhibited over a nontrivial subset of the parameters space.
First, a classical sensitivity analysis is performed by employing the method discussed in [27]: the state variables ODE’s are coupled to the equations of sensitivity.
This allows us to compute a numerical solution to the whole set of equations, thus simultaneously obtaining both the state variables and the associated sensitivity coefficients.
Unfortunately, the above methods can efficiently evaluate changes in steady state concentrations only for low-dimensional parameter variations. In [20], a global sensitivity method, named bioSDP, is proposed to evaluate the effect of multiple (simultaneous) parameter variations on the system’s dynamics. More specifically, this method can be used to compute some bounds on the maximum variation of the equilibrium points induced by changes of the parameter values. The computation is based on the solution of a dual problem (see Methods for more details): given a subset of the state space, say Ω, the bioSDP algorithm is able to certificate that, for any admissible realization of the parameter vector p, Ω does not contain any equilibrium point.
The bioSDP algorithm takes as inputs the admissible range of variation of the parameters, defined as a box ${\mathcal{\mathcal{B}}}_{p}$ in the parameter space, and an initial outer approximation, in the form of a box ${\stackrel{~}{\mathcal{S}}}^{0}$, of the subset ${\mathcal{X}}_{e}$ of all the admissible equilibrium points of system (1) subject to $p\in {\mathcal{\mathcal{B}}}_{p}$.
Then, these outer boundaries are iteratively narrowed by applying a bisection algorithm. As a result, the state space is partitioned in one or more subsets containing all the equilibrium points that fall inside the initial search subspace ${\stackrel{~}{\mathcal{S}}}^{0}$. In fact, due to the computational burden, the bisection algorithm can only be applied to systems of dimension less or equal than three (see for example Figure 5 in [21], where a two-dimensional system is analyzed). For higher-order systems like ours, the algorithm resorts to a box shrinkage procedure, i.e., it just tries to reduce the size of the initial box as far as possible.
Due to the above limitations, in our case the method proposed in [20] would not be able to distinguish the two distinct steady state subsets. To overcome this issue, we devise a strategy that leverages the bioSDP algorithm but, instead of computing one set containing all the equilibrium points, aims to separately compute two distinct robust steady state subsets ${\stackrel{~}{\mathcal{S}}}_{1}$ and ${\stackrel{~}{\mathcal{S}}}_{2}$, which define the boundaries for the variation of x_{e 1} and x_{e 2}, respectively. Thus, we need two initial outer approximation subsets, let us denote them by ${\stackrel{~}{\mathcal{S}}}_{1}^{0}$ and ${\stackrel{~}{\mathcal{S}}}_{2}^{0}$, respectively.
Guessing two good initial outer approximations would in general turn out to be a daunting task for systems of nontrivial dimension. In our case, exploiting the previous analysis and by virtue of continuity arguments, we surmise that, for small-enough variations of the parameter values, the DA’s represent good initial guesses. Thus, we let ${\stackrel{~}{\mathcal{S}}}_{i}^{0}={\rho}_{i}{\stackrel{~}{\mathcal{D}}}_{i}$, i=1,2, where ρ_{ i }>0 is a scaling factor and apply the bioSDP algorithm separately to these two initial boxes, with ρ_{ i }=1; if the algorithm does not find a solution, it is re-applied iteratively with different values of ρ_{ i } in a predefined interval [ρ_{min},ρ_{max}], until a solution is found.
Note that, setting ${\stackrel{~}{\mathcal{S}}}_{i}^{0}$ as the initial search space for the bioSDP algorithm we are focusing the analysis on those equilibrium points that belong to a neighborhood of x_{ e i }, instead of searching for all the equilibrium points. Performing this analysis separately, first on x_{e 1} and then on x_{e 2}, enables us to ascertain whether the bistability is preserved against parametric perturbations: the answer is affirmative if we are able to compute two disjoint robust steady state subsets, i.e., ${\stackrel{~}{\mathcal{S}}}_{1}\cap {\stackrel{~}{\mathcal{S}}}_{2}=\varnothing $. If this problem is feasible for an assigned parameter box ${\mathcal{\mathcal{B}}}_{p}$, then we are guaranteed that bistability is preserved for all p belonging to ${\mathcal{\mathcal{B}}}_{p}$.
To alleviate the computational burden of the procedure, the multi-parametric sensitivity analysis has been limited to the parameters subset Θ:={k_{1},k_{2},k_{5},k_{7},k_{8},k_{9},μ_{13},μ_{16},μ_{17}}, which, according to the local sensitivity analysis, have a major influence on the system dynamics (see Figure 4). The robustness has been evaluated against increasingly larger ranges of parameter variations, corresponding to ±2%, ±5%, ±10%, ±20%, ±30% and ±50%, with respect to the nominal values given in Table 2. Figure 2 displays the computed robust steady state boxes for the various cases. The bistable behavior of the GAL regulatory network is guaranteed for parametric variations up to ±20% with respect to the nominal parameters value. For such uncertainty values, indeed, the computed subsets ${\stackrel{~}{\mathcal{S}}}_{1}$ and ${\stackrel{~}{\mathcal{S}}}_{2}$ are still disjoint, since the intervals of G _{3a} and G _{80,3a} are not overlapping. For parametric variations of ±30% or more, the intersection of the two subsets is no longer empty (see Figures 2, panels g and h); in the latter case, it is no longer possible to guarantee that the system preserves bistability for all admissible parameter values.
Discussion and conclusions
Robustness, intended as the capability to cope with fluctuations of the molecular species concentrations, caused by endogenous and exogenous perturbations, and to preserve biological functions despite inter-individual variability of kinetic parameters, is a key feature of biological systems. This fundamental feature poses an important challenge when trying to describe biological phenomena by means of mechanistic mathematical models: a set of differential/algebraic equations which, for some value of the parameters, interpolates experimental data, cannot be considered a plausible model if it does not possess the aforementioned robustness properties. Recognizing the power of these arguments as tools for testing novel biomodels, and with the aim of supporting the validity of our proposed model of the GAL regulatory network, we have devised an analytical procedure which can be exploited to investigate the robustness properties of biomodels of bistable biological systems. The procedure exploits several complementary methods for the analysis of nonlinear quadratic systems (i.e., mass-action models), devised both by our group and by other authors, and its effectiveness has been demonstrated by applying it to thoroughly characterize the robustness of bistability for a new model of the galactose metabolism regulatory system.
The procedure consists of three phases: in the first phase, the properties of the nominal system (i.e., parameters values are assumed to be certain) are investigated, since the first requirement is that the reaction network is structurally compatible with the existence of multiple equilibrium points. This can be ascertained through the use of CRNT, which also allows the computation of a candidate set of parameter values. Subsequently, the second stage of the procedure focuses on the analysis of the DA’s of the equilibrium points, using the method devised in [13], since the DA can be regarded as a robustness measure against perturbations that push the system away from its steady state operative condition. The third phase of the procedure consists in the analysis of the robustness of the system’s bistability with respect to parametric uncertainty. Traditionally, this analysis is based on sensitivity and bifurcations analysis; however, these tools are rather limited, due to their inability to take into account multiple simultaneous parameter variations. To overcome these limitations, we have proposed a multi-parametric robustness analysis strategy: by opportunely leveraging a global sensitivity analysis method, and combining it with the information provided by our DA’s analysis technique, we were able to certify the persistence of bistability in the face of multiple variations of the uncertain parameters.
Beyond its specific application for validation of the proposed model of the GAL regulatory network, the overall procedure provides a powerful approach for the analysis and validation of any biochemical network model which is required to robustly reproduce bistable dynamics, underlying persistent memory, molecular switches and cell differentiation phenomena, without recourse to large-scale numerical simulations.
Methods
Chemical reaction network theory
Given a reaction network, the capability of the associated model to exhibit two or more equilibrium points depends on the mathematical form of the reaction rates and on the specific values of the kinetic parameters.
While the characterization of multistability for a generic nonlinear system requires an ad hoc mathematical treatment, in the case of mass-action systems it can be performed through a powerful analytical tool, namely Chemical Reaction Network Theory (CRNT) [23, 24]. CRNT links the structure of a biochemical network, endowed with mass-action kinetics, to the capability of the network to admit multiple positive steady-states. The advantage of CRNT is that it provides an immediate way to analyze the type of dynamical behavior that one can expect from an arbitrarily complex reaction network, just by inspection of the topology of the associated graph. More specifically, CRNT enables us to establish the conditions for the existence, multiplicity and stability of fixed points for the associated ODE system, without even the need to write down the kinetic equations nor to assign values to the kinetic parameters. This point makes CRNT especially suitable for dealing with biomolecular systems, whose parameters are often unknown or subject to significant variability among different individuals.
For a complete description of the CRNT, the interested reader is referred to the original articles [23, 24], or to [28] for an introductory overview of the main results. Despite the complexity of its theoretical foundations, the application of the main CRNT results is straightforward through the use of the CRNT algorithm, which is implemented in the CRNT Toolbox^{1}. Given the reaction network that we want to study, the algorithm establishes whether the associated mass-action dynamical system can admit multiple positive steady states for some values of the kinetic parameters. In the affirmative case, the algorithm also computes a set of values of the kinetic parameters for which the system is multistable.
Analysis of the domains of attraction
The exact computation of the DA’s for a nonlinear system is generally a very hard problem to solve, especially for systems of medium/high order. The problem of computing estimates of the DA’s has been studied for a long time and several methods, based on Lyapunov functions, were originally proposed, e.g. in [29, 30]. More recently, Chesi has devised novel results concerning DA analysis based on Sum of Squares (SOS) representation of polynomials and Semi-Definite Programming (SDP), [31, 32]. Moreover, effective examples of the usefulness of SOS/SDP-based approaches to elucidate the properties of biological systems are provided in [33, 34].
Topcu et al. in [35, 36] have dealt with the topic of estimating a robust DA in the case of uncertain parameters. Compared to the method used in this work, their results can deal with a larger class of systems, namely polynomial dynamical systems; however, they cast the problem in the form of Bilinear Matrix Inequalities (BMI’s), whose solution is much more demanding than LMI’s (used by our approach) and, thus, its practical applicability is limited to systems of low order with few optimization variables..
When dealing with nonlinear quadratic systems, an alternative approach to DA analysis, proposed by Amato and coworkers in [13, 14], can be adopted: this method allows one to check whether an assigned box (or, more in general, a polytope) in the state space belongs to the domain of attraction of a given equilibrium point. Such problem can be cast as a LMI feasibility problem [37], which can be effectively tackled through effective off-the-shelf numerical tools. In what follows we provide a brief overview of the main result used in the present work.
where $A,{F}_{i}\in {\mathbb{R}}^{9\times 9}$. We can now precisely state the problem to be solved. Note that, for the sake of brevity, the statement of the problem refers to the zero equilibrium point, i.e. the origin of the state space. Nevertheless, it is easy to generalize the definition and the entire procedure to non-zero equilibrium points, via a change of state variables, as shown in [13].
Problem 1
Assume that each eigenvalue of matrix A in (5) has strictly negative real part (i.e., the origin is an asymptotically stable equilibrium point); then, given a box $\mathcal{S}$, with the origin of the state space lying in the interior of $\mathcal{S}$, establish whether $\mathcal{S}$ belongs to the DA of the zero equilibrium point.
The following Theorem provides sufficient conditions to solve Problem 1.
Theorem 1
For a fixed γ, conditions (7) constitute a set of LMI’s, which can be easily solved through off-the-shelf efficient numerical tools (e.g., the LMILAB provided in the MATLAB Robust Control Toolbox [38]).
In order to find the largest possible estimate of the DA, Theorem 1 can be applied iteratively, starting from a small initial box ${\mathcal{P}}_{0}$, surrounding the equilibrium point, and then stretching the box at each iteration along the different dimensions of the state space, until conditions (7) become unfeasible.
Local sensitivity analysis
The set of differential equations that constitutes the dynamical system is coupled to the equations of the sensitivity coefficients. This allows computing a numerical solution to the whole set of equations, thus simultaneously obtaining both the state variables and the associated sensitivity coefficients.
It is worth recalling that traditional sensitivity analysis methods are only valid locally with respect to a particular point in the model’s parameter space, i.e., in the neighborhood of a certain parameter set. Another significant limitation consists in their capability to consider the sensitivity of the model with respect to the variations of just one parameter at a time; indeed, a model might display low sensitivity to single parameter variations, while being extremely sensitive to simultaneous multiple parameter changes.
Bifurcations analysis
A bifurcation occurs at values of p such that small changes of the parameters can dramatically alter the number or types of attractors of system (9) [39].
The changes in the map of equilibrium points can be effectively visualized by using a ’bifurcation diagram’, in which the steady state value of one state variable is plotted as a function of a bifurcation parameter. The calculation of bifurcations can be performed through continuation softwares, like the MatCont package [40], which we have used to perform the analysis illustrated in Figures 5 and 6.
Bifurcation diagrams are powerful tools in order to investigate the robustness of nonlinear biomodels in the face of parametric uncertainty. However, it is necessary to take into account that a) analytical solutions for bifurcations are only available for low-dimensional models, and b) that bifurcation diagrams are practically applicable only to study the effect of one or two parameters variation at a time.
Global sensitivity analysis via infeasibility certificates
In view of the limitations of the approaches presented in the previous two sections, while it is possible to employ them to obtain preliminary information on the parametric robustness of a given model, particular care must be taken in drawing any conclusions about global properties of the system under investigation.
To overcome these limitations, we have exploited a global sensitivity analysis technique for biochemical networks proposed in [20]. Given the admissible parameters variation box, the approach proposed by Waldherr et al. allows one to compute an outer approximation, $\stackrel{~}{\mathcal{S}}$, of the region of the state space that contains all the equilibrium points, denoted by ${\mathcal{X}}_{e}$. The problem can be formalized as follows.
Problem 2
Note that, apart from trivial cases, the calculation of an analytical form of ${\mathcal{X}}_{e}$ is practically impossible. Moreover, computational brute-force approaches are applicable only to very low-order systems. Monte-Carlo techniques can be applied in the other cases, although they may require a large computational effort and guarantee only probabilistic results.
where p_{j,min}, p_{j,max}, j=1,…,m, define the admissible parameter box ${\mathcal{\mathcal{B}}}_{p}$, and x_{i,min}, x_{i,max}, i=1,…,n, are the extremal values of the box $\stackrel{~}{\mathcal{S}}$ of the state space to be tested as a candidate solution to Problem 2.
The optimization problem (11) is not easy to deal with from the computational point of view. However, it can be tackled by solving its dual version, that is the problem of computing regions of the state space that are guaranteed not to contain any steady state for any parameter value in ${\mathcal{\mathcal{B}}}_{p}$. The latter can be relaxed to become a SDP problem [41], and solved by means of computationally efficient convex optimization tools. For a detailed description of this procedure, the reader is referred to [20, 21]. In this works, the computation of a solution to problem 2 constitutes the core of an iterative procedure, implemented by the bioSDP algorithm: starting from an initial large region of the state space, the algorithm tries to compute one or more partitions containing ${\mathcal{X}}_{e}$. The procedure is very effective for low-order systems (n≤3), since in this case a bisection algorithm can be used for the partitioning. For systems of higher order, a box shrinkage procedure is employed, which can only return one partition $\stackrel{~}{\mathcal{S}}$ and, therefore, is not useful for analyzing the persistence of bistability.
In order to solve this problem, we have devised an alternative strategy, which combines the results of the DA’s analysis with the bioSDP algorithm and has proven to be effective in the analysis of our case study. The details of this approach have been already reported in the results Section ‘Local and global analysis confirm robustness of the bistable equilibria to parametric uncertainty’.
Endnotes
^{1}The CRNT algorithm is implemented in the CRNT toolbox, which is freely available at http://www.che.eng.ohio-state.edu/~feinberg/crnt/
Notes
Declarations
Acknowledgements
DGB acknowledges support from EPSRC Research Grant EP/I017445/1.
Authors’ Affiliations
References
- Bhat PJ: Galactose Regulon of Yeast: From Genetics to Systems Biology. 2008, Berlin: Spring-VerlagView ArticleGoogle Scholar
- Acar M, Becskei A, van Oudenaarden A: Enhancement of cellular memory by reducing stochastic transitions. Nature. 2005, 435 (7039): 228-232. 10.1038/nature03524.PubMedView ArticleGoogle Scholar
- de Atauri P, Orrell D, Ramsey S, Bolouri H: Evolution of ‘design’ principles in biochemical networks. Syst Biology, IEE Proc. 2004, 1: 28-40. 10.1049/sb:20045013.View ArticleGoogle Scholar
- Yang R, Lenaghan SC, Wikswo JP, Zhang M: External control of the GAL network in S. cerevisiae: A view from control theory. PLoS ONE. 2011, 6 (4): e19353-10.1371/journal.pone.0019353.PubMedPubMed CentralView ArticleGoogle Scholar
- Smidtas S, Schäcter V, Képès F: The adaptive filter of the yeast galactose pathway. J Theor Biol. 2006, 242 (2): 372-381. 10.1016/j.jtbi.2006.03.005.PubMedView ArticleGoogle Scholar
- Kulkarni VV, Kareenhalli V, Malakar PP, Pao LY, Safonov MG, Viswanathan A Ganesh: Stability analysis of the GAL regulatory network in Saccharomyces cerevisiae and Kluyveromyces lactis. BMC Bioinformatics. 2010, 11 (Suppl 1): S43-10.1186/1471-2105-11-S1-S43.PubMedPubMed CentralView ArticleGoogle Scholar
- Cosentino C, Salerno L, Passanti A, Merola A, Bates DG, Amato F: Structural bistability of the GAL regulatory network and characterization of its domains of attraction. J Comput Biol. 2012, 19 (2): 148-162. 10.1089/cmb.2011.0251.PubMedView ArticleGoogle Scholar
- Mélykúti B, August E, Papachristodoulou A, El-Samad H: Discriminating between rival biochemical network models: three approaches to optimal experiment design. BMC Syst Biol. 2010, 4: 38-10.1186/1752-0509-4-38.PubMedPubMed CentralView ArticleGoogle Scholar
- Wagner A: Robustness and Evolvability in Living Systems. 2005, Princeton: Princeton University PressGoogle Scholar
- Morohashi M, Winn AE, Borisuk MT, Bolouri H, Doyle J, Kitano H: Robustness as a measure of plausibility in models of biochemical networks. J Theor Biol. 2002, 216: 19-30. 10.1006/jtbi.2002.2537.PubMedView ArticleGoogle Scholar
- Bates DG, Cosentino C: Validation and invalidation of systems biology models using robustness analysis. IET Syst Biol. 2011, 5 (4): 229-244. 10.1049/iet-syb.2010.0072.PubMedView ArticleGoogle Scholar
- Anderson J, Papachristodoulou A: On validation and invalidation of biological models. BMC Bioinformatics. 2009, 10: 132-10.1186/1471-2105-10-132.PubMedPubMed CentralView ArticleGoogle Scholar
- Amato F, Cosentino C, Merola A: On the region of attraction of nonlinear quadratic systems. Automatica. 2007, 43 (12): 2119-2123. 10.1016/j.automatica.2007.03.022.View ArticleGoogle Scholar
- Amato F, Calabrese F, Cosentino C, Merola A: Stability analysis of nonlinear quadratic systems via polyhedral lyapunov functions. Automatica. 2011, 47 (3): 614-617. 10.1016/j.automatica.2010.12.005.View ArticleGoogle Scholar
- Stelling J, Sauer U, Szallasi Z, Doyle III FJ, Doyle J: Robustness of cellular functions. Cell. 2004, 118 (6): 675-685. 10.1016/j.cell.2004.09.008.PubMedView ArticleGoogle Scholar
- Eißing T, Waldherr S, Allgöwer F, Scheurich P, Bullinger E: Steady state and (bi-) stability evaluation of simple protease signalling networks. BioSystems. 2007, 90 (3): 591-601. 10.1016/j.biosystems.2007.01.003.PubMedView ArticleGoogle Scholar
- Kim J, Bates DG, Postlethwaite I, Ma L, Iglesias PA: Robustness analysis of biochemical network models. IEE Proc: Syst Biol. 2006, 153 (3): 96-104. 10.1049/ip-syb:20050024.Google Scholar
- Blüthgen N, Herzel H: How robust are switches in intracellular signaling cascades?. J Theor Biol. 2003, 225 (3): 293-300. 10.1016/S0022-5193(03)00247-9.PubMedView ArticleGoogle Scholar
- Eißing T, Allgöwer F, Bullinger E: Robustness properties of apoptosis models with respect to parameter variations and intrinsic noise. Syst Biol, IEE Proc. 2005, 152 (4): 221-228. 10.1049/ip-syb:20050046.View ArticleGoogle Scholar
- Waldherr S, Findeisen R, Allgower F: Global sensitivity analysis of biochemical reaction networks via semidefinite programming. Proc. of 17th World Congress, The International Federation of Automatic Control: 6-11 July 2008; Seoul, Korea. 2008, 9701-9706.Google Scholar
- Hasenauer J, Rumschinski P, Waldherr S, Borchers S, Allgöwer F, Findeisen R: Guaranteed steady state bounds for uncertain (bio-)chemical processes using infeasibility certificates. J Process Control. 2010, 20 (9): 1076-1083. 10.1016/j.jprocont.2010.06.004.View ArticleGoogle Scholar
- Angeli D, Ferrell JEJ, Sontag ED: Detection of multistability, bifurcations, and hysteresis in a large class of biological positive-feedback systems. Proc Nat Acad Sci USA. 2004, 101 (7): 1822-1827. 10.1073/pnas.0308265100.PubMedPubMed CentralView ArticleGoogle Scholar
- Feinberg M: Chemical reaction network structure and the stability of complex isothermal reactors-I. The deficiency zero and deficiency one theorems. Chem Eng Sci. 1987, 42 (10): 2229-2268. 10.1016/0009-2509(87)80099-4.View ArticleGoogle Scholar
- Feinberg M: Chemical reaction network structure and the stability of complex isothermal reactors-II. Multiple steady states for networks of deficiency one. Chem Eng Sci. 1988, 43: 1-25. 10.1016/0009-2509(88)87122-7.View ArticleGoogle Scholar
- Ruhela A, Verma M, Edwards JS, Bhat PJ, Bhartiya S, Venkatesh KV: Autoregulation of regulatory proteins is key for dynamic operation of GAL switch in Saccharomyces cerevisiæ. FEBS Lett. 2004, 576 (12): 119-126.PubMedView ArticleGoogle Scholar
- Pannala VR, Hazarika SJ, Bhat PJ, Bhartiya S, Venkatesh KV: Growth-related model of the GAL system in saccharomyces cerevisiae predicts behaviour of several mutant strains. IET Syst Biol. 2012, 6 (2): 44-53. 10.1049/iet-syb.2010.0060.PubMedView ArticleGoogle Scholar
- Ingalls BP, Sauro HM: Sensitivity analysis of stoichiometric networks: An extension of metabolic control analysis to non-steady state trajectories. J Theor Biol. 2003, 222: 23-36. 10.1016/S0022-5193(03)00011-0.PubMedView ArticleGoogle Scholar
- Cosentino C, Bates DG: Feedback Control in Systems Biology. 2011, Boca Raton: CRC PressGoogle Scholar
- Genesio R, Tartaglia M, Vicino A: On the estimation of asymptotic stability regions: State of the art and new proposals. IEEE Trans Autom Control. 1985, 30 (8): 747-755. 10.1109/TAC.1985.1104057.View ArticleGoogle Scholar
- Vannelli A, Vidyasagar M: Maximal Lyapunov functions and domains of attraction for autonomous nonlinear systems. Automatica. 1985, 21: 69-80. 10.1016/0005-1098(85)90099-8.View ArticleGoogle Scholar
- Chesi G: Estimating the domain of attraction for uncertain polynomial systems. Automatica. 2004, 40 (11): 1981-1986. 10.1016/j.automatica.2004.06.014.View ArticleGoogle Scholar
- Chesi G: Domain of Attraction. Analysis and control Via SOS Programming. 2011, London: Springer-VerlagView ArticleGoogle Scholar
- El-Samad H, Prajna S, Papachristodoulou A, Doyle J, Khammash M: Advanced methods and algorithms for biological networks analysis. Proc IEEE. 2006, 94 (4): 832-853.View ArticleGoogle Scholar
- August E, Papachristodoulou A: A new computational tool for establishing model parameter identifiability. J Comput Biol. 2009, 16 (6): 875-885. 10.1089/cmb.2008.0211.PubMedView ArticleGoogle Scholar
- Topcu U, Packard A: Local stability analysis for uncertain nonlinear systems. IEEE Trans Autom Control. 2009, 54: 1042-1047.View ArticleGoogle Scholar
- Topcu U, Packard A, Seiler P, Balas G: Robust region-of-attraction estimation. IEEE Trans Autom Control. 2010, 55: 137-141.View ArticleGoogle Scholar
- Boyd S, El Ghaoui L, Feron E, Balakrishnan V: Linear Matrix Inequalities in System and Control Theory, Volume 15 of Studies in Applied Mathematics. 1994, Philadelphia: SIAMView ArticleGoogle Scholar
- MathWorks T: MATLAB Robust Control Toolbox. 2011, Natick: MathWorks Inc.Google Scholar
- Kuznetsov Y: Elements of Applied Bifurcation Theory. No. v. 112 in Applied Mathematical Sciences. 2004, SpringerView ArticleGoogle Scholar
- Dhooge A, Govaerts W, Kuznetsov YA: MATCONT: A MATLAB package for numerical bifurcation analysis of ODEs. ACM Trans Math Software. 2003, 29 (2): 141-164. 10.1145/779359.779362.View ArticleGoogle Scholar
- Vandenberghe L, Boyd S: Semidefinite programming. SIAM Rev. 1996, 38: 49-95. 10.1137/1038003.View ArticleGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.