Skip to main content

Proportionality between variances in gene expression induced by noise and mutation: consequence of evolutionary robustness



Characterization of robustness and plasticity of phenotypes is a basic issue in evolutionary and developmental biology. The robustness and plasticity are concerned with changeability of a biological system against external perturbations. The perturbations are either genetic, i.e., due to mutations in genes in the population, or epigenetic, i.e., due to noise during development or environmental variations. Thus, the variances of phenotypes due to genetic and epigenetic perturbations provide quantitative measures for such changeability during evolution and development, respectively.


Using numerical models simulating the evolutionary changes in the gene regulation network required to achieve a particular expression pattern, we first confirmed that gene expression dynamics robust to mutation evolved in the presence of a sufficient level of transcriptional noise. Under such conditions, the two types of variances in the gene expression levels, i.e. those due to mutations to the gene regulation network and those due to noise in gene expression dynamics were found to be proportional over a number of genes. The fraction of such genes with a common proportionality coefficient increased with an increase in the robustness of the evolved network. This proportionality was generally confirmed, also under the presence of environmental fluctuations and sexual recombination in diploids, and was explained from an evolutionary robustness hypothesis, in which an evolved robust system suppresses the so-called error catastrophe - the destabilization of the single-peaked distribution in gene expression levels. Experimental evidences for the proportionality of the variances over genes are also discussed.


The proportionality between the genetic and epigenetic variances of phenotypes implies the correlation between the robustness (or plasticity) against genetic changes and against noise in development, and also suggests that phenotypic traits that are more variable epigenetically have a higher evolutionary potential.


Plasticity and robustness are basic concepts in evolutionary and developmental biology. Plasticity refers to the changeability of phenotypes in response to external environmental perturbations. Indeed many important concepts in biology are concerned with the changeability in the system. This changeability depends on each phenotype: some phenotypes are more variable than others. How is such degree of changeability characterized quantitatively?

On the other hand, robustness is another basic concept in evolutionary and developmental biology. Here, phenotypic robustness is defined as the ability of the system to continue to function despite perturbations to it [17]. Phenotypes important for survival are expected to be robust, at least to some degree, to enable organisms to survive under such perturbations.

For both plasticity and robustness, there are epigenetic and genetic sources of perturbations to a biological system, which act in different time scales. Epigenetic perturbation works at a faster scale. Phenotypes are changed through noise in gene expression and developmental dynamics. Environmental variation gives another source for variability. Genetic variation, on the other hand, works at a longer time scale through evolution. Now, is there any relationship between the changes of genetic and epigenetic origins? If a phenotype is changeable easily epigenetically through development or environment, is it also more feasible to change genetically? Similarly, if a phenotype is robust to developmental perturbations, is it also robust to genetic variations through evolution? When we consider generally any dynamical systems, such relationship would not be expected. However, as a biological system is a result of evolution, existence of some relationship between the genetic and epigenetic robustness may be expected.

The relationship between evolution and robustness has been long debated since the pioneering studies by Schmalhausen and Waddington [8, 9]. Waddington, then coined the term "genetic assimilation", in which phenotypic changes induced environmentally are then assimilated to genetic changes through evolution. Although important, these pioneering studies have mostly emphasized the qualitative aspects of the relationship between robustness and evolution. However, advances in quantitative studies on cell biology have facilitated the quantitative assessment of this relationship. In particular, fluctuations of phenotypes (e.g., gene expression levels) that have been measured extensively through the fluorescent-protein techniques [1013] can provide a tool to explore such relationship.

In quantitative terms, robustness can be considered as a measure of the insensitivity of phenotypes to external disturbances and plasticity as a measure of changeability of phenotypes. On the other hand, fluctuation is the degree of phenotypic variation induced by perturbations. Hence, the phenotypic fluctuation (variance) increases with a decrease in robustness, and vice versa. Thus, the variance can serve as an index (inverse) for robustness, and also for plasticity. Now, the question concerning robustness and evolution can also be posed in terms of phenotypic variances. How is the evolution speed correlated with the variance? Does the variance increase or decrease through evolution?

Indeed, findings of previous studies involving artificial selection experiments with bacteria suggest that the rate of evolution, i.e., the increase in the fitness per generation, is proportional to the phenotypic variance among isogenic individuals [14, 15]. This relationship, originally defined on the basis of a macroscopic distribution theory, was also confirmed in in-silico experiments by using gene regulation networks (GRNs) and metabolic networks [16, 17].

This observed relationship is noteworthy because evolvability is characterized by non-genetic variation of phenotypes. Even if the rate of genetic change (mutation or recombination) is identical, the rate of evolution can differ according to this variation. To elucidate this point, recall again that there are 2 sources in phenotypic variances, genetic and epigenetic. Quantitatively, the former is characterized by the phenotypic variance in a heterogenic population and is due to genetic modifications, as, denoted as V g , whereas the latter, denoted here as V ip , is the phenotypic variance in an isogenic population due to noise during the developmental process. The former reflects the structural robustness of the phenotype, i.e., the rigidity of the phenotype against changes induced by genetic mutations, whereas the latter reflects the robustness of the phenotype against the stochasticity encountered during the developmental process or that induced by environmental changes. (Phenotypic variance of non-genetic origin is traditionally discussed as environmental variance V e . Here, we are concerned with the variance due to fluctuation during developmental process, and thus adopt this notation, but one could regard this variance as a component of V e [18, 19].)

It is obvious that evolution speed, which indicates the change in phenotype due to genetic changes, is correlated with V g , as was demonstrated by Fisher [18, 20]. Thus, the proportionality between the evolution speed and V ip suggests the proportionality between V g and V ip throughout the course of evolution. Indeed, this relationship was confirmed by the evolutionary stability theory and numerical simulations [14, 16], which imply that robustness to noise and to mutation are correlated.

