Extraction of kinetics from equilibrium distributions of states using the Metropolis Monte Carlo method

Sergei F. Chekmarev Institute of Thermophysics, SB RAS, 630090 Novosibirsk, Russia Department of Physics, Novosibirsk State University, 630090 Novosibirsk, Russia [email protected]

Abstract

The Metropolis Monte Carlo (MC) method is used to extract reaction kinetics from a given equilibrium distribution of states of a complex system. The approach is illustrated by the folding/unfolding reaction for two proteins - a model $\beta$ -hairpin and $\alpha$ -helical protein $\alpha_{3}$ D. For $\beta$ -hairpin, the free energy surfaces (FESs) and free energy profiles (FEPs) are employed as the equilibrium distributions of states, playing a role of the potentials of mean force to determine the acceptance probabilities of new states in the MC simulations. Based on the FESs and PESs for a set of temperatures that were simulated with the molecular dynamics (MD) method, the MC simulations are performed to extract folding/unfolding rates. It has been found that the rate constants and first-passage time (FPT) distributions obtained in the MC simulations change with temperature in good agreement with those from the MD simulations. For $\alpha_{3}$ D protein, whose equilibrium folding/unfolding was studied with the single-molecule FRET method (Chung et al., J. Phys. Chem. A, 115, 2011, 3642), the FRET-efficiency histograms at different denaturant concentrations were used as the equilibrium distribution of protein states. It has been found that the rate constants for folding and unfolding obtained in the MC simulations change with denaturant concentration in reasonable agreement with the constants that were extracted from the photon trajectories on the basis of theoretical models. The simulated FPT distributions are single-exponential, which is consistent with the assumption of two-state kinetics that was made in the theoretical models. The promising feature of the present approach is that it is based solely on the equilibrium distributions of states, without introducing any additional parameters to perform simulations, which suggests its applicability to other complex systems.

^†^†preprint: AIP/123-QED

I Introduction

The equilibrium distribution of states and the kinetics of their interconversion are two inherent characteristics of any complex system that are important to know for understanding its behavior. At the same time, the former are generally much easier to determine than the latter, both in experiment and in simulations. One representative example is the protein folding/unfolding reaction. Single-molecule experiments, such as the Förster resonance energy transfer (FRET), atomic-force microscopy (AFM) and optical tweezers (OT) techniques, make it possible to determine the equilibrium distributions of states, either in the form of the FRET-efficiency histograms (FRET schuler2002probing ; huang2009direct ; chung2011extracting ; liu2012exploring ; schuler2013single ; eaton2021modern , folding/unfolding) or free energy profiles (AFM yu2012energy ; batchelor2020protein and OTT bustamante2020single , unfolding). However, in order to extract kinetics from these distributions, theoretical models are needed (FRET mckinney2006analysis ; gopich2009decoding ; gopich2010fret ; ramanathan2015method ; osadko2016amplitude , and AFM and OTT hummer2001free ; hummer2005free ; hyeon2007measuring ; dudko2008theory ). A similar difference in the possibility of obtaining the equilibrium distribution of states and determining the kinetics exists for the molecular dynamics (MD) simulations of folding of relatively large proteins (of the order of one hundred of residues and larger). While brute force all-atom simulations, which are able to give a detailed insight into the folding kinetics, are very time-consuming shaw2010atomic ; lindorff-larsen2011how , several methods exist that allow very efficient calculation of the equilibrium distributions of protein states camilloni2018advanced , such as the replica exchange molecular dynamics sugita1999replica ; rhee2003multiplexed ; rao2003replica and the weighted histogram analysis method ferrenberg1989 in combination with sampling from biased molecular dynamics simulations shea2001from . Therefore, it is tempting to find out whether the knowledge of the equilibrium distribution of states of a system allows us to obtain information about its kinetics, and if so, to what extent.

