Evolution of Robustness and Plasticity under Environmental Fluctuation: Formulation in terms of Phenotypic Variances
aa r X i v : . [ q - b i o . P E ] M a y Noname manuscript No. (will be inserted by the editor)
Kunihiko Kaneko
Evolution of Robustness and Plasticityunder Environmental Fluctuation:Formulation in terms of PhenotypicVariances
Received: date / Accepted: date
Abstract
The characterization of plasticity, robustness, and evolvability, animportant issue in biology, is studied in terms of phenotypic fluctuations. Bynumerically evolving gene regulatory networks, the proportionality betweenthe phenotypic variances of epigenetic and genetic origins is confirmed. Theformer is given by the variance of the phenotypic fluctuation due to noisein the developmental process; and the latter, by the variance of the pheno-typic fluctuation due to genetic mutation. The relationship suggests a linkbetween robustness to noise and to mutation, since robustness can be definedby the sharpness of the distribution of the phenotype. Next, the proportion-ality between the variances is demonstrated to also hold over expressionsof different genes (phenotypic traits) when the system acquires robustnessthrough the evolution. Then, evolution under environmental variation is nu-merically investigated and it is found that both the adaptability to a novelenvironment and the robustness are made compatible when a certain degreeof phenotypic fluctuations exists due to noise. The highest adaptability isachieved at a certain noise level at which the gene expression dynamics arenear the critical state to lose the robustness. Based on our results, we revisitWaddington’s canalization and genetic assimilation with regard to the twotypes of phenotypic fluctuations.
Keywords
Robustness, Fluctuation-Response Relationship, Evolution,Genetic Variance
In evolutionary biology, plasticity and robustness are considered the basiccharacteristics of phenotypes; these characteristics have been widely dis-
Center for Complex Systems Biology and Department of Basic Science, Univ. ofTokyo, 3-8-1 Komaba, Tokyo 153-8902, Japan cussed for decades. In general, phenotypes are shaped from genotypes asa result of developmental dynamics , under an environmental condition. Al-though these dynamics are determined by genes, they can be stochastic innature, owing to some noise in the developmental process; further, these dy-namics also depend on the environmental condition.Plasticity refers to the changeability of the phenotype against the en-vironmental change. Through developmental dynamics, the influence of theenvironment is amplified or reduced[1-6]. In other words, plasticity concernshow the developmental dynamics are affected by the environmental change.Of course, the phenotype depends on the genotype. This changeabilityagainst genetic change is called evolvability and is related to the sensitivity ofdevelopmental dynamics against genetic change. In this sense, both plasticityand evolvability represent the responsiveness against external perturbationsand, in simple terms, a sort of ”susceptibility” in statistical physics.Robustness, on the other hand, is defined as the ability to function againstpossible disturbances in the system [7-13]. Such disturbances have two dis-tinct origins: non-genetic and genetic. The former concerns the robustnessagainst the stochasticity that can arise during the developmental process,while the latter concerns the structural robustness of the phenotype, i.e., itsrigidity against the genetic changes produced by mutations. If the varianceof phenotype owing to disturbances such as noise in gene expression dynam-ics is smaller, then the robustness is increased. In this sense, the varianceof phenotypic fluctuations serves as a measure of robustness, as has alreadybeen discussed[14].Now, there exist certain basic questions associated with robustness andevolution[1-8]. Does robustness increase or decrease through evolution? If itincreases, the rigidity of the phenotype against perturbations also increases.Consequently, the plasticity, as well as evolvability, may decrease with evolu-tion. If that is the case, then the question that arises is: How can the plasticityneeded to cope with a novel environment be sustained?The decrease in plasticity and increase in robustness with evolution wasactually observed in laboratory experiments under fixed environmental andfitness conditions and was also confirmed through numerical experiments.In a simulated evolution of catalytic reaction and gene regulatory networksunder a given single fitness condition, the fluctuation decreases through thecourse of evolution[14,15]. Such networks evolve to reduce the fluctuation bynoise, for a sufficient noise level. Through this evolution, robustness againstnoise increases, leading to the decrease in phenotypic plasticity. However, asa system is more and more adapted to one environment, the phenotype fittedfor it would lose the potentiality to adapt to a novel environmental condition.In the present organisms, however, neither the evolutionary potential northe phenotypic fluctuation vanishes. Even after evolution, the phenotype inquestion is not necessarily fixed at its optimal value, but its variance oftenremains rather large. There can be several sources for the deviation from the Here, ”development” refers to a dynamic process that shapes the phenotypeand is not restricted to multicellular organisms. Developmental dynamics are alsoobserved in unicellular organisms, for instance, gene (protein) expression dynamicsby regulatory networks. optimal state, which is neglected in the idealized numerical and laboratoryexperiments with a single fitness condition. One of the most typical factorsfor this deviation is environmental variation. With environmental change, thephenotype for producing the fittest state is not fixed but may vary over gen-erations. Then, we are faced with the following questions: First, under envi-ronmental variation, can robustness and plasticity coexist through the courseof evolution? Second, are the phenotypic fluctuations sustained, to maintainthe adaptability to environmental changes? Third, how are the opposing fea-tures of robustness and adaptability to the new environment compromised?Finally, is there a noise level optimal for achieving both adaptability androbustness? We address these questions in the present paper.
As for the fluctuation, there is an established relationship between the evolu-tionary speed and the variance of the fitness, namely, the so-called fundamen-tal theorem of natural selection proposed by Fisher[19,20,21]; this theoremstates that the evolution speed is proportional to V g , the variance of fit-ness due to genetic variation. Besides the fitness itself, any phenotypes aregenerally changed by the genetic variation. Now, relationship between theevolution speed of each phenotype with its variance due to genetic change isestablished as Price equation, formulated through the covariance between aphenotype and fitness[22,23]. Statistical-mechanical interpretation of Fisher’stheorem in terms of the fluctuation-response relationship[24] or fluctuationtheorem[25] was also discussed.Besides the genetic change, there are sources of phenotypic fluctuationseven among isogenic individuals, and under a fixed environment. In fact,recent observations in cell biology show that there are relatively large fluctu-ations even among isogenic individuals. The protein abundances over isogenicindividual cells exhibit a rather broad distribution, i.e., the concentration ofmolecules exhibits quite a large variance over isogenic cells [26,27,28,29]. Thisvariance is due to stochastic gene expressions and to other external pertur-bations. Furthermore some of such phenotypic fluctuations are indeed tightlycorrelated with the fitness, and hence the fitness also shows (relatively large)isogenic fluctuations. For example, large fluctuations in the growth speed (ordivision time) are observed among isogenic bacterial cells[30,31].Hence, there are two sources of variations for the fitness or phenotype: ge-netic and epigenetic. The former concerns structural change in developmentaldynamics with genetic change, and the latter concerns the noise during thegene expression dynamics. As mentioned, the relationship between the for-mer variance with the evolution speed was established as Fisher’s theorem(for fitness) and Price equation (for general phenotype). Then, is there anyrelationship between the latter variance with the evolution speed?At a glance, the latter variance might not seem to be relevant to evo-lution, since this epigenetic change itself due to noise is not inherited tothe offspring, in general. This is not the case. Here, it should be noted thatthe degree of variance itself is a nature of developmental dynamics to shape the phenotype, and thus depends on genotype, and can be inherited. Hence,there may exist some correlation between evolution speed and this isogenicphenotypic variance, which is denoted as V ip here .Indeed, from evolution experiments for increasing the fluorescence in aninserted protein in bacteria[17] and also from numerical experiments forevolving reaction networks or gene regulatory networks[15,14] to increasea given fitness, we have observed that V ip ∝ Evolution Speed,where the evolution speed is defined as the increase of the fitness. In theexperiment, for each mutant bacteria, average fluorescence of isogenic cellswas measured and that with highest average fluorescence was selected. Si-multaneously, the variance of isogenic bacterial population was measuredto obtain V ip . The same procedure was applied to measure the evolutionspeed and the fitness variance in numerical evolution to increase a givenfitness. Interestingly, these experiments support the above relationship, atleast approximately. Some ’phenomenological explanation’ was proposed byassuming a Gaussian-type distribution P ( x ; a ) of the fitness (phenotype) x as parameterized by a ”genetic” parameter, a . and a linear change in thepeak position of x against a [16,17,18]. It should be noted, however, that theargument based on the distribution P ( x ; a ) is not a ”derivation” but rathera phenomenological description. Indeed, the description by P ( x ; a ) itself isan assumption: for example, whether a genotype is represented by a scalarparameter a is an assumption.Now, considering the established Fisher’s relationship, the above relation-ship suggests the proportionality between V ip and V g through an evolutionarycourse. Indeed, this relationship was confirmed from several simulations ofmodels. Again, this relationship is not derived from established relationshipsin population genetics. Indeed, the proportionality between the two is notobserved in the first few generations, but observed after robust evolutionpreserving a single peak in the fitness is progressed. Considering this obser-vation, we previously discussed the relationship, by postulating evolutionarystability of the distribution over phenotype and genotype[14,15,35].The relationship between V ip and V g also suggests a possible link betweendevelopmental robustness to noise and evolutionary robustness to geneticchanges (mutation). For this link, we first note that the two types of variances V g and V ip lead to two kinds of robustness: rigidity of the fitness (phenotype)against genetic changes produced by mutations and robustness against thestochasticity that can arise in an environment or during the developmentalprocess. When the fitness (phenotype) is robust to noise in developmentalprocess, it is rather insensitive to the noise, and therefore, its distributionis sharper. Hence, the (inverse of the) variance of isogenic distribution, V ip ,gives an index for robustness to noise in developmental dynamics[14,15]. Onthe other hand, if V g is smaller, the phenotype is rather insensitive to ge- This is not standard terminology. Phenotypic variance by non-genetic originsis often termed as environmental variation V e . Here, however, the source of thevariance is not necessarily environmental change but, primarily, noise in the devel-opmental process; hence, we use a different term. Note that the total phenotypicvariance is V ip + V g under a certain ideal condition[20,21], if they are added inde-pendently. netic changes, implying higher genetic (or mutational) robustness. Hence thecorrelation between V ip and V g implies correlation between the two typesof robustness. Indeed, previous simulations suggest that robustness to noisefosters robustness to mutation. Congruence between evolutionary and de-velopmental robustness was also discussed by Ancel and Fontana for theevolution of RNA[5].To close this section, a brief remark on the measurement of V g is givenhere. Because the fitness (phenotype) distribution exists even in isogenicindividuals, the variance of its distribution over a heterogenic populationincludes both the variance among isogenic individuals and that due to geneticvariation. To distinguish between the two contributions, we first measure theaverage fitness (phenotype) over isogenic individuals and then compute thevariance of this average over a heterogenic population. This variance will onlybe attributed to genetic heterogeneity. This variance is denoted by V g . F itness = f ( Genotype ). Here the evolu-tion is discussed as a hill-climbing process through random change in geno-type space and selection. Although this viewpoint has been important inevolutionary studies, another facet in evolution has to be also considered,to discuss the phenotypic plasticity and evolution-development relationship,that is genotype-phenotype mapping shaped by the developmental dynamics.Here we focus on this evolution-development relationship.Note that the fitness is a function of (some) phenotypes (say a set ofprotein abundances). The phenotypes are determined by developmental dy-namics (say gene expression dynamics), whose rule (e.g., set of equations orparameters therein) is governed by genes. Thus this evolution-developmentscheme is represented as follows:(i)
Fitness= F(phenotype) (ii)
Phenotype determined as a result of developmental dynamics (iii)
Rule of developmental dynamics given by genotype According to (ii) and (iii), genotype-phenotype mapping is shaped, whichis not necessarily deterministic. As discussed, the phenotypes from isogenicdistribution are distributed, since the dynamics (ii) involve stochasticity (ingene expression). Indeed, with this stochasticty, reached attractors by thedynamics are not necessarily identical, and accordingly the fitness may varyeven among the isogenic individuals, which lead to the variance V ip . If thedynamics to shape the phenotype is very stable, the variance V ip is smaller. According to the conventional terminology in population genetics, this varianceis referred to as ”additive” genetic variance; The term ”additive” is included, toremove the variance due to sexual recombination. Here we do not discuss the influ-ence of recombination and, thus, do not need to distinguish genetic variance fromthe additive one.
Further, the degree of the change in phenotype by the genetic change dependsboth on (ii) and (iii), and V g depends on the sensitivity of the phenotype tothe change in the rule of the dynamics. With the evolutionary process thegenotype, i.e., the rule of the dynamics, changes, so that V ip and V g changewith the evolution.The developmental dynamics, in general, involve a large number of vari-ables (e.g., proteins expressed by genes), and are complex. Accordingly, genotype-phenotype mapping is generally complex, and the mapping from genotypeto the fitness is not simple, even if the function in (i) is simple. In severalstudies with adaptive fitness landscape, complexity in the fitness landscapehas been taken into account, while the developmental dynamics (ii) are not.We take a simple fitness function (say, the number of expressed genes), butinstead take complex developmental dynamics into account, to discuss thephenotypic plasticity in developmental and evolutionary dynamics.If we are concerned only with the fitness landscape (i), we can introducea single ’energy’-like function for it and formulate the adaptive evolution interms of standard statistical mechanics. Here, however, we need to considerthe dynamics to shape the phenotype (ii) [5,13,14,15]. If we adopt statistical-mechanical formulation, we need two ’energy’-like functions, one for fitnessand the other for Hamiltonian for the development dynamics, as is formulatedby two-temperature statistical physics[33,34].So far, such study on evolution-development (so called ’evo-devo’) lacksmathematically sophisticated formulation as compared with the celebratedpopulation genetics, and thus we sometimes have to resort to heuristic studiesbased on numerical evolution experiments, from which we extract ’generic’properties. We then provide some plausible arguments (or phenomenologi-cal theory), while mathematical or statistical-physical formulation has to bepursued in future.3.2 Gene expression network modelFollowing the argument in the lase section and to discuss the issue of plas-ticity and robustness in genotype-phenotype mapping, we adopted a simplemodel[14,39,40] for gene expression dynamics with a sigmoid input-outputbehavior[36,37,38]; it should be noted, however, that several other simula-tions in the form of other biological networks will give essentially the sameresult. In this model, the dynamics of a given gene expression level, x i , isdescribed by the following: dx i /dt = γ { f ( M X j J ij x j ) − x i } + ση i ( t ) , (1)where J ij = − , ,
0. The noise term η i ( t ) is due to stochasticity in geneexpression. For simplicity, it is taken to be Gaussian white noise given by < η i ( t ) η j ( t ′ ) > = δ i,j δ ( t − t ′ ), while the (qualitative) results to be discussed below is independent of the choice . The amplitude of noise strength is givenby σ , which determines stochasticity in gene expression. M denotes the totalnumber of genes; and k , the number of target genes that determine fitness.Here, the function f ( x ) represents a threshold function for gene expression.Previously, we adopted the model (1) with f ( x ) = tanh( βx ) , (2)where the gene expression is ”on” if x > x < f ( x ) = 1 / ( exp ( − β ( x − θ i ) + 1) + δ (3)Here, x is positive and is scaled so that the maximal level expression is ∼ δ takes a small positive value corresponding to a spontaneous expres-sion level. If the input sum from other genes P Mj J ij x j to gene i exceeds thethreshold θ i , the gene (protein) i is expressed so that x i ∼
1, as f approachesa step function for sufficiently large β ( > x i ’s are initially smaller than the threshold θ i , they remains so if δ < θ i . Hence, to generate some expression patterns,some input term is needed. Here, we introduce ”input” genes j = 1 , ...k inp ,in which x j is set at x j = I j > j = 1 , , ..., k inp ) , (4)where I j is an external input to a set of ”input genes.” These inputs areneeded to express genes; otherwise, the expression levels remain sub-threshold.In this case with eq. (3), the gene is expressed if x i > θ i , The initial conditionis given by a state where none of the genes are expressed, i.e., x i ∼ F , is determined bywhether the expressions of the given ”target” genes are expressed after a suf-ficient time. This fitness condition is given such that k target genes are ”on”(expressed), i.e., x i > θ i for M − k < i ≤ M . Because the model includes anoise component, the fitness can fluctuate at each run, which leads to a distri-bution in F and x i , even among individuals sharing the same gene regulatorynetwork. For each network, we compute the average fitness F over L runs aswell as the variance. This variance over L identical individuals (having theidentical networks J ij ) due to noise is V ip , as it represents the fluctuation ofthe fitness over isogenic individuals (isogenic phenotypic variance). V g , to bediscussed below, is the variance of F over N heterogenic individuals (havingdifferent networks J ij ). Indeed, multiplicative noise depending on x i , as well as a stochastic reactionmodel simulated with the use of Gillespie algorithm gives a qualitatively samebehavior[41]. Of course the magnitude of the phenotype variance as well as thethreshold noise level for error catastrophe, to be discussed below, depends on theform of noise. However, relationship of such threshold noise level with V g /V ip , aswell as the proportionality between V ip and V g does not depend on the specificform of noise. Here V ip is not the variance over time in a single run, but the varianceover runs for identical individuals. With developmental process by gene ex-pression dynamics from a given initial condition, the gene expression pattern( x ( j ); j = 1 , .., M ) reaches an attractor. Reached attractors can be differ-ent by each run, due to the stochasticity in expression dynamics during thetransient time steps to reach the attractor. The fitness is computed on thisattractor, which is also distributed as a result of noise. Once the attractoris reached, the fluctuation of the fitness around it is negligible, since the fit-ness is not given directly by x i but defined through a threshold function of x i − θ i . The variance is smaller as more runs result in the same attractor (orattractors of the same fitness) under noise. If gene expression dynamics withsuch global attraction to the same-fitness attractor are shaped through theevolution, the variance V ip is smaller.Now, at each generation, there are N individuals with slightly differentgene regulatory networks { J ij } , and their average fitness F may differ byeach. Among the networks, we select the ones with higher fitness values.From each of the selected networks, J ij is ”mutated,” i.e., J ij for a certainpair i, j selected randomly with a certain fraction switches to one of the valuesamong ± , N s ( < N ) networks with higher F values are selected, each ofwhich produces N/N s mutants. We repeat this selection-mutation processover generations. (For example, we choose N = L = 500 or 700 for mostsimulations, and N s = N/
4; the conclusion, to be shown below, does notchange as long as these values are sufficiently large. We use β = 7, γ = . M = 64, and k = 8 and initially choose J ij randomly. The probability istaken to be 1/4 for J ij = ± V ip , V g is defined as thevariance of F over the N individuals having different genes (gene regulatorynetworks). This shows the variance due to genetic change. As V g is decreased,the fitness becomes insensitive to the genetic change, i.e., the robustness tomutation is increased. On the other hand, as V ip is decreased, the robustnessto noise is increased, since V ip measures the variance due to noise in a de-velopmental process. (It should be remarked that the absolute value of thefitness F is not important, since the model behavior is invariant against thetransformation F → c × F . Accordingly the absolute values of V ip and V g are not important, but the ratio between the two or relative change of thevariances over generations is important. On the other hand, each expression x ( i ) is already scaled so that the maximal expression is ∼ σ c beyond which the evolutionof robustness progresses, so that both V g and V ip decrease. Here, most of the -0.06-0.05-0.04-0.03-0.02-0.01 0 0 20 40 60 80 100 120 140 A v e r age f i t ne ss Generation σ =.06 σ =.03 σ =.005 Fig. 1
Evolutionary time course of the average fitness < F > . The average ofmean fitness < F > over all N = 500 individuals that have different genotypes(i.e., networks J ij ) in each generation is plotted. The mean fitness of each genotypeis computed from L = 500 runs. The plotted points are for different values ofnoise strength, σ =0.005, 0.03, 0.06, with different colors. For Figs.1-3, we choose M = 64, k = 8, and k inp = 8 with I j = 1, while θ i is distributed in [0.1,.0.3]. individuals take the highest fitness value. In contrast, for a lower noise level σ < σ c , mutants that have very low fitness values always remain. Severalindividuals take the highest fitness value, whereas the fraction of individualswith much lower fitness values does not decrease; hence, V g remains large(see Figs.1 and 2).(2) At around the threshold noise level, V g approaches V ip . For σ < σ c , V g ∼ V ip holds, whereas for σ > σ c , V ip > V g is satisfied. For robust evolutionto progress, this inequality is satisfied.(3) When the noise is larger than this threshold, the two variances de-crease, while V g ∝ V ip is maintained through the evolution course. Hence,the proportionality between the two variances is confirmed.Why does the system not maintain the highest fitness state under a smallphenotypic noise level with σ < σ c ? Indeed, the dynamics of the top-fitnessnetworks that evolved under such low noise levels have distinguishable fea-tures from those that evolved under high noise levels. It was found thatfor networks evolved under σ > σ c , a large portion of the initial conditionsreached attractors that give the highest fitness values, whereas for networksevolved under σ < σ c , only a tiny fraction (i.e., in the vicinity of the all-offstates) reached such attractors.In other words, for σ > σ c , the ”developmental” dynamics that give afunctional phenotype have a global, smooth attraction to the target. In fact,such types of developmental dynamics with global attraction are known to beubiquitous in protein folding dynamics [42,44], gene expression dynamics[43], V g Vip σ =.06 σ =.05 σ =.03 σ =.005Vg=Vip Fig. 2
Relationship between V g and V ip , the variances of the fitness. V g is com-puted from P ( F ), the distribution of mean fitness F and V ip is computed from theisogenic variance of the fitness over L = 500 runs, for each genotype, and then thesevalues are averaged over all existing individuals. (We also confirmed the overall re-lationship by using the variance for a gene regulatory network that gives the peakfitness value in P ( F )). The plotted points are over 140 generations. σ = 0 .
005 ( (cid:3) ),0.03 ( ∗ ), 0.05 ( × ) and 0.06 (+). For σ > σ c ≈ .
02, both the variances decreasewith generatons, so that the right-top is the first generation and left-bottom is the140th generation. After initial few generations both decrease roughly in proportion,while at σ = 0 .
03, the deviation from the proportionality is larger possibly becausethe noise level is near the critical point to lose the robust evolutionary process. For σ = 0 . and so forth. On the other hand, the landscape evolved at σ < σ c is rugged.Except for the vicinity of the given initial conditions, the expression dynamicsdo not reach the target pattern.The observed proportionality between V ip and V g is not self-evident. In-deed, if a random network is considered for a gene regulatory network, suchproportionality is not observed. In the present simulation, after a few gen-erations of evolution, both the variances decrease, following proportionality,if σ > σ c . Although there is no complete derivation for this relationship, itis suggested that this proportionality as well as the relationship V ip > V g is a consequence of evolutionary stability to keep a single-peakedness in thedistribution P ( x = phenotype, a = genotype ), under conditions of strongselective pressure and low mutation rate[15,35].4.2 Proportionality between the two variances across genesAs mentioned above, the gene expression dynamics evolved under a sufficientlevel of noise have a characteristic property; the attractor providing the phe- notype of the highest fitness has a large basin volume and is, hence, attractedglobally by a developmental process under noise.Note that the expression level x j of non-target genes j could be either onor off, because there is no selection pressure directed at fixing their expressionlevel. Still, each expression level x j can have some correlation with the fitnessin general. Hence it is also interesting to study the variance of each expressionlevel and discuss its evolutionary changes.Similar to the variances for the fitness, the phenotypic variance V ip ( i ) foreach gene i in an isogenic population is defined on the basis of the varianceof the expression of each gene i , with each X i = Sign ( x i − θ i ), in an isogenicpopulation. Accordingly, the variance computed by using the distribution of X i in this heterogenic population gives V g ( i ) for each gene i .Following Price equation[22], the rate of change in each expression levelbetween generations is expected to be correlated with V g ( i ), as it has director indirect influence to the fitness, through the gene expression dynamics (1).In contrast, here, we are interested in the variance of isogenic phenotypic fluc-tuations of each expression level V ip ( i ). Indeed, this variance decreases overgenerations for most genes. As the evolution progresses, these expressionsalso start to be rigidly fixed so that their variances decrease over most genes.In Fig.3(a), we have plotted V g ( i ) versus V ip ( i ) over generations for severalgenes i . We can see that they decrease (roughly) in proportion, over gen-erations. Furthermore, the proportion coefficient V g ( i ) /V ip ( i ) seems to takeclose values across many genes.In Fig.3(b) we have plotted ( V ip ( i ), V g ( i )), across all expressed genes,after evolution reached the genotype with the highest fitness. As shown,the proportionality (or strong correlation) between V g ( i ) and V ip ( i ) holdsacross many (expressed) genes for a system through evolution. The ratio ρ = V g ( i ) /V ip ( i ) increases with the decrease in σ , and at around σ c , it approaches ∼ P ( x i , a )for each gene expression x i and by further assuming that the distributionmaintains a single peak up to a common mutation rate (i.e., a common mu-tation rate for the error catastrophe threshold) over genes[40]. The latterhypothesis may be justified as a highly robust system: In such a system withincreased error-threshold mutation rates, once the error occurs, it propagatesand percolates to many genes, so that a common error threshold value is ex-pected. Note that there are some preliminary experimental supports on thisproportionality of the two variances over genes or phenotypic traits[18,46,47,48], although future studies are required for the confirmation. Throughout the present paper, V g ( i ) and V ip ( i ) with ( i ) denote the variancesof each phenotype, expression level i , while, V g and V ip denote the variances of thefitness.2 V g ( i ) Vip(i) i=16i=21i=22i=31i=52Vg=Vip 1e-06 1e-05 0.0001 0.001 0.01 0.1 1 0.001 0.01 0.1 V g Vip σ =.06 σ =.03 σ =.005Vg=Vip Fig. 3 (a) Plot of ( V ip ( i ) , V g ( i )) for the genes i = 16 , , , ,
52 for the gener-ations 5-40. Each of the variances decreases over generations, so that variance foreach gene changes from right (upper) to left (lower) in the figure. As describedin the text, V ip ( i ) was computed as the variance of the distribution of Sign ( x i )over L = 500 runs for an identical genotype, while V g ( i ) was computed as a vari-ance of the distribution of ( Sign ( x i )) over N = 500 individuals, where Sign ( x i )refers to the mean over 500 runs. With generations, both the variances decreaseroughly with proportion, with a trend of common proportion coefficient. (b) Plotof ( V ip ( i ) , V g ( i )) across all expressed genes i (i.e. for such genes i that x ( i ) > θ i ),after evolution is completed (for the generations 25-60), for three values of noiselevels σ = 0 .
005 (*), 0.03 ( × ), and 0.06 (+).3 / environmental con-dition, both the fluctuations and the rate of evolution decrease. The systemloses plasticity against the change caused by external noise or external mu-tation. Nevertheless, in nature, neither the fluctuations nor the evolutionpotential vanish. How are phenotypic plasticity, fluctuations, and evolution-ary potential sustained in nature?One possible origin for the preservation of plasticity may be environmen-tal fluctuation[51], as has also been studied in terms of statistical physics[52,53,54]. The plasticity of a biological system is relevant for coping with the en-vironmental change that may alter phenotypic dynamics in order to achievea higher fitness. In the present modeling, there can be two ways to includesuch an environmental change.One is a direct method, in which an input term given in eq.(4) is changedwith each generation, while the fitness condition is maintained. The othermethod is indirect, in which environmental change is introduced as the changein the fitness condition, while preserving the dynamics itself. Here, we discussthe simulation result of the former procedure first and will consider the resultfrom the other procedure later in § I j (1 ≤ j ≤ k inp ) fromthe generation at which the system had already adapted to the environmentand decreased the phenotypic variances. An example is plotted in Fig.4. Here, I j initially takes I j = 1 for 1 ≤ j ≤ k inp / I j = 0 for k inp / < j ≤ k inp ,before switching to 1 − I j after the 29th generation. By switching the envi-ronment, the fitness first decreases and later adapts to the new environment(Fig.4a).To determine the evolution of phenotypic plasticity, we computed thevariances of the fitness, V ip and V g , over successive generations (see Fig.4b).After the switch of the fitness condition, both V ip and V g first increase to arelatively high level and continuously increasing further over a few genera-tions. At later generations, both V ip and V g again decrease, maintaining theproportionality. The proportionality law between the genetic and epigeneticvariances is satisfied with both increase and decrease in plasticity throughthe evolution.With the increase in V ip , the fitness is more variable with noise, which alsoleads to higher changeability against environmental conditions. The gene ex-pression dynamics regain plasticity, which allows for the switch of the targetgenes after further generations. Then, with the increase in V g , the change-ability against genetic change increases, thus increasing the evolvability.Next, we explore the change in the variances of the expression of eachgene. As shown in Fig.5, both variances V ip ( i ) and V g ( i ) increase, as a resultof environmental change. Expressions of all genes are more variable againstnoise and mutation. Most gene expressions gain higher plasticity by increas- -2-1.5-1-0.5 0 0.5 1 0 10 20 30 40 50 60 A v e r age f i t ne ss , V i p , V g GenerationAverage fitnessV ip V g
20 30 40 50 60 70 80Vip V g Switch after 29th Vg=Vip 0.0001 0.001 0.01 0.1 1 1e-05 0.0001 0.001 0.01 0.1
Fig. 4
The time course in the fitness and the variance of the fitness over genera-tions. First, the evolution under an environment I j = 1 for 1 ≤ j ≤
4, and I j = 0for 5 ≤ j ≤ k inp = 8 is simulated to progress up to 30 generations. After the 29thgeneration, we switch the fitness condition to I j = 0 for 1 ≤ j ≤
4, and I j = 1 for5 ≤ j ≤
8, which is maintained at later generations. M = 64, L = N = 700, and θ i is distributed uniformly in [0.1, 0.3]. The switch initially causes a decrease in thefitness, but after a few dozens of generations, almost all networks evolve to adaptto the new fitness condition. The noise level is set at σ = 0 . > σ c . Top: Thetime course of the average fitness and the variance V ip throughout the evolution.Bottom: The plot of the variances of the fitness, V g versus V ip for each generationafter the switch of the environmental condition. The generations (up to 60) arerepresented by different colors. Both the variances increase in correlation, after theswitch, and later, they decrease in proportion, to adapt to the new condition.5 V g Fig. 5
Change in the variances V ip ( i ) and V g ( i ) /V ip ( i ) at after the switch in en-vironmental condition after the 29th generation, as given in Fig. 4. The variance V ip ( i )) is plotted up to generation 60, where the switch is given at the 30th genera-tion. Plotted for the genes i =13,15,28,40,47,51. For the genes 13, 28, and 47, V g ( i )before the switch after 29th generation is smaller than 10 − . The color represents V g ( i ) /V ip ( i ), as shown in the right bar. ing the variance in their expression. Besides the increase in the variances,the ratio ρ i = V g ( i ) /V ip ( i ) also increase at some generation, to approach V ip ∼ V g (see the color change in Fig.5, where the red color shows higherratio ρ i . This increase in V g ( i ) /V ip ( i ) suggests that the system is closer tothe error catastrophe point, where the stability condition in the distributionfunction P ( x i = phenotype, a = genotype ) is lost. This leads to the increasein plasticity, with which the adaptation to a novel condition is achieved.Once the fitted phenotypes are generated by this adaptation, the variancesdecrease, with restoring the proportionality between V g ( i ) and V ip ( i ).To sum up, adaptation to novel environment is characterized by the phe-notypic variances as follows. With the increase in V ip ( i ), the sensitivity of thephenotype to noise and environmental change is increased, thereby increasingplasticity. With the increase in V g ( i ), changeability of the phenotype by mu-tation is increased, thereby accelerating evolutionary change of each expres-sion level. With the increase in V g ( i ) /V ip ( i ), the sensitivity to genetic changeis further increased, thus facilitating the evolution. With these trends—theincrease in V ip ( i ), V g ( i ), ρ i = V g ( i ) /V ip ( i ), the adaptation to a novel environ-ment is fostered. With this increase in plasticity, gene expression dynamicsfor adapting to a novel condition are explored. Once these fitted dynamicsare shaped, the variances decrease, leading to a decrease in plasticity andincrease in robustness. -1.2-1-0.8-0.6-0.4-0.2 0 0 50 100 150 200 250 300 350 400 A v e r age f i t ne ss Generation σ =.1 σ =.05 σ =.01 σ =.001 Fig. 6
The average of the mean fitness < F > plotted for each generation, undercontinual environmental variation, as described in the text. M = 64, k = 8, and k inp = 4, while I j are changed randomly within [0,.8]. The average of the meanfitness, F , of each individual (over L = 100 runs) is computed over the totalpopulation ( N = 100) at each generation. The noise level, σ , is 0.1 (red), 0.05( ∼ σ c ; green), 0.01(blue), and 0.001 (pink). At around σ ∼ .
05, the average fitnessreaches the highest level. I i within [0 ,
1] pergeneration.When environmental changes are continuously repeated, the decrease andincrease in the variances V ip and V g are repeated. Note that it takes moregenerations to adapt to a new fitness condition, if the phenotypic varianceshave been smaller. In our model, if the noise level in development is larger,the phenotypic variances already take a small value during the adaptationto satisfy the fitness condition. Hence, in this case, it takes more generationsto adapt to a new fitness condition. On the other hand, if the noise level σ is smaller than σ c ∼ .
05, robust evolution does not progress. Hence, forcontinuous environmental change, there will be an optimal noise level toboth adapt sufficiently fast to a new environment and evolve the robustnessof fitness for each environmental condition. In Fig.6, we have plotted thetime course of the average fitness in population. If the noise level is large,the system cannot follow the frequent environmental change and the averagefitness cannot increase sufficiently. On the other hand, if the noise level issmall, the fitness increases; however, if it is too small, the fitness of someindividuals remains rather low. Indeed, there is an optimal noise level atwhich the average fitness is maximal, as shown in Fig.7, where the average - A v e r age F i t ne ss , V i p , V g σ -Average FitnessVgVip Fig. 7
The average fitness and the variances, V ip and V g , through the course of theevolution under environmental variation as in Fig.6. Each value is further averagedover 200-700 generations. The overall temporal averages are plotted against thenoise level σ . Instead of the fitness itself, its sign inversion − < F > is plotted. Ataround σ ∼ .
05, the average fitness takes a maximal value, and at σ slightly belowit, V g approaches V ip . V g Vip σ =.1 σ =.05 σ =.01 σ =.001Vg=Vip Fig. 8
The variances of the fitness ( V ip , V g ) through the course of the evolutionunder environmental variation as in Fig.6. Each point is a result of one generation,and the plot is taken over 200–700 generations. The noise level σ is 0.1 (red), 0.005( ∼ σ c ; green), 0.01 (blue), and 0.001(pink).8 fitness over generations is plotted against the noise level σ . This optimal noiselevel is close to the value of the robustness transition σ c .Next, we plotted the variances V ip and V g over generations (see Fig.8).When σ < σ c , then V g > V ip and both the variances remain rather large,demonstrating that robustness has not evolved at all. For σ ≫ σ c , V g < V ip and the variances remain small. The robustness has evolved, but the systemcannot adapt to an environmental change as the variances have become toosmall. In contrast, for σ ∼ σ c , V ip and V g vary between low and high valuesover generations, maintaining the proportionality between the two variances,with V ip slightly larger than V g .5.3 Adaptation against switches of the fitness condition -4.5-4-3.5-3-2.5-2-1.5-1-0.5 0 400 450 500 550 600 F i t ne ss Generation
Fig. 9
The average of fitness plotted per generation, where the fitness conditionin the model (1)-(2) is switched for every 10 generations between + + + + + + ++and + + + + − − −− . The average of the mean fitness F of each individual (over L = 200 runs) is computed over the total population ( N = 200) at each generation.The noise level σ is 0.1 (red), 0.008 ( ∼ σ c ; green) and 0.001 (blue). For confirmation of the result in the last section, we also carried outnumerical experiments by adopting a separate procedure, i.e., by switchingthe fitness condition. As a specific example, we carried out the simulation bytaking the model (1)-(2), without including input terms (4). After the geneexpression dynamics are evolved with the fitness to prefer x i > i = 1 , , ...k (= 8) adopted already, then at a certain generation,we change the fitness condition so that the genes i = 1 , , .., k/ − − −− ,instead of + + + + + + ++: In this model, the gene is off if x i < x i > - A v e r age F i t ne ss , V i p , V g σ -Average FitnessVipVg Fig. 10
The temporal average of the average fitness and the variances V ip and V g in the evolution, with the change in the fitness condition for every 10 generationsas in Fig.9. The overall temporal averages over population and over generationsare plotted against the noise level σ . Instead of the fitness itself, its sign inversion − < F > is plotted. The temporal average is taken over 500–1000 generations. Ataround σ ∼ . σ slightly belowit, V g exceeds V ip . V g Vip σ =.1 σ =.002 σ =.008 σ =.001Vg=Vip Fig. 11
Variances of the fitness, ( V ip , V g ), are plotted over generations through thecourse of the evolution with fitness change after every 10 generations, as describedin Fig.9. Each point is a result of one generation, and the plot is taken over 500–1000 generations. The noise level σ is 0.1 (red), 0.02(green), .008 ( ∼ σ c ;blue), and0.001 (pink). In this case as well, the variances of the fitness, V ip and V g , first increasein proportion to adapt to a new fitness condition. Later, they decrease in pro-portion to gain robustness to noise and mutation. Next, we again computed V ip ( i ) and V g ( i ) successively through the course of the evolution. Immedi-ately after the switch in the fitness, the variances V ip ( i ) and V g ( i ) increase as well as the ratio V g ( i ) /V ip ( i ). These variances in gene expression levelsfacilitate plasticity, and adaptation to a new environment.When environmental changes are continuously repeated, the decrease andincrease processes of the variances V ip and V g are repeated. In Fig.9, we haveplotted the time course of the average fitness in population, when the fitnesscondition is switched after every 10 generations. In this case again, whenthe noise level is near σ c , fast adaptation to a new environment and theincrease of robustness at later generations are compatible. Dependence ofthe average fitness on the noise level is shown in Fig.10, which also showsan optimal noise level near the robustness transition σ c (which is lower thanthe case in § V g exceeds V ip and therobustness to mutation is lost. The plot of the variances V ip and V g overgenerations in Fig.11 shows that at σ ∼ σ c , they go up and down, maintainingan approximate proportionality between the two, with V ip slightly less than V g . Overall, the behavior here under the fitness switch agrees well with thatunder the environmental variation in § In the present paper, we have studied biological robustness and plasticity,in terms of phenotypic fluctuations. The results are summarized into threepoints.(1) Confirmation of earlier results on the phenotypic variances due tonoise and due to genetic variation: The two variances, V ip (due to noise)and V g (due to mutation) decrease in proportion through the course of ro-bust evolution under a fixed environmental condition. After the evolution toachieve robustness to noise and mutation is completed, V ip ( i ) and V g ( i ) areproportional across expressions of most genes (or different phenotypic traits).In short, plasticity (changeability) of phenotype ∝ V ip ∝ V g ∝ evolutionspeed through the course of the evolution and across phenotypic traits (expressionsof genes).(2) Increase in the phenotypic variances and recovery of plasticity: Whenrobustness is increased under a given environmental condition, the systemloses plasticity to adapt to a novel environment. When the environmentalcondition is switched, both the phenotypic variances due to noise and due togenetic variation increase to gain plasticity and thus to adapt to the novelenvironment. This increase is observed both for the variance of fitness andof each gene expression level.(3) Optimal noise level to achieve both robustness and plasticity undercontinuous environmental change. There is generally a threshold level of noisein gene expression (or in developmental dynamics), beyond which robustnessto noise and mutation evolves. If the noise level is larger, however, the systemloses plasticity to adapt to environmental changes, whereas if it is much lower,a robust, fitted phenotype is not generated. At around the noise level for the”robustness transition,” the system can adapt to environmental changes andachieve a higher fitness. There, the phenotypic variances V ip and V g increase and decrease, roughly maintaining the proportionality between the two, whilesustaining V ip > ∼ V g .From a statistical-physicist viewpoint, the relationship between respon-siveness and fluctuation is expected as a proper extension of the fluctuation-response relationship. In this context, it is rather natural that the response ra-tio to environmental change, i.e., plasticity, is proportional to V ip , the isogenicphenotypic fluctuation. On the other hand, as Fisher stated, the variancedue to genetic change, V g , is proportional to evolution speed, i.e., responseof the fitness against mutation and selection. Interestingly, our simulationsand evolutionary stability argument suggest the proportionality between V ip and V g . This implies the proportionality between environmental plasticityand evolvability.In fact, Waddington[49,50] coined the term genetic assimilation, in whichphenotypic changes induced by environmental changes foster later geneticevolution. Since then, positive roles of phenotypic plasticity in evolution havebeen extensively discussed[3,4,5]. Our study gives a quantitative representa-tion of such relationship in terms of fluctuations.Existence of the threshold noise level below which robustness is lost isreminiscent of a glass transition in physics: For a higher noise level, dynam-ical systems for global attraction to a functional phenotype are generatedthrough evolution, whereas for a lower noise level, the dynamics follow mo-tion in a rugged landscape, where perturbation to it leads to a failure inthe shaping of the functional phenotype. In fact, Sakata et al. considered aspin-glass model whose interaction matrix evolves to generate a high-fitnessthermodynamic state. The transition to lose robustness was found by lower-ing the temperature. Interestingly, this transition is identified as the replicasymmetry breaking transition from a replica symmetric phase[33,34].Under environmental fluctuation, the evolution to achieve both plasticityto a new environment and robustness of a fitted state is possible near thistransition, for losing the robustness. In other words, one may regard that abiological systems favors the ”edge-of-glass” state. Acknowledgements
I would like to thank T. Yomo, C. Furusawa, M. Tachikawa, A. Sakata,K. Hukushima, L. Lafuerza and S. Ishihara for stimulating discussion. Thiswork was supported by a Grant-in-Aid for Scientific Research (No 21120004)on Innovative Areas “The study on the neural dynamics for understandingcommunication in terms of complex hetero systems (No.4103)” of MEXT,Japan.
References
1. de Visser JA, et al. (2003) Evolution and detection of genetic robustness.
Evolution :1959-1972.2. Callahan HS, Pigliucci M., and Schlichting CD (1997) Developmental phe-notypic plasticity: where ecology and evolution meet molecular biology. Bioessays :519-525.23. West-Eberhard MJ (2003) Developmental Plasticity and Evolution (OxfordUniversity Press).4. Kirschner MW, and Gerhart JC (2005),
The Plausibility of Life (Yale Uni-versity Press).5. Ancel LW, and Fontana W (2002) Plasticity, evolvability, and modularity inRNA.
J. Exp. Zool. :242-283.6. Frank, S.A., (2012)
Natural selection. II. Developmental variability and evo-lutionary rate* , Journal of Evolutionary Biology , , 2310–2320,7. Wagner A (2002) Robustness against mutations in genetic networks of yeast. Nature Genetics :355-361.8. Wagner A. 2005, Robustness and evolvability in living systems , PrincetonUniversity Press, Princeton NJ.9. Wagner GP, Booth G, Bagheri-Chaichian H, 1997, A population genetic the-ory of canalization.
Evolution :329-347.10. Barkai N and Leibler S (1997) Robustness in simple biochemical networks. Nature , :913-917.11. Alon U, Surette MG, Barkai N, Leibler S (1999) Robustness in bacterialchemotaxis. Nature , 1999, :168-171.12. Siegal ML and Bergman A (2002) Waddington’s canalization revisited: Devel-opmental stability and evolution.
Proc Nat. Acad. Sci. U S A :10528-10532.13. Ciliberti S, Martin OC, Wagner A (2007) Robustness can evolve gradually incomplex regulatory gene networks with varying topology. PLoS Comp. Biology :e15.14. Kaneko K (2007) Evolution of robustness to noise and mutation in gene ex-pression dynamics. PLoS ONE :e434.15. Kaneko K and Furusawa C (2006) An evolutionary relationship between ge-netic variation and phenotypic fluctuation. J. Theo. Biol. :78-86.16. Kaneko K (2006)
Life: An Introduction to Complex Systems Biology (Springer, Heidelberg and New York).17. Sato K, Ito Y, Yomo T, Kaneko K (2003) On the relation between fluctuationand response in biological systems.
Proc. Nat. Acad. Sci. U S A :14086-14090.18. Lehner B. and Kaneko K., 2011, A macroscopic relationship between fluc-tuation and response in biology.
Cellular and Molecular Life Sciences , ,1005-101019. Fisher RA (1930) The Genetical Theory of Natural Selection , (Oxford Uni-versity Press).20. Futuyma DJ (1986)
Evolutionary Biology (Second edition), (Sinauer Asso-ciates Inc., Sunderland).21. Hartl DL and Clark AG (2007)
Principles of Population Genetics (SinauerAssoc. Inc., Sunderland) 4th ed.22. Price, G. R. (1970). Selection and covariance. Nature :520-521.23. Lande R., and Arnold S.J., (1983)
The Measurement of Selection on Corre-lated Characters
Evolution 37: 1210-122624. Ao P., (2005)
Laws in Darwinian evolutionary theory , Physics of Life Reviews, : 117-15625. Mustonen V., and Lassig M.,(2010) Fitness flux and ubiquity of adaptiveevolution, Proc. Nat. Acad. Sci. U S A , 107, 4248-425326. Elowitz MB, Levine AJ, Siggia ED, Swain PS (2002) Stochastic gene expres-sion in a single cell.
Science :1183-1187.27. Bar-Even A, et al. (2006) Noise in protein expression scales with naturalprotein abundance.
Nature Genetics :636-643.28. Kaern M, Elston TC, Blake WJ, Collins JJ (2005) Stochasticity in gene ex-pression: From theories to phenotypes. Nat. Rev. Genet. :451-464.29. Furusawa C, Suzuki T, Kashiwagi A, Yomo T, Kaneko K (2005) Ubiquity oflog-normal distributions in intra-cellular reaction dynamics. Biophysics :25-31.30. Tsuru S., Ichinose J., Sakurai T.,Kashiwagi A., Ying B-W., Kaneko K., andYomo T., ”Noisy cell growth rate leads to fluctuating protein concentrationin bacteria” Physical Biology 6 (2009) 036015331. Wakamoto Y. Ramsden J. and Yasuda K., (2005) Single-cell growth and di-vision dynamics showing epigenetic correlations Analyst, , 311-31732. Wright, S., (1932) The roles of mutation, inbreeding, crossbreeding and selec-tion in evolution , Proceedings of the sixth international congress on genetics , , 356–366,33. Sakata A, Hukushima K, Kaneko K (2009) Funnel landscape and mutationalrobustness as a result of evolution under thermal noise. Phys. Rev. Lett. :148101.34. Sakata A, Hukushima K, Kaneko K (2011), Replica symmetry breaking in anadiabatic spin-glass model of adaptive evolution, arXiv1111.5770v1, submit-ted to Europhys. Lett.35. Kaneko K (2009) Relationship among phenotypic plasticity, genetic and epi-genetic fluctuations, robustness, and evolvability.
J. BioSci . :529-542.36. Glass L and Kauffman SA (1973) The logical analysis of continuous, non-linearbiochemical control networks. J. Theor. Biol. :103-129.37. Mjolsness E, Sharp DH, Reisnitz J (1991) A connectionist model of develop-ment. J. Theor. Biol. :429-45338. Salazar-Ciudad I, Garcia-Fernandez J, and Sole, RV (2000) Gene networkscapable of pattern formation: from induction to reaction–diffusion.
J. Theor.Biol. : 587-603.39. Kaneko K. (2008) Shaping robust system through evolution.
Chaos :026112.40. Kaneko K. (2011) Proportionality between variances in gene expression in-duced by noise and mutation: consequence of evolutionary robustness. BMCEvol Biol :2741. Lafuerza, L.F. and Kaneko K., in preparation.42. Onuchic J.N., Wolynes P.G., Luthey-Schulten Z., and Socci N.D. (1995) “To-ward an Outline of the Topography of a Realistic Protein-Folding Funnel”, Proc. Nat. Acad. Sci. USA , 3626.43. Li, F., Long, T., Lu, Y., Ouyang, Q., and Tang, C. (2004) ”The yeast cell-cycle network is robustly designed’, Proc. Nat. Acad. Sci. U S A , , pp.10040–1004644. Abe H. and Go. N. (1980) ” Noninteracting local-structure model of foldingand unfolding transition in globular proteins. I. Formulation”, Biopolymers , 1013.45. Eigen M. and Schuster P. (1979) The Hypercycle , (Springer, Heidelberg).46. Landry CR, Lemos B, Rifkin, SA, Dickinson WJ, and Hartl DL (2007) GeneticProperties Influencing the Evolvability of Gene Expression.
Science : 118.47. Lehner B (2010) Genes Confer Similar Robustness to Environmental, Stochas-tic, and Genetic Perturbations in Yeast.
PLoS One e9035.48. Stearns SC, Kaiser M, Kawecki TJ (1995) The differential genetic and envi-ronmental canalization of fitness components in Drosophila melanogaster.
J.Evol. Biol. :539-557.49. Waddington CH (1957) The Strategy of the Genes (Allen & Unwin, London).50. Schmalhausen II (1949, reprinted 1986)
Factors of Evolution: The Theory ofStabilizing Selection (University of Chicago Press, Chicago).51. Slatkin, M. and Lande, R. 1976,
Niche width in a fluctuating environment-density independent model , American Naturalist, 31-55,52. Leibler S. and Kussel E., (2010)
Individual histories and selection in hetero-geneous populations Proc. Nat. Acad. Sci. U S A
The value of information for populationsin varying environments , Journal of Statistical Physics, , 1124–116654. Mustonen V., and Lassig M.,(2007) Adaptations to fluctuating selection in