Evolution enhances mutational robustness and suppresses the emergence of a new phenotype: A new computational approach for studying evolution.

The aim of this paper is two-fold. First, we propose a new computational method to investigate the particularities of evolution. Second, we apply this method to a model of gene regulatory networks (GRNs) and explore the evolution of mutational robustness and bistability. Living systems have developed their functions through evolutionary processes. To understand the particularities of this process theoretically, evolutionary simulation (ES) alone is insufficient because the outcomes of ES depend on evolutionary pathways. We need a reference system for comparison. An appropriate reference system for this purpose is an ensemble of the randomly sampled genotypes. However, generating high-fitness genotypes by simple random sampling is difficult because such genotypes are rare. In this study, we used the multicanonical Monte Carlo method developed in statistical physics to construct a reference ensemble of GRNs and compared it with the outcomes of ES. We obtained the following results. First, mutational robustness was significantly higher in ES than in the reference ensemble at the same fitness level. Second, the emergence of a new phenotype, bistability, was delayed in evolution. Third, the bistable group of GRNs contains many mutationally fragile GRNs compared with those in the non-bistable group. This suggests that the delayed emergence of bistability is a consequence of the mutation-selection mechanism.

Introduction

As living systems have developed through the long history of evolution, their present forms reflect their evolutionary histories. Thus, some properties of living systems should be consequences of the particularity inherent in evolution. In contrast, some properties may be more commonly observed, so that they do not depend on the evolutionary pathway. In this respect, the common properties and particularities of evolution are of specific interest. However, experimental studies on this topic are limited because the existing living organisms are products of evolution; thus, evolutionary experiments can only provide outcomes that reflect evolutionary histories. Therefore, computational methods are indispensable for obtaining information on evolutionary processes.

Although evolutionary simulation (ES) is a powerful method for studying the evolutionary process, its outcomes depend strongly on evolutionary pathways; therefore, ES alone is insufficient for our purposes and we need a reference system for comparison. In this study, we investigated the evolution of gene regulatory networks (GRNs). The reference system that we consider appropriate for this case is an ensemble of randomly sampled GRNs. If some properties are commonly observed in the ensemble, they should be realized irrespective of the evolutionary pathway. If there are some differences between the results of the ES and reference ensemble, they are manifestations of the particularity of evolution.

The aim of this paper is two-fold. First, we propose a research method for determining the properties of the evolutionary process by comparing a randomly sampled ensemble of genotypes obtained using a method based on statistical mechanics with genotypes obtained using the ES. Second, we apply it to a model of GRNs and investigate the evolution of mutational robustness and emergence of bistability.

Let us discuss our methodology concretely in the context of mutational robustness and bistability of GRNs. Living systems possess many types of robustness including that against environmental and internal noises as well as that against genomic mutations [1-3]. Among the various types of robustness, the most important is mutational robustness, which enables a living system to retain its function and continue to exist despite genome mutation.

Mutational robustness has been demonstrated experimentally. For example, comprehensive single-gene knockout experiments for bacteria and yeasts have revealed that most gene knockouts do not affect viability [4-6]. An artificial rewiring experiment for the GRN of E. coli showed that the bacterium remains viable after most artificial additions of regulatory links [7]. Such mutational robustness must have been acquired through an evolutionary process. Mutationally robust genotypes have been selected because of the mutation-selection mechanism, which is not directly related to any function. In this respect, enhancement of mutational robustness through evolution is called “second-order selection” [2].

Random sampling is suitable for identifying properties that do not depend on evolutionary history. Ciliberti et al. conducted random sampling of GRNs to investigate mutational robustness [8, 9]. The fitness of their model had only two values—viable and non-viable. They investigated the interrelations of viable GRNs and argued that most GRNs belong to a large cluster connected by neutral mutations, similar to the neutral networks found in the RNA sequence space [10]. However, such a simple random sampling method is not useful for systems with a more complex fitness landscape, because highly fit GRNs are rare. Burda et al. and Zagorski et al. employed the Markov chain Monte Carlo method to sample highly fit GRNs [11, 12]. They found that GRNs exhibiting multistability contained a common network motif.

The abovementioned methods are insufficient for sampling GRNs with a wide range of fitness levels. Saito and Kikuchi proposed the use of the multicanonical Monte Carlo (McMC) method to investigate the mutational robustness of GRNs [13]. McMC was originally developed in statistical physics for sampling configurations within a wide range of energies [14, 15]. However, this method has also been found to be useful for sampling nonphysical systems [16]. It is particularly effective for generating very rare states and estimating the probabilities of their appearance [16-20].

Nagata and Kikuchi investigated a GRN model using McMC [21]; they regarded fitness as the “energy” of GRNs and sampled GRNs with very low fitness to very high fitness uniformly and randomly. They considered a neural network-like model of GRNs, with one input gene and one output gene, and set the fitness in a manner that fitness was high if the GRN responded sensitively to changes in input. As each gene in their model did not respond ultrasensitively [22, 23], and the network structures were restricted, simple network motifs did not give rise to bistability [24]. Despite this, they found that highly fit GRNs always exhibited bistability. Therefore, bistability has emerged as a consequence of the cooperation of many genes. According to their results, a new phenotype of bistability appears regardless of the evolutionary path. They also found that mutationally robust GRNs were not rare among highly fit GRNs.

Bistable or multistable responses of GRNs are widely observed in living systems. The best-known example is the toggle switch for lysogenic-lytic transition in phage λ, which has been extensively studied both experimentally and theoretically [25-28]. Another example is the cdk1 activation system in Xenopus eggs [29-31]. The bistable switches of GRNs are also utilized in cell-fate decisions; a well-known example is the bistability of the MAPK cascade which regulates maturation of the Xenopus oocyte [32]. While the roles of small motifs have been the focus of research of such systems, the importance of cooperativity among many genes has been stressed theoretically [33].

In this study, we constructed a reference ensemble using McMC and compared it with the results of ES. We extended the research of Ref. [21] to more general network structures and explored both the enhancement of mutational robustness by evolution and how evolution affects the emergence of bistability.

Model

Genes encoded by DNA in cells are read by RNA polymerase and transcribed to mRNA, which is then used to assemble proteins. A category of proteins called transcription factors acts as activators or repressors of other genes. Many genes regulate each other in this manner and form a complex network called a GRN. GRNs are used to alter the cell state to adapt to environmental changes or to control the cell cycle or cell differentiation. One of the best-known GRN mathematical models is the Boolean network model proposed by Kauffman, in which each fixed point of the dynamical system is considered to represent a cell state [34-36].

In this study, we considered regulatory relations and ignored the details of gene expression. Such connectionist-type modeling has been widely used in theoretical studies [21, 33, 37–44]. We represented GRNs as directed graphs, with nodes as genes and edges as regulatory interactions. For simplicity, we considered GRNs with one input node and one output node. In contrast to Ref. [21], wherein several restrictions were imposed on the network structure, we allowed any network as long as the number of edges from one node to another was at most one. In the following sections, we have restricted our discussion to networks with N = 32 nodes and K = 80 edges.

