Skip to main content
  • Research article
  • Open access
  • Published:

Threshold-dominated regulation hides genetic variation in gene expression networks



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.


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.


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.


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 [411]. There is also extensive experimental documentation of sigmoidal or binary responses in gene regulation [5, 1116] as well as theoretical justifications based on classical principles from physical chemistry [1720].

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 [1820]. 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 [2327] 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

The above basic assumptions of Gene Regulatory Networks lead to the following generic model for the time course of gene product concentrations:

y ˙ j = κ j R j ( Z ) λ j Q j ( Z ) y j , MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGacaGaaiaabeqaaeqabiWaaaGcbaGafmyEaKNbaiaadaWgaaWcbaGaemOAaOgabeaakiabg2da9GGaciab=P7aRnaaBaaaleaacqWGQbGAaeqaaOGaemOuai1aaSbaaSqaaiabdQgaQbqabaGccqGGOaakcqWGAbGwcqGGPaqkcqGHsislcqWF7oaBdaWgaaWcbaGaemOAaOgabeaakiabdgfarnaaBaaaleaacqWGQbGAaeqaaOGaeiikaGIaemOwaOLaeiykaKIaemyEaK3aaSbaaSqaaiabdQgaQbqabaGccqGGSaalaaa@4703@

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].

In almost all cases, regulation of the degradation is disregarded, thus we assume Q j = 1. We let S be a Hill function, Z j = S(y j , θ j , p) = y j p/(y j p+ θ j p), with the property that when p → ∞, then Z j approaches the Heaviside step function with threshold θ j , and put all p j = p (Fig. 1). Of course, in real systems sigmoids most likely do not have the same steepness. Our justification for taking all p j equal is to investigate systematically whether there is a threshold robustness effect for varying steepnesses. Once this effect is established, one might take one step further to investigate necessary and sufficient conditions on the steepnesses for robustness in more realistic models.

Figure 1
figure 1

Regulatory functions used in the simulations. (a) The Hill function describes the dose-responce relationship between the amount of a regulator x1 and the relative production rate of the regulated gene. The threshold parameter θ1 = 2 is the same for all four curves. (b) Regulatory function for two regulators x1 and x2. Two Hill functions are combined by the algebraic equivalent of the Boolean AND function (see Table 1). The parameter values used in this panel are θ1 = θ2 = 2 and p = 10.

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.

An illustrative example is provided by a popular model for a negatively autoregulated gene, which are very common in e.g. E. coli

m ˙ = r + α ( 1 Z ) γ m , y ˙ = κ m λ y , MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGacaGaaiaabeqaaeqabiWaaaGcbaqbaeWabiqaaaqaaiqbd2gaTzaacaGaeyypa0JaemOCaiNaey4kaSccciGae8xSdeMaeiikaGIaeGymaeJaeyOeI0IaemOwaOLaeiykaKIaeyOeI0Iae83SdCMaemyBa0MaeiilaWcabaGafmyEaKNbaiaacqGH9aqpcqWF6oWAcqWGTbqBcqGHsislcqWF7oaBcqWG5bqEcqGGSaalaaaaaa@46C6@

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 y ˙ MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGacaGaaiaabeqaaeqabiWaaaGcbaGafmyEaKNbaiaaaaa@2D57@ = (κ/γ)[r + α(1 - Z)] - λy, which is of the generic type Eq. (1).

After elimination of m, the equilibrium equation for y can be solved graphically (Fig. 2). There are three qualitatively different situations to consider: the solution can be close to the threshold and approach the threshold as steepness increases (red line), or the gene can be almost off (green line) or almost constitutively on (magenta line), the differences being accentuated for higher steepness. In the first case the stationary point is a SSP, in the two other cases it is regular. In the singular case we find the solution in the step function limit by putting y* = θ and solve for m*, getting m* = λθ/κ, Z* = 1 + r/α - g, where g = γλθ/(ακ). From the requirement Z* 0, 1 it follows that this solution exists when r/α <g < 1 + r/α. Thus, when g lies in this interval, the protein concentration y* is very close to the threshold, independent of parameter values, and y* = θ in the step function limit.

Figure 2
figure 2

Graphical solution of the equilibrium condition of Eq. (2) for varying steepness of the response function and varying relative degradation rates. The blue curves are graphs of r + α(1 - Z), the straight lines graph the degradation term γ y. For intermediate values of γ the equilibrium concentration is close to θ and approaches θ when the response function steepness increases (red line), and there is active regulation. When γ is small, the basal transcription rate is sufficient to balance the degradation, and y* gets large (green line). When γ is large, degradation is so rapid that the protein concentration never reaches the level where it regulates. Maximal production is necessary to balance the degradation (magenta line).

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 (Z1, Z2) - γ j x j ,   j = 1,...,3,

