Transition in relaxation paths in allosteric molecules : Enzymatic kinetically constrained model

A hierarchy of timescales is ubiquitous in biological systems, where enzymatic reactions play an important role because they can hasten the relaxation to equilibrium. We introduced a statistical physics model of interacting spins that also incorporates enzymatic reactions to extend the classic model for allosteric regulation. Through Monte Carlo simulations, we found that the relaxation dynamics are much slower than the elementary reactions and are logarithmic in time with several plateaus, as is commonly observed for glasses. This is because of the kinetic constraints from the cooperativity via the competition for an enzyme, which has different affinity for molecules with different structures. Our model showed symmetry breaking in the relaxation trajectories that led to inherently kinetic transitions without any correspondence to the equilibrium state. In this paper, we discuss the relevance of these results for diverse responses in biology. ∗ hatakeyama@complex.c.u-tokyo.ac.jp 1 ar X iv :1 80 8. 05 03 1v 3 [ ph ys ic s. bi oph ] 2 0 Ju n 20 19

B iological systems are known to have a hierarchy of timescales (1). Ordinarily, the timescale of biochemical reactions is of the subsecond order, that of organisms' behaviors is of the order of seconds to hours, and that of lifespans is of the order of years to a century. How organisms fill the gaps between such timescales is one of the most important problems in biology and physics; mechanisms that fill the gap between the timescales of biochemical reactions and organisms' behaviors have yet to be unveiled.
A biochemical process often consists of an enzymatic reaction with a complex protein, which has multiple modification sites (2)(3)(4). The multisite modification of such molecules has recently been reported to be essential for regulating biological timescales (5)(6)(7)(8)(9). A prominent example is epigenetic regulation by histone and DNA modification; silencing gene expression by epigenetic modification takes several tens of hours, and the erasure of such memory takes a much longer time (5). Moreover, the timescale of silencing and erasure is broadly distributed among clonal cells.
The regulation of biological timescales by multisite modification has also been observed in signaling pathways. In an Erk/Elk-1 signaling pathway, the timescales of the phosphorylation reactions of multiple phosphorylation sites are broadly distributed (6). Mutants in which some phosphorylation sites are replaced by an amino acid mimicking the unphosphorylated state have been used to demonstrate that the phosphorylation speed of each site depends on both the site itself and the competition for Erk, which works for an enzyme, between multiple modification sites.
In these examples, the modification speed of proteins with multisite modification can be estimated neither from the turnover rate of an enzyme nor from the affinity between the substrate and the enzyme, as in ordinal Michaelis-Menten kinetics. Extensive theoretical studies have shown that multisite modification and the competition for an enzyme can change the kinetics as well as the steady-state modification level (7,9,10). Notably, sequential multisite modification has been reported to generate a variety of timescales, some of which are much slower than the enzymatic turnover rate (10). Such regulation of the modification kinetics, including the slow dynamics, is considered to result from intermolecular cooperativity. A dynamical-system model with chemical kinetics has previously been proposed for the average relaxation process. To reveal the intermolecular cooperativity and the fluctuations in the slow relaxation process, a model and analysis that go beyond dynamical systems are required.
Concepts from statistical physics may be useful for investigating the slow stochastic biochemical dynamics and fluctuation. Such slow dynamics have been extensively and intensively studied with regard to the physics of glasses (11). There are two theoretical phenomenological approaches for such glassy behavior. One is the thermodynamic approach, which has been studied for the spin-glass model (12). The energy landscape is rugged, and the state is trapped in many local energy minima; thus, relaxation to the global minimum . .