So far the proportionality between V ip and V g of a given fitness through the course of evolution (over generations) was confirmed. Now, let us come back to the original question on a possible relationship between genetic and epigenetic robustness (or changeability) over many phenotypic traits. To discuss such problem, we need to study the relationship of the variances V ip and V g over many phenotypes or expressions of many genes, for a given individual (not over the evolutionary course). This is the focus of the present paper.

The phenotypes or gene expression levels are generally associated with several genes, as known as epistasis. Even if a fitness value may be directly related to a single phenotype or the expression of a single gene, many genes may modify the expression of each other. These expressions are interrelated through a complex gene regulation network; therefore, the nature of their correlation may change during the course of evolution. This may give some evolutionary constraint on the changes of expression levels of genes. Then, does a gene with a higher fluctuation in its expression level have a higher potentiality in evolution than others? Is there correlation between phenotypic changes by epigenetic noise in gene expression dynamics and by genetic variation? We will demonstrate such proportionality over expressions of many genes, rather than the previously studied proportionality of the variances through the evolutionary course.

Using a numerical model in GRNs and simulating their evolutionary changes required to increase a given level of fitness, we first found that the rate of evolution of the expression levels of several genes was highly correlated with (or roughly proportional to) the respective variances of these genes. Next, we proceeded to present evidence for the proportionality between the 2 types of variances of gene expressions, i.e., of genetic and epigenetic origins over many genes. This proportionality was achieved after a selection process under a given fitness condition and is true whenever the phenotype of evolved system is robust to transcriptional noise in gene expression dynamics and genetic mutation. Further, the generality of this relationship was confirmed by studying a variety of models and evolutionary stability theory of multivariate distribution. We will also discuss some experimental evidences for this relationship.


Modeling Strategy

Here, we considered the evolution of a simple model for "development." For this purpose, we postulated the following conditions for development and evolution.

(i) The set of variables x i (i = 1, , M ) represents the expression levels of M genes. These variables take continuous values, which we set such that if gene i is expressed, then x i > 0 and if not, x i < 0. (Choice of a threshold at zero is a matter of convenience. Indeed, we also carried out simulations of a model in which the expression x i takes only non-negative values with certain positive threshold values for activation. Results to be discussed were not altered by this change).

(ii) Gene expressions mutually activate or inhibit each other and are regulated by the GRN. The temporal evolution of the gene expression level x i is generally not simple because the value of M is large. The phenotype of each individual is defined by the set of gene expression level x i after evolution for "developmental time," which starts from the time point with a given initial gene expression level (set as x i < 0 for all i, i.e., none of the genes is expressed). The developmental time for the system is set to a large period to allow gene expression patterns to settle into an attractor state.

(iii) Gene expression dynamics are noisy. Owing to their dependence on chemical reactions involving molecular collisions, gene expression dynamics are also stochastic in nature. In particular, since the number of molecules (e.g., mRNA and proteins) involved in gene expression is not extremely large, a deviation from the average rate equation for the reaction is possible. This deviation is represented as the noise applied to the average gene regulation expression dynamics. The amplitude of noise is denoted as σ.

(iv) Genotype: Depending on the genotype, the structure of GRN changes, i.e., network paths that determine the genes responsible for activating (or repressing) a given gene. This interaction between genes is represented by a connection matrix J ij , which takes a value of 1 (if activating), -1 (if repressing), or 0 (if there is no influence).

(v) Population: A population of individuals with different genotypes, i.e., with each individual having a slightly different GRN-the matrix J ij . This gives a genotypic distribution.

(vi) Fitness: The pattern of expression in x i is determined on the basis of gene expression dynamics. Fitness is a function of the expression pattern of a subset of x i , i.e., the expression of a given set of "target" genes i = 1, , k. If the pattern of x i is closer to the prescribed pattern of genes, the fitness is higher. We take the prescribed pattern as "all on" for the target genes, unless otherwise mentioned. If all components of x i (i = 1, , k) are positive, the fitness is set at 0, and it is decreased by 1 when the number of negative x i is increased by 1.

(vii) Mutation and selection: Offspring are produced according to the fitness: individuals with low fitness cannot produce offspring. The selection process according to the fitness thus progresses, while mutation introduces a change in the GRN, represented as 1,-1, or 0, in a few elements of the matrix J ij . The fraction of the elements altered is given by the mutation rate.

We studied numerical models based on the abovementioned characteristics. (For details see Method). Note that due to the complex nature of gene expression dynamics, GRN with a higher fitness value is generally rare. Further, since the expression dynamics are subjected to perturbations by noise, the "goal" of reaching the highest fitness may not be easily achieved. However, through the evolutionary processes, including natural selection and mutations, a GRN with the highest fitness value is generated. Individuals with the highest fitness value in a population evolve over some generations to reach the highest fitness value, i.e., x i > 0, for the target genes i = 1, 2, ..., k.

We previously reported that the fitness distribution in a heterogenic population undergoes a transition when the noise level is increased [17, 21]. Before we discuss the main results of the present paper from the next subsection, we briefly summarize the earlier numerical result in this subsection.

When the noise level σ was below a given threshold (σ < σ c ), some individuals within the population might take lower fitness values, thereby inducing a considerable difference in the distribution of fitness over different genotypes. Some mutants derived from individuals with genotypes having high fitness values might take much lower fitness values, even after many generations of evolution. However, when the noise level was greater than the threshold (σ c < σ), the distribution of fitness was sharp, with concentration of mutants with high fitness values. Thus, low-fitness mutants were eliminated through the evolution. (If the noise level was too high, the expression levels were not fixed in time and increased or decreased over time. We did not examine such cases).