A variable x ∈ [0, 1] that represents the expression level of a gene is assigned to each node, where i indicates the node number. x obeys the following discrete-time dynamics [38, 39]: where t denotes the time step. J represents regulation from the j-th node to the i-th node. For simplicity, we assume that J takes one of the three values −0, ±1; + 1 indicates activation, −1 indicates repression, and 0 indicates the absence of regulation. The 0th node is the input gene, and I ∈ [0, 1] is the strength of the input signal [44]. The response of the genes is given by the following sigmoidal function: where β represents the steepness of the function, and μ is the threshold. This function is widely used in theoretical studies [33, 40, 44–46]. In the present study, we set β = 2 and μ = 0. These parameters are the same as those in the neural network model used by Hopfield and Tank [47], and provide a gradual increase to the response function, which reflects the stochastic nature of gene expression [33]. The response function R(x) with μ = 0 is not ultrasensitive because a single gene with this response function has a single fixed point even when the auto-activation loop is attached. Therefore, for the emergence of multistability, many genes must act cooperatively. An example of the parameters of R(x) that a single gene with an auto-activation loop exhibits bistability is β = 6 and μ = 0.5. Although spontaneous expression R(0) = 0.5 is rather high, we do not believe that it caused any problems in the present investigation.

Following the definition of fitness presented by Ref. [21], we set the fitness such that it was larger when the difference in the expression levels of the output gene between I = 0 and 1 (“off” and “on”) was large. was the fixed-point value of the output node for fixed I; as the initial condition, we set the expression of all genes as 0.5, the spontaneous expression. If x(t) behaved oscillatorily over time instead of reaching the fixed point, we used the temporal average for ; however, the system reached a fixed point in most cases. Fitness f was defined as follows:

Fitness takes a value in [0, 1] by definition.

We sampled the GRNs using McMC to construct the reference ensemble. For this purpose, we divided the entire range of fitness into 100 bins. McMC enabled us to sample GRNs such that the numbers of GRNs in all bins were almost the same. The GRNs within each bin were randomly sampled in principle. Hereafter, we call this method “random sampling”. We performed two types of ES: in one, half of the GRNs in the population were preserved at each generation, and in the other type, 90% of the GRNs were preserved. We referred to these as Evo50 and Evo90, respectively. In the following sections, we mainly describe the results of Evo50, unless otherwise stated. For each ES, we prepared a random initial population comprising of a thousand GRNs. The details of the computational methods are summarized in the Methods section.

Results

Genotypic entropy and speed of evolution

The blue line in Fig 1a shows the base 10 logarithm of the appearance probability Ω(f) for each bin of fitness f obtained by random sampling. As the logarithm of probability is entropy, we refer it to as the “genotypic entropy.” The sum of Ω(f) was normalized to 1. The probability for f ≥ 0.99 was ∼10−16. As we could count the total possible number of GRNs as , highly fit GRNs were numerous but rare. The fitness dependence of genotypic entropy is divided roughly into three regions. The majority of GRNs concentrate near f ∼ 0. The number of GRNs then decreases exponentially with f, and, for a very high f, they decrease faster than the exponential rate. For comparison, we have shown the genotypic entropy for a steeper response function, β = 4 and μ = 0, in Fig 1b. It also comprises three regions, namely, the majority around f ∼ 0, an exponential decrease, and a faster-than-exponential decrease, as in the case of β = 2 (the jaggy outline is not because of statistical error).

Fig 1

Genotypic entropy and evolution of fitness.

The fitness f is divided into 100 bins. The blue dashed lines (left axis) show the base 10 logarithms of the appearance probability Ω(f) of each bin, obtained by random sampling. The orange solid lines (right axis) represent the average fitness of each generation calculated for the lineages obtained using Evo50. Averages were taken over 100,000 lineages. The vertical lines indicate the fitness at which Ω(f) starts to decrease faster than the exponential rate. (a) β = 2 and μ = 0. (b) β = 4 and μ = 0.

The orange lines represent the average fitness of each generation obtained with Evo50. The vertical axis represents the generation in the downward direction. Evolution progresses almost linearly in the early stage and slows down drastically for a large f. The values of f for which the increase in fitness starts to slow down roughly coincides with the values for which the faster-than-exponential decrease in Ω(f) begins for both β = 2 and 4. This may be because the number of possible destinations that a GRN can transit to by the mutation is restricted by Ω(f). In other words, when the number of GRNs with higher fitness levels decreases drastically, the possibility that the fitness increases by chance also decreases. A comparison of Evo50 and Evo90 for β = 2 is given in S1 Fig. Although the evolution was slower for Evo90, the overall tendency was similar. This result suggests that evolutionary speed depends in a large part on genotypic entropy.

Evolutionary enhancement of mutational robustness

To discuss mutational robustness, we introduce a measure of robustness. In Ref. [21], a single-edge deletion was considered as a mutation. It was found that the edges split into two classes, neutral and essential, for highly fit GRNs, and that the essential edges were minor among them; the essential edge means that the deletion of such an edge leads to fitness close to zero. We considered a single-edge deletion as a mutation in the present model as well, and obtained similar results. We then define the following quantity, r, as the robustness measure: where is the fitness after the ith edge is deleted, and the sum is taken for all edges. Because r should increase with f, comparing the r of GRNs with different f is not meaningful. However, by comparing this quantity for GRNs obtained by random sampling and ES at the same f, we can investigate how evolution affects mutational robustness.

Fig 2 shows the average of r against f for random sampling, Evo50, and Evo90. For the ES, we classified GRNs in the obtained lineages by fitness into 100 bins, as with random sampling. The average was taken over all the GRNs in the corresponding bin. For the highest fitness of the ES, only data in the range [0.990, 0.991) were used. 〈r〉 obtained by ES increased monotonically and coincided with those obtained by random sampling up to f ≃ 0.5. For a larger f, the value of evolution departed upward from that of random sampling, and the difference became increasingly significant as f increased. Evo50 and Evo90 behaved almost similarly, except for the highest fitness, where Evo90 exhibited a slight decrease. The reason for this decrease is not clear, but the standard errors are very small, and this decrease is not caused by a statistical error.

Fig 2

Average of the robustness measure r against fitness.

The average was taken over all the samples in each bin. The blue line represents random sampling. The orange and green lines represent Evo50 and Evo90, respectively. GRNs obtained by ES were classified according to f into 100 bins, similar to random sampling. The slight decrease at the highest fitness for Evo90 was not caused by a statistical error.

To scrutinize the difference, we show the probability distributions of r for f ∈ [0.5, 0.51), [0.8, 0.81), and [0.99, 1.0] in Fig 3a–3c. The data for the ES are taken from Evo50. Considering that the distribution of f within each bin differs for random sampling and ES, we divided each bin into ten sub-bins and reweighed the distribution obtained by ES, so that the distribution of f coincided with that of random sampling. While both distributions roughly agreed for f ∈ [0.5, 0.51), we observed a deviation for f ∈ [0.8, 0.81). The two distributions exhibited distinct differences for f ≥ 0.99, and the distribution obtained by ES was biased to a large r compared with that of random sampling. Therefore, evolution was found to enhance mutational robustness.

Fig 3

Probability distributions of the robustness measure r.

(a)f ∈ [0.5, 0.51) (b)f ∈ [0.8, 0.81) (c)f ∈ [0.99, 1.0]. The orange solid lines represent Evo50 and the blue dashed lines represent random sampling. r was divided into 80 bins because otherwise, unnecessary oscillation appeared in the distribution, owing to the discreteness of the number of essential edges. The distributions for the ES were corrected using the reweighting method described in the text.

