Global Relationships in Fluctuation and Response in Adaptive Evolution

Cells generally change their internal state to adapt to an environmental change, and accordingly evolve in response to the new conditions. This process involves phenotypic changes that occur over several different time scales, ranging from faster environmental adaptation without a corresponding change in the genomic sequence to slower evolutionary dynamics involving genetic mutations and subsequent selection. In this regard, a question arises as to whether there are any relationships between such phenotypic changes over the different time scales at which adaptive evolution occurs. In this study, we analyzed simulated adaptive evolution in a simple cell model, and found that proportionality between concentration changes in adaptation and evolution over all components, and the proportion coefficients were closely linked to the change in the growth rate of a cell. Furthermore, we demonstrated that the phenotypic variances in component concentrations due to (non-genetic) noise and genomic alternations are proportional across all components. These global relationships in cellular states were also supported by phenomenological theory and transcriptome analysis of laboratory evolution in {\it Escherichia coli}. These findings provide a basis for the development of a quantitative theory of plasticity and robustness, and to determine the general restriction of phenotypic changes imposed by evolution.


I. INTRODUCTION
When an environmental condition is changed, biological systems change their state to adapt and evolve to the environmental change. Despite the recognized importance to characterize the capacity of adaptation and evolution, discussions on evolvability and plasticity have thus far remained at a qualitative, rather than quantitative, level. On the other hand, the cellular internal state can now be quantitatively determined by measuring the abundances of a variety of components, including proteins and metabolites. Recent advances in high-throughput experimental analysis enable quantification of changes within such a high-dimensional state space [1].
After an environmental change, cells may first respond by changing the abundances of cellular components without changing the genome sequence. The typical time scale of such environmental adaptation is generally shorter than several generations. On the other hand, over the long-term, i.e., over many generations, the internal state is gradually changed by evolutionary dynamics, in which the genome sequence is altered by mutations and individuals with higher fitness are generally selected. Indeed, experimental data of both changes in phenotype, reflecting changes in gene expression profiles, and changes in the genomic sequence throughout the course of evolution are now available; for example, much data are available from results obtained from the experimental evolution of E. coli [2][3][4].
Therefore, an important question arises: is there a general relationship between shortterm phenotypic changes in adaptation and long-term phenotypic changes in evolution? Of course, the phenotypic changes that occur over different time scales are generally caused by different mechanisms, and thus the existence of any relationship between them would be nontrivial. However, it should be noted that the essence of cellular dynamics is reproduction, in which the abundance of each cellular component is roughly doubled, and this constraint imposed by cellular reproduction imposes a restriction on the time development in highdimensional cellular state space. That is, it is possible that a (glivingh) cellular state would be restricted to a sub-space of the high-dimensional state space, described by a relatively smaller number of variables. Such restriction to low-dimensional dynamics can provide a non-trivial link between the phenotypic changes occurring in adaptation and the long-term changes occurring over the course of evolution [5]. In fact, some studies have suggested that there is a common trend over the thousands of gene expression changes in adaptation and evolution, in which genes whose expressions exhibit a larger response to environmental change tend to also show a larger response in their expression at the evolutionary scale [6,7]. Furthermore, global cellular behavior is represented by few macroscopic variables such as the growth rate and fitness of the system, which govern the entire (high-dimensional) dynamics in a cell. Therefore, it is important to uncover the possible relationships in the high-dimensional cellular dynamics between the expression of the thousands of proteins and metabolites and a macroscopic variable such as the growth rate, throughout the process of adaptation and evolution.
In addition to the phenotypic changes that occur after environmental changes during adaptation and evolution, the cellular state also generally exhibits fluctuations even under a constant environment and without genomic alternations, which originate from the stochastic nature of intra-cellular chemical reactions [8,9]. The possible relationship between such non-genetic fluctuations and adaptive responses has also garnered much attention recently.
The proportionality between such fluctuations and the evolutionary rate of fitness and phenotype has been demonstrated in bacterial experimental evolution and in simulations of toy cell models, which are supported by phenomenological theory [12][13][14][15], analogous to the proportionality between the fluctuation and response in statistical physics that has been well established since Einstein [10,11]. In the present context, evolution is considered the response to genetic change, and thus the proportionality between fluctuation and evolutionary rate means that components that are more variable by noise are also more variable by genetic change.
Considering the suggested proportionality between the response of cellular states to the environmental change and to genetic change, and also between the response and fluctuations, one may expect the existence of proportionality among two-by-two quantities, namely, fluctuations and responses induced by environmental (non-genetic) perturbations (noise) and by genetic changes (mutation). This grand relationship, if confirmed, is of critical importance to evolutionary biology, as it would provide a theoretical basis for the quantitative study on the plasticity and robustness underlying adaptive evolution, and could also set a general restriction as to the extent of phenotypic changes possible through (future) evolution. Indeed, the quantities constituting the relationship can now be measured over high-dimensional cellular dynamics, reflected as changes in the expression of thousands of genes. However, thus far, experimental confirmation of this grand relationship remains premature, and is in need of further scrutiny. Accordingly, at this stage, it is important to examine such a relationship by adopting an integrative approach combining in silico evolution of a cell model consisting of thousands of chemical species, laboratory evolution of bacteria under environmental stress, and phenomenological theory for time development in a high-dimensional state space.
In the present study, we aimed to uncover such statistical laws underlying the fluctuation and response of high-dimensional state variables occurring through adaptive evolution, and to connect them with changes in growth rate or fitness.

