Details of simulations and comparison with BW analysis

To study the dynamics of this model we used the Monte Carlo algorithm with local moves, which include end moves, corner flips and crankshaft moves. The excluded volume constraint was enforced by not allowing any multiple occupancy. (Further details can be found in a previous work.[10]) To calculate the mean first passage time ( ) as a function of temperature the following procedure was used. For each temperature many simulations were performed, each starting from a different random unfolded initial condition. The simulations were run until the chain found the folded state or until a maximum number of time steps had elapsed. The mean first passage time is simply the average of the folding times for the different runs. For an intermediate range of temperatures all of the runs would find the folded state within steps. However for very low and high temperatures only some fraction would find the folded state in the allowed time. For these temperatures the mean calculated is a lower bound to the actual mean first passage time. This is because for trials that do not fold within steps, we average in which is less than the actually folding time for that run. This is done for practical reasons since we have a finite amount of computer time for the simulations. The key point is that is much greater than for a broad range of temperatures; which in fact, it is. Figure 2 shows the mean first passage time simulated for a range of temperatures.

At low temperatures, grows rapidly. This slow down is caused by trapping in local minima. It is well known that heteropolymers undergo an ideal glass transition at low enough temperature. The random energy model estimate of this temperature for the present system is . The measured energetic fluctuations 's are 51 at and 22 at T=1.1 (slightly above the glass transition). The respective configurational entropies for the collapsed states, , are around 20 and 12. Both cases provide a . This agrees with the critical temperature at which replica symmetry breaking sets in. This is determined by finding the temperature at which becomes of order 1 (see figure 8).

  
Figure 8: Qualitative replica symmetry breaking value plot at three temperatures ( , and ). Note that although at the global glass temperature ( ) discrete states are apparent even for small degrees of nativeness, at the folding temperature the discrete intermediates are largely native-like.

is the Boltzmann occupation of a microstate. Y is not the perfect measure of replica symmetry breaking since it assumes the basin of attraction of a microstate can be identified with the state itself. The corresponding neglect of correlation between these high Q states, however, is able to provide a reasonable estimate of for a finite mesoscopic system.

These estimates of are very similar to those of the kinetic glass transition which is the temperature at which the folding time slows down substantially. We define the kinetic glass transition temperature ( ) to be the temperature at which the folding time is sufficiently slower than the fastest folding time observed. There are several ways of choosing this point which will give roughly the same answer. We choose , to be the temperature at which . For this sequence, .

In addition to the kinetic runs to calculate , we performed a number of thermodynamic runs to determine the free energy profiles. Again most of the details can be found in a previous work.[11] We used the Monte Carlo histogram technique which allows us to calculate extensive quantities like the free energy. The basic idea is to store a histogram as a function of Q and E. These histograms measure the probability that a given (Q,E) pair will occur. From these we can calculate the density of states n(Q,E) up to a multiplicative constant (which can be determined since we know the ). We can then calculate any thermodynamic quantity; in particular, we determined the free energy as a function of Q (F[Q]) for several temperatures. The free energy is shown in figure 5.

  
Figure: Autocorrelation function for Q (see eq. 9) plotted on a semi-log scale. . The correlations have an initial rapid decay followed by a slower decay. By fitting a sum of exponentials (usually 3 or 4) the autocorrelation time can be determined.

For a broad range of temperatures the free energy has a rough double well profile, with one minimum at the folded state (Q=28) and another for the molten globule ( -7). The transition state varies from at T=2 to at the lower temperatures. At first it might appear the transition state is at since this is the highest local maximum in the plot at most temperatures.[20] However, this turns out not to be the correct transition point. One must be careful in interpreting free energy plots of this type. There are Monte Carlo moves that connect states with Q=23 directly to the ground state short circuiting this barrier. In fact any barrier that is smaller in width than 5 can also be by passed by certain Monte Carlo moves. Note, the correct barriers from 22 to 16 are much broader. In fact it is better to think in terms of a transition ``region'' rather than a unique transition point. The way in which we identified this region was to look at the time trajectories (figure 4). Since we are in a diffusive regime, a trajectory that reaches the transition region should have some probability of diffusing back down the barrier rather than crossing over. For values between 16-22 we see such an effect. However, by the time is chosen there would be virtually no trajectories that cross back without first reaching the folded state. We have also previously used a 2 dimensional analysis of the trajectories and free energy surfaces and again we find the transition region to be between 16-22 and not at 25.[3] Since we know that the transition state is somewhere between 16 and 22 we define a specific transition point, for use in calculating the barrier for the Kramer's rate equation, as the value of Q that is a maximum in this region.

