# Evolutionary Shaping of Low-Dimensional Path Facilitates Robust and Plastic
Switching Between Phenotypes

Ayaka Sakata<EMAIL_ADDRESS>Department of Statistical Inference &
Mathematics, The Institute of Statistical Mathematics, Tachikawa, Tokyo
190-8562, Japan  Department of Statistical Science, The Graduate University
for Advanced Science (SOKENDAI), Hayama-cho, Kanagawa 240-0193, Japan
Kunihiko Kaneko  Center for Complex Systems Biology, Universal Biology
Institute, University of Tokyo, Komaba, Tokyo 153-8902, Japan  The Niels Bohr
Institute, University of Copenhagen, Blegdamsvej 17, Copenhagen 2100-DK,
Denmark

###### Abstract

Biological systems must be robust for stable function against perturbations,
but robustness alone is not sufficient. The ability to switch between
appropriate states (phenotypes) in response to different conditions is
essential for biological functions, as observed in allosteric enzymes and
motor proteins. How are robustness and plasticity simultaneously acquired
through evolution? In this study, we examine the evolution of genotypes that
realize plastic switching between two endpoint phenotypes upon external
inputs, as well as stationary expressions of phenotypes. Here, we introduce a
statistical physics model consisting of spins, with active sites and
regulatory sites, which are distinct from each other. In our model, we
represent the phenotype and genotype as spin configurations and the spin-spin
interactions, respectively. The fitness for selection is given by the spin
configuration, whose behavior is governed by the genotypes. Specifically, the
fitness for selection is given so that it takes a higher value as more of the
active sites take two requested spin configurations depending on the states of
the regulatory sites. The remaining spins do not directly affect the fitness,
but they interact with other spins. We numerically evolve the matrices of
spin-spin interactions (genotypes) by changing them with mutations and
selection of those with higher fitness. Our numerical simulations show that
characteristic genotypes with higher fitness evolve slightly above the phase
transition temperature between replica symmetric and replica symmetry breaking
phase in spin-glass theory. These genotypes shape two spin configurations
separately depending on the regulation, where the two phenotypes are
dominantly represented by the genotypes’ first and second eigenmodes, where
smooth switching of two phenotypes are achieved by following a one-
dimensional, quarter-circle path connecting the two phenotypes. Upon changes
in regulations, spin configurations are attracted to this path, which allows
for robust and plastic switching between the two phenotypes. The statistical-
physics analysis based on the two eigenmodes show that the free energy
landscape has a valley along the quarter-circle one-dimensional switching
path. Robust attraction to the path is achieved through the evolution of
interaction within non-active and non-regulatory spin sites, which themselves
do not contribute to fitness. Our finding indicates that the compatibility of
the robustness and plasticity is acquired by the evolution of the low-
dimensionality in the phenotype space, which will be relevant to the
understanding of the robust function of protein as well as the material
design.

##  I Introduction

Biological systems are inherently complex, comprising numerous elements.
Despite such complexity, they function robustly under environmental and
stochastic perturbations. The biological function, in general, is given as a
result of phenotypes, which are generated via dynamics based on genetic
information. As a consequence, the function-related phenotypes need to be
robustly shaped through the dynamics. However, a single robust phenotype or
fitted state is insufficient for a biological system to function under varying
conditions. Phenotypes must exhibit plasticity, shifting to appropriate
patterns in response to relevant signals or inputs [1, 2]. For instance, the
active sites of enzyme proteins can change between two conformations known as
tense and relaxed states, induced by allosteric regulation [3, 4, 5, 6]. Motor
proteins, such as the myosin, kinesin, and dynein families, exhibit large-
scale conformational changes in response to binding events [7, 8].
Phosphorylation of substrates in the mitogen-activated protein kinase cascades
can switch between two states depending on modification by phosphatase or
diphosphatase [9]. Gene expression pattern switches in response to signals are
also necessary for cell survival. Thus, the ability to switch between
appropriate phenotypes in response to different conditions is essential for
biological functions. Accordingly, the presence of multiple phenotypes and
transitions among them in response to inputs must be shaped through evolution.
Considering such changes in phenotypes, then, plasticity to external
conditions also needs to be required for the switching to different
phenotypes, in addition to the robust expression of each phenotypes. In
general, how robustness and plasticity are compatible remains a fundamental
question in biology [10, 11].

To study such plastic responses in biological systems, it is essential to
understand the nature of switching pathways, in addition to the multiple
phenotypes corresponding to endpoint structures. An understanding of the
switching pathways can aid in development of engineering techniques, such as
drug design, that target the intermediate states of the switching pathways
[12, 13]. Despite advances in structural biology in recent decades, molecular-
level characterization of switching remains a challenge due to limitations in
macromolecular X-ray crystallography, nuclear magnetic resonance, and small-
angle X-ray scattering [14]. Hence, theoretical or numerical approaches are
necessary to understand general characteristics of large-scale conformational
switching [15, 16]. For instance, the plastic network model, an extension of
the elastic network model [17, 18, 19], was utilized to generate
conformational switching pathways that are consistent with experimental data
of the intermediate structures in Escherichia coli adenylate kinase [7]. The
resulting pathways resemble combinations of low-energy normal modes obtained
for the endpoint structures [7]. It has then been suggested that such
preferred directionality may contribute to catalysis in many enzymes,
achieving extraordinary rate acceleration and specificity [20]. For Src
kinase, such switching paths were explored by using a coarse-grained, two-
state Go model, characterized by a two-dimensional free energy landscape [21].

In general, theoretical and numerical methods to explore conformational
changes assume the existence of probable switching paths, which minimize
energy, free energy, or action [22]. The existence of a probable path implies
that possible transient changes are constrained along the path. Further, low-
dimensional approximations using principle component analysis have often been
adopted to simplify the numerically or experimentally obtained switching paths
[16]. These studies suggest the importance of understanding how low-
dimensional switching paths are shaped and evolved in the phenotypic spaces.

As for the stationary states, recent experimental and numerical observations
have shown that evolved phenotypes are often constrained within a low-
dimensional manifold despite the high dimensionality of the phenotype space.
For example, changes in (logarithmic) concentrations of mRNAs or proteins have
been found to be correlated [23, 24, 25, 26] or proportional [27, 28] across
all components under various environmental stresses. Numerical simulations of
cell models with catalytic reaction networks have also demonstrated that
evolved phenotypic changes caused by environmental and mutational changes are
constrained within a low-dimensional manifold [29]. This reduction in
dimensionality from high-dimensional phenotypes has also been observed in the
structural changes of proteins, as a result of data analysis [30].
Additionally, such dimensional reduction is suggested to be a result of the
robustness of phenotypes shaped by evolution. However, such studies are
limited to phenotypes around the endpoint structures, i.e., the stationary
conditions. In this study, we examine the evolution of the switching path from
the viewpoint of dimensional reduction.

In particular, we address the following questions:

  * •

Under what conditions and how are multiple endpoint phenotypes shaped
depending on external inputs and stabilized through evolution?

  * •

Are low-dimensional constraints of switching paths shaped through evolution?

  * •

What are the characteristics of switching paths between endpoints?

  * •

What are the characteristics of evolved genotypes that allow robust switching
paths?

To address these questions, we extend a spin-statistical physics model
introduced previously [31]. In this model, the spin variables $\bm{S}$ and
their interaction variables $\bm{J}$ represent phenotype and genotype [32, 33,
34], respectively, and fitness is provided by certain spin configurations. We
consider two endpoint structures, corresponding to those under regulation and
without regulation. We introduce active and regulatory sites in the spin
system to represent the effect of external regulation applied to the
regulatory sites. The fitness of selective evolution depends on the
appropriate expression of configurations. Fitted interactions can provide two
configurations of active spins, corresponding to regulated and non-regulated
cases.

Numerical evolution allows us to examine how the robustness of each phenotype,
as well as its plasticity to switch between the two configurations, is shaped
by regulation. Our result shows that, as a result of evolution, the
dimensional reduction to a two-dimensional phenotype space appears, under a
certain range of temperatures, while a one-dimensional path is shaped for the
switch between the two phenotypes in the regulated and non-regulated cases.
The shaped switching path is robust to thermal noise and genetic mutation. In
terms of statistical physics, the robustness of the fitted phenotype is
achieved in the replica-symmetric phase. In contrast, the plasticity of the
switch increases as the temperature approaches the replica-symmetry breaking
(RSB) transition. We then will show that robust response is achieved near the
RSB transition.

##  II Model

Figure 1: (a) Schematic representation of the conformational change induced by
regulation. A, R, and S denote active sites, regulatory sites, and substrates,
respectively. The molecule denoted by L is a ligand that regulates the protein
through regulatory sites. (b) The spin model for conformational switching that
has active and regulatory sites. (c) A landscape picture of the conformational
changes induced by regulation discussed in this study. (d) Free energy
landscape of the multi-pattern embedding in the associative memory model such
as Hopfield networks.

Here, we introduce an abstract model of interacting spins, as a simplified
representation of proteins whose active sites conformation is regulated by
regulatory sites. Fig. 1 (a) gives a simplified picture of the regulation and
related conformational changes adopted in this study. The protein shown as
grey in the figure has an active (‘A’) and a regulatory site (‘R’), which are
parts of the protein, consisting of amino-acid residues. In general allosteric
regulation, the active and regulatory sites are located sufficiently apart and
do not interact directly. As shown in Fig. 1 (a), binding of the ligand to the
regulatory site leads to the conformational change in the active site, via
interaction with sites other than the active sites and regulatory sites. In
contrast, without the binding of the activator to the regulatory site, such
conformational changes in the active site do not occur and remain in its
original conformation.

