Structural Bistability of the GAL Regulatory Network and Characterization of its Domains of Attraction

Abstract

Bistability is a system-level property, exploited by many biomolecular interaction networks as a key mechanism to accomplish different cellular functions (e.g., differentiation, cell cycle, switch-like response to external stimuli). Bistability has also been experimentally found to occur in the regulatory network of the galactose metabolic pathway in the model organism Saccharomyces cerevisiae. In this yeast, bistability generates a persistent memory of the type of carbon source available in the extracellular medium: under the same experimental conditions, cells previously grown with different nutrients generate different responses and get stably locked into two distinct steady states. The molecular interactions of the GAL regulatory network have been thoroughly dissected through wet-lab experiments; thus, this system provides a formidable benchmark to our ability to characterize and reproduce in silico the behavior of bistable biological systems. To this aim, a number of models have been proposed in the literature; however, we found that they are not able to replicate the persistent memory behavior observed in (Acar et al., 2005). The present study proposes a novel model of the GAL regulatory network, which, in addition to reproducing in silico the experimental findings, can be formally analyzed for structural multistability of the network, using chemical reaction network theory (CRNT), and allows the characterization of the domains of attraction (DA). This work provides further insights into the GAL system and proposes an easily generalizable approach to the study of bistability-associated behaviors in biological systems.

1. Introduction

Aprimary goal of biologists in the postgenomic era is to understand how characteristic cellular phenotypes arise from the interactions between the basic elements of life, namely genes, proteins, and metabolites. By now, it is utterly clear that the properties of complex biomolecular interaction networks cannot be solely unveiled by studying the single components, but rather they have to be investigated at the system level, exploiting the methodological tools of bioinformatics, and systems and control theory (Cosentino et al., 2011). Among these properties, an especially important one is the capability of several signaling networks to exhibit a switch-like biochemical response, which enables a cell to toggle between two discrete, alternative, equilibrium conditions following an external stimulus and to remain locked into the new steady state even when the stimulus is removed. Well-known examples of this class of biological systems are the mitogen-activated protein kinase (MAPK) cascades in animal cells (Bagowski et al., 2001; Ferrell et al., 1998) and the cell cycle regulatory network in Xenopus and Saccharomyces cerevisiae (Sha et al., 2003; Pannala et al., 2010). In S. cerevisiae, one of the most extensively studied genetic switch is that implemented by the regulatory network of the galactose metabolic pathway: S. cerevisiae normally uses glucose as the carbon and energy source. However, in the absence of glucose, it can metabolize galactose as an alternative source, through Leloir metabolic pathway (Bhat, 2008). Such pathway is controlled by a set of enzymes, commonly referred to as the GAL system, whose molecular interactions have been thoroughly dissected in the last years through the classical reductionist molecular biology experimental approach. A number of mathematical models of the GAL system, in conjunction with guided experiments, have investigated system-level properties such as ultrasensitivity, memory, noise attenuation, rapid response, and sensitive response (Acar et al., 2005). In particular, the capability of the GAL system to exhibit persistent memory of the past exposure to different carbon sources, has been experimentally found in Acar et al. (2005), thus demonstrating the bistability of the underlying regulatory network. However, notwithstanding this experimental evidence, the mathematical models proposed so far are not able to reproduce this particular property. This fact has been recently pointed out in Kulkarni et al. (2010), where the authors test the dynamical nonlinear quadratic model proposed in Smidtas et al. (2006) against the aforementioned experimental findings and conclude that it does not exhibit bistability. Kulkarni et al. (2010) extend the Smidtas' model and derive some conclusions about cellular memory by applying a method devised by our group (Amato et al., 2007) to analyze the domains of attraction (DA) of an equilibrium point in nonlinear quadratic systems. On the hand, however, Kulkarni et al. (2010) do not provide a theoretical analysis nor numerical simulations to prove whether the proposed model actually exhibits a bistable behavior.

In this work, we propose a novel model that maintains the quadratic form of the most recently proposed models, derived by applying mass-action kinetics. Indeed, such class of dynamical systems can be easily tested for structural multistability, applying chemical reaction network theory (CRNT) (Feinberg, 1987, 1988). CRNT is a powerful theoretical tool to establish whether an assigned reaction network, endowed with mass action kinetics, can exhibit multiple steady states, independently of the particular choice of the kinetic parameter values. Applying CRNT, we were able to conclude that both Smidtas' and Kulkarni's models cannot exhibit bistability (see Appendix); instead, our model verifies the CRNT conditions for bistability, and in fact, it can reproduce the persistent memory behavior observed in Acar et al. (2005). Subsequently, we have also performed an analysis of the DA, applying our method devised in Amato et al. (2007), to characterize the persistence of cellular memory in correspondence of both equilibrium points.

This article is structured as follows: Section 2 provides a brief overview of the biochemical network of the GAL regulatory system. In Section 3, we describe the new mathematical model devised to replicate the persistent memory property of the GAL system. The theoretical analysis tools used in the presented study are discussed in Section 4, and their application is illustrated in the Section 5, along with the results obtained through numerical simulations. Finally, in Section 6, we discuss the results and provide some concluding remarks. Some additional materials, including detailed models and parameter values, are given in the Appendix.

2. Gal Regulatory Network in S. Cerevisiæ

S. cerevisiae exhibits the capability to metabolize galactose by activating the coordinated transcription of both structural and regulatory genes. This set of genes enables yeast cells to use galactose as an alternative energy and carbon source to glucose, when the latter is absent (Bhat, 2008; Lohr, et al., 1995). The entire pathway, which leads to the consumption of galactose comprises two main subsystems: (i) a metabolic pathway, from galactose cellular internalization to the production of glucose, and (ii) a regulatory pathway, which produces all the proteins that activate the transcription of galactose metabolizing enzymes (Pannala et al., 2010).