where x j is the scaled version of y j , α3 0, 1 and all γ j 0, 1, and Z j = x j p/(x j p+ 1). The prime denotes differentiation with respect to a scaled time. The Boolean functions R j are chosen so that the system has a unique stable point in which both x1 and x2 are singular (close to threshold at equilibrium). Altogether 14 models satsify this requirement (Fig. 3, see also the Methods section).

Figure 3
figure 3

Connectivity diagrams for the 14 network models in the simulation study. Genes 1 and 2 are represented by circles, the downstream gene 3 being omitted for clarity. The sign of an arrow indicates whether the type of regulation is activation (+), in which case the input variable is Z i , or inhibition (-), in which case the input variable is 1 - Z i . When a gene has two regulators, the individual signals are combined with a logic block, represented by a rectangle, merging the two signals into one by the continuous analogue of the Boolean functions AND or OR. (See the Methods section for explanations of the Boolean variables and functions.)

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 ( p ) MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGacaGaaiaabeqaaeqabiWaaaGcbaGaemiEaG3aa0baaSqaaiabdQgaQbqaaiabgEHiQaaakiabcIcaOiabdchaWjabcMcaPaaa@32EA@ 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 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGacaGaaiaabeqaaeqabiWaaaGcbaGaemiEaG3aa0baaSqaaiabdMgaPbqaaiabgEHiQaaaaaa@2FC3@ for each of the 14 models we generated 81 parameter sets, giving a total of 1134 different systems with a unique SSP for which x1 and x2 are singular and x3 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 CVuni ≈ 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.


Models 1, 9, 12, and 14 represent four different classes among the 14 models: Model 1 represents models with a negative feedback loop between the two genes plus autoregulation, Model 9 has just a pure negative feedback loop and no autoregulation, Model 12 has interaction but no feedback loop between genes 1 and 2, and Model 14 has no interaction at all between genes 1 and 2 (Fig. 3). The comparison (Fig. 4) of the distributions of CV1kand CV2kover the 81 parameter sets to the distribution of CV3kshows that there is a marked difference between the CVs of the singular variables x1 and x2 and the downstream, regular variable x3. This difference is most marked for high Hill coefficient, but is present even under Michaelis-Menten conditions. While in almost all cases the variation in x 3 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGacaGaaiaabeqaaeqabiWaaaGcbaGaemiEaG3aa0baaSqaaiabiodaZaqaaiabgEHiQaaaaaa@2F5C@ is larger than CVuni, the variations in x 1 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGacaGaaiaabeqaaeqabiWaaaGcbaGaemiEaG3aa0baaSqaaiabigdaXaqaaiabgEHiQaaaaaa@2F58@ and x 2 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGacaGaaiaabeqaaeqabiWaaaGcbaGaemiEaG3aa0baaSqaaiabikdaYaqaaiabgEHiQaaaaaa@2F5A@ are considerably smaller for most parameter sets, in Model 14 for all parameter sets. The result for Model 14 is not surprising, as it is well known that negative autoregulation leads to a high degree of robustness [30, 33, 34].

Figure 4
figure 4

Variation in steady state values for Models 1, 9, 12, and 14. Boxplots show the distributions of the 81 coefficients of variance for all three genes using seven different Hill coefficients ranging from p = 1 to p = 100. For each Hill coefficients the three plots show from left to right the coefficient of variation for gene 1 (red), 2 (green), and 3 (blue), respectively. The boxes show the quartiles and the median. The black vertical lines extend to the largest observation and the smallest. The long black horisontal line shows the coefficient of variation 0.288 of the perturbed production rates α j .

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 x3 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.