We introduce an abstract statistical-physics model with interacting spins
representing conformation, as shown in Fig.1(b). The model consists of spin
variable $\bm{S}=\\{S_{1},\cdots,S_{N}\\}\in\\{-1,+1\\}^{N}$ which represents
the conformational change in each amino-acid residue, whereas the coupling
$\bm{J}$ between spins represents the interaction among residues. These are
respectively denoted by nodes and edges in Fig. 1(b). Here, we set $\bm{J}$ as
$N\times N$ symmetric matrix. Its elements are given by
$J_{ii}=0~{}(i=1,\cdots,N)$, and $J_{ij}\in\Omega_{J}$ for $i\neq j$ where
$\Omega_{J}=\\{-1/\sqrt{N},0,1/\sqrt{N}\\}$. Active (‘A’) and regulatory (‘R’)
sites are represented by $N_{A}$ and $N_{R}$ spins among the $N$ spins. The
label sets of active and regulatory spins are denoted by ${\cal A}$ and ${\cal
R}$, respectively. We set $J_{ij}=0$ for $i\in{\cal A}$ and $j\in{\cal R}$ or
$i\in{\cal R}$ and $j\in{\cal A}$, to prohibit the direct interaction between
regulatory and active sites. The spin variables other than those at the active
and regulatory sites are called free sites, as shown in Fig. 1 (b).

For the dynamics of spin-variables under given $\bm{J}$, we adopt the
transition rule of spins from $\bm{S}$ to $\bm{S}^{\prime}$ under given
Hamiltonian with the interaction matrix $\bm{J}$ and temperature $T$ as

$\displaystyle\mathrm{Pr}[\bm{S}\to\bm{S}^{\prime}|\bm{J}]=\min\\{\exp(-\beta\Delta_{H}(\bm{S},\bm{S}^{\prime}|\bm{J})),1\\},$
(1)

where $\beta=T^{-1}$ is the inverse of temperature, and
$\Delta_{H}(\bm{S},\bm{S}^{\prime}|\bm{J})\equiv
H(\bm{S}^{\prime}|\bm{J})-H(\bm{S}|\bm{J})$. Here, we set the Hamiltonian as

$\displaystyle H(\bm{S}|\bm{J})=-\frac{1}{2}\bm{S}^{\top}\bm{J}\bm{S},$ (2)

where $\top$ denotes the matrix transpose.

Note that in this statistical-physics model, we adopt the spin variables
$\\{-1,1\\}$, instead of continuous conformational variables in residues. This
is a highly simplified model by nature (see [35] for examples of spin models
for protein dynamics). Here, we aim to elucidate how certain stochastic
dynamics for generating functional phenotypes are shaped through evolution. To
this end, the present model capture the essence of such dynamics and genotype-
phenotype mapping, in which spin variables $\bm{S}$ corresponding to the
phenotypes are shaped by high-dimensional dynamics under genetic rules given
by the interaction matrix $\bm{J}$, whereas regulation is referred to as
change in a part of “regulatory” spins, as defined below.

Next, the functional change in the active sites is postulated by the
appropriate change in the configuration of active spins $\bm{S}_{\cal
A}=\\{S_{i}|i\in{\cal A}\\}$, depending on the configurations of the
regulatory site $\bm{S}_{\cal R}=\\{S_{i}|i\in{\cal R}\\}$. Here, instead of
introducing the binding of ligands to regulatory sites as external variables,
we assume that the configuration of the spins is set at $\bm{S}_{\cal R}^{+}$
upon the binding. That is, among $2^{N_{R}}$ possible configurations of the
regulatory spins, the regulatory spins only take the configuration in
$\bm{S}_{\cal R}^{+}$ when the ligand binding occurs. Further, we consider
that $\bm{S}_{\cal R}^{+}$ cannot appear without the ligand binding.
Accordingly, the equilibrium distribution upon the regulation and non-
regulation is given by

$\displaystyle P_{\beta}^{+}(\bm{S}|\bm{J})=\frac{1}{Z_{\beta}^{+}}\exp(-\beta
H(\bm{S}|\bm{J})),~{}~{}Z_{\beta}^{+}(\bm{J})=\sum_{\bm{S}|\bm{S}_{\cal
R}\in\bm{S}_{\cal R}^{+}}\exp(-\beta H),$ (3) $\displaystyle
P_{\beta}^{-}(\bm{S}|\bm{J})=\frac{1}{Z_{\beta}^{-}}\exp(-\beta
H(\bm{S}|\bm{J})),~{}~{}~{}Z_{\beta}^{-}(\bm{J})=\sum_{\bm{S}|\bm{S}_{\cal
R}\notin\bm{S}_{\cal R}^{+}}\exp(-\beta H),$ (4)

where $\bm{S}|\bm{S}_{\cal R}\in\bm{S}_{\cal R}^{+}$ and $\bm{S}|\bm{S}_{\cal
R}\notin\bm{S}_{\cal R}^{+}$ indicate the set of possible configurations for
regulated and non-regulated states, respectively.

Next, the functional change in configurations of the active spins in response
to the regulation is given by the change in regulatory spins from
$\bm{S}_{\cal A}^{-}$ to $\bm{S}_{\cal A}^{+}$: Thus, the conformational
change induced by regulation is modeled as follows: if the configuration of
regulatory spins is set at $\bm{S}_{\cal R}^{+}$, the configuration of the
active spins turns into $\bm{S}_{\cal A}^{+}$; otherwise, the configuration of
the active sites stays at $\bm{S}_{\cal A}^{-}$.

The function of the present system to express the target spin pattern
$\bm{S}_{\cal A}^{\pm}$ appropriately can be measured by the magnetization
$m_{A}^{\pm}$ defined as the overlap of the spins in the active sites with the
corresponding target spin patterns as

$\displaystyle m_{A}^{+}$ $\displaystyle=\frac{1}{N_{A}}\sum_{i\in{\cal
A}}S_{i}S_{i}^{+}$ (5) $\displaystyle m_{A}^{-}$
$\displaystyle=\frac{1}{N_{A}}\sum_{i\in{\cal A}}S_{i}S_{i}^{-}.$ (6)

Finally, the overall fitness that measures the functionality of the present
system is given by the sum of the expectations of $m_{A}^{\pm}$ as

$\displaystyle\psi(\bm{J})=\frac{1}{2}\left\\{\langle|m_{A}^{+}|\rangle_{+}+\langle|m_{A}^{-}|\rangle_{-}\right\\},$
(7)

where $\langle\cdot\rangle_{+}$ and $\langle\cdot\rangle_{-}$ are the
expectation values according to the equilibrium distributions for regulated
and non-regulated cases eq.(3) and eq.(4),

The evolution of genotypes $\bm{J}$ is then based on the above fitness
$\psi(\bm{J})$. Higher fitness genotypes are selected under given selective
pressure: At generation $g$, the evolutionary change in $\bm{J}$ to increase
the fitness is given by

$\displaystyle\mathrm{Pr}[\bm{J}^{(g)}\to\bm{J}^{(g+1)}]=\min\\{\exp(\beta_{J}\Delta\psi),1\\},$
(8)

where $\Delta\psi=\psi(\bm{J}^{(g+1)})-\psi(\bm{J}^{(g)})$. The parameter
$\beta_{J}=T_{J}^{-1}$ expresses the selection pressure, and the genotypes are
selected uniformly at high temperatures $T_{J}\to\infty$, whereas at low
$T_{J}$, genotypes with higher fitness values are preferred.

Remark: The celebrated Hopfield neural network model can be used for embedding
several patterns in spin models. In this case, as schematically shown in Fig.
1 (d), multiple patterns with different spin configurations were reached
depending on the initial condition given by Hamiltonian dynamics. In contrast,
in our case, by external inputs to regulatory spins (i.e., with inputs or with
different boundary conditions), different spin configurations are reached
depending on if the regulatory spins are regulated or not, from the same
initial conditions for the two, as in Fig.1 (c), which has been introduced in
the study of the reshaping of the energy landscape of a protein by allostery
[36]. (In the context of the neural network model, this corresponds to the
associative memory model upon external inputs [37].)

##  III Numerical simulation

Without loss of generality, we set the indices of the regulatory sites and
active sites as ${\cal R}=\\{N-N_{R}+1,\ldots,N\\}$ and ${\cal
A}=\\{1,\ldots,N_{\mathrm{A}}\\}$, respecively. Further, we set the
configuration $\bm{S}_{\cal R}^{+}$ as $\bm{S}_{\cal
R}^{+}=\\{\\{+1,\cdots,+1\\},\\{-1,\cdots,-1\\}\\}$. For the desirable
configurations of the active sites, we set
$\bm{S}_{A}^{+}=\\{\\{+1,\cdots,+1\\},\\{-1,\cdots,-1\\}\\}$ and
$\bm{S}_{A}^{-}=\\{\\{+1,-1,\cdots\\},\\{-1,+1,\cdots\\}\\}$. In the genotype
evolution process, we induce a 10-point mutation at each generation to
generate the candidate of the next generation $\bm{J}^{\prime}$ from $\bm{J}$,
maintaining the symmetry ${\bm{J}^{\prime}}^{\top}=\bm{J}^{\prime}$. Here, we
mainly show the results for $N=100$, $N_{A}=5$ and $N_{R}=10$, and the free
sites consist of $N-N_{A}-N_{R}=85$ spin variables. We update $\bm{J}$ at a
sufficiently large value as $\beta_{J}=100$, and discuss $T$-dependencies.

###  III.1 Fitness, rugged landscape, and separation of two patterns

Figure 2: Evolutional dynamics of $\langle m_{A}^{\pm}\rangle_{\pm}$
associated with the evolution of a genotype $\bm{J}$. An example at $T=0.91$
is shown in (a), and examples at $T=0.67$ are shown in (b)-(d).