The regulation pathway (Fig. 1) is made up of the following autoregulated factors: a transcriptional activator protein Gal4p, a repressor protein Gal80p, and a signal transducer protein Gal3p (Timson et al., 2002). Gal4p activates transcription of genes GAL1, GAL2, GAL7, and GAL10, but is unable to activate the transcription in the absence of galactose. The protein encoded by gene GAL2 acts as a mediator of galactose transport into the cell. In the absence of external galactose, Gal80p binds to the activation domain of Gal4p, thus inhibiting the expression of the GAL genes (Melcher et al., 2001). In the presence of galactose, instead, the inducer Gal3p is activated to form the complex Gal3p*:Gal80p, which promotes the shuttling of Gal80p from the nucleus to the cytoplasm. This decreases the concentration of Gal80p, reducing the fraction of Gal80p-bound Gal4p in the nucleus (Peng et al., 2002). Thus, galactose relieves the inactivation of Gal4p and promotes transcription of the GAL genes. The presence of glucose (even when galactose is available) inhibits the synthesis of Gal4p through Mig1p-mediated repression of GAL genes (Klein et al., 1998).

FIG. 1.

Schematic diagram of the GAL regulatory network in Saccharomyces cerevisiæ. External galactose is internalized into the cell, then it shuttles between the cytoplasm and the nucleus. The signal transducer Gal3p in galactose bound-stage is highlighted to differentiate it from unbound form. The pointed arrows indicate activation, whereas the blunt arrows indicate inhibition.

The persistence of cellular memory exhibited by the GAL regulatory network is a system-level property, resulting from the interactions of several species, forming multiple nested loops. Two coupled positive feedback loops include the galactose permease Gal2p and the signaling protein Gal3p, respectively, and a negative feedback loop involves the inhibitor Gal80p. In Angeli et al. (2004), it has been noted that the existence of a positive-feedback loop, or a mutually inhibitory, double-negative-feedback loop, is a necessary condition for the occurrence of multistability. History-dependent experiments on galactose regulatory gene activity revealed the capability of this molecular circuitry to toggle between two discrete, alternative stable steady states, based on the galactose concentration. Exposing two cell cultures, one grown in the absence and the other in the presence (2%) of galactose, to extreme (i.e., either high or low) galactose concentrations, they reach steady states, independently of the past history (absence of memory). For intermediate concentrations, instead, the two cultures show significantly different expression of the GAL genes (persistent memory). Signal transducer Gal3p was identified as the most responsible in the generation of two stable expression states depending on the previous galactose consumption states. The other positive feedback loop mediated by the Gal2p is only able to increase the difference among the two levels of expression. A higher level of inhibitor Gal80p results in destabilization of persistent memory (Acar et al., 2005).

3. Mathematical Model of the Gal Regulatory Network in S. Cerevisiæ

To evaluate the dynamics of the GAL regulatory pathway from external galactose uptake to its downstream elements, several models have already been proposed in the literature, including both simplified reduced-order models (Smidtas et al., 2006; Kulkarni et al., 2010) and more detailed ones (de Atauri et al., 2004). We have tested these models against the experimental findings of Acar et al. (2005) to ascertain whether they can reproduce the persistent memory behavior; however, both numerical simulations and theoretical analysis performed on such models have provided a negative answer (see Appendix).

On the basis of these findings, we have devised a novel mathematical model, based on the schematic diagram of interactions reported in Figure 1. The complete set of reactions is reported in Appendix, and the state variables (concentrations of each species) of the dynamical model are given in Table 1. The main goal of our work is to achieve a simple model, derived by applying mass-action kinetics, which can be characterized with respect to bistability by means of a theoretical analysis tool, namely CRNT. Therefore, the model does not aim to achieve a detailed description of the biochemical reaction, which would require much more complex nonlinear kinetics, with the drawback of rendering the stability analysis a daunting task. Additionally, the formulation of the model in the framework of dynamical quadratic systems enables us to quantitatively characterize the DA of the equilibrium points, as will be illustrated in the following.

Table 1.

Description of the State Variables of System

State Variable	Description
G₃	Gal3p protein concentration
G_i	Internalized galactose concentration
G_3a	Active Gal3p protein concentration
G₄	Gal4p protein concentration
G₈₀	Gal80p protein concentration
G_4,80	Gal4p:Gal80p complex concentration
G_3a,80	Active Gal3:Gal80p complex concentration
G_ex	Extracellular galactose concentration
G₂	Gal2p protein concentration

Note that, in order to focus our study on the bistable behavior, the following simplifying assumptions have been made: (i) the GAL gene inhibition mechanisms, which are activated in the presence of glucose, are neglected, since we always assume a zero concentration of glucose; (ii) the galactose uptake into the cell occurs only through the protein Gal2 (i.e., we neglect the intrinsic Gal2-indipendent transport); (iii) the nucleoplasmic shuttling is not modeled, and the cytoplasm and nucleus are treated as a single cellular compartment; (iv) binding/unbinding of Gal4 to/from DNA is neglected; (v) no dimerization of proteins is modeled; and (vi) only one binding site for Gal4 is modeled for each species.

The pathway starts from the external galactose uptake in the cell membrane and is composed of 16 reaction steps. External galactose is internalized into the yeast cell via the action of Gal2, as described by the reaction \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*} G_e + G_2 \rightarrow G_i + G_2 + G_e \tag{1} \end{align*} \end{document}

For each internalized molecule, G_i, a new molecule of external galactose, G_e, is added. This is a technical assumption, which is required to model a medium that is kept at constant galactose concentration, as in the experiments described in Acar et al. (2005). The transcription and translation of Gal2, Gal3, and Gal80, promoted by Gal4p, is modeled as \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*} G_4 \rightarrow G_y + G_4 , \ y = 2 , 3 , 80 \tag{2} \end{align*} \end{document}