Significance Statement
In biological systems, slow dynamics with an orders-ofmagnitude difference in timescales are ubiquitous. In the physics of glasses, as has been extensively investigated, such slow dynamics are understood to be generated by a kinetic constraint with multibody interactions, even without a complex, rugged energy landscape. Here, we introduce the "enzymatic kinetically-constrained model" to bridge the two fields, which fits well with the relaxation dynamics in allosteric biomolecules. The relaxation dynamics therein are found to be drastically slowed-down with several plateaus. Moreover, the relaxation trajectories to equilibrium show the first-and second-order-like transitions against changes in the enzyme concentration and temperature, respectively, which will provide a novel concept in statistical physics and a basis for understanding diverse kinetic responses in cells.
T.S.H. and K.K. designed research; T.S.H. performed research; and T.S.H. and K.K. wrote the paper.
The authors have no conflicts of interest to declare.

D R A F T
is slowed down. The other is a kinetic constraint, in which no multiple thermodynamic metastable states exist in the energy landscape but relaxation to equilibrium is kinetically suppressed (13). The spin-glass model has sometimes been adopted for neural (14) and gene-regulation (15) networks. On the other hand, a promising alternative mechanism exists for kinetically slowed-down processes without resorting to the existence of multiple metastable states, namely, controlling the enzyme abundance. Because the reaction rate is controlled enzymatically, the lack of an enzyme may suppress the corresponding reactions. Despite the possible relevance of the kinetic constraint concept to biochemical processes, however, it has not been fully explored in biology.
In this study, to uncover a relationship between the slow dynamics in biology and the kinetic constraint, we adopted the Monod-Wyman-Changeux (MWC) model for multiple modifications of the protein state (Fig. 1). This is the classic established model for concerted allosteric regulation. The original MWC model was introduced in 1965 to explain the allosteric effect of hemoglobin against changes in the O 2 concentration (16). The model includes both the binding and unbinding of ligands and a large change in shape between the relaxed (R) and tense (T) states. When no ligands are bound, the molecule tends to be in the T state, which has a low affinity to the ligands. As ligand binding progresses, the molecule tends to be in the R state, which has higher affinity to the ligands. Thus, ligand binding is accelerated as the ligand concentration increases, which results in the allosteric effect. We call this nonlinear change in binding the (1) energetic effect, which has been considered in statistical-physics models (17).
Because a modification reaction (e.g., phosphorylation, methylation, and acetylation) is generally catalyzed by an enzyme, we extended the original MWC model to incorporate kinetic facilitation by the enzyme. When the modification progresses and a molecule transitions to the R state, binding with the enzyme is accelerated, and the reaction proceeds faster. We call this the (2) enzymatic (kinetic) effect; this effect has not been considered in previous studies.
We used statistical physics to analyze the model by mapping the state of the MWC molecule according to the 'spin." Each molecule has a single spin representing the structural change and multiple spins representing the modification of each site. The energy of each molecule depends on the conformation of spins, and a spin flip requires the enzyme, so that both the energetic and enzymatic effects are included. Although both effects accelerate reactions, the model counterintuitively showed a slow relaxation to the equilibrium state. A typical time course for relaxation showed multiple plateaus, where the modification progress was transiently frozen far from equilibrium. The relaxation was logarithmic in time rather than exponential. Furthermore, the relaxation time did not follow the standard Arrhenius form, an exponential of the inverse temperature. The origin of such slow relaxation is explained as the decrease in the kinetic constants of reactions, which is similar to the kinetically constrained model (KCM) for glass (13). In contrast to the KCMs studied thus far, however, the present enzymatic kinetically constrained model (eKCM) showed symmetry breaking of the trajectories to equilibrium. This suggests the existence of a phase transition in the relaxation trajectories rather than in the equilibrium thermodynamic quantities. We discuss the biological relevance of these results in this paper.