Evolutionary simulations of a simple cell model
We employed our previously established mutually catalytic reaction network model, as this model is capable of capturing the basic characteristic of cells such as the power-law abundances, log-normal fluctuations, fluctuation-response relationship of fitness, adaptation with fold-change detection, and so forth, in spite of its simplicity [9,13,16,17]. In the model, the cellular state is represented by a set of molecule numbers (N 1 , N 2 , · · · , N K ), where N i is the number of molecules of the chemical species i, which ranges from i = 1 to K. For the internal chemical reaction dynamics, we chose a catalytic network among these K chemical species, where each reaction from some chemical i to some other chemical j is catalyzed by a third chemical ℓ. Some resources (nutrients) are supplied from the environment by transportation through the cell membrane with the aid of some other chemicals that are termed 'transporters'. The environmental condition is given by the concentrations of nutrient chemicals. Through catalytic reactions, these nutrients are transformed into cell-component chemicals, and a cell divides when the amount of component chemicals reaches a certain threshold. Here, to achieve a higher growth rate, the synthesis of the cell components, transporters, and chemicals that catalyze the synthesis of those components need to be harmonized with the nutrient uptake. We allowed the above toy-cell model consisting of catalytic reaction networks to evolve by rewiring the network paths with a given mutation rate and selecting the pathways with a certain fraction of cells that showed a higher growth rate (See Methods for details). For a given environmental condition, evolution progresses so that the cell growth rate, i.e., the inverse of the average division time, is increased ( Fig. 1). To study the response to environmental change, we then switched the nutrient condition after evolution under a fixed condition for 3000 generations (denoted by the arrow in Fig. 1). The growth rate initially decreased following this environmental change, and then recovered through genetic evolution over generations. Next, we explored the phenotypic state changes in response to the environmental and evolutionary changes in order to evaluate the relationship between non-genetic and genetic responses upon environmental change.
As phenotypic state variables for the cell, we computed the abundances of each chemical N i at the division event. Here, it is convenient to choose X i = log N i as a phenotypic variable, since the abundance generally increases exponentially over time through cellular growth, and perturbation in a network is also generally amplified exponentially. Indeed, this choice of logarithmic abundances is also relevant to the theoretical argument presented below, as well as to transcriptome analysis of gene expression. Note also that the abundances N i are distributed by cells, even for those sharing the same reaction network, due to stochasticity in reaction dynamics; thus, the average abundance over all cells is required to study the mean response of cells, denoted by · · · .
After the change in nutrient condition, the abundances of all the components change.
Let us denote the average change of these abundances by: δX Env In other words, evolution shows a common tendency to abolish the changes in components introduced by the environmental change. This common proportionality across all chemicals suggests that the proportion coefficient r(m) is a "global variable" over a huge number of chemical species. A reasonable candidate for such a global variable is the cell growth rate µ. Hence, it is natural to compare the coefficient r(m) with the growth rate. Toward this end, we again computed the change in the growth rate δµ Env = µ(1) − µ(0)(< 0) and δµ Gen (m) = µ(m) − µ(0) at the m th generation. The ratio δµ Gen (m)/δµ Env gives an index for the recovery in this growth rate from the decrease caused by the environmental change, with 0 and 1 representing full and null recovery, respectively. In Fig. 3, the proportion coefficient r(m) is plotted against this growth rate recovery δµ Gen (m)/δµ Env . The proportionality between the two is clearly discernible.
Recalling the possible relationship between fluctuation and response, as is typical in statistical physics, we then evaluated whether there also exists a common relationship among the variances of all the components. Here, as previously reported [9], the distribution of each N i follows an approximately log-normal distribution, as confirmed experimentally in the protein abundances in the present cells. Hence, it is again relevant to adopt X i = log N i as a phenotype variable, so that the distribution of X i follows a roughly Gaussian distribution [12]. The phenotypic variance V ip (i) for each component i is defined as the variance of X i in an isogenic population. On the other hand, the variance due to genetic change V g (i) is defined as the variance of mean X i over a heterogenic distribution, where the mean is computed across clones of a given genotype (i.e., a network), while the heterogenic distribution is related to different genotypes (networks) that exist at a given generation.
The results of simulations are given in Fig. 4, which shows proportionality between V ip (i) and V g (i) across the components i for the evolved population. As the mutation rate is increased over evolution, V g (i) increases as the genotype distribution is broadened, whereas V ip (i) remains at the same level, so that the ratio V g (i)/V ip (i) is increased while maintaining the proportionality. This observed proportionality means that the components that are more variable owing to noise in the reaction dynamics are also more variable owing to mutation.
So far, we have confirmed the existence of common proportionality between non-genetic and genetic variances, as well as between the environmental and evolutionary responses. As the proportionality between the fluctuation and response is a natural outcome in statistical physics, we compared the response and fluctuations in more detail. However, direct comparison of the phenotypic variance V ip (i) with the environmental response did not show a clearly discernable proportionality. This is probably due to the discrepancy in the definitions of the two quantities: the variance originates from very high-dimensional dynamics without any specific directional change, while in the environmental response, only one specific environmental change considering only a few nutrients is applied. To make a more direct comparison, we then sampled the environmental responses against a variety of external changes introduced by different nutrient conditions to define the average environmental response R env (i) = (δX env ) 2 with · · · over 10 4 environmental conditions (see Methods).
As shown in Fig. S1, this average environmental response showed clear proportionality with V ip (i) and V g (i), respectively. Hence, the proportionality relationships among two-by-two quantities, i.e., genetic and non-genetic responses and fluctuations, hold over all components (see Fig. 5).
Thus, the degree of plasticity required to achieve an adaptive response to a new environment is characterized by fluctuations V ip (i), i.e., those that do not consider environmental or genetic changes. On the other hand, when cells are exposed to a novel environment, the potential of adaptation is expected to increase. Therefore, when placed in a novel condition, it is expected that the phenotypic fluctuations would increase to allow the cells to adapt to the new environment. In Fig. S2, the variances (V ip (i), V g (i)) are plotted before and after the environmental change (arrow in Fig. 1). In this case, all of the variances increase while roughly maintaining their proportionality. After this increase, the variances decrease over generations under a fixed environmental condition, while the proportionality between V ip (i) and V g (i) is maintained.