What caused this difference in distribution? As stated previously, the edges of randomly sampled GRNs for f ≥ 0.99 are split into two classes: neutral and essential. Thus, when we delete one edge, f′ is either close to f or almost zero, and other types of edges are scarce. Therefore, we investigated the number distribution of essential edges for f ≥ 0.99 by stating that an edge is essential if f′ < 0.8. Although this definition is arbitrary, it hardly affects the results. Fig 4 shows the probability distribution of the number of essential edges n. GRNs obtained by ES are significantly biased toward a small number of sides compared to random sampling. The highest probability for ES was at n = 3, and the distribution was narrow; further, more than 2% of GRNs had no essential edge. In contrast, the highest probability for random sampling was at n = 15, and the distribution was much broader. The number of GRNs lacking an essential edge was only 0.4%. Therefore, the small number of essential edges is the cause of enhanced mutational robustness by evolution.

Fig 4

Probability distribution of the number of essential edges n.

The data for f ∈ [0.99, 1.0] are shown. An edge is regarded as essential if the fitness f′ becomes less than 0.8 after the edge is deleted. The orange solid line represents Evo50, and the blue dashed line represents random sampling. The distribution of the ES was corrected using the reweighting method described in the text.

One possible explanation for this is that fewer nodes affect the output in the evolutionarily obtained GRNs. To check this, we counted the number of nodes n that had at least one path to the output node. Fig 5 shows the distributions, and the distribution of random networks is also plotted for comparison. The peak of the distribution was at n = 28, which is slightly less than that for the peak of random networks; however, the distributions were indistinguishable between ES and random sampling. Therefore, the number of effective nodes does not cause a difference in the number of essential edges.

Fig 5

Distribution of the number of effective nodes n.

The data for f ∈ [0.99, 1.0] are shown. If a node has at least one path to the output node, the node is regarded as “effective.” The orange solid line represents Evo50, the blue dashed line represents random sampling, and the grey dotted line represents random networks as a reference. The distribution for the ES was corrected using the reweighting method described in the text.

Delayed emergence of bistability

The model in Ref. [21] exhibits bistability as f approaches its maximum value. In other words, when I is changed continuously, two saddle-node bifurcations occur, in contrast to the case of a small f, wherein a single fixed point moves to follow the change in I. We call the latter type of GRNs monostable. We found that the present model behaves similarly despite the fact that the network structures are significantly different owing to the different restrictions imposed on network construction.

Bistable GRNs are classified into three categories. The first is the toggle switch, in which two saddle-node bifurcations occur within I ∈ [0, 1]. The toggle switch is found, for instance, in phage λ and utilized for adaptation to environmental change [25-27]. The second is the one-way switch [48]. In this case, only one saddle-node bifurcation point is found in the range I ∈ [0, 1], and another bifurcation point is present outside this range. One-way switches give rise to cell maturation or cell differentiation; a typical example of such switches is the MAPK cascade in the maturation of Xenopus oocytes [32]. In the last category, both bifurcation points are outside the range of I. As this type may not work as a switch as long as the effect of noise is not considered, we called it “unswitchable.” In this study, we did not distinguish among these three and treated them equally as bistable GRNs. It is straightforward to change the definition of bistability to deal only with, for example, toggle switches.

First, we investigated bistability using a strict criterion. Bistability was checked as follows: Starting from the steady state at I = 0, I was increased by 0.001, and the dynamics were run until the steady state was reached. This procedure was repeated for up to I = 1. Then, the inverse process, from I = 1 to 0, was performed. If a difference in larger than 10−6 was observed between these two processes in a range of I, the GRN was considered to be bistable. We employed such a strict criterion because the monostable GRNs and bistable GRNs were mathematically different as dynamical systems, in that, the number of fixed points was different. Thus, classifying them as strictly as possible is meaningful. In this respect, we may regard monostability and bistability as different phenotypes.

Fig 6a shows the fraction of bistable GRNs P2(f) against fitness. The blue line is the result of random sampling. The bistable GRNs began to appear at f ≃ 0.5 and increased rapidly until all GRNs became bistable for f → 1. Therefore, a new phenotype of bistability emerged as fitness increased, and all GRNs converged to such a phenotype as fitness approached its maximum value. The orange and green lines are the results of Evo50 and Evo90, respectively. For all data, the standard errors were smaller than the mark. Fig 6a shows that P2(f) of the ES was substantially lower compared to that of random sampling at the same f. In other words, evolution delays the rapid increase in P2(f). Nonetheless, the eventual emergence of a bistable phenotype is inevitable as fitness increases, because all GRNs are bistable for f → 1. Evo50 and Evo90 behave similarly except for the early stage of increase; Evo90 initially coincides with random sampling and soon becomes confluent with Evo50. This indicates that the emergence of bistability in evolution depends, partly, on the evolutionary speed.

Fig 6

Fitness dependence of the fraction of bistable GRNs.

The blue line represents random sampling. The orange and green lines represent Evo50 and Evo90, respectively. The standard errors are smaller than the mark. Bistability was checked as follows. First, starting from the steady state at I = 0, I increased by ΔI, and the dynamics were run until a steady state was reached. This procedure was repeated for up to I = 1. Then, the inverse process, from I = 1 to 0, was performed. If the difference in larger than the threshold x was observed between these two processes in a range of I, the GRN was regarded as bistable. GRNs obtained by ES were classified according to f into 100 bins, similar to random sampling. (a) Strict criterion, ΔI = 0.001 and x = 10−6. (b) Loose criterion, ΔI = 0.01 and x = 0.5.

Although the strict criterion of bistability employed above is mathematically meaningful, very weak bistability may be biologically irrelevant, because it may not be distinguishable from monostability in living systems. Thus, we also investigated bistability using a looser criterion; the checking interval of I was set at 0.01, and the criterion of bistability was set as a difference in larger than 0.5. The results are shown in Fig 6b. While the rapid increase in the bistable GRNs moves to a higher f compared to that in Fig 6a, the tendency that the emergence of bistability is delayed in evolution compared to that of random sampling remains unchanged. Therefore, we consider that the delayed emergence of bistability is a biologically relevant phenomenon.

Relation between bistability and mutational robustness

So far, we observed the enhancement of mutational robustness and the delayed emergence of bistability by evolution. Thus, we expected that there would be a relationship between them. Fig 7a and 7b show the number distributions of fitness f and robustness r for the GRNs obtained by random sampling for the monostable and bistable GRNs, respectively. Two distributions show pronounced difference. For monostable GRNs, r increases with f and is distributed within a narrow range. In contrast, bistable GRNs include those with very low robustness, irrespective of their fitness values. As a result, the robustness of bistable GRNs is widespread. This difference suggests that if relatively robust GRNs for mutation are favored by evolution, bistable GRNs tend to be avoided.

Fig 7

Number distribution of fitness f and robustness measure r obtained via random sampling.

(a) the monostable GRNs (b) the bistable GRNs. Both f and r are divided into bins of width 0.01. Color of each bin indicates the base 10 logarithm of the number of GRNs. No GRN was found in white bins. The strict criterion was used to determine bistability. Note that the total number of GRNs in each bin of f is about 50000 irrespective of the value of f.

Motif analysis