Activation of Gal3p by internalized galactose is described as \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*} G_3 + G_i \rightarrow G_{3a} \tag{3} \end{align*} \end{document}

The two steps of inactivation of the inhibitor and subsequent release of Gal4p are modeled as a single reaction: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*} G_{4 , 80} + G_{3a} \rightarrow G_{3a , 80} + G_4 \tag{4} \end{align*} \end{document}

Finally, each species is subject to degradation, which is modeled by first-order kinetics (i.e., each species is transformed into a null product): \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $$G_x \rightarrow \emptyset$$ \end{document} .

Assuming mass action kinetics for each reaction rate, we obtain a set of nonlinear ordinary differential equations (ODEs), containing constant, linear, and quadratic terms in the state variables. The model, which describes the concentration changes over time of the considered species, can be conveniently written in matrix form as \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*} \dot{x} ( t ) = Ax ( t ) + B ( x ) + Nu ( t ) , \tag{5} \end{align*} \end{document}

where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*} x ( t ) = [ G_3 ( t ) \ G_i ( t ) \ G_{3a} ( t ) \ G_4 ( t ) \ G_{80} ( t ) \ G_{4 , 80} ( t ) \ G_{3a , 80} ( t ) \ G_2 ( t ) ] ^T , \tag{6} \end{align*} \end{document}

for t ≥ 0 and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*} B ( x ) = \left( \begin{matrix}x^TB_1 \\ \vdots \\ x^TB_n\end{matrix} \right) x , \tag{7} \end{align*} \end{document}

with A, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $$B_i \in {\mathbb R}^{n \times n} , i = 1 , \ldots 8.$$ \end{document}

Looking at the topologies of previous models (Fig. 2), we can spot some basic differences with respect to the one proposed in this article. In particular, in Smidtas' detailed model (Fig. 2a), Gal3p, Gal80p, and galactose bind together, decreasing the binding of Gal80p to Gal4p. Subsequently, the decreasing of Gal80p in the nucleus activates Gal4p, triggering transcription of GAL3 and GAL80 genes and closing the feedback loops, as newly produced Gal3p and Gal80p switch the equilibrium back towards Gal4p inactivation (Gal4/80p). In the Smidtas' reduced model (Fig. 2b), an oversimplification of these feedback loops is adopted, preserving the qualitative dynamics of the detailed model. However, the two feedback loops involving GAL3 (positive) and GAL80 (negative) are lumped together in a single negative feedback loop. In the absence of galactose, Gal3/80p binds the Receptor (R), thus inhibiting its transcriptional activity of Gal4p. The presence of galactose causes the release of Gal4p and consequential activation of GAL3 and GAL80 transcription. Thus, newly produced Gal3/80p closes the feedback towards Gal4p inactivation. Starting from Smidtas' simplified model, Kulkarni et al. (2010) have proposed an extended model, by including further molecular reactions and, hence, more state variables (Fig. 2c). However, their model inherits the simplification regarding the single negative feedback loop of Gal3/80p. Our model (Fig. 2d) includes more state variables, since we have also taken into account the Gal2p-mediated uptake of galactose and the formation of Gal3/80p and Gal4/80p complexes from Gal3p, Gal4p, and Gal80p as reactants, which in the other models are involved in bound stage only. Structurally, our system is more adherent to the circuit depicted in Acar et al. (2005), consisting of two positive and one negative feedback loops, which affect the uptake of galactose, the nucleoplasmic shuttling of regulator proteins, and the transcription of GAL genes. Therefore, we can conclude that the presence of these distinct positive and negative feedback loops is the key ingredient that determines the occurrence of bistability in the GAL system.

FIG. 2.

Graphical representation of the GAL core regulatory pathway. Boxes, species; circles, reactions. Arrows point to reaction products; reagents are connected to the reaction by lines (without arrows). A dashed line indicates that the linked species is both a reagent and a product. (a) Smidtas' detailed model. Gal, galactose; bound receptor (BR), the Gal3/80p/galactose complex; Gal4_80p and Gal3_80p denote, respectively, the Gal4p and Gal80p and the Gal3p and Gal80p proteins in bound stage. (b) Smidtas' reduced model. Receptor (R), Gal4/3/80p complex, in which Gal4p is sequestered and cannot promote its transcriptional activity, whereas G4 denotes the transcriptional activator protein Gal4p. (c) Kulkarni's model, in which G4a indicates the activated Gal4p protein and G4 denotes the Gal4/3/80p complex. (d) Our model. It is possible to note that more state variables and molecular reactions are considered. Here, Gx (where x = 2, 3, 4, and 80) denotes the Galxp protein; Gex and Gic, the extra- and the intra-cellular galactose; G3att, the signaling protein Gal3p in active form; G80G3att and G4G80, the complexes formed by Gal80p and Gal3p (in active form), and by Gal4p and Gal80p proteins, respectively.

In the next section, we will introduce two theoretical tools to investigate the properties of a dynamical nonlinear quadratic system in the general form (eq. 5).

4. Structural Multistability and da of Nonlinear Quadratic Systems

4.1. Multistability analysis: CRNT

We recall that the main goal of our study is to ascertain whether the modeled system exhibits bistability (or, more generally, multistability). Note that, given a certain reaction network, the capability of the system to exhibit multiple equilibrium points depends on the associated reaction rates and on the values of the kinetic parameters. Therefore, it is in general not possible to draw any conclusion without restricting the scope of the analysis to a more specific class of systems. In our case, having assumed mass-action kinetics, the structure of the nonlinear ODE system is well defined and corresponds to the class of dynamical quadratic systems with nonnegative-valued parameters. Even for models of this class, however, the existence of multiple steady states is in general dependent on the particular choice of the parameters.