For all 14 models a decrease in p reduces robustness of both singular variables, defined as the minima of CV1kand CV2kover all 81 parameter sets (Fig. 5). But even at p = 1 there is less variation in the equilibrium values than in the perturbed parameter values. The cases of highest robustness for each of the 14 models and all steepness values show large variation among the models (Fig. 5). But for all models, even for the least robust Models 1 and 2, robustness increases with increasing Hill exponent, and is always smaller than CVuni. In Fig. 5c the differences among the models are accentuated for high Hill exponents. For each model the number CVminmax plotted along the vertical axis is computed as follows: for each of the 81 parameter sets first compute C V max k MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGacaGaaiaabeqaaeqabiWaaaGcbaGaem4qamKaemOvay1aa0baaSqaaiGbc2gaTjabcggaHjabcIha4bqaaiabdUgaRbaaaaa@33C9@ = max{CV1k, CV2k}, then find CVminmax = min{ C V max k MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGacaGaaiaabeqaaeqabiWaaaGcbaGaem4qamKaemOvay1aa0baaSqaaiGbc2gaTjabcggaHjabcIha4bqaaiabdUgaRbaaaaa@33C9@ }. For each model this procedure selects the parameter set for which both CV1 and CV2 are small, giving CVminmax as a measure of the highest robustness when both genes 1 and 2 are taken into account.

Figure 5
figure 5

The coefficient of variation as function of the Hill coefficient for the most robust parameter sets for gene 1 and 2 across all 14 models. (a)-(b): Minimum of CV1kand CV2k, respectively, over all 81 parameter sets. (c): CVminmax as function of the Hill coefficient for each of the 14 models. An explanation for the remarkably high values for Models 1, 2, and 11 for high Hill coefficient values is given in the text. (d) The ratio of minimum of CV2kto minimum of CV1kas function of the Hill coefficient for each of the 14 models.

Models 1, 2, and 11 have distinctly larger values of CVminmax than the other models. This behaviour can be explained as a consequence of the shape of the parameter space domain ΩSSP in (μ1, μ2)-space, where μ j = α j /γ j , j = 1, 2. This is the parameter domain for which both x 1 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGacaGaaiaabeqaaeqabiWaaaGcbaGaemiEaG3aa0baaSqaaiabigdaXaqaaiabgEHiQaaaaaa@2F58@ and x 2 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGacaGaaiaabeqaaeqabiWaaaGcbaGaemiEaG3aa0baaSqaaiabikdaYaqaaiabgEHiQaaaaaa@2F5A@ attain threshold values in the step function limit, in other words, the domain in which there is threshold robustness in both variables. For Models 1, 2, and 11 ΩSSP is concave, while for all the other models it is either wedge-shaped or rectangular (Fig. 6). In the concave domains of Models 1, 2, and 11 there is no point giving threshold robustness in both variables for the full parameter perturbation range of 50% up and down (see the Methods section, which also contains the derivation of this result), contrary to all the other models where threshold robustness is obtained for all parameter values sufficiently far from the boundary. With this result in mind it is reasonable to expect a drastically reduced robustness for these three models compared to the rest.

Figure 6
figure 6

The shaded areas are the robustness domains ΩSSP in the (μ1, μ2)-plane for (a) Models 1, 2, 11, (b) Models 3, 13, (c) Model 6, (d) Models 4, 5, 7–10, 12, 14. For parameter values in ΩSSP both x 1 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGacaGaaiaabeqaaeqabiWaaaGcbaGaemiEaG3aa0baaSqaaiabigdaXaqaaiabgEHiQaaaaaa@2F58@ and x 2 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGacaGaaiaabeqaaeqabiWaaaGcbaGaemiEaG3aa0baaSqaaiabikdaYaqaaiabgEHiQaaaaaa@2F5A@ are singular variables and approach their threshold values in the step function limit. Then they exhibit threshold robustness for all parameter perturbations which leave the perturbed values inside ΩSSP.

A comparison of CV1 and CV2 for all models show a distinct difference in robustness of x1 and x2 for Model 2 and Model 11 (Fig. 5d). We can explain this difference by how the character of the stationary point varies over the parameter space (Fig. 7 for Model 2). One can see that x1 is singular for all μ2 > 1 independent of μ1, while the domain in which x2 is singular is much smaller and with a strongly narrowing band. In this band, all points are close to the boundary, and robustness in x2 is very easily lost. Accordingly, the probability of having a perturbed point in parameter space in which the singular state is preserved is much less for x2 than for x1, just as seen in Fig. 5. For Model 11 the situation is similar.

Figure 7
figure 7