The present paper examines the possibility of extracting kinetics from the equilibrium distributions of protein states with the use of the Metropolis’s Monte Carlo (MC) algorithm metropolis1953equation ; frenkel2004introduction , assuming that the acceptance probability of a new state can be determined by the equilibrium state distribution known from simulations or experiment. The typical examples of such distributions are free energy surfaces (FESs) and free energy profiles (FEPs), although, in general, any equilibrium distribution of states along one or several collective variables can be used for this purpose. In contrast to the nonequilibrium distributions of protein states that determine the probability fluxes of folding chekmarev2008hydrodynamic ; kalgin2014first , the equilibrium distributions represent potentials of mean-force (PMF) roux1995calculation ; trzesniak2007comparison ; shea2001from , and thus they can be assumed to govern kinetics. A known shortcoming of the MC methods is that a generated MC trajectory represents a Markovian sequence of states rather than the change of these states in time. Therefore, only the relative values of the reaction rates for different reaction conditions, e.g., for temperature, are available. In order to have absolute value of rates, a fit to the time-scale of the original data is required huitema1999can ; rutkai2010dynamic ; sanz2010dynamic ; bal2014time . At the same time, in contrast to other possible approaches which could employ the PMFs to calculate kinetics, e.g., the use of Fokker-Planck or Langevin equations kikuchi1991metropolis ; tiana2007use , no additional parameters are required to perform the MC simulations for a given equilibrium distribution of states. It is a considerable advantage of the present approach.

The approach is illustrated for two systems. One is a model $\beta$ -hairpin protein. In order to have baseline results, the equilibrium MD simulations were performed. The obtained FESs and FEPs were employed as the input distributions of states for subsequent MC simulations. It has been shown that the rate constants of folding and unfolding calculated from the FESs and FEPs change with temperature in good agreement with those from the MD simulations. The other system is an $\alpha$ -helical protein $\alpha_{3}$ D, for which folding and unfolding rate constants for different denaturant concentrations were extracted from single-molecule FRET trajectories chung2011extracting using theoretical models gopich2009decoding ; gopich2010fret . In this case, the FRET-efficiency histograms (FEHs) obtained in the experiment were employed as the equilibrium distributions of states for the MC simulations. It has been found that the folding and unfolding rate constants obtained in the MC simulations vary with the denaturant concentration in reasonable agreement with those extracted from the FRET trajectories. For both proteins, the first-passage time (FPT) distributions in the MC simulations have been found to be single-exponential, which agree with those obtained in the MD simulations ( $\beta$ -hairpin protein) and correspond to the assumption of two-state kinetics employed in the theoretical models to extract rate constants ( $\alpha_{3}$ D protein).

The paper is organized as follows. Section II describes the MD (Sect. II.1) and MC (Sect. II.2) simulations for the $\beta$ -hairpin protein, Sect. III present the results of the MC simulations for the $\alpha_{3}$ D protein, and Sect. IV briefly summarizes the results of the work.

II $\beta$ -hairpin Protein: System and Simulation

II.1 System and Molecular Dynamics Simulations

To perform MD simulations, a coarse-grained (C_α-bead) protein model and a Gō-type interaction potential go1983theoretical were used. Briefly, the approach is as follows (for details, see Ref. chekmarev2013protein, ). The C_α-model was constructed on the basis of the solution NMR data for 12-residue HP7 protein (2evq.pdb) andersen2006minimization . Specifically, the first conformation in the NMR ensemble of the protein structures was employed for this purpose, which was then considered as a reference native structure of the protein. The Gō-type potential accounted for the rigidity of the backbone and the contributions of native and non-native contacts hoang2000molecular . Two C_α-beads were considered to be in native contact if they were not the nearest neighbors along the protein chain and had the interbead distance not longer than $d_{\mathrm{cut}}=7.5$ Å. In this case, the number of native contacts is $N_{\mathrm{nat}}=N_{\mathrm{nat}}^{\mathrm{NAT}}=27$ . The simulations were performed with a constant-temperature MD based on the coupled set of Langevin equations biswas1986simulated . The time-step was $\Delta t=0.0125\tau$ , where $\tau$ is the characteristic time. At the length scale $l=7.5\mathrm{\AA}$ Å and the attractive energy $\epsilon=2.2$ kcal/mol miyazawa1996residue , $\tau=(Ml^{2}/\epsilon)^{1/2}\approx 2.6$ ps, where $M=110$ Da is the average mass of the residue. The friction constant $\gamma$ in the Langevin equations was $\gamma=10M/\tau$ . In what follows, the length, in particular, the radius of gyration, is measured in angstroms. Unless otherwise noted, the other quantities (energy, temperature and time) are dimensionless; specifically, the energy and temperature are measured in units $\epsilon$ (in the latter case with the Boltzmann constant set to unity), and the time is measured in units of $\tau$ .

