UV intensity may affect autoimmune disorder.

Phenotype of biological systems needs to be robust against mutation in order to sustain themselves between generations. On the other hand, phenotype of an individual also needs to be robust against fluctuations of both internal and external origins that are encountered during growth and development. Is there a relationship between these two types of robustness, one during a single generation and the other during evolution? Could stochasticity in gene expression have any relevance to the evolution of these types of robustness? Robustness can be defined by the sharpness of the distribution of phenotype; the variance of phenotype distribution due to genetic variation gives a measure of 'genetic robustness', while that of isogenic individuals gives a measure of 'developmental robustness'. Through simulations of a simple stochastic gene expression network that undergoes mutation and selection, we show that in order for the network to acquire both types of robustness, the phenotypic variance induced by mutations must be smaller than that observed in an isogenic population. As the latter originates from noise in gene expression, this signifies that the genetic robustness evolves only when the noise strength in gene expression is larger than some threshold. In such a case, the two variances decrease throughout the evolutionary time course, indicating increase in robustness. The results reveal how noise that cells encounter during growth and development shapes networks' robustness to stochasticity in gene expression, which in turn shapes networks' robustness to mutation. The necessary condition for evolution of robustness, as well as the relationship between genetic and developmental robustness, is derived quantitatively through the variance of phenotypic fluctuations, which are directly measurable experimentally.


INTRODUCTION
Robustness is ability to function against changes in the parameter of a system [1,2,3,4,5]. In a biological system, the changes have two distinct origins, genetic and epigenetic. The former concerns with genetic robustness, i.e., rigidity of phenotype against mutation, which is necessary to maintain a high fitness state. The latter concerns with fluctuation in number of molecules and external environment.
Indeed, phenotype of isogenic individual organisms is not necessarily identical. Chemotaxis [6], enzyme activities, and protein abundance [7,8,9,10,11] differ even among those sharing the same genotype. Recent studies on stochastic gene expression elucidated the sources of fluctuations [7]. The question most often asked is how some biological functions are robust to phenotypic noise [11,12], while there may also be positive roles of fluctuations in cell differentiation, pattern formation, and adaptation [13,14,15,16].
Noise, in general, can be an obstacle in tuning a system to the fittest state and maintaining it there. Phenotype of an organism is often reproducible even under fluctuating environment or under molecular fluctuations [2]. Therefore, phenotype that is concerned with fitness is expected to keep some robustness against such stochasticity in gene expression, i.e., robustness in 'developmental' dynamics to noise. Phenotype having a higher fitness is maintained under noise. How is such ''developmental robustness'' achieved through evolution? In the evolutionary context, on the other hand, another type of robustness, robustness to mutation need to be considered. When genetic changes occur, gene expression dynamics are perturbed so that phenotype with a high fitness may no longer be maintained. The ''genetic robustness'' concerns with the stability of a high-fitness state against mutation.
Whether these two types of robustness emerge under natural selection have long been debated in the context of developmental dynamics and evolution theory [3,5,17,18,19,20], since the proposition of stabilization selection by Schmalhausen [21] and canalization by Waddington [22,23,24]. Are developmental robustness to noise and genetic robustness to mutation related? Is phenotypic noise relevant to attain robustness to mutation? In the present paper, we answer these questions quantitatively with the help of dynamical network model of gene expression.
Under the presence of noise in gene expression, phenotype as well as fitness, of isogenic organisms is distributed, usually following a bell-shaped probability function. When the phenotype is less robust to noise, this distribution is broader. Hence, the variance of this distribution, i.e., variance of isogenic phenotypic fluctuation denoted as V ip , gives an index for robustness to noise in developmental dynamics. On the other hand, robustness to mutation is measured from the fitness distribution over individuals with different genotypes. An index for it is given by variance of phenotypic fluctuation arising from diversity of genotypes in a population [25,26,27], denoted here as V g . This variance V g increases as the fraction of low-fitness mutants increases.
Here we show that evolution to increase both types of robustness is possible only when the inequality V ip $V g is satisfied.
Since the isogenic phenotypic fluctuation V ip increases with noise, this means that evolution of robustness is possible only when the amplitude of phenotypic noise is larger than some critical value as derived by V ip $V g , implying a positive role of noise to evolution. We demonstrate that both the two variances V ip and V g decrease in the course of evolution, while keeping the proportionality between the two. This proportionality is consistent with an observation in a bacterial evolution experiment [16,17,18].
We explain the origin of the critical noise strength, by noting that smooth dynamical behavior free from a rugged potential landscape evolves as a result of phenotypic noise. When the noise amplitude is smaller than the threshold, we observe that low-fitness mutants are accumulated, so that robustness to mutation is not achieved. Generality and relevance of our results to biological evolution are briefly discussed.

