# Evolution enhances the mutational robustness and suppresses emergence of a
new phenotype

Tadamune Kaneko  Department of Physics, Osaka University  Cybermedia Center,
Osaka University  Macoto Kikuchi<EMAIL_ADDRESS>Cybermedia Center,
Osaka University  Department of Physics, Osaka University  Graduate School of
Frontier Bioscience, Osaka University

###### Abstract

Mutational robustness is an indispensable trait that living systems must
possess. This property developed through a long history of evolution. In this
paper, we discuss the evolutionary enhancement of mutational robustness for a
toy model of the gene regulatory network (GRN). For exploring the
particularity of evolution, we need a reference system. Such a reference
system that we consider appropriate is a set of randomly sampled GRNs.
However, since GRNs with high fitness are rare, a simple random sampling
method is of no use for this purpose. Recently, Nagata and Kikuchi (PLoS
Comput Biol 16 (2020) e1007969) proposed an efficient procedure for randomly
sampling GRNs in a wide range of fitness based on the multicanonical Monte
Carlo method. We employ that method to construct the reference set and compare
it with the lineages obtained by evolutionary simulations. We consider a
neural-network-like model of GRN, which has one input gene and one output
gene. The fitness is defined so that it is large when the difference in the
expressions of the output gene between two input values (“on” and “off”)
becomes large. By comparing the evolution of fitness and the fitness
distribution of the reference set, we find that evolution slows down
significantly near the value of fitness that the number of GRNs becomes
extremely rare. It implies that the evolutionary speed is determined mainly by
the number of available GRNs, or in other words,“entropy”. Random sampling
reveals that the highly-fit GRNs in this model inevitably exhibit bistability.
Accordingly, the evolutionary simulations produce bistable GRNs as fitness
becomes high. But the increase in the fraction of bistable GRNs occurs at
higher fitness than the reference set. Thus, the emergence of a new phenotype,
bistability, is delayed in evolution. We also study an indicator for
mutational robustness. While in the early stage of evolution, it coincides
with that of the reference set, in the late stage, it becomes significantly
higher than that of the reference set. In short, evolution enhances mutational
robustness. The reason is that GRNs with fewer lethal edges, that is, GRNs
having fewer such regulatory interactions that the fitness drops largely when
cut, are selected by evolution. We also find that local network motifs are
relevant to function but irrelevant to mutational robustness. Thus mutational
robustness is a property related to a global structure of GRNs. While those
results are intuitively natural, we show them quantitatively by our method.

##  I Introduction

Living systems have developed their specific functions through a long history
of evolution. Not only adaptation to the environment, but they have also
acquired many kinds of robustness, such as the one against environmental
noises as well as internal noises and the one against mutation of
genomes.Kitano (2004); Wagner (2005); Masel and Siegal (2009) Among them, the
most important is the mutational robustness, by which a living system does not
lose its function and can continue its life even when the genome is mutated.

Mutational robustness has been demonstrated experimentally. For example,
comprehensive single-gene knock-out experiments for bacterium and yeasts have
revealed that most of the knock-outs do not affect viability.Giaever  _et al._
(2002); Baba  _et al._ (2006); Kim  _et al._ (2010) An artificial rewiring
experiment for the gene regulatory network of E. Coli also has shown that the
bacteria is robust for most of the artificial additions of regulatory
links.Isaran  _et al._ (2008) Since living systems are subjected constantly to
genetic mutations, if they easily lost their functions, such systems could not
have produced offsprings during the evolutionary process. Thus, the mutational
robustness must have also evolved through the evolutionary process along with
the function. The central question we pose in this paper is whether the
mutational robustness evolves.

However, as Wagner pointed out,Wagner (2005) experimental studies on this
issue would face difficulty more or less, because the existing living
organisms are already the products of evolution, and evolutionary experiments
can provide only outcomes reflecting the evolutionary pathways. Therefore,
numerical methods play indispensable roles in obtaining information. In this
study, we conducted a numerical investigation on the evolution of the
mutational robustness using a toy model of the gene regulatory networks.

Genes coded on DNA in cells are read by RNA polymerase and transcribed to
mRNA, whose information is then used to assemble proteins. A category of
proteins called the transcription factors (TF) act as activators or repressors
for other genes. Many genes regulate each other in this way and form a complex
network called the gene regulatory network (GRN). GRNs are utilized to change
the cell states to adapt to environmental change or to control the cell cycles
or cell differentiation. One of the best known mathematical models of GRN is
the boolean network model proposed by Kauffman, in which each fixed point of
the dynamical system are considered to represent a cell state.Kauffman (1969);
Kauffman and Levin (1987); Kauffman (1993).

To investigate whether the mutational robustness evolves or not, we need some
reference system independent of evolutionary pathways. In the case of GRN, the
reference system that we consider appropriate is the set of randomly-sampled
GRNs. Random sampling is suitable to find universal properties that do not
depend on histories of evolution. In this direction of research, Ciliberti et
al. conducted a random sampling of GRNs to discuss mutational
robustness.Ciliberti  _et al._ (2007a, b) The fitness of their model was only
two-valued: viable and non-viable. They investigated the interrelations of
viable GRNs and argued that a majority of them belong to a large cluster
connected by neutral mutations. However, such a simple random sampling method
is not useful for systems having a more complex fitness landscape because
highly-fit GRNs are very rare. Then, Burda et al. and Zagorski et al. employed
the Markov-Chain Monte Carlo method to sample highly-fit GRNs.Burda  _et al._
(2011); Zagorski  _et al._ (2013) They found that GRNs exhibiting multi-
stability contain some common network motif. These researches focused on the
universal properties of highly-fit GRNs.