Fortunately, under certain conditions, the characterization of multistability of mass-action systems can be performed through a powerful analysis tool, CRNT. The advantage of CRNT is that it provides a straightforward way to analyze the type of dynamical behavior that one can expect from an arbitrarily complex network of chemical reactions, just by inspection of the topology of the associated graph. More specifically, CRNT enables us to establish whether an assigned reaction network can exhibit one or multiple equilibrium points, without even the need to write down the kinetic equations and 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 the sake of brevity, it is not possible to provide here the nuts and bolts of CRNT, which is quite involved and requires a considerable number of preliminary definitions. The interested reader is referred to the original articles (Feinberg, 1987, 1988) for a detailed treatment or to Cosentino et al. (2011) for an introductory overview of the main results. Despite the complexity of its theoretical foundations, a great advantage of CRNT is that it requires very little background for practical application; indeed, the CRNT analysis conditions have been coded in the CRNT Toolbox (http://www.chbmeng.ohio-state.edu/∼feinberg/crnt), which is freely available. Using this toolbox, it is sufficient to fill in the network's reactions and run the algorithm to get a comprehensive report elucidating all the properties of the network of interest. In particular, the CRNT algorithm indicates either that the associated dynamical system admits multiple positive steady states for some values of the kinetic parameters, or else that no such combination of the parameter values exist. In the affirmative case, the algorithm will also provide a set of values of the kinetic parameters for which the system is multistable.

Finally, it is worth remarking a key point of this approach to the analysis of multistability of chemical reaction systems: CRNT enables us to characterize multistability as a structural property of the given reaction network, regardless of any specific choice of the model parameters and of the complexity of the network. Instead, classical methods for the study of equilibrium points (Murray, 2001) can provide parametric multistability analysis only for low-order systems, since they are based on the analytical computation of the equilibrium points.

In the Results section, we will apply CRNT to invalidate the models proposed in Smidtas et al. (2006) and Kulkarni et al. (2010) and to validate the model proposed in Section 3.

4.2. Estimation of the DA

Determining the number of equilibrium points and their stability is a fundamental step in the analysis of biochemical reaction models. However, a complete characterization of the system's behavior requires the computation of the DA of each asymptotically stable equilibrium point. The DA of an asymptotically stable equilibrium point, x_e, is a subset of the state space, Ω, such that every state trajectory starting in Ω converges to x_e.

DAs are of paramount importance in biological systems; indeed, when a system that is operating in the neighborhood of an equilibrium point leaves the boundaries of its region of attraction, it is usually abruptly led to a new operating condition (corresponding to another equilibrium point or to oscillations). This phenomenon is at the basis of many on-off regulatory mechanisms of biological functions and periodic behaviors at the molecular, cellular, or population level and is closely interweaved with an important issue in systems biology, namely the relation between cellular decision making and biological noise. Cellular decision making is the process by which cells can assume different and heritable fates modifying their genome or the environment in which they reside (Balazsi et al., 2011). Some examples of regulatory molecular switch concern virus, bacteria, yeast, lower metazoans, and mammals. In lambda phage, a gene regulatory network controls the lysis/lysogeny decision in which two repressors (CI and Cro) operate as a mutual toggle switch (lambda switch) with bistable dynamics (Gardner et al., 2000). Another example is represented by the Nanog-Oct4 gene regulatory network, which controls embryonic stem cell differentiation. In the latter system, regulatory interactions, involving positive and negative feedback, generate bistability from noise-induced transitions (Chickarmane et al., 2006).

Typical information that one can derive from the knowledge of the DA is the extent of the perturbations that the system can tolerate without changing its operative condition. In other words, the DA is an indicator of the capability of the system to adapt to changing environmental conditions, which is recognized in biology as an essential feature for every living organism.

Unfortunately, the exact computation of the DA for a generic nonlinear system is a very difficult problem to solve, especially for systems of nontrivial dimension. A more affordable problem is that of determining an estimate of (and possibly close to) the actual DA. Similarly to what has been discussed for the multistability analysis, the problem of DA estimation is greatly simplified when dealing with the class of nonlinear quadratic systems. In particular, below we recall the main result of Amato et al. (2007), which allows to easily and efficiently solve the aforementioned problem. Since the adopted method yields a polytopic estimate of the DA, let us first recall the following definition:

Definition 1

A polytope (or box) \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $${ \mathscr P} \subset {\mathbb R}^n$$ \end{document} can be described as follows: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*} & {\mathscr P} = { \rm conv} \{ x_{ ( i ) } , x_{ ( 2 ) } , \ldots , x_{ ( p ) } \} \tag{ 8{\rm a} } \\ & = \{ x \in {\mathbb R}^n :{a_k}^T x \leq 1 , k = 1 , 2 , \ldots , q \} \tag{8{\rm b} } \end{align*} \end{document}

where p and q are suitable integers numbers, x_(i) denotes the ith vertex of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $${ \mathscr P}$$ \end{document} , and conv[·] indicates the operation of taking the convex hull of the argument.

Next, we state the problem that is solved by the proposed method. Note that, for the sake of simplicity, the statement of the problem refers to the equilibrium point in the origin. However, it is easy to generalize the definition and the procedure to equilibrium points different than zero, via a change of state variables, as shown in Amato et al. (2007).

Problem 2

Assume that every eigenvalue of A in (5) has strictly negative real part; then, given a polytope \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $${ \mathscr P}$$ \end{document} , with the origin of the state space lying in the interior of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $${ \mathscr P}$$ \end{document} , establish whether \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $${ \mathscr P}$$ \end{document} belongs to the DA of system (5).

We can now recall the following sufficient conditions under which Problem 1 admits a solution (Amato et al., 2007).

Theorem 3