Model
The original MWC model has two properties: modification and structural change.
(1) Modification (or ligand attachment) at several sites. Each molecule has several sites (M ) that represent multiple modification states in a single protein.
When modification progresses, a protein changes its structure. If this structural change produces an allosteric effect, the protein shows a large deformation, and modification becomes more feasible. In the original MWC model, these two structures are called the T and R states (16).
As noted earlier, we represent each modification and each structural state as two types of spins; that is, each molecule has M modification spins and a single state spin taking a down or up state. Thus, each molecule has 2 M +1 states. The modification spin si flips from 0 to 1 when each site is modified, where i is 0 ≤ i ≤ M , and the total number of up spins is denoted as m. The structural spin σ flips from down to up with a structural change from the T (σ = 0) state to R (σ = 1) state. The internal energy of a single molecule is defined as the summation of all modification spins. Thus, the microscopic Hamiltonian of the single molecule is given as where ǫ 1 σ (h) and ǫ 0 σ are the energies of the up and down spins, respectively, of the σ state protein. ǫ 1 σ (h) is a function of h, which is the "chemical field" derived from the concentration of the coenzyme required for the transfer of functional groups, and it is set as constant without losing generality. {si} is a set of modification spins. The numbers of molecules and enzymes in the present system are fixed to N and n, respectively. Therefore, the partition function is given as

D R A F T
Corresponding to the magnetization in the ordinal spin models, we introduce two quantitative measures: the fraction of unmodified monomers U and the T-state molecule ratio T . These are defined as 1 − ΣN m/(N M ) and N0/N = 1 − N1/N , respectively, where N0 and N1 are the numbers of T-state σ = 0 and R-state σ = 1 molecules, respectively. Such quantities in the equilibrium state are easily derived from the partition function because there is no interaction term among the molecules in the Hamiltonian, whereas we need to specify the kinetics to study the relaxation process. When investigating the kinetics, the allosteric effects cannot be neglected. When modification progresses, a molecule tends to flip its state, and the modification reaction is accelerated both energetically and enzymatically. Although the enzyme works as a catalyst for modification and does not change the detailed balance condition, competition for the enzyme among molecules introduces a kinetic effect rather than an energetic effect. By assuming that the timescale of enzyme binding is much faster than that of modification and state change, the binding reaction can be eliminated adiabatically. Therefore, the kinetics of the modification-spin flip are governed by the product of the binding probability (P b σ ) and the probability to traverse across the energy barrier for activation, whereas the kinetics of the state-spin flip are governed only by the latter. The energetic effect is represented by the microscopic energy of each state, whereas the enzymatic effect is represented by the changes in P b 0 and P b 1 following the state-spin flip. The two types of allosteric effects are formulated as follows. (

1) Energetic effect
If a molecule is modified at many modification sites, such a protein tends to change from the R state to the T state. The modification sites of R-state proteins are easier to modify than those of T-state proteins. Therefore, the energy of each modification spin has to satisfy the following inequality: For simplicity, we set the difference of each energy as unity. By considering the detailed balance condition, the transition probability of each protein state is given as , · · · , sM }|σ, {0, · · · , sM }) and p(σ, {0, · · · , sM }|σ, {1, · · · , sM }) are the transition probabilities for the modification and unmodification of the 1st modification spin of a σ-state molecule, respectively. The same transition probabilities are adopted for the other modification spins. p (1, {si}|0, {si}) and p(0, {si}|1, {si}) are the transition probabilities for the flipping of a molecule with modification state {si} from σ = 0 to σ = 1 and σ = 1 to σ = 0, respectively.
Owing to the allosteric effect, the state tends to change to up as modification progresses; therefore, the modification spins tend to be up. When m is small, the state spin tends to be σ = 0 because ǫ 0 1 > ǫ 0 0 , whereas the state spin tends to be σ = 0 because ǫ 1 0 > ǫ 1 1 for large m. The transition probabilities from σ = 0 to 1 and σ = 1 to 0 are identical for m = M/2. Here, we assume that the state-flip reaction can always occur within the characteristic time of the state flip when the microscopic energy decreases. Therefore, we set the maximum value of the flip probability of the state spin as unity under the detailed balance conditions given in Eq. (5).
For the modification-spin flip, the activation energy is set to be equal to ǫ 0 1 for all modification reactions for simplicity. The modification reaction of the R-state molecule has no energy barrier and can always occur with the characteristic timescale of the modification-spin flip, whereas that of the T-state molecule has an energy barrier.
The state change between R and T increases the affinity with the enzyme; that is, it increases the probability of binding with the enzyme. As modification progresses, the fraction of R-state molecules increases, which tends to bind the enzyme more and is modified faster than the T state. We represent such an enzymatic allosteric effect by using equilibrium statistical physics. We assume that only a single enzyme can bind to the molecule. The T-and R-state molecules have different binding energies with the enzyme of ǫ b 0 and ǫ b 1 , respectively, where ǫ b 0 is lower than ǫ b 1 for the positive allosteric effect. Thus, the transition probability for the modification reaction under the detailed balance condition (Eq. (4)) (see the Supplemental Information for a derivation): , [8] where nσ is the average number of enzymes that bind to the σ-state molecule. Fig. S1 shows the binding probabilities calculated thus far. We set the timescales of the modification flip and state flip as τm and τs, respectively, where τm is longer than τs. We set the initial condition such that all molecules are in the σ = 0 and m = 0 states, and we investigated the relaxation dynamics to the equilibrium state using the Monte Carlo method.