Theoretical Argument
By using a simple cell model, we have confirmed the common proportionality over thousands of components for genetic and non-genetic variances, and environmental and genetic responses. The results suggest the existence of a global variable that governs adaptive evolution. Here, the growth rate µ of a cell is a candidate for such a variable, since, for a cell to maintain its composition, every component has to be synthesized in conjunction with the growth rate. Indeed, in [5] we considered the dynamics of gene expression where µx i gives the dilution of the concentration by the increase in cell volume V , and x i is the concentration of the component i, , the original stationary state is given by Now, with the change in environmental condition E and genetic change G, the expression X i is shifted to X * i + δX i , and µ is shifted to µ + δµ. Assuming that the change in logarithmic concentration is not so large, and taking only the linear part of the changes and using the where γ E i ≡ ∂F i ∂E and γ G i ≡ ∂F i ∂G , respectively. Here, G is a coordinate introduced to represent the genetic change. It is not evident that the genetic change is represented by only a single variable. However, considering that under this scenario evolution progresses under a stressed environmental condition, one could project high-dimensional genetic change in the direction required to increase fitness (growth rate) under the condition, indicating that a single variable G can be introduced; indeed, several studies conducted to date support this assumption [12,13,15]. Accordingly, the variable G has the same dimensions as E, and can be scaled so that G and E induce the same degree of change in expression. The genetic evolution following the initial stress δE is expected to diminish the environmental stresses, so that evolution occurs in the direction δG < 0, if the environmental change δE is positive. Considering that evolution occurs through the projected direction in δE, it is natural to assume γ E i = γ G i (although this might be a crude approximation). Under the linear conditions of interest, the change in µ is proportional to δE or δG, with δµ(δE, δG) = α(δE + δG). Note that, again, the direction of δG is opposite to δE. Thus, we obtain, Then, over the course of evolution δG = 0 to δG(m), under a given environmental condition In other words, all the expression changes are proportional, as confirmed in the present simulations. To check the validity of the theory, we compared the proportion coefficient in the expression change (LHS of eq.5) with the change in growth rate (RHS) numerically through the course of the evolution simulation. As shown in Fig. 3, the relationship of eq. 5 holds rather well. Note that if there is deviation from γ E i = γ G i over i, the proportionality over all genes will deviate. In other words, the deviation from δX j (δE, δG) ∝ δX j (δE, 0) across genes i reflects the deviation between γ E j and γ G j . The relationship in the variances V ip (i) and V g (i) is considered in a similar manner.
Consider that the fluctuation in δE and δG induce fluctuation in each expression i, induced by either noise or genetic variation. This fluctuation induces variation in the growth rate, according to δµ = αδE or αδG, so that where δΥ is either δE or δG, i.e., phenotypic change induced by variation in the environment (i.e., noise) or by genetic change (e.g., by mutation), and · · · is the average over the distribution induced by the phenotypic noise or genetic variation. The variance V ip (j) and V g (j) are (δX j (δE)) 2 and (δX j (δG)) 2 , respectively, so that Thus, the ratio of the two variances takes on the same value independent of j, which is determined by the ratio of variances in growth rate fluctuations induced by noise to those induced by genetic variation. This relationship was again confirmed in our simulated evolution model (see Fig. S3) (see also [14] for an alternative derivation of the common proportionality between V ip (i)/V g (i) assuming the common error of catastrophe in the phenotype distribution).
The above theoretical interpretations on both the proportionality in responses and in variances suggest the importance of changes in the growth rate. Since the growth rate globally governs all of the concentrations through dilution, its dominance over each component is a reasonable assumption. To further evaluate the relationship between environmental and evolutionary dynamics, however, we need to also assume that evolution progresses as to assimilate environmental change, as Waddington proposed [18]. In our theory, this genetic assimilation is formulated by the introduction of the variable G that has a similar effect with the environment (or compensates for the environmental stress), so that