In the case of very steep sigmoids, the (μ1, μ2) space of Model 2 is divided into 5 domains, each domain comprising the parameter values giving a particular type of SSP. For example, in the domain denoted (Z1, 1), x1 is at its threshold and is singular, thus Z1 ≠ 0, 1, while x2 is above its threshold, thus Z2 = 1. Only in the shaded domain are both variables singular and actively regulating. For steep, but not infinitely steep response functions the relations are approximately true.

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 x1 and x2 as active regulators. We consider x j an active regulator if Z j MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGacaGaaiaabeqaaeqabiWaaaGcbaGaemOwaO1aa0baaSqaaiabdQgaQbqaaiabgEHiQaaaaaa@2F89@ 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 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGacaGaaiaabeqaaeqabiWaaaGcbaGaemiEaG3aa0baaSqaaiabigdaXaqaaiabgEHiQaaaaaa@2F58@ and x 2 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGacaGaaiaabeqaaeqabiWaaaGcbaGaemiEaG3aa0baaSqaaiabikdaYaqaaiabgEHiQaaaaaa@2F5A@ .


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 [3941] 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 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGacaGaaiaabeqaaeqabiWaaaGcbaGaemiEaG3aa0baaSqaaiabdQgaQbqaaiabgEHiQaaaaaa@2FC5@ 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 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGacaGaaiaabeqaaeqabiWaaaGcbaGaemOwaO1aa0baaSqaaiabdQgaQbqaaiabgEHiQaaaaaa@2F89@ 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 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGacaGaaiaabeqaaeqabiWaaaGcbaGaemiEaG3aa0baaSqaaiabdQgaQbqaaiabgEHiQaaaaaa@2FC5@ 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 Rn, F is a differentiable function F : RnRn, and F(x0) = g0, can only have a differentiable solution x = G(g), x0 = G(g0), 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].


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.


Model equations and regulatory functions

To simplify the simulations our model system

y ˙ j = κ j R j ( Z 1 , Z 2 ) λ j y j , j = 1 , ... , 3 , MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGacaGaaiaabeqaaeqabiWaaaGcbaqbaeqabeGaaaqaaiqbdMha5zaacaWaaSbaaSqaaiabdQgaQbqabaGccqGH9aqpiiGacqWF6oWAdaWgaaWcbaGaemOAaOgabeaakiabdkfasnaaBaaaleaacqWGQbGAaeqaaOGaeiikaGIaemOwaO1aaSbaaSqaaiabigdaXaqabaGccqGGSaalcqWGAbGwdaWgaaWcbaGaeGOmaidabeaakiabcMcaPiabgkHiTiab=T7aSnaaBaaaleaacqWGQbGAaeqaaOGaemyEaK3aaSbaaSqaaiabdQgaQbqabaGccqGGSaalaeaacqWGQbGAcqGH9aqpcqaIXaqmcqGGSaalcqGGUaGlcqGGUaGlcqGGUaGlcqGGSaalcqaIZaWmcqGGSaalaaaaaa@4F61@

was scaled to give the scaled equations (3)

x j ' = α j R j (Z1, Z2) - γ 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 regulatory functions R1 and R2 were chosen as algebraic equivalents of Boolean functions with two inputs, subject to the requirements ∂R j /∂Z j ≤ 0, j = 1, 2, and (∂R1/∂Z2)(∂R2/∂Z1) ≤ 0 to ensure that there is no positive loop in the system and accordingly no multistationarity, and a globally attracting stationary point. Excluding TRUE and FALSE which imply no regulation, and two more functions that do not fulfil the above mono-stationarity conditions, we were left with the 12 regulatory functions listed in Table 1. The algebraic equivalents are computed by the rules given in [25], X1 and X2 are Boolean variables, X ¯ MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGacaGaaiaabeqaaeqabiWaaaGcbaGafeiwaGLbaebaaaa@2D22@ i = NOT X i ; Z1 and Z2 are the corresponding sigmoidal functions as defined above. The third function R3(Z1, Z2) = Z1Z2 in all cases.

Table 1 The 12 Boolean functions used and their algebraic equivalents

Disregarding cases that were equivalent due to the x1-x2 symmetry of Eq. (3), we were left with 14 different models (Fig. 3, only x1 and x2 are shown). For example, the scaled equation system for Model 6 is