The above methods are insufficient to sample GRNs in a wide range of fitness.
Saito and Kikuchi proposed to utilize the Multicanonical Monte Carlo method
(McMC) to investigate the mutational robustness of GRNs.Saito and Kikuchi
(2013) McMC has been developed originally in the field of statistical physics
for sampling configurations within a wide range of energy evenly.Berg and
Neuhaus (1991, 1992) But It has been realized to be useful also in sampling
for nonphysical systems. It is particularly effective for generating very rare
states and estimating the appearance probabilities of such rare states.Saito
_et al._ (2010); Saito and Iba (2011); Kitajima and Iba (2011); Iba  _et al._
(2014); Kitajima and Kikuchi (2015)

By using McMC, Nagata and Kikuchi investigated a model of GRN.Nagata and
Kikuchi (2020) They regarded fitness as the “energy” of GRNs and sampled GRNs
randomly from very low fitness to very high fitness, and then classified them
by fitness to explore their universal properties that are independent of the
evolutionary pathway. Put concretely, they considered a neural-network-like
model of GRN with one input gene and one output gene and set the fitness so
that it is high if the GRN responded to change of input sensitively. Since
each gene in their model did not respond ultrasensitively,Goldbeter and
Koshland (1981); Ferrell and Ha (2014a) and the network structures were
restricted, simple network motifs did not give rise to bistability.Alon (2020)
Despite that setting, they found that highly-fit GRNs always exhibited
bistability. Thus, this bistability emerges as a consequence of the
cooperation of many genes. According to their results, the new phenotype that
exhibits bistability appears whatever evolutionary path the system follows. In
this case, the fitness function constrains the possible phenotype. Adding to
that, they found that mutationally robust GRNs are not rare among such highly-
fit GRNs. Thus the intuitive expectation that the optimization process would
lead to mutational fragility is not necessarily correct for evolution of this
model. Then, does evolution further enhance the mutational robustness?

To answer this question, we conducted evolutionary simulations. By comparing
these results with the reference set of GRNs obtained by the method of Ref.
Nagata and Kikuchi (2020), we discuss the universality and the characteristics
of evolution.

##  II Model and method

To model GRN, we take into account only regulatory relations ignoring the
details of gene expression. Such connectionist-type modeling has been used
widely in theoretical studies.Mjolsness  _et al._ (1991); Wagner (1994, 1996);
Siegal and Bergman (2002); Masel (2004); Kaneko (2007); Espinosa-Soto  _et
al._ (2011); Inoue and Kaneko (2013, 2018); Nagata and Kikuchi (2020) We
represent a GRN as a directional graph by identifying nodes as genes and edges
as regulatory interactions. We consider GRNs with one input node and one
output node. In contrast to Ref. Nagata and Kikuchi (2020) that several
restrictions were imposed to the network structure, we allow any network as
long as the number of edges from one node to another is at most one. Thus,
GRNs can contain auto-regulation and mutual regulation. Genes and edges that
do not contribute to output are also permitted. As a result, the effective
network size becomes smaller in general than the prepared size. In the
following sections, we restrict ourselves to the networks of the number of
nodes $N=32$ and the number of edges $K=80$.

A variable $x_{i}\in[0,1]$ that represents the expression of a gene is
assigned to each node with $i$ the number of the gene. $x_{i}$ obeys the
following discrete-time dynamics: Wagner (1994, 1996)

$x_{i}(t+1)=R\left(I\delta_{i,0}+\sum_{j}J_{ij}x_{j}(t)\right),$ (1)

where $t$ is the time step. $J_{ij}$ represents the regulation from the $j$-th
gene to the $i$-th gene. For simplicity, we assume that $J_{ij}$ takes one of
the three values, $0,\pm 1$; $+1$ indicates activation, $-1$ indicates
repression, and 0 means the absence of regulation. The $0$th gene is the input
gene, and $I\in[0,1]$ is the strength of the input signal.Inoue and Kaneko
(2013) The response of the genes is given by the following sigmoidal function:

$R(x)=\frac{1}{1+e^{-\beta(x-\mu)}},$ (2)

where $\beta$ represents the steepness of the function, and $\mu$ is the
threshold. This function is used widely in theoretical studies.Siegal and
Bergman (2002); Furusawa and Kaneko (2008); Pujato  _et al._ (2013); Inoue and
Kaneko (2013, 2018) In the present study, we set $\mu=0$ and $\beta=2$. These
parameters are the same as those used for the neural network model by Hopfield
and TankHopfield and Tank (1985) and give a rather gradual increase to the
response function, which reflect the stochastic nature of gene
expression.Inoue and Kaneko (2018) The response function $R(x)$ is not
ultrasensitive in the sense that a single gene with this response function has
a single fixed point even when the auto-activation loop is attached.
Therefore, for multistability to emerge, many genes are necessary to act
cooperatively. Although the spontaneous expression $R(0)=0.5$ is rather large,
we do not consider it causes any problem for the present purpose.

We define fitness for this model. Following Ref. Nagata and Kikuchi (2020), we
set the fitness so that it becomes large when the difference in the expression
of the output gene between $I=0$ and $1$ (“off” and “on”) is large. Let
$\bar{x}_{out}(I)$ be the corresponding output to input value $I$.
$\bar{x}_{out}(I)$ is the fixed-point value of the output gene for fixed $I$;
as the initial condition, we set the expression of all genes as 0.5, the
spontaneous expression. If $x_{out}(t)$ behaves oscillatory in time instead of
reaching the fixed point, we use the temporal average for $\bar{x}_{out}(I)$;
the system reaches a fixed point in most cases. The fitness $f$ is defined as
follows:

$f=|\bar{x}_{out}(0)-\bar{x}_{out}(1)|.$ (3)

The fitness takes values in $[0,1]$ by definition.

We studied this model using two different methods. One is random sampling
using the multicanonical Monte Carlo method (more precisely, entropic
sampling,Lee (1993) which is one of the variations of the McMC). Details of
the method are described in Ref. Nagata and Kikuchi (2020). We divided fitness
into 100 bins and determined the weight for each bin so that the GRNs appeared
evenly using the Wang-Landau method.Wang and Landau (2001a, b) One update of
McMC consisted of the following successive two processes: Deleting a randomly
selected edge and putting a new edge to an unlinked gene pair that is also
chosen randomly. Whether accepting this change or not was decided by the
Metropolis method. One Monte Carlo step (MCS) consisted of $K$ such updates,
and we recorded $f$ at every MCS. We sampled GRNs at every 20 MCSs to reduce
the correlation between the samples. We conducted ten independent runs each of
which consisting of $10^{7}$ MCSs. As a result, we obtained 50000 samples on
average for each bin. Ideally, we can gather GRNs randomly within each bin.
Although inter-sample correlations remain to some extent, we call this method
“the random sampling” in this paper. The set of these randomly sampled GRNs
provides the reference set.

Another method is the evolutionary simulation. As the initial population, we
generated 1000 GRNs randomly. In one step of evolution, we selected the top
500 GRNs according to the order of fitness and made one copy for each; then,
these copies were subjected to mutation. The mutation operation for each GRN
was the same as the procedure in the McMC; an edge was deleted randomly from
the network, and then a new edge was added between a randomly selected
unlinked gene pair. As these procedures repeated, the average fitness of the
population increased. After 150 generations, we picked the GRN having the
highest fitness among the population and traced its ancestors to obtain a
single lineage. The above population size and the number of generations were
arbitrarily determined. We repeated this evolutionary simulation 10000 times
independently, and as a result, we collected 10000 lineages.

##  III Results

###  III.1 Fitness landscape and speed of evolution

Figure 1: Fitness landscape and evolution of fitness. Fitness $f$ is divided
into 100 bins. The blue dashed lines (left axis) show the appearance
probability of each bin $\Omega(f)$, which was obtained by McMC. The orange
solid lines (right axis) is the average fitness of each generation by the
evolutionary simulation. The averages are taken over 10000 lineages. (a)
$\beta=2$ (b) $\beta=4$.

The blue lines in Fig. 1(a) show the appearance probability $\Omega(f)$ of
each bin of fitness $f$. It can be called the fitness landscape. The sum of
$\Omega(f)$ is normalized to 1. The probability for $f\geq 0.99$ is $\sim
10^{-16}$. Since we can count the total possible number of GRNs as
$\binom{N^{2}}{K}2^{K}\simeq 10^{145}$, the highly-fit GRNs are numerous but
rare. Majority of GRNs concentrate near $f\sim 0$. GRNs decrease exponentially
with $f$, and for very high $f$, they decrease faster than exponential. For
comparison, we show the fitness landscape for steeper response function,
$\beta=4$, in Fig. 1(b). Although the fitness landscape is rough (This
roughness is not a statistical error), it is divided into three regions as in
the case of $\beta=2$.

The orange lines are the average fitness of each generation obtained by the
evolutionary simulation. The vertical axis represents the generation in the
downward direction. The evolution progresses almost linearly in the early
stage and slows down drastically for large $f$. The values of $f$ that the
increase in fitness starts to slow down roughly coincide with the value that
the faster-than-exponential decrease in $\Omega(f)$ begins for both $\beta=2$
and $4$. The reason is intuitively understandable: The number of the possible
destination that a GRN can transit by the mutation is restricted by
$\Omega(f)$. In other words, when the number of GRNs with higher fitness
decreases drastically, the possibility that the fitness raises by chance also
decreases. Since $\log\Omega(f)$ can be called “entropy”, this result suggests
that the evolutionary speed depends in large part on the entropy.

###  III.2 Delayed emergence of bistability

The model in Ref. Nagata and Kikuchi (2020) exhibits bistability as $f$
approaches its maximum value. Namely, the dynamical system reaches different
fixed points for $I=0$ and $1$. When $I$ is changed continuously, two saddle-
node bifurcations take place,Tyson  _et al._ (2003) in contrast to the case of
small $f$ that a single fixed point moves to follow the change of $I$. We
found that the present model behaves similarly.

The bistable GRNs are classified into three categories. The first is the
toggle switch that two saddle-node bifurcations occur within $I\in[0,1]$. The
toggle switch is found, for instance, in phage $\lambda$ and utilized for
adaptation to environmental change.Ptashne and Gann (2002); Ptashne (2004);
Zong  _et al._ (2010); Watson  _et al._ (2013) The second is the one-way
switch.Tyson  _et al._ (2003) In this case, only one saddle-node bifurcation
point is found in the range $I\in[0,1]$, and another bifurcation point is out
of this range. One-way switches give rise to cell maturation or cell
differentiation; a typical example is the MAPK cascade in the maturation of
the Xenopus oocyte.Ferrell and Machleder (1998) In the last category, both of
two bifurcation points are out of the range of $I$. Since this type may not
work as a switch as long as the effect of noise is not taken into account, we
call it “unswitchable”. In this paper, we do not distinguish these three and
treat them equally as the bistable GRNs. It is straightforward to change the
definition of bistability to deal only with, for example, toggle switches.