In Fig.2, we show examples of evolutionary dynamics of $\bm{J}$ through the
evolutionary changes of $\langle|m_{A}^{\pm}|\rangle_{\pm}^{(g)}$, where
$\langle\cdot\rangle^{(g)}$ denotes the expectation according to the
distribution $P_{\beta}^{\pm}(\bm{S}|\bm{J}^{(g)})$ with a $g$-th generation
genotype $\bm{J}^{(g)}$. These quantities measure the tendency to exhibit
desirable patterns depending on the regulatory site, whereas fitness is given
by their mean, as eq.(7). Fig.2 (a) shows an example of the generational
changes of $\langle|m_{A}^{\pm}|\rangle_{\pm}^{(g)}$ at $T=0.91$. They show a
negative correlation; when $\langle|m_{A}^{+}|\rangle_{+}^{(g)}$ increases,
$\langle|m_{A}^{-}|\rangle_{-}^{(g)}$ decreases, and vice versa. For the
genotypes that show this behavior, which are the most evolved genotypes around
$T=0.91$, the simultaneous expression of $\bm{S}_{A}^{+}$ and
$\bm{S}_{A}^{-}$, depending on the regulatory sites, is difficult. When the
active sites take one of the configurations in $\bm{S}_{A}^{+}$ or
$\bm{S}_{A}^{-}$, irrespective of the regulatory sites, we obtain
$|m_{A}^{+}|=1$ and $|m_{A}^{-}|=0.2$ or $|m_{A}^{+}|=0.2$ and
$|m_{A}^{-}|=1$, respectively. Therefore, the fitness value of the genotype
that can express only one desirable pattern among $\bm{S}_{A}^{\pm}$ can reach
a highest value of 0.6. Meanwhile, at a lower temperature $T=0.67$, both
$\langle|m_{A}^{\pm}|\rangle_{\pm}^{(g)}$ increase simultaneously after the
evolution, as shown in Fig.2(b)-(d), where the fitness reach around 0.9. There
are three evolutionary courses; $\langle|m_{A}^{+}|\rangle_{+}^{(g)}$ or
$\langle|m_{A}^{-}|\rangle_{-}^{(g)}$ increases first (Fig.2(b) or (c)), or
they increase simultaneously (Fig.2(d)). Among 100 samples, 21, 33, and 46
samples follow each course, respectively.

Figure 3: $T$-dependence of (a) fitness, (b) the fraction of genotypes on
which the BP algorithm does not converge, (c) $\langle m_{A}^{+}\rangle_{+}$
and $\langle m_{A}^{-}\rangle_{-}$, and (d) Similarity between regulated and
non-regulated states. Each data point is averaged over 100 samples of evolved
$\bm{J}$. The vertical dotted line, two-dot chain line, and one-dot chain line
denote $T_{0}$, $T_{1}$, and $T_{2}$, respectively.

For each $T$, we obtain 100 samples of the evolved $\bm{J}$ updated for
$g=10^{5}$ generations, and denote the temperature-dependence ensemble of
evolved genotypes as ${\cal J}(T)$. Fig.3 (a) shows the $T$-dependence of the
mean of the fitness among ${\cal J}(T)$. With a decrease in $T$, the fitness
value $\Psi$ increases from 0.4, which is a trivial value given by the uniform
distribution of the phenotype $\bm{S}$. $T_{0}$ is defined as the transition
temperature characterized by fitness below, in which the fitness value
increases as $T$ decreases. To be precise, it is defined by the point where
the second derivative of $\Psi$ shows discontinuity. We term the phase
$T>T_{0}$ as the paramagnetic phase. Next, the energy landscape on
$\bm{J}\in{\cal J}(T)$ governing the phenotype expression dynamics changes at
$T=T_{2}$: We term the phases at $T_{0}>T>T_{2}$ and $T<T_{2}$ as replica-
symmetric (RS) phase and replica symmetry breaking (RSB) phase, respectively.
The difference between the two phases can be detected by the belief
propagation (BP) algorithm [38, 39, 40]. In the fully connected spin-glass
system, the stability condition of the BP algorithm agrees with the validity
of the RS assumption in the replica analysis, which is known as de Almeida-
Thouless (AT) instability [41, 42]; hence, when the BP algorithm converges,
the system on $\bm{J}$ corresponds to the RS phase, otherwise the RSB phase.
The RSB indicates the rugged landscape with exponential orders of metastable
states, and the phenotype expression dynamics is not robust to thermal
fluctuation [32]111 While the AT instability cannot be derived from not fully
connected systems, including our evolutionary model, the BP algorithm is often
used as a numerical method to identify the phase of not fully connected
systems.. At $T>T_{2}$, most of the evolved genotypes in ${\cal J}(T)$ exhibit
rapid convergence of the BP algorithm. Meanwhile, the BP algorithm cannot
converge for most evolved genotypes in ${\cal J}(T)$ when $T<T_{2}$. In
Fig.3(b), we present the fraction of evolved genotypes for which the BP
algorithm does not appear to converge within $10^{5}$ steps, which increases
as $T$ is decreased below $T_{2}$.

The existence of these transitions from the paramagnetic phase to the RS
phase, and then to the RSB phase, is common with the evolving spin-glass model
to express one specific phenotype [32, 33, 34, 31]. In the present model,
however, another transition appeared at $T_{1}$, with respect to the
achievability of two patterns. In Fig.3(c), we show the temperature dependence
of the overlaps $\langle|m_{A}^{+}|\rangle_{+}$ and
$\langle|m_{A}^{-}|\rangle_{-}$, whose mean corresponds to fitness. At
$T_{1}<T<T_{0}$, $\langle|m_{A}^{+}|\rangle_{+}$ contributes more to fitness,
and only the phenotype expression with regulation is preferentially shaped. At
$T<T_{1}$, both the increase of $\langle|m_{A}^{-}|\rangle_{+}$ and
$\langle|m_{A}^{-}|\rangle_{-}$ are achieved depending on the regulatory site.
We term the phases $T_{0}>T>T_{1}$ and $T_{1}>T>T_{2}$ as RS1 and RS2,
respectively. Here, we note a negative correlation between $\langle
m_{A}^{\pm}\rangle_{\pm}^{(g)}$ is observed in the paramagnetic and RS1
phases, as shown in Fig.2(a). In the RS2 and RSB phases, the increase of both
$\langle m_{A}^{\pm}\rangle_{\pm}^{(g)}$ is achieved after a sufficient
update, as shown in Fig.2(b)-(d).

To study the transition at $T=T_{1}$, we examine the probability distributions
of spin configurations, $p_{\beta}^{+}(\bm{S}|\bm{J})$ and
$p_{\beta}^{-}(\bm{S}|\bm{J})$, with and without regulations, by means of the
component-wise expected phenotype for each $i=1,\ldots,N$ defined as

$\displaystyle\mu_{i}^{\pm}=\left\langle\mathrm{sign}\left(\sum_{i=1}^{N_{A}}S_{i}\right)S_{i}\right\rangle_{\pm},$
(9)

where the term $\mathrm{sign}(\sum_{i=1}^{N_{A}}S_{i})$ is introduced to break
the Z2 symmetry. In Fig.3(d), we show $T$-dependence of the similarities
between two mean phenotypes measured by
$\sum_{i=1}^{N}\mu_{i}^{+}\mu_{i}^{-}/N$. As shown in Fig.3(d), the overlap
shows a peak at $T=T_{1}$, and it decreases as $T$ decreases below $T_{1}$.
According to the decrease in the overlap, the transition between the phenotype
with and without regulation involves large conformational changes.

###  III.2 Two-dimensional structure in the phenotype space

Figure 4: Comparison between the eigenvectors ($\bm{\xi}_{1}$ and
$\bm{\xi}_{2}$) and $\bm{\mu}^{\pm}$ for one sample of evolved $\bm{J}\in{\cal
J}(T)$ at (a) $T=0.833$ and (b) $T=0.667$. The solid lines show the
relationship $y=\tanh(\beta\sqrt{N}x)$.

We investigate how the two patterns shaped by evolution are separated in the
RS2 phase. To compare $N$-dimensional mean phenotypes $\bm{\mu}^{+}$ under and
$\bm{\mu}^{-}$ without regulation, it is convenient to determine a reference
coordinate system. We adopte the eigenvectors of the evolved genotypes as the
axes to represent mean phenotypes. Using the eigenvectors and corresponding
eigenvalues, the genotype $\bm{J}$ is decomposed as

$\displaystyle\bm{J}=\sum_{r=1}^{N}\lambda_{r}\bm{\xi}_{r}\bm{\xi}_{r}^{\top},$
(10)

where $\bm{\xi}_{i}$ and $\lambda_{i}$ are $i$-th eigenvector and $i$-th
eivenvalue. We set the indices of the eigenmodes to be
$\lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{N}$.

In Fig. 4, we show the scatter plots between $\mu_{i}^{\pm}$ against
eigenvectors $\xi_{1i}$ and $\xi_{2i}$ for $i=1,\ldots,N$ at (a) $T=0.833$
(RS1 phase) and (b) $T=0.667$ (RS2 phase) under one realization of
$\bm{J}\in{\cal J}(T)$. In the RS1 phase, the mean phenotypes
$\bm{\mu}^{\pm}$, in particular $\bm{\mu}^{-}$, are highly correlated with
$\bm{\xi}_{1}$, as described by $y=\tanh(\beta\sqrt{N}x)$ (see Fig. 4(a)).
Here, the function $\tanh$ is consistent with the mean-field form of the
magnetization $\mu_{i}^{\pm}=\tanh(\beta\sum_{j\neq i}J_{ij}\mu^{\pm}_{j})$.
Meanwhile, in the RS2 phase, the regulated $\bm{\mu}^{+}$ and non-regulated
$\bm{\mu}^{-}$ states exhibit correlations with $\bm{\xi}_{1}$ and
$\bm{\xi}_{2}$, respectively, as shown in Fig.4(b). In both phases, the
correlations between $\bm{\mu}^{\pm}$ and $\bm{\xi}_{r}$ ($r\geq 3$) are
negligible.