Experimental Verification
The  [19,20]. In this experiment, after cultivation of approximately 1,000 generations (2,500 hours) under 5% ethanol stress, 6 independent ethanol-tolerant strains were obtained, which exhibited an approximately 2-fold increase in specific growth rates in comparison to the ancestor. For all independent culture series, mRNA samples were extracted from approximately 10 8 cells at 6 different time points, and the absolute expression levels were quantified by using microarray analysis. All mRNA samples were obtained from the cells in exponential growth phase, which means that the changes in cellular state over the time scale of several generations were negligible, and each expression level represented cells in a steady-growth state (see [20] for details of materials and methods).
Using the time-series expression data of bacterial adaptive evolution, we analyzed the common proportionality in expression changes. The environmental response of the i-th gene δX Env i is defined by the log-transformed ratio of the expression level of the i-th gene obtained 24 hours after exposure to the stress condition to that obtained under the no-stress condition. Similarly, the evolutionary response at n hours after the exposure to the stress δX Gen i (n) is defined by the log-transformed ratio of the expression level at n hours to that of the non-stress condition. We found a common trend between the environmental and genetic responses over all genes, as shown in Fig. 6(a). Furthermore, we also found that the proportion coefficient r(n) for δX Gen i (n)/δX Env i is roughly proportional to the growth recovery ratio δµ Gen (n)/δµ Env , as shown in Fig. 6(b), where δµ Gen (n) and δµ Env are the growth rate differences of n hours and 24 hours after the exposure to stress, respectively. The results demonstrated that the evolutionary dynamics with growth recovery were accompanied by gene expression changes to eliminate the phenotypic changes introduced by the new environment, and agreed well with the simulation results of the simple cell model shown in Fig.   2 as well as the theoretical argument presented above.
Furthermore, there is some indirect experimental support for the proposed relationships between the variances. Stearns and colleagues measured the isogenic variance V ip of five lifehistory traits (such as body weight, lifespan, etc.) in Drosophila melanogaster, as well as the genetic variance V g between different genetic lines observed during laboratory evolution to increase the traits, and observed proportionality between the two [21]. The correlation between the isogenic variances in trait expression and variance due to mutation (but without selection) was measured across a few thousand genes in Saccharomyces cerevisiae. Correlation between the two variances was observed [22], whereas proportionality was not so clear. This is possibly because evolution without selection was applied in the experiment, and therefore only the variance resulting from random mutation was measured.