Figure 2: Fitness dependence of the fraction of bistable GRNs. The orange
solid line is for evolutionary simulation, and the blue dashed line is for
random sampling. Bistability is checked as follows: First, 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 up to $I=1$.
Then, the inverse process from $I=1$ to $0$ was performed. If the difference
in $\bar{x}_{out}$ larger than $10^{-6}$ was observed between these two
processes in some range of $I$, the GRN was regarded as bistable. GRNs
obtained by evolutionary simulations were classified according to $f$ into 100
bins similarly to random sampling. $P_{2}(f)$ was calculated for every other
bin.

Fig.2 shows the fraction of the bistable GRNs $P_{2}(f)$ against fitness. The
blue line is the result of random sampling. The bistable GRNs starts to appear
at $f\simeq 0.5$ and increases rapidly until all the GRNs become bistable for
$f\rightarrow 1$. Thus, the new phenotype that exhibits bistability emerges as
the fitness increases, and all the GRNs converge to such a phenotype as the
fitness approaches its maximum value. The orange line is $P_{2}(f)$ obtained
by evolutionary simulation. We classified GRNs in the obtained 10000 lineages
by fitness as in the case of the random sampling. The figure clearly shows
that $P_{2}(f)$ of the evolutionary simulation is significantly low compared
to that of the random sampling at the same $f$. In other words, evolution
delays the rapid increase of $P_{2}(f)$. It suggests that evolution tends to
suppress the emergence of a new phenotype. The reason is that a lineage
reflects its evolutionary history; nonetheless, the eventual emergence of the
bistable phenotype is inevitable as the fitness increases because all the GRNs
are bistable for $f\rightarrow 1$.

###  III.3 Evolution of mutational robustness

To discuss mutational robustness, we introduce a measure for robustness. In
Ref. Nagata and Kikuchi (2020), the single-edge deletion is considered as a
mutation; it was found that the edges split into two classes, neutral and
lethal, for highly-fit GRNs, and the lethal edges are minor among them. We
consider the single-edge deletion as the mutation also for the present model
and obtained the same results. Then, we defined the following quantity $r$ as
the robustness measure:

$r\equiv\frac{1}{K}\sum_{i=1}^{K}f^{\prime}_{i},$ (4)

where $f^{\prime}_{i}$ is fitness after the $i$th edge is deleted, and the sum
is taken for all the edges. Since $f^{\prime}$ should have a tendency to
increase with $f$, comparing $r$ of GRNs with different $f$ is not meaningful.
However, by comparing this quantity for GRNs obtained by random sampling and
evolutionary simulation at the same $f$, we can discuss the particularity of
evolution.

Figure 3: Average of the robustness measure $r$ against fitness. The average
was taken over all the samples in each bin. The orange solid line is for
evolutionary simulation, and the blue dashed line is for random sampling. The
averages for evolutionary simulation were calculated by the reweighting method
described in the text.

Fig3 shows the average of $r$ against $f$. The average was taken over all the
GRNs in the corresponding bin. $\langle r\rangle$ obtained by evolutionary
simulation increases monotonously and coincides with that obtained by random
sampling up to $f\simeq 0.5$. For larger $f$, The value of evolution departs
upward from that of random sampling, and the difference becomes more and more
significant as $f$ increases.

Figure 4: Probability distribution of the robustness measure $r$. The orange
solid line is for evolutionary simulation, and the blue dashed line is for
random sampling. The distribution for evolutionary simulation were corrected
by the reweighting method described in the text. (a)$f\in[0.5,0.51)$
(b)$f\in[0.8,0.81)$ (c)$f\in[0.99,1.0]$

To scrutinize the difference, we show the probability distribution of
$f^{\prime}$ for $f=0.5,0.8$, and $0.99$ in Figs. 4(a)-(c), respectively.
Here, the value of $f$ indicates the lower bound of the bin; for example,
$f=0.99$ means $f\in[0.99,1.0]$. Taking into account that the distribution of
$f$ in each bin differs for random sampling and evolution, we divided each bin
into ten sub-bins and reweighed the distribution obtained by evolutions so
that the distribution of $f$ coincides with that of random sampling (The
values plotted in Fig.3 shown before were also calculated based on this
reweighing method). We divided $f^{\prime}$ into 80 bins because otherwise,
unnecessary oscillation appears in the distribution due to the discreteness of
the number of lethal edges. While both distributions roughly agree for
$f=0.5$, we can see deviation for $f=0.8$. The two distributions exhibit
distinct difference for $f=0.99$; the distribution obtained by evolution is
biased to large $r$ compared to random sampling. Therefore, evolution enhances
mutational robustness.

What causes this difference? As was stated before, the edges of randomly
sampled GRNs for $f\geq 0.99$ split into neutral and lethal. Namely, when we
cut one edge, $f^{\prime}$ is either close to $f$ or almost 0, and other types
of edges are scarce. Based on this observation, we investigated the number
distribution of the lethal edges by defining that an edge is lethal if
$f^{\prime}<0.8$ . Although this definition is arbitrary, it hardly affects
the result. Fig.5 shows the probability distribution of the number of the
lethal edges $n_{L}$ for $f\geq 0.99$. GRNs obtained by evolutionary
simulation are biased significantly to the small number side compared to
random sampling. The highest probability is at $n_{L}=3$, and the distribution
is narrow; moreover, more than 2% of GRNs have no lethal edge. In contrast,
the highest probability for random sampling is at $n_{L}=15$, and the
distribution is much wider; GRNs lacking lethal edge are only 0.4%. Thus, it
is clear that the smallness of the number of lethal edges is the reason for
the enhancement of mutational robustness by evolution.