Next, we investigated the network motifs. In the model presented in Ref. [21], coherent feedforward loops (FFL+), and positive feedback loops (FBL+) were significantly abundant in highly fit GRNs. As the present model allows both auto- and mutual-regulations, we also explored patterns other than triangular patterns. We counted the number of auto-regulations, mutual regulations between node pairs, triangles, and mutual activation or mutual repression of two nodes accompanied by auto-activation of both nodes. Because the motifs are defined as network patterns that are abundant compared to random networks [24, 49], we also counted them for random networks.

As a result, the following patterns were greater in number than those in the random networks: auto-activation, mutual activation, mutual repression, FFL+, FBL+, mutual activation accompanied by auto-activations of both nodes, and mutual repression accompanied by auto-activation of both nodes. Although their abundances were not remarkable, we called them motifs. The number distribution of auto-activation is shown in S2a Fig, and those of the other motifs are shown in S3 Fig. Other patterns, such as auto-repression, incoherent feedforward loop, and negative feedback loop, were fewer in number compared to those in random networks. In S2 Fig, we compare the number distributions of auto-activation and auto-repression; the former is a motif, whereas the latter is not. Although the number of auto-regulations is low, there are very few GRNs that do not undergo auto-activation. In contrast, almost half of the GRNs do not exhibit auto-repression. Thus, auto-activation is favored over auto-repression, but is not indispensable.

Overall, the distributions of these motifs were almost the same for both ES and random sampling. Therefore, whether or not these highly fit GRNs are products of evolution is not reflected in the distributions of these local motifs.

Path distribution

As a characteristic of the global structure, we counted the number of paths n connecting the input and output nodes. Fig 8a shows the distribution of paths starting from the input node and reaching the output node without passing the same node more than once. The data for f ∈ [0.99, 1.0] are shown for random sampling whereas the data for f ∈ [0.99, 0.991) are shown for ES. Two distributions exhibited a distinct difference; while the probability of GRNs with only one path reached 4% for random sampling, it was 0.1% for ES. In addition, evolutionarily obtained GRNs with less than approximately 100 paths were fewer than those obtained through random sampling. This suggests that the global structures of GRNs obtained using these two methods differ significantly. However, the difference in path distribution does not explain everything. Fig 8b shows a scatter plot of the number of paths and the number of essential edges. The figure shows that the number of essential edges is lower for ES, irrespective of the number of paths. Even for n = 1, the number of essential edges is distributed broadly. This means that the locations of the essential edges were not limited to the “on-path” locations between the input and output nodes.

Fig 8

Number of paths n starting from the input node and reaching the output node without passing the same node twice.

Orange indicates ES, and blue indicates random sampling. (a) Probability distribution of n. The data for n ≤ 400 are shown. The probability of GRNs with only one path reached 4% for random sampling. (b) Scatter plot of n and the number of essential edges n. The data for n ≤ 1000 are shown. The data for f ∈ [0.99, 1.0] are shown for random sampling whereas the data for f ∈ [0.99, 0.991) are shown for ES.

Steady-state evolution of the mutational robustness

Evolutionary simulations would eventually reach a steady state, and the robustness distribution in the steady state was expected to differ from that of random sampling, considering the results shown in Fig 3. To investigate the steady state, we conducted very long simulations of Evo90, and found that the fitness of all preserved GRNs exceeded 0.9999999 at the two-millionth generation. Such extremely high fitness values should be an artifact of the deterministic nature of GRN dynamics and should not be considered biologically relevant. We thus expect that noise is important for high-fitness GRNs.

Instead of introducing noise, we imposed the upper limit f ≤ 0.99 on fitness and conducted simulations of Evo90. In the simulations, fitness values greater than 0.99, were regarded as 0.99. Fig 9a shows the distributions of the number of essential edges n for GRNs of f = 0.99 at the generation where the fitness of all preserved GRNs first reached 0.99. The results of ten independent runs are shown. The number distributions differed from run to run, and while some populations exhibited a small number of essential edges, the overall tendency was that n was distributed broadly. As all preserved GRNs have the same f, the fitness-driven evolution should have ceased at the generations shown in the figure, after which neutral evolution would continue.

Fig 9

Distribution of the number of essential edges n for GRNs with f = 0.99.

We introduced the upper limit of 0.99 for fitness, and conducted Evo90. In other words, a fitness larger than 0.99 was regarded as f = 0.99. The data for all GRNs with f = 0.99 for each population were used. (a) Distribution of n at the generation where the fitness of all the preserved GRNs first reached 0.99. The results of ten independent runs are shown. (b) Average distribution of n in the steady state. After 2000 generations, we collected GRNs at every 2000 generations up to one million generations and took their average. The results of six independent runs are shown.

We found that the distribution became steady after approximately 2000 generations. We then conducted runs of one million generations and collected all the GRNs at every 2000 generations. Fig 9b shows the average distributions of n for GRNs with f = 0.99. The results of six independent runs are shown. Only two distinct distributions were observed; six runs were classified into four and two runs. We considered populations with the same essential edge distribution to be genetically similar. These distributions are biased to low n compared to Fig 9a, and the peaks of n are 2 and 4 for the two distributions. Moreover, the ratio of GRNs without an essential edge reached approximately 4% of each population. These results indicate that after the fitness distribution reached the maximum, the selection driven by mutational robustness progressed until a steady state was reached. As a result, only a limited number of genotypic groups remained in the steady state; in these six runs, we observed that they converged to only two distinct groups.

Summary and discussions

In this paper, we proposed a new computational method for studying the properties of evolutionary processes. By generating a reference ensemble via random sampling using the multicanonical Monte Carlo method and comparing it with the outcomes of evolutionary simulations, we can quantitatively explore the commonly observed properties and particularities of evolution. This method is both powerful and general. Using this method, we investigated the evolution of a gene regulatory network model focusing on mutational robustness and bistability.

We found that mutational robustness was markedly enhanced as fitness increased, compared to that in the randomly sampled ensemble with the same fitness, even though the selection is imposed only on fitness. The mechanism for this enhancement can be explained by “second-order selection” [2]. The mutation we considered comprises two successive procedures. First, a randomly selected edge is deleted. Next, a new edge is added between a randomly selected node pair. For fitness values higher than some intermediate values, the edges start to be divided roughly into two types, as observed in Ref. [21]: almost neutral ones, and those causing substantially decreased fitness when deleted. Here, we name the latter type as “essential.” If the deleted edge is essential the possibility of fitness to recover via the random addition of a new edge is very low. Therefore, the more essential edges a GRN has, the harder it is for its copy to survive.

The condition for mutational robustness to evolve has been discussed theoretically in the case of neutral evolution, based on population dynamics. According to this theory, the product of population size P, the number of edges K, and mutation rate μ should be sufficiently large [2, 50]. As μ is comparatively large in our ES, we consider that this condition is satisfied. We note a difference: the evolution was not neutral in our ES. The present results showed that mutational robustness was enhanced with increasing fitness in evolution. As mutations that increase fitness are considered rare, fitness increases intermittently. In contrast, mutational robustness can evolve even when mutations are almost neutral because deleterious mutations are expelled by selection.

Based on a GRN model with two-valued fitness, Ciliberti et al. reported that mutational robustness is enhanced significantly in GRNs that experienced natural selection compared to that in GRNs randomly selected from viable ones [8]. Several different selection pressures have also been reported to result in to mutationally robust GRNs [51]. Our results were consistent with these findings.