III. DISCUSSION
We have shown proportionality in the change in the concentrations of most intra-cellular components as a result of adaptive evolution, which was confirmed in simulated evolution of catalytic reaction network models of cells, laboratory experiments of bacterial evolution, and phenomenological theory. As the theoretical argument, albeit phenomenological, is rather general, we expect that the observed relationships obtained from the simulation and laboratory experiments represent a general phenomenon, independent of the specific models or organisms considered. This proportionality across thousands of components implies that there is a strong constraint in phenotypic evolution. In particular, the expression of different components cannot evolve independently, but rather change together, for the most part, along a one-dimensional path provided by eq.(5). Phenotypic change in adaptive evolution under a fixed environmental condition is highly constrained; thus, we here quantitatively formulate the general restriction or feasibility of the direction of phenotype changes in future evolution.
Our result also implies that the changes in the concentrations of most components that are induced by the environmental change become relaxed through the evolutionary process.
This suggests a phenomenon of strong homeostasis, that is, a tendency to restore to the original, adapted, intra-cellular states, via genetic change. In some sense, this homoeostasis is similar to the Le Chatelier principle in thermodynamics, in that changes introduced by external perturbations are relaxed by subsequent temporal evolution.
There could be a huge variety of genetic changes that yield the phenotypic changes required for adaptation. Indeed, in our simulations, there were a variety of network structures that could achieve phenotypic adaptation. When the simulation was run again with a different seed of random numbers for mutations, the resulting network (i.e., genotypes) was different in each run, but the change in concentrations (phenotypes) followed the proportionality given by eq.(5), independently of the specific genetic changes occurring during evolution.
Furthermore, in bacterial evolution experiments, the results from different strains tended to follow the same proportionality law described by eq.(5). It is interesting to note that such correlated change in expression levels by genetic changes is also suggested in several experiments [6,7]. It will be important to further confirm the relationship eq.(5) in more laboratory evolution experiments, and to also unveil the underlying genotype-phenotype map that achieves the common, restricted change in expression levels observed in the experimental data.
We have also found proportionality in the fluctuations in expression levels across components. As expected from Fisher's fundamental theorem of natural selection [23], the higher the genetic variance, the higher the evolutionary rate. Hence, the proportionality between V ip (i) ∝ V g (i) suggests that a higher isogenic variance of a given expression level due to noise would be accompanied by a higher rate in the change in the expression level due to evolution. Hence, our results suggest that the direction of evolutionary change in phenotypic space is likely to be predetermined by the isogenic variance of expression level due to noise.
According to our theoretical framework, the responses and fluctuations in expression levels are represented by the macroscopic growth rate and its fluctuation. Therefore, the relationship between the response and fluctuations, analogous to thermodynamics, is represented by the landscape of the growth rate as a function of phenotype (expression level) and the environment, in contrast to the established fitness landscape represented in genetic space proposed by Sewall Wright [24]. We hope that the present study will provide a basis for the development of a future macroscopic theory for phenotypic evolution.