Figure 5: Probability distribution of the number of lethal edges $n_{L}$ for
$f\in[0.99,1.0]$. The orange solid line is for evolutionary simulation, and
the blue dashed line is for random sampling. An edge is regarded as lethal if
fitness after that edge is cut is less than $0.8$. The distribution for
evolutionary simulation were corrected by the reweighting method described in
the text.

One possibility is that the nodes that affect the output are fewer in the
evolutionary obtained GRNs. To check this, we counted the number of nodes
$n_{N}$ that have at least one path to the output node. Fig.6 shows the
distributions; that for the random networks is also plotted for comparison.
While the peak of the distribution is at $n_{N}=28$, which is slightly fewer
than the random networks, the distributions are indistinguishable between
evolution and random sampling. Therefore, the number of effective nodes is not
the reason for the difference in the number of lethal edges.

Figure 6: Distribution of effective node number $n_{N}$ for $f\in[0.99,1.0]$.
If a node has at least one path to the output node, that node is regarded as
“effective”. The orange solid line is for evolutionary simulation, and the
blue dashed line is for random sampling. The distribution for the random
networks is also shown as a reference by the gray dotted line. The
distribution for evolutionary simulation were corrected by the reweighting
method described in the text.

###  III.4 Motif analysis

Next, we investigate network motifs. In the model of Ref. Nagata and Kikuchi
(2020), the coherent feedforward loops and the positive feedback loops were
significantly abundant in highly-fit GRNs. Since the present model allows
auto- and mutual-regulations, we also explored other patterns adding to the
triangular patterns. Put concretely, we counted the number of auto-
regulations, mutual-regulations between gene pairs, triangles, and mutual
activation or repression of two nodes combined with auto-activation of both
nodes. The last pattern frequently appears in multistable GRNs.Ptashne and
Gann (2002); Ptashne (2004); Huang  _et al._ (2007); Burda  _et al._ (2011);
Jia  _et al._ (2017); Huang  _et al._ (2017) Since the motifs are defined as
the network patterns that are abundant compared to the random networks,Shen-
Orr  _et al._ (2002); Alon (2020) we also counted them for the random
networks.

Figure 7: Probability distributions of the number of motifs for
$f\in[0.99,1.0]$. The orange solid lines are for evolutionary simulation, and
the blue dashed lines are for random sampling. The distributions for the
random networks are also shown as references by the gray dotted lines. The
distributions for evolutionary simulation were corrected by the reweighting
method described in the text.(a)auto-activation (b)mutual activation (c)mutual
repression (d)coherent feed-forward loop (e)positive feedback loop (f)mutual
activation or mutual repression accompanied by auto-activation of both genes

As a result, the following patterns are more in number than the random
networks: auto-activation, mutual-activation, mutual-repression, coherent
feedforward loops (FFL+), positive feedback loops (FBL+), and mutual
activation or repression plus auto-activation of two genes. We show their
number distributions in Fig.7(a)-(f). Although their abundances are not so
remarkable, we may call them motifs. Other patterns such as auto-repression,
incoherent feedforward loop, and negative feedback loop are less than the
random networks. The notable point here is that the distributions of these
motifs are almost the same for both evolution and random sampling. Thus,
whether or not these highly-fit GRNs are products of evolution is not
reflected in distribution of these local motifs. Therefore, the difference
between evolution and random sampling should lie in global network structures.

###  III.5 Path distribution

Figure 8: Probability distribution of the number of paths $n_{path}$ starting
from the input node and reaching the output node without passing the same node
twice. $n_{path}\leq 400$ is shown. The orange line is for evolutionary
simulation, and the blue line is for random sampling. For random sampling,
$P(1)$ reaches 4%. Figure 9: Scatter plot of the number of paths $n_{path}$
and the number of the lethal edges $n_{L}$. $n_{path}\leq 1000$ is shown. The
orange dots are for evolutionary simulation, and the blue dots are for random
sampling. While the data for $f\in[0.99,1.0]$ are shown for random sampling,
the data for $f\in[0.99,0.991)$ are shown for the evolutionary simulations.

As a characteristic of the global structure, we counted the number of paths
$n_{path}$ connecting the input and the output. Fig. 8 shows the distribution
of paths starting from the input node and reaching the output node without
passing the same node twice. Two distributions exhibit a distinct difference;
while the probability of GRNs with only one path reaches 4% for random
sampling, it is 0.1% for evolution, and evolutionary obtained GRNs with the
number of paths less than about 100 are fewer than those of random sampling.
It suggests that the global structures of GRNs obtained by these two methods
differ significantly. However, the difference in the path distribution does
not explain all. Fig. 9 is a scatter plot of the number of paths and the
number of the lethal edges. We took data for $f\in[0.99,0.991)$ from the
evolutionary obtained GRNs. The figure clearly shows that the number of the
lethal edges is less for evolution irrespective of the number of paths. It
should be noted that even for $n_{path}=1$, the number of the lethal edge
distributes broadly. It means that the locations of the lethal edges are not
limited to the “on-path” between the input node and the output node.

##  IV Summary and discussions