Next, we discuss bistability. Random sampling revealed that almost all (probably all) GRNs become bistable as fitness approaches its maximum value. Thus, bistability is an inevitable consequence of increased fitness, irrespective of the evolutionary pathway. We confirmed this using an ES. As bistability is not explicitly considered in fitness, it is an emergent property. We may thus regard bistability as a “new phenotype” and the present results indicate that this new phenotype would always appear even if evolution was rewound and restarted. Nagata and Kikuchi obtained the same results for their GRN model [21]. In contrast to our model, auto-regulation and mutual regulation were prohibited in their model. As a result, the distributions of network motifs were different for the two models. Nevertheless, both models exhibited similar bistabilities. Thus, bistability is a common property, irrespective of the network structure. This observation suggests that the possible phenotype is constrained by the form of fitness function.

By comparing the outcomes of ES with random sampling, we found that the appearance of bistability was significantly delayed in evolution. Random sampling revealed that the bistable group of GRNs contained many mutationally fragile GRNs compared to those in the monostable group of GRNs. In other words, bistable GRNs and monostable GRNs behave differently in terms of mutational robustness. This is a nontrivial finding that our methodology made possible. This result suggests that the delay in the emergence of bistability may arise from the tendency that mutationally robust GRNs are selected by evolution. Whether this scenario applies to other phenotypes when different fitness functions are considered is of interest for future studies.

A set of GRNs with high fitness composes the neutral space. Thus, we collected members of the neutral space using our method of random sampling by McMC. Ciliberti et al. analyzed the structure of the neutral space for the above-mentioned model in which simple random sampling was valid and found that high-fitness GRNs belong to a large neutral space [8]. Unfortunately, a similar analysis was difficult for GRNs obtained using McMC. Instead, we set the maximum value f = 0.99 and regarded all fitness values greater than this as 0.99 for conducting long ES. After all the preserved GRNs reached this maximum value, neutral evolution continued and eventually, a steady state was reached. By investigating the steady state, we found that the neutral space is divided into only a small number of parts. This is consistent with the results reported in Ref. [8]. The situation is similar to that studied in the population dynamics mentioned above and is consistent with it, given PKμ ≫ 1 [50]. A detailed analysis of the neutral space is required in future research.

Next, we discuss network structures. Characteristic motifs were found for the highly fit GRNs. Among them, coherent feedforward loops are frequently observed in actual GRNs [24, 49, 52]. The following structures are also identified as motifs, which are known to exhibit bistability if the response of each gene is ultrasensitive [23, 32, 53–58]: auto-activation, mutual-activation, mutual-repression, and mutual-activation/repression accompanied by auto-activation of both genes. The last motif is widely observed in multistable GRNs [11, 25, 26, 59–61]. However, these are not relevant to mutational robustness. This suggests that mutational robustness is related to the global structure of GRNs. We found that the number of paths connecting the input and output nodes differed between randomly sampled GRNs and evolutionarily obtained ones. Although it is understandable that GRNs with many such paths are robust because of redundancy, this does not fully explain the origin of mutational robustness.

Finally, evolutionary speed was roughly determined by the number of available GRNs or, in other words, “genotypic entropy.” This is partly because of our definition of a single-valued fitness function, which can be computed from the dynamics of a given GRN for a predetermined single task and has the maximum value. This setup is somewhat artificial, and fitness is not a simple function in reality. However, some experimental studies have addressed the situation discussed in this study. For example, Sato et al. reported a similar evolutionary study for a protein [62]. They found that evolutionary speed decreased as fitness increased. We consider that the variation in evolutionary speed in their experiment is explained almost fully by the entropic effect. In a natural situation departing from experimental conditions, evolution is not restricted to proceed in only one direction. When evolution in one direction becomes difficult, the direction changes. We consider that the concept of genotypic entropy affecting evolutionary speed can also apply to such a natural situation.

In summary, we have shown using a new computational method that mutational robustness is enhanced during evolution. We have pointed out the possibility that it affects the emergence of a new phenotype of bistability in GRNs. The research method proposed in this paper can be applied widely to evolution-related phenomena, not restricted to the evolution of GRNs. Further, the basic idea of generating a reference ensemble using McMC can be extended to other fields, such as the learning process of machine learning.

Methods

Random sampling

Random sampling was realized using the McMC method (more precisely, entropic sampling [63], which is one of the variations of McMC). The details of this method are described in Ref. [21]. We divided fitness into 100 bins and determined the weight for each bin such that the GRNs appeared evenly, using the Wang-Landau method [64, 65]. One McMC update comprised the following two successive processes: deleting a randomly selected edge and adding a new edge to an unlinked node pair that was also chosen randomly. Whether to accept this change was determined using the Metropolis method. One Monte Carlo step (MCS) comprised K such updates, and we recorded f at every MCS. We sampled GRNs at every 20 MCSs to reduce the correlation between samples. We conducted ten independent runs, each consisting of 107 MCSs. Therefore, we obtained an average of 50,000 samples for each bin. Although inter-sample correlations should have remained to some extent, we named this method “random sampling” in this study. The ensemble of these randomly sampled GRNs was considered as the reference ensemble.

Evolutionary simulation

We conducted two types of evolutionary simulations: Evo50 and Evo90. For both simulations, we prepared an initial population comprising 1000 randomly generated GRNs. The population size remained unchanged during the simulation. In each generation of Evo50, we selected the top 500 GRNs according to the fitness level and made one copy for each. In the case of Evo90, 90% of the population, namely, 900 GRNs, from the highest fitness were selected for preservation, and the remaining 10% of GRNs were discarded in each generation. We randomly selected 100 GRNs from the 900 preserved GRNs and made one copy for each. Then, these copies were subjected to mutation for both simulations. The mutation procedure for each GRN was the same as that in the McMC procedure; an edge was deleted randomly from the network, and a new edge was then added between a randomly selected unlinked node pair. As these procedures were repeated, the average fitness of the population increased. After 150 generations for Evo50 and 200 generations for Evo90, we selected the GRN with the highest fitness in the population and traced its ancestors to obtain a single lineage. We repeated this evolutionary simulation 100,000 times and 55,000 times independently for Evo50 and Evo90, respectively, and, as a result, we collected 100,000 lineages and 55,000 lineages, respectively.

Evolution of fitness for Evo50 and Evo90.

The orange and green lines represent the average fitness of each generation, calculated for the lineages obtained by evolutionary simulations, Evo50 and Evo90, respectively. Averages were taken over 100,000 and 55,000 lineages, respectively. The vertical line indicates the fitness at which Ω(f) starts to decrease faster than the exponential rate.

(PDF)

Click here for additional data file.

Probability distribution of the number of auto-regulations.

(a) Number of auto-activations n. (b) Number of auto-repressions n. The data for f ∈ [0.99, 1.0] are shown. The orange solid line and blue dashed line indicate ES and random sampling, respectively. The black dotted line represents the random networks as a reference. The distribution for the ES was corrected using the reweighting method described in the text. Auto-activation is a motif, whereas auto-repression is not a motif.

(PDF)

Click here for additional data file.

Probability distributions of the number of motifs.