The polytope \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $${ \mathscr P}$$ \end{document} defined in (8) belongs to the DA of system (5) if there exist scalars γ,c and a symmetric matrix P such that \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*} & 0 < \gamma < 1 , \tag{ 9{\rm a} } \\ & c > 0 , \tag{9{\rm b} } \\ & P > 0 , \tag{9{\rm c} } \\ & \bigg ( \begin{matrix}1 & \gamma {a_k}^T \\ \gamma a_k & P / c\end{matrix} \bigg ) \geq 0 , k = 1 , 2 , \ldots , 2n \tag{9{\rm d} } \\ & {x_{ ( i ) }}^T ( P / c ) x_{ ( i ) } \leq 1 , \quad i = 1 , 2 , \ldots , 2^n \tag{9{\rm e} } \\ & \gamma ( A^T P + PA ) + \left( \begin{matrix}{x_{ ( i ) }}^T B_1 \\ {x_{ ( i ) }}^T B_2 \\ \vdots \\ {x_{ ( i ) }}^T B_n\end{matrix} \right) P \\ & + { ( B_1}^T x_{ ( i ) } \ {B_2}^T x_{ ( i ) } \cdots {B_n}^T x_{ ( i ) } ) P < 0 , i = 1 , 2 , \ldots 2^n. \tag{9{\rm f} } \end{align*} \end{document}

The inequality conditions given in Theorem 3 can be easily solved through off-the-shelf convex optimization algorithms, for example, using the LMI solvers provided in the MATLAB Robust Control Toolbox (MathWorks, 2011).

In the Results section, after proving bistability through CRNT, we will apply the above result to characterize the DA of the equilibrium points of the GAL system model devised in Section 3.

5. Results

5.1. Structural bistability of the GAL regulatory network

Given the reaction diagram for the GAL regulatory network devised in Section 3 and assuming mass-action kinetics, we have applied the CRNT algorithm to establish whether this regulatory network can admit multiple equilibrium points. Recall that CRNT enables us to answer this question regardless of the kinetic parameters. Once the network is defined in the CRNT Toolbox, we applied the Higher Deficiency Analysis (Ellison, 1998) and Mass Injectivity Analysis (Craciun et al., 2005) routines. These analyses confirm that the network can admit two steady states and provide us with a set of kinetic parameter values (Table 2) for which the following steady states, x_e₁, x_e₂, exist: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*} & x_{e_1} = [ 63.86 \ 63.86 \ 1.0 \ 1.0 \ 1.12 \ 8.64 \ 21.91 \ 1.0 ] ^T \ [ \mu M ] , \\ & x_{e_2} = [ 126.01 \ 126.01 \ 6.67 \ 2.25 \ 1.12 \ 7.39 \ 62.86 \ 2.25 ] ^T \ [ \mu M ] . \end{align*} \end{document}

Table 2.

Kinetic Parameter Values of System

Parameter	Value	Parameter	Value
k₁	7.353 × 10⁻³ [μM⁻¹h⁻¹]	k₁₁	86.79 × 10⁻³ [μM⁻¹h⁻¹]
k₂	7.078 [h⁻¹]	μ₁	1.0 [h⁻¹]
k₃	28.282 [μM⁻¹h⁻¹]	μ₂	1.0 [h⁻¹]
k₄	0.1158 [h⁻¹]	μ₃	1.0 [h⁻¹]
k₅	12.03 [μM⁻¹h⁻¹]	μ₄	1.0 [h⁻¹]
k₆	3.741 [μM⁻¹h⁻¹]	μ₅	1.0 [h⁻¹]
k₇	31.67 [μh⁻¹]	μ₆	1.0 [h⁻¹]
k₈	1.0 [μh⁻¹]	μ₇	1.0 [h⁻¹]
k₉	86.79 [μh⁻¹]	μ₈	1.0 [h⁻¹]
k₁₀	9.639 [μM × h⁻¹]

Note that the numerical values provided by the CRNT algorithm have been scaled to the micromolar concentrations in order to be consistent with previously published experimental results and numerical simulations.

With these values of the kinetic parameters, we have run in silico experiments to further validate the capability of our model to replicate the persistent memory feature evidenced by Acar et al. (2005). The results of the numerical simulations are shown in Figure 3. The preliminary phase where the cells are grown on different media, in the presence and in the absence of galactose, is taken into account by setting different initial concentrations of the molecular species in the GAL regulatory network. In the simulated experiments, we set all the initial concentrations at either one of the two equilibrium values \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $$x_{e_1}$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $$x_{e_2}$$ \end{document} , and then we perturb the initial concentration of one species, to check whether this perturbation leads the system to a different equilibrium point.

FIG. 3.

Evolution over time of the concentration of activated Gal3p starting from different initial conditions of Gal80p (a) and of internalized galactose (b). The simulations show that trajectories funnel into either one of the two stable steady states (denoted by dashed and solid lines for the low and the high equilibrium, respectively), depending on the initial conditions. The system converges at the low steady state, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $$x_{e_1}$$ \end{document} , regardless of the initial concentration of Gal80p, whereas different initial concentrations of internalized galactose lead the system to different equilibrium conditions.

In Figure 3, we have reported the perturbation experiments on Gal80p and on the internalized galactose, to illustrate how different species can have significantly different effects on the destabilization of the current operative condition. Figure 3a shows that the changes in the initial concentration of Gal80p do not make Gal3p active to switch to the high equilibrium condition; on the contrary, the initial concentration of internalized galactose strongly affects the response of the system. Indeed, Figure 3b demonstrates that the response of the system is bistable since the trajectories funnel into either one of the two steady states, depending on the initial value of the state variables, thus featuring a persistent memory of the past history of the system.

5.2. Characterization of the DA

In the previous section, we have established that the proposed model can exhibit two distinct steady states. Here, we want to characterize the DA of each of the equilibrium points, applying the analysis result devised in Amato et al. (2007) and recalled in Section 4.2.