x 1 = α 1 ( 1 Z 1 ) Z 2 γ 1 x 1 , x 2 = α 2 ( 1 Z 1 ) γ 2 x 2 , x 3 = α 3 Z 1 Z 2 γ 3 x 3 . MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGacaGaaiaabeqaaeqabiWaaaGcbaqbaeWabmqaaaqaaiqbdIha4zaafaWaaSbaaSqaaiabigdaXaqabaGccqGH9aqpiiGacqWFXoqydaWgaaWcbaGaeGymaedabeaakiabcIcaOiabigdaXiabgkHiTiabdQfaAnaaBaaaleaacqaIXaqmaeqaaOGaeiykaKIaemOwaO1aaSbaaSqaaiabikdaYaqabaGccqGHsislcqWFZoWzdaWgaaWcbaGaeGymaedabeaakiabdIha4naaBaaaleaacqaIXaqmaeqaaOGaeiilaWcabaGafmiEaGNbauaadaWgaaWcbaGaeGOmaidabeaakiabg2da9iab=f7aHnaaBaaaleaacqaIYaGmaeqaaOGaeiikaGIaeGymaeJaeyOeI0IaemOwaO1aaSbaaSqaaiabigdaXaqabaGccqGGPaqkcqGHsislcqWFZoWzdaWgaaWcbaGaeGOmaidabeaakiabdIha4naaBaaaleaacqaIYaGmaeqaaOGaeiilaWcabaGafmiEaGNbauaadaWgaaWcbaGaeG4mamdabeaakiabg2da9iab=f7aHnaaBaaaleaacqaIZaWmaeqaaOGaemOwaO1aaSbaaSqaaiabigdaXaqabaGccqWGAbGwdaWgaaWcbaGaeGOmaidabeaakiabgkHiTiab=n7aNnaaBaaaleaacqaIZaWmaeqaaOGaemiEaG3aaSbaaSqaaiabiodaZaqabaGccqGGUaGlaaaaaa@685A@


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 x3 = (K/κ3)y3.

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.


This subsection describes the generation of parameter sets for Eq. (3). The parameters γ1, γ2, γ3 were sampled from independent uniform distributions over 0, 1. Since we were interested in parameter sets that exhibit threshold regulation in the step-function limit, we did not sample the production rates directly but rather worked our way back from the values of the equilibrium values Z 1 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGacaGaaiaabeqaaeqabiWaaaGcbaGaemOwaO1aa0baaSqaaiabigdaXaqaaiabgEHiQaaaaaa@2F1C@ and Z 2 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGacaGaaiaabeqaaeqabiWaaaGcbaGaemOwaO1aa0baaSqaaiabikdaYaqaaiabgEHiQaaaaaa@2F1E@ in the limit p → ∞. Granted that these lie in 0, 1, the steady state values of x1 and x2 are at their thresholds θ j = 1. From the steady state conditions of Eq. (3), the maximal production rates for gene 1 and 2 are then given by

α j = γ j R j ( Z 1 , Z 2 ) , j = 1 , 2. MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGacaGaaiaabeqaaeqabiWaaaGcbaqbaeqabeGaaaqaaGGaciab=f7aHnaaBaaaleaacqWGQbGAaeqaaOGaeyypa0tcfa4aaSaaaeaacqWFZoWzdaWgaaqaaiabdQgaQbqabaaabaGaemOuai1aaSbaaeaacqWGQbGAaeqaaiabcIcaOiabdQfaAnaaDaaabaGaeGymaedabaGaey4fIOcaaiabcYcaSiabdQfaAnaaDaaabaGaeGOmaidabaGaey4fIOcaaiabcMcaPaaakiabcYcaSaqaaiabdQgaQjabg2da9iabigdaXiabcYcaSiabikdaYiabc6caUaaaaaa@46D9@

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 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGacaGaaiaabeqaaeqabiWaaaGcbaGaemOwaO1aa0baaSqaaiabigdaXaqaaiabgEHiQaaaaaa@2F1C@ , Z 2 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGacaGaaiaabeqaaeqabiWaaaGcbaGaemOwaO1aa0baaSqaaiabikdaYaqaaiabgEHiQaaaaaa@2F1E@ ) 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 φkwe generated a set of 50 perturbed parameter values φ ˜ k MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGacaGaaiaabeqaaeqabiWaaaGcbaacciGaf8NXdyMbaGaadaqhaaWcbaGaeS4eHWgabaGaem4AaSgaaaaa@305F@ , ℓ = 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 k MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGacaGaaiaabeqaaeqabiWaaaGcbaGaemiEaG3aa0baaSqaaiabdQgaQjabloriSbqaaiabgEHiQiabdUgaRbaaaaa@3255@ of the three gene products in Eq. (3).

Measuring robustness towards pertubation

The robustness properties of the the three gene products in Eq. (3) were evaluated by comparing the coefficient of variation (CV) for parameters to the coefficient of variation for the steady state values. For each k {1,...,81} the robustness of the steady state values of x j against perturbations in the production rates α j is analysed using the coefficient of variation