We studied a toy model of gene regulatory networks and showed that evolution
enhances mutational robustness and delays the emergence of a new phenotype.
Both consequences are intuitively natural, and we were able to show them
quantitatively by comparing the evolutionary simulation with randomly sampled
GRNs as a reference set. The multicanonical ensemble Monte Carlo method was
effective in realizing random sampling of GRNs in the whole range of fitness.

The answer to the central question of this study, “does the mutational
robustness evolve?”, is “yes”. In our model, evolution progresses in two
stages; in the early stage, the measure of mutational robustness took the same
values as random sampling. At the later stage, in contrast, mutational
robustness is markedly enhanced as the fitness increases, compared to the
randomly sampled set.

The mechanism for the enhancement of mutational robustness is intuitively
understandable. The mutation we considered consists of two successive
procedures: First, an edge selected randomly is deleted. Next, a new edge is
added between a randomly selected node pair. As was shown in Ref. Nagata and
Kikuchi (2020), the edges start to be divided roughly into two types at an
intermediate value of fitness: neutral ones and ones that the fitness
decreases substantially when deleted. Here we call the latter type “lethal”,
although fitness does not necessarily drop close to zero unless fitness before
the cut is large enough. If the deleted edge is lethal, fitness drops largely,
and the possibility for fitness to recover by a randomly added new edge is
very low. Therefore, the more lethal edges a GRN has, the harder its copy
survives. In contrast, when a neutral edge is deleted, whether or not fitness
increases is determined randomly. In this scenario, the mutational robustness
evolves because GRNs containing many lethal edges are eliminated easily from
the population. This mechanism is called “the second-order selection”
according to Wagner.Wagner (2005) Details of the selection process require
further investigation.

The new phenotype, namely the GRNs exhibiting bistability, inevitably emerge
as fitness increases. It is an example that fitness constrains phenotypes. We
showed that even under such a situation, evolution behaves conservatively; the
appearance of bistability is delayed significantly compared to random
sampling. In other words, evolution proceeds so as not to allow easily the
change in phenotype. As a result, the phenotypic change occurs rapidly at the
last stage of evolution. This phenomenon may explain the punctuated
equilibrium.

The emergence of the cooperative response of GRNs consisting of many genes
having non-cooperative response function has been pointed out based on a
numerical study.Inoue and Kaneko (2018) Our result is consistent with that
finding. From the experimental side, an example has been reported that a motif
with non-cooperative response gives rise to the cooperative bistability when
built in a larger GRN, although the mechanism of cooperativity is different
from the present case.Tan  _et al._ (2009) In contrast, when the response of
each gene is ultrasensitive, small motifs can be enough to exhibit
bistability. The global structure of the network is expected to be different
largely in such cases, and thus we consider that the present results apply
mainly to large GRNs.

Next, we discuss network motifs. Some characteristic motifs are found for
highly-fit GRNs, although not prominent compared to the model in Ref Nagata
and Kikuchi (2020); the reason for this difference is attributed to the
difference in the allowed network structure. Among these motifs, coherent
feedforward loops are a frequently observed motif in real GRNs.Shen-Orr  _et
al._ (2002); Mangan and Alon (2003); Alon (2020) The following structures are
also identified as motifs, which are known to exhibit bistability if the
response of each gene is ultrasensitive: Ferrell and Machleder (1998); Ferrell
and Xiong (2001); Ferrell (2002); Angeli  _et al._ (2003); Pfeuty and Kaneko
(2009); Tiwari and Igoshin (2012); Ferrell and Ha (2014a, b) Auto-activation,
mutual activation, mutual repression, and mutual activation/repression along
with auto-activation. The last one is observed widely in multistable
GRNs.Ptashne and Gann (2002); Ptashne (2004); Huang  _et al._ (2007); Burda
_et al._ (2011); Jia  _et al._ (2017); Huang  _et al._ (2017) Since the
response function in our model is not ultrasensitive, the bistability is
realized as a cooperative phenomenon of many genes. Nevertheless, these motifs
are observed in our model, though not indispensable. They, however, are not
relevant to mutational robustness. It implies that mutational robustness is
related to the global structures of GRNs. We found that the number of paths
connecting the input node and the output node differs between randomly sampled
GRNs and evolutionary obtained ones. Although it is intuitively understandable
that GRNs having many such paths are robust because of redundancy, this
picture does not explain the origin of mutational robustness fully.

Evolutionary speed is roughly determined by the number of available GRNs or,
in other words, “entropy”. It certainly is a consequence of the fact that we
defined 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 such a simple
function in reality. However, some evolutionary experiments correspond to the
situation discussed in this paper. For example, Sato et al. reported such an
experimental study.Sato  _et al._ (2003) They put DNA coding the green
fluorescent protein (GFP) along with a random sequence to the genome of E.Coli
and conducted an artificial evolution with the selection according to the
fluorescence intensity. In this experiment, an observable single-valued
fitness was set, which certainly has the maximum value. They found that the
evolution speed slows down as the fluorescence intensity increases. It implies
a decrease in the available genotype as fitness increases. They discussed that
phenotypic divergence is relevant to evolutionary speed. Although we consider
that the variation of evolutionary speed in their experiment is explained
almost fully by the entropic effect, the relationship between phenotypic
divergence and genotypic entropy requires further exploration. In the natural
situation departed from the experimental condition, evolution is not
restricted to proceed in one direction; when evolution to one direction
becomes difficult, it will change direction. We consider that the concept of
entropic effect on evolutionary speed can apply also to such a natural
situation.