Note that, in principle, one could try to compute the DAs through numerical experiments, varying the initial concentration of each species and simulating the system response to check which steady state each trajectory converges to. However, it is clear that this method can be successful only for very low-order systems, because the number of simulations increases combinatorially with the number of state variables. Moreover, the computational approach can only be based on a gridding of the state space, thus generally yielding very rough results.

Note that the two equilibrium points x_e₁, x_e₂ ≠ 0; therefore, we have used two changes of variables, before computing each DA estimate, as described in Amato et al. (2007). Then, we have defined the polytopes surrounding each point and have applied the conditions of Theorem 3 to solve Problem 2 for each polytope. In order to find the best possible estimates of the DAs, we have adopted a trial and error strategy, by iteratively enlarging the polytopes along different dimensions, stopping the iterations when conditions (9) could no longer be satisfied. The larger estimates obtained through this procedure are the boxes \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*} {\mathscr P}_1 = [ 53.86 , 73.86 ] \times [ 62.86 , 64.86 ] \times [ 0.0 , 4.0 ] \times [ 0.9 , 4.0 ] \times \\ \times [ 1.02 , 29.12 ] \times [ 5.64 , 8.74 ] \times [ 21.61 , 137.0 ] \times [ 0.0 , 1.5 ] , \end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*} {\mathscr P}_2 & = [ 84.01 , 130.01 ] \times [ 84.01 , 126.31 ] \times [ 6.57 , 356.67 ] \times [ 2.24 , 4.15 ] \times \\ & \times [ 0.62 , 10.12 ] \times [ 4.39 , 7.49 ] \times [ 62.76 , 62.96 ] \times [ 2.0 , 3.25 ] \end{align*} \end{document}

for the low and high equilibrium point, respectively. The LMI conditions (9) have been implemented and solved by means of the MATLAB Robust Control Toolbox (MathWorks, 2011).

To validate the DA estimates, we have run simulated experiments: as shown in Figure 4, initial conditions inside the box \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $${ \mathscr P}_1$$ \end{document} (resp. \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $${ \mathscr P}_2$$ \end{document} ) generate responses that converge to the equilibrium point x_e₁ (resp. x_e₂). Recall that the result of Theorem 3 is only sufficient; therefore, the computed polytopic regions are conservative estimates of the actual DAs. Thus, if the initial condition is set to a point that is outside the two boxes, we cannot predict where the trajectory will converge. Nevertheless, we have performed this converse test, by making some numerical simulations with initial conditions outside, but close to the frontiers of the boxes. These simulations, which in Figure 4 correspond to the dashed lines, show that each of these trajectories escapes away from the nearby box and converges to the other equilibrium point. Although further investigation would be needed to confirm this point, this seems to suggest that the estimated domains are close to the actual DAs.

FIG. 4.

State response of the concentrations of Gal4p (a) and Gal2p (b) from perturbed initial conditions belonging to the boxes \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $${ \mathscr P}_1$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $${ \mathscr P}_2$$ \end{document} . Solid (dashed) lines denote trajectories with initial condition included in \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $${ \mathscr P}_1$$ \end{document} (resp. \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $${ \mathscr P}_2$$ \end{document} ), converging to \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $$x_{e_1}$$ \end{document} (resp. \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $$x_{e_2}$$ \end{document} ). Star (circle) markers denote trajectories with initial condition outside \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\mathscr P}_1$$ \end{document} (resp. \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\mathscr P}_2$$ \end{document} ), converging to \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $$x_{e_2}$$ \end{document} (resp. \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $$x_{e_1}$$ \end{document} ).

6. Discussion

The dynamical behavior of the GAL system has been analyzed both numerically and analytically by several models already proposed in the literature. However, the mathematical models devised so far are not able to reproduce persistent memory of the past exposure to galactose, a key property of the GAL regulatory network, which has been experimentally evidenced in Acar et al. (2005). By exploiting CRNT to analyze the models based on mass-action kinetics, we have found that the model cannot exhibit the persistent memory feature because the structure of the underlying reaction networks is such that they cannot admit multiple steady states, regardless of the values assigned to the kinetic parameters.

Motivated by these findings, we have devised a novel model, still based on mass-action kinetics, of the GAL regulatory network, under certain simplifying assumptions. Using CRNT again, we have performed a structural bistability analysis of the proposed model, concluding that it can exhibit two steady states for some values of the kinetic parameters. Using the values obtained through the CRNT toolbox, the model has been simulated to replicate in silico the persistent-memory experiments conducted in Acar et al. (2005), confirming the capability to exhibit such features. Additionally, we have tested alternative structures of the reaction network, looking for simpler models; however, such tests have evidenced that further simplification of the model entails significant changes in the global dynamical behavior and, in particular, loss of bistability.

As determination of the multiple steady states is not sufficient to completely characterize the behavior of a dynamical nonlinear system, we have subsequently tackled the problem of estimating the DAs of the computed equilibrium points. Such a problem is hardly solvable for general nonlinear systems; however, when dealing with the class of quadratic systems, it is possible to exploit the results devised by our group in Amato et al. (2007). These conditions have enabled us to effectively compute two polytopic estimates of the DAs of the two distinct equilibrium points found by means of CRNT. The computation of the DAs is a key issue for a quantitative characterization of the behavior of a biological system; indeed, the DA of an equilibrium point provides valuable information about the capability of the system to adapt to exogenous perturbations or transcriptional noise.

The results reported in the present article provide further insights into the functioning of the GAL regulatory system, and an interesting development would be the experimental validation of the DA estimated on the basis of the proposed mathematical model.