C V j k = std ( { x j k } = 1 , ... , 50 ) mean ( { x j k } = 1 , ... , 50 ) . MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGacaGaaiaabeqaaeqabiWaaaGcbaGaem4qamKaemOvay1aa0baaSqaaiabdQgaQbqaaiabdUgaRbaakiabg2da9KqbaoaalaaabaGaee4CamNaeeiDaqNaeeizaqMaeiikaGIaei4EaSNaemiEaG3aa0baaeaacqWGQbGAcqWItecBaeaacqGHxiIkcqWGRbWAaaGaeiyFa03aaSbaaeaacqWItecBcqGH9aqpcqaIXaqmcqGGSaalcqGGUaGlcqGGUaGlcqGGUaGlcqGGSaalcqaI1aqncqaIWaamaeqaaiabcMcaPaqaaiabb2gaTjabbwgaLjabbggaHjabb6gaUjabcIcaOiabcUha7jabdIha4naaDaaabaGaemOAaOMaeS4eHWgabaGaey4fIOIaem4AaSgaaiabc2ha9naaBaaabaGaeS4eHWMaeyypa0JaeGymaeJaeiilaWIaeiOla4IaeiOla4IaeiOla4IaeiilaWIaeGynauJaeGimaadabeaacqGGPaqkaaGaeiOla4caaa@66EF@

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 CVuni = 1/ 12 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGacaGaaiaabeqaaeqabiWaaaGcbaWaaOaaaeaacqaIXaqmcqaIYaGmaSqabaaaaa@2DD0@ ≈ 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 x1 and x2 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.

From Fig. 3 it follows that the scaled equations of motion for x1 and x2 in Model 1 are

x 1 = α 1 [ 1 Z 1 Z 2 ] γ 1 x 1 , x 2 = α 2 [ 1 Z 2 + Z 1 Z 2 ] γ 2 x 2 . MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGacaGaaiaabeqaaeqabiWaaaGcbaqbaeWabiqaaaqaaiqbdIha4zaafaWaaSbaaSqaaiabigdaXaqabaGccqGH9aqpiiGacqWFXoqydaWgaaWcbaGaeGymaedabeaakiabcUfaBjabigdaXiabgkHiTiabdQfaAnaaBaaaleaacqaIXaqmaeqaaOGaemOwaO1aaSbaaSqaaiabikdaYaqabaGccqGGDbqxcqGHsislcqWFZoWzdaWgaaWcbaGaeGymaedabeaakiabdIha4naaBaaaleaacqaIXaqmaeqaaOGaeiilaWcabaGafmiEaGNbauaadaWgaaWcbaGaeGOmaidabeaakiabg2da9iab=f7aHnaaBaaaleaacqaIYaGmaeqaaOGaei4waSLaeGymaeJaeyOeI0IaemOwaO1aaSbaaSqaaiabikdaYaqabaGccqGHRaWkcqWGAbGwdaWgaaWcbaGaeGymaedabeaakiabdQfaAnaaBaaaleaacqaIYaGmaeqaaOGaeiyxa0LaeyOeI0Iae83SdC2aaSbaaSqaaiabikdaYaqabaGccqWG4baEdaWgaaWcbaGaeGOmaidabeaakiabc6caUaaaaaa@5D23@

Assume the stationary point is a SSP in which both variables are singular, and that p is very large. Then at equilibrium, x 1 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGacaGaaiaabeqaaeqabiWaaaGcbaGaemiEaG3aa0baaSqaaiabigdaXaqaaiabgEHiQaaaaaa@2F58@ θ1 = 1, x 2 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGacaGaaiaabeqaaeqabiWaaaGcbaGaemiEaG3aa0baaSqaaiabikdaYaqaaiabgEHiQaaaaaa@2F5A@ θ2 = 1, and 0 < Z i MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGacaGaaiaabeqaaeqabiWaaaGcbaGaemOwaO1aa0baaSqaaiabdMgaPbqaaiabgEHiQaaaaaa@2F87@ < 1. Solving the equilibrium conditions with respect to μ j , j = 1, 2 under these assumptions we find that μ j > 1 and