In addition to the free energy barrier we need to measure the diffusion coefficient in this system in order to calculate the folding time. We compute the diffusion coefficient by measuring the autocorrelation time of Q. Specifically we calculate:

 

Figure 9 shows a semi-log plot of the auto correlation function at T=1.509, the folding temperature. The figure is multi-exponential with a short fast decay followed by a much longer, slower decay. The final drop-off is mostly due to errors in sampling at very long times. We are interested in the long time decay and calculated the exponent associated with this by fitting a function of several exponentials (usually 3 or 4) to this plot. We repeated this calculation at several temperatures (note that we can not use the histogram method here do the large amount of histogram information that would have to be stored). Figure 10 is a log-log plot of the autocorrelation function for several temperatures. As the temperature is decreased the autocorrelation time increases.

  
Figure 10: Q Autocorrelation function plotted on a log-log scale for several temperatures. As the temperature decreases the correlation time increases.

Finally the diffusion coefficient was estimated as .

With the diffusion coefficients and free energy surfaces in hand, it is now possible to test the analytical prediction given by eq. 2. Using the discrete form of the full double integral (with a constant diffusion coefficient):

 

we obtain the results presented in figure 11.

  
Figure: Comparison of mean first passage times from Monte Carlo simulations with folding times calculated both the discrete from of the double integral (eqs. 2 and 10) and the simple rate formula which assumes a constant diffusion coefficient and that the curvature at the top and bottom of the free energy are the same.

Note the sum is not taken to Q=28 since we do not expect the diffusion coefficient to remain constant in that region (between Q=23-28). However, the time it takes to go from Q=23 to the folded state (Q=28) is small compared to (as seen in figure 6) so the error from truncating the sum will also be small. Overcoming the barrier ( ) is the limiting step. Additionally, we can also use a simplified form of eq. 2 which assumes that the free energy surface is well approximated by a parabolic potential around the molten globule states that extends all the way to the barrier (eq. 3). If we further assume that and that the two curvatures are the same then the folding time is given by (which is also shown in figure 11). This simple result is useful because it depends only on the correlation time in the molten globule and the free energy barrier. The agreement with the full formula is remarkable. For both formulas, the fastest time at the optimal kinetic folding temperature is obtained to well within a factor of two. Also, an extremely good description of the behavior for the slow down anticipated for low and high temperatures is obtained. Thus we see the analytical theory based on the actual molten globule dynamics and funnel free energy profile is not simply qualitatively correct but can be used for quantitative predictions of the folding time over a wide thermodynamic range.

Up to this point, the values for the diffusion coefficients have been obtained directly from the simulation data for the time correlation functions around the molten globule. It is interesting now to compare the obtained results with the existing models of configurational diffusion, particularly the BW theory. Since our results fall in the range between and , the intermediate temperature formula given by eq. 7 would be seem to be the appropriate comparison. Figure 12 shows a plot of the diffusion coefficients as a function of .

  
Figure 12: Diffusion coefficient plotted versus

The slope obtained from the semi-log plot provides a fitted value for of . This is a little bit smaller than the 's obtained directly from simulations at this temperature range. This is expected because correlations in the landscape cause only part of the total energetic ruggedness to come into the activation barriers for escaping from traps. Qualitatively agreement between the model and simulations suggests that while the REM results are a first approximation, the correlations in the energy landscape are significant and the effective fluctuations for escape from a trap involves only a fraction of the chain.[18, 21] We would expect these correlation effects to become more important as the effective chain length of the protein grows, however we should bear in mind that the same qualitative behavior still applies. In the temperature range studied, we note the data could be equally well fitted by an Arrhenius temperature plot (see figure 13).

  
Figure 13: Diffusion coefficient plotted versus and

The effective activation energy is ( ). This very large number implies that the effective moves of Q require substantial organization of the protein since the stability gap is only . The simple Ferry law fit is also adequate providing a (see eq. 6 and figure 13).