In Fig.5 (a) and (b), we show the temperature dependence of the correlation
between the eigenvectors and $\bm{\mu}^{\pm}$ by introducing the correlation
coefficient between $\\{\xi_{ri}\\}$ and $\\{\mathrm{atanh}(\mu_{i}^{\pm})\\}$
for $r=1,2,3$. Here, the function $\mathrm{atanh}$ is introduced by
considering the $\tanh$-form dependencies of $\bm{\mu}^{\pm}$ on
$\bm{\xi}_{1}$ or $\bm{\xi}_{2}$, as shown in Fig.4. We denote the vector
consisting of $\mathrm{atanh}(\mu_{i}^{\pm})~{}(i=1,\cdots,N)$ as
$\mathrm{atanh}(\bm{\mu}^{\pm})\in\mathbb{R}^{N}$. As shown in Fig.5 (a), the
correlation coefficient between the first eigenvector $\bm{\xi}_{1}$ and the
mean phenotype with regulation $\bm{\mu}^{+}$ increases at $T<T_{0}$, namely
in the RS1 phase. As the temperature is lowered further below $T_{1}$ (towards
the RS2 phase), the correlation between the regulated state and the second
eigenvector increases to be larger than that between the first eigenvector and
the regulated state. Meanwhile, as shown in Fig.5(b), the correlation between
the first eigenvector and the non-regulated state $\bm{\mu}^{-}$ is always
higher than that of other eigenvectors at $T<T_{0}$. For the higher order-
eigenvectors than the second-order, their correlation between the regulated
and non-regulated states is small, as with $\bm{\xi}_{3}$ shown in Fig.5(a)
and (b).

To summarize, typical phenotypes evolved in the RS1 phase are concentrated on
the direction of the first eigenvector, for both with and without regulation.
Meanwhile, in the RS2 phase, the typical phenotypes with and without
regulation are distinctively along the second eigenvector and the first
eigenvector of genotype, respectively. Thus, the typical phenotypes generated
by two distributions $p_{\beta}^{\pm}(\bm{S}|\bm{J})$ are almost orthogonal to
each other.

The contributions of the first and second eigenmodes to the mean phenotypes
$\bm{\mu}^{\pm}$ are given by the magnitudes of their corresponding
eigenvalues. In Fig.5(c), we show the $T$-dependence of the expected value of
the first, second, and third eigenvalues of the evolved genotypes. The
horizontal lines denote their expected values for the symmetric matrices whose
components independently and identically obey the uniform distribution over
$\Omega_{J}$. For $T<T_{0}$, the first eigenvalue shows distinct increases
from the expected value, whereas for $T\simeq T_{1}$, the second eigenvalue
increases. Meanwhile, the third and higher-order eigenvalues show slight
changes. Therefore, the two desirable phenotypes are achieved by the
contribution of the first and second order eigenmodes.

Figure 5: The correlation coefficient between the first, second, and third
eigenvectors with (a) regulated state $\bm{\mu}^{+}$ and (b) non-regulated
state $\bm{\mu}^{-}$. The $T$-dependence of the first, second, and third
eigenvalues of the evolved genotypes is shown in (c), where the three
horizontal lines represent the expected eigenvalues for randomly generated
$\bm{J}$. The vertical dashed line, one-dot-chain line, and two-dot-chain line
denote $T_{0}$, $T_{2}$ and $T_{1}$, respectively. Each point is averaged over
100 samples of the evolved $\bm{J}$.

Following these observations, we map the mean phenotypes with and without
regulation onto the two-dimensional space spanned by the first and second
eigenvectors, $\bm{\xi}_{1}$ and $\bm{\xi}_{2}$, of the evolved genotypes. In
the RS1 and RS2 phases, characteristic mappings are observed as shown in Fig.6
(a), where the mean phenotypes with and without regulation are denoted by
$\star$ and $\bullet$, respectively. In the RS1 phase, the first eigenvector
is dominant to express both mean phenotypes with and without regulation for
most of the evolved genotypes. We term this case as the overlapped phenotypes
(Fig.6 (a) left). Meanwhile, in the RS2 phase, phenotypes are shaped by the
first and second eigenvectors of the evolved genotypes. These genotypes can
satisfy the required fitness conditions both without and with regulation,
respectively. We term the case as separated phenotypes as shown in Fig.6 (a)
right. Hereafter, we term the genotype $\bm{J}$ that gives overlapped and
separable phenotypes as type J1 and type J2, respectively.

Fig.6(b) shows the temperature dependence of the fraction of the type J1 and
J2 genotypes among the ensemble of evolved genotypes ${\cal J}(T)$. At
sufficiently large $T$, their fractions are equal to 0.25, which is indicated
by horizontal lines. The value of 0.25 is the expected value of the fraction
of type J1 and J2 for the randomly generated $\bm{J}$s, as there are two other
cases of mapping; the case that $\star$ and $\bullet$ located along
$\bm{\xi}_{2}$, and that $\star$ and $\bullet$ are along $\bm{\xi}_{1}$ and
$\bm{\xi}_{2}$, respectively. As $T$ decreased toward the RS1 phase, the
fraction of genotypes of type J1 increases up to 0.8. By lowering the
temperature further in the RS2 phase, the dominant genotype is replaced by
type J2. For lower $T<T_{1}$, the dominancy of type J2 decreases as $T$
decreases, and the fraction of types J1 and J2 approaches 0.25.

Figure 6:  (a) Characteristic mapping of the mean phenotypes with regulation
($\star$) and without regulation ($\bullet$), where the diagonal dashed line
with a 45-degree slope is a guide for the eyes. (b) Fraction of genotypes of
types J2 (cyan, left axis) and J1 (yellow, right axis). The horizontal lines
represent 0.25, which is the trivial value for randomly distributed genotypes.
The vertical dashed line, one-dot chain, and two-dot chain line represent
$T_{0}$, $T_{2}$, and $T_{1}$, respectively.

###  III.3 Why the separable phenotype appears at $T<T_{1}$?

Here, we discuss why the type J1 and J2 genotypes are dominant at
$T_{1}<T<T_{0}$ and $T_{2}<T<T_{1}$, respectively. To answer this question, we
observe the fitness of the evolved genotypes ${\cal J}(T)$ under a trial
temperature $T_{\mathrm{tr}}$. The evolutionary process in our model selected
genotypes among possible $\bm{J}$s; hence, $\bm{J}\in{\cal J}(T)$ can be a
candidate for genotypes in ${\cal J}(T_{\mathrm{tr}})$ ($T\neq
T_{\mathrm{tr}}$), in principle. By evaluating the fitness of $\bm{J}\in{\cal
J}(T)$ at a different temperature $T_{\mathrm{tr}}$, we discuss the reason why
$\bm{J}\in{\cal J}(T)$ cannot be selected at different temperatures.

Fig.7 shows the $T_{\mathrm{tr}}$ dependence of the fitness $\Psi$ on
$\bm{J}\in{\cal J}(T)$ for $T=0.91$ (type J1; RS1) and $T=0.63$ (type J2;
RS2). At sufficiently large $T_{\mathrm{tr}}$, J1 and J2 fitness values are
not much different. Type J1 and J2 are subject to one- and two-dimensional
constraints, respectively, hence, the possible configurations of type J1 are
larger than that of type J2. From the thermodynamic perspective, the dominance
of the type J1 in the RS1 phase is caused by this entropic effect. In the RS2
phase, the fitness of type J2 is sufficiently large to overcome the entropic
effect, and they can be dominant in this phase. This observation indicates
that the changes in the ensemble of ${\cal J}(T)$ can be regarded as a phase
transition with respect to genotypes between the type J1 and J2.

Figure 7: Trial temperature $T_{\mathrm{tr}}$ dependence of the fitness for
$\bm{J}\in{\cal J}(T)$ at $T=0.91$ (denoted by J1) and $T=0.63$ (denoted by
J2). The shaded region indicates the difference between $\langle
m_{A}^{-}\rangle_{-}$ and $\langle m_{A}^{+}\rangle_{+}$.

###  III.4 Evolutional dynamics of genotypes on the two-dimensional plane

For the understanding of the evolutionary construction of the separated
phenotypes, we simplify the evolutionary dynamics using the two-dimensional
space spanned by the first and second eigenvectors, although the two-
dimensional approximation was not necessarily accurate in the early stages of
evolution, even in the RS2 phase. In Fig.8, we show evolutionary change of the
mean phenotypes in RS2 phase corresponding to the series shown in Fig.2 (b),
where $\langle m_{A}^{-}\rangle$ increased before $\langle m_{A}^{+}\rangle$.
The panels of Fig.8 show the time evolution of the mean phenotypes
$\bm{\mu}^{+}(\bm{J}^{(g)})$ ($\bigcirc$) and $\bm{\mu}^{-}(\bm{J}^{(g)})$
($\square$) mapped onto the two-dimensional space spanned by the first and
second eigenvectors of the genotype at each generation denoted in the panels.
The localization of the mean phenotype without regulation appears on the first
eigenvector 141-500 generations before that of the regulation case. From
generations 641-800, the contribution of the second eigenvector to the mean
phenotype increases with regulation. After the reorganization of the
distributions at generations 801-980, the characteristic phenotype mapping for
the type J2 genotype appeares.

When $\langle|m_{A}^{+}|\rangle_{+}$ increases before the increase of
$\langle|m_{A}^{+}|\rangle_{-}$, the localization of $\bm{\mu}^{+}$ on the
second eigenvector appears in the early stage of evolution. Additionally, the
localization of $\bm{\mu}^{-}$ follows with the reorganization of
$\bm{\mu}^{+}$. When both $\langle|m_{A}^{\pm}|\rangle_{\pm}$ increased
simultaneously, $\bm{\mu}^{\pm}$ are localized almost simultaneously (see
supplement material).

Figure 8: The evolution of the mean phenotypes. Corresponding to the evolution
of Fig.2 (b), this figure shows the evolution of the mean phenotype in the
two-dimensional space spanned by the first and second eigenvectors of the
genotype for the evolution generations $[1-140]$, $[141-500]$,
$[501-640],\cdots,[981-1100]$. The results with and without regulations are
plotted by $\circ$ and $\square$, respectively.

##  IV Switching trajectory

