- Email: [email protected]
Ab initio Folding Simulation of the Trp-cage Mini-protein Approaches NMR Resolution
doi:10.1016/S0022-2836(03)00177-3
J. Mol. Biol. (2003) 327, 711–717
Ab initio Folding Simulation of the Trp-cage Mini-protein Approaches NMR Resolution Shibasish Chowdhury, Mathew C. Lee, Guoming Xiong and Yong Duan* Department of Chemistry and Biochemistry, Center of Biomedical Research Excellence in Structural and Functional Genomics, University of Delaware, Newark, DE 19716 USA
Here, we report a 100 ns molecular dynamics simulation of the folding process of a recently designed autonomous-folding mini-protein designated as tc5b with a new AMBER force field parameter set developed based on condensed-phase quantum mechanical calculations and a Generalized Born continuum solvent model. Starting from its fully extended conformation, our simulation has produced a final structure resembling ˚ main-chain root mean square that of NMR native structure to within 1 A deviation. Remarkably, the simulated structure stayed in the native state for most part of the simulation after it reached the state. Of greater significance is that our simulation has not only reached the correct main-chain conformation, but also a very high degree of accuracy in side-chain packing conformation. This feat has traditionally been a challenge for ab initio simulation studies. In addition to characterization of the trajectory, comparison of our results to experimental data is also presented. Analysis of the trajectory suggests that the rate-limiting step of folding of this miniprotein is the packing of the Trp side-chain. q 2003 Elsevier Science Ltd. All rights reserved
*Corresponding author
Keywords: ab initio protein folding; Trp-cage; molecular dynamics; protein design; AMBER
Introduction Despite decades of intensive research, the protein folding problem remains a major challenge to both experimentalists and theoreticians. However, in recent years, exciting advances in this field have been made by both experimental and computational biologists. On the experimental front, artificially designed autonomous-folding mini-proteins have been achieved in several laboratories.1 – 4 Examples include the 28-residue bba motif by Dahiyat & Mayo,1 the 23-residue bba mini-protein by Imperiali and co-workers,3 and the latest entrant of a 20-residue protein that folds into a novel Trp-cage fold.5 These experiments have helped us to address the fundamental questions regarding the minimum requirements for a protein sequence to fold; they have generated model systems to elucidate the general principles governing protein folding. On the computational front, increase in computer Abbreviations used: RMSD, root mean square deviation. E-mail address of the corresponding author: [email protected]
speed and improvements in force field along with more efficient computation algorithms have brought realistic computer simulation of the folding/unfolding process to within the realm of possibility. Some notable examples are the folding/unfolding simulation of chymotrypsin inhibitor 2 by Daggett and co-workers,6 the folding simulation of villin headpiece sub-domain by Duan & Kollman7 and the recent work by Pande and co-workers.8 These computer simulations have given us a microscopic view of the important molecular events in atomic details and have complemented the macroscopic views of the folding process based on experiments and low-resolution simulations. In this short article, we report our ab initio folding simulation of the aforementioned 20-residue Trp-cage mini-protein, the smallest protein-like model system presently known. In the paragraphs below, we briefly describe the structural features which make this protein a unique model system, followed by characterization of our simulation and comparisons to experimental results. Using Trp-cage sub-domain of Exe-4 as a starting point, Andersen and co-workers performed a series of mutagenesis and truncation experiments
0022-2836/03/$ - see front matter q 2003 Elsevier Science Ltd. All rights reserved
712
Ab initio Folding of the Trp-cage Mini-protein
Figure 1. (A) Native NMR structure of tc5b. (B) Representative simulated structure of tc5b. Backbone is shown in red ribbon. The Trp-cage forming residues, Trp25 (yellow), Pro31 (green), Pro37 (magenta) and Pro38 (blue) are shown in CPK model. All other side-chains are shown in stick.
that finally led to a stable mini-protein designated as tc5b5 whose Trp side-chain was shown to be 100% buried in aqueous solution by NMR measurements at 280 K. We will follow the numbering scheme of Anderson and co-workers in which the primary sequence of tc5b is numbered as N20LYIQW25LKDGG30PSSGR35PPPS39. NMR structure (1L2Y) of this Trp-cage motif shows important features of a small protein (Figure 1(A)); it has multiple secondary structural elements, side-chain – side-chain packing interactions, well defined tertiary contacts and a compact overall structure. Starting from the Nterminal, the sequence folds into one a-helix spanning residues 21 – 27 followed by a 310-helix (residues 30 – 33) and a short stretch of coil rigidified by three consecutive proline residues extending to the C-terminal. In addition to the secondary structural elements, hydrogen bonding between NH of Gly30 and the backbone carbonyl oxygen of Trp25, NH11 of the Trp25 indole ring to the backbone carbonyl of Arg35 are both important stabilizing factors that characterize the NMR structure. Side-chain rotamer restriction observed in the ensemble of NMR structure is a strong evidence that tc5b is a highly compact globular protein with little conformation flexibility. The hydrophobic core that makes up the Trp-cage (two proline residues on each side of the Trp indole ring plus Tyr22) is thought to be the driving force for the folding process of tc5b. Taken together, the small size of the protein, the architectural simplicity and the limited conformational flexibility make this mini-protein an ideal model system for computational study.
Results and Discussions Started from a straight-chain conformation, the protein rapidly collapsed into a meta-stable state and the main-chain root mean square deviation ˚ (RMSD) decreased from an initial value of , 10 A ˚ within 10 ps of simulation time. After this to , 4 A initial transition, the main-chain RMSD persisted ˚ for about 23 ns. In the mean time, at around 4 A the heavy atom RMSD decreased in a staircase ˚ . At around 23 ns of simufashion to around 4.5 A lation time, the meta-stable structure underwent a series of major conformational changes to reach a stable state that closely resembles the NMR structure at around 29 ns. At this point the main-chain ˚ and the heavy-atom RMSD was around 1.0 A ˚ RMSD reached 2 A (Figure 2(A)). Once the structure fell into the native basin, it stayed there for the remainder of the simulation. This is remarkable. The simulation clearly showed that the native state is the most stable state sampled during the simulation. In comparison to the ensemble of NMR structures, the heavy atom RMSD for the last 60 ns of ˚ with a standard deviation the simulation was , 2 A ˚ , whereas the pair-wise RMSD of the 38 of 0.2 A ˚ to models in the NMR ensemble ranges from 1.6 A ˚ ˚ 2.8 A with a standard deviation of 0.3 A. These statistics place the quality of the simulated structures to be well within the range of experimental uncertainties, further highlighting the excellent agreement between our simulation and experimental observations. Closer examination of the simulated structure revealed that the same cage pattern was observed
Ab initio Folding of the Trp-cage Mini-protein
Figure 2. (A) RMSD profile of main-chain atoms (blue) and all heavy atoms (red) with respect to corresponding atoms in the representative NMR structure. Residues from 3 to 18 are considered for the RMSD calculation. If we consider full protein the average RMSD is increased ˚ . (B) Fraction of NOE violation (PDB code 1L2Y) by 0.5 A is plotted during the entire folding simulation. (C) Number of native atom to atom contacts of Trp25 side-chain with all other atoms was plotted. Atoms in the Trp25 side-chain are considered to be in contact with other ˚. atoms if the distance between two atoms is within 3 A (D) Fractional residual native contact is plotted. Time in nanoseconds is plotted along x-axis in all four plots.
in our simulation as that of the NMR structure in the Trp-cage region; the landmark residue Trp25 was securely sandwiched by two proline rings on either side of its indole ring, while Tyr22 stood nearby to close up the cage (Figure 1(B)). Furthermore, we also observed the same characteristic hydrogen bonding pattern between NH of Gly30 and the backbone carbonyl oxygen of Trp25. These comparisons reinforce our conclusion that our simulation has indeed correctly folded tc5b from straight chain to its experimentally demonstrated native conformation. The experimental native structure ensemble of this mini-protein was derived from 169 NOE distances. As a rigorous benchmark of our results, we compared our simulated structures with these NOE distances (obtained from the PDB) in a timeseries plot of NOE violations (Figure 2(B)). The simulated structures showed an initial 60% NOE violation in its initial extended state which decreased steadily and reached 48% at around 2 ns. After a relatively sharp decrease in NOE violations between 2 ns and 4 ns, we observed a slow phase which lasted for about 20 ns (up to 23 ns). The NOE violations then decreased again sharply and fell below 20% by 30 ns of simulation time. In the final 1 ns of simulation, only , 14% NOE violations were observed. Detailed analysis revealed that NOE violations between Pro37 and Trp25 sidechains contribute the major part of the overall NOE violation (, 6%). NOE violations between Trp25 and long flexible side-chains of Leu26 as well as Arg35 also contribute a notable fraction (, 3%)
713
of overall NOE violation. However, if we increase ˚ , most of the NOE distance restraints by 0.4 A these NOE violations disappear and the average NOE violation fell down to 3– 7%. In contrast, the ensemble of NMR models possesses an average of 3.6% NOE distance violations. Thus, the simulated structure closely resembles the native NMR structures. Considering the dynamic nature, this result should be viewed as being physically reasonable, and lends a high degree of rigor to our simulated structures in terms of benchmarking simulation against experimental observation. The percentage of native contacts over the simulation is shown in Figure 2(D). Not surprisingly, the development of the native contacts correlated closely to the drop of NOE distance violation. In our analysis, residues were considered in contact ˚, if any of the atom pairs were closer than 4.0 A starting from the residues between i and i þ 3. With this definition of native contact, we identified on average of about 61(^ 5) native contacts within the NMR structure ensemble. The average fraction of native contacts in the simulation rose quickly to 45% within the first nanosecond. After that there was a slow phase of native contact formation (1 – 15 ns), during which the fraction of native contact decreased slightly. Subsequently, fractional native contact increased sharply to around , 85% within the next 15 ns of simulation with a brief slow phase during , 18 –22 ns. An average of 90% native contact was observed in the last 20 ns of simulation, indicating that the folding simulation generated structures with almost all the native contacts. The high percentage of native contacts reinforces the notion that our simulated trajectory captured some of the main aspects of the folding events. Hydrogen bonding is arguably one of the most important weak interactions that determine the formation of secondary structural elements in proteins. In our simulation, we observed three to four a-helical i, i þ 4 main-chain hydrogen bonds within the short a-helical fragment in the last 20 ns of simulation; this bonding pattern is a characteristic feature observed in the ensemble of the NMR structures. Formation of the first hydrogen bond in the a-helical region of the simulated protein was complete by 1.7 ns. On average, a total of , 1.5 main-chain hydrogen bonds were formed between 11 ns and 12 ns. By 25 ns of simulation time almost all main-chain 4 ! 1 hydrogen bonds were formed within the a-helical segment. The formation of hydrogen bond analysis of the ahelical segment closely resembled that of helicity analysis based on F and C angles (data not shown). Around 0.3 hydrogen bonds formed at 25 ns within the short 310-helix segment. By 29 ns of simulation, formation of the main-chain hydrogen bond within the 310-helical region was completed. In short, the chronological order of hydrogen bond formation showed that formation of the 310-helix initiated after the formation of ahelical segment and the formation of all helices
714
was completed around 29 ns, which was roughly the folding time of the entire protein (judged by the RMSD value with respect to native structure). The ability to form both a- and 310-helices is noteworthy given the on-going debate on the relative significance of a- versus 310-helices,11 indicating a reasonable balance between these two conformations in the force field. Close inspection of the folding events indicates that formation of the 310helix was preceded by a short a-helix. Thus, contrary to the earlier suggestions,12 the 310-helix is not a transition state of a-helix. Perhaps one of the most interesting aspects of molecular dynamics simulation is that one can track the time evolution of structural changes at high temporal resolutions, a feat difficult to do with current experimental methods. In our analysis, particular attention was paid to the changes in the native contacts made by residue Trp25 as it constitutes the landmark feature of the mini-protein tc5b where its side-chain is buried by three nearby proline residues and one tyrosine. Visualization of the folding movie showed that the Trp side-chain interacted with Pro36, 37 and 38 in the early stage of simulation (within 1 ns of simulation) which stabilized the native topology of the protein. The hydrophobic contacts between Trp25 and Leu21 and Ile23 were also noted in the early stage of simulation. Thus, the first stage of folding consisted of the formation of native-like backbone topology, stabilized by the clustering of hydrophobic residues around residue Trp25. However, contacts between Trp25 and Pro31 as well as between Asp28 and Arg35 (potential salt bridge) were not observed in the early stage of simulation. The second milestone in the folding process was the formation of secondary structures. During this stage, the a-helical contact between Trp25 and Tyr22 began to build up after 1 ns of simulation, signifying that the a-helix was initiated from the N terminus. This point was further supported by the early contacts made between Gln24 and Leu21. In contrast, the contacts between Gln24 and Lys27 were established at around 15 ns, completing the formation of the C terminus end of the a-helix. The stable salt bridge between Asp28 and Arg35 was also observed after 15 ns of simulation. The fact that salt bridge formation occurred at this later stage of folding alludes to our initial hypothesis that the hydrophobic force was the main driving force for the native protein topology, whereas the salt bridges provided additional stability to the native structure. Different (structural) states in the folding pathway were examined by clustering analysis. This analysis would inform us about the different meta-stable states within the initial straight chain conformation and final native state. It would also tell us about the relative population of each detectable state which is a measure of the relative stability of each state. A total of 19 clusters were identified by the hierarchical clustering approach, however, only nine clusters were populated with
Ab initio Folding of the Trp-cage Mini-protein
more than 50 snapshots. Visualization of these nine clusters further subdivided the clusters into four intermediate states and one native state, which correspond closely to the development in main-chain RMSD. The first intermediate state represents the initial transient state which was formed by the initial hydrophobic collapse of the straight chain structure. It was stable for no more than 1 ns. The second intermediate state represents structures that exhibited partial helical turn with Leu2-carbonyl to Trp6-NH and Tyr3-carbonyl to Leu7-NH hydrogen bonds, another evidence that helix formation was initiated from the N terminus. About 13% of all saved structures fell into this state. Although the overall topology of this intermediate state differed significantly from the native structure, the hydrophobic contacts between Trp25 sidechain and Pro36-38 were present. The third intermediate state, having a population of 10% of the total saved snapshots, also contained structures with partial helix. However, the overall topology of these structures was much closer to the native structure. The fourth intermediate state, consisting of 4% of the saved snapshots, is characterized by the presence of the completed a-helical segment, although the small 310-helical segment was still absent. Finally, the native state structure, with ahelical and short 310-helical segment completed, is the most populated state with about 72% of all saved structures. Superimposition of five structures randomly taken from the native state cluster to that of the NMR ensemble is shown in Figure 3. One can clearly see from the stereographic Figure that the two ensembles of structures coincided very closely with each other. Most notably, the side-chain packing in this state is almost identical to the native NMR structure. Our simulation results indicate that the rate-limiting step is not the initiation or formation of ahelical segment. We observed that the helix could form its first turn rapidly within nanosecond time scale and completed within 15 ns. There is at least 20 ns time lag between the initiation of helix and completion of protein folding. In our opinion, disparity between these two important events rules out the possibility that the helix initiation played any significant role in the folding of tc5b. Our data suggest that final sharp drop of the RMSD values which brings tc5b structure into the vicinity of the native structure coincide with the sharp increase in native contacts or packing of the Trp25 side-chain. The time series of heavy atom RMSD plot (Figure 2(A)) clearly showed four stages of conformational changes as there are four plateaus in this time series. In contrast, there are only three plateaus in the atom to atom native contact time series plot of Trp25 (Figure 2(C)). The first, third and the fourth major transitions in heavy atom RMSD all correlated to a corresponding change in Trp25 native contacts. The final descent in RMSD value (at , 29 ns) correlated with a sharp increase in native contacts of Trp25 which
Ab initio Folding of the Trp-cage Mini-protein
715
Figure 3. Stereo view of the five representative NMR structures (red) superimposed on the five simulated snap shots (blue). Simulated snap shots are taken from the last 5 ns of simulation at 1 ns interval. Only heavy atoms are shown for clarity.
represents the last step of conformational rearrangement involving correct packing of Trp25 side-chain that led to the formation of native tc5b structure. This can be visualized from the folding movie. Further detailed contact analysis between side-chain atoms of Trp25 and individual atom of each other residue revealed that contacts between Trp25 side-chain and side-chains of Tyr22, Gly30, Pro31, Pro37 and Pro38 are mainly responsible for the final sharp drop in heavy atom RMSD. Interestingly these are the residues that essentially entrap the Trp25 side-chain. Thus, the last step of tc5b folding could be the correct packing of Trp side-chain. Although unequivocal assignment of the ratelimiting step requires more detailed thermodynamic analysis, based on our observation, we believe that the correct packing of Trp25 side-chain is most likely the rate-limiting step where upon packing, complete folding of the entire protein immediately followed. Indeed, in our subsequent simulations, packing of Trp side-chains was seen as the ratelimiting step, despite the quick formation of both secondary and tertiary structures (detail will be reported elsewhere). Nevertheless, one possible experimental validation could be measurement of helix formation in comparison to the packing of fluorescent Trp residues. Our simulation suggests that helix formation precedes completion of the folding. Identification of the rate-limiting step and the nature of transition state has been the focus in the field of protein folding mechanism studies. Presently, an important theory attributes the transition state to the chain entropy.13 On the other hand, Fvalue analyses14 have suggested a much more complex picture. We think our results agree with the latter better simply because formation of the native-like main-chain structure was early in the process and did not appear to be the rate-limiting step. In our opinion, this is not necessarily contrary to the conclusion of Plaxco et al., instead it is comp-
lementary. In the study of Plaxco et al. folding rates of multiple small proteins were examined. Given the diversity of the protein structures, a common denominator would the contribution of the chain entropy and the specific interactions that play roles in the folding of individual proteins would appear as the “noise” and uncertainties. Therefore, when one asks the common factor affecting folding of a diverse pool of proteins, chain-entropy would play an important role. However, in the folding of individual protein, specific interactions, such as those identified by F-value analysis and those described in this study, would be needed to describe the process. It should be noted that our present study focuses on the folding events observed in one simulation trajectory. Despite its high degree of realism, we must keep in mind that alternative folding processes are also likely to occur given the complexity of the free energy landscape. Therefore, we should focus on the qualitative aspects of the simulation results, for example, the sequence of events. The time scales of the events should be obtained by multiple simulations. This is particularly true since we employed a continuum solvent model in our simulation which is known to accelerate the process due to lack of solvent viscosity. Thus our observed time scales are likely to be much faster than what happen in reality. Nevertheless, by focusing on the qualitative aspects and the sequence of events, a good deal can be learned, though much of these remain to be tested by large number of simulations. Simmerling et al. recently studied the same protein with a similar approach15 in which a standard AMBER force field was used with a set of adjusted main-chain torsion parameters. In their simulation, the heavy atom RMSD also showed a step-wise reduction. The ability to fold the same protein with different force fields clearly demonstrates the robustness of all-atom molecular dynamics simulations. Interestingly, their observed folding time
716
was around 7 ns, much shorter than what we observed in our simulation. Now, we have gone a step further. Our simulation suggests that the ratelimiting step could be packing of the Trp25 sidechain and formation of the main-chain tertiary contacts could be an early event in the folding process. Qiu et al. recently studied folding of the protein16 by laser temperature jump experiments. They found that the extrapolated folding time is about 4 ms by monitoring the fluorescence of Trp25 sidechain which measures its packing. In our opinion, this is consistent with our observation that packing of Trp25 side-chain is the rate-limiting step. The difference between the simulated and experimental folding times is likely due to the absence of solvent viscosity in the simulation and the time was extrapolated from single simulation. Our observed formation time of helical secondary structure was on the order of 25 ns, close to experimental observations on isolated helices. Thus, the a-helix may form early. Our simulation also showed that tertiary contacts can form quickly. When tertiary contacts form, however, the Trp25 formed a hydrophobic cluster outside the cage and the cage was in a closed state, both of which prevented the Trp from entering the cage quickly. Thus, breaking the hydrophobic cluster, opening up the cage, and encaging Trp are all parts of the process that eventually leads to the packing of the Trp side-chain. These steps could be much slower than what we observed due to the lack of solvent viscosity in the simulation.
Concluding Remarks Present folding simulation of tc5b explores the folding pathway of 20-residue mini-protein. The main-chain as well as side-chain configuration of the simulated structure is in excellent agreement with the NMR ensembles. The formation of Trpcage was also closely monitored, which reveals that hydrophobic interaction between Trp and Pro residues initiates the folding process by forming the correct protein topology. Initial data indicate that main-chain conformation and secondary structure formation is not the rate-limiting step; rather, correct packing of Trp25 side-chain may be the rate-limiting step of the protein folding. Additional folding as well as unfolding simulations (work in progress) may give a more accurate account of the folding kinetics of this protein.
Method In our folding simulation with AMBER package, we began from a fully extended conformation of tc5b. Using the Generalized Born implicit solvent model9 to treat solvation effects and a new all-atom point-charge force field recently developed in our laboratory, we performed 100 ns of molecular dynamics simulation to fold the protein ab initio. The detail of the force field will be published elsewhere. In summary, the charges were
Ab initio Folding of the Trp-cage Mini-protein
obtained by fitting to the electrostatic potential calculated quantum mechanically in organic solvent ð1 ¼ 4Þ: The main-chain torsion parameters were obtained by fitting to the energy surface of alanine dipeptide. The energy surface was also calculated quantum mechanically in organic solvent ð1 ¼ 4Þ: Other parameters, including bonds, bond angles, torsion parameters (except the main-chain torsion angles), and van der Waals parameters were identical to those described by Wang et al.10 The initial extended chain conformation was constructed from the primary sequence using leap module in AMBER package. After 500 steps of energy minimization, random velocities were assigned at 100 K. The simulation was allowed to continue for 100 ns with an integration time step of 2.0 fs. The temperature was controlled at 300 K by loosely coupling the system to a heat bath with a coupling constant of 2.0 ps. Non-bonded interactions were calculated without cutoff (in both MD and minimization). The trajectory was saved at 10 ps intervals for further analysis. The simulation was done on a dual processor (Pentium III 1 GHz) PC work station which gave a throughput of about 5 ns a day. The entire simulation took roughly 20 days to complete.
Acknowledgements Computer time was provided by Pittsburgh Supercomputer Center. This work was supported by research grants from NIH (RR15588, PI Lenhoff and GM64458 to Y.D.), the state of Delaware, and University of Delaware Research Fund. Usage of VMD, UCSF Midas, WebLab viewer graphics packages are gratefully acknowledged.
References 1. Dahiyat, B. I. & Mayo, S. L. (1997). De novo protein design: fully automated sequence selection. Science, 278, 82 – 87. 2. Hill, R. B. & DeGrado, W. F. (1998). Solutions structure of a2D, a nativelike de novo designed protein. J. Am. Chem. Soc. 120, 1138– 1145. 3. Struthers, M. D., Cheng, R. P. & Imperiali, B. (1996). Design of a monomeric 23-residue polypeptide with defined tertiary structure. Science, 271, 342– 345. 4. Kortemme, T., R, A. M. & Serrano, L. (1998). Design of a 20-amino acid, three-stranded b-sheet protein. Science, 281, 253– 256. 5. Neidigh, J. W., Fesinmeyer, R. M. & Andersen, N. H. (2002). Designing a 20-residue protein. Nature Struct. Biol. 9, 425–430. 6. Li, A. & Daggett, V. (1996). Identification and characterization of the unfolding transition state of chymotrypsin inhibitor 2 by molecular dynamics simulations. J. Mol. Biol. 257, 412–429. 7. Duan, Y. & Kollman, P. A. (1998). Pathways to a protein folding intermediate observed in a 1 ms simulation in aqueous solution. Science, 282, 740– 744. 8. Snow, C. D., Nguyen, H., Pande, V. S. & Gruebele, M. (2002). Absolute comparison of simulated and experimental protein-folding dynamics computationally. Nature, in press.
Ab initio Folding of the Trp-cage Mini-protein
9. Tsui, V. & Case, D. A. (2000). Molecular dynamics simulations of nucleic acids with a generalized born solvation model. J. Am. Chem. Soc. 122, 2489– 2498. 10. Wang, J. M., Cieplak, P. & Kollman, P. A. (2000). How well does a restrained electrostatic potential (RESP) model perform in calculating conformational energies of organic and biological molecules? J. Comp. Chem. 21, 1049– 1074. 11. Tran, T. T., Zeng, J., Treutlein, H. & Burgess, A. W. (2002). Effects of thioamide substitutions on the conformation and stability of a- and 310-helices. J. Am. Chem. Soc. 124, 5222– 5230. 12. Basu, G., Kitao, A., Hirata, F. & Go, N. (1994). A collective motion description of the 310-/a-helix transition—implications for a natural reaction coordinate. J. Am. Chem. Soc. 116, 6307– 6315.
717
13. Plaxco, K. W., Simons, K. T. & Baker, D. (1998). Contact order, transition state placement and the refolding rates of single domain proteins. J. Mol. Biol. 277, 985 –994. 14. Fersht, A. R. (1995). Mapping the structures of transition-states and intermediates in folding—delineation of pathways at high-resolution. Phil. Trans. Roy. Soc. London B-Biol. Sci. 348, 11 –15. 15. Simmerling, C., Strockbine, B. & Roitberg, A. E. (2002). All-atom structure prediction and folding simulations of a stable protein. J. Am. Chem. Soc. 124, 11258 –11259. 16. Qiu, L., Pabit, S. A., Roitberg, A. E. & Hagen, S. J. (2002). Smaller and faster: the 20-residue Trp-Cage protein folds in 4 ms. J. Am. Chem. Soc. 124, 12952 – 12953.
Edited by M. Levitt (Received 10 September 2002; received in revised form 26 November 2002; accepted 23 January 2003)
J. Mol. Biol. (2003) 327, 711–717
Ab initio Folding Simulation of the Trp-cage Mini-protein Approaches NMR Resolution Shibasish Chowdhury, Mathew C. Lee, Guoming Xiong and Yong Duan* Department of Chemistry and Biochemistry, Center of Biomedical Research Excellence in Structural and Functional Genomics, University of Delaware, Newark, DE 19716 USA
Here, we report a 100 ns molecular dynamics simulation of the folding process of a recently designed autonomous-folding mini-protein designated as tc5b with a new AMBER force field parameter set developed based on condensed-phase quantum mechanical calculations and a Generalized Born continuum solvent model. Starting from its fully extended conformation, our simulation has produced a final structure resembling ˚ main-chain root mean square that of NMR native structure to within 1 A deviation. Remarkably, the simulated structure stayed in the native state for most part of the simulation after it reached the state. Of greater significance is that our simulation has not only reached the correct main-chain conformation, but also a very high degree of accuracy in side-chain packing conformation. This feat has traditionally been a challenge for ab initio simulation studies. In addition to characterization of the trajectory, comparison of our results to experimental data is also presented. Analysis of the trajectory suggests that the rate-limiting step of folding of this miniprotein is the packing of the Trp side-chain. q 2003 Elsevier Science Ltd. All rights reserved
*Corresponding author
Keywords: ab initio protein folding; Trp-cage; molecular dynamics; protein design; AMBER
Introduction Despite decades of intensive research, the protein folding problem remains a major challenge to both experimentalists and theoreticians. However, in recent years, exciting advances in this field have been made by both experimental and computational biologists. On the experimental front, artificially designed autonomous-folding mini-proteins have been achieved in several laboratories.1 – 4 Examples include the 28-residue bba motif by Dahiyat & Mayo,1 the 23-residue bba mini-protein by Imperiali and co-workers,3 and the latest entrant of a 20-residue protein that folds into a novel Trp-cage fold.5 These experiments have helped us to address the fundamental questions regarding the minimum requirements for a protein sequence to fold; they have generated model systems to elucidate the general principles governing protein folding. On the computational front, increase in computer Abbreviations used: RMSD, root mean square deviation. E-mail address of the corresponding author: [email protected]
speed and improvements in force field along with more efficient computation algorithms have brought realistic computer simulation of the folding/unfolding process to within the realm of possibility. Some notable examples are the folding/unfolding simulation of chymotrypsin inhibitor 2 by Daggett and co-workers,6 the folding simulation of villin headpiece sub-domain by Duan & Kollman7 and the recent work by Pande and co-workers.8 These computer simulations have given us a microscopic view of the important molecular events in atomic details and have complemented the macroscopic views of the folding process based on experiments and low-resolution simulations. In this short article, we report our ab initio folding simulation of the aforementioned 20-residue Trp-cage mini-protein, the smallest protein-like model system presently known. In the paragraphs below, we briefly describe the structural features which make this protein a unique model system, followed by characterization of our simulation and comparisons to experimental results. Using Trp-cage sub-domain of Exe-4 as a starting point, Andersen and co-workers performed a series of mutagenesis and truncation experiments
0022-2836/03/$ - see front matter q 2003 Elsevier Science Ltd. All rights reserved
712
Ab initio Folding of the Trp-cage Mini-protein
Figure 1. (A) Native NMR structure of tc5b. (B) Representative simulated structure of tc5b. Backbone is shown in red ribbon. The Trp-cage forming residues, Trp25 (yellow), Pro31 (green), Pro37 (magenta) and Pro38 (blue) are shown in CPK model. All other side-chains are shown in stick.
that finally led to a stable mini-protein designated as tc5b5 whose Trp side-chain was shown to be 100% buried in aqueous solution by NMR measurements at 280 K. We will follow the numbering scheme of Anderson and co-workers in which the primary sequence of tc5b is numbered as N20LYIQW25LKDGG30PSSGR35PPPS39. NMR structure (1L2Y) of this Trp-cage motif shows important features of a small protein (Figure 1(A)); it has multiple secondary structural elements, side-chain – side-chain packing interactions, well defined tertiary contacts and a compact overall structure. Starting from the Nterminal, the sequence folds into one a-helix spanning residues 21 – 27 followed by a 310-helix (residues 30 – 33) and a short stretch of coil rigidified by three consecutive proline residues extending to the C-terminal. In addition to the secondary structural elements, hydrogen bonding between NH of Gly30 and the backbone carbonyl oxygen of Trp25, NH11 of the Trp25 indole ring to the backbone carbonyl of Arg35 are both important stabilizing factors that characterize the NMR structure. Side-chain rotamer restriction observed in the ensemble of NMR structure is a strong evidence that tc5b is a highly compact globular protein with little conformation flexibility. The hydrophobic core that makes up the Trp-cage (two proline residues on each side of the Trp indole ring plus Tyr22) is thought to be the driving force for the folding process of tc5b. Taken together, the small size of the protein, the architectural simplicity and the limited conformational flexibility make this mini-protein an ideal model system for computational study.
Results and Discussions Started from a straight-chain conformation, the protein rapidly collapsed into a meta-stable state and the main-chain root mean square deviation ˚ (RMSD) decreased from an initial value of , 10 A ˚ within 10 ps of simulation time. After this to , 4 A initial transition, the main-chain RMSD persisted ˚ for about 23 ns. In the mean time, at around 4 A the heavy atom RMSD decreased in a staircase ˚ . At around 23 ns of simufashion to around 4.5 A lation time, the meta-stable structure underwent a series of major conformational changes to reach a stable state that closely resembles the NMR structure at around 29 ns. At this point the main-chain ˚ and the heavy-atom RMSD was around 1.0 A ˚ RMSD reached 2 A (Figure 2(A)). Once the structure fell into the native basin, it stayed there for the remainder of the simulation. This is remarkable. The simulation clearly showed that the native state is the most stable state sampled during the simulation. In comparison to the ensemble of NMR structures, the heavy atom RMSD for the last 60 ns of ˚ with a standard deviation the simulation was , 2 A ˚ , whereas the pair-wise RMSD of the 38 of 0.2 A ˚ to models in the NMR ensemble ranges from 1.6 A ˚ ˚ 2.8 A with a standard deviation of 0.3 A. These statistics place the quality of the simulated structures to be well within the range of experimental uncertainties, further highlighting the excellent agreement between our simulation and experimental observations. Closer examination of the simulated structure revealed that the same cage pattern was observed
Ab initio Folding of the Trp-cage Mini-protein
Figure 2. (A) RMSD profile of main-chain atoms (blue) and all heavy atoms (red) with respect to corresponding atoms in the representative NMR structure. Residues from 3 to 18 are considered for the RMSD calculation. If we consider full protein the average RMSD is increased ˚ . (B) Fraction of NOE violation (PDB code 1L2Y) by 0.5 A is plotted during the entire folding simulation. (C) Number of native atom to atom contacts of Trp25 side-chain with all other atoms was plotted. Atoms in the Trp25 side-chain are considered to be in contact with other ˚. atoms if the distance between two atoms is within 3 A (D) Fractional residual native contact is plotted. Time in nanoseconds is plotted along x-axis in all four plots.
in our simulation as that of the NMR structure in the Trp-cage region; the landmark residue Trp25 was securely sandwiched by two proline rings on either side of its indole ring, while Tyr22 stood nearby to close up the cage (Figure 1(B)). Furthermore, we also observed the same characteristic hydrogen bonding pattern between NH of Gly30 and the backbone carbonyl oxygen of Trp25. These comparisons reinforce our conclusion that our simulation has indeed correctly folded tc5b from straight chain to its experimentally demonstrated native conformation. The experimental native structure ensemble of this mini-protein was derived from 169 NOE distances. As a rigorous benchmark of our results, we compared our simulated structures with these NOE distances (obtained from the PDB) in a timeseries plot of NOE violations (Figure 2(B)). The simulated structures showed an initial 60% NOE violation in its initial extended state which decreased steadily and reached 48% at around 2 ns. After a relatively sharp decrease in NOE violations between 2 ns and 4 ns, we observed a slow phase which lasted for about 20 ns (up to 23 ns). The NOE violations then decreased again sharply and fell below 20% by 30 ns of simulation time. In the final 1 ns of simulation, only , 14% NOE violations were observed. Detailed analysis revealed that NOE violations between Pro37 and Trp25 sidechains contribute the major part of the overall NOE violation (, 6%). NOE violations between Trp25 and long flexible side-chains of Leu26 as well as Arg35 also contribute a notable fraction (, 3%)
713
of overall NOE violation. However, if we increase ˚ , most of the NOE distance restraints by 0.4 A these NOE violations disappear and the average NOE violation fell down to 3– 7%. In contrast, the ensemble of NMR models possesses an average of 3.6% NOE distance violations. Thus, the simulated structure closely resembles the native NMR structures. Considering the dynamic nature, this result should be viewed as being physically reasonable, and lends a high degree of rigor to our simulated structures in terms of benchmarking simulation against experimental observation. The percentage of native contacts over the simulation is shown in Figure 2(D). Not surprisingly, the development of the native contacts correlated closely to the drop of NOE distance violation. In our analysis, residues were considered in contact ˚, if any of the atom pairs were closer than 4.0 A starting from the residues between i and i þ 3. With this definition of native contact, we identified on average of about 61(^ 5) native contacts within the NMR structure ensemble. The average fraction of native contacts in the simulation rose quickly to 45% within the first nanosecond. After that there was a slow phase of native contact formation (1 – 15 ns), during which the fraction of native contact decreased slightly. Subsequently, fractional native contact increased sharply to around , 85% within the next 15 ns of simulation with a brief slow phase during , 18 –22 ns. An average of 90% native contact was observed in the last 20 ns of simulation, indicating that the folding simulation generated structures with almost all the native contacts. The high percentage of native contacts reinforces the notion that our simulated trajectory captured some of the main aspects of the folding events. Hydrogen bonding is arguably one of the most important weak interactions that determine the formation of secondary structural elements in proteins. In our simulation, we observed three to four a-helical i, i þ 4 main-chain hydrogen bonds within the short a-helical fragment in the last 20 ns of simulation; this bonding pattern is a characteristic feature observed in the ensemble of the NMR structures. Formation of the first hydrogen bond in the a-helical region of the simulated protein was complete by 1.7 ns. On average, a total of , 1.5 main-chain hydrogen bonds were formed between 11 ns and 12 ns. By 25 ns of simulation time almost all main-chain 4 ! 1 hydrogen bonds were formed within the a-helical segment. The formation of hydrogen bond analysis of the ahelical segment closely resembled that of helicity analysis based on F and C angles (data not shown). Around 0.3 hydrogen bonds formed at 25 ns within the short 310-helix segment. By 29 ns of simulation, formation of the main-chain hydrogen bond within the 310-helical region was completed. In short, the chronological order of hydrogen bond formation showed that formation of the 310-helix initiated after the formation of ahelical segment and the formation of all helices
714
was completed around 29 ns, which was roughly the folding time of the entire protein (judged by the RMSD value with respect to native structure). The ability to form both a- and 310-helices is noteworthy given the on-going debate on the relative significance of a- versus 310-helices,11 indicating a reasonable balance between these two conformations in the force field. Close inspection of the folding events indicates that formation of the 310helix was preceded by a short a-helix. Thus, contrary to the earlier suggestions,12 the 310-helix is not a transition state of a-helix. Perhaps one of the most interesting aspects of molecular dynamics simulation is that one can track the time evolution of structural changes at high temporal resolutions, a feat difficult to do with current experimental methods. In our analysis, particular attention was paid to the changes in the native contacts made by residue Trp25 as it constitutes the landmark feature of the mini-protein tc5b where its side-chain is buried by three nearby proline residues and one tyrosine. Visualization of the folding movie showed that the Trp side-chain interacted with Pro36, 37 and 38 in the early stage of simulation (within 1 ns of simulation) which stabilized the native topology of the protein. The hydrophobic contacts between Trp25 and Leu21 and Ile23 were also noted in the early stage of simulation. Thus, the first stage of folding consisted of the formation of native-like backbone topology, stabilized by the clustering of hydrophobic residues around residue Trp25. However, contacts between Trp25 and Pro31 as well as between Asp28 and Arg35 (potential salt bridge) were not observed in the early stage of simulation. The second milestone in the folding process was the formation of secondary structures. During this stage, the a-helical contact between Trp25 and Tyr22 began to build up after 1 ns of simulation, signifying that the a-helix was initiated from the N terminus. This point was further supported by the early contacts made between Gln24 and Leu21. In contrast, the contacts between Gln24 and Lys27 were established at around 15 ns, completing the formation of the C terminus end of the a-helix. The stable salt bridge between Asp28 and Arg35 was also observed after 15 ns of simulation. The fact that salt bridge formation occurred at this later stage of folding alludes to our initial hypothesis that the hydrophobic force was the main driving force for the native protein topology, whereas the salt bridges provided additional stability to the native structure. Different (structural) states in the folding pathway were examined by clustering analysis. This analysis would inform us about the different meta-stable states within the initial straight chain conformation and final native state. It would also tell us about the relative population of each detectable state which is a measure of the relative stability of each state. A total of 19 clusters were identified by the hierarchical clustering approach, however, only nine clusters were populated with
Ab initio Folding of the Trp-cage Mini-protein
more than 50 snapshots. Visualization of these nine clusters further subdivided the clusters into four intermediate states and one native state, which correspond closely to the development in main-chain RMSD. The first intermediate state represents the initial transient state which was formed by the initial hydrophobic collapse of the straight chain structure. It was stable for no more than 1 ns. The second intermediate state represents structures that exhibited partial helical turn with Leu2-carbonyl to Trp6-NH and Tyr3-carbonyl to Leu7-NH hydrogen bonds, another evidence that helix formation was initiated from the N terminus. About 13% of all saved structures fell into this state. Although the overall topology of this intermediate state differed significantly from the native structure, the hydrophobic contacts between Trp25 sidechain and Pro36-38 were present. The third intermediate state, having a population of 10% of the total saved snapshots, also contained structures with partial helix. However, the overall topology of these structures was much closer to the native structure. The fourth intermediate state, consisting of 4% of the saved snapshots, is characterized by the presence of the completed a-helical segment, although the small 310-helical segment was still absent. Finally, the native state structure, with ahelical and short 310-helical segment completed, is the most populated state with about 72% of all saved structures. Superimposition of five structures randomly taken from the native state cluster to that of the NMR ensemble is shown in Figure 3. One can clearly see from the stereographic Figure that the two ensembles of structures coincided very closely with each other. Most notably, the side-chain packing in this state is almost identical to the native NMR structure. Our simulation results indicate that the rate-limiting step is not the initiation or formation of ahelical segment. We observed that the helix could form its first turn rapidly within nanosecond time scale and completed within 15 ns. There is at least 20 ns time lag between the initiation of helix and completion of protein folding. In our opinion, disparity between these two important events rules out the possibility that the helix initiation played any significant role in the folding of tc5b. Our data suggest that final sharp drop of the RMSD values which brings tc5b structure into the vicinity of the native structure coincide with the sharp increase in native contacts or packing of the Trp25 side-chain. The time series of heavy atom RMSD plot (Figure 2(A)) clearly showed four stages of conformational changes as there are four plateaus in this time series. In contrast, there are only three plateaus in the atom to atom native contact time series plot of Trp25 (Figure 2(C)). The first, third and the fourth major transitions in heavy atom RMSD all correlated to a corresponding change in Trp25 native contacts. The final descent in RMSD value (at , 29 ns) correlated with a sharp increase in native contacts of Trp25 which
Ab initio Folding of the Trp-cage Mini-protein
715
Figure 3. Stereo view of the five representative NMR structures (red) superimposed on the five simulated snap shots (blue). Simulated snap shots are taken from the last 5 ns of simulation at 1 ns interval. Only heavy atoms are shown for clarity.
represents the last step of conformational rearrangement involving correct packing of Trp25 side-chain that led to the formation of native tc5b structure. This can be visualized from the folding movie. Further detailed contact analysis between side-chain atoms of Trp25 and individual atom of each other residue revealed that contacts between Trp25 side-chain and side-chains of Tyr22, Gly30, Pro31, Pro37 and Pro38 are mainly responsible for the final sharp drop in heavy atom RMSD. Interestingly these are the residues that essentially entrap the Trp25 side-chain. Thus, the last step of tc5b folding could be the correct packing of Trp side-chain. Although unequivocal assignment of the ratelimiting step requires more detailed thermodynamic analysis, based on our observation, we believe that the correct packing of Trp25 side-chain is most likely the rate-limiting step where upon packing, complete folding of the entire protein immediately followed. Indeed, in our subsequent simulations, packing of Trp side-chains was seen as the ratelimiting step, despite the quick formation of both secondary and tertiary structures (detail will be reported elsewhere). Nevertheless, one possible experimental validation could be measurement of helix formation in comparison to the packing of fluorescent Trp residues. Our simulation suggests that helix formation precedes completion of the folding. Identification of the rate-limiting step and the nature of transition state has been the focus in the field of protein folding mechanism studies. Presently, an important theory attributes the transition state to the chain entropy.13 On the other hand, Fvalue analyses14 have suggested a much more complex picture. We think our results agree with the latter better simply because formation of the native-like main-chain structure was early in the process and did not appear to be the rate-limiting step. In our opinion, this is not necessarily contrary to the conclusion of Plaxco et al., instead it is comp-
lementary. In the study of Plaxco et al. folding rates of multiple small proteins were examined. Given the diversity of the protein structures, a common denominator would the contribution of the chain entropy and the specific interactions that play roles in the folding of individual proteins would appear as the “noise” and uncertainties. Therefore, when one asks the common factor affecting folding of a diverse pool of proteins, chain-entropy would play an important role. However, in the folding of individual protein, specific interactions, such as those identified by F-value analysis and those described in this study, would be needed to describe the process. It should be noted that our present study focuses on the folding events observed in one simulation trajectory. Despite its high degree of realism, we must keep in mind that alternative folding processes are also likely to occur given the complexity of the free energy landscape. Therefore, we should focus on the qualitative aspects of the simulation results, for example, the sequence of events. The time scales of the events should be obtained by multiple simulations. This is particularly true since we employed a continuum solvent model in our simulation which is known to accelerate the process due to lack of solvent viscosity. Thus our observed time scales are likely to be much faster than what happen in reality. Nevertheless, by focusing on the qualitative aspects and the sequence of events, a good deal can be learned, though much of these remain to be tested by large number of simulations. Simmerling et al. recently studied the same protein with a similar approach15 in which a standard AMBER force field was used with a set of adjusted main-chain torsion parameters. In their simulation, the heavy atom RMSD also showed a step-wise reduction. The ability to fold the same protein with different force fields clearly demonstrates the robustness of all-atom molecular dynamics simulations. Interestingly, their observed folding time
716
was around 7 ns, much shorter than what we observed in our simulation. Now, we have gone a step further. Our simulation suggests that the ratelimiting step could be packing of the Trp25 sidechain and formation of the main-chain tertiary contacts could be an early event in the folding process. Qiu et al. recently studied folding of the protein16 by laser temperature jump experiments. They found that the extrapolated folding time is about 4 ms by monitoring the fluorescence of Trp25 sidechain which measures its packing. In our opinion, this is consistent with our observation that packing of Trp25 side-chain is the rate-limiting step. The difference between the simulated and experimental folding times is likely due to the absence of solvent viscosity in the simulation and the time was extrapolated from single simulation. Our observed formation time of helical secondary structure was on the order of 25 ns, close to experimental observations on isolated helices. Thus, the a-helix may form early. Our simulation also showed that tertiary contacts can form quickly. When tertiary contacts form, however, the Trp25 formed a hydrophobic cluster outside the cage and the cage was in a closed state, both of which prevented the Trp from entering the cage quickly. Thus, breaking the hydrophobic cluster, opening up the cage, and encaging Trp are all parts of the process that eventually leads to the packing of the Trp side-chain. These steps could be much slower than what we observed due to the lack of solvent viscosity in the simulation.
Concluding Remarks Present folding simulation of tc5b explores the folding pathway of 20-residue mini-protein. The main-chain as well as side-chain configuration of the simulated structure is in excellent agreement with the NMR ensembles. The formation of Trpcage was also closely monitored, which reveals that hydrophobic interaction between Trp and Pro residues initiates the folding process by forming the correct protein topology. Initial data indicate that main-chain conformation and secondary structure formation is not the rate-limiting step; rather, correct packing of Trp25 side-chain may be the rate-limiting step of the protein folding. Additional folding as well as unfolding simulations (work in progress) may give a more accurate account of the folding kinetics of this protein.
Method In our folding simulation with AMBER package, we began from a fully extended conformation of tc5b. Using the Generalized Born implicit solvent model9 to treat solvation effects and a new all-atom point-charge force field recently developed in our laboratory, we performed 100 ns of molecular dynamics simulation to fold the protein ab initio. The detail of the force field will be published elsewhere. In summary, the charges were
Ab initio Folding of the Trp-cage Mini-protein
obtained by fitting to the electrostatic potential calculated quantum mechanically in organic solvent ð1 ¼ 4Þ: The main-chain torsion parameters were obtained by fitting to the energy surface of alanine dipeptide. The energy surface was also calculated quantum mechanically in organic solvent ð1 ¼ 4Þ: Other parameters, including bonds, bond angles, torsion parameters (except the main-chain torsion angles), and van der Waals parameters were identical to those described by Wang et al.10 The initial extended chain conformation was constructed from the primary sequence using leap module in AMBER package. After 500 steps of energy minimization, random velocities were assigned at 100 K. The simulation was allowed to continue for 100 ns with an integration time step of 2.0 fs. The temperature was controlled at 300 K by loosely coupling the system to a heat bath with a coupling constant of 2.0 ps. Non-bonded interactions were calculated without cutoff (in both MD and minimization). The trajectory was saved at 10 ps intervals for further analysis. The simulation was done on a dual processor (Pentium III 1 GHz) PC work station which gave a throughput of about 5 ns a day. The entire simulation took roughly 20 days to complete.
Acknowledgements Computer time was provided by Pittsburgh Supercomputer Center. This work was supported by research grants from NIH (RR15588, PI Lenhoff and GM64458 to Y.D.), the state of Delaware, and University of Delaware Research Fund. Usage of VMD, UCSF Midas, WebLab viewer graphics packages are gratefully acknowledged.
References 1. Dahiyat, B. I. & Mayo, S. L. (1997). De novo protein design: fully automated sequence selection. Science, 278, 82 – 87. 2. Hill, R. B. & DeGrado, W. F. (1998). Solutions structure of a2D, a nativelike de novo designed protein. J. Am. Chem. Soc. 120, 1138– 1145. 3. Struthers, M. D., Cheng, R. P. & Imperiali, B. (1996). Design of a monomeric 23-residue polypeptide with defined tertiary structure. Science, 271, 342– 345. 4. Kortemme, T., R, A. M. & Serrano, L. (1998). Design of a 20-amino acid, three-stranded b-sheet protein. Science, 281, 253– 256. 5. Neidigh, J. W., Fesinmeyer, R. M. & Andersen, N. H. (2002). Designing a 20-residue protein. Nature Struct. Biol. 9, 425–430. 6. Li, A. & Daggett, V. (1996). Identification and characterization of the unfolding transition state of chymotrypsin inhibitor 2 by molecular dynamics simulations. J. Mol. Biol. 257, 412–429. 7. Duan, Y. & Kollman, P. A. (1998). Pathways to a protein folding intermediate observed in a 1 ms simulation in aqueous solution. Science, 282, 740– 744. 8. Snow, C. D., Nguyen, H., Pande, V. S. & Gruebele, M. (2002). Absolute comparison of simulated and experimental protein-folding dynamics computationally. Nature, in press.
Ab initio Folding of the Trp-cage Mini-protein
9. Tsui, V. & Case, D. A. (2000). Molecular dynamics simulations of nucleic acids with a generalized born solvation model. J. Am. Chem. Soc. 122, 2489– 2498. 10. Wang, J. M., Cieplak, P. & Kollman, P. A. (2000). How well does a restrained electrostatic potential (RESP) model perform in calculating conformational energies of organic and biological molecules? J. Comp. Chem. 21, 1049– 1074. 11. Tran, T. T., Zeng, J., Treutlein, H. & Burgess, A. W. (2002). Effects of thioamide substitutions on the conformation and stability of a- and 310-helices. J. Am. Chem. Soc. 124, 5222– 5230. 12. Basu, G., Kitao, A., Hirata, F. & Go, N. (1994). A collective motion description of the 310-/a-helix transition—implications for a natural reaction coordinate. J. Am. Chem. Soc. 116, 6307– 6315.
717
13. Plaxco, K. W., Simons, K. T. & Baker, D. (1998). Contact order, transition state placement and the refolding rates of single domain proteins. J. Mol. Biol. 277, 985 –994. 14. Fersht, A. R. (1995). Mapping the structures of transition-states and intermediates in folding—delineation of pathways at high-resolution. Phil. Trans. Roy. Soc. London B-Biol. Sci. 348, 11 –15. 15. Simmerling, C., Strockbine, B. & Roitberg, A. E. (2002). All-atom structure prediction and folding simulations of a stable protein. J. Am. Chem. Soc. 124, 11258 –11259. 16. Qiu, L., Pabit, S. A., Roitberg, A. E. & Hagen, S. J. (2002). Smaller and faster: the 20-residue Trp-Cage protein folds in 4 ms. J. Am. Chem. Soc. 124, 12952 – 12953.
Edited by M. Levitt (Received 10 September 2002; received in revised form 26 November 2002; accepted 23 January 2003)