Results
eKCM shows a slow relaxation to equilibrium and multiple plateaus. First, we calculated the relaxation dynamics of U and T . Although < U >ens and < T >ens, where <>ens is the ensemble average, finally relaxed to the equilibrium values Ueq and Teq, respectively, their time courses varied depending on the temperature 1/β and < n >, which is the number of enzymes binding to a substrate for a fixed number of molecules N (Fig. 2). When the temperature was high, the relaxation of < U >ens decreased exponentially with time with no plateau. As the temperature decreased, the relaxation slowed down and decreased logarithmically with time. Two plateaus appeared as the temperature decreased further (see Fig. 2a). The two plateaus were clearly discernible when < n > was reduced below < n > /N = 0.95.
The value of < U >ens at the plateaus increased with decreasing < n > (see Figs. 2b and c). The relaxation of < T >ens also showed a similar dependence on < n > and the temperature, as shown in Fig. S2. Such slow dynamics with plateaus have often been observed in glasses (11). Anomalous parameter dependence on the relaxation time.

D R A F T
Glasses are known to exhibit a nonstandard temperature dependence on the relaxation time that is distinct from the normal exp(Eactβ) form derived by the standard thermal activation process with the activation energy Eact. To characterize the glassy phase, we computed the dependencies of the relaxation time for equilibration on < n > and the temperature. When the temperature was fixed and < n > was varied, the transition probabilities for the number of modified spins changed in proportion with the number of the complex < n > in accordance with the binding probability. Indeed, for lower < n >, the relaxation time was proportional to the inverse of < n >. However, it was further prolonged beyond < n > −1 as < n > decreased to approach 1 (Fig. 3a). Such an < n > dependence cannot be explained by the simple < n > dependence of the binding probability of the enzyme.
The relaxation time also showed an anomalous temperature dependence. Because the energy barrier for the modification of the T-state molecule was set to unity while the modification of the R-state molecule was temperature-independent, the relaxation time would normally be expected to be proportional to exp(β). Indeed, at high temperatures, the relaxation time approximately followed exp(β). However, as the temperature decreased, the rate of increase in the relaxation time against the temperature was enhanced and reached its maximum at around β = 1. The relaxation time showed a bending point around β ∼ 1, and for low temperatures, it approximately obeyed exp(6β) (see Fig. 3b). The fact that the activation energy was six times the energy barrier for each modification (see Fig. 1) suggests that the transition from the (m = 0, σ = 0) state to the (m = 0, σ = 1) state is rate-limiting.
The two salient features, namely, the plateaus and anomalous parameter sensitivity of the relaxation time, suggest that the present model has characteristics similar to those of glasses.