IV. METHODS: MODEL SIMULATIONS
The cellular state can be represented by a set of numbers (N 1 , N 2 , · · · , N K ), where N i is the number of molecules of the chemical species i with i ranging from i = 1 to K. For the internal chemical reaction dynamics, we chose a catalytic network among these k chemical species, where each reaction from some chemical i to some other chemical j is assumed to be catalyzed by a third chemical ℓ, i.e., (i + ℓ → j + ℓ). A catalytic network is chosen randomly such that the probability that any two chemicals, i and j, are connected is given We studied the evolution of the replication dynamics by generating slightly modified networks and selecting those that grew faster. First, n parent cells were generated, and the connecting paths of catalytic networks were chosen randomly with connection rate ρ.
From each of the n parent cells, L mutant cells were generated by randomly replacing mρK 2 reaction paths, where ρK 2 is the total number of reactions and m is the mutation rate per reaction per generation. Then, reaction dynamics were simulated for each of the nL cells to determine the rate of growth of each cell; that is, the inverse of the time required for division. Within the cell population, n cells with faster growth rates were selected to be the parent cells of the next generation, from which nL mutant cells were again generated in the same manner. Throughout the simulation, the parameters were set as n = 1000 and L = 5.
In the simulations shown in Fig. 1, the mutation rate m was set to 1 × 10 −3 .
The environmental change is given by changing the nutrient concentration ratio in the en-vironment. In the evolutionary simulation shown in Fig. 1, there are two nutrient chemicals, each associated with one transporter chemical. Before adding the new environmental condition (generation≤0), the concentrations of these two nutrients in the environment (c 1 , c 2 ) were set to (0.5, 0.5), while after the environmental change (generation>0), they were set to (0.9, 0.1). In the result shown in Fig. S1, to add a variety of environmental changes, we randomly selected a nutrient chemical and a transporter chemical for this nutrient among K total chemical species. Then, the concentrations of the new nutrient c N ew and the original nutrients were set to (c 1 , c 2 , c N ew ) = (0.45, 0.45, 0.1). We iterated the random addition of a nutrient 10 4 times to obtain the average environmental response R env (i).
Research on Innovative Areas [25128715 and 26119719 to C.F.] from MEXT, Japan.
Author contribution: C.F. and K.K. designed the study. C.F. performed the simulations and data analysis. K.K. performed the theoretical analysis. C.F. and K.K. wrote the manuscript.
Competing financial interests: The authors declare no competing financial interests.        The relationship between V ip (i) and V g (i). The variances were computed by using the network and environment at generation 0 (before the environmental change) shown in Fig.   S1 with various mutation rates. V ip (i) and V g (i) were calculated based on the simulation results of randomly generated 10 5 networks. The solid line is y = x for reference.  for n = 384, 744, 1224, 1824, and 2496 hours. The growth recovery rate δµ Gen (n)/δµ Env was calculated based on the experimental measurements (see [20] for details). Among the 6 independent culture lines in [20], the results of 5 culture lines without genome duplication are plotted. The black line is y = x for reference.   and V g (µ)/V i p(µ). The variance ratio V g (j)/V ip (j) is calculated by the least square method for all components. The data points are obtained with m = 10 −6 × 2 ℓ for ℓ = 1, 2, · · · , 8.