Figure 9: (a) MC steps required for the switching from regulated to non-
regulated state ($\bigcirc$) and non-regulated to regulated state ($\square$).
The dashed vertical line, two-dotted chain line, and one-dotted chain line
denote $T_{0}$, $T_{1}$, and $T_{2}$, respectively. The inset magnifies the
difference between the RS1 and RS2 phase. The transition time to shift the
active sites from the state without regulation to that with regulation is
evaluated as follows. After the sufficient time updates of $\bm{S}$ under the
non-regulated condition $\bm{S}_{\cal R}\notin\bm{S}_{\cal R}^{+}$, the
regulatory sites are changed to $\bm{S}_{\cal R}\in\bm{S}_{\cal R}^{+}$, and
then $\bm{S}$ (except the regulatory region) is updated according to (1). We
compute the target magnetization $|\sum_{i\in{\cal A}}S_{i}/N_{A}|$ at each MC
step to obtain the step where $|\sum_{i\in{\cal A}}S_{i}/N_{A}|$ first reaches
the value $\langle|m_{A}^{+}|\rangle_{+}$, which is defined as the transition
time. (b) Switching trajectories of local magnetizations projected to the
first and the second eigenvectors defined on an evolved genotype at $T=0.67$
(RS2). The component projected onto the third order eigenvector is denoted by
dashed lines.

Figure 10: Heat maps on the two-dimensional space for switching trajectories
from regulated state to non-regulated state, defined on an evolved genotype at
$T=0.67$ (RS2). Here, the two-dimensional space is meshed by 0.01, and
$\log_{10}$-frequencies of the trajectories during the given steps are
plotted. $\star$ and $\bullet$ denote the regulated state and non-regulated
state projected onto the two-dimensional space, respectively. The lines below
these points represent the first and second eigenmodes of fluctuation around
these points. Here, the length of the lines is magnified to be discernible.
However, the ratio of the lines is proportional to the root of the ratio of
the eigenvalues.

Under the $\bm{J}$s of type J2, the shift between the regulated and non-
regulated states involves a large conformational change. We employ the MCMC
method according to (1) for simulating the transition dynamics from regulated
to non-regulated cases, or from non-regulated to regulated cases, and compute
the MC steps required for the shift between the two cases.

Fig. 9(a) shows the transition time calculated by the MCMC method from the
regulated to non-regulated states ($\bigcirc$), and from the regulated to non-
regulated states ($\square$), respectively. Here, the upper limit of the MC
step is set at $10^{5}$. In the RS1 phase, there is little change in
phenotypes with and without regulation, and the transition time is within 20
steps. Compared with the RS1 phase, the transition time required in the RS2
phase increases. This increase in relaxation time is associated with the large
conformational change of phenotype under the J2 genotypes. However, the large
conformational change does not qualitatively change the relaxation time. As in
the RS1 phase, the relaxation time in the RS2 phase is in the order of
$10^{2}$. In the RSB phase, the MC steps required for switching diverge as $T$
decreases. This phenomenon in the RSB phase is consistent with the property of
the RSB phase where the metastable states hamper relaxation.

The trajectories shifting between two states lie in $2^{N}$-dimensional space.
However, particularly in the RS2 phase, the two-dimensional space spanned by
the first and second eigenvectors of the evolved genotype is sufficient to
describe the switching trajectories. This low-dimensional constraint is
already observed as the equilibrium property in the RS2 phase, as shown in
Fig.5 (a) and (b). Fig. 9 (b) shows the trajectories of the components
projected onto the first ($\circ$), second ($\square$), and third (dashed
line) eigenvectors defined on an evolved $\bm{J}$ of type J2 at $T=0.67$
(RS2). In the regulate-to-non-regulate switching, the change in the first
component is much larger, and in the non-regulate-to-regulate switching, the
change in the second component is much larger. Meanwhile, the third-order (and
higher) components are nearly constant during regulated-to-non-regulated or
non-regulated-to-regulated switching.

We generate 1000 switching trajectories on a certain $\bm{J}\in{\cal J}(T)$,
and map them onto the two-dimensional space spanned by the first and second
eigenvectors of the evolved genotypes. Fig.10 shows the heat map on the two-
dimensional space for the switching trajectories defined on an evolved
genotype of type J2 at $T=0.68$ (RS2) from non-regulated to regulated states.
The mean phenotypes with and without regulation, $\bm{\mu}^{+}$ and
$\bm{\mu}^{-}$, after sufficient time steps of updating are denoted by $\star$
and $\bullet$, respectively. Additionally, the direction of the fluctuation of
these points is indicated by two lines below the points. The switching
trajectory when regulation is removed is shown in the supplementary material.
For both cases of switching, most of the trajectories follow a quarter-circle
path. This quarter-circle path is restricted to a one-dimensional path within
the two-dimensional space. With this restriction, the transition time between
the two states remains small, even though the two phenotypes are far apart, as
shown in Fig. 9(a).

The quarter-circle path on the two-dimensional plane restricts the
trajectories of the convergence from arbitrary initial conditions to the
phenotypes with and without regulation. The heat map for the relaxation
dynamics on a type J2 genotype evolved at $T=0.67$ from arbitrary initial
conditions is shown in Fig. 11 for the regulated case. Most of the
trajectories are attracted once to the quarter-circle line where the switching
paths are concentrated, and then approach the regulated state. (The relaxation
dynamics of the non-regulated phenotype are shown in the supplementary
material.) The quarter-circle path is attractive in the sense that any state
tends towards the regulated or non-regulated state through this path.

Figure 11: Heat maps on the two-dimensional space for relaxation trajectories
from an initial condition to the regulated state defined on the evolved
genotype of type J2 at $T=0.67$ (RS2), the same value as adapted in Fig.10.
$\star$, $\bullet$, and the orthogonal lines below these points are the same
as in Figs.10.

##  V Two dimensional approximation of free energy landscape

To understand the characteristic switching path in the two-dimensional space,
we examine the free energy landscape. The free energies for the regulated and
non-regulated cases, denoted by $f_{+}$ and $f_{-}$, are defined as

$\displaystyle f_{+}$
$\displaystyle=-\frac{1}{N\beta}\ln\sum_{\bm{S}|\bm{S}_{{\cal
R}}\in\bm{S}_{{\cal R}}^{+}}\exp(-\beta H),$ (11) $\displaystyle f_{-}$
$\displaystyle=-\frac{1}{N\beta}\ln\sum_{\bm{S}|\bm{S}_{{\cal
R}}\notin\bm{S}_{{\cal R}}^{+}}\exp(-\beta H).$ (12)

Following the result of the numerical simulations, we consider the two-rank
approximation of the evolved $\bm{J}$ as
$\bm{J}\simeq\lambda_{1}\bm{\xi}^{1}{\bm{\xi}^{1}}^{\top}+\lambda_{2}\bm{\xi}^{2}{\bm{\xi}^{2}}^{\top}$.
Under the two-rank approximation, the Hamiltonian is given by

$\displaystyle
H=-\sum_{k=1}^{2}\frac{\lambda_{k}}{2}\left\\{\left(\sum_{i=1}^{N}\xi_{i}^{k}S_{i}\right)^{2}-1\right\\}.$
(13)

For the two-rank approximation form, one can represent the free energy as a
function of $m_{1}$ and $m_{2}$ defined by

$\displaystyle m_{1}^{\pm}$
$\displaystyle=\frac{1}{\sqrt{N}}{\bm{\xi}^{1}}^{\top}\bm{\mu}^{\pm},$ (14)
$\displaystyle m_{2}^{\pm}$
$\displaystyle=\frac{1}{\sqrt{N}}{\bm{\xi}^{2}}^{\top}\bm{\mu}^{\pm},$ (15)

where $m_{1}^{+}$ and $m_{2}^{+}$ correspond to the projection of the local
magnetization with regulation onto the first and second eigenvectors,
respectively, whereas $m_{i}^{-}~{}(i=1,2)$ are those without regulation.
Following the calculation shown in the Appendix, the free energies are given
by

$\displaystyle f_{+}$
$\displaystyle=\\!\\!\sum_{k=1}^{2}\frac{\lambda_{k}{m_{k}^{+}}^{2}}{2}-\frac{1}{N\beta}\left\\{\sum_{i=1}^{N}\ln\\!\left(2\cosh(\beta
h_{i}^{+})\right)\\!+\\!\log(p_{+}^{+}+p_{-}^{+})\\!\right\\}+\frac{1}{N}\sum_{i<j,i,j,\in{\cal
R}}J_{ij}$ (16) $\displaystyle f_{-}$
$\displaystyle=\\!\sum_{k=1}^{2}\frac{\lambda_{k}{m_{k}^{-}}^{2}}{2}-\frac{1}{N\beta}\left\\{\sum_{i=1}^{N}\ln\left(2\cosh(\beta
h_{i}^{-})\right)+\ln\left(1-\left(p_{+}^{-}+p_{-}^{-}\right)\right)\right\\}$
(17)

where
$h_{i}^{\pm}=\lambda_{1}m_{1}^{\pm}\sqrt{N}\xi_{i}^{1}+\lambda_{2}m_{2}^{\pm}\sqrt{N}\xi_{i}^{2}$
and

$\displaystyle p_{\pm}^{+}=\prod_{i\in{\cal R}}\frac{\exp(\pm\beta
h_{i}^{+})}{2\cosh(\beta h_{i}^{+})},~{}~{}~{}p_{\pm}^{-}=\prod_{i\in{\cal
R}}\frac{\exp(\pm\beta h_{i}^{-})}{2\cosh(\beta h_{i}^{-})}.$ (18)

The saddle point equations for $m_{k}^{\pm}$ are given by