Symmetry breaking of the relaxation trajectories in the eKCM.
To reveal the mechanism of the anomalous parameter dependence, we analyzed the relaxation-time distribution of the samples. In the region where the relaxation time showed anomalous parameter dependence, the relaxation-time distribution changed from unimodal to multimodal (Fig. 4). The multiple peaks that emerged are named the first, second, and third peaks in ascending order of relaxation time.
When < n > was varied at a fixed temperature, the positions of the peaks changed in proportion to < n > −1 (see Fig. 4a). Thus, we studied the distribution of the relaxation time normalized by < n > −1 . When < n > /N was close to 1, two peaks were observed. Then, as < n > decreased, the first peak disappeared and was replaced by the third one. Finally, the second peak disappeared completely at < n > /N = 0.4. This change in the distribution is similar to the first-order phase transition in equilibrium thermodynamics. The transition here, however, was with regard to the relaxation time rather than the thermodynamic quantities. It should be recalled that the same equilibrium state was reached over all samples independent of the relaxation courses. This multimodal distribution of the relaxation time was concerned with the symmetry breaking of the trajectories and not with the equilibrium state.
The temperature dependence of the relaxation-time distribution when < n > /N was fixed close to 1 differed from its < n > dependence (see Fig. 4b). In this case, the position of the first peak changed in proportion to exp(4β), as explained later. At high temperatures, the relaxation-time distribution was unimodal. This distribution broadened with decreasing temperature at around β = 1.0, where it became bimodal, and the distance between these two peaks increased with decreasing temperature. This change in the relaxation-time distribution is similar to the second-order phase transition in equilibrium.
As shown previously, the average relaxation time was proportional to exp(6β) for low temperatures (Fig. 3b). In fact, the second peak changed with exp(6β), and the position of the first peak followed exp(4β). This suggests that the samples showing different peaks experienced different rate-limiting steps during relaxation. As described before, the exp(6β) dependency was due to the transition from (m = 0, σ = 0) to
As the relaxation progressed and the enzyme was occupied by the R-state molecules from the positive allosteric effect, it became more difficult for modification reactions of the T-state molecules to occur (see Fig. 5a), which provided a kinetically constrained condition. Thus, the transition from the T state to the R state was rate-limiting. Under the present conditions, three possible steps are candidates for the rate-limiting transition: from (m = 0, σ = 0) to (m = 0, σ = 1), from (m = 1, σ = 0) to (m = 1, σ = 1), and from (m = 3, σ = 0) to (m = 3, σ = 1).
The step that is rate-limiting depends on the distribution of modifications at the time when the R-state molecule increases and occupies the enzyme. In the parameter region close to the transition point, some samples included m = 0 state molecules, whereas others did not. Therefore, there were two groups with regard to the relaxation time: one following exp(6β) and the other following exp(4β). Indeed, only the samples included in the larger peak value had m = 0 state molecules and showed three plateaus (Fig. S3). Note that the trimodal distribution shown in Fig. 4a was due to the transition from (m = 3, σ = 0) to (m = 3, σ = 1). With this mechanism, symmetry breaking appeared in the relaxation trajectory.
In summary, if m = 2-, 1-, and 0-state molecules coexist at the start of the kinetically constrained condition, they relax following exp(2β), exp(4β), and exp(6β), respectively. Based on the above considerations, the multimodal distribution of the logarithmic relaxation time is an important indicator of the transition to the glass phase. Indeed, the variance of the logarithmic relaxation time attained a high value at the upper boundary of the glass phase (Fig. 4d). It is noted that when we changed the temperature with a fixed small < n > value, the relaxation-time distribution did not show a transition from unimodal to bimodal (Fig. 4c). This suggests that such a change is similar to crossover rather than the transition (see also Fig. 5b). In fact, there was no change in the appearance of the plateaus at low < n > values, as shown in Fig. 2a.

