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

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. How are robustness and plasticity simultaneously acquired through evolution? We examine the evolution of genotypes that realize plastic switching between two phenotypes upon external inputs, as well as stationary expressions of phenotypes. We introduce a statistical physics model consisting of spins, with active and regulatory sites, which are distinct from each other. We represent the phenotype and genotype as spin configurations and the spin-spin interactions, respectively. 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 regulation. We numerically evolve the interaction matrix by changing them with mutations and selection of those with higher fitness. Our numerical simulations show that characteristic genotypes evolve slightly above the transition temperature between replica symmetric and replica symmetry breaking phase. 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, and smooth switching of two phenotypes are achieved by following a one-dimensional 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 show that the free energy landscape has a valley along the switching path. Our finding indicates that the compatibility of the robustness and plasticity is acquired by the evolution of the low-dimensionality in the phenotype space.


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, molecularlevel 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 confor-mational 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 S and their interaction variables 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 replicasymmetric 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
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  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 S = {S 1 , • • • , S N } ∈ {−1, +1} N which represents the conformational change in each amino-acid residue, whereas the coupling J between spins represents the interaction among residues.These are respectively denoted by nodes and edges in Fig. 1(b).Here, we set J as N × N symmetric matrix.
Its elements are given by J ii = 0 (i = 1, • • • , N ), and J ij ∈ Ω J for i = j where Ω J = {−1/ √ N , 0, 1/ √ 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 A and R, respectively.We set J ij = 0 for i ∈ A and j ∈ R or i ∈ R and j ∈ 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 J , we adopt the transition rule of spins from S to S under given Hamiltonian with the interaction matrix J and temperature T as where β = T −1 is the inverse of temperature, and ∆ H (S, S |J ) ≡ H(S |J ) − H(S|J ).Here, we set the Hamiltonian as where 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 genotypephenotype mapping, in which spin variables S corresponding to the phenotypes are shaped by high-dimensional dynamics under genetic rules given by the interaction matrix 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 S A = {S i |i ∈ A}, depending on the configurations of the regulatory site S R = {S i |i ∈ 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 S + R upon the binding.That is, among 2 N R possible configurations of the regulatory spins, the regulatory spins only take the configuration in S + R when the ligand binding occurs.Further, we consider that S + R cannot appear without the ligand binding.Accordingly, the equilibrium distribution upon the regulation and non-regulation is given by where S|S R ∈ S + R and S|S R / ∈ S + 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 S − A to S + A : Thus, the conformational change induced by regulation is modeled as follows: if the configuration of regulatory spins is set at S + R , the configuration of the active spins turns into S + A ; otherwise, the configuration of the active sites stays at S − A .The function of the present system to express the target spin pattern S ± A appropriately can be measured by the magnetization m ± A defined as the overlap of the spins in the active sites with the corresponding target spin patterns as Finally, the overall fitness that measures the functionality of the present system is given by the sum of the expectations of m ± A as where • + and • − are the expectation values according to the equilibrium distributions for regulated and non-regulated cases eq.( 3) and eq.( 4), The evolution of genotypes J is then based on the above fitness ψ(J ).Higher fitness genotypes are selected under given selective pressure: At generation g, the evolutionary change in J to increase the fitness is given by where ∆ψ = ψ(J (g+1) ) − ψ(J (g) ).The parameter β J = T −1 J expresses the selection pressure, and the genotypes are selected uniformly at high temperatures T J → ∞, 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 R = {N −N R +1, . . ., N } and A = {1, . . ., N A }, respecively.Further, we set the configuration For the desirable configurations of the active sites, we set In the genotype evolution process, we induce a 10-point mutation at each generation to generate the candidate of the next generation J from J , maintaining the symmetry J = J .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 J at a sufficiently large value as β J = 100, and discuss T -dependencies.

A. Fitness, rugged landscape, and separation of two patterns
In Fig. 2, we show examples of evolutionary dynamics of J through the evolutionary ± , where • (g) denotes the expectation according to the distribution P ± β (S|J (g) ) with a g-th generation genotype 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).and vice versa.For the genotypes that show this behavior, which are the most evolved genotypes around T = 0.91, the simultaneous expression of S + A and S − A , depending on the regulatory sites, is difficult.When the active sites take one of the configurations in S + A or S − A , irrespective of the regulatory sites, we obtain |m

and |m −
A | = 1, respectively.Therefore, the fitness value of the genotype that can express only one desirable pattern among S ± A can reach a highest value of 0.6.Meanwhile, at a lower temperature ± increase simultaneously after the evolution, as shown in Fig. 2(b)-(d), where the fitness reach around 0.9.There are three evolutionary courses; − 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.
For each T , we obtain 100 samples of the evolved J updated for g = 10 5 generations, and denote the temperature-dependence ensemble of evolved genotypes as J (T ).Fig. 3 (a) shows the T -dependence of the mean of the fitness among J (T ).With a decrease in T , the fitness value Ψ increases from 0.4, which is a trivial value given by the uniform distribution of the phenotype 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 Ψ shows discontinuity.We term the phase T > T 0 as the paramagnetic phase.Next, the energy landscape on J ∈ 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 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] [43].At T > T 2 , most of the evolved genotypes in J (T ) exhibit rapid convergence of the BP algorithm.
Meanwhile, the BP algorithm cannot converge for most evolved genotypes in 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 [31][32][33][34].In the present model, however, another transition appeared at T 1 , with respect to the achievability of two patterns.In Fig. 3(c T 1 > T > T 2 as RS1 and RS2, respectively.Here, we note a negative correlation between ± is observed in the paramagnetic and RS1 phases, as shown in Fig. 2(a).In the RS2 and RSB phases, the increase of both m ± A (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 + β (S|J ) and p − β (S|J ), with and without regulations, by means of the component- wise expected phenotype for each i = 1, . . ., N defined as where the term sign( N A i=1 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 N i=1 µ + i µ − 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.

B. Two-dimensional structure in the phenotype space
We investigate how the two patterns shaped by evolution are separated in the RS2 phase.
To compare N -dimensional mean phenotypes µ + under and µ − 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 J is decomposed as where ξ i and λ i are i-th eigenvector and i-th eivenvalue.We set the indices of the eigenmodes to be In Fig. 4, we show the scatter plots between µ ± i against eigenvectors ξ 1i and ξ 2i for i = 1, . . ., N at (a) T = 0.833 (RS1 phase) and (b) T = 0.667 (RS2 phase) under one realization of J ∈ J (T ).In the RS1 phase, the mean phenotypes µ ± , in particular µ − , are highly correlated with ξ 1 , as described by y = tanh(β √ N x) (see Fig. 4(a)).Here, the function tanh is consistent with the mean-field form of the magnetization µ ± i = tanh(β j =i J ij µ ± j ).Meanwhile, in the RS2 phase, the regulated µ + and non-regulated µ − states exhibit correlations with ξ 1 and ξ 2 , respectively, as shown in Fig. 4(b).In both phases, the correlations between µ ± and ξ r (r ≥ 3) are negligible.
In Fig. 5 (a) and (b), we show the temperature dependence of the correlation between the eigenvectors and µ ± by introducing the correlation coefficient between {ξ ri } and {atanh(µ ± i )} for r = 1, 2, 3. Here, the function atanh is introduced by considering the tanh-form dependencies of µ ± on ξ 1 or ξ 2 , as shown in Fig. 4. We denote the vector consisting of As shown in Fig. 5 (a), the correlation coefficient between the first eigenvector ξ 1 and the mean phenotype with regulation µ + 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 µ − 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 ξ 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 ± β (S|J ) are almost orthogonal to each other.The contributions of the first and second eigenmodes to the mean phenotypes µ ± are given by the magnitudes of their corresponding eigenvalues.In Fig. 5(c independently and identically obey the uniform distribution over Ω J .For T < T 0 , the first eigenvalue shows distinct increases from the expected value, whereas for T 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.
Following these observations, we map the mean phenotypes with and without regulation onto the two-dimensional space spanned by the first and second eigenvectors, ξ 1 and ξ 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 and •, 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 J that gives overlapped and separable phenotypes as type J1 and type J2, respectively.C. 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 J (T ) under a trial temperature T tr .The evolutionary process in our model selected genotypes among possible J s; hence, J ∈ J (T ) can be a candidate for genotypes in J (T tr ) (T = T tr ), in principle.By evaluating the fitness of J ∈ J (T ) at a different temperature T tr , we discuss the reason why J ∈ J (T ) cannot be selected at different temperatures.
Fig. 7 shows the T tr dependence of the fitness Ψ on J ∈ J (T ) for T = 0.91 (type J1; RS1) and T = 0.63 (type J2; RS2).At sufficiently large T 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 J (T ) can be regarded as a phase transition with respect to genotypes between the type J1 and J2.

D. 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 m − A increased before m + A .The panels of Fig. 8 show the time evolution of the mean phenotypes µ + (J (g) ) ( ) and µ − (J (g) ) ( ) 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 |m +
A | + increases before the increase of |m + A | − , the localization of µ + on the second eigenvector appears in the early stage of evolution.Additionally, the localization of µ − follows with the reorganization of µ + .When both |m ± A | ± increased simultaneously, µ ± are localized almost simultaneously (see supplement material).

IV. SWITCHING TRAJECTORY
Under the 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 ( ), and from the regulated to non-regulated states ( ), 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 (•), second ( ), and third (dashed line) eigenvectors defined on an evolved 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 J ∈ 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, µ + and µ − , after sufficient time steps of updating are denoted by and •, 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 FIG. 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.and • 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.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.

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 Following the result of the numerical simulations, we consider the two-rank approximation of the evolved J as J λ 1 ξ 1 ξ 1 + λ 2 ξ 2 ξ 2 .Under the two-rank approximation, the Hamiltonian is given by For the two-rank approximation form, one can represent the free energy as a function of m 1 and m 2 defined by 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 where The saddle point equations for m ± k are given by Fig. 12 (a) and (b) show the landscape of f + and f − , respectively, plotted on the twodimensional 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 and •; 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.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 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 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 app (assumption A3).A particular difference between f ± and f app is that the valleys around the mean of the phenotype with and without regulation (Fig. 12) cannot be described by f 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 app would be invalid.
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,31,44,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 geno-type, 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 to the components of the eigenvectors can be replaced with the integral according to the Gaussian distribution.as By introducing the saddle point method to the integrals of m 1 and m 2 , we obtaine the approximated free energy f app .

FIG. 1 .
FIG. 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.

FIG. 2 .
FIG. 2. Evolutional dynamics of m ±A ± associated with the evolution of a genotype J .An example at T = 0.91 is shown in (a), and examples at T = 0.67 are shown in (b)-(d).

8 TFIG. 3 .
FIG. 3. T -dependence of (a) fitness, (b) the fraction of genotypes on which the BP algorithm does not converge, (c) m + A + and m − A − , and (d) Similarity between regulated and non-regulated states.Each data point is averaged over 100 samples of evolved J .The vertical dotted line, two-dot chain line, and one-dot chain line denote T 0 , T 1 , and T 2 , respectively.
), we show the temperature dependence of the overlaps |m + A | + and |m − A | − , whose mean corresponds to fitness.At T 1 < T < T 0 , |m + A | + contributes more to fitness, and only the phenotype expression with regulation is preferentially shaped.At T < T 1 , both the increase of |m − A | + and |m − A | − are achieved depending on the regulatory site.We term the phases T 0 > T > T 1 and

FIG. 5 .
FIG. 5.The correlation coefficient between the first, second, and third eigenvectors with (a) regulated state µ + and (b) non-regulated state µ − .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 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 J .

FIG. 6 .
FIG. 6.(a) Characteristic mapping of the mean phenotypes with regulation ( ) and without regulation (•), 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.

Fig. 6 (
Fig.6(b)  shows the temperature dependence of the fraction of the type J1 and J2 genotypes among the ensemble of evolved genotypes 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 J s, as there are two other cases of mapping; the case that and • located along ξ 2 , and that and • are along ξ 1 and ξ 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.

FIG. 7 .
FIG. 7. Trial temperature T tr dependence of the fitness for J ∈ J (T ) at T = 0.91 (denoted by J1) and T = 0.63 (denoted by J2).The shaded region indicates the difference between m − A − and m + A + .

8 FIG. 8 .
FIG. 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], • • • , [981 − 1100].The results with and without regulations are plotted by • and , respectively.

8 FIG. 9 .
FIG. 9. (a) MC steps required for the switching from regulated to non-regulated state ( ) and nonregulated to regulated state ( ).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 S under the non-regulated condition S R / ∈ S + R , the regulatory sites are changed to S R ∈ S + R , and then S (except the regulatory region) is updated according to(1).We compute the target magnetization| i∈A S i /N A |at each MC step to obtain the step where | i∈A S i /N A | first reaches the value |m + A | + , 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.

5 FIG. 11 .
FIG.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., •, and the orthogonal lines below these points are the same as in Figs.10.

3 FIG. 13 .
FIG. 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.