$\displaystyle m_{k}^{+}$ $\displaystyle=\frac{1}{\sqrt{N}}\sum_{i\notin{\cal
R}}\xi_{i}^{k}\tanh(\beta h_{i}^{+})+\frac{1}{N}\sum_{j\in{\cal
R}}\xi_{j}^{k}\frac{p_{+}^{+}-p_{-}^{+}}{p_{+}^{+}+p_{-}^{+}}$ (19)
$\displaystyle m_{k}^{-}$ $\displaystyle=\frac{1}{\sqrt{N}}\sum_{i\notin{\cal
R}}\xi_{i}^{k}\tanh(\beta h_{i}^{-})+\frac{1}{\sqrt{N}}\sum_{i\in{\cal
R}}\xi_{i}^{k}\frac{\tanh(\beta
h_{i}^{-})-(p_{+}^{-}-p_{-}^{-})}{1-\left(p_{+}^{-}+p_{-}^{-}\right)}.$ (20)

Figure 12: Two-dimensional approximation of the free energy landscape of one
genotype with separable phenotype space for (a) regulated case and (b) non-
regulated case at $p_{A}=0.05$ and $p_{R}=0.1$. $\star$ and $\bullet$ show the
projection of the regulated state and non-regulated state.

Fig.12 (a) and (b) show the landscape of $f^{+}$ and $f^{-}$, respectively,
plotted on the two-dimensional space of one evolved J2 genotype in the RS2
phase at $T=0.67$. The minima of the free energies are consistent with the
numerically observed phenotypes with and without regulation, which are
indicated by $\star$ and $\bullet$; hence, the two-dimensional approximation
of the free energy is valid. As shown in Fig.12, along the quarter-circle
shape that connects the regulated and non-regulated states, the free energy
remains small. The trajectories shown in Fig. 10 are restricted to this
quarter-circle, wherein free energy is small.

###  V.1 Rough approximation of free energy

What property of the evolved genotype gives the quarter-circle shape of the
free energy landscape. We introduce the following assumptions.

A1

The difference between the first and second eigenvalues is negligible.

A2

The components of the eigenvectors $\bm{\xi}^{1}$ and $\bm{\xi}^{2}$ are
independently and identically distributed, according to the Gaussian
distribution ${\cal N}(0,1/N)$.

A3

The active and regulatory sites are negligible.

Following assumption A1, we replaced the first and second eigenvalues
$\lambda_{1}$ and $\lambda_{2}$ with their mean $\overline{\lambda}$. Under
these assumptions, we obtained the following form of free energy, as explained
in the Appendix;

$\displaystyle
f_{\mathrm{app}}=\frac{\overline{\lambda}}{2}(m_{1}^{2}+m_{2}^{2})-\frac{1}{\beta}\int
Dz\ln\cosh\left(\beta\overline{\lambda}\sqrt{m_{1}^{2}+m_{2}^{2}}z\right),$
(21)

where $Dz=\frac{dz}{\sqrt{2\pi}}\exp\left(-\frac{z^{2}}{2}\right)$. The form
of (21) indicates that the free energy under the approximations A1-A3 depends
on $m_{1}$ and $m_{2}$ through $m_{S}\equiv\sqrt{m_{1}^{2}+m_{2}^{2}}$. Hence,
the approximated free energy has the same value according to $m_{S}$, even
when the individual values of $m_{1}$ and $m_{2}$ are different. The saddle
point of $m_{S}$ is given by

$\displaystyle m_{S}=\int Dz~{}z\tanh(\beta\overline{\lambda}m_{S}z).$ (22)

Figure 13:  Free energy landscape on the two-dimensional plane under the
assumptions A1-A3 at $T=0.67$ (RS2). The solid line represents the minimum
free energy and the star and circle represent the mean of phenotypes with and
without regulation, respectively.

Fig. 13 shows the free energy landscape under the assumptions A1-A3 defined on
a J2 genotype evolved at $T=0.67$ (RS2). As expected from the form of eq.
(21), the approximated free energy shows a quarter-circle landscape. The
quarter-circle curve represents the minimum of the free energy
$f_{\mathrm{app}}$, whereas the equilibrium states with and without regulation
are denoted by stars and circles, respectively, which are located near the
extremum line of $f_{\mathrm{app}}$. Therefore, the one-dimensional and
quarter-circle switching path is considered to be provided by the free sites,
as the active and regulatory sites are ignored in deriving $f_{\mathrm{app}}$
(assumption A3). A particular difference between $f_{\pm}$ and
$f_{\mathrm{app}}$ is that the valleys around the mean of the phenotype with
and without regulation (Fig. 12) cannot be described by $f_{\mathrm{app}}$.
For the description of these valleys, it is necessary to consider the active
and regulatory sites. Thus, free energy consists of quarter-circle switching
paths provided by free sites and valleys around the mean phenotypes provided
by active and regulatory sites. Further, the assumption A2 suggests that
randomness in the embedded pattern in the free sites is significant for the
description of the quarter-circle path. Therefore, as the number of free sites
decreases, or equivalently as the number of active and regulatory sites
increases, the description by $f_{\mathrm{app}}$ would be invalid.

##  VI Summary and Conclusion

In this study, we investigated the evolution of a spin model to generate two
specific configurations of active sites depending on the regulation. A fitness
function was designed to increase when the appropriate spin configurations
(phenotypes) with and without regulation appears with high probability. Our
analysis revealed three transition points, $T_{0}$, $T_{1}$ and $T_{2}$: The
fitness was increased from the trivial value for $T<T_{0}$. For
$T_{2}<T<T_{0}$, the evolved system belonged to the RS phase. The RS phase was
further divided into two regions at $T=T_{1}$, the RS1 ($T_{1}<T<T_{0}$) and
RS2 ($T_{2}<T<T_{1}$) phases, with dominant genotypes differing in these
regions, type J1 for RS1 and J2 for RS2 phases. For $T_{1}<T<T_{0}$, the
phenotypes, i.e., spin configurations, other than active sites barely depended
on the regulation. In contrast, for $T_{2}<T<T_{1}$, the two phenotypes with
and without regulation, showed a large difference, contrasting the small
difference in the RS1 phase.

In the RS2 phase, the two phenotypes were provided by using the first and
second eigenmodes to express non-regulated and regulated phenotypes,
respectively, where the switching path between the two phenotypes can be
described by the first and second eigenmodes of the two endpoint phenotypes. A
one-dimensional quarter-half shape switching path connected the two endpoint
phenotypes in the two-dimensional space spanned by the first and second
eigenvectors of the J2 genotype. This switching path was robust to
perturbations, in the sense that any trajectories deviating from the path were
attracted to the path. Evolutionary construction of this one-dimensional path
met the requirements of plasticity against regulatory changes and robustness
in phenotypes. Further, the low-dimensionality of the switching path allowed
for quick switching between two stable phenotypes depending on the regulation.

To understand the evolutionary origin of the one-dimensional switching path,
we applied a two-dimensional approximation to the free energy landscape for
the evolved genotype in the RS2 phase. By only considering randomness in the
free sites of two endpoint phenotypes, it was found that the free energy takes
a minimum along a quarter-circle shape in two dimensions. The two endpoint
phenotypes were located near the quarter-circle path, and the switching
trajectories followed the valley of the free energies connecting the two
endpoints. In this case, the minima relate to the sites that were active and
regulated. The cooperative evolution of the active, regulatory, and free sites
provided stable expression of the endpoint phenotypes and robust switching
paths.

Our findings suggest that low dimensionality plays a crucial role in achieving
both stable expressions of two phenotypes and large conformational changes
over a stable path. This leads to the acquisition of both robustness and
plasticity. Constraints on adaptive changes in phenotypes upon environmental
and evolutionary changes have recently gathered much attention [29, 44, 31,
45]. The constraint attracts a low-dimensional subspace within the high-
dimensional space, supporting the robustness. Here, we demonstrated that the
state change relevant to function is facilitated by the one-dimensionally
constrained path on the two-dimensional plane, which allows large-amplitude
plastic motion that is advantageous for functional changes. Notably, this
constrained path is already “prepared” as a relaxation path during the
evolution course (Fig.8).

Previous studies have demonstrated that genotypes providing a single function
by expressing a specific phenotype can evolve in the RS phase [31]. In our
study, we found the transition that occurs in the RS phase for two-functional
phenotypes. The genotypes that achieve switching between two functional
phenotypes depending on the regulation were dominant in the RS2 phase, i.e.,
in the temperature region closer to the RSB within the RS phase. For the
evolution to achieve more functions, further transitions within the RS phase
can be expected. With such successive transitions, the genotype will approach
the RSB transition point, where further plasticity will be achieved. This may
be consistent with the observation of critical behavior in protein dynamics
[30], wherein plasticity and robustness are compatible.

Here, we did not impose any driving force to create the one-dimensional
switching path; rather, the evolution under the fitness defined by the two
endpoint phenotypes resulted in genotypes that provide not only stable
expression of the phenotypes but also robust and plastic switching. This
observation presents the possibility of the evolutionary construction of
proteins [46] with allosteric effects based on the binding ability of the
active site, under conditions characterized by the RS phase, in addition to
synthetic approaches [47]. Further analysis of interacting spin systems that
achieve robust multiple functions is essential for the evolution of proteins
and material design [48, 49].

Investigation of the microscopic properties of evolved genotypes is an
important future research direction. However, the focus of this study was on
the extraction of macroscopic low-dimensional structures. Frustration is a
potential measure to characterize the genotype, which captures consistency in
interactions, and the increase in frustration can indicate a rugged landscape.
Generally, as the number of embedded patterns increases, the level of
frustration in the interactions increases [50, 51]. We observed an increase in
frustration in our model in comparison to the one-desirable phenotype case
(see Supplement). In actual proteins, steric frustration can be utilized by
multisubstrate enzymes to facilitate the rate-limiting product-release step
[52]. Understanding the relationship between frustration and the number of
embedded patterns may provide insights into the properties of real
biomolecules.