Discussion
We propose a statistical physics model for allosteric biomolecules and enzymes that catalyze the modification reaction of the allosteric molecules. The present model showed a slow relaxation to equilibrium with several plateaus, in which the relaxation time depended anomalously on the temperature. This behavior is known in glass theory. In glasses, such slow dynamics are considered to be caused by the cooperativity among molecules. In our model, there are two types of cooperativity. One is the intramolecular cooperativity, which is caused by the cooperative flip of the molecular structure of all M monomers. The maximum value of the variance among samples approaches M , reflecting such intramolecular cooperativity (see Fig. S4). The other is intermolecular cooperativity, which is generated by the interaction of allosteric molecules with each other via the enzyme. The slow relaxation in the present study originates not from a complex energy landscape but from a kinetic constraint by enzyme limitation. In our model, allosteric molecules interact with each other kinetically without energetic interaction, as shown in Eq. (1) in the "Methods" section. Such an intermolecular cooperativity in kinetics is similar to that studied in the KCM for glasses. Our eKCM showed dynamical heterogeneity, which is important in glass theory (18). In the slow relaxation of the eKCM, the relaxation of modification states of the T-state molecules was frozen, whereas R-state molecules were partially equilibrated. The number of frozen molecules depended on the history and differed for each sample (Figs. S4 and S5).
The eKCM demonstrated the symmetry breaking of the relaxation paths to equilibrium depending on the temperature. Moreover, the symmetry breaking of paths was also observed when the concentration of the enzyme was close to that of the allosteric molecule. Following such symmetry breaking in trajectories, slow relaxation and plateaus appeared, similar to the glass transition. In contrast, most of the KCMs studied D R A F T thus far do not show the phase transition of paths against the change in temperature (13). In the eKCM, the strength of the cooperativity in kinetics depends on the enzyme abundance and temperature, unlike in the case of many KCMs. When the amount of enzyme is greater than that of the substrate, there is no competition, and glassy behavior does not appear. The competition for a limited amount of the enzyme leads to interactions among molecules and a transition in the paths.
The symmetry breaking discovered in this study is considered to be a new nonequilibrium phase transition, in contrast to those concerned with the steady state, as previously reported (19)(20)(21). Here, we found a phase transition in the relaxation process to equilibrium. Although such a phase transition has been discussed in the KCMs of glasses, such discussions have been limited to the transition against the "activity" (22), which is a fictitious parameter introduced for large-deviation analysis and not a real parameter that can be controlled, as in our model. We expect that further analysis of the eKCM will lead to deeper understanding of the nonequilibrium phase transition with glass behavior.
In recent decades, the slow relaxation process has also attracted much interest in the field of biology. In bacterial chemotaxis, chemoreceptors, which are often described by an MWC-type model (23,24), are known to form clusters with each other. A recent experiment demonstrated that such chemoreceptor clusters show a slow logarithmic change in their structure with time in response to the addition and removal of a ligand (25). This slow relaxation is a result of the interactions among the chemoreceptors that modify their structure without using a specific adaptation by signaling enzymes (cheR and cheB). Therefore, the model in our study fits well with this problem of relaxation of activated chemoreceptors. A plateau has not yet been observed in such experiments, although such logarithmic relaxation without plateaus was also observed in our model around the critical temperature. To find further correspondence between theory and experiments on the chemoreceptor cluster, a two-dimensional space should be considered in the eKCM.
Another biological example that corresponds with our theory is Ca 2+ /calmodulin-dependent protein kinase II (CaMKII), which is known to be important for maintaining long-term potentiation in a synapse (26). Here, the structural changes in CaMKII that accompany phosphorylation are sustained over a few minutes against a stimulus in the hippocampal synapse (27). On the other hand, a recent study has suggested that the phosphorylation of CaMKII in the dorsal raphe nucleus regulates the sleep-awake cycle, particularly the length of sleeping hours (i.e., typically 6 h) (28). These results suggest that the phosphorylation of CaMKII can show various timescales during relaxation ranging from minutes to hours, which was the focus of our study.
Our study also demonstrated that a small fluctuation in molecules can be amplified to a large variation in the relaxation time. This provides a way for the fluctuations in the relaxation time in biology to be studied. Although the noise in chemical concentrations has recently attracted much attention from many physicists and biologists (29), there has been little study on the fluctuations in the relaxation paths in biology. Our results suggest that cells with the same genetic and environmental backgrounds can have quite different relaxation times and different fates. Further theoretical and experimental research will unveil the biological significance of the relaxation-time fluctuations.
Possible candidates for confirming the glassy relaxation and symmetry breaking in the relaxation courses are measuring the time courses of the modification levels, excitation energy, and concentration of enzymes. The critical temperature for the transition in the relaxation behavior that we studied depends on the difference in the energy of each state of a monomer in our model. This difference can be altered by external conditions such as the concentration of a coenzyme, e.g., ATP as a phosphate donor for phosphorylation, or a mutation to mimic a (un)modified form. We expect that experimental studies on the relaxation process along with theoretical studies will deepen our understanding of the relationship between the regulation of biological timescales and glass theory in physics. We assume that only a single enzyme can bind to the molecule. The Tand R-state molecules have different binding energies with the enzyme of ǫ b 0 and ǫ b 1 , respectively, where ǫ b 0 is lower than ǫ b 1 for the positive allosteric effect. Therefore, the grand partition function and the average number of binding enzymes < n > in equilibrium are represented in a similar manner as Langmuir's adsorption isotherm [1] and are given as