On the other hand, we believe that the main contribution of this work consists in the introduction of a novel approach to the analysis of some system-level properties in biological systems such as persistent memory, switch-like biochemical response, and irreversible cell differentiation, which are based on the capability of the underlying biomolecular regulatory networks to exhibit multiple steady states. We surmise that the proposed approach, consisting of (a) the analysis for structural multistability of the reaction network and (b) the estimation of the DAs of the equilibrium points, which has been successfully applied to the analysis of the well-known GAL system in yeast, can be effectively generalized to other regulatory, signaling, and metabolic networks.

7. Appendix

7.1. Model developed in the present work

The GAL regulatory network model proposed in the present work comprises the following set of reactions:

\documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*} & G_3 + G_i \rightleftharpoons G_{3a} \rightarrow \emptyset \\ & G_{3a , 80} \rightarrow \emptyset \\ & G_{4 , 80} + G_{3a} \rightleftharpoons G_{3a , 80} + G_4 \\ & G_e + G_2 \rightarrow G_i + G_2 + G_3 \end{align*} \end{document}

From these reactions, assuming mass action kinetics, we can derive the complete dynamical model equations as follows: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*} \dot{G}_3 = k_9 G_4 - k_1 G_3 G_i + k_2 G_{3a} - \mu_1 G_3 ( { \rm A.1a} ) \end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*} & \dot{G}_i = k_{11}G_e G_2 - k_1 G_3 G_i + k_2 G_{3a} - \mu_8 G_i \tag{ {\rm A.1b} } \\ & \dot{G}_{3a} = k_1 G_3 G_i - k_2 G_{3a} - k_5 G_{4 , 80}G_{3a} + k_6 G_{3a , 80}G_4 - \mu_3 G_{3a} \tag{ {\rm A.1c} } \\ & \dot{G}_4 = k_{10} - k_3 G_4 G_{80} + k_4 G_{4 , 80} + k_5 G_{4 , 80} G_{3a} - k_6 G_{3a , 80} G_4 - \mu_4 G_4 \tag{ {\rm A.1d} } \\ & \dot{G}_{80} = - k_3 G_4 G_{80} + k_4 G_{4 , 80} + k_7 G_4 - \mu_2 G_{80} \tag{ {\rm A.1e} } \\ & \dot{G}_{4 , 80} = k_3 G_4 G_{80} - k_4 G_{4 , 80} - k_5 G_{4 , 80} G_{3a} + k_6 G_{3a , 80} G_4 - \mu_6 G_{4 , 80} \tag{ {\rm A.1f} } \\ & \dot{G}_{3a , 80} = k_5 G_{4 , 80} G_{3a} - k_6 G_{3a , 80} G_4 - \mu_7 G_{3a , 80} \tag{ {\rm A.1g} } \\ & \dot{G}_2 = k_8 G_4 - \mu_5 G_2 \tag{ {\rm A.1h} } \\ & G_e = { \rm const} \tag{ {\rm A.1i} } \end{align*} \end{document}

7.2. Smidtas' and Kulkarni's models

The detailed model proposed in Smidtas et al. (2006) includes the following equations: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*} & c \times Gal4p / DNA \rightarrow c \times Gal4p / DNA + Gal3p + Gal80p \tag{ {\rm A.2a} } \\ & Gal4p / DNA + Gal80p \rightleftharpoons Gal80p / Gal4p / DNA \tag{ {\rm A.2b} } \\ & Gal3p + Gal80p \rightleftharpoons Gal3p / Gal80p \tag{ { \rm A.2c} } \\ & Gal3p / Gal80p + Galactose \rightarrow Glactose / Gal3p / Gal80p \tag{ {\rm A.2d} } \\ & Gal3p \rightarrow \emptyset \tag{{\rm A.2e} } \\ & Gal80p \rightarrow \emptyset \tag{ {\rm A.2f} } \end{align*} \end{document}

In the presence of galactose, Gal3p, Gal80p binds to galactose, thus decreasing the binding of Gal80p to Gal4p. Reduction of Gal80p activates Gal4p and induces transcription of the GAL3 and GAL80 genes. This closes the feedback loop, and Gal3p and Gal80p shift the equilibrium back towards Gal4p inactivation (Gal4p/80p). By means of CRNT, we can obtain that, for arbitrary kinetics (subject to very weak constraints), the corresponding differential equations cannot admit a steady state at which all species concentrations are positive, nor can they admit a cyclic composition trajectory that passes through a point at which all species concentrations are positive. In particular, there is no set of rate constants for which the resulting mass action system can admit two distinct strictly positive steady states that are stoichiometrically compatible or a degenerate steady state.

Kulkarni et al. (2010) tested the dynamical nonlinear quadratic model proposed in Smidtas et al. (2006) and derived an ODE model, including more chemical reactions: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*} & G4 \rightleftharpoons G4^* + GAL3 / 80 \tag{ {\rm A.3a} } \\ & R + Gal \rightleftharpoons G4^* + BR \tag{ {\rm A.3b} } \\ & G4* + GAL3 / 80 \rightleftharpoons R \tag{ {\rm A.3c} } \\ & G4^* \rightarrow \emptyset \tag{ {\rm A.3d} } \\ & GAL3 / 80 \rightarrow \emptyset \tag{ {\rm A.3e} } \\ & R \rightarrow \emptyset \tag{ {\rm A.3f} } \\ & BR \rightarrow \emptyset { {\rm A.3g} } \\ & Gal \rightarrow \emptyset \tag{ {\rm A.3h} } \end{align*} \end{document}

where GAL3/80 is the complex in which are lumped Gal3p and Gal80p; R is the inactive complex formed by Gal3p, Gal80p, and Gal4p; BR is the bound receptor comprising Gal3p, Gal80p, and galactose; and G4* represents the activated transcriptional activator Gal4p protein. The total Gal4p concentration is assumed to be constant during the GAL response.

