- Research article
- Open Access
Threshold-dominated regulation hides genetic variation in gene expression networks
- Arne B Gjuvsland^{1, 3}Email author,
- Erik Plahte^{2, 3} and
- Stig W Omholt^{1, 3}
https://doi.org/10.1186/1752-0509-1-57
© Gjuvsland et al; licensee BioMed Central Ltd. 2007
- Received: 29 June 2007
- Accepted: 06 December 2007
- Published: 06 December 2007
Abstract
Background
In dynamical models with feedback and sigmoidal response functions, some or all variables have thresholds around which they regulate themselves or other variables. A mathematical analysis has shown that when the dose-response functions approach binary or on/off responses, any variable with an equilibrium value close to one of its thresholds is very robust to parameter perturbations of a homeostatic state. We denote this threshold robustness. To check the empirical relevance of this phenomenon with response function steepnesses ranging from a near on/off response down to Michaelis-Menten conditions, we have performed a simulation study to investigate the degree of threshold robustness in models for a three-gene system with one downstream gene, using several logical input gates, but excluding models with positive feedback to avoid multistationarity. Varying parameter values representing functional genetic variation, we have analysed the coefficient of variation (CV) of the gene product concentrations in the stable state for the regulating genes in absolute terms and compared to the CV for the unregulating downstream gene. The sigmoidal or binary dose-response functions in these models can be considered as phenomenological models of the aggregated effects on protein or mRNA expression rates of all cellular reactions involved in gene expression.
Results
For all the models, threshold robustness increases with increasing response steepness. The CV s of the regulating genes are significantly smaller than for the unregulating gene, in particular for steep responses. The effect becomes less prominent as steepnesses approach Michaelis-Menten conditions. If the parameter perturbation shifts the equilibrium value too far away from threshold, the gene product is no longer an effective regulator and robustness is lost. Threshold robustness arises when a variable is an active regulator around its threshold, and this function is maintained by the feedback loop that the regulator necessarily takes part in and also is regulated by. In the present study all feedback loops are negative, and our results suggest that threshold robustness is maintained by negative feedback which necessarily exists in the homeostatic state.
Conclusion
Threshold robustness of a variable can be seen as its ability to maintain an active regulation around its threshold in a homeostatic state despite external perturbations. The feedback loop that the system necessarily possesses in this state, ensures that the robust variable is itself regulated and kept close to its threshold. Our results suggest that threshold regulation is a generic phenomenon in feedback-regulated networks with sigmoidal response functions, at least when there is no positive feedback. Threshold robustness in gene regulatory networks illustrates that hidden genetic variation can be explained by systemic properties of the genotype-phenotype map.
Keywords
- Gene Regulatory Network
- Parameter Perturbation
- Robustness Property
- Singular Component
- Singular Variable
Background
Historical perspective
In the early 1970s, Leon Glass, Stuart Kauffman, and René Thomas started their pioneering efforts in exploring the possibility of modelling what was then called "Genetic Control Circuits" (Thomas) and "Biochemical Control Networks" (Glass and Kauffman) by using concepts and ideas from mathematical logic. Combining these ideas with earlier ideas from Monod and others on allostery and cooperativity which suggested sigmoidal rate dependences of key metabolites, Glass and Kauffman [1, 2] and Thomas [3] proposed that gene transcription could be modelled by sigmoidal response functions depending on transcription factor concentrations. In the case of several transcription factors acting on a gene, they assumed the effect could be expressed by Boolean combinations of the separate response functions, and proposed a simple framework of ordinary differential equations for modelling of gene regulatory networks based on these principles. Glass and Kauffman observed that the behaviour of these regulatory networks was remarkably insensitive to the steepness of the sigmoids, and suggested to use Heaviside or step function in stead of sigmoids as dose-response functions to simplify the models and their analysis.
From these early attempts, phenomenologic frameworks for the modelling of Gene Regulatory Networks (GRNs) have been developed, based on a few fundamental premises: (i) genes are controlled by transcription factors (TFs) which combine into logical input functions, and these can be described by Boolean logic; (ii) the effect of a transcription factor on the transcription rate of a gene (the response function) can be described by a sigmoidal function of its concentration with a pronounced threshold behaviour (graded response) or by a Heaviside step function (binary response); (iii) this can be modelled in a discrete way in which transcription factors are either absent of present, and proteins are either transcribed or not, or in a continuous way by means of ordinary differential equations; (iv) proteins act as transcription factors, so that networks become closed with feedback loops; (v) posttranscriptional, translational and posttranslational regulation, transport processes, metabolic processes etc. can be phenomenologically encompassed by the sigmoidal or binary response functions.
Recent experimental evidence seem to confirm the validity and usefulness of the basic assumption on which GRNs are resting. A number of studies have shown that gene networks have cis-regulatory elements governed by a Boolean-like logic [4–11]. There is also extensive experimental documentation of sigmoidal or binary responses in gene regulation [5, 11–16] as well as theoretical justifications based on classical principles from physical chemistry [17–20].
How common is steep transcription response? Analyses based on classical methods from statistical physical chemistry show that a steep transcription response curve could be due to cooperativity in the transcription factor binding [18–20]. It has been shown that transcriptional and signalling cascades do in fact lead to graded or binary responses [21, 22]. There is also extensive evidence that transcription response in single cells is binary (see references in [14]), and that individual cells responds in an on/off way to varying external inputs [13]. Thus, there are good reasons to expect that high gain regulation is quite common in gene regulatory networks.
The method developed in [23–27] to deal with models with steep sigmoidal response functions works for quite general models, also with other nonlinearities in addition to the steep sigmoidal functions. In the course of this work it was discovered that when the responses functions are very steep, equilibrium values for actively regulating variables show a remarkable robustness towards changes in all parameters except the level of the threshold around which the active regulation occurs. We call this phenomenon threshold robustness. To be precise, this is a mathematical result valid in the limit when the sigmoid function approaches a step function (Heaviside function), but for continuity reason it is also valid for steep sigmoid functions. But to what extent is it found in models with more empirically sound threshold functions? We have investigated this question by a simulation study of a wide class of 3-dimensional regulatory systems where the regulatory dose-response relationships are varied from a hyperbolic Michaelis-Menten situation to an extremely steep sigmoidal situation.
If conserved when the steepness of the sigmoidal interactions is slackened to realistic values, insensitivity or robustness to functional genetic polymorphism may be a deep generic property of some of the loci in a wide range of regulatory networks. When present, threshold robustness adds significant and characteristic phenomena to the genotype-phenotype map. This implies for example that the functional mutational changes in network which shows threshold robustness will only results in small phenotypic variations in the homeostatic values of the protein products.
Analytical foundation
where y_{ j }is the concentration of gene product number j, j = 1,...,n, Z_{ j }= S(y_{ j }, θ_{ j }, p_{ j }) is a sigmoid or binary function with threshold θ_{ j }and steepness parameter p_{ j }, and y and Z are the vectors with y_{ j }and Z_{ j }as components. The functions R_{ j }∈ [0, 1] and Q_{ j }∈ ⟨0, 1] are regulatory functions, frequently taken to be algebraic equivalents of Boolean functions [25], describing the regulation of production and decay, respectively, while the positive parameters κ_{ j }and λ_{ j }represent the maximal production and decay rates.
Eq. (1) could be justified in at least two ways. It could be considered a model of transcription regulation with the y_{ j }still representing protein concentrations. This model could be derived from a larger model for protein and mRNA concentrations where transcription of mRNA is regulated by protein concentrations, and the conversion from mRNA to protein is described by linear equations. If all mRNA degradation rates are much larger than all protein degradation rates, we can apply a quasi-stationary hypothesis to the mRNA concentrations, leading to Eq. (1). This procedure can be justified mathematically as well as biologically. A simple example is presented below, and the case n = 2 is studied in [28].
Alternatively, taken as a model of gene regulation, Eq. (1) is a generic phenomenological model of protein concentration dynamics, not a mechanistic description of gene regulation. The threshold functions model the aggregated effect of all the processes involved in the real cellular regulatory networks [29]: transcription, translation, intracellular transport, post-translational modifications, protein-protein interactions, metabolic processes, and signalling cascades. Such drastic simplification is hard to justify theoretically, but models based on the generic Eq. (1) have been applied successfully to many real systems [5, 30, 31]. Considered in this way, Eq. (1) is a generic, phenomenological framework assumed to catch the essential features of a wide range of regulatory systems, where the regulatory control may be at the level of transcription, mRNA stability, translation, or post-translation, and where the state variables may for example be concentrations of proteins, hormones, mRNA, and intracellular ions [29].
Models encompassed by Eq. (1) have been extensively investigated mathematically, in particular in the special case of S_{ j }being step functions. (See [32] for an extensive list of references.) An efficient way to analyse models of the generic type Eq. (1) with steep sigmoids is presented in [26, 27].
Frequently models encompassed by the generic type Eq. (1) have stationary points lying close to the thresholds of one or several variables when the sigmoids are steep, i.e. when p is large. A stationary point y*(p) is called a singular stationary point (SSP) for Eq. (1) if at least one of its components approaches its threshold when p → ∞. These components are called singular, the others regular. It has been proved that if a SSP Y* exists in the step function limit, then there exists a stationary point y*(p) for sufficiently large p with the property that y*(p) → Y* when p → ∞ [26]. Mathematically, SSPs of Eq. (1) have remarkable generic robustness properties. The key issue is that the singular components of Y* are locked to their thresholds, while the regular components vary with the other parameters in the model. For y*(p) this implies that when p is large, the singular components are highly insensitive to changes in all parameters except the thresholds of the singular variables. Biologically this insensitivity means that the expression levels of genes corresponding to singular components will be kept virtually constant despite stochastic or mutational variation in the expression process rates.
where m and y are the concentrations of mRNA and protein, respectively, Z = S(y, θ, p) is the sigmoidal response function, r is the basal transcription rate, θ is the regulation threshold, p is the steepness parameter, and the remaining four parameters are production and degradation rate constants. Incidentally, asssuming quasistationarity of mRNA concentration as explained above, leads to $\dot{y}$ = (κ/γ)[r + α(1 - Z)] - λy, which is of the generic type Eq. (1).
However, if g is perturbed outside this interval, y* is no longer locked to the threshold, and robustness is lost. If p → ∞ and g <r/α, then Z* = 1 (the magenta case), m* = r/γ, y* = κr/(γλ), while if g > 1 + r/α, then Z* = 0 (the green case), m* = (α + r)/γ, y* = κ(α + r)/(γλ). We see that in both these cases, y* vary with the other parameters.
Model system
For our simulations we chose a particular realisation of Eq. (1) which after a scaling is given by the dimensionless equations
x_{ j }' = α_{ j }R_{ j }(Z_{1}, Z_{2}) - γ_{ j }x_{ j }, j = 1,...,3,
The key question is whether the robustness of SSPs is still generic and preserved when the Hill exponents atttain smaller, more realistic values. Considering instances of Eq. (1) as models for gene regulatory networks, we checked this out for a large set of 3-dimensional particular realisations of Eq. (3). We took the set of stable equilibrium values ${x}_{j}^{\ast}(p)$ as the phenotype and the set {α_{ j }, γ_{ j }, θ_{ j }, p, R_{ j }}, j = 1,...,3, of parameter values and regulatory functions as the genotype. The equilibrium conditions for Eq. (3) then define the genotype-phenotype map for this system. Our interest is to investigate the robustness properties of the phenotype under mutations, i.e. under perturbations of the genotype. For the 14 models with a unique SSP in the step function limit we investigated the robustness properties of the singular and the regular components of the SSP for a range of parameter perturbations and for varying Hill exponent of the response functions from p = 1 (Michaelis-Menten conditions) to p = 100, which for all practical purposes is very close to a step function.
The coefficient of variation CV for a distribution is defined as the ratio of the standard deviation to the mean. Being a dimensionless number and scaled by the mean, it is suitable for comparison of the variation of distributions with large differences in mean values. To measure and compare the CV s for the equilibrium values ${x}_{i}^{\ast}$ for each of the 14 models we generated 81 parameter sets, giving a total of 1134 different systems with a unique SSP for which x_{1} and x_{2} are singular and x_{3} regular. For each data set we sampled 50 random pertubations of each production parameter α_{ j }from the uniform distribution U(α_{ j }/2, 3α_{ j }/2) with corresponding coefficient of variation CV_{uni} ≈ 0.288. We then computed the coefficients of variation CV_{ j }^{ k }, j = 1, 2, 3, k = 1,...,81 for the steady state levels of all three variables in all the 14 networks separately. Details are described in the Methods section. We use the minimum coefficients of variation as robustness measure in order to be able to compare the robustness of most favourable parameter sets.
Results
What is more important is that in almost all cases, all the coefficients of variation for genes 1 and 2 are significantly smaller than those for gene 3. The reduced sensitivity of gene 3 when p decreases towards 1 can be seen as a consequence of the fact that when the sigmoids in the rate equation for x_{3} are slackened, a certain variation in the inputs gives smaller variation in the response. However, for all the models except Model 14 there are parameter sets for which the system is not robust, which shows that robustness is generally only present for a certain range of parameter values. The decisive factor is how easily a perturbation shifts the position of the stable point away from the switching domain. If this happens, the gene's status in the regulatory system is changed: it is no longer an active regulator and its robustness is lost. Accordingly, the gene that was regulated is now either off or constitutively on unless it is still effectively regulated by the other gene.
To further illustrate the difference between the least robust Models 1, 2, and 11 and the rest of the models, we computed for each realisation the number of perturbations N for which the stable point is a SSP with both x_{1} and x_{2} as active regulators. We consider x_{ j }an active regulator if ${Z}_{j}^{\ast}$ lies in the interval ⟨0.05, 0.95⟩. Among all 81 parameter sets the highest observed N was 37, 37, and 35 for Models 1, 2, and 11 respectively. For all other models one can always find a parameter set with N = 50, i.e. for which all perturbations render a SSP with a high degree of robustness in both ${x}_{1}^{\ast}$ and ${x}_{2}^{\ast}$.
Discussion
A number of different sources of robustness in cellular function and biochemical networks are discussed in the literature (see e.g. [35, 36]). Considered as as a systemic property of a developmental or functional unit in an organism, robustness has been explained as a consequence of both negative and positive feedback and several other network properties [30, 37, 38], as a consequence of network topology and connectivity [39–41] or of modularity or redundancy of the network [42, 43]. Conversely, complexity has been seen as a consequence of selection for robustness rather than the other way round [44, 45].
Distinguishing the phenomenon of robustness from homeostasis and stability which concern the system's ability to maintain a stable state, Kitano defines robustness as "a property that allows a system to maintain its functions against external perturbations [37, 46]." The function maintained by threshold robustness is the gene's or gene product's ability to act as an active regulator of itself or another gene despite parameter perturbations, which can be seen as consequences of external noise.
A variable x_{ j }exhibits threshold robustness when it is a singular (also called switching) variable of a SSP, in other words, when the stationary value ${x}_{j}^{\ast}$ approaches the threshold θ_{ j }as the steepness of the associated sigmoid function tends to infinity. In that case the stationary value is locked to the threshold whatever the value of the other parameters (as long they are not perturbed so strongly that x_{ j }is no longer a singular variable in the perturbed system). Then the equilibrium value of the response function ${Z}_{j}^{\ast}$ is neither close to 0 nor to 1, these cases corresponding to the gene in question being either constitutively off or constitutively on. We illustrated this by our analysis of Model 2 given above. For all parameter perturbations which maintain the parameter values within the shaded area in Fig. 7, both variables are singular and locked to their threshold, and robustness is preserved. Only if the equilibrium point is perturbed outside the shaded area, one or both ${x}_{j}^{\ast}$ slip away from threshold onto one of the flatter part of the response curve where corresponding Z_{ j }-value approaches 0 or 1 and is rather insensitive to variations in the x_{ j }-value. Then the x-variable is no longer an active regulator, the previously regulated gene now being either constantly off or constantly on, and its robustness is gone. This analysis pertains to the limiting case of very steep responses functions, but is approximately valid when the response is more gradual.
Thus, threshold robustness of a transcription factor or gene product stems from the fact that it is actively regulating itself or other genes in the stable equilibrium state. There is a general theorem stating that the system of equations F(x) = g, where x ∈ R^{ n }, F is a differentiable function F : R^{ n }→ R^{ n }, and F(x_{0}) = g_{0}, can only have a differentiable solution x = G(g), x_{0} = G(g_{0}), if there is a feedback loop involving all n variables [27, 47]. When this is applied to a SSP, it follows that for the SSP to exist there must be a feedback loop among all the singular variables, mediated by the sigmoidal terms Z_{ j } [27]. Thus, at least for the models investigated in the present paper, if there is a sigmoid-mediated feedback loop among a subset of the variables, and there is a SSP in which these variables are singular, and the sigmoids involved are sufficiently steep, the system will exhibit threshold robustness in this SSP for all the singular variables. Mathematically, threshold robustness is not restricted to any particular type of feedback system. Rather, it is a generic feature of GRNs with steep response functions. Due to the generality of the above-mentioned theorem and the concept and properties of SSPs [25, 26], we conjecture that threshold robustness is a general property of singular stationary points. Our findings suggest that this feature is generic in a wider sense, not being dependent on response functions being steep, but it becomes weaker as the system approaches Michaelis-Menten conditions. Also, we may conclude that threshold robustness is a systemic property effectuated by the feedback loop between the sigmoidal interaction terms of the variables at threshold.
The fact that an actively regulating variable necessarily is part of a feedback loop explains how it can maintain its regulation ability. Due to the feedback and because the point is stable, a perturbation that shifts a singular variable away from its threshold, will eventually be counteracted. Thus, a singular variable, i.e. a variable with threshold robustness, is actively regulating and also itself being actively regulated. On the other hand, a feedback loop, even between sigmoidal terms, do not necessarily imply that a SSP exists. It may for certain parameter value combinations but not for others. Thus, it would not be fair to say that threshold robustness is a necessary consequence of feedback. Rather, it is a property of a SSP, and the SSP can only exist if there is a sigmoid-mediated feedback loop among its singular variables. René Thomas has suggested the reasonable conjecture that this must be a negative feedback loop [48]. Our results seem to support this conjecture and suggest that threshold robustness is maintained by negative feedback.
In our numerical simulations, we chose both thresholds fixed at θ_{ j }= 1. However, this was just a matter of convenience and not a model limitation. Without this choice, all parameter variations would still have been expressed by the same two parameters, but with the modification that μ_{ j }= α_{ j }/(γ_{ j }θ_{ j }). Thus, variation in the thresholds would also result in a variation of μ_{ j }, the only additional effect being a shift in the position of the SSP as long as the perturbed parameters would not fall outside the domain in parameter space where both variables are singular. The system would still exhibit threshold robustness, but now around the shifted threshold values. For the majority of the models this domain covers the major part of the parameter space (c.f. Fig. 6). For this reason threshold robustness might also be called adaptive robustness.
Seen as genotype-phenotype maps the studied gene regulatory networks share some interesting features related to the genetic phenomena hidden genetic variation and neutrality. Hidden (or cryptic) genetic variation is defined as "standing genetic variation that does not contribute to the normal range of phenotypes observed in a population, but that is available to modify a phenotype that arises after environmental change or the introduction of novel alleles [49]." It is known that negative autoregulation hides variation in the copy numbers of genes [33]. Our results show a more general connection between negative feedback, threshold regulation and hidden genetic variation. In our simulations we fixed a certain amount of parameter variation, corresponding to genetic variation, for a range of different gene regulatory models, and showed that this variation does not result in a corresponding variation in the phenotype. Thus, in a network which exhibits threshold robustness, functional mutations are hidden for phenotypic selection. Our results imply that mutations causing changes in the maximum production rate or the relative decay rate, but keeping the threshold of regulation intact, may have almost no phenotypic signature if the regulatory dose-response relationships are steep enough. Such mutations are neutral in the sense of Wagner's definition: "A neutral mutation does not change a well defined aspect of a biological system's function in a specific environment and genetic background [36]." It is implicit in this definition that neutral mutations may aquire a phenotypic signature if the system conditions change. Threshold robustness also offers an explanation of how genotypic variation that is hidden under one condition may be released by for instance a mutation causing a functional change in the regulatory machinery. For instance, the hidden genetic variation could be released by a change in the regulatory structure beyond the limit inside which robustness for one or several variables is conserved. That would turn singular components into regular components without threshold robustness and susceptible to parameter variations. Thus, a single key mutation in a regulatory structure may release a substantial amount of hidden, potential variation, as illustrated in a simulation study by Bergman and Siegal [50].
Conclusion
This paper presents a simulation study of a class of gene regulatory models in which regulation is modelled with sigmoidal response functions combined by operators mimicking Boolean functions. From a mathematical analysis it is known that when the sigmoidal response functions are very steep, the equilibrium values of the regulating agents are locked to the thresholds, thus are very insensitive to perturbations in all parameters except the threshold levels. This implies that they retain their active regulating power despite parameter perturbations. Our simulations show that this threshold robustness is preserved also for more gentle responses, and is qualitatively present even under Michaelis-Menten conditions. Even though the models investigated are simple, there are reasons to believe that they give a phenomenological description of a large number of different system in which the aggregated effects of a series of transcriptional, translational, and post-translational processes as well as protein-protein interactions and metabolic processes can be described by threshold-dominated response functions.
According to Kitano, robustness is an ability to maintain a function under noise-like perturbations. Threshold robustness is the ability of a protein or transcription factor to maintain an active regulation of a gene in homeostasis under external perturbations. The feedback loop that the system necessarily possesses in the homeostatic state, ensures that the robust members of the loop are themselves regulated and kept close to their threshold values.
Some of the 14 models investigated show a much lower degree of robustness than the rest, a fact that we have explained by a specific analysis of the models in which we compared the shape and extension of the parameter space domain in which the robustness property is preserved under parameter perturbations. For the models with lower degree of robustness, this domain is smaller and concavely shaped such that the perturbed values very easily fall outside the robustness domain. We have also seen that one variable may be considerably less robust than the other, despite all response functions being equally steep, and have offered an explanation of this phenomenon.
Threshold robustness may offer increased insight into genetic phenomena such as maintenance and release of genetic variation in evolution, but a closer investigation of these matters was beyond the scope of the present paper.
Methods
Model equations and regulatory functions
was scaled to give the scaled equations (3)
x_{ j }' = α_{ j }R_{ j }(Z_{1}, Z_{2}) - γ_{ j }x_{ j }, j = 1,...,3,
where α_{3} ∈ ⟨0, 1⟩ and all γ_{ j }∈ ⟨0, 1⟩, and Z_{ j }= x_{ j }^{ p }/(x_{ j }^{ p }+ 1). Details are given below.
The 12 Boolean functions used and their algebraic equivalents
Boolean function | Algebraic equivalent | |
---|---|---|
1 | X_{1} AND X_{2} | Z _{1} Z _{2} |
2 | X_{1} AND $\overline{\text{X}}$_{2} | Z_{1}(1 - Z_{2}) |
3 | $\overline{\text{X}}$_{1} AND X_{2} | (1 - Z_{1})Z_{2} |
4 | $\overline{\text{X}}$_{1} AND $\overline{\text{X}}$_{2} | (1 - Z_{1})(1 - Z_{2}) |
5 | X_{1} | Z _{1} |
6 | X_{2} | Z _{2} |
7 | $\overline{\text{X}}$ _{1} | 1 - Z_{1} |
8 | $\overline{\text{X}}$ _{2} | 1 - Z_{2} |
9 | X_{1} OR X_{2} | Z_{1} + Z_{2} - Z_{1}Z_{2} |
10 | $\overline{\text{X}}$_{1} OR $\overline{\text{X}}$_{2} | 1 - Z_{1}Z_{2} |
11 | X_{1} OR $\overline{\text{X}}$_{2} | 1 - Z_{2} + Z_{1}Z_{2} |
12 | $\overline{\text{X}}$_{1} OR X_{2} | 1 - Z_{1} + Z_{1}Z_{2} |
Scaling
Here we outline the scaling steps to non-dimensionalise the model equations (4) and standardise the parameter ranges.
1. For j = 1, 2, introduce x_{ j }= y_{ j }/θ_{ j }to scale the thresholds for x_{ j }to θ_{ j }= 1.
2. Assume λ_{ j }∈ ⟨0, Λ⟩, j = 1, 2, 3, and scale the time t by introducing τ = Λt.
3. Assume κ_{3} ∈ ⟨0, K⟩ and introduce x_{3} = (K/κ_{3})y_{3}.
This leads to the scaled equations (3). The new parameters are α_{ j }= κ_{ j }/(Λθ_{ j }), j = 1, 2, α_{3} = κ_{3}/K, and γ_{ j }= λ_{ j }/Λ, j = 1, 2, 3, with the properties α_{3} ∈ ⟨0, 1⟩ and all γ_{ j }∈ ⟨0, 1⟩. Below we use μ_{ j }= α_{ j }/γ_{ j }, j = 1, 2.
Sampling
For each of the 14 models we generated 81 initial parameter sets φ^{ k }= {α_{1k}, γ_{1k}, α_{2k}, γ_{2k}, α_{3k}, γ_{3k}} by taking the equilibrium values (${Z}_{1}^{\ast}$, ${Z}_{2}^{\ast}$) from the square lattice {0.1, 0.2,...,0.9} × {0.1, 0.2,...,0.9}.
Having drawn a parameter set we introduced variation in the parameters. For each parameter set φ^{ k }we generated a set of 50 perturbed parameter values ${\tilde{\phi}}_{\ell}^{k}$, ℓ = 1,...,50 by keeping γ_{ j }fixed and sampling α_{ j }of all three genes uniformly in the range ⟨α_{ j }/2, 3α_{ j }/2⟩. This is justified by the fact that the stable state values only depend on the ratios μ_{ j }. For each set of perturbed parameter values we used Matlab's odesolver to find the steady state ${x}_{j\ell}^{\ast k}$ of the three gene products in Eq. (3).
Measuring robustness towards pertubation
As the coefficient of variation is invariant under scaling, the parameter CV s will be close to the CV of the uniform distribution U(0.5, 1.5) which is CV_{uni} = 1/$\sqrt{12}$ ≈ 0.288. Since we wanted to study how the robustness depends on the steepness of the dose-response function, we carried out this procedure for p = 1, 2, 5, 10, 20, 50, 100.
Explaining the differences in robustness among the models
The reason for the reduced robustness of Models 1, 2, and 11 compared to the remaining models can be explained by investigating the shape and size of the robustness domain Ω_{SSP} in the (μ_{1}, μ_{2})-plane for which both x_{1} and x_{2} are active regulators in the stable point in the step function limit. We find that for Models 1, 2, and 11, Ω_{SSP} is a narrow, concave stripe, while for the remaining models it is a convex domain (Fig. 6). Below we show how to derive this result for Model 1. The analyses for the rest of the models are similar.
Then Ω_{SSP} is the domain covered by this family of curves when ${Z}_{2}^{\ast}$ vary between 0 and 1. One easily finds Ω_{SSP} to be the domain between the three curves μ_{1} = 1, μ_{2} = 1, and μ_{2} = μ_{1}/(μ_{1} - 1) (Fig. 6a). As both μ_{1} and μ_{2} are perturbed 50% up and down in our simulations, there is no point in Ω_{SSP} in which both perturbed values are bound to stay inside Ω_{SSP}.
Declarations
Acknowledgements
This work was supported in part by the National Programme for Research in Functional Genomics in Norway (FUGE) in the Research Council of Norway (grant no. NFR153302).
Authors’ Affiliations
References
- Glass L, Kauffman SA: Co-operative components, spatial localization and oscillatory cellular dynamics. Journal of Theoretical Biology. 1972, 34: 219-237. 10.1016/0022-5193(72)90157-9PubMedView ArticleGoogle Scholar
- Glass L, Kauffman SA: The logical analysis of continuous, non-linear biochemical control networks. Journal of Theoretical Biology. 1973, 39: 103-129. 10.1016/0022-5193(73)90208-7PubMedView ArticleGoogle Scholar
- Thomas R: Boolean Formalization of Genetic-Control Circuits. Journal of Theoretical Biology. 1973, 42 (3): 563-585. 10.1016/0022-5193(73)90247-6PubMedView ArticleGoogle Scholar
- Yuh CH, Bolouri H, Davidson EH: Genomic cis-Regulatory Logic: Experimental and Computational Analysis of a Sea Urchin Gene. Science. 1998, 279 (5358): 1896-1902. 10.1126/science.279.5358.1896PubMedView ArticleGoogle Scholar
- Gardner T, Cantor C, Collins J: Construction of a genetic toggle switch in Escherichia coli. Nature. 2000, 403 (6767): 339-342. 10.1038/35002131PubMedView ArticleGoogle Scholar
- Yuh CH, Bolouri H, Davidson EH: cis-regulatory logic in the endo16 gene: switching from a specification to a differentiation mode of control. Development. 2001, 128 (5): 617-629.PubMedGoogle Scholar
- Guet CC, Elowitz MB, Hsing WH, Leibler S: Combinatorial synthesis of genetic networks. Science. 2002, 296 (5572): 1466-1470. 10.1126/science.1067407PubMedView ArticleGoogle Scholar
- Setty Y, Mayo AE, Surette MG, Alon U: Detailed map of a cis-regulatory input function. PNAS. 2003, 100 (13): 7702-7707. 10.1073/pnas.1230759100PubMed CentralPubMedView ArticleGoogle Scholar
- Kramer BP, Fischer C, Fussenegger M: BioLogic gates enable logical transcription control in mammalian cells. Biotechnology and Bioengineering. 2004, 87 (4): 478-484. 10.1002/bit.20142PubMedView ArticleGoogle Scholar
- Istrail S, Davidson EH: Logic functions of the genomic cis-regulatory code. PNAS. 2005, 102 (14): 4954-4959. 10.1073/pnas.0409624102PubMed CentralPubMedView ArticleGoogle Scholar
- Mayo AE, Setty Y, Shavit S, Zaslaver A, Alon U: Plasticity of the cis-Regulatory Input Function of a Gene. PLoS Biology. 2006, 4 (4): e45- 10.1371/journal.pbio.0040045PubMed CentralPubMedView ArticleGoogle Scholar
- Goldbeter A, Koshland DE: An Amplified Sensitivity Arising from Covalent Modification in Biological Systems. PNAS. 1981, 78 (11): 6840-6844. 10.1073/pnas.78.11.6840PubMed CentralPubMedView ArticleGoogle Scholar
- Rossi FMV, Kringstein AM, Spicher A, Guicherit OM, Blau HM: Transcriptional control: Rheostat converted to on/off switch. Molecular Cell. 2000, 6 (3): 723-728. 10.1016/S1097-2765(00)00070-8PubMedView ArticleGoogle Scholar
- Biggar SR, Crabtree GR: Cell signaling can direct either binary or graded transcriptional responses. The EMBO Journal. 2001, 20: 3167-3176. 10.1093/emboj/20.12.3167PubMed CentralPubMedView ArticleGoogle Scholar
- Rosenfeld N, Young JW, Alon U, Swain PS, Elowitz MB: Gene Regulation at the Single-Cell Level. Science. 2005, 307 (5717): 1962-1965. 10.1126/science.1106914PubMedView ArticleGoogle Scholar
- Kim J, Bates DG, Postlethwaite I, Ma L, Iglesias PA: Robustness analysis of biochemical network models. IEE Proceedings – Systems Biology. 2006, 153 (3): 96-104. 10.1049/ip-syb:20050024PubMedView ArticleGoogle Scholar
- Wolf DM, Eeckman FH: On the Relationship Between Genomic Regulatory Element Organization and Gene Regulatory Dynamics. Journal of Theoretical Biology. 1998, 195 (2): 167-186. 10.1006/jtbi.1998.0790PubMedView ArticleGoogle Scholar
- Buchler NE, Gerland U, Hwa T: Nonlinear protein degradation and the function of genetic circuits. PNAS. 2005, 102 (27): 9559-9564. 10.1073/pnas.0409553102PubMed CentralPubMedView ArticleGoogle Scholar
- Bintu L, Buchler NE, Garcia HG, Gerland U, Hwa T, Kondev J, Phillips R: Transcriptional regulation by the numbers: models. Current Opinion in Genetics & Development. 2005, 15 (2): 116-124. 10.1016/j.gde.2005.02.007View ArticleGoogle Scholar
- Bintu L, Buchler NE, Garcia HG, Gerland U, Hwa T, Kondev J, Kuhlman T, Phillips R: Transcriptional regulation by the numbers: applications. Current Opinion in Genetics & Development. 2005, 15 (2): 125-135. 10.1016/j.gde.2005.02.006View ArticleGoogle Scholar
- Thattai M, van Oudenaarden A: Attenuation of noise in ultrasensitive signaling cascades. Biophysical Journal. 2002, 82 (6): 2943-2950.PubMed CentralPubMedView ArticleGoogle Scholar
- Hooshangi S, Thiberge S, Weiss R: Ultrasensitivity and noise propagation in a synthetic transcriptional cascade. PNAS. 2005, 102 (10): 3581-3586. 10.1073/pnas.0408507102PubMed CentralPubMedView ArticleGoogle Scholar
- Plahte E, Mestl T, Omholt SW: Global Analysis of Steady Points for Systems of Differential-Equations with Sigmoid Interactions. Dynamics and Stability of Systems. 1994, 9 (4): 275-291.View ArticleGoogle Scholar
- Mestl T, Plahte E, Omholt S: A mathematical framework for describing and analyzing gene regulatory networks. Journal of Theoretical Biology. 1995, 176 (2): 291-300. 10.1006/jtbi.1995.0199PubMedView ArticleGoogle Scholar
- Plahte E, Mestl T, Omholt SW: A methodological basis for description and analysis of systems with complex switch-like interactions. Journal of Mathematical Biology. 1998, 36 (4): 321-348. 10.1007/s002850050103PubMedView ArticleGoogle Scholar
- Plahte E, Kjøglum S: Analysis and generic properties of gene regulatory networks with graded response functions. Physica D: Nonlinear Phenomena. 2005, 201 (1–2): 150-176. 10.1016/j.physd.2004.11.014.View ArticleGoogle Scholar
- Veflingstad SR, Plahte E: Analysis of gene regulatory network models with graded and binary transcriptional responses. Biosystems. 2007, 90: 323-339. 10.1016/j.biosystems.2006.09.036PubMedView ArticleGoogle Scholar
- Cherry J, Adler F: How to make a biological switch. Journal of Theoretical Biology. 2000, 203: 117-133. 10.1006/jtbi.2000.1068PubMedView ArticleGoogle Scholar
- Brazhnik P, de la Fuente A, Mendes P: Gene networks: how to put the function in genomics. Trends in Biotechnology. 2002, 20 (11): 467-472. 10.1016/S0167-7799(02)02053-XPubMedView ArticleGoogle Scholar
- Becskei A, Serrano L: Engineering stability in gene networks by autoregulation. Nature. 2000, 405 (6786): 590-593. 10.1038/35014651PubMedView ArticleGoogle Scholar
- Rosenfeld N, Elowitz MB, Alon U: Negative autoregulation speeds the response times of transcription networks. Journal of Molecular Biology. 2002, 323: 785-93. 10.1016/S0022-2836(02)00994-4PubMedView ArticleGoogle Scholar
- de Jong H: Modeling and Simulation of Genetic Regulatory Systems: A Literature Review. Journal of Computational Biology. 2002, 9: 67-104. 10.1089/10665270252833208PubMedView ArticleGoogle Scholar
- Thomas R, D'Ari R: Biological Feedback. 1990, Boca Raton, USA: CRC PressGoogle Scholar
- Alon U: An Introduction to Systems Biology. Design Principles of Biological Circuits. 2007, Boca Raton, London, New York: Chapman & Hall (CRC)Google Scholar
- Stelling J, Sauer U, Szallasi Z, Doyle FJ, Doyle J: Robustness of Cellular Functions. Cell. 2004, 118 (6): 675-685. 10.1016/j.cell.2004.09.008PubMedView ArticleGoogle Scholar
- Wagner A: Robustness and Evolvability in Living Systems. 2005, Princeton: Princeton University PressGoogle Scholar
- Kitano H: Biological robustness. Nature Reviews Genetics. 2004, 5 (11): 826-837. 10.1038/nrg1471PubMedView ArticleGoogle Scholar
- Venkatesh KV, Bhartiya S, Ruhela A: Multiple feedback loops are key to a robust dynamic performance of tryptophan regulation in Escherichia coli. FEBS Letters. 2004, 563 (1–3): 234-240. 10.1016/S0014-5793(04)00310-2PubMedView ArticleGoogle Scholar
- Bhartiya S, Chaudhary N, Venkatesh KV, Doyle FJ: Multiple feedback loop design in the tryptophan regulatory network of Escherichia coli suggests a paradigm for robust regulation of processes in series. Journal of the Royal Society Interface. 2006, 3 (8): 383-391. 10.1098/rsif.2005.0103PubMed CentralView ArticleGoogle Scholar
- Ma W, Lai L, Ouyang Q, Tang C: Robustness and modular design of the Drosophila segment polarity network. Molecular Systems Biology. 2006, 2:Google Scholar
- Siegal ML, Bergman A: Waddington's canalization revisited: Developmental stability and evolution. PNAS. 2002, 99 (16): 10528-10532. 10.1073/pnas.102303999PubMed CentralPubMedView ArticleGoogle Scholar
- von Dassow G, Meir E, Munro E, Odell G: The segment polarity network is a robust development module. Nature. 2000, 406 (6792): 188-192. 10.1038/35018085PubMedView ArticleGoogle Scholar
- Albert R, Jeong H, Barabasi AL: Error and attack tolerance of complex networks. Nature. 2000, 406 (6794): 378-382. 10.1038/35019019PubMedView ArticleGoogle Scholar
- Hartwell LH, Hopfield JJ, Leibler S, Murray AW: From molecular to modular cell biology. Nature. 1999, 402 (6761): C47-C52. 10.1038/35011540PubMedView ArticleGoogle Scholar
- Lauffenburger DA: Cell signaling pathways as control modules: Complexity for simplicity?. PNAS. 2000, 97 (10): 5031-5033. 10.1073/pnas.97.10.5031PubMed CentralPubMedView ArticleGoogle Scholar
- Kitano H: Towards a theory of biological robustness. Mol Syst Biol. 2007, 3: 137- 10.1038/msb4100179PubMed CentralPubMedView ArticleGoogle Scholar
- Thomas R, Kaufman M: Multistationarity, the basis of cell differentiation and memory. I. Structural conditions of multistationarity and other nontrivial behavior. Chaos. 2001, 11: 170-179. 10.1063/1.1350439PubMedView ArticleGoogle Scholar
- Thomas R: Circular causality. IEE Proceedings Systems Biology. 2006, 153 (4): 140-153. 10.1049/ip-syb:20050101PubMedView ArticleGoogle Scholar
- Gibson G, Dworkin I: Uncovering cryptic genetic variation. Nat Rev Genet. 2004, 5 (9): 681-90. 10.1038/nrg1426PubMedView ArticleGoogle Scholar
- Bergman A, Siegal ML: Evolutionary capacitance as a general feature of complex gene networks. Nature. 2003, 424 (6948): 549-552. 10.1038/nature01765PubMedView 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.