The data for f ∈ [0.99, 1.0] are shown. The orange solid lines represent Evo50, and the blue dashed lines represent random sampling. The distributions for the random networks are also shown as references with gray dotted lines. The distributions for the ES were corrected using the reweighting method described in the text. (a) Mutual activation (b) Mutual repression (c) Coherent feed-forward loop (d) Positive feedback loop (e) Mutual activation accompanied by auto-activation of both genes (f) Mutual repression accompanied by auto-activation of both genes.

(PDF)

Click here for additional data file.

30 Mar 2021

Dear Dr. Kikuchi,

Thank you very much for submitting your manuscript "Evolution enhances mutational robustness and suppresses the emergence of a new phenotype" for consideration at PLOS Computational Biology.

As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by several independent reviewers. In light of the reviews (below this email), we would like to invite the resubmission of a significantly-revised version that takes into account the reviewers' comments.

We cannot make any decision about publication until we have seen the revised manuscript and your response to the reviewers' comments. Your revised manuscript is also likely to be sent to reviewers for further evaluation.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to the review comments and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Please prepare and submit your revised manuscript within 60 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email. Please note that revised manuscripts received after the 60-day due date may require evaluation and peer review similar to newly submitted manuscripts.

Thank you again for your submission. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Alexandre V. Morozov, Ph.D.

Associate Editor

PLOS Computational Biology

Sushmita Roy

Deputy Editor

PLOS Computational Biology

***********************

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: General overview

The authors consider in silico gene regulatory networks (GRNs) based on discrete time parallel updates of continuous (transcriptional) expression levels along with (i) an associated fitness function and (ii) evolutionary (mutational) dynamics. In that framework they compare properties of GRNs produced within two ensembles: the "uniform" ensemble of all GRNs of a given fitness (referred to random sampling ensemble) and the GRNs resulting from evolutionary dynamics given an initial starting point. The work's main conclusions are as follows. 1: At given fitness, mutational robustness is higher in the second ensemble; 2: the GRNs of high fitness typically exhibit bistability (and thus hysteresis), but less so in the second ensemble, which the authors interpret as a "delayed phenotype" effect; 3: evolution slows down as one reaches high fitness values.

Main recommendations of the referee

1: The author's use of a "random sampling" is both technically a major accomplishment and a useful tool for obtaining insights into what aspects are due to evolutionary dynamics and not just fitness. The authors use a sampling method (a Wang-Landau multicanonical MC) that produces a quantitative characterisation of the "fitness landscape" which allows the authors to refer to what they call "genotypic entropy". That is very good work and I think some of your readers would benefit further if you put a bit more accent on this point. The authors could have used simulations in the evolutionary ensemble to obtain a characterisation of the same type, leading to genotypic entropy in the evolutionary ensemble too. In physics terminology, you want to answer the question "What is the distribution of fitness in that ensemble", obtainable by direct simulation in the steady state. Theoreticians would be interested to compare that distribution to the Boltzmann one (adjusting the temperature to get the peak at the same fitness position), as far as I know one has little theoretical control over that distribution. Note that this requested computation uses the steady state distribution in the evolutionary dynamics, distribution that depends on the way mutations are produced and accepted from one generation to the next. A "noise level" (related to an effective temperature) is simply the fraction of top fitness individuals produced that are passed on to the next generation in your evolutionary simulation; if you accepted more than 500 individuals, the effective temperature would be higher. I recommend you implement this suggestion so that you no longer use non-equilibrium effects produced by transients and have thus a well defined (steady state) ensemble framework. Changing to this way of computing in your paper would be little work (in fact the authors perhaps already have the data for the noise level they used) but such a change would significantly enrichen the theoretical and conceptual messages of the paper. Furthermore I am quite convinced that the 3 main results put forward by the authors would hold in the steady state of the evolutionary dynamics ensemble. Lastly, note that an advantage of using steady states vs transients is that there then is no dependence on generation number. So to conclude this first point, please do the steady state measurements rather than transient ones and implement the "temperature" parameter in the evolutionary dynamics so that you can scan the average fitness, thereby allowing a much more systematic comparison of the two ensembles.

2: The first of the 3 conclusions given by the authors in their abstract, namely that the second ensemble enhances mutational robustness, is completely expected. Indeed, previous studies have shown that the ensemble based on evolutionary dynamics leads to GRNs with greater mutational robustness than random sampling. That property is very simple to justify qualitatively but it also has a deeper mathematical proof in the case of "0-1" fitness landscapes as given in the paper "Neutral evolution of mutational robustness" (PNAS 1999). So your work can be considered to be a generalization of the results from "0-1" fitness landscapes to arbitrary ones. It is a useful contribution but it is also completely expected. The following analogy from quantum mechanics may help you qualitatively justify the result in your case of continuous fitness landscapes: if one considers a square well potential in 1 dimension, the ground state will have low density at the well's border, corresponding to higher mutational robustness in that ground state compared to the mutational robustness in the trial wave function set to a constant in the well.

3: The second of the 3 conclusions given by the authors in their abstract, namely that bistability is delayed in the evolutionary ensemble, is an observation based on the authors' simulations but remains unexplained. Since the authors established it in their "non equilibrium" framework, we don't know if it is due to the transient dynamics or is a feature of the (steady state) evolutionary ensemble. This gives another reason to follow the request I made of using the steady state ensemble in point 1. Probably the result of bistability holds in the steady state ensemble too. Then a possible interpretation and justification for your result would arise if the "bistability" phenotype was anticorrelated with the mutational robustness. (I have to admit that this seems a bit counter intuitive to me but in fact I have very little basis for having any intuition for this subtle point.) Could you provide a scatter plot and test for such a negative correlation? If the correlation did come out to be negative, you would have a justification for the result, thus a more palatable message. Whether my intuition is right or wrong, I think the paper would be much improved if you could provide new material to put the claim on more solid ground and include a justification. Without a good result there via a revision, I fear your work will have little impact.

4: The third of the 3 conclusions given by the authors in their abstract, namely that evolution slows down as one reaches high fitness, is in a sense inevitable or even a tautology: if fitness cannot improve its rate of change goes to zero... Nevertheless it is just a qualitative observation whereas in a research paper one expects quantitative measures. Furthermore, the current approach of the authors which is to consider the dependence of fitness on generation time has 2 major drawbacks: it is based on the outcome of non-equilibrium dynamics and it is logically inevitable since as one approaches the maximum fitness it can no longer increase. If the authors want to keep this point, I recommend they use the same definition of "slowness" as in Monte Carlo via an autocorrelation time. In the steady state, you can measure that quantity. I would then expect it to grow (slower dynamics) as the selection pressure in the population is increased (or as the mean fitness in the population increases). Allowing such an increase requires that the evolutionary dynamics include a "temperature" parameter, for instance the fraction of kept children over the population size.

Minor comments

1: In the introduction, the authors write "Thus, the central question we pose in this study is whether mutational robustness really evolves". Since the stationary probabilities in your reference ensemble and that produced by evolutionary dynamics inevitably are different, there will necessarily be evolution of mutational robustness in your study. I don't think it is good to say your paper is mainly about this point.

2: You write "this response function has a single fixed point even when the auto-activation loop is attached". Please mention that this feature depends on the values of beta and mu (and possibly at which values auto-activation allows for multistability at the single gene level).

3: Since the Jij space is discrete and your response function R has a maximum that is strictly less than 1, there has to be a drop in GRN entropy at f approaching that limit which is less than 1. This effect must be very small since it is not visible in Fig 1. But it also means that the two curves in that figure have to become vertical before f=1, close or not to the point where the dependence on f becomes super exponential.