By means of CRNT, we can still obtain that, for arbitrary kinetics (subject to very weak constraints), the corresponding system cannot admit a steady state at which all species concentrations are positive, nor can they admit a cyclic composition trajectory that passes through a point at which all species concentrations are positive. In particular, there is no set of rate constants for which the resulting mass action system can admit two distinct strictly positive steady states that are stoichiometrically compatible or a degenerate steady state.

Footnotes

Disclosure Statement

No competing financial interests exist.

References

Acar

, Becskei

, van Oudenaarden

2005. Enhancement of cellular memory by reducing stochastic transitions. Nature, 435:228–232.

Amato

2006. Robust Control of Linear Systems Subject to Uncertain Time-Varying Parameters. Springer Verlag: Berlin.

Amato

, Cosentino

, Merola

2007. On the region of attraction of nonlinear quadratic systems. Automatica, 43:2119–2123.

Angeli

, Ferrel

J.E.

, Sontag

E.D.

2004. Detection of multistability, bifurcation, and hysteresis in a large class of biological positive-feedback systems. Proc. Natl. Acad. Sci. USA, 101:1822–1827.

Bagowski

C.P.

, Ferrell

J.E.

2001. Bistability in the JNK cascade. Curr. Biol., 11:1176–1182.

Balázsi

, van Oudernaarden

, Collins

J.J.

2011. Cellular decision making and biological noise: from microbes to mammals. Cell, 144:910–925.

Bhat

P.J.

2008. Galactose Regulon of Yeast: From Genetics to Systems Biology. Springer-Verlag: Berlin.

Boyd

S.P.

, El Ghaoui

, Feron

et al. 1994. Linear Matrix Inequalities in System and Control Theory. SIAM Press: Philadelphia.

Chickarmane

, Troein

, Nuber

U.A.

et al. 2006. Transcriptional dynamics of the embryonic stem cell switch. PLoS Comput. Biol., 2:e123.

10.

Cosentino

, Bates

D.G.

2011. Feedback Control in Systems Biology. CRC Press: Boca Raton, FL.

11.

Craciun

, Feinberg

2005. Multiple equilibria in complex chemical reaction networks. I. The injectivity property. SIAM J. Appl. Math., 65:1526–1546.

12.

de Atauri

, Orrell

, Ramsey

et al. 2004. Evolution of design principles in biochemical networks. IEEE Proc. Syst. Biol., 1:28–40.

13.

Ellison

1998. The advanced deficiency algorithm and its applications to mechanism discrimination [Ph.D. dissertation] University of Rochester: Rochester, NY.

14.

Feinberg

1988. Chemical reaction network structure and the stabiliyt of complex isothermal reactors. II. Multiple steady states for network of deficiency one. Chem. Eng. Sci., 43:1–25.

15.

Feinberg

1987. Chemical reaction network structure and the stability of complex isothermal reactors. I. The deficiency zero and deficiency one theorems. Chem. Eng. Sci., 42:2229–2268.

16.

Ferrell

J.E.

Jr. , Machleder

E.M.

1998. The biochemical basis of an all-or-none cell fate switch in Xenopus oocytes. Science, 280:895–898.

17.

Gardner

T.S.

, Cantor

C.R

, Collins

J.J.

2000. Construction of a genetic toggle switch in Escherichia coli. Nature, 403:339–342.

18.

Johnston

1987. A model fungal gene regulatory mechanism: the GAL genes of Saccharomyces cereviasiae. Microbiol. Rev., 51:458–476.

19.

Khalil

M.K.

1991. Nonlinear Systems. Macmillan: New York.

20.

Klein

C.J.L.

, Olsson

, Nielsen

1998. Glucose control in Saccharomyces cerevisiae: the role of MIG1 in metabolic functions. Microbiology, 144:13–24.

21.

Kulkarni

V.V.

, Kareenhalli

, Malakar

et al. 2010. Stability analysis of the GAL regulatory network in Saccharomyces cerevisiae and Kluyveromyces lactis. BMC Bioinform., 11:S43.

22.

La Salle

J.P.

1960. Some extensions of Lyapunov's second method. IRE Trans. Circuit Theor, CT-7:520–527.

23.

Lohr

, Venkov

, Zlatanova

1995.

24.

Transcriptional regulation in the yeast GAL gene family: a complex genetic network. FASEB J., 9:777–786.

25.

MathWorks. 2011. MATLAB Robust Control Toolbox. MathWorks, Inc.: Natick, MA.

26.

Melcher

, Xu

H.E.

2001. Gal80-Gal80 interaction on adjacent Gal4p binding sites is required for complete GAL gene repression. EMBO J., 20:841–851.

27.

Murray

J.D.

2001. Mathematical Biology: I. An Introduction. Springer: New York.

28.

Pannala

V.R.

, Bhat

P.J.

, Bhartiya

et al. 2010. System biology of GAL regulon in Saccharomyces cerevisiae. WIREs Syst. Biol. Med., 2:98–106.

29.

Peng

, Hopper

J.E.

2002. Gene activation by interaction of an inhibitor with a cytoplasmic signaling protein. Proc. Natl. Acad. Sci. USA, 99:8548–8553.

30.

Pomenering

J.R.

, Sontag

E.D.

, Ferrell

J.E.

Jr.

2003. Building a cell cycle oscillator: hysteresis and bistability in the activation of Cdc2. Nat. Cell Biol., 5:346–351.

31.

Sha

, Moore

, Chen

et al. 2003. Hysteresis drives cell-cycle transitions in Xenopous laevis egg extracts. Proc. Natl. Acad. Sci. USA., 100:975–980.

32.

Smidtas

, Schachter

, Képès

2006. The adaptive filter of the yeast galactose pathway. J. Theor. Biol., 242:372–381.

33.

Timson

D.J.

, Ross

H.C.

, Reece

R.J.

2002. Gal3p and Gal1p interact with the transcriptional repressor Gal80p to form a complex of 1:1 stoichiometry. Biochem. J., 363:515–520.