μ 1 = 1 1 Z 1 Z 2 , μ 2 = 1 1 Z 2 + Z 1 Z 2 . MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGacaGaaiaabeqaaeqabiWaaaGcbaqbaeWabiqaaaqaaGGaciab=X7aTnaaBaaaleaacqaIXaqmaeqaaOGaeyypa0tcfa4aaSaaaeaacqaIXaqmaeaacqaIXaqmcqGHsislcqWGAbGwdaqhaaqaaiabigdaXaqaaiabgEHiQaaacqWGAbGwdaqhaaqaaiabikdaYaqaaiabgEHiQaaaaaGccqGGSaalaeaacqWF8oqBdaWgaaWcbaGaeGOmaidabeaakiabg2da9KqbaoaalaaabaGaeGymaedabaGaeGymaeJaeyOeI0IaemOwaO1aa0baaeaacqaIYaGmaeaacqGHxiIkaaGaey4kaSIaemOwaO1aa0baaeaacqaIXaqmaeaacqGHxiIkaaGaemOwaO1aa0baaeaacqaIYaGmaeaacqGHxiIkaaaaaOGaeiOla4caaaaa@4DC4@

Upon elimination of Z 1 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGacaGaaiaabeqaaeqabiWaaaGcbaGaemOwaO1aa0baaSqaaiabigdaXaqaaiabgEHiQaaaaaa@2F1C@ this leads to

μ 2 = μ 1 ( 2 Z 2 ) μ 1 1 . MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGacaGaaiaabeqaaeqabiWaaaGcbaacciGae8hVd02aaSbaaSqaaiabikdaYaqabaGccqGH9aqpjuaGdaWcaaqaaiab=X7aTnaaBaaabaGaeGymaedabeaaaeaacqGGOaakcqaIYaGmcqGHsislcqWGAbGwdaqhaaqaaiabikdaYaqaaiabgEHiQaaacqGGPaqkcqWF8oqBdaWgaaqaaiabigdaXaqabaGaeyOeI0IaeGymaedaaOGaeiOla4caaa@3FCA@