4: "In other words, the dynamical system reaches different fixed points for I = 0 and 1." You can improve this sentence so that the readers will understand your point better. The fixed points at I=0 and I=1 are different even at small f. What happens at larger f values is that bistability emerges: in a range of intermediate I values, one has 2 fixed points, having respectively low and high outputs. So although individual genes do not have bistability (for the beta and mu parameters you chose), the GRN does.

5: The high fitness GRNs correspond to amplifying the differential input of node 0. An engineer would be able to design a set of Jijs to do that "simply". Are you able to give such a "design" interpretation to your GRNs of high fitness?

6: Out of curiosity, have you understood whether the autoactivation plays any role, i.e. would the results be almost the same without Jii=+1? I assume that all the high fitness GRN have Jii=0 or +1, that is Jii=-1 is "not good" for fitness. Any results there will be welcome.

7: "This implies that the difference between evolution and random sampling should lie in global network structures." I assume you are comparing the GRNs at a given fitness. Then if it is fitness that drives motifs one does not expect much differences between the two ensembles. The motifs should appear as fitness rises, that is probably the relevant correlation, while the "evolutionary dynamics" ensemble is mainly relevant for mutational robustness.

8: In the Discussion, you write "The new GRN phenotype, bistability, inevitably emerges as fitness increases". I would be more careful about the claim, although the observation holds in your system, it may not be "inevitable". For instance if the parameter mu were adjustable for each gene, one could probably just use a multi-layer fan out architecture to amplify the deviation from 0.5 of each gene's input signal and there would be no bistability. As a result, I am not so sure that your interpretation of "delayed" phenotypic change is generic. Specifically, I am not sure that mutational robustness (enhanced under evolution) is antagonistic with having a particular phenotype (here multistability). This connects to my point 3 in "Main recommendations of the referee" where I think you can improve your study to really reach a conclusive result.

Reviewer #2: In this computational study, the authors simulate a simple model of a genetic regulatory network to investigate how the ensemble of genotypes emerging through an evolutionary process is atypical relative to the ensemble of all high-fitness genomes. For this, they capitalize on a method described in a recent paper sharing one of the co-authors, which employed the multicanonical Monte Carlo technique to sample genotypes in a given fitness bin in an evolution-independent way. This approach allows comparing evolved high-fitness genomes to randomly sampled high-fitness genomes. The authors show that the evolved genomes are more mutationally robust than typical.

I would not necessarily call this reported result groundbreaking; it confirms the intuition many people would probably say they have, but it does so in a clear quantitative way, and the technique used here will likely be useful to other studies. (I, for one, was pleased to learn about the capabilities of McMC.) The paper is well-written and, barring a few comments below, clearly explained. I though the introduction could probably do a better job integrating into the context of the existing literature: currently, it focuses rather narrowly on mutational robustness & GRN, but the idea that evolved genomes are atypical in potentially predictable ways recurs in the broader literature, and to me it seems that the authors' approach is relevant to this broader conversation, which might help this work gain broader visibility. Some examples I'm familiar with are Orr, Genetics (2003); Rice Good Desai, Genetics (2015); Tikhonov Kachru Fisher, PNAS (2020); the work of Erik van Nimwegen. I'm not saying these specific references should be added, merely that these papers, or references therein, might inspire the authors to add a broader paragraph to their discussion section, saying where a similar McMC methodology could be useful.

A few more specific comments:

Major:

1. The "model" section is completely silent on the details of how the evolutionary simulations were performed. Was a single lineage simulated at a time, or a population (how large?) How were mutations implemented? (This is mentioned very late in the paper, and only because it helps explain some finding.) The evolutionary procedure must be clearly described in the methods section.

2. "Evolution behaves conservatively" - I don't know what this means. Evolved genomes are atypical; various observables are skewed in one or another direction, and one can invoke intuitive arguments to rationalize these differences if it helps explain the observations. But interpreting this as "evolution behaves conservatively" / "evolution prevents easy phenotypic changes" does not strike me as justified.

3. In the concluding paragraph, the authors say that one of their "three main results" is the "relationship between evolutionary speed and genetic entropy". The value of this "genotypic entropy" discussion is unclear to me. As you approach a fitness peak, the further increase slows down. This is neither new nor surprising, and if the authors choose to dedicate a paragraph of the discussion to this (e.g. lines 301-321), I think it should be made more clear what specific point they are trying to make. As it is, the discussion is rather meandering and strangely singles out one particular study by Sato et al.

4. "If a difference in x_out larger than 10−6 was observed between these two processes in some range of I, the GRN was regarded as bistable." <-- The value 10^-6 seems so incredibly low that I'm wondering if this is a typo. Surely that level of "bistability" would be completely irrelevant for any biological system? I would imagine (or certainly hope) that the findings in the paper should not change if the threshold is set to something more biologically reasonable. Any value here is an arbitrary choice and I do of course understand that this is not meant as a practical model of any real GRN, but still - any value below e.g. 1% seems like a very odd choice. If the findings DO depend on the threshold value, then I would find this concerning and would certainly like to see some discussion of this.

Minor:

Line 113 - I think some reference or mention that this curve was generated by the McMC protocol is necessary also in the text, not just the figure legend.

Line 117 "GRNs decrease exponentially with f" --> number of GRNs?

"Fitness landscape is rough" -> Is this referring to the wiggles in the curve? A rough/rugged landscape usually means something quite specific, and this doesn't seem to be the same meaning as here.

Line 120 "it was divided into three regions" -> unclear what this is referring to

Are there any error bars in Fig 2? Otherwise how can we tell that it is substantially different?

Line 178: Increased monotonically

**********

Have all data underlying the figures and results presented in the manuscript been provided?

Large-scale datasets should be made available via a public repository as described in the PLOS Computational Biology data availability policy, and numerical data that underlies graphs or summary statistics should be provided in spreadsheet form as supporting information.

Reviewer #1: Yes

Reviewer #2: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: Yes: Olivier C. Martin

Reviewer #2: No

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, . PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at .

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example in PLOS Biology see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, we recommend that you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. Additionally, PLOS ONE offers an option to publish peer-reviewed clinical study protocols. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols

16 Aug 2021

Submitted filename: reply_to_reviewers.pdf

Click here for additional data file.

23 Sep 2021

Dear Dr. Kikuchi,

Thank you very much for submitting your manuscript "Evolution enhances mutational robustness and suppresses the emergence of a new phenotype" for consideration at PLOS Computational Biology.

As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by several independent reviewers. In light of the reviews (below this email), we would like to invite the resubmission of a significantly-revised version that takes into account the reviewers' comments.

We cannot make any decision about publication until we have seen the revised manuscript and your response to the reviewers' comments. Your revised manuscript is also likely to be sent to reviewers for further evaluation.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to the review comments and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Please prepare and submit your revised manuscript within 60 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email. Please note that revised manuscripts received after the 60-day due date may require evaluation and peer review similar to newly submitted manuscripts.

Thank you again for your submission. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Alexandre V. Morozov, Ph.D.

Associate Editor

PLOS Computational Biology

Sushmita Roy

Deputy Editor

PLOS Computational Biology

***********************

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: General overview

The authors have worked hard to address the criticisms of both referees and the revised version is definitely of quality.