The evolutionary spin model considered in this study is rather simple and
abstract. There is room to consider more realistic settings and discuss the
generality of the results. For instance, several biological molecules have
multiple regulatory or active sites, and their phenotype expression is more
complicated. G protein-coupled receptors show dual ligand binding events where
the binding of one ligand enhances that of the other [53, 54]. Thiamine
diphosphate in the two active sites of pyruvate dehydrogenase complex can
communicate with each other over a distance of 20 angstroms using a proton to
switch the conformation [55]. The contribution of these kinds of cooperation
to the evolution of robustness and plasticity needs to be revealed. In
contrast to the global (all-to-all) coupling model, the study of models with
spatially localized interactions is also important [56, 57, 58].

To conclude, we have shown that the stable expression and switching of
phenotypes takes advantage of evolutionary constructed low-dimensional
phenotypic constraints, with which robustness and plasticity are compatible.
Our finding indicates that the evolution of low-dimensionality can be a
unified view for the understanding of evolutional phenomena.

##  Appendix A Derivation of the free energy density

We introduce the equality

$\displaystyle 1$ $\displaystyle=\int
dm_{k}\delta\left(m_{k}-\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\xi_{i}^{k}S_{i}\right)$

for $k=1,2$ to the partition function with regulation as

$\displaystyle Z_{+}$ $\displaystyle=\int
dm_{1}^{+}dm_{2}^{+}\sum_{\bm{S}|\bm{S}_{\cal R}^{+}\in\bm{S}_{\cal
R}}\prod_{k=1}^{2}\delta\left(m_{k}-\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\xi_{i}^{k}S_{i}\right)\exp\left[\sum_{k=1}^{2}\frac{\beta\lambda_{k}}{2}\left(m_{k}^{2}-1\right)\right]$
$\displaystyle=\int
dm_{1}^{+}dm_{2}^{+}d\hat{m}_{1}^{+}d\hat{m}_{2}^{+}\sum_{\bm{S}|\bm{S}_{\cal
R}^{+}\in\bm{S}_{\cal
R}}\prod_{k=1}^{2}\exp\left(-Nm_{k}\hat{m}_{k}+\sqrt{N}\hat{m}_{k}\sum_{i=1}^{N}\xi_{i}^{k}S_{i}+\frac{\beta\lambda_{k}}{2}(m_{k}^{2}-1)\right)$
$\displaystyle=\int
dm_{1}^{+}dm_{2}^{+}d\hat{m}_{1}^{+}d\hat{m}_{2}^{+}\prod_{k=1}^{2}\exp\left(-Nm_{k}\hat{m}_{k}+\frac{\beta\lambda_{k}}{2}(m_{k}^{2}-1)\right)$
$\displaystyle\hskip
28.45274pt\times\prod_{i=1}^{N}\left\\{2\cosh(\sqrt{N}\sum_{k=1}^{2}\hat{m}_{k}\xi_{i}^{k})\right\\}\times\left\\{p_{+}^{+}+p_{-}^{+}\right\\}.$
(23)

The integrals are implemented by utilizing the saddle point method.

For the without-regulation case, we obtain

$\displaystyle Z_{-}$ $\displaystyle=\int
dm_{1}^{-}dm_{2}^{-}d\hat{m}_{1}^{-}d\hat{m}_{2}^{-}\sum_{\bm{S}|\bm{S}_{\cal
R}\notin\bm{S}_{\cal
R}^{+}}\prod_{k=1}^{2}\exp\left(-Nm_{k}\hat{m}_{k}+\sqrt{N}\hat{m}_{k}\sum_{i=1}^{N}\xi_{i}^{k}S_{i}+\frac{\beta\lambda_{k}}{2}(m_{k}^{2}-1)\right)$
$\displaystyle=\int
dm_{1}^{-}dm_{2}^{-}d\hat{m}_{1}^{-}d\hat{m}_{2}^{-}\prod_{k=1}^{2}\exp\left(-Nm_{k}\hat{m}_{k}+\frac{\beta\lambda_{k}}{2}(m_{k}^{2}-1)\right)$
$\displaystyle\hskip
28.45274pt\times\prod_{i=1}^{N}\left\\{2\cosh(\sqrt{N}\sum_{k=1}^{2}\hat{m}_{k}\xi_{i}^{k})\right\\}\times\left\\{1-(p_{+}^{-}+p_{-}^{-})\right\\}.$
(24)

###  A.1 Free energy under assumptions A1-A3

Under the assumptions A1 and A3, $Z^{+}=Z^{-}$ holds, and we denote it as
$Z_{\pm}$ given by

$\displaystyle Z_{\pm}$
$\displaystyle=\sum_{\bm{S}}\exp\left(\beta\overline{\lambda}\sum_{i<j}(\xi_{i}^{1}\xi_{j}^{1}+\xi_{i}^{2}\xi_{j}^{2})S_{i}S_{j}\right)$
$\displaystyle\sim\sum_{\bm{S}}\int
dm_{1}dm_{2}\delta\left(m_{1}-\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\xi_{i}^{1}S_{i}\right)\delta\left(m_{2}-\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\xi_{i}^{2}S_{i}\right)\exp\left(\displaystyle\frac{\beta\overline{\lambda}N}{2}\left(m_{1}^{2}+m_{2}^{2}\right)\right)$
$\displaystyle=\int
dm_{1}dm_{2}d\hat{m}_{1}d\hat{m}_{2}\exp(-N\hat{m}_{1}m_{1}-N\hat{m}_{2}m_{2})\exp\left(\displaystyle\frac{\beta\overline{\lambda}N}{2}\left(m_{1}^{2}+m_{2}^{2}\right)\right)$
$\displaystyle\hskip
28.45274pt\times\prod_{i}2\cosh(\hat{m}_{1}\sqrt{N}\xi_{i}^{1}+\hat{m}_{2}\sqrt{N}\xi_{i}^{2})$
$\displaystyle=\int
dm_{1}dm_{2}\exp\left(-\displaystyle\frac{\beta\overline{\lambda}N}{2}\left(m_{1}^{2}+m_{2}^{2}\right)+\sum_{i}\ln\cosh(\beta\overline{\lambda}\sqrt{N}(m_{1}\xi_{i}^{1}+m_{2}\xi_{i}^{2}))\right),$
(25)

Here, we introduced the saddle point methods for the integral with respect to
$\hat{m}_{1}$ and $\hat{m}_{2}$. Under the assumption A2 at sufficiently large
system size, the summation with respect to the components of the eigenvectors
can be replaced with the integral according to the Gaussian distribution. as

$\displaystyle Z_{\pm}$ $\displaystyle=\int
dm_{1}dm_{2}\exp\left(-\displaystyle\frac{\beta\overline{\lambda}N}{2}\left(m_{1}^{2}+m_{2}^{2}\right)+N\int
Dz\ln\cosh\left(\beta\overline{\lambda}\sqrt{m_{1}^{2}+m_{2}^{2}}z\right)\right).$
(26)

By introducing the saddle point method to the integrals of $m_{1}$ and
$m_{2}$, we obtaine the approximated free energy $f_{\mathrm{app}}$.

###### Acknowledgements.

The authors thank Qian Yuan Tang and Tuan Pham for helpful comments and
discussions. This work is partially supported by Grant in-Aid for Scientific
Research (A) (20H00123) from the Ministry of Education, Culture, Sports,
Science, and Technology (MEXT) of Japan. KK is also supported by the Novo
Nordisk Foundation.