D R A F T
For both in vivo and in vitro reactions, N and n are almost constant and are the control parameters of the system. Here, we assume that < n > is the same as n and that the chemical potential µ is a function of < n >, N 0 , and N 1 = N − N 0 . This assumption is justified when N and n are sufficiently large, i.e., at the thermodynamic limit. Then, exp(−βµ) is given as where A = (< n > −N + N 0 ) exp(−βǫ b 0 ) + (< n > −N 1 ) exp(−βǫ b 1 ). In the limit < n >→ N , µ approaches −∞, and all molecules bind to the enzyme. The arguments so far can be summarized to give the binding probability where n σ is the average number of enzymes that bind to the σ-state molecule. It is noted that in the thermodynamic limit, the binding reaction is described by the Michaelis-Menten equation for the timescale separation. The Michaelis-Menten equation for multiple substrates with the conservation of the enzyme concentration is given as gradually with decreasing temperature (see Fig. S4a). This corresponded to the gradual appearance of the plateau in the averaged relaxation dynamics (Fig. 2a in the main text). As < n > decreased, the change in the variance became more drastic. When < n > /N was 1, there was only one peak at an early step (see the purple line in Fig. S4b). At < n > /N = 0.9, a second plateau appeared suddenly in the average relaxation dynamics (see Fig. 2b in the main text). Accordingly, a second peak in the variance appeared suddenly (see the blue line in Fig. S4b). The variance increased at the beginning of the plateaus because the number of modified spins varied with each sample, and the induced variance froze at the plateaus (see Fig.  S5). Therefore, the variance among samples provided a good indicator of the plateaus.
To clarify the region where multiple plateaus appear, we plotted the maximum variance after the first peak by subtracting the equilibrium value of the variance (Fig. S6). The diagram shows that the variance was high under low-< n > and low-temperature conditions, which indicates the existence of multiple plateaus similar to glassy dynamics. Accordingly, we call this region the "glass phase."   This was calculated as the maximum variance among the ensemble in the time course after 20 MCSs. It was normalized by dividing it by the average < U > ens and was subtracted from the analytically calculated value < (U − U eq ) 2 > eq /U eq .