Minor recommendations of the referee

(1) In your abstract you write "This implies that the delayed emergence of bistability is a consequence of the mutation-selection mechanism of evolution". Although I am happy you were able to confirm my intuition on this issue, it is not a proof. You should replace "implies" in this sentence by "strongly suggests".

(2) The term "universal" in physics is strong and I recommend you remove it and tone down the "universality" claims you make in this paper.

(3) The term "randomly sampled GRNs" will be too often interpreted by readers as uniform with no constraint on fitness. I think for your purposes you could replace it by "GRNs sampled uniformly at given fitness" and possibly refer to this as your "reference ensemble" of GRNs. (I prefer "reference ensemble" to "reference set" as it is more general and is directly taken from statistical physics.

(4) "E. Coli" -> "E. coli"

(5) You would benefit from having a native English speaker go through the text.

(6) I found the "Summary and discussions" section rather long.

Other comments

To follow up on my previous report, perhaps you did not understand my sugestions concerning temperature. Your use of different fractions of the population to duplicate corresponds to different selection pressures. I was implicitly suggesting to use the more standard framework whereby the fitness of an individual gives a probability of reproduction. Generally fitness is written as

f = exp(g(genotype)) and the probability to reproduce is f up to a factor depending on the population average of f. This standard method corresponds to Fisher-Wright and allows one to keep the population size constant. My suggestion to introduce a temperature just corresponds to changing that rule to

f = exp(g(genotype)/T). A given T produces a steady state ensemble that will be unique, thereby avoiding (at least theoretically) the disconnected components of a neutral network that you find with your approach.

Reviewer #2: Please see attachment

**********

Have the authors made all data and (if applicable) computational code underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data and code underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data and code should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data or code —e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, . PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at .

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example in PLOS Biology see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, we recommend that you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. Additionally, PLOS ONE offers an option to publish peer-reviewed clinical study protocols. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols

Submitted filename: Review post-revisions.pdf

Click here for additional data file.

2 Nov 2021

Submitted filename: reply_to_reviewers.pdf

Click here for additional data file.

4 Dec 2021

Dear Dr. Kikuchi,

Thank you very much for submitting your manuscript "Evolution enhances mutational robustness and suppresses the emergence of a new phenotype: A new numerical approach for studying evolution" for consideration at PLOS Computational Biology. As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by several independent reviewers. The reviewers appreciated the attention to an important topic. Based on the reviews, we are likely to accept this manuscript for publication, providing that you modify the manuscript according to the review recommendations.

Please prepare and submit your revised manuscript within 30 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to all review comments, and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Thank you again for your submission to our journal. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Alexandre V. Morozov, Ph.D.

Associate Editor

PLOS Computational Biology

Sushmita Roy

Deputy Editor

PLOS Computational Biology

***********************

A link appears below if there are any accompanying review attachments. If you believe any reviews to be missing, please contact ploscompbiol@plos.org immediately:

[LINK]

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: General overview

In my previous report I no longer had major scientific recommendations. In their revision, the authors handled properly my recommendations and they furthermore improved the logical flow of the paper and streamlined it, taking it from a meadering if not branched and dispersed presentation to a coherent linear demonstration.

I also went through their response to the second referee. Part of his/her reservations overlapped with mine. Since this revision is now readable and follows a clear logic, in my opinion the authors have also responded satisfactorily to the recommendations of this second referee, though clearly he/she had the same reading effort to do as me.

Minor comments

1: In the title and elsewhere, you could replace "numerical" by "computational".

2: I would not use the term "lethal edge": it is the deletion of the edge that is lethal...

3: Figure 7 is a bit misleading because the orange points are displayed after the blue ones and so hide them. It might be better to use crosses for blue rather than the circles so that the overlaps are not a problem.

4: Although Editage was useful, they did not do as good a job on the abstract and introduction (which is overly "wordy"). If you could improve the English there, you'll have more readers going to the end of your paper.

Reviewer #2: I thank the authors for the revisions. They have addressed my concerns. Hopefully the authors find that the changes did indeed benefited clarity.

**********

Have the authors made all data and (if applicable) computational code underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data and code underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data and code should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data or code —e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: Yes: Olivier C. Martin

Reviewer #2: No

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org.

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example in PLOS Biology see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, we recommend that you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. Additionally, PLOS ONE offers an option to publish peer-reviewed clinical study protocols. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols

References:

Review your reference list to ensure that it is complete and correct. If you have cited papers that have been retracted, please include the rationale for doing so in the manuscript text, or remove these references and replace them with relevant current references. Any changes to the reference list should be mentioned in the rebuttal letter that accompanies your revised manuscript.

If you need to cite a retracted article, indicate the article’s retracted status in the References list and also include a citation and full reference for the retraction notice.

16 Dec 2021

Submitted filename: Reply_to_reviewers.pdf

Click here for additional data file.

27 Dec 2021

Dear Dr. Kikuchi,

We are pleased to inform you that your manuscript 'Evolution enhances mutational robustness and suppresses the emergence of a new phenotype: A new computational approach for studying evolution' has been provisionally accepted for publication in PLOS Computational Biology.

Before your manuscript can be formally accepted you will need to complete some formatting changes, which you will receive in a follow up email. A member of our team will be in touch with a set of requests.

Please note that your manuscript will not be scheduled for publication until you have made the required changes, so a swift response is appreciated.

IMPORTANT: The editorial review process is now complete. PLOS will only permit corrections to spelling, formatting or significant scientific errors from this point onwards. Requests for major changes, or any which affect the scientific understanding of your work, will cause delays to the publication date of your manuscript.

Should you, your institution's press office or the journal office choose to press release your paper, you will automatically be opted out of early publication. We ask that you notify us now if you or your institution is planning to press release the article. All press must be co-ordinated with PLOS.

Thank you again for supporting Open Access publishing; we are looking forward to publishing your work in PLOS Computational Biology.

Best regards,

Alexandre V. Morozov, Ph.D.

Associate Editor

PLOS Computational Biology

Sushmita Roy

Deputy Editor

PLOS Computational Biology

***********************************************************

13 Jan 2022

PCOMPBIOL-D-21-00210R3

Evolution enhances mutational robustness and suppresses the emergence of a new phenotype: A new computational approach for studying evolution

Dear Dr Kikuchi,

I am pleased to inform you that your manuscript has been formally accepted for publication in PLOS Computational Biology. Your manuscript is now with our production department and you will be notified of the publication date in due course.

The corresponding author will soon be receiving a typeset proof for review, to ensure errors have not been introduced during production. Please review the PDF proof of your manuscript carefully, as this is the last chance to correct any errors. Please note that major changes, or those which affect the scientific understanding of the work, will likely cause delays to the publication date of your manuscript.

Soon after your final files are uploaded, unless you have opted out, the early version of your manuscript will be published online. The date of the early version will be your article's publication date. The final article will be published to the same URL, and all versions of the paper will be accessible to readers.

Thank you again for supporting PLOS Computational Biology and open-access publishing. We are looking forward to publishing your work!

With kind regards,

Olena Szabo

PLOS Computational Biology | Carlyle House, Carlyle Road, Cambridge CB4 3DN | United Kingdom ploscompbiol@plos.org | Phone +44 (0) 1223-442824 | ploscompbiol.org | @PLOSCompBiol