The equilibrium simulations were carried out for a range of temperatures at which both folding and unfolding events occurred quite frequently, specifically, at $T=$ 0.35, 0,375, 0.4, 0.425, and 0.45. The regime of folding varied from downhill (low temperatures) to uphill (larger temperatures) folding through a marginal two-state folding (intermediate temperatures). For each temperature, the MD trajectory was run until $10^{3}$ folding/unfolding events took place. The protein was considered to be unfolded if $N_{\mathrm{nat}}\leq 5$ , and to be folded if the root-mean-square-deviation (RMSD) from the reference native structure was less than 1Å chekmarev2021first . Based on the simulated trajectories, the FESs and FEPs were calculated. The FES was constructed using the number of native contacts ( $N_{\mathrm{nat}}$ ) and the radius of gyration ( $R_{\mathrm{g}}$ ) as collective variables best2013native ; shea2001from

F(N_{\mathrm{nat}},R_{\mathrm{g}})=-T\ln P(N_{\mathrm{nat}},R_{\mathrm{g}})

(1)

where $P(N_{\mathrm{nat}},R_{\mathrm{g}})$ is the probability to find the protein in the point $(N_{\mathrm{nat}},R_{\mathrm{g}})$ . The FEP was calculated as

F(N_{\mathrm{nat}})=-T\ln P(N_{\mathrm{nat}})

(2)

where $P(N_{\mathrm{nat}})$ is the probability that the protein has $N_{\mathrm{nat}}$ native contacts. The resulting FESs and FEPs for the lower, intermediate and upper temperatures are shown in Fig. 1a-c and Fig. 2, respectively. In structure, they are common for folding of $\beta$ -hairpins (experiment munoz1997folding , and simulation studies in explicit zhou2001free ; gao2015correct and implicit dinner1999understanding ; zagrovic2001beta solvents).

Refer to caption — Figure 1: $\beta$ -hairpin, free energy surfaces: (a) $T=0.35$ , (b) $T=0.4$ , and (c) $T=0.45$ . The radius of gyration is given in angstroms.

Both folding and unfolding reactions were considered. For this, the equilibrium MD trajectory was divided into segments of the trajectory corresponding to the transition from the unfolded to the native state (folding) and to the backward transition (unfolding) chekmarev2021first . Based on these segments of the trajectory, the MD rate constants of folding and unfolding, as well as the corresponding FPT distributions, were calculated (see below).

II.2 Monte Carlo Simulations

Having the FESs and FEPs, the MC simulations were performed. According to the Metropolis algorithm metropolis1953equation ; frenkel2004introduction , the acceptance probability of the transition from state A to B is determined as