Finally, we should be careful about the model dependence. We believe that the
three main results, that is, the evolution of mutational robustness, delay in
the emergence of a new phenotype, and relationship between evolutionary speed
and genetic entropy, apply universally to other systems. Similar studies for
different models and other objects are necessary to confirm.

###### Acknowledgements.

The authors thank Koich Fujimoto, Katsuyoshi Matsushita, and Hajime Yoshino
for fruitful discussions and comments.

## References

  * Kitano (2004) H. Kitano, Nat Rev Genet 5, 828 (2004).
  * Wagner (2005) A. Wagner,  _Robustness and Evolvability in Living Systems_ (Princeton Univ Press, 2005).
  * Masel and Siegal (2009) J. Masel and M. L. Siegal, Trends Genet 25, 395 (2009).
  * Giaever  _et al._ (2002) G. Giaever, A. M. Chu, L. Ni, C. Connelly, L. Riles, S. Veronneau, S. Dow, A. Lucau-Danila, K. Anderson, B. Andre, A. P. Arkin, A. Astromoff, M. E. Bakkoury, R. Bangham, R. Benito, S. Brachat, S. Campanaro, M. Curtiss, K. Davis, A. Deutschbauer, K. D. Entian, P. Flaherty, F. Foury, D. J. Garfinkel, M. Gerstein, D. Gotte, U. Gludener, J. H. Hegemann, S. Hempel, Z. Herman, D. F. Jaramillo, D. E. Kelly, S. L. Kelly, P. Kotter, D. LaBonte, D. C. Lamb, N. Lan, H. Liang, H. Liao, L. Liu, C. Luo, M. Lussier, R. Mao, P. Menard, S. L. Ooi, J. L. Revuelta, C. J. Roberts, M. Rose, P. Ross-Macdonald, B. Scherens, G. Schimmack, B. Shafer, D. D. Shoemaker, S. Sookhai-Mahadeo, R. K. Storms, J. N. Strathern, G. Valle1, M. Voet, G. Volckaert, C. Y. Wang, T. R. Ward, J. Wilhelmy, E. A. Winzeler, Y. Yang, G. Yen, E. Youngman, K. Yu, H. Bussey, J. D. Boeke, M. Snyder, P. Philippsen, R. W. Davis, and M. Johnston, Nature 418, 387 (2002).
  * Baba  _et al._ (2006) T. Baba, T. Ara, M. Hasegawa, Y. Takai, Y. Okumura, M. Baba, K. A. Datsenko, M. Tomita, B. L. Wanner, and H. Mori, Mol Syst Biol 2, 2006.0008 (2006).
  * Kim  _et al._ (2010) D. U. Kim, J. Hayles, D. Kim, V. Wood, H. O. Park, M. Won, H. S. Yoo, T. Duhig, M. Nam, G. Palmer, S. Han, L. Jeffery, S. T. Baek, H. Lee, Y. S. Shim, M. Lee, L. K. and. K. S. Heo, E. J. Noh, A. R. Lee, Y. J. Jang, K. S. Chung, S. J. Choi, J. Y. Park, Y. Park, H. M. Kim, S. K. Park, H. J. Park, E. J. Kang, H. B. Kim, H. S. Kang, H. M. Park, K. Kim, K. Song, K. B. Song, P. Nurse, and K. L. Hoe, Nat Biotechnol 28, 617 (2010).
  * Isaran  _et al._ (2008) M. Isaran, C. Lemerle, K. Michalodimitrakis, C. Horn, P. Beltrao, E. Raineri, M. Garriga-Canut, and L. Serrano, Nature 452, 840 (2008).
  * Kauffman (1969) S. A. Kauffman, J Theor Biol 22, 437 (1969).
  * Kauffman and Levin (1987) S. A. Kauffman and S. Levin, J Theor Biol 128, 11 (1987).
  * Kauffman (1993) S. A. Kauffman,  _The Origins of Order: Self-Organization and Selection in Evolution_ (Oxford Univ Press, 1993).
  * Ciliberti  _et al._ (2007a) S. Ciliberti, O. C. Martin, and A. Wagner, PLoS Comput Biol 3, e15 (2007a).
  * Ciliberti  _et al._ (2007b) S. Ciliberti, O. C. Martin, and A. Wagner, Proc Natl Acad Sci USA 104, 13591 (2007b).
  * Burda  _et al._ (2011) Z. Burda, A. Krzywicki, O. C. Martin, and M. Zagorski, Proc Natl Acad Sci USA 108, 17263 (2011).
  * Zagorski  _et al._ (2013) M. Zagorski, A. Krzywicki, and O. C. Martin, Phys Rev E 87, 012727 (2013).
  * Saito and Kikuchi (2013) N. Saito and M. Kikuchi, New J Phys 15, 053037 (2013).
  * Berg and Neuhaus (1991) B. A. Berg and T. Neuhaus, Phys Lett B267, 249– (1991).
  * Berg and Neuhaus (1992) B. A. Berg and T. Neuhaus, Phys Rev Lett 68, 9 (1992).
  * Saito  _et al._ (2010) N. Saito, Y. Iba, and K. Hukushima, Phys Rev E 82, 031142 (2010).
  * Saito and Iba (2011) N. Saito and Y. Iba, Comput Phys Commun 182, 223–225 (2011).
  * Kitajima and Iba (2011) A. Kitajima and Y. Iba, Comput Phys Commun 182, 251 (2011).
  * Iba  _et al._ (2014) Y. Iba, N. Saito, and A. Kitajima, Ann Inst Stat Math 66, 611 (2014).
  * Kitajima and Kikuchi (2015) A. Kitajima and M. Kikuchi, PLoS ONE 10, e0 125062 (2015).
  * Nagata and Kikuchi (2020) S. Nagata and M. Kikuchi, PLoS Comput Biol 16, e1007969 (2020).
  * Goldbeter and Koshland (1981) A. Goldbeter and D. E. Koshland, Proc Natl Acad Sci USA 78, 6840 (1981).
  * Ferrell and Ha (2014a) J. E. Ferrell and S. H. Ha, Trends Biochem Sci 39, 496 (2014a).
  * Alon (2020) U. Alon,  _An Introduction to Systems Biology (second edition)_ (CRC Press, 2020).
  * Mjolsness  _et al._ (1991) E. Mjolsness, D. H. Sharp, and J. Reinitz, J Theor Biol 152, 429 (1991).
  * Wagner (1994) A. Wagner, Proc Natl Acad Sci USA 91, 4387 (1994).
  * Wagner (1996) A. Wagner, Evolution 50, 1008 (1996).
  * Siegal and Bergman (2002) M. L. Siegal and A. Bergman, Proc Natl Acad Sci USA 99, 10528 (2002).
  * Masel (2004) J. Masel, J Evol Biol 17, 1106 (2004).
  * Kaneko (2007) K. Kaneko, PLoS one 2, e434 (2007).
  * Espinosa-Soto  _et al._ (2011) C. Espinosa-Soto, O. C. Martin, and A. Wagner, J Evol Biol 24, 1284 (2011).
  * Inoue and Kaneko (2013) M. Inoue and K. Kaneko, PLoS Comput Biol 9, e1003001 (2013).
  * Inoue and Kaneko (2018) M. Inoue and K. Kaneko, Eutophys Lett 124, 38002 (2018).
  * Furusawa and Kaneko (2008) C. Furusawa and K. Kaneko, PLoS Comput Biol 4, e3 (2008).
  * Pujato  _et al._ (2013) M. Pujato, T. MacCarthy, A. Fiser, and A. Bergman, PLoS Comput Biol 9, e1002865 (2013).
  * Hopfield and Tank (1985) J. J. Hopfield and D. W. Tank, Biol Cybern 52, 141 (1985).
  * Lee (1993) J. Lee, Phys Rev Lett 71, 211–214 (1993).
  * Wang and Landau (2001a) F. Wang and D. P. Landau, Phys Rev Lett 86, 2050– (2001a).
  * Wang and Landau (2001b) F. Wang and D. P. Landau, Phys Rev E 64, 056101 (2001b).
  * Tyson  _et al._ (2003) J. J. Tyson, K. C. Chen, and B. Novak, Curr Opin Cell Biol 15, 221 (2003).
  * Ptashne and Gann (2002) M. Ptashne and A. Gann,  _Genes & Signals_ (Cold Spring Harbor Laboratory Press, 2002).
  * Ptashne (2004) M. Ptashne,  _A Genetic Switch (third edition)_ (Cold Spring Harbor Laboratory Press, 2004).
  * Zong  _et al._ (2010) C. Zong, L. h. So, L. A. Sepuiveda, S. O. Skinner, and I. Golding, Mol Syst Biol 6, 440 (2010).
  * Watson  _et al._ (2013) J. D. Watson, T. A. Baker, S. P. Bell, A. Gann, M. Levine, and R. Losick,  _Molecular Biology of the Gene (7th edition)_ (Pearson Education, 2013).
  * Ferrell and Machleder (1998) J. E. Ferrell and E. M. Machleder, Science 280, 895 (1998).
  * Huang  _et al._ (2007) S. Huang, Y.-P. Guo, G. May, and T. Enver, Dev Biol 305, 695 (2007).
  * Jia  _et al._ (2017) D. Jia, M. K. Jolly, W. Harrison, M. Boareto, E. Ben-Jacob, and H. Levine, Phys Biol 14, 035007 (2017).
  * Huang  _et al._ (2017) B. Huang, M. Lu, D. Jia, E. Ben-Jacob, H. Levine, and J. N. Onuchic, PLoS Comput Biol 13, e1005456 (2017).
  * Shen-Orr  _et al._ (2002) S. S. Shen-Orr, R. Milo, S. Mangan, and U. Alon, Nat Genet 31, 64 (2002).
  * Mangan and Alon (2003) S. Mangan and U. Alon, Proc Natl Acad Sci USA 100, 11980 (2003).
  * Ferrell and Xiong (2001) J. E. Ferrell and W. Xiong, Chaos 11, 227 (2001).
  * Ferrell (2002) J. E. Ferrell, Curr Open Chem Biol 6, 140 (2002).
  * Angeli  _et al._ (2003) D. Angeli, J. E. Ferrell, and E. D. Sontag, Proc Natl Acad Sci USA 101, 1822 (2003).
  * Pfeuty and Kaneko (2009) B. Pfeuty and K. Kaneko, Phys Biol 6, 046013 (2009).
  * Tiwari and Igoshin (2012) A. Tiwari and O. A. Igoshin, Phys Biol 9, 055003 (2012).
  * Ferrell and Ha (2014b) J. E. Ferrell and S. H. Ha, Trends Biochem Sci 39, 612 (2014b).
  * Tan  _et al._ (2009) C. Tan, P. Marguet, and L. You, Nat Chem Biol 5, 842 (2009).
  * Sato  _et al._ (2003) K. Sato, Y. Ito, T. Yomo, and K. Kaneko, Proc Natl Acad Sci USA 100, 14086 (2003).