## References

  * Stout  _et al._ [2005] R. D. Stout, C. Jiang, B. Matta, I. Tietzel, S. K. Watkins, and J. Suttles, Macrophages sequentially change their functional phenotype in response to changes in microenvironmental influences, The Journal of Immunology 175, 342 (2005).
  * Pisarchik and Feudel [2014] A. N. Pisarchik and U. Feudel, Control of multistability, Physics Reports 540, 167 (2014).
  * Changeux [1961] J. P. Changeux, The feedback control mechanisms of biosynthetic-threonine deaminase by l-isoleucine, in  _Cold Spring Harbor symposia on quantitative biology_ (Cold Spring Harbor Laboratory Press, 1961) pp. 313–318.
  * Monod  _et al._ [1965] J. Monod, J. Wyman, and J. P. Changeux, On the nature of allosteric transitions: a plausible model, J. Mol. Biol. 12, 88 (1965).
  * Changeux [2013] J.-P. Changeux, 50 years of allosteric interactions: the twists and turns of the models, Nature reviews Molecular cell biology 14, 819 (2013).
  * Fischer  _et al._ [2011] S. Fischer, K. W. Olsen, K. Nam, and M. Karplus, Unsuspected pathway of the allosteric transition in hemoglobin, Proceedings of the National Academy of Sciences 108, 5608 (2011), https://www.pnas.org/doi/pdf/10.1073/pnas.1011995108 .
  * Maragakis and Karplus [2005] P. Maragakis and M. Karplus, Large amplitude conformational change in proteins explored with a plastic network model: adenylate kinase, Journal of molecular biology 352, 807 (2005).
  * Karplus and Gao [2004] M. Karplus and Y. Q. Gao, Biomolecular motors: the f1-atpase paradigm, Current opinion in structural biology 14, 250 (2004).
  * Markevich  _et al._ [2004] N. I. Markevich, J. B. Hoek, and B. N. Kholodenko, Signaling switches and bistability arising from multisite phosphorylation in protein kinase cascades , Journal of Cell Biology 164, 353 (2004), https://rupress.org/jcb/article-pdf/164/3/353/1312675/jcb1643353.pdf .
  * Hatakeyama and Kaneko [2020] T. S. Hatakeyama and K. Kaneko, Transition in relaxation paths in allosteric molecules: Enzymatic kinetically constrained model, Physical Review Research 2, 012005(R) (2020).
  * Stern  _et al._ [2020] M. Stern, M.B. Pinson, and A. Murugan, Continual learning of multiple memories in mechanical networks, Physical Review X 10, 031044 (2020).
  * Laine  _et al._ [2010] E. Laine, C. Goncalves, J. C. Karst, A. Lesnard, S. Rault, W.-J. Tang, T. E. Malliavin, D. Ladant, and A. Blondel, Use of allostery to identify inhibitors of calmodulin-induced activation of bacillus anthracis edema factor, Proceedings of the National Academy of Sciences 107, 11277 (2010).
  * Laine  _et al._ [2012] É. Laine, L. Martínez, D. Ladant, T. Malliavin, and A. Blondel, Molecular motions as a drug target: mechanistic simulations of anthrax toxin edema factor function led to the discovery of novel allosteric inhibitors, Toxins 4, 580 (2012).
  * Panjkovich and Svergun [2016] A. Panjkovich and D. I. Svergun, Deciphering conformational transitions of proteins by small angle x-ray scattering and normal mode analysis, Physical Chemistry Chemical Physics 18, 5707 (2016).
  * van der Vaart [2006] A. van der Vaart, Simulation of conformational transitions, Theoretical Chemistry Accounts 116, 183 (2006).
  * Orellana [2019] L. Orellana, Large-scale conformational changes and protein function: breaking the in silico barrier, Frontiers in molecular biosciences 6, 117 (2019).
  * Tirion [1996] M. M. Tirion, Large amplitude elastic motions in proteins from a single-parameter, atomic analysis, Phys. Rev. Lett. 77, 1905 (1996).
  * Bahar  _et al._ [1997] I. Bahar, A. R. Atilgan, and B. Erman, Direct evaluation of thermal fluctuations in proteins using a single-parameter harmonic potential, Folding and Design 2, 173 (1997).
  * Haliloglu  _et al._ [1997] T. Haliloglu, I. Bahar, and B. Erman, Gaussian dynamics of folded proteins, Physical review letters 79, 3090 (1997).
  * Henzler-Wildman  _et al._ [2007] K. A. Henzler-Wildman, V. Thai, M. Lei, M. Ott, M. Wolf-Watz, T. Fenn, E. Pozharski, M. A. Wilson, G. A. Petsko, M. Karplus, _et al._ , Intrinsic motions along an enzymatic reaction trajectory, Nature 450, 838 (2007).
  * Yang and Roux [2008] S. Yang and B. Roux, Src kinase conformational activation: thermodynamics, pathways, and mechanisms, PLoS computational biology 4, e1000047 (2008).
  * Delarue  _et al._ [2017] M. Delarue, P. Koehl, and H. Orland, Ab initio sampling of transition paths by conditioned langevin dynamics, The Journal of chemical physics 147, 152703 (2017).
  * Bergmann  _et al._ [2003] S. Bergmann, J. Ihmels, and N. Barkai, Similarities and differences in genome-wide expression data of six organisms, PLOS Biology 2, E9 (2003).
  * Gunasekera  _et al._ [2008] T. S. Gunasekera, L. N. Csonka, and O. Paliy, Genome-wide transcriptional responses of escherichia coli k-12 to continuous osmotic and heat stresses, Journal of Bacteriol 190, 3712 (2008).
  * Marguerat  _et al._ [2012] S. Marguerat, A. Schmidt, S. Codlin, W. Chen, R. Aebersold, and J. Bähler, Quantitative analysis of fission yeast transcriptomes and proteomes in proliferating and quiescent cells, Cell 151, 671 (2012).
  * Matsumoto  _et al._ [2013] Y. Matsumoto, Y. Murakami, S. Tsuru, B. Y. Ying, and T. Yomo, Growth rate-coordinated transcriptome reorganization in bacteria, BMC Genomics 14, 808 (2013).
  * Kaneko  _et al._ [2015] K. Kaneko, C. Furusawa, and T. Yomo, Universal relationship in gene expression changes for cells in steady-growth state, Physical Review X 5, 011014 (2015).
  * Schmidt  _et al._ [2016] A. Schmidt, K. Kochanowski, S. Vedelaar, E. Ahrné, B. Volkmer, L. Callipo, K. Knoops, M. Bauer, R. Aebersold, and M. Heinemann, The quantitative and condition-dependent escherichia coli proteome, Nat. Biotech. 34, 104 (2016).
  * Furusawa and Kaneko [2018] C. Furusawa and K. Kaneko, Formation of dominant mode by evolution in biological systems, Physical Review E 97, 042410 (2018).
  * Tang and Kaneko [2020] Q.-Y. Tang and K. Kaneko, Long-range correlation in protein dynamics: Confirmation by structural data and normal mode analysis, PLoS computational biology 16, e1007670 (2020).
  * Sakata and Kaneko [2020] A. Sakata and K. Kaneko, Dimensional reduction in evolving spin-glass model: Correlation of phenotypic responses to environmental and mutational changes, Phys. Rev. Lett. 124, 218101 (2020).
  * Sakata  _et al._ [2009a] A. Sakata, K. Hukushima, and K. Kaneko, Funnel landscape and mutational robustness as a result of evolution under thermal noise, Phys. Rev. Lett. 102, 148101 (2009a).
  * Sakata  _et al._ [2009b] A. Sakata, K. Hukushima, and K. Kaneko, Statistical-mechanical study of evolution of robustness in noisy environments, Phys. Rev. E 80, 051919 (2009b).
  * Sakata  _et al._ [2012] A. Sakata, K. Hukushima, and K. Kaneko, Replica symmetry breaking in an adiabatic spin-glass model of adaptive evolution, EPL (Europhysics Letters) 99, 68004 (2012).
  * Saito  _et al._ [1997] S. Saito, M. Sasai, and T. Yomo, Evolution of the folding ability of proteins through functional selection, Proceedings of the National Academy of Sciences 94, 11324 (1997).
  * Swain and Gierasch [2006] J. F. Swain and L. M. Gierasch, The changing landscape of protein allostery, Current opinion in structural biology 16, 102 (2006).
  * Kurikawa and Kaneko [2021] T. Kurikawa and K. Kaneko, Multiple-timescale neural networks: generation of history-dependent sequences and inference through autonomous bifurcations, Frontiers in Computational Neuroscience 15 (2021).
  * Pearl [1988] J. Pearl,  _Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference_ (Morgan Kaufmann, 1988).
  * Koller and Friedman [2009] D. Koller and N. Friedman,  _Probabilistic Graphical Models: Principles and Techniques_ (The MIT Press, 2009).
  * Mezard and Montanari [2009] M. Mezard and A. Montanari,  _Information, physics, and computation_ (Oxford University Press, 2009).
  * De Almeida and Thouless [1978] J. De Almeida and D. J. Thouless, Stability of the sherrington-kirkpatrick solution of a spin glass model, Journal of Physics A: Mathematical and General 11, 983 (1978).
  * Kabashima [2003] Y. Kabashima, A cdma multiuser detection algorithm on the basis of belief propagation, Journal of Physics A: Mathematical and General 36, 11111 (2003).
  * Note [1] While the AT instability cannot be derived from not fully connected systems, including our evolutionary model, the BP algorithm is often used as a numerical method to identify the phase of not fully connected systems.
  * Sato and Kaneko [2020] T. U. Sato and K. Kaneko, Evolutionary dimension reduction in phenotypic space, Phys. Rev. Res. 2, 013197 (2020).
  * Tang and Kaneko [2021] Q.-Y. Tang and K. Kaneko, Dynamics-evolution correspondence in protein structures, Phys. Rev. Lett. 127, 098103 (2021).
  * Togashi and Mikhailov [2007] Y. Togashi and A. S. Mikhailov, Nonlinear relaxation dynamics in elastic networks and design principles of molecular machines, Proceedings of the National Academy of Sciences 104, 8697 (2007).
  * Raman  _et al._ [2014] S. Raman, N. Taylor, N. Genuth, S. Fields, and G. M. Church, Engineering allostery, Trends in Genetics 30, 521 (2014).
  * Bornscheuer and Pohl [2001] U. T. Bornscheuer and M. Pohl, Improved biocatalysts by directed evolution and rational protein design, Current opinion in chemical biology 5, 137 (2001).
  * Bahar  _et al._ [2017] I. Bahar, R. L. Jernigan, and K. A. Dill,  _Protein actions: Principles and modeling_ (Garland Science, 2017).
  * Hopfield [1982] J. J. Hopfield, Neural networks and physical systems with emergent collective computational abilities., Proceedings of the national academy of sciences 79, 2554 (1982).
  * Amit  _et al._ [1985] D. J. Amit, H. Gutfreund, and H. Sompolinsky, Spin-glass models of neural networks, Physical Review A 32, 1007 (1985).
  * Li  _et al._ [2019] W. Li, J. Wang, J. Zhang, S. Takada, and W. Wang, Overcoming the bottleneck of the enzymatic cycle by steric frustration, Phys. Rev. Lett. 122, 238102 (2019).
  * Wootten  _et al._ [2013] D. Wootten, A. Christopoulos, and P. M. Sexton, Emerging paradigms in gpcr allostery: implications for drug discovery, Nature reviews Drug discovery 12, 630 (2013).
  * Lu  _et al._ [2017] J. Lu, N. Byrne, J. Wang, G. Bricogne, F. K. Brown, H. R. Chobanian, S. L. Colletti, J. Di Salvo, B. Thomas-Fowlkes, Y. Guo, _et al._ , Structural basis for the cooperative allosteric activation of the free fatty acid receptor gpr40, Nature structural & molecular biology 24, 570 (2017).
  * Frank  _et al._ [2004] R. A. Frank, C. M. Titman, J. V. Pratap, B. F. Luisi, and R. N. Perham, A molecular switch and proton wire synchronize the active sites in thiamine enzymes, Science 306, 872 (2004).
  * Hemery and Rivoire [2015] M. Hemery and O. Rivoire, Evolution of sparsity and modularity in a model of protein allostery, Physical review E 91, 042704 (2015).
  * Tlusty  _et al._ [2017] T. Tlusty, A. Libchaber, and J. P. Eckmann, Physical model of the genotype-to-phenotype map of proteins, Phys. Rev. X 7, 021037 (2017).
  * Rivoire [2019] O. Rivoire, Parsimonious evolutionary scenario for the origin of allostery and coevolution patterns in proteins, Phys. Rev. E 100, 032411 (2019).