W(\mathrm{A}\rightarrow\mathrm{B})=\begin{cases}1&{\text{if}}\,\,\,\,P(\mathrm{B}/P(\mathrm{A})\geq 1\\ P(\mathrm{B})/P(\mathrm{A})&{\text{if}}\,\,\,\,P(\mathrm{B})/P(\mathrm{A})<1\end{cases}

(3)

where $P(\mathrm{A})$ and $P(\mathrm{B})$ are the equilibrium probabilities of the old ( $\mathrm{A}$ ) and new ( $\mathrm{B}$ ) states, respectively. A trial move from $\mathrm{A}$ to $\mathrm{B}$ is accepted if $P(\mathrm{B})/P(\mathrm{A})\geq\alpha$ , where $\alpha$ is a random number uniformly distributed between 0 and 1. In application to the free energy landscapes of Figs. 1 and 2, $P(\mathrm{B})/P(\mathrm{A})=\exp\{-[F(\mathrm{B})-F(\mathrm{A})]/T\}$ . For FESs (Fig. 1), $\mathrm{A}$ and $\mathrm{B}$ represent the points in the $(N_{\mathrm{nat}},R_{\mathrm{g}})$ space, and for FEPs (Fig. 2), the points in the $N_{\mathrm{nat}}$ space. The corresponding free energies $F(\mathrm{A})$ and $F(\mathrm{B})$ are determined by Eqs. (1) and (2).

In the case of FESs, the MC simulation space represented a rectangle with $0\leq N_{\mathrm{nat}}\leq 28$ and $0\leq R_{\mathrm{g}}\leq 15$ Å. It was divided into $28\times 30$ segments ( $dN_{\mathrm{nat}}=1$ and $dR_{\mathrm{g}}=0.5$ ), forming a mesh of $29\times 31$ discrete states. The trial moves consisted of random shifts of the point along the $N_{\mathrm{nat}}$ and $R_{\mathrm{g}}$ axes by -1 or +1 mesh step so that the 8 nearest points surrounding the current point were sampled with equal probability. To verify that detailed balance is satisfied, a long ("equilibrium") MC trajectory with many folding/unfolding events was run. It was found that the balance was approximately maintained for each point-to-point transition and was exactly fulfilled for the transitions within the simulation cell, i.e. the number of transitions from the central point of the cell to the eight surrounding points was equal to the number of the backward transitions. Accordingly, the equilibrium MC simulations successfully reproduced the original equilibrium distributions on which they were based (one example is shown in Fig. 2). To simulate the process of folding, the trajectories were initiated at point (5, 11.8), representing the unfolded state, and terminated at point (27, 5.78), representing the native state, (Fig. 1a-c). For the unfolding trajectories, these points served as the terminal and initial states, respectively.

In the case of FEPs, the MC simulations were conducted in the linear $N_{\mathrm{nat}}$ space with the boundaries $N_{\mathrm{nat}}=0$ and $N_{\mathrm{nat}}=27$ . The trial moves represented random changes of the current value of $N_{\mathrm{nat}}$ by -1 or +1, which were taken with equal probability. In equilibrium, detailed balance was strictly maintained. The folding trajectories started at $N_{\mathrm{nat}}=5$ and terminated at $N_{\mathrm{nat}}=27$ , and the unfolding trajectories used these points as terminal and initial state, respectively.

In each case (the FES or FEP based simulations, the protein folding or unfolding, and a given temperature), $10^{3}$ trajectories were run. For the FES based simulations, the acceptance probability increased from $\approx 0.57$ ( $T=0.35$ ) to $\approx 0.67$ ( $T=0.45$ ) for folding trajectories, and from $\approx 0.57$ ( $T=0.35$ ) to $\approx 0.62$ ( $T=0.45$ ) for unfolding trajectories, and for the FEP based simulations, it increased from $\approx 0.87$ ( $T=0.35$ ) to $\approx 0.9$ ( $T=0.45$ ) for both folding and unfolding trajectories.

In order to fit the MC steps to the MD time steps, the average ratios btween the MC and MD rate constants over the whole temperature range were calculated. In particular, for the FES based folding trajectories, the time step was determined as $(\Delta t)^{\mathrm{F}}_{\mathrm{FES}}=(1/N)\sum_{i=1}^{i=N}(k_{\mathrm{F}})_{\mathrm{MC},i}/(k_{\mathrm{F}})_{\mathrm{MD},i}$ , where $N=5$ is the number of temperatures, and $(k_{\mathrm{F}})_{\mathrm{MC}}$ and $(k_{\mathrm{F}})_{\mathrm{MD}}$ are the rate constants for the MC and MD folding trajectories, respectively. The rate constants were calculated as the inverses of the corresponding mean first-passage times (MFPTs) in units of number of steps (the MC simulations) and $\tau$ (the MD simulations). Similarly, the time steps for the FEP based and unfolding MC simulations were determined. It was found that $(\Delta t)^{\mathrm{F}}_{\mathrm{FES}}\approx 3.4\times 10^{-3}\tau$ , $(\Delta t)^{\mathrm{F}}_{\mathrm{FEP}}\approx 6.3\times 10^{-3}\tau$ , and $(\Delta t)^{\mathrm{U}}_{\mathrm{FES}}\approx(\Delta t)^{\mathrm{U}}_{\mathrm{FEP}}\approx 0.18\tau$ . These figures show that the MC time scales generally depend on both the direction of the process (folding or unfolding) and the dimensionality of the state distributions for the MC simulations (FES or FEP).

Figures 3 and 4 compare the results of the MC and MD simulations. Figure 3 presents the rate constants obtained in the MC simulations in comparison to those in the MD simulations. The MC "time” scales were adjusted to the MD time scales using the above values of the MC steps in $\tau$ units, e.g., for the MC folding simulations based on the FES, the adjusted rate constant at $T=0.35$ is obtained to be $(k_{\mathrm{F}})_{\mathrm{MC}}=(\tilde{k}_{\mathrm{F}})_{\mathrm{MC}}/(\Delta t)^{\mathrm{F}}_{\mathrm{FES}}\approx 2.6\times 10^{-3}$ , where $(\tilde{k}_{\mathrm{F}})_{\mathrm{MC}}\approx 8.75\times 10^{-6}$ is the original rate constant in the inverse number of the MC steps, and $(\Delta t)^{\mathrm{F}}_{\mathrm{FES}}\approx 3.4\times 10^{-3}\tau$ . Figure 3 shows that the rate constants from the MC simulations are in good agreement with those from the MD simulations through all temperature range. The same result would evidently be obtained if the MC step in $\tau$ units was determined by equating the MC rate constant at one temperature with the corresponding MD rate constant. The above values of the time steps were also used to adjust the time scales for the MC FPT distributions. Figure 4 presents these distributions in comparison with the distributions obtained in the MD simulations and shows that, similar to the rate constants, they are also in agreement, including a single-exponential decay. Importantly, the good agreement of the MC rate constants with the MD constants through all temperature range evidences that the time step in the MC simulations is not affected by temperature.

III $\alpha_{3}$ D Protein: Experiment and Monte Carlo Simulations

The equilibrium folding/unfolding of this 73-residue three-helix bundle protein (2a3d.pdb) at room temperature and different denaturant (GdmCl) concentrations (1.5 M to 3 M) was studied using single-molecule FRET chung2011extracting . Two models assuming two-state kinetics were employed to analyze the photon trajectories. One is a Maximum Likelihood Estimation (MLE), in which the photon trajectory were fitted to the sequences of transitions between unfolded and folded states gopich2009decoding , and the other is an approximation of the FRET-efficiency histogram (FEH) by a sum of Gaussians whose weights depended on the probabilities of surviving the folded and unfolded states for the bin time (3G) gopich2010fret . Using these approaches, the rate constants of folding and unfolding were determined, which were found in good agreement, except for low values of the denaturant concentration, specifically for 1.5 M and, partly, for 2.0 M (cf. the corresponding lines of Table 1 of Ref. chung2011extracting, ). Figure 5 shows the FEHs for freely

diffusing protein molecules. The time bin to count donor and acceptor photons is 2 ms. Two peaks are significant - at intermediate and higher values of FRET-efficiency ( $E$ ), which represent the unfolded (U) and folded (F) states, as indicated in Fig. 5. The peaks at $E\sim 0$ originate from molecules with inactive acceptors and are not included in consideration. It was found that the FEHs are well approximated by a sum of three Gaussians (3G model). Two Gaussians describe the unfolded and folded states, and the third one accounts for the intermediate values of $E$ that appear due to the transitions between the folded and unfolded states.

It should be noted that the present FEHs may not represent the true PMFs, as it would be if they obeyed the Förster "distance-efficiency" equation $E=1/[1+(r/R_{0})^{6}]$ ( $r$ is the donor-to-acceptor separation distance, and $R_{0}$ is the Förster radius). First, the number of photons in time bin ( $\sim 10^{2}$ ) is not sufficiently large to avoid the effect of shot noise on the distribution of protein states. In this case, the broadening of the folded and unfolded peaks may mainly be a results of short noise rather than of the presence of folded- and unfolded-like conformations of the protein. Second, the binning time (2 ms) is comparable to the relaxation time ( $\sim 1$ ms), so that the signal from different protein conformations, which typically interconvert on a much shorter time scale ( $\sim 1\mu$ s), may be placed in the same bin. In particular, the intermediate values of $E$ between the unfolded and folded peaks may represent a mixture of photons that come from the unfolded and folded states rather than the transitional conformations. However, if detailed balance is fulfilled, it is more important that the intermediate values of $E$ adequately represent the probabilities of the intermediate protein conformations rather than the conformations themselves.

The MC moves were determined by Eq. (3), in which the values of FRET-efficiency represented the protein states, and the numbers of bursts in the FEHs served as (non-normalized) probabilities of these states, $P(E)$ . According to the FEHs of Fig. 5, the whole range of FRET-efficiencies, from 0 to 1, was divided into 31 segments. To take into account that the peak at $E\sim 0$ corresponds to the inactive acceptor, the populations of states at $0<E<E_{\mathrm{thr}}$ , where $E_{\mathrm{thr}}$ is the value of $E$ at which $P(E)$ has a minimum between the inactive and unfolded peaks, were set to zero, $P(E)=0$ . Each MC step was taken to be an equally probable random move to a neighboring segment. Detailed balance was strictly satisfied similar to the MC simulations for $\beta$ -hairpin (Sect. II.2). The unfolded and folded states were chosen to correspond to the maxima of the unfolded and folded peaks (Fig. 5), with the values of FRET-efficiency $E_{\mathrm{U}}=0.61$ and $E_{\mathrm{F}}=0.93$ , respectively. To simulate folding, the MC trajectories were initiated at $E=E_{\mathrm{U}}$ and terminated at $E=E_{\mathrm{F}}$ , and to simulate unfolding, the states $E_{\mathrm{F}}$ and $E_{\mathrm{U}}$ were used as the initial and terminal states, respectively. In all cases, i.e., at each denaturant concentration for folding and unfolding, $10^{3}$ MC trajectories were run. For both folding and unfolding reactions, the acceptance probability varied insignificantly, from $\approx 0.87$ to $\approx 0.92$ . Based on the simulated folding and unfolding trajectories, the rate constants for folding ( $k_{\mathrm{F}}$ ) and unfolding ( $k_{\mathrm{U}}$ ) were determined, and the corresponding FPT distributions were constructed. Similar to $\beta$ -hairpin, the MC “time” scale was adjusted to the experimental time scale by calculating the average ratios between simulated and experimental rate constants. Specifically, the MC time steps for folding and unfolding were determined, respectively, as $(\Delta t)^{\mathrm{F}}=(1/N)\sum_{i=1}^{i=N}(k_{\mathrm{F}})_{\mathrm{MC},i}/(k_{\mathrm{F}})_{\mathrm{exp},i}$ and $(\Delta t)^{\mathrm{U}}=(1/N)\sum_{i=1}^{i=N}(k_{\mathrm{U}})_{\mathrm{MC},i}/(k_{\mathrm{U}})_{\mathrm{exp},i}$ , where $N=3$ is the number of denaturant concentrations ranging from 2.0 M to 2.5 M (i.e., 2.0 M, 2.25 M, and 2.5 M), at which the relative populations of both folded and unfolded states were larger than $1\%$ . It was obtained $(\Delta t)^{\mathrm{F}}\approx 7.5\times 10^{-3}$ ms and $(\Delta t)^{\mathrm{U}}\approx 8.7\times 10^{-3}$ ms. Given the MC time steps, the rate constants for folding and unfolding were recalculated as $(k_{\mathrm{F}})_{\mathrm{MC}}=(\tilde{k}_{\mathrm{F}})_{\mathrm{MC}}/(\Delta t)^{\mathrm{F}}$ and $(k_{\mathrm{U}})_{\mathrm{MC}}=(\tilde{k}_{\mathrm{U}})_{\mathrm{MC}}/(\Delta t)^{\mathrm{U}}$ , where the rate constants denoted by tilde are the original MC constants in the inverse number of the MC steps.

Figure 6a-b compares the adjusted MC rate constants with those extracted from the FRET-trajectories with the 3G model (Table 1 of Ref. chung2011extracting, ). It shows that the MC simulations based on the FRET-histograms reasonably predict the change of the rate constants with denaturant concentration, for both folding and unfolding reactions. The only exception are the concentrations at which the normalized population of the target state is as small as $\approx 1\%$ , as at 3M GdmCl in the case of folding, where $P(0.93)\approx 1.2\times 10^{-2}$ , and at 1.5M GdmCl in the case of unfolding, where $P(0.61)\approx 9\times 10^{-3}$ . More generally, the MC simulations are expected to give satisfactory results when the peaks for the target states are resolved. Figure 7a-b also presents the simulated FPT distributions for folding and unfolding. The distributions are single-exponential, which is consistent with the assumption of two-state kinetics that was made in the theoretical models for extracting the rate constants from the FRET data chung2011extracting .

IV Conclusions

It has been shown that the Metropolis MC simulations based on the equilibrium distributions of states can be successfully used to predict how characteristic reaction times in a complex system change with the determining conditions. As an example, the process of folding and unfolding for two proteins was considered. The first protein, a $\beta$ -hairpin protein, whose coarse-grained (C_α-bead) model was constructed on the basis of the solution NMR data for 12-residue HP7 protein andersen2006minimization , was employed to develop and test the approach. Performing MD simulations, the baseline characteristics of the folding/unfolding process were determined for a set of temperatures - the equilibrium free energy surfaces (FESs) and free energy profiles (FEPs), the rate constants of folding and unfolding, and the first-passage time (FPT) distributions for these reactions. Using the FESs and FEPs as the potentials of mean-force (PMFs), the MC simulations were carried out to obtain relative rate constants of folding and unfolding and the FPT distributions. It has been found that they change with temperature in good agreement with those from the MD simulations and, after fitting the MC steps to the MD time steps, they agree well in absolute value. The other protein is a three-helical $\alpha_{3}$ D protein, whose folding/unfolding reaction was studied at different different denaturant concentrations using single-molecule FRET method chung2011extracting . In this case, the measured FRET-efficiency histograms (FEHs) were employed as the equilibrium distributions of states. Because of the uncertainty posed by experimental limitations (a relatively small number of photons in time bins, and a long binning time in comparison to the characteristic times of interconversion of protein molecules), it was not quite clear whether these FEHs could be used as the PMFs. However, it has been found that for the denaturant concentrations at which the FEH peak for the target folded (unfolded) state is resolved, the rates constants of folding (unfolding) change with the denaturant concentration in satisfactory agreement with those extracted from photon trajectories with the use of theoretic models gopich2009decoding ; gopich2010fret . Also, it has been found the FPT distributions are single-exponential, which is consistent with the assumption of two-state kinetics in the theoretical models gopich2009decoding ; gopich2010fret for extracting the rate constants from the FRET data. To obtain absolute values of the rate constants, a time-dependent quantity is required, which, on the one hand, would be known in the experiment (physical or computational) and, on the other, could be obtained in the MC simulations based on the given equilibrium distribution of states, e.g., the self-diffusion coefficient, if this distribution were in physical space huitema1999can . The promising feature of the present approach is that it is based solely on the equilibrium distributions of states, without introducing any additional parameters to perform simulations, which suggests its applicability to other complex systems.

Acknowledgements.

I thank Irina Gopich for valuable discussions of the results of the FRET-experiments for folding of

\alpha_{3}

D protein, and Hoi Sung Chung for sharing the corresponding FRET-efficiency histograms. The work was supported by a contract with the IT SB RAS.

Author Declarations

Conflict of Interest
The authors have no conflicts to disclose.

Data Availability Statement

AVAILABILITY OF DATA	STATEMENT OF DATA AVAILABILITY
Data available in article or on request from the author	The data that support the findings of this study are available within the article or from the author upon reasonable request.

References

(1) B. Schuler, E. A. Lipman, and W. A. Eaton, Nature 419, 743 (2002).
(2) F. Huang, L. Ying, and A. R. Fersht, Proc. Natl. Acad. Sci. USA 106, 16239 (2009).
(3) H. S. Chung, I. V. Gopich, K. McHale, T. Cellmer, J. M. Louis, and W. A. Eaton, J. Phys. Chem. A 115, 3642 (2011).
(4) J. Liu, L. A. Camposa, M. Cerminara, X. Wang, R. Ramanathan, D. S. Englisha, and V. Muñoz, Proc. Natl. Acad. Sci. USA 109, 179 (2012).
(5) B. Schuler and H. Hofmann, Curr. Opin. Struct. Biol. 23, 36 (2013).
(6) W. A. Eaton, J. Phys. Chem. B 125, 3452 (2021).
(7) H. Yu, A. N. Gupta, X. Liu, K. Neupane, A. M. Brigley, I. Sosova, and M. T. Woodside, Proc. Natl. Acad. Sci. USA 109, 14452 (2012).
(8) M. Batchelor, K. Papachristos, M. Stofella, Z. T. Yew, and E. Paci, Biochim. Biophys. Acta 1864, 129613 (2020).
(9) C. Bustamante, L. Alexander, K. Maciuba, and C. M. Kaiser, Annu. Rev. Biochem. 89, 443 (2020).
(10) S. A. McKinney, C. Joo, and T. Ha, Biophys. J. 91, 1941 (2006).
(11) I. V. Gopich and A. Szabo, J. Phys. Chem. B 113, 10965 (2009).
(12) I. V. Gopich and A. Szabo, J. Phys. Chem. B 114, 15221 (2010).
(13) R. Ramanathan and V. Muñoz, J. Phys. Chem. B 119, 7944 (2015).
(14) I. S. Osad’ko, J. Phys. Chem. C 120, 697 (2016).
(15) G. Hummer and A. Szabo, Proc. Natl. Acad. Sci. USA 98, 3658 (2001).
(16) G. Hummer and A. Szabo, Acc. Chem. Res. 38, 504 (2005).
(17) C. Hyeon and D. Thirumalai, J. Phys.: Condens. Matter 19, 113101 (2007).
(18) O. K. Dudko, G. Hummer, and A. Szabo, Proc. Natl. Acad. Sci. USA 105, 15755 (2008).
(19) D. E. Shaw, P. Maragakis, K. Lindorff-Larsen, S. Piana, R. O. Dror, M. P. Eastwood, J. A. Bank, J. M. Jumper, J. K. Salmon, Y. Shan, and W. Wriggers, Science 330, 341 (2010).
(20) K. Lindorff-Larsen, S. Piana, R. O. Dror, and D. E. Shaw, Science 334, 517 (2011).
(21) C. Camilloni and F. Pietrucci, Advances in Physics: X 3, 1477531 (2018).
(22) Y. Sugita and Y. Okamoto, Chem. Phys. Lett. 314, 141 (1999).
(23) Y. M. Rhee and V. S. Pande, Biophys. J. 84, 775 (2003).
(24) F. Rao and A. Caflisch, J. Chem. Phys. 119, 4035 (2003).
(25) A. M. Ferrenberg and R. H. Swendsen , Phys. Rev. Lett. 63, 1195 (1989).
(26) J.-E. Shea and C. L. Brooks III, Annu. Rev. Phys. Chem. 52, 499 (2001).
(27) N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. N. Teller, and E. Teller, J. Chem. Phys 21, 1087 (1953).
(28) D. Frenkel, Introduction to Monte Carlo Methods, Computational Soft Matter: From Synthetic Polymers to Proteins, Lecture Notes, John von Neumann Institute for Computing, Jülich, NIC Series, N. Attig, K. Binder, H. Grubmüller, K. Kremer (Eds.), 2004, 23, 29-60.
(29) S. F. Chekmarev, A. Yu. Palyanov, and M. Karplus, Phys. Rev. Lett. 100, 018107 (2008).
(30) I. V. Kalgin, S. F. Chekmarev, and M. Karplus, J. Phys. Chem. B 118, 4287 (2014).
(31) B. Roux, Comput. Phys. Comm. 91, 275 (1995).
(32) D. Trzesniak, A.-P. E. Kunz, and W. F. van Gunsteren, ChemPhysChem 8, 162 (2007).
(33) H. E. A. Huitema and J. P. van der Eerden, J. Chem. Phys. 110, 3267 (1999).
(34) G. Rutkai and T. Kristóf, J. Chem. Phys. 132, 104107 (2010).
(35) E. Sanz and D. Marenduzzo, J. Chem. Phys. 132, 194102 (2010).
(36) K. M. Bal and E. C. Neyts, J. Chem. Phys. 141, 204104 (2014).
(37) K. Kikuchi, M. Yoshida, T. Maekawa, and H. Watanabe, Chem. Phys. Lett. 185, 335 (1991).
(38) G. Tiana, L. Sutto, and R. A. Broglia, Physica A: Stat. Mech. Appl. 380, 241 (2007).
(39) N. Gō, Annu. Rev. Biophys. Bioeng. 12, 183 (1983).
(40) S. F. Chekmarev, J. Chem. Phys. 139, 145103 (2013).
(41) N. H. Andersen, K. A. Olsen, R. M. Fesinmeyer, X. Tan, F. M. Hudson, L. A. Eidenschink, and S. R. Farazi, J. Am. Chem. Soc. 128, 6101 (2006).
(42) T. X. Hoang and M. Cieplak, J. Chem. Phys. 112, 6851 (2000).
(43) R. Biswas and D. R. Hamann, Phys. Rev. B 34, 895 (1986).
(44) S. Miyazawa and R. L. Jernigan, J. Mol. Biol. 256, 623 (1996).
(45) S. F. Chekmarev, Phys. Chem. Chem. Phys. 23, 17856 (2021).
(46) R. B. Best, G. Hummer, and W. A. Eaton, Natl. Acad. Sci. USA 110, 17874 (2013).
(47) V, Muñoz, P. A. Thompson, J. Hofrichter, and W. A. Eaton, Nature 390, 196 (1997).
(48) R. Zhou, B. J. Berne, and R. Germain, Proc. Natl. Acad. Sci. USA 98, 14931 (2001).
(49) Y. Gao, Y. Li, L. Mou, B. Lin, J. Z. H. Zhang, and Y Mei, Sci. Rep. 5, 10359 (2015).
(50) A. R. Dinner, T. Lazaridis, and M. Karplus, Natl. Acad. Sci. USA 96, 9068 (1999).
(51) B. Zagrovic, E. J. Sorin, and V. Pande, J. Mol. Biol. 313, 151 (2001).