Atomic-detailed milestones along the folding trajectory of protein G

The high computational cost of carrying out molecular dynamics simulations of even small–size proteins is a major obstacle in the study, at atomic detail and in explicit solvent, of the physical mechanism which is at the basis of the folding of proteins. Making use of a biasing algorithm, based on the principle of the ratchet–and–pawl, we have been able to calculate eight folding trajectories (to an RMSD between 1.2Å and 2.5Å) of the B1 domain of protein G in explicit solvent without the need of high–performance computing. The simulations show that in the denatured state there is a complex network of cause-effect relationships among contacts, which results in a rather hierarchical folding mechanism. The network displays few local and nonlocal native contacts which are cause of most of the others, in agreement with the NOE signals obtained in mildly-denatured conditions. Also nonnative contacts play an active role in the folding kinetics. The set of conformations corresponding to the transition state display φ–values with a correlation coefficient of 0.69 with the experimental ones. They are structurally quite homogeneous and topologically nativelike, although some of the side chains and most of the hydrogen bonds are not in place. 1 ar X iv :0 90 5. 28 75 v1 [ qbi o. B M ] 1 8 M ay 2 00 9 Molecular–dynamics simulations in explicit solvent can be a very useful complement to experimental studies of protein folding, in keeping with the fact that they provide insight into the time evolution of the process with atomic detail, under fully controlled conditions [1]. On the other hand, they are computationally very demanding, even in the case of small proteins. Among the most massive folding simulations ever realized is a 10μs molecular dynamics (MD) folding trajectory of the 38–residue WW domain, lasting for about 3 months on 329 cores and reaching conformations which are ∼ 50% similar to the native conformations in terms of number of contacts [2]. To be statistically sound, Pande and coworkers carried out 410 simulations of the folding of the 35–residue Villin Headpiece, the average duration being 863 ns. The calculation lasted for 54 machine years on a distributed computer, and eighteen of these trajectories reached the native conformation [3]. The intrinsic and unavoidable computational problem in carrying out folding simulations with realistic protein models is the wide range of time scales involved: the time step of the simulation must be tuned to femtoseconds, corresponding to the time scale of atomic vibrations, while the overall folding process spans over interval of time ranging from milliseconds to seconds. In an attempt to overcome this difficulty, a number of investigations focused on the study of unfolding simulations at high temperature [4, 5]. A decade ago Marchi and Ballone developed an adiabatic bias molecular dynamics (ABMD) method [6], to generate MD trajectories between pairs of points in the conformational space of complex systems. It was applied for the first time to protein unfolding by Paci and Karplus [7]. The method is based on the introduction of a biasing potential which is zero when the system is moving towards the desired arrival point and which damps the fluctuations when the system attempts at moving in the opposite direction. As in the case of the ratchet and pawl system, propelled by thermal motion of the solvent molecules, the biasing potential does not exert work on the system. Consequently, the resulting trajectories are physically correct. On the other hand, the algorithm cannot provide the statistical weight of the visited states nor the time scales associated with the trajectory. In the present work we report on results of the application of the ABMD algorithm to the study, with the help of the Amber force field [8] in explicit solvent and without recurring to high–performance computing, of the folding of the 56–residues B1 domain of protein G starting from 16 thermally–unfolded conformations. From the eight trajectories which reached a RMSD lower than 2.5Å we have extracted the conditional probabilities of contact