Thermodynamics and Folding Pathways of Trpzip2: An Accelerated

Department of Chemistry, Texas A&M University, College Station, Texas 77843 ..... The physical meaning of short trajectories shot from the barrier for...
0 downloads 0 Views 964KB Size
Download PDF
J. Phys. Chem. B 2009, 113, 803–808

803

Thermodynamics and Folding Pathways of Trpzip2: An Accelerated Molecular Dynamics Simulation Study Lijiang Yang, Qiang Shao, and Yi Qin Gao† Department of Chemistry, Texas A&M UniVersity, College Station, Texas 77843 ReceiVed: April 11, 2008; ReVised Manuscript ReceiVed: NoVember 21, 2008

In this paper, we apply an enhanced sampling method introduced earlier to study the folding mechanism of a β-hairpin, trpzip2, using an all-atom potential for the protein and an implicit model for the solvent. The enhanced sampling method allows us to obtain multiple protein folding and unfolding trajectories in relatively short simulations. The sufficient sampling of folding and unfolding events of trpzip2 makes possible a more detailed investigation of its folding landscape and thermodynamics, leading to the identification of folding pathways. The analysis of the thermodynamics involved in the folding of trpzip2 showed that this polypeptide folds by two stages: a downhill hydrophobic collapse followed by formation of native hydrogen bonds. During the hydrogen bond formation, a transition state can be identified with only one native hydrogen bond being formed, which is more consistent with a ‘zip-out’ mechanism. To address the dynamics in a more reliable way, explicit solvent was used for the estimation of the diffusion constant, which was used in the Kramer’s theory, together with the free energy profile calculated using the enhanced sampling and the implicit solvent, to calculate the folding rate. The calculated folding time agrees well with the experimental value of 2.5 µs. Introduction Understanding of the protein folding mechanism is one of the most challenging problems in biological science. Due to the complexity of the protein folding process, great efforts have been paid in both experimental and theoretical studies during the past several decades.1-9 Over the recent 10 years, numerous computer simulations at the all-atom level have been performed to investigate protein folding processes from the statistical point of view as a result of the great improvement of the computational power. These simulations include, for example, the microsecond time-scale molecular dynamics folding simulation on the villin headpiece10 and the submillisecond atomistic protein folding simulation using a network of tens of thousands of PCs around the world.11,12 At the same time, numerous unfolding simulations have been performed13-16 to understand the folding mechanism since (1) such simulations, starting from a crystal or NMR structure, improve the sampling of the experimentally relevant regions of conformational space, and (2) by introducing high temperatures, energetic barriers are easily overcome in unfolding simulations. However, it was argued that the basis of unfolding simulations, “the principle of microscopic reversibility”, is a statement made for equilibrium systems and might not hold under the strongly nonequilibrium conditions of unfolding simulations.17 Recently, we folded the mini protein Trp-cage successfully using the accelerated molecular dynamics (MD) simulation method.18 Among 24 independent 10 ns accelerated MD simulations, 22 reached the folded states, and the best structure obtained in the simulations had a very low backbone/heavyatom rmsd (0.62/1.39 Å) compared with the NMR structure. In the accelerated MD simulation, a bias potential is added to control the enhanced sampling in a desired energy range without under-sampling the low energy or oversampling the high energy configurations, which therefore leads to a more uniform † To whom correspondence [email protected].

should

be

addressed.

E-mail:

distribution as a function of energy.19 As a result, the sampling of the conformational space can be very efficient in accelerated MD simulations. When a suitable bias potential is used, the concomitance of folding and unfolding events can be observed in the same trajectory since barriers between folded, intermediate, and unfolded states are decreased so that transitions in both folding and unfolding directions happen more easily. This feature makes the accelerated MD method a useful tool for thermodynamics calculations, since important regions of conformational space near both folding and unfolding basins can be sampled. In this study, we apply the accelerated MD simulation method to the folding of a stable β-hairpin, tryptophan zipper 2 (trpzip2, SWTWENGKWTWK). Trpzip2 is designed by Cochran et al.20 and has a strong energy bias toward its native structure. Trpzip2 has been proven to be a good, simple system to study the folding mechanism of β-hairpins and has been investigated in several computational studies.21-24 In the current study, trpzip2 was folded and unfolded iteratively through accelerated MD simulations, which allowed more detailed and quantitative investigation of its folding landscape and thermodynamics, leading to the identification of folding pathways to the native state, the estimation of the folding rate, and the understanding of the general principles in the folding of β-hairpins. Methods Accelerated MD Simulation. In the accelerated MD simulation (umbrella sampling), a bias potential is added to the original potential to enhance the sampling of high energy ranges and thus allows more frequent transitions. A straightforward way of adding a bias potential to the potential energy function of a complex system to enhance sampling of the important but rarely sampled high energy states has already been suggested in the original paper by Torrie and Valleau.25 Recently, in a series of publications, Hamelberg et al.26-29 proposed and applied a method to alter the potential energy landscape, leading to an

10.1021/jp803160f CCC: $40.75  2009 American Chemical Society Published on Web 12/29/2008

804 J. Phys. Chem. B, Vol. 113, No. 3, 2009

Yang et al.

efficient sampling of the conformational space. In this method, the sampling over the high energy range was increased at the cost of a largely reduced sampling over the states with lower energies. More recently, we proposed an approach19 to accelerate MD simulations in which a bias potential was added to control the enhanced sampling in a desired range of energy (e.g., near the energy barrier and in the entire important energy range, which includes the energies of states that are highly populated at the desired temperature as well as energies that are high enough for the reaction to occur) without under-sampling the low energy or over-sampling the high energy range, leading to a more uniform distribution as a function of energy. The bias potential can be written as a sum of Gaussian terms:18 n

f(V) )

∑ aie-[(V - V )/σ ] i

i

2

(1)

i

where Vi is the energy of a state that is to be sampled with enhancement, ai is a negative constant that defines how much the enhancement will be. The parameters in eq 1 are chosen in a way that the magnitude of f(V) is small for energies that are heavily sampled in the low temperature normal MD simulation and for energies that are unimportant (e.g., very high or low energies) and is large otherwise. First, several short normal MD simulations at different temperatures (e.g., one at room temperature, one at a relatively lower temperature, and one at a relatively higher temperature 400 K) would be run to find out a reasonable energy range to be explored. Then, an iterative procedure is used to determine a suitable bias potential18 that would enable good samplings of the whole predetermined energy range. (The enhanced sampling in the energy space was equivalently achieved using a temperature-based method.30,31) An accelerated MD simulation with the biased potential function V˜(r) ) V + f(V) at equilibrium yields a distribution function, over the configuration space, ˜ ˜ F˜(r) ) e-βV(r) /Q

(2)

where β ) 1/kBT, kB being the Boltzmann constant, T being ˜ ) ∫e- β˜ V(r) dr. Becuase V˜(r) is a known the temperature, and Q function of V(r), the distribution function with original (unbiased) potential F(r) can be recovered from F˜ (r) as: ˜

F(r) ) e-βV(r) /Q ) e-βV(r) × ˜

˜

˜ /Q e-β[V(r)-V(r)] /Q ) F˜ (r)e-β[V(r)-V(r)]Q

(3)

where Q ) ∫e- βV(r) dr. Other thermodynamic properties can also be obtained in a similar fashion. Simulation Setup and Analysis. All simulations were carried out by using the AMBER 9.0 simulation package. The protein was modeled with the AMBER FF9632 all-atom force field and the generalized Born/solvent accessible surface area (GB/SA) implicit solvent model.33,34 The salt concentration was set to 0.2 M, and the default surface tension was 0.005 kcal/mol/Å2. The SHAKE algorithm35 with a relative geometric tolerance of 10-5 was used to constrain all bonds including hydrogen. The system was maintained at the objective temperature 300 K using the weak-coupling algorithm36 with a coupling constant of 5.0 ps-1. No nonbonded cutoff was used in simulations. First, short simulations were run at the desired temperature 300 K to yield information on energy distributions. On the basis of the

Figure 1. The sampled potential energy in the accelerated and normal MD simulations. (a) Potential energy visited in the 10 ns normal (black) and accelerated (red) MD simulations; (b) potential energy probability distribution function calculated from all four accelerated MD simulations.

knowledge of energy distributions, a bias potential in the form of eq 2 (n ) 5, a1 ) -1.5, V1 ) -400, σ12 ) 1000; a2 ) -0.8, V2 ) -380, σ22 ) 300; a3 ) -0.032, V3 ) -345, σ32 ) 200; a4 ) -1.6, V4 ) -300, σ42 ) 300, a5 ) -3.0, V5 ) -290, σ52 ) 600; a, V, σ are all in units of kcal/mol) was added to broaden the energy distribution so that the high energies that are only significantly populated at high temperatures will also be efficiently sampled in a simulation at the lower temperature. Then four independent accelerated MD simulations using the bias potential were performed and each one of them was extended to 500 ns. The root-mean-square deviation from experiment structures was calculated by using NMR structure 1LE1 as the reference. Besides rmsd, radius of gyration, Rg, and number of native backbone hydrogen bond were also analyzed to calculate the potential of mean force of trpzip2. Only CR atoms were counted when calculating radius of gyration. To characterize the folding mechanism, hierarchical clustering analysis was applied with pairwise rmsd of 3 Å. Results Folding Landscape and Pathway. Starting from the fully extended structure of the peptide chain trpzip2, four independent sets of 500 ns accelerated MD simulations (see Methods section) were conducted at 300 K. With the introduction of the bias potential, the energy range explored in accelerated MD simulations became considerably broader (between -390 and -280 kcal/mol) than that explored in normal MD simulations (between -380 and -330 kcal/mol). Figure 1a shows the potential energy distributions in the 10 ns accelerated and normal MD simulations. The potential energy probability distribution function calculated from all four accelerated MD simulations is shown in Figure 2b, which indicates that all four simulations have good sampling of both high and low energy regions. The sampling of high energy region enables the escaping from local energy minima, and the sampling of low energy regions helps to identify folded states whenever they are approached, so that the searching of the global minimum is accelerated.18 Another effect of the sampling of high energy regions is that both folding and unfolding events could be easily observed in accelerated MD simulations. The enhanced sampling of high energy regions easily initiates the unfolding from the native basin, as can be seen from Figure 2 (the evolutions of CR-rmsd in the four accelerated MD

Thermodynamics and Folding Pathways of Trpzip2

J. Phys. Chem. B, Vol. 113, No. 3, 2009 805

Figure 2. The CR-rmsd distributions in the four 500 ns accelerated MD simulations. Multiple folding/unfolding events are observed. Figure 4. The natural logarithm of the visiting probability as a function of CR-rmsd and potential energy. (The interval between contour lines is 1.)

Figure 3. The best folded structure (blue) from accelerated MD simulations with 0.53 Å CR-rmsd and 1.05 Å heavy-atom rmsd compared with the NMR structure 1LE1 (red). The four tryptophan residues which form the hydrophobic core of trpzip2 are shown with the stick representation.

simulations). If we define the state with CR-rmsd less than 2.0 Å as the folded state, then a number of folding and unfolding events were observed in each 500 ns accelerated MD simulation. The best folded structure is shown in Figure 3. It has a 0.53 Å CR-rmsd and a 1.05 Å heavy-atom rmsd. In addition to a high degree of resemblance between the backbone structures, the packing patterns of the hydrophobic core are almost identical: the four tryptophan residues Trp2, Trp4, Trp9, and Trp11 are tightly packed against each other and in nearly the exact same patterns as those in the NMR structure. Figure 4 shows the logarithm of the visiting probability as a function of the CRrmsd and the potential energy. Apparently, the folded states locate mainly in the low energy region (-360 to -350 kcal/ mol). This observation agrees with the experimental finding: trpzip2 is strongly biased toward its native structure.20 From another facet, it indicates that the force field we used, AMBER FF96,32 is suitable for the folding simulation of β-hairpin trpzip2. To detect the populated conformations from the four 500 ns trajectories, hierarchical clustering analyses37,38 were conducted on the ensemble of 20 000 conformations (1 snapshot/100 ps). Clustering was based on the CR-rmsd between structures, and the cluster radius was selected as 3 Å. Four most-populated clusters were identified from the clustering analysis. They represent four different basins of attractions: the native state N, a partially folded state P, an intermediate compact state C, and a misfolded state M. Most of the other clusters can be combined to form an unfolded state U. To map all states onto the folding landscape, the radius of gyration of the hydrophobic core

Figure 5. The free energy landscape as a function of the number of native backbone hydrogen bonds nHB and the radius of gyration of the hydrophobic core Rcore g . The compact state, the misfolded state, the partially folded state, and native state are labeled on the surface as C, M, P, and N, respectively. (The interval between contour lines is 2kBT.)

Rgcore(i.e., the four Trp residues) and the number of native backbone hydrogen bonds nHB were calculated for the four states to classify the different states quantitatively. The conformations belonging to state N have a number of native backbone hydrogen bonds nHB ) 5 and an average of radius j gcore ) 4.1 Å; the state P, of gyration of the hydrophobic core R j gcore ) 4.1 Å; is defined within the range of 2 e nHB e 4 and R core j the state C, nHB e 1 and Rg ) 5.7 Å; and the state M, nHB j gcore ) 4.6 Å. e 1 and R The folding landscape constructed as a function of nHB and is shown in Figure 5, and the four major clusters are labeled Rcore g on the free energy surface. According to Figure 5, a possible folding pathway along reaction coordinates nHB and Rgcore can be deduced. Starting from the extended structure, the peptide first collapses to the state C with Rcore g decreasing from more than 12 Å to around 6 Å. (The ordering of events is based on the progression along reaction coordinates, not a progression in time.) All hydrophobic residues prefer a more compact structure for stronger interactions. This observation indicates that the essential driving force of the early folding process is the hydrophobic interactions. It was also observed that the hydrogen bond is not the main driving force at this stage since almost none of the native backbone hydrogen bonds are formed. The peptide remains to be at a compact structure during subsequent folding processes, which is evidently shown by the L-shape

806 J. Phys. Chem. B, Vol. 113, No. 3, 2009

Yang et al.

Figure 6. The free energy landscape as a function of CR-rmsd and (a) the number of native backbone hydrogen bonds nHB and (b) the radius of gyration of the hydrophobic core Rcore g . (The interval between contour lines is 2kBT.)

feature of the folding landscape (Rgcore is centered on 4 Å after trpzip2 folded to the state P). In the latter stage, the hydrophobic interactions no longer drive the folding and hydrogen bond interactions become the major driving force: the fluctuation of Rgcore remains rather small while the number of native backbone hydrogen bonds increases gradually. To show the folding mechanism more clearly, the free energy surfaces as a function of CR-rmsd and Rgcore/nHB are drawn in Figure 6. The similar folding pattern is observed. As shown in Figure 6a, at the early stage, the CR-rmsd decreases quickly along with the decrease of Rcore g , indicating that the hydrophobic collapse is the dominant event in the folding. However, after Rgcore reaches 6 Å, the decrease of CR-rmsd slows down, which indicates that the protein has already reached a rather compact structure and that the hydrophobic interactions mainly play a role of stabilizing the compact structure rather than driving the peptide to fold further. On the contrary, Figure 6b shows that the hydrogen bond interaction has little contribution in the early folding stage but it becomes the driving force in the later stage. Apparently, the folding of trpzip2 is not a downhill process. There is a free energy barrier of ∼2-3 kcal/mol on the folding free energy landscape, which corresponds to the formation of the first backbone hydrogen bond. The unfolding of trpzip2 is simply the reverse process of the folding according to Figures 5 and 6. First, the breaking of native hydrogen bonds initiates the unfolding process. The compact structures are then disassembled because of the loss of essential hydrophobic interactions. The unfolding process proposed here agrees with the finding in the previous replica exchange simulations and experiments by Gruebele et al.21 As mentioned before, a misfolded state was also observed in the simulations, although it consists only a very small portion of the ensemble (4.2%, reweighted from the direct population 3.7%. The reweighting procedure is described in an earlier paper19). The misfolded state has a Rgcore value close to that of the state N, but there is no native backbone hydrogen bonds formed. In Figure 5, we can see that there is no apparent free energy barrier between the state C and the misfolded state M, so that the peptide would easily reach this state. By applying the hydrogen bond analysis to all snapshots that are in states C and M, it is found that the state M has 43% more non-native backbone hydrogen bonds than the state C. (These results indicate that the peptide could be easily trapped in the state M if normal MD simulations are employed. Strong non-native backbone hydrogen bonds would restrict the peptide from

Figure 7. The representative structures for the unfolded state (U), the compact state (C), the misfolded state (M), the partially folded state (P), and the native state (N). Arrow lines show the possible folding pathways.

escaping the state M. However, in accelerated MD simulations, the peptide would not trap in this basin since the bias potential easily breaks the non-native hydrogen bonds interactions.) The representative structures of states C, P, N, and M and the outline of a typical folding pathway based on the thermodynamics of the simulations are shown in Figure 7. From the extended structure U, the peptide first collapses to the compact state C mainly through hydrophobic interactions. Then the peptide takes the path to either the state M or the state P. And there is a path with low barrier from the state M to the state P. Finally, the formation of the native backbone hydrogen bonds drives the peptide to fold to the native state N. On the basis of the analysis of the folding landscape and several important states, we thus conclude that the folding is initiated by the hydrophobic collapse. The following processes are the formation of the native backbone hydrogen bonds. However, a question remains as in what order hydrogen bonds are formed? To answer such a question, the conformations in the transition state (TS) ensemble are investigated. According to the free energy surfaces in Figures 5 and 6, the transition region is defined as CR-rmsd > 2.0 Å and nHB ) 2 or 3. By monitoring the trajectories of the accelerated MD simulations, we search for the snapshots that satisfy the criteria above. There are 1157 snapshots located in the transition region. Through a hierarchical clustering analysis, the 1157 snapshots are divided into two clusters: the larger cluster TS1 consists of 1112 snapshots, which is 95.4% of the transition state ensemble (the states have been reweighted to the Boltzmann distribution of the original system without the bias potential), and the smaller one (TS2) only contains 45 snapshots, consisting 4.6% of the transition state ensemble. The representative structures of TS1 and TS2 are shown in Figure 8. For the snapshots in the TS1 ensemble, the two inner native backbone hydrogen bonds HB1 and HB2 are formed (see Figure

Thermodynamics and Folding Pathways of Trpzip2

J. Phys. Chem. B, Vol. 113, No. 3, 2009 807

Figure 8. The representative structures for the transition state ensembles: TS1 and TS2. The native folded structure with 5 native backbone hydrogen bonds labeled is also shown for comparison.

8 for the definitions of HB1-HB5). The HB3 forms in about half of the snapshots in the TS1, whereas no HB4 and HB5 backbone hydrogen bonds are found. This intermediate structure clearly demonstrates the early formation of the turn and indicates a hydrogen bond-centric “zip-out” folding mechanism.22,39 As Munoz et al. mentioned in their work,40 zip-out is just the most probable way to form β-hairpin structures, and other mechanisms could also exist. In the transition state ensemble TS2, the two inner hydrogen bonds HB1 and HB2 are absent, whereas the outer hydrogen bonds HB4 and HB5 are kept intact. This would be consistent with the “zip-in” mechanism,41 in which backbone hydrogen bonds fold from the terminal. Although structures that are consistent with both zip-out and zip-in (less favorable) mechanisms have been observed, detailed kinetic studies are needed to quantitatively determine their importance in the folding process. Estimation of Folding Rate. On the basis of the accelerated MD simulations of folding/unfolding of trpzip2, thermodynamics properties are easily obtained. However, the kinetics can not be obtained directly from accelerated MD simulations. A combined use of enhanced sampling simulations (for fast free energy calculations) and normal MD simulations (for dynamics properties) is a feasible way to get kinetics of the system. In this study, the protein folding time can be estimated from the reaction rate theory.42 Kramers’ theory of unimolecular reaction rates in solution43-45 assumes that the dynamics can be described by one-dimensional diffusion along a reaction coordinate in which both the reactant well (in this case, the unfolded well) and the barrier top are parabolic. The folding time is given by eq 4,42

τfolding )

( )

( )

2πkBT 2πkBT ∆G* ∆G* exp ≈ 2 exp ωminωmaxDmax kBT kBT ωminDmin (4)

where ωmin and ωmax are frequencies that characterize the curvature of the free energy profile at the unfolded well and the barrier top, respectively; Dmin and Dmax are the diffusion constants at the unfolded well and the barrier top; ∆G* is the height of the free energy barrier; kB is Boltzmann’s constant; and T is the absolute temperature. To estimate the folding time of trpzip2, two reaction coordinates, CR-rmsd and number of native contact, are used. The calculated free energy profiles as a function of CR-rmsd or number of native contact based on the four 500 ns accelerated MD simulations are shown in Figure 9. In Figure 9, the height of the free energy barrier ∆G* along CR-rmsd and number of native contact is 1.9 and 1.4 kcal/mol,

Figure 9. The free energy profile of trpzip2 using CR-rmsd and number of native contacts as the reaction coordinates.

respectively. The frequency at the unfolded well ωmin is 0.21 s-1 and 0.08 s-1, respectively. To obtain the diffusion constant at the unfolded well, normal MD simulations starting from the unfolded well need to be considered. Because the absence of solvent friction in the implicit solvent model may overestimate the diffusion constant, a 5 ns normal MD simulation with explicit solvent was conducted. First the structure at the unfolded well obtained in accelerated MD simulations was solvated with TIP3P water molecules. Then the system was heated up to 300 K. Finally, a 5 ns normal MD simulation was run for the diffusion constant calculation. Starting from the unfolded wells of Figure 9, panels a and b, we obtained the diffusion constant Dmin ) 5.0 × 10-12 m2/s and 8.1 × 10-12 m2/s. Using eq 1, the folding time is estimated to be 3.0 and 5.0 µs respectively. Both values are close to the experiment result of 2.5 µs.39 Discussions and Conclusions In this paper, we studied the folding of a β-hairpin trpzip2 by accelerated MD simulations. Multiple folding and unfolding events can be observed in the accelerated MD simulations with a total length of 2 µs. The acquirement of both folding and unfolding trajectories makes the mapping of the thorough free energy landscape more efficiently since important regions of conformational space involved in the folding and unfolding processes can be reached. By constructing two-dimensional profiles from the simulations using the number of native backbone hydrogen bonds, nHB, the radius of gyration of the hydrophobic core, Rgcore, and the CR-rmsd as the reaction coordinates, the folding landscape and a favorable folding/ unfolding pathway were obtained. It was found that both the hydrophobic interactions and hydrogen bonds play important roles in the folding process. The comparison between two transition state ensembles suggested that the dominant formation order of backbone hydrogen bonds is from inner to terminal in the hairpin folding process, which indicates a zip-out folding mechanism. We further show by this study that the accelerated MD simulation method makes the quick sampling of the energy and configuration spaces feasible, so that it is very helpful in mapping the folding landscape and calculating thermodynamic properties. A series of accelerated MD simulations for small β-hairpin trpzip4, GB1, β-sheet DPDP, R-helix Villin headpiece, and R,β-BBA5 proteins have all yielded successful folding/ unfolding trajectories. Further simulations are being performed to collect enough data for the analysis of folding mechanism as well as calculation of thermodynamics for these systems. These studies strongly show the bright prospect of the method.

808 J. Phys. Chem. B, Vol. 113, No. 3, 2009 Similar to many other enhanced sampling methods, one should note that kinetics is lost in this method. One way to obtain kinetics, as shown in the paper, is to use Kramer’s theory based on both the free energy profiles calculated from accelerated MD simulations and the diffusion constant calculated from normal MD simulations. We used the Kramer’s equation instead of shooting short trajectories46-48 because of the following consideration: the pre-exponential factor in the rate expression of this protein folding process is on the order of a fraction of microseconds (instead of the sub-picoseconds for bond-making and -breaking chemical reactions), and this small pre-exponential factor corresponds to a long “molecular time scale” (again on the order of microsecond) as defined in ref 49. The successful probability of a short trajectory that is in the picoseconds or even nanoseconds range is thus not guaranteed to follow the exponential distribution when the molecular time scale is in the sub-microseconds range and would likely contain limited kinetic information for protein folding, although picoseconds trajectories are known to be good for chemical reactions.49-51 The difference between typical chemical reactions and the more diffusive protein folding process can also be seen from the barrier height, which is typically much lower for the latter (e.g., it is only a couple of kcal/mol for the protein under study). As a result, even if starting from the transition state, since one gains only a very small factor by overcoming the barrier (for Trpzip2, a factor of exp(3.7) ≈ 50), it takes a rather long time for the trajectory to reach the reactant or product states. The physical meaning of short trajectories shot from the barrier for the folding of Trpzip2 would thus be debatable. Acknowledgment. This work is partially supported by the ACS-PRF and the Camille and Henry Dreyfus Foundation. Y.Q.G. is a 2006 Searle Scholar. We thank all the reviewers for their excellent advice. References and Notes (1) Anfinsen, C. B. Science 1973, 181, 223. (2) Levintha, C. J. Chim. Phys. Phys., Chim. Biol 1968, 65, 44. (3) Karplus, M.; Weaver, D. L. Nature 1976, 260, 404. (4) Dobson, C. M.; Sali, A.; Karplus, M. Angew Chem. Int. Ed. 1998, 37, 868. (5) Bryngelson, J. D.; Wolynes, P. G Proc. Natl. Acad. Sci. USA 1987, 84, 7524. (6) Karplus, M.; Sali, A. Curr. Opin. Struc. Biol. 1995, 5, 58. (7) Dinner, A. R.; Sali, A.; Smith, L. J.; Dobson, C. M.; Karplus, M. Trends Biochem. Sci. 2000, 25, 331. (8) Krivov, S. V.; Karplus, M. Proc. Natl. Acad. Sci. USA 2004, 101, 14766. (9) Dinner, A. R.; Lazaridis, T.; Karplus, M. Proc. Natl. Acad. Sci. USA 1999, 96, 9068. (10) Duan, Y.; Kollman, P. A. Science 1998, 282, 740. (11) Zagrovic, B.; Snow, C. D.; Shirts, M. R.; Pande, V. S. J. Mol. Biol. 2002, 324, 1051. (12) Pande, V. S.; Baker, I.; Chapman, J.; Elmer, S. P.; Khaliq, S.; Larson, S. M.; Rhee, Y. M.; Shirts, M. R.; Snow, C. D.; Sorin, E. J.; Zagrovic, B. Biopolymers 2003, 68, 91. (13) Mark, A. E.; Vangunsteren, W. F. Biochemistry 1992, 31, 7745.

Yang et al. (14) Tiradorives, J.; Jorgensen, W. L. Biochemistry 1993, 32, 4175. (15) Daggett, V.; Levitt, M. Curr. Opin. Struc. Biol. 1994, 4, 291. (16) Daggett, V.; Li, A. J.; Itzhaki, L. S.; Otzen, D. E.; Fersht, A. R. J. Mol. Biol. 1996, 257, 430. (17) Shea, J. E.; Brooks, C. L. Annu. ReV. Phys. Chem. 2001, 52, 499. (18) Yang, L. J.; Grubb, M. P.; Gao, Y. Q. J. Chem. Phys. 2007, 126, 125102. (19) Gao, Y. Q.; Yang, L. J. J. Chem. Phys. 2006, 125, 114403. (20) Cochran, A. G.; Skelton, N. J.; Starovasnik, M. A. Proc. Natl. Acad. Sci. USA 2002, 99, 9081. (21) Yang, W. Y.; Pitera, J. W.; Swope, W. C.; Gruebele, M. J. Mol. Biol. 2004, 336, 241. (22) Zhang, J.; Qin, M.; Wang, W. Proteins 2006, 62, 672. (23) Roitberg, A. E.; Okur, A.; Simmerling, C. J. Phys. Chem. B 2007, 111, 2415. (24) Pitera, J. W.; Haque, I.; Swope, W. C. J. Chem. Phys. 2006, 124, 141102. (25) Torrie, G. M.; Valleau, J. P. J. Comput. Phys. 1977, 23, 187. (26) Hamelberg, D.; Mongan, J.; McCammon, J. A. J. Chem. Phys. 2004, 120, 11919. (27) Hamelberg, D.; Shen, T.; McCammon, J. A. J. Chem. Phys. 2005, 122, 241103. (28) Hamelberg, D.; Shen, T.; McCammon, J. A. J. Am. Chem. Soc. 2005, 127, 1969. (29) Shen, T. Y.; Hamelberg, D.; McCammon, J. A. Phys. ReV. E 2006, 73, 041908. (30) Gao, Y. Q. J. Chem. Phys. 2008, 128, 064105. (31) Gao, Y. Q. J. Chem. Phys. 2008, 128, 134111. (32) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 5179. (33) Still, W. C.; Tempczyk, A.; Hawley, R. C.; Hendrickson, T. J. Am. Chem. Soc. 1990, 112, 6127. (34) Qiu, D.; Shenkin, P. S.; Hollinger, F. P.; Still, W. C. J. Phys. Chem. A 1997, 101, 3005. (35) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327. (36) Berendsen, H. J. C.; Postma, J. P. M.; Vangunsteren, W. F.; Dinola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684. (37) Snow, C. D.; Qiu, L. L.; Du, D. G.; Gai, F.; Hagen, S. J.; Pande, V. S Proc. Natl. Acad. Sci. USA 2004, 101, 4077. (38) Feig, M.; Karanicolas, J.; Brooks, C. L. J. Mol. Graphics Modell. 2004, 22, 377. (39) Du, D. G.; Zhu, Y. J.; Huang, C. Y.; Gai, F. Proc. Natl. Acad. Sci. USA 2004, 101, 15915. (40) Munoz, V.; Thompson, P. A.; Hofrichter, J.; Eaton, W. A. Nature 1997, 390, 196. (41) Munoz, V.; Ghirlando, R.; Blanco, F. J.; Jas, G. S.; Hofrichter, J.; Eaton, W. A. Biochemistry 2006, 45, 7023. (42) Kubelka, J.; Hofrichter, J.; Eaton, W. A. Curr. Opin. Struc. Biol. 2004, 14, 76. (43) Kramers, H. A. Physica 1940, 7, 284. (44) Berne, B. J.; Borkovec, M.; Straub, J. E. J. Phys. Chem. 1988, 92, 3711. (45) Socci, N. D.; Onuchic, J. N.; Wolynes, P. G. J. Chem. Phys. 1996, 104, 5860. (46) Bolhuis, P. G. Proc. Natl. Acad. Sci. USA 2003, 100, 12129. (47) Juraszek, J.; Bolhuis, P. G. Proc. Natl. Acad. Sci. USA 2006, 103, 15859. (48) Snow, C. D.; Rhee, Y. M.; Pande, V. S Biophys. J. 2006, 91, 14. (49) Chandler, D. Introduction to Modern Statistical Mechanics; Oxford University Press: New York, 1987. (50) Dellago, C.; Bolhuis, P. G.; Csajka, F. S.; Chandler, D. J. Chem. Phys. 1998, 108, 1964. (51) Fu, X. B.; Yang, L. J.; Gao, Y. Q. J. Chem. Phys. 2007, 127, 154106.

JP803160F