Theoretical Framework on Genetic-Phenotypic Relationship
In natural population, both the phenotype and genotype differ among individuals. Let us consider population distribution P(x, a) where x is a variable characterizing a phenotype and a is that for the corresponding genotype [18]. Here the phenotype x is responsible for the fitness of an individual, and the selection depending on x is considered as an evolutionary process. Since the phenotype differs even among isogenic individuals, the distribution P(x; a = a 0 ) for a fixed genotype a 0 has some variance. This isogenic phenotypic variance V ip , defined as the variance over clones, is written as V ip (a)~Ð (x{x(a)) 2 P(x,a)dx, where x(a) is the average phenotype of a clonal population sharing the genotype a, namely x(a)~Ð P(x,a)xdx. This variation of phenotype is a result of noise through the developmental process to shape the phenotype. If this variance is smaller, the phenotype is less influenced by noise, and thus V ip works as a measure of robustness of the phenotype against noise.
On the other hand, the standard evolutionary genetics [25,26,27] mainly studies the phenotypic variance due to genetic variation. It measures phenotypic variability due to diversity in genotypes in a population. This phenotypic variance by genetic variance, which is termed V g here, is then defined as the variance of the average x(a), over genetically heterogeneous individuals. It is given by V g~Ð (x(a){vxw) 2 P(a)da, where P(a) is the distribution of the genotype a and ,x _ . is the average of x(a) over genotypes a. While V ip is defined as variance over clones, i.e., individuals with the same genotype, V g comes from those with different genotypes. As V g is smaller, the phenotypic change by genetic variation is smaller. Hence V g gives a measure of robustness of the phenotype against mutation.
We have previously derived an inequality V ip .V g between the two variances, by assuming evolutionary stability of the population distribution P(x, a), that is preservation of single-peakedness through the course of evolution [18] (see Supporting Text S1). Indeed the single-peaked distribution collapses as V ip approaches V g , where the distribution is extended to very low values of x (fitness). In other words, error catastrophe occurs at V g <V ip ; (Here error catastrophe means accumulation of low-fitness mutants in the population after generations, and the term is used here by extending its original meaning by Eigen [28]). For each course of evolution under a fixed mutation rate, the proportionality between V g and V ip is derived, since the genetic variance increases roughly proportionally to the mutation rate [18].
Note, however, that the derivation of these relationships (V ip $V g , error catastrophe at V g <V ip , and proportionality between V g and V ip for a given course of evolution) is based on the existence of two-variable distribution function P(x = phenotype, a = gene), and the postulate that single-peaked distribution is maintained throughout evolution, which is not trivial. Hence the above relationships need to be examined by some models for evolution. In addition, why does the population distribution extend to low-fitness values when the phenotypic fluctuation V ip is smaller than V g ? Or, put it another way, why do systems with small phenotypic noise run into ''error catastrophe''? In fact, the emergence of error catastrophe as a result of decreasing isogenic phenotypic fluctuation below V g may look rather counter-intuitive, since in general one expects fluctuation to perturb a system from the fittest state. The necessity of fluctuation for evolution to increase robustness to noise and to mutation needs theoretical examination.

Model
To study the proposed relationships, we need to consider seriously how the phenotype is shaped through complex ''developmental process''. In the present paper, we use the term 'development', in a broad sense, including a process in uni-cellular organisms to reach cell division. It is a dynamical process to shape a phenotype at a 'matured' state (where fitness is defined) from a given initial state. In general, this dynamic process is complex so that the process may not reach the identical phenotype due to the noise through this developmental process. This leads to the isogenic variance of the phenotype V ip . On the other hand, the equation governing the developmental process is varied as a result of mutation. The phenotype variance over a population with distributed genotypes gives V g .
We consider a simple model to satisfy the requirement on 'development' above. It consists of a complex dynamic process to reach a target phenotype under a noise which may alter the final phenotypic state. We do not choose a biologically realistic model that describes a specific developmental process, but instead take a model as simple as possible, to satisfy a minimal requirement for our study. Here we take a simplified model, borrowed from a gene regulatory network, where expression of a gene activates or inhibits expression of other genes under noise. These interactions between genes are determined by the network. The expression profile changes in time, and eventually reaches a stationary pattern. This gene expression pattern determines fitness. Selection occurs after introduction of mutation at each generation in the gene network. Among the mutated networks, we select a network with a higher fitness value. Since there is a noise term in the gene expression dynamics, fitness fluctuates even among the individuals with an identical gene network, which leads to the isogenic fluctuation V ip . On the other hand, the expression pattern varies by mutation in the network, and gives rise to variation in the average fitness, resulting in V g .
This simplified gene expression follows a typical switch-like dynamics with a sigmoid input-output behavior [29,30,31,32,33] widely applied in models of signal transduction [34] and neural networks [35] (For a related evolution model with discrete states, see e.g., [24]). The dynamics of a gene expression level x i is described by where J ij = 21,1,0, and g(t) is Gaussian white noise given by ,g(t)g(t9). = d(t2t9). M is the total number of genes, and k is the number of output genes that are responsible for fitness to be determined. The value of s represents noise strength that determines stochasticity in gene expression (For simplicity we mainly consider the case that the noise amplitude is independent of x i , while inclusion of such x-dependence of noise amplitude does not alter the conclusion to be discussed). By following a sigmoid function tanh, x i has a tendency to approach either 1 or 21, which is regarded as ''on'' or ''off'' of gene expression. Even though x is defined over [2', '], it is attracted to the range [21,1] (or slightly above or below the range due to the noise term). We consider a developmental process leading to a matured phenotype from a fixed initial state, which is given by (21,21,…,21); i.e., all genes are off, unless noted otherwise. (This specific choice of initial condition is not important). Let us define a fitness function so that gene expression levels (x i ) for genes i = 1,2,…,k(,M) would reach an ''on'' state, i.e., x i .0. The fitness is maximum if all k genes are on after a transient time span T ini , and minimum if all are off. To be specific, we define the fitness function by where S(x) = 1 for x.0, and 0 otherwise, […] temp is time average between t = T ini and t = T f (The time average here is not important, because the gene expressions x i are fixed after some time, in most cases). Adoption of the value (S(x j )21) after initial time T ini leads to the same result (. The fitness function takes the maximum value F = 0 when the selected pattern of gene expression (x i ; i = 1,2,…,k) is always ''on'' and takes the minimum (F = 2k) when all k genes are always off. Note that fitness is calculated only after time T ini , which is chosen sufficiently large so that the temporal average can be computed after the gene expression dynamics has fallen on an attractor. This initial time can be considered as the time required for developmental dynamics.
As the model contains a noise term, fitness fluctuates at each run, which leads to the distribution in F, even for the same network. Hence we obtain the distribution p(F; g), for a given network ''g'', whose variance gives isogenic phenotypic fluctuation. At each generation, we compute the fitness F over L runs, to obtain the average fitness value F _ of a given network. Now we consider the evolutionary process of the network. Since the network is governed by J ij which determines the 'rule' of the dynamics, it is natural to treat J ij as a measure of genotype. Individuals with different genotype have a different set of J ij At each generation there are N individuals with different sets of J ij For each individual network, we compute the average fitness F _ . Then we select the top N s (,N) networks that have higher fitness values. (The value N/N s corresponds to the selective pressure. As it is larger, the evolution speed increases. However, specific choice of this value itself is not important to the result to be discussed).
At each generation, mutation changes the network, i.e., changes J ij at a given mutation rate m. We rewire the network at a given rate so that changes in J ij produce N new networks. (In most simulations, only a single path, i.e., a single pair of i, j is changed. The mutation rate can be lowered by changing a path only for some probability. Although it is important to the evolution speed and the error catastrophe point to be discussed, the conclusion to be discussed is not altered by specific choice of m.) Here we make N/N s mutants from each of the top N s networks, so that there will be N networks again for the next generation. From this population of networks we repeat the process of the developmental dynamics, fitness calculation, selection, and mutation (Instead of this simple genetic algorithm, we can also assume that the number of offspring increases with the fitness. This choice does not alter the conclusion to be presented).
Simulations start from a population of random networks with a given fraction of paths (for example, 50% of J ij are nonzero). At each generation, the N individuals have slightly different networks J ij , so that the values of F _ are different. We denote the fitness distribution over individuals with different genotype as P(F _ ). On the other hand, the fitness distribution for an identical network (''g'') is computed to obtain p(F; g).

Remark:
Developmental dynamics and selection process in our model are too simple. Still, this model is relevant to examine general statement on phenotypic fluctuations, as the model at least captures complex dynamics giving rise to a phenotype, stochasticity in dynamics, mutation, and selection according to a given phenotype. Indeed, real gene expression dynamics depend on environmental conditions, and the fitness is defined as expression patterns to adapt each environmental condition. We have also carried out some simulations by imposing such fitness but the results to be discussed (with regards to V g and V ip ) are invariant.

RESULTS
Let us first see how the evolutionary process changes as a function of the noise strength s. After generations, the peak position in P(F _ ) increases, so that the top of F _ in the population approaches the highest value 0. Indeed, in all cases, the top group quickly evolves to the highest fitness state F _ = 0 (see Fig. 1a; even for s = 0.2, the highest fittest value approaches 0 after a few hundred more generations.). The time necessary for the system to reach this state becomes shorter as the phenotypic noise decreases (see Fig. 2). On the other hand, the time evolution of the distribution function P(F _ ) depends drastically on the noise strength s. When s is small, the distribution is broad and the existing individual with the lowest F _ remains at the low fitness state, while for large s, even the individuals with the lowest fitness approach F _ = 0 (see Fig. 1b and Fig. 3). There is a threshold noise s c (<0.02), below which the distribution P(F _ ) is broadened, as is discernible in the data of the variance of the distribution, V g in Fig. 2. Here, the top individuals reach the highest fitness, leaving others at the very low fitness state. As a result, the average fitness over all individuals, vF w~Ð FP(F )dF is low. ,F _ . and the lowest fitness over individuals F _ min , after a sufficiently large number of generations, are plotted against s in Fig. 2. The abrupt decrease in fitness suggests threshold noise s c , below which low-fitness mutants always remain in the distribution. For s.s c , the distribution P(F _ ) takes a sharp peak at F _ ,0, where the variance is rather small. Distribution below and above s c are displayed in Fig. 3. (This type of transition is also observed by increasing the mutation rate, while fixing the noise strength at s.s c ).
Let us study the relationship between V g and V ip Here V ip is defined as variance from the distribution p(F; genotype), i.e., over individuals with the same genotype. As the distribution p depends on each individual with different genotype, the variance changes accordingly. Naturally, the top individual has a smaller variance, and the individual with lower fitness has a larger variance. As a measure of V ip , we used either the average of the variance over all individuals or the variance of phenotype from a gene network that is located closest to the peak in the distribution P(F _ ). Both estimates of V ip do not differ much, and the following conclusion is drawn in both cases. V g , on the other hand, is simply the variance of the distribution P(F _ ), i.e., over individuals having different genotypes present.
The relationship between V g and V ip thus evaluated is plotted in Fig. 4. We find that both the variances decrease through the evolutionary time course when s.s c , where we note: (ii) V g !V ip during the evolutionary time course under a fixed noise strength s(.s c ) and a fixed mutation rate. As s is lowered toward s c , V g increases so that it approaches V ip .
(iii) V g <V ip at s<s c , where error catastrophe occurs. In other words, the fittest networks maintaining a sharp distribution around the peak dominate only when V ip .V g is satisfied. As s is decreased to s c , V ip approaches V g , error catastrophe occurs and a considerable fraction of low-fitness mutants accumulates. Hence, the relationships proposed theoretically assuming evolutionary stability of a two-variable distribution function P(x = phenotype, a = genotype) is confirmed. Here, without introducing phenomenological assumptions, the three relationships are observed in a general stochastic gene-network model. Why didn't the system maintain the highest fitness state under small phenotypic noise s,s c ? To study the difference in dynamics evolved for different values of s, we choose the top individual (network) that has evolved at s = s 0 , and place it under a different noise strength s = s9. In Fig. 5, we have plotted the fraction of runs giving rise to F = 0 under such circumstance. As shown, the successful fraction decreases when s9 goes beyond s 0 . In other words, the network evolved under a given noise strength successfully reaches the target gene expression up to that critical noise level, while it begins to fail doing so at a higher noise strength. Accordingly, the distribution p(F; gene) extends to lower values in fitness, when a network evolved under small phenotypic    noise is developed under a higher noise level. On the other hand, the network evolved under high level noise maintains a high fitness value, even when the noise is lowered.
Next we study the basin structure of attractors in the present system. Note that an orbit of the network with the highest fitness, starting from the prescribed initial condition, is within the basin of attraction for an attractor corresponding to the target state (x i .0 for i = 1,…,k). Hence the basin of attraction for this target attractor is expected to be larger for the dynamics evolved under higher level noise. We have simulated the dynamics (1) for the evolved fittest network under zero noise, starting from a variety of initial conditions over the entire phase space, and measured the distribution of F at each attractor. The distribution is shown in Fig. 6 (Due to the symmetry against x j = 1«x j = 21 in the model, the distribution is symmetric around F = 2k/2 when all initial conditions are taken. In fact, by starting from x i = 1 for all i, the orbit reaches an attractor x j ,0 for j = 1,…,k, resulting in F = 2k). For the network evolved under s.s c , the distribution has a sharp peak at F = 0 (and F = 2k due to the symmetry), with more than 40% attraction to each. On the other hand, for the networks evolved under s,s c , the peak height at F,0 is very small, i.e., the basin for the attractor with F = 0 is tiny. There exist many small peaks corresponding to attractors with 2k,F,0, having similar sizes of basin volumes. In fact, the basin volumes for attractors with 2k,F,0 grow as s is decreased, and are dominant for s,s c.

Dynamic Origin of Robust Evolution
The difference in the basin structure suggested by Fig. 6 is schematically displayed in Fig. 7. For the network evolved under s.s c there is a large, smooth attraction to the target state, while for the dynamics evolved under s,s c , the phase space is split into small basins. Let us consider the distance between the basin boundaries from a ''target orbit'' starting from 21,…,21 and reaching x i .0 (for 1#i#k) which is defined by D here. The distance D remains small for the dynamics evolved under a low noise strength s,s c , and increases for those evolved under a higher noise. It is interesting to note that evolution influences the basin structure globally over the phase space, although the fitness condition is applied locally to an orbit starting from a specific initial condition.
The results in Fig. 5 and Fig. 6 imply that the gene regulation networks that operate and evolve under noisy environment exhibit qualitatively different dynamics compared to those subjected to low level noise. In our model, the fitness of an individual changes when its gene expression x j for j = 1,…,k changes its sign. Recall that the fixed point solution x i = tanh(S j J ij x j ) changes its sign when S j J ij x j in the sigmoid function changes its sign. This change may occur during the developmental dynamics by noise, and we call  such points in the phase space 'turning points'. When an orbit of eq.(1) passes over turning points, x j takes a negative value for some j (for 1#i#k) at the attractor (see Fig. 8 for schematic representation). Since there are many variables for gene expression and the values of J ij are distributed over 21 and 1, the term tanh(S j J ij x j ) changes its sign at several points in the phase space {x j } generally. Hence there can be many turning points in the phase space. The fittest network with F _ <0 chooses such orbits having no turning points within the noise range from the original orbit. An orbit from the fittest individual evolved under low-level noise  Figure 6. Distribution of the fitness value when the initial condition for x j is not fixed at 21, but is distributed over [21,1]. We choose the evolved network as in Fig. 5, and for each network we take 10000 initial conditions, and simulated the dynamics (1) without noise to measure the fitness value F after the system reached an attractor (as the temporal average 400,t,500). The histogram is plotted with a bin size 0.1. doi:10.1371/journal.pone.0000434.g006 encounters many turning points when subjected to noisy environment. The average distance between the turning points and an orbit that has reached the target gene expression pattern is estimated by the distance D defined above. Recall that the distance D is small for the dynamics evolved under a low noise strength. Such dynamics, if perturbed by a higher level of noise, are easily caught in the turning points, which explains the behavior shown in Fig. 5.
Let us now discuss the relationship between V g and V ip . Noise disturbs an orbit so that it may go across the basin boundary (turning points) with some probability. We denote the standard deviation of the location of the orbit due to noise as d p , which is proportional to the noise strength s. Since the distance between the orbit and the basin boundary is deviated by d p , and the fitness value drops when the orbit crosses the basin boundary, the variance V ip is estimated to be proportional to (d p /D) 2 .
Next, we discuss how the mutation in the network influences the dynamics. When the network is altered, i.e., a path is added or removed as a result of mutation in J ij , there exists a variation of the order of O(1= ffiffiffiffi ffi N p ) in the threshold function term in eq.(1). This leads to a deviation of the location of turning points (or basin boundary). We denote this deviation as d g , which increases with the mutation rate. The variance V g is estimated to be proportional to (d g /D) 2 , with the same proportion coefficient as that between V ip and (d g /D) 2 .
Under the presence of noise, there is a selective pressure to avoid the turning points (basin boundaries) that exist within the distance d p from the ''target'' orbit. This leads to an increase in D. However, if d p is larger than d p , the memory of this distance between the target and the boundaries will not be propagated to the next generation, due to large perturbation to the original network by the mutation. Hence increase in D (i.e., increase in robustness) is expected only if d p .d g . Since d p and d g increase with the noise strength s and the mutation rate m respectively, there exists a critical noise strength s c beyond which this inequality is satisfied. From the relationship between d p,g and V ip,g , the condition for robust evolution is rewritten as V ip .V g .
When the condition V ip .V g (i.e., s.s c ) is satisfied, the system increases D during evolution. We have computed the temporal evolution of basin distribution. With generations, the distribution evolves from the pattern at a low level noise in Fig. 7, to that at large s characterized by enhanced peak at F = 0. Accordingly D increases with generations. Recall that V ip ,(d p /D) 2 , and V g ,(d g / D) 2 , both variances decrease with generations, while V ip /V g is kept constant.

DISCUSSION
We have demonstrated the inequality and proportionality between V g and V ip , through numerical evolution experiment of a gene network. As phenotypic noise is decreased and the inequality V ip .V g is broken, low-fitness mutants are no longer eliminated. This is because the mutants fail to reach the target gene expression pattern, by crossing the boundary of the basin of attraction to the target. When the amplitude of the noise is larger, on the other hand, the networks of the dynamics with a large basin volume for the target attractor are selected and thus mutants with lower fitness are removed successively through the selection process. Hence noise increases developmental robustness through evolution, together with genetic robustness.
Although we used a specific example to demonstrate the relationship between V ip and V g and the error catastrophe, we expect this relationship to be generally applicable to systems satisfying the following conditions: (i) Fitness is determined through developmental dynamics.
(ii) Developmental dynamics is sufficiently complex so that its orbit, when deviated by noise, may fail to reach the state with the highest fitness.
(iii) There is effective equivalence between mutation and noise in the developmental dynamics with regards to phenotype change.
Note that the present system as well as the previous cell model [18] satisfies these conditions. The condition (i) is straightforward in our model, and the condition (ii) is satisfied because of the complex dynamics having many turning points in the phase space. Noise in developmental dynamics sometimes perturbs an orbit to cross the basin boundary so as to escape from attraction to the target gene expression pattern, while a mutation in the network may also induce such failure, by shifting basin boundaries. Hence the condition (iii) is satisfied.
When developmental process fails due to phenotypic noise, the fitness function takes a low value. Evolution under noise acts to prevent such failure within the range of noise. On the other hand, due to the condition (iii), mutation may also lead to such lethality. When the effect of mutation goes beyond the range given by the phenotypic noise, mutants with very low fitness values begin to accumulate. Hence there appears a threshold level of phenotypic noise below which low-fitness (or deleterious) mutants accumulate (or threshold mutation rate beyond which such mutants accumulate). In this sense, we expect that for robust evolution, the inequality V g ,V ip must be satisfied in order for the low-fitness mutants to be eliminated. Violation of the inequality leads to accumulation of low-fitness (or deleterious) mutants, a phenomenon known as error catastrophe [28]. Only under the presence of noise in the developmental process, systems acquire robustness through evolution. In other words, developmental robustness to stochasticity in gene expression implies genetic robustness to mutation. Quantitative analyses of stochasticity in protein abundance during the laboratory evolution of bacteria are possible [17,36]. By carefully measuring the variation V g of given phenotype in mutants, and comparing it with that of isogenic bacteria, V ip , one can examine the validity of our conclusion between V g and V ip .
It is worthwhile to mention that in a class of theoretical models, fitness landscape is given as an explicit continuous function of a gene sequence (e.g., energy function in a spin glass [37]), where a minute change in sequence does not lead to a drastic change in fitness. On the other hand, in a system satisfying (i) and (ii), a small change in genotype (e.g., a single change in the network path) may result in a large drop in fitness, since fitness is determined after the developmental dynamics. Indeed, there may appear mutants with very low fitness values from an individual with a high fitness value, only by a single change of a path in the network. Such deleterious mutations are also observed in nature [27].
It is interesting to note that a larger basin of attraction to a target attractor (with the highest fitness value) is formed through a mutation and selection process. As a result, dynamics over the entire phase space are simplified to those having only a few attractors, even though the fitness function is given locally without scanning over the entire phase space. When the time-course is represented as a motion along a potential surface, our results suggest that the potential landscape becomes smoother and simpler through evolution, and loses ruggedness after generations. Indeed, existence of such global attraction in an actual gene network has recently been reported in yeast cell-cycle [38].
Such smooth landscape was also studied in protein folding [39,40]. Saito et al. [41] observed an evolutionary process from a rugged to the so-called funnel-like landscape in an interacting spin system abstracting protein folding dynamics. Under a general framework of statistical mechanics [42], a relationship between the degree of variance in coupling coefficients J ij between spins (corresponding to V g ) and the temperature (i.e., phenotypic noise for spin x i , corresponding to V ip ) is formulated. Such relationship may be relevant to understand the relationships between V g and V ip in our study.
According to established Fisher's theorem on natural selection, evolution speed of phenotype is proportional to the phenotypic variance by genetic variation, V g [25,26,27]. The demonstrated proportionality between V ip and V g then suggests that the evolution speed is proportional to the isogenic phenotypic fluctuation, as is also supported by an experiment on bacterial evolution in a laboratory [17] and confirmed by simulations of a reaction network model of a growing cell [18].
Isogenic phenotypic fluctuation is related to phenotypic plasticity, which is a degree of phenotype change in a different environment. Positive roles of phenotypic plasticity in evolution have been discussed [20,43,44,45]. Since susceptibility to the environmental change and the phenotypic fluctuation are positively correlated according to the fluctuation-response relationship [16,46.47], our present results on the relationship between phenotypic fluctuations and evolution imply, inevitably, a relationship between phenotypic plasticity and evolution akin to genetic assimilation proposed by Waddington [22].

SUPPORTING INFORMATION
Text S1 Derivation of General Relationship on Fluctuations. Mathematical derivation on general relationships on phenotypic variances is presented. Found at: doi:10.1371/journal.pone.0000434.s001 (0.05 MB PDF)