Accordingly, an increase in the noise strength lead to a transition in the robustness to mutation. GRNs that evolve under a low noise level did not have robustness to mutation, whereby some mutants could not sustain the high fitness value. In contrast, robustness to mutation was achieved for GRNs evolved only under a high noise level σ > σ.

This counterintuitive robustness to mutation for σ > σ c can be explained as follows. According to the dynamics of GRN that evolves under higher noise level, a large portion of the initial conditions reach target attractors that give the highest fitness values, thereby achieving robustness against noise, while for those evolved under σ < σ c , only a tiny fraction reaches target attractors. The developmental landscape for σ > σ c gives a global, smooth attraction to the target, whereas the landscape evolved at σ < σ c is rugged. Now, consider mutation to a network to slightly change gene expression dynamics. In the smooth landscape with global attraction, such a perturbation will cause little change to the final expression pattern, while under the dynamics with rugged developmental landscape, it often destroys the attraction to the target attractor. In other words, robustness to mutation evolves only under robustness to noise during development.

For σ > σ c , robustness to both noise and mutation evolved. We computed robustness to both factors using the 2 types of variances, V ip and V g , where the former was the noise-induced variance in the distribution of the fitness (i.e., the number of "on" target genes), in a population of isogenic individuals and V g is the variance in the distribution of fitness in a heterogenic population. The latter was obtained by first computing the mean fitness value for each genotype, and then by measuring the variance of these mean values for a heterogenic population. The 2 variances of the fitness decreased proportionally, throughout the course of evolution, for σ > σ c , in accordance with the principle of evolution of robustness (see [17, 21]).

V g - V ip relationship over genes

Apart from the fitness level, the expression level x i of each gene i was also distributed, even for isogenic distribution with the same J ij ; this was because of the stochasticity in each gene expression dynamics. Similar to the variances for the fitness, the phenotypic variance V ip (i) for each gene i in an isogenic population is defined on the basis of the variance of the expression of each gene i, with each X i = Sign(x i ), in an isogenic population. On the other hand, the mean expression level X i ¯ over the isogenic population depended on each genotype (i.e., the matrix J ij ). The variance computed using the distribution of X i ¯ in this heterogenic population, then, gives the genetic variance V g (i) for each gene i.

As mentioned above, our model also accounted for many non-target genes that do not contribute to fitness. The expression level x i of such non-target genes i could be either positive or negative because there was no selection pressure directed at fixing their expression level. However, we found that the expression levels of many non-target genes become fixed to positive or negative values over the course of evolution when σ > σ c . To achieve robust fitness, the expression levels of some of the non-target genes were also fixed to either the on or off status, with consistency across individuals in the gene expression status. Hence, the average expression level < X i ¯ > in a heterogenic population increased or decreased to either positive or negative values ( ... ¯ is the mean over an isogenic population, while < ... > is the average over a heterogenic population).

We computed the rate of evolution for each gene expression over a given number of generations, as the rate of either increase or decrease of the average expression level < X i ¯ > in a heterogenic population. It was computed after evolution for some generation to achieve an increase in fitness. We then examined the validity of the relationship between the evolution speed and the fluctuations over genes, by plotting | < Δ X i ¯ > | against V ip (i), where < Δ X i ¯ > is the change in the average expression level < X i ¯ > of a given gene over some generations (20, in Figure 1). The proportionality was found to be valid over different genes, both target and non-target, as shown in Figure 1. Genes with higher variances would have a higher potentiality to evolve (In Figure 1, the target genes had less variance and evolution speeds over these generations. This is because the fitness had already increased in this generation, fixing the expression levels to positive values, and there was little room for further increase).

Figure 1
figure 1

Relationship between evolution speed and isogenic phenotypic fluctuations. Following the Method, we computed the mean expression level X ( i ) ¯ = S i g n ( x i ) ¯ over 200 runs for each gene i. The ordinate shows the change of its average < | X ( i ) ¯ | > over all (= 200) individuals between the generations 80 and 90, representing the evolution speed of each gene expression i. The abscissa showed V ip (i), the variance of X(i) over 200 runs for the identical genotype (J ij ). Computed as the variance at generation 80. The evolution progressed at around these generations and was completed before 200 generations. The noise level was σ = 0.08. See Method for details.

Considering the possible generalization of Fisher's theorem [18, 20], which states that the evolution speed is proportional to V g , and applying it to the expression levels over genes, we may expect that the evolution speed of each gene expression level is proportional to V g (i). Then, the proportionality between V g (i) and V ip (i) over the genes can be expected from the proportionality between the evolution speed and V ip (i). Considering this, we examined the relationship between V g (i) and V ip (i) for several values of noise levels (Figure 2). From the plots, we obtained the following findings:

Figure 2
figure 2

Relationship between V g ( i ) and V ip ( i ). As described in the Method section of the text, V ip (i) was computed as the variance of the distribution of Sign(x i ) over 200 runs for an identical genotype, while V g (i) was computed as a variance of the distribution of ( S i g n ( x i ) ) ¯ over 200 individuals, where S i g n ( x i ) ¯ was the mean over 200 runs. σ = 0.1 (a), 0.05 (b), 0.01 (c), and 0.001 (d). The plot of (V g (i) and V ip (i)) for all genes i at the 300th generation: Target genes as red , and nontarget genes as blue +.

(i) Proportionality between V ip (i) and V g (i) was satisfied over many genes, i.e., r i V ip (i)/V g (i) took a common value ρ(< 1) for many genes i, when the evolution of robustness progressed at σ > σ c .

(ii) Target genes always lay on the proportional line, with relatively low values of V ip (i) (and accordingly, V g (i)), while the variances of many non-target genes also lay on the same proportionality line r i ~ ρ. Variances of only a few non-target genes did not exhibit the abovementioned proportional relationship, and the ratios r i for such genes were scattered between ρ and 1. (See also Additional file 1, Figure S1).

(iii) The fraction of such genes that fitted on the single proportional line increased with the noise strength σ. As the noise level was lowered, an increase was noted in the fraction of genes showing expression variances that deviate from the abovementioned proportional relationship. At around the threshold noise level σ c , most genes approached r i ~ 1. (See Additional file 1, Figure S1, for noise dependence of the variances V ip (i) and V g (i)).

We then plotted the histogram of r i over all genes i, sampled over a few sets of generations (see Figure 3a, plotted in log-scale for r i ). The figure showed that the peak at ρ was more prominent with an increase in σ. On the other hand, the broader distribution ranging between r i ~ ρ and 1 became more prominent as the noise level decreased, until the distribution around r i ~ 1 dominated at σ ~ σ c .

Figure 3
figure 3

Histogram of the ratio V g ( i ) /V ip ( i ). (a)The variances V g (i) and V ip (i) and their ratio were computed for each gene i, in the same manner as that adopted in Figure 2. Histogram of the l o g 10 V g ( i ) V i p ( i ) was computed and sampled for all genes i from 280 to 300 generations. σ = .001 (red +), .02 (green x), .03(blue *), .05(pink ), and .08(sky blue *) and .1 (black ). (b) Temporal evolution of the histogram for the ratio V g (i)/V ip (i) for σ = 0.1. The histogram was computed over all networks and for every 20 generations. We used the logarithmic scales for both axes, so that the values for abscissa shows l o g 10 V g ( i ) V i p ( i ) : The histograms over generations 1-20, 21-40, 41-60,101-120, 181-200, and 281-300 were plotted with a different color for each network. With evolution, a peak was formed at V g (i)/V ip (i) = ρ ~ 10-2.3. See Method for the choice of the parameters.

The proportionality between V ip (i) and V g (i) was not a property of every gene expression dynamics but was evident only after the system achieved robustness through evolution. The fraction of genes showing a proportional relationship increased during the course of evolution. Indeed, for σ > σ c , the peak at r i ~ ρ increased over generations until it approached the distribution shown in Figure 3a (see Figure 3b). Summing up, the evolution of robustness was characterized by the formation of the peak at ρ < 1, in the distribution of r i = V ip (i)/V g (i).

The proportionality between the 2 variances implies the existence of a correlation between the noise- and mutation-induced changes in the gene expression statuses (see also Additional file 1, Figure S2 for the correlation in variances). Such a correlation was observed by computing the frequency of errors, i.e., changes in the on/off status of gene expression due to noise (without a change in the network) and as a result of mutation to the network (without adding the noise). The frequency of these 2 errors was highly correlated over genes for the GRN evolved at σ > σ c (see Figure 4). In other words, genes that were switched on or off more frequently by noise were also switched more frequently by mutation. This was in strong contrast with the GRN evolved at σ < σ c where no such correlation was observed (see Additional file 1, Figure S3). To sum up, the changeability of each gene expression level by noise and mutation was correlated, for a robust evolved system.

Figure 4
figure 4

Correlation between errors by noise and mutation over all genes. The frequency of errors by mutation versus noise plotted over all genes for 3 networks evolved at low noise σ = .1 > σ c . The abscissa showed the number of events for which the expression of the gene is switched on or off by noise. Computation was carried out as follows: Taking an evolved network (with highest fitness), we simulated the gene regulation network(GRN) dynamics without noise for a sufficiently long duration and obtained the original value of x(i) for each gene i. (i) The GRN dynamics were then simulated in the same manner under the noise level σ = .1 for over 105 different runs, and the number of switch events (i.e., the runs in which x(i) changes its sign) was computed. (ii) We generated 105 networks by changing 50 paths randomly chosen from the original network, and simulated these GRNs without noise and computed the number of switch events for each gene i. Step (i) gives the number of error by noise, and step (ii) gives the error by mutation, for each gene i. The relations between these 2 errors were plotted for 3 original networks with a different color for each network.


We confirmed that the proportionality between phenotypic variances of genetic and epigenetic origins held true for a system with evolved robustness, by simulating our model and its extended versions over several conditions. For all the conditions below, we confirmed (a) transition to robust evolution with the increase in noise level (b) proportionality between V ip (i) and V g (i) throughout the course of evolution and over many genes i.

(i) Against the change in k (target set) and M (number of genes): With an increase in the fraction of target genes, evolution to the fittest state became increasingly dicult, and the noise level for robust evolution σ c was slightly increased, but the proportionality was valid for σ > σ c (see Additional file 1, Figure S4a).

(ii) Considering that the density in the connection paths in the actual GRN is rather low, the validity of the results was verified by decreasing the path rates in the model. The conclusion remained unchanged as long as the network was percolated. The fraction of genes deviating from the proportionality slightly increased as the path rate in the network decreased (e.g., to .05 per gene) (see Additional file 1, Figure S4b).

(iii) Even if the noise level depended on each gene i, the conclusion was valid (see Additional file 1, Figure S5). After evolution, the variance in the expression level of each gene was not correlated with the noise level of each gene, implying that the fluctuation of gene expression was mainly controlled by (evolved) gene-to-gene regulation J ij .

(iv) To consider the influence of extrinsic rather than intrinsic noise, we also introduced the same level of noise to all gene expressions. Again, the robustness and the proportionality between the variances persisted as long as the level of the intrinsic noise was larger than σ c . Although the "extrinsic" (common) noise also contributed to the evolution of robustness, it played only a minor role in this respect. (See Additional file 1, Figure S6).

(v) By the variation in the environmental condition, the fitness condition for the target genes varied accordingly. By switching the condition for the target genes from 'on' to 'off' after some generations, we verified whether the evolution of GRN copes with this environmental variation. When the condition was switched, both the variances of epigenetic and genetic origins, as well as r i = V g (i)/V ip (i), increased to adapt novel environmental condition. Once the adaptation was achieved, the variances as well as r i decreased, to regain robustness (Figure 5).

Figure 5
figure 5

Change in the variances after the switching of the fitness condition. After evolution under the fitness condition to favor x i >0 ("on") for all the target genes i = 1, 2, ... k(= 8) as studied already, the fitness condition was switched at a certain generation to favor x i >0 for i = 1, 2, .., k/2 and x i <0 for i = k/2 + 1, .., k (i.e., the fittest gene expression pattern was given by + + + + − − −−, instead of + + + + + + ++). The switching was applied after a sufficiently large number of generations when the fittest networks are evolved (i.e., with x i >0 for the target genes). The switch initially caused a decrease in the fitness, but after a few dozens of generations, almost all networks evolved to adapt to the new fitness condition if σ > σ c . The values of the parameter and the procedure for computing the variances were identical with those in the previous cases. (i) The plot of the variances of the fitness, V g versus V ip per generation after the switching of the fitness condition. The color represented the generation number from the switching. There was a correlation between the increase in both the variances after the switch, and then, there was a proportional decrease as adaptation to the new fitness condition progressed. (b)(c) The plot of (V g (i) and V ip (i)) over all genes i at the generation 10, 20, and 30 for (b), and 60, 80, and 120 for (c). After the switch V g (i) and V ip (i) increased up to 30-60 generations, while the ratio V g (i)/V ip (i)approached unity for many genes. For generations >60, the variances decreased while the proportionality between V g (i) and V ip (i) was regained.

When the target phenotype was changed periodically, we observed that the increase and decrease of the variances were consistently repeated, when the noise level was near the transition value (σ c ), where rapid adaptation to new environment and robustness in phenotype were compatible.

(vi) We also extended our GRN model to account for diploids with sexual recombination. Here, each individual had a pair of matrices J i j 1 and J i j 2 , and the gene expression dynamics were given as a result of summation of the two matrices instead of the equation in the Method. By considering recombination of two matrices from a parent, we evolved GRN to achieve a higher fitness. The proportionality of the 2 variances was again confirmed, while another noteworthy finding was that in the case of heterozygotes, the robustness was further enhanced (suppresses the variances of expression).

Further, the proportionality between the variances was not confined to the present model. Indeed, in a model of catalytic-reaction-network, such a proportionality between the variances of fitness [16] and over chemicals evolved, whereas robustness transition with the increase of noise was confirmed in an abstract spin model for protein folding [22]. We expect that the relationship holds true as long as the fitness is determined through complex developmental dynamics with noise and the high-fitness states are not easily achieved so that error catastrophe appears with the increase in the mutation rate.

A phenomenological distribution theory

Considering the generality of the proportionality between the phenotypic variances of genetic and epigenetic origins, we provide a distribution theory for it without going into detailed setups of the model. We adopt the evolutionary stability argument first introduced for the proportionality between V ip and V g (of the fitness) through the course of evolution [16, 23].

We consider a multivariate distribution function with regard to gene expression level x i and the genotype a. Considering that multiple genes are involved, we assume that the genotype is represented by a scalar parameter a (e.g., by a Hamming distance from the fittest genetic sequence). Now, we assume evolutionary stability, in which the distribution maintains a single peak through evolution. Then, by a suitable transformation of variables, this peak position is taken to be 0; further, the form is approximated by Gaussian distribution around the origin to give the following equation

P ( x i , a ) = N 0 e x p ( x i 2 2 α i + C i x i a a 2 2 μ ) .

with N0, a normalization constant so that ∫P (x : a)dx = 1. Here, α i V ip (i) is the variance of the gene expression level, while μ is the mutation rate that determines the variance in the genotypes. Only a linear change of x i with regard to a is considered by neglecting higher order terms. Eq. (1) is rewritten as

P ( x i , a ) = e x p ( ( x i C i a α i ) 2 2 α i 1 2 ( 1 μ C i 2 α i ) a 2 )

Now, recall the stability of the distribution P (x i , a), i.e., whether it has a peak in the space with x i and a. This condition is given by 1 / μ C i 2 α i > 0 , i.e., μ < μ m a x i 1 C i 2 α i . For the mutation rate larger than μ m a x i , the distribution is flattened. In this case, the peaked distribution concentrated at a certain gene expression level is no longer sustained. This can be interpreted as a kind of error catastrophe originally introduced by Eigen [24], i.e., collapse of the sustenance of the localized distribution of functional genes. The critical mutation rate for the error catastrophe is given by μ = μ m a x j , which can take independent values by any gene i. However, GRNs that have achieved robustness to noise and mutation through evolution may have some constraints among expressions of different genes.

To have higher robustness, the error threshold for the fitness should be postponed to a higher mutation rate. When the system has achieved robustness to noise and mutation through evolution, the fitness level changes only to a small degree (i.e., remains almost neutral) against a considerable amount of change in the GRNs due to mutation [21]. Until the occurrence of this number of mutations, the genes rarely undergo any change in their expression statuses (on or off). This introduces a constraint on the change in gene expression against mutations.

If each of non-target genes were switched on or off independent of each other, the error in the expression of target genes that could be influenced by each switch would occur frequently. Indeed, the evolved robust GRN has some constraint on the errors by noise and mutation, as plotted in Figure 4. To achieve higher robustness, there needs to be some correlation between the changes in the expression statuses of genes. By suitable mutual interaction (J ij ≠ 0) among genes, the error frequency in the target gene expression can be reduced. This reduction works up to the mutation rate μ m a x i , while for μ ~ μ m a x i , errors in the expression status of one gene can be propagated to the expression of many other genes, which in turn will induce changes in the expression statuses of other genes. Hence, for a robust network having higher error threshold mutation rate, many genes may be switched on or off simultaneously within a GRN, once an error catastrophe in one gene expression occurs. Accordingly, many genes i are expected to share a common critical mutation rate. Hence, μ m a x i is roughly equal for many genes, when robustness is evolved, i.e., at σ > σ c . In fact, we computed the expression pattern of GRNs by adding a larger number of mutations to the evolved GRN, to obtain the variances V ip (i) and V g (i). When the original network was evolved under high noise conditions, the variances touched with the line V ip (i) = V g (i) simultaneously over many genes, as the mutation rate was increased (see Figure 6). The flattening of the distribution occurred at similar mutation rates over many genes.

Figure 6
figure 6

Plot of ( V ip ( i ), V g ( i )) over all genes, after 200 generations, computed in the same manner as that described in Figure 2. Instead of a single mutation per network adopted thus far, we used much a larger number of mutations. Each color represented the variances as a result of 20, 75, 100, and 140 mutations per network, which gave additions or deletions of paths chosen randomly in the gene regulation network (GRN). After mutation number ~ 100, many of genes touched the line V g (i) = V ip (i), where error catastrophe was percolated over many genes.

Following this argument, we may expect that

μ m a x i = ( C i 2 α i ) 1 = i n d e p e n d e n t o f m a n y g e n e s

when robustness is evolved (i.e., at σ > σ c ). Note that x i ¯ , the mean of x i for a given genotype a, is given by C i α i a according to eq. (2). The variance of x i ¯ due to this "genetic change" is given by the distribution of a. Thus, we get

V g ( i ) = < ( δ x i ) 2 > = C i 2 α i 2 < ( δ a ) 2 > .

Since V ip (i) = α i , we get

r i = V g ( i ) / V i p ( i ) = C i 2 α i < ( δ a ) 2 > ,

but this value is independent of gene i, according to eq.(3). Hence, the proportionality over genes is explained by a common error catastrophe threshold value over different gene expression levels.

As stated earlier, this proportionality over genes is not a general property of (gene expression) dynamical systems, but emerges only when evolution occurs under a sufficient level of noise (see Figure 3b). The use of the distribution function and assumption of stability implies that we are concerned with the stationary distribution that is attained after the progress in evolution. As the noise level reduces resulting in a decrease in robustness, the fraction of genes sharing a common value of V g (i)/V ip (i) = ρ decreases. This decrease is interpreted as a decrease in the number of genes that have a common error threshold value.


The relationship between the variances over genes, rather than the evolutionary course, will be easier to confirm experimentally because variances in this case need not be traced across many generations. By measuring directly the isogenic phenotypic variance and mutational variance over many genes (proteins), the correlation between the 2 can be examined. Although direct experimental support is not yet available, recent studies conducted by Laundry et al. [25] on such variances in yeast suggest the existence of a correlation between the 2 types of variances (see also [26]). In fact, they measured "expression noise" for each gene as the variance from its expression in isogenic organisms, and "mutational variance" as the variance of the change in the expression levels of the genes after the occurrence of mutations. The former corresponds to V ip (i), while the latter correlates with V g (i) since both measure the variations in phenotype (gene expression) induced by genetic changes. Although experimental data are scattered, a positive correlation is noted between the expression noise and mutational variance, as is consistent with the inference of the proportionality between V ip (i) and V g (i). Note that this proportionality holds true only for a set of genes whose expression levels are mutually related and directly or indirectly related with the fitness, whereas the experimental data cover all genes. In fact, by choosing a set of genes that have a stronger mutual relationship, the correlation between the variances is increased [27].

According to the theoretical argument we presented here, the phenotype variable is not restricted to gene expression level, but can represent any trait. Stearns et al. carried out selection experiments of Drosophila melanogaster on several fitness conditions such as age at eclosion, weight, and lifespan. They measured the variance of these phenotype traits within lines (i.e., corresponding to V ip (i)) and among lines (corresponding to V g (i)). Interestingly, the 2 types of variances (even after being normalized by the mean) showed remarkable proportionality over different phenotypic traits (see Figure 2 of [28]), in agreement with our theory. Note that they used the population selected on the basis of certain phenotypic traits, as in our model and theory. This must explain why the observed proportionality between the 2 variances is clearer than that in [25].


The characterization of robustness, evolvability, and plasticity [1, 2931] is an important issue in the field of evolutionary and developmental biology; however, studies on this issue are often qualitative. In the present study, we have demonstrated that the phenotypic fluctuations provide quantitative measures for these. Consider a population of organisms evolved under a single fitness condition, where the phenotypes that directly or indirectly influence the fitness are given as a result of (gene expression) dynamics under noise (determined by transcriptional networks). Through selection and mutation, the rules for the dynamics (i.e., the transcriptional networks) evolve leading to the achievement of a higher fitness level. Previously, we defined V ip as the variance of the fitness within an isogenic population, and V g as the variance of the average fitness within a heterogenic population, and obtained V ip V g evolution speed through the course of evolution [16, 17, 23].

In the present paper we defined the variances at each gene expression level i due to noise as V ip (i) and that due to mutation as V g (i). Then, the conclusion of the present paper is summarized as follows:

(1) For a population of organisms at a given generation evolved after some generations, V ip (i) V g (i) for most genes i. In other words, r i = V g (i)/V ip (i) took a common value (ρ < 1) over many genes, and the number of such genes increased as the robustness of fitness increased. (Total phenotype variance is given by V ip (i) + V g (i) if the 2 variances are added independently. In this case the heritability [18, 19], defined as the ratio of V g (i) to the total phenotypic variance is given by r i /(1 + r i ). Hence, the heritability takes a common value for mutually correlated traits, for evolved population under a fixed, single fitness condition). The previous relationship through the evolutionary course, implies that organisms with larger phenotypic variances have higher rates of evolution. The relationship (1) we found here implies that genes (or phenotypic traits) that have larger fluctuations have higher evolution speed.

(2) As the fraction of genes sharing a common value r i = ρ < 1 increased, there was a decrease in the degree of freedom to change the gene expressions independently by mutation. This increase in correlated change lead to an increase in robustness of the fitness to mutation and noise. On the other hand, the expression of genes with larger r i ~ 1 was easily switched by noise or mutation, and provided plasticity to environmental as well as mutational change (see also Figure 5).

The generality of our results was confirmed by several extensions of the model including environmental fluctuations, gene-dependent noise amplitudes, diploid with recombination, and so forth, as well as a catalytic reaction network model. Here, we should note that a correspondence between the change induced by noise in development and that induced by mutation was a source of correlation between V ip and V g in our study. Indeed, Waddington coined the term "genetic assimilation" as a process in which environment-induced phenotypic changes are subsequently embedded into genes [9]. The proportionality among phenotypic plasticity, V ip (i), and V g (i) is regarded as a quantitative expression of this genetic assimilation.


A simplified gene expression dynamics with a sigmoid input-output behavior [32, 33] is adopted here, although several simulations in the form of biological networks will give essentially the same result. In this model, the dynamics of a given gene expression level, x i , is described by the following:

d x i / d t = γ { tanh [ β j > k M J i j x j ] x i } + σ η i ( t ) ,

where J ij = 1, 1, 0, and η i (t) is a Gaussian white noise given by < η i (t)η j (t') > = δ i, j δ(t t'). M is the total number of genes, and k is the number of target genes that determine fitness. The summation only for j > k is introduced to eliminate possible influences from the target genes, which might also fix other gene expressions. Without this restriction and just by the summation over all genes, however, conclusions of the present numerical results are invariant. Of course, the matrix J ij is generally asymmetric. The amplitude of noise strength is given by σ that determines stochasticity in gene expression. The initial condition is given by (-1,-1,...,-1); i.e., all genes are off. The fitness F is determined by how many of the "target" genes are on after a sufficient time, i.e., the number of i such that x i > 0 for i = 1, 2, ..k < M. Because the model includes a noise component, the fitness can fluctuate at each run, which leads to a distribution in the fitness F and x i , even among individuals sharing the same gene regulation network. For each network, we compute the average fitness F ¯ over L runs.

At each generation, there are N individuals with slightly different J ij . Among the networks, we select those with higher fitness values. From the selected networks, J ij is "mutated," i.e., J ij for a certain pair i, j selected randomly with a certain fraction is changed among ± 1,0. For example, each of the N s (< N ) networks with higher values F ¯ produce N/N s mutants. (We also used the selection procedure such that the offspring number is proportional to the fitness, with normalization of the total population to N, but the same results presented here were obtained). We repeat this selection-mutation process over generations.

The variance of fitness or gene expression Sign(x i ) of identical networks over L runs gives V ip or V ip (i), while the variances of their mean over different networks J ij give V g or V g (i). we chose N = L = 200, and N s = N/4, while the conclusion to be shown below does not change as long as these values are sufficiently large. We use β = 7, γ = .1, M = 64 and k = 8, and initially chose J ij randomly with equal probability for ± 1,0, unless otherwise mentioned.


  1. de Visser JA, et al: Evolution and detection of genetic robustness. Evolution. 2003, 57: 1959-1972. 10.1554/02-750R.

    PubMed  Google Scholar 

  2. Wagner A: Robustness against mutations in genetic networks of yeast. Nature Genetics. 2002, 24: 355-361. 10.1038/74174.

    Article  Google Scholar 

  3. Ciliberti S, Martin OC, Wagner A: Robustness can evolve gradually in complex regulatory gene networks with varying topology. PLOS Comp Biology. 2007, 3: e15-10.1371/journal.pcbi.0030015.

    Article  Google Scholar 

  4. Wagner GP, Booth G, Bagheri-Chaichian H: A population Genetic Theory of Canalization. Evolution. 1997, 51: 329-347. 10.2307/2411105.

    Article  Google Scholar 

  5. Siegal ML, Bergman A: Waddington's canalization revisited: Developmental stability and evolution. Proc Nat Acad Sci USA. 2002, 99: 10528-10532. 10.1073/pnas.102303999.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Barkai N, Leibler S: Robustness in simple biochemical networks. Nature. 1997, 387: 913-917. 10.1038/43199.

    Article  CAS  PubMed  Google Scholar 

  7. Alon U, Surette MG, Barkai N, Leibler S: Robustness in bacterial chemotaxis. Nature. 1999, 397: 168-171. 10.1038/16483.

    Article  CAS  PubMed  Google Scholar 

  8. Schmalhausen II: (1949, reprinted 1986) Factors of Evolution: The Theory of Stabilizing Selection (University of Chicago Press, Chicago).

  9. Waddington CH: The Strategy of the Genes. 1957, (Allen & Unwin, London)

    Google Scholar 

  10. Elowitz MB, Levine AJ, Siggia ED, Swain PS: Stochastic gene expression in a single cell. Science. 2002, 297: 1183-1187. 10.1126/science.1070919.

    Article  CAS  PubMed  Google Scholar 

  11. Furusawa C, Suzuki T, Kashiwagi A, Yomo T, Kaneko K: Ubiquity of log-normal distributions in intra-cellular reaction dynamics. Biophysics. 2005, 1: 25-31. 10.2142/biophysics.1.25.

    Article  CAS  Google Scholar 

  12. Bar-Even A, et al: Noise in protein expression scales with natural protein abundance. Nature Genetics. 2006, 38: 636-643. 10.1038/ng1807.

    Article  CAS  PubMed  Google Scholar 

  13. Kaern M, Elston TC, Blake WJ, Collins JJ: Stochasticity in gene expression: From theories to phenotypes. Nat Rev Genet. 2005, 6: 451-464. 10.1038/nrg1615.

    Article  CAS  PubMed  Google Scholar 

  14. Kaneko K: Life: An Introduction to Complex Systems Biology. 2006, (Springer, Heidelberg and New York)

    Google Scholar 

  15. Sato K, Ito Y, Yomo T, Kaneko K: On the relation between fluctuation and response in biological systems. Proc Nat Acad Sci USA. 2003, 100: 14086-14090. 10.1073/pnas.2334996100.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Kaneko K, Furusawa C: An evolutionary relationship between genetic variation and phenotypic fluctuation. J Theo Biol. 2006, 240: 78-86. 10.1016/j.jtbi.2005.08.029.

    Article  Google Scholar 

  17. Kaneko K: Evolution of robustness to noise and mutation in gene expression dynamics. PLoS ONE. 2007, 2: e434-10.1371/journal.pone.0000434.

    Article  PubMed  PubMed Central  Google Scholar 

  18. Futuyma DJ: Evolutionary Biology. 1986, (Sinauer Associates Inc., Sunderland), Second

    Google Scholar 

  19. Hartl DL, Clark AG: Principles of Population Genetics. 2007, (Sinauer Assoc. Inc., Sunderland), 4

    Google Scholar 

  20. Fisher RA: The Genetical Theory of Natural Selection. 1930, (Oxford University Press)

    Book  Google Scholar 

  21. Kaneko K: Shaping Robust system through evolution. Chaos. 2008, 18: 026112-10.1063/1.2912458.

    Article  PubMed  Google Scholar 

  22. Sakata A, Hukushima K, Kaneko K: Funnel landscape and mutational robustness as a result of evolution under thermal noise. Phys Rev Lett. 2009, 102: 148101-10.1103/PhysRevLett.102.148101.

    Article  PubMed  Google Scholar 

  23. Kaneko K: Relationship among Phenotypic Plasticity, Genetic and Epigenetic Fluctuations, Robustness, and Evolovability;. J BioSci. 2009, 34: 529-542. 10.1007/s12038-009-0072-9.

    Article  PubMed  Google Scholar 

  24. Eigen M, Schuster P: The Hypercycle. 1979, (Springer, Heidelberg)

    Book  Google Scholar 

  25. Landry CR, Lemos B, Dickinson WJ, Hartl DL: Genetic Properties Influencing the Evolvability of Gene Expression. Science. 2007, 317: 118-10.1126/science.1140247.

    Article  CAS  PubMed  Google Scholar 

  26. Lehner B: Selection to minimize noise in living systems and its implications for the evolution of gene expression. Mol Syst Biol. 2008, 4: 170-10.1038/msb.2008.11.

    Article  PubMed  PubMed Central  Google Scholar 

  27. Lehner B: Genes Confer Similar Robustness to Environmental, Stochastic, and Genetic Perturbations in Yeast PLoS One. 2010, e9035-

    Google Scholar 

  28. Stearns SC, Kawecki TJ: The differential genetic and environmental canalization of fitness components in Drosophila melanogaster. J Evol Biol. 1995, 8: 539-557. 10.1046/j.1420-9101.1995.8050539.x.

    Article  Google Scholar 

  29. Kirschner MW, Gerhart JC: The Plausibility of Life. 2005, (Yale University Press)

    Google Scholar 

  30. Ancel LW, Fontana W: Plasticity, evolvability, and modularity in RNA. J Exp Zool. 2002, 288: 242-283. 10.1002/1097-010X(20001015)288:3<242::AID-JEZ5>3.0.CO;2-O.

    Article  Google Scholar 

  31. West-Eberhard MJ: Developmental Plasticity and Evolution. 2003, (Oxford University Press)

    Google Scholar 

  32. Mjolsness E, Sharp DH, Reisnitz J: A connectionist model of development. J Theor Biol. 1991, 152: 429-453. 10.1016/S0022-5193(05)80391-1.

    Article  CAS  PubMed  Google Scholar 

  33. Salazar-Ciudad I, Garcia-Fernandez J, Sole RV: Gene networks capable of pattern formation: from induction to reaction-diffusion. J Theor Biol. 2000, 205: 587-603. 10.1006/jtbi.2000.2092.

    Article  CAS  PubMed  Google Scholar 

Download references


I would like to thank B. Lehner, T. Yomo, C. Furusawa, M. Tachikawa, and S. Ishihara, for stimulating discussion. This research was supported by the people of Japan, via a Grant-in-Aid for scientific research(No.21120004) from MEXT Japan.

Author information

Authors and Affiliations


Corresponding author

Correspondence to Kunihiko Kaneko.

Electronic supplementary material


Additional file 1: Figure S1 Dependence of the variances V ip ( i ) and V g ( i ) on the noise strength. Figure S2 Correlation between gene expressions. Figure S3 Correlation between errors by noise and mutation over all genes. Figure S4 Relationship between V (i) and V ip (i) for evolved GRNs with a larger fraction of target genes, and a smaller fraction of nonzero genes. Figure S5 Relationship between V g (i) and V ip (i) for the gene expression dynamics whose noise level depends on each gene. Figure S6 Relationship between V g (i) and V ip (i) for a model with "extrinsic noise." (PDF 66 KB)

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

Kaneko, K. Proportionality between variances in gene expression induced by noise and mutation: consequence of evolutionary robustness. BMC Evol Biol 11, 27 (2011).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI:


  • Gene Expression Level
  • Gene Regulation Network
  • High Fitness
  • Noise Strength
  • Evolution Speed