Then ΩSSP is the domain covered by this family of curves when Z 2 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGacaGaaiaabeqaaeqabiWaaaGcbaGaemOwaO1aa0baaSqaaiabikdaYaqaaiabgEHiQaaaaaa@2F1E@ 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.


  1. 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-9

    Article  CAS  PubMed  Google Scholar 

  2. 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-7

    Article  CAS  PubMed  Google Scholar 

  3. Thomas R: Boolean Formalization of Genetic-Control Circuits. Journal of Theoretical Biology. 1973, 42 (3): 563-585. 10.1016/0022-5193(73)90247-6

    Article  CAS  PubMed  Google Scholar 

  4. 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.1896

    Article  CAS  PubMed  Google Scholar 

  5. Gardner T, Cantor C, Collins J: Construction of a genetic toggle switch in Escherichia coli. Nature. 2000, 403 (6767): 339-342. 10.1038/35002131

    Article  CAS  PubMed  Google Scholar 

  6. 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.

    CAS  PubMed  Google Scholar 

  7. Guet CC, Elowitz MB, Hsing WH, Leibler S: Combinatorial synthesis of genetic networks. Science. 2002, 296 (5572): 1466-1470. 10.1126/science.1067407

    Article  CAS  PubMed  Google Scholar 

  8. 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.1230759100

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  9. 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.20142

    Article  CAS  PubMed  Google Scholar 

  10. Istrail S, Davidson EH: Logic functions of the genomic cis-regulatory code. PNAS. 2005, 102 (14): 4954-4959. 10.1073/pnas.0409624102

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  11. 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.0040045

    Article  PubMed Central  PubMed  Google Scholar 

  12. 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.6840

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  13. 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-8

    Article  CAS  PubMed  Google Scholar 

  14. 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.3167

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  15. 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.1106914

    Article  CAS  PubMed  Google Scholar 

  16. 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:20050024

    Article  CAS  PubMed  Google Scholar 

  17. 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.0790

    Article  CAS  PubMed  Google Scholar 

  18. Buchler NE, Gerland U, Hwa T: Nonlinear protein degradation and the function of genetic circuits. PNAS. 2005, 102 (27): 9559-9564. 10.1073/pnas.0409553102

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  19. 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.007

    Article  CAS  Google Scholar 

  20. 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.006

    Article  CAS  Google Scholar 

  21. Thattai M, van Oudenaarden A: Attenuation of noise in ultrasensitive signaling cascades. Biophysical Journal. 2002, 82 (6): 2943-2950.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  22. Hooshangi S, Thiberge S, Weiss R: Ultrasensitivity and noise propagation in a synthetic transcriptional cascade. PNAS. 2005, 102 (10): 3581-3586. 10.1073/pnas.0408507102

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  23. 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.

    Article  Google Scholar 

  24. 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.0199

    Article  CAS  PubMed  Google Scholar 

  25. 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/s002850050103

    Article  CAS  PubMed  Google Scholar 

  26. 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.

    Article  CAS  Google Scholar 

  27. 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.036

    Article  CAS  PubMed  Google Scholar 

  28. Cherry J, Adler F: How to make a biological switch. Journal of Theoretical Biology. 2000, 203: 117-133. 10.1006/jtbi.2000.1068

    Article  CAS  PubMed  Google Scholar 

  29. 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-X

    Article  CAS  PubMed  Google Scholar 

  30. Becskei A, Serrano L: Engineering stability in gene networks by autoregulation. Nature. 2000, 405 (6786): 590-593. 10.1038/35014651

    Article  CAS  PubMed  Google Scholar 

  31. 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-4

    Article  CAS  PubMed  Google Scholar 

  32. de Jong H: Modeling and Simulation of Genetic Regulatory Systems: A Literature Review. Journal of Computational Biology. 2002, 9: 67-104. 10.1089/10665270252833208

    Article  CAS  PubMed  Google Scholar 

  33. Thomas R, D'Ari R: Biological Feedback. 1990, Boca Raton, USA: CRC Press

    Google Scholar 

  34. Alon U: An Introduction to Systems Biology. Design Principles of Biological Circuits. 2007, Boca Raton, London, New York: Chapman & Hall (CRC)

    Google Scholar 

  35. 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.008

    Article  CAS  PubMed  Google Scholar 

  36. Wagner A: Robustness and Evolvability in Living Systems. 2005, Princeton: Princeton University Press

    Google Scholar 

  37. Kitano H: Biological robustness. Nature Reviews Genetics. 2004, 5 (11): 826-837. 10.1038/nrg1471

    Article  CAS  PubMed  Google Scholar 

  38. 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-2

    Article  CAS  PubMed  Google Scholar 

  39. 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.0103

    Article  CAS  PubMed Central  Google Scholar 

  40. 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 

  41. Siegal ML, Bergman A: Waddington's canalization revisited: Developmental stability and evolution. PNAS. 2002, 99 (16): 10528-10532. 10.1073/pnas.102303999

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  42. 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/35018085

    Article  CAS  PubMed  Google Scholar 

  43. Albert R, Jeong H, Barabasi AL: Error and attack tolerance of complex networks. Nature. 2000, 406 (6794): 378-382. 10.1038/35019019

    Article  CAS  PubMed  Google Scholar 

  44. Hartwell LH, Hopfield JJ, Leibler S, Murray AW: From molecular to modular cell biology. Nature. 1999, 402 (6761): C47-C52. 10.1038/35011540

    Article  CAS  PubMed  Google Scholar 

  45. Lauffenburger DA: Cell signaling pathways as control modules: Complexity for simplicity?. PNAS. 2000, 97 (10): 5031-5033. 10.1073/pnas.97.10.5031

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  46. Kitano H: Towards a theory of biological robustness. Mol Syst Biol. 2007, 3: 137- 10.1038/msb4100179

    Article  PubMed Central  PubMed  Google Scholar 

  47. 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.1350439

    Article  PubMed  Google Scholar 

  48. Thomas R: Circular causality. IEE Proceedings Systems Biology. 2006, 153 (4): 140-153. 10.1049/ip-syb:20050101

    Article  CAS  PubMed  Google Scholar 

  49. Gibson G, Dworkin I: Uncovering cryptic genetic variation. Nat Rev Genet. 2004, 5 (9): 681-90. 10.1038/nrg1426

    Article  CAS  PubMed  Google Scholar 

  50. Bergman A, Siegal ML: Evolutionary capacitance as a general feature of complex gene networks. Nature. 2003, 424 (6948): 549-552. 10.1038/nature01765

    Article  CAS  PubMed  Google Scholar 

Download references


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).

Author information

Authors and Affiliations


Corresponding author

Correspondence to Arne B Gjuvsland.

Additional information

Authors' contributions

SWO and EP initiated the work, ABG and EP devised the experimental design, ABG ran all the simulations, carried out the data analysis and designed the figures, all authors contributed to and approved the final manuscript.

Authors’ original submitted files for images

Rights and permissions

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 (, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions

About this article

Cite this article

Gjuvsland, A.B., Plahte, E. & Omholt, S.W. Threshold-dominated regulation hides genetic variation in gene expression networks. BMC Syst Biol 1, 57 (2007).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: