Mu, Y.

Yang, Y.

Juelich Shared Electronic Resources

John von Neumann Institute for ComputingModeling of Water Explicitly in theReplica-Exchange Simulation Method forProtein FoldingMu Yuguang, Yang Yepublished inNIC Workshop 2006,From Computational Biophysics to Systems Biology,Jan Meinke, Olav Zimmermann,Sandipan Mohanty, Ulrich H.E. Hansmann (Editors)John von Neumann Institute for Computing, Ju¨lich,NIC Series, Vol. 34, ISBN-10: 3-9810843-0-6,ISBN-13: 978-3-9810843-0-6, pp. 119-224, 2006.c© 2006 by John von Neumann Institute for ComputingPermission to make digital or hard copies of portions of this work forpersonal or classroom use is granted provided that the copies are notmade or distributed for profit or commercial advantage and that copiesbear this notice and the full citation on the first page. To copy otherwiserequires prior specific permission by the publisher mentioned above.http://www.fz-juelich.de/nic-series/volume34Modeling of Water Explicitly in the Replica-ExchangeSimulation Method for Protein FoldingMu Yuguang∗ and Yang YeSchool of Biological Sciences,Nanyang Technological University, Singapore 637551∗Correspondent Author: E-mail: ygmu@ntu.edu.sgConsidering the replica-exchange simulation of protein in explicit water to be a two-Hamiltonian system with one Hamiltonian for conformational sampling in molecular dynamicssimulation and another for controlling the exchanges between difference temperature replicasin Monte Carlo jumps, we introduce an approximation on the latter Hamiltonian using a con-tinuum solvent model, surface-generalized Born model. Such replica exchange simulation withhybrid Hamiltonians method is applied to fold the C-terminus (residue 41-56) of protein G fromextended structures. Promising results show that not only the total number of replica needed islargely reduced but also the folding efficiency is greatly enhanced. Combined with recently in-vented dihedral principle component analysis a general framework for ab initio folding a smallprotein merely from sequence knowledge is emerging.PACS numbers: 87.15.Cc, 87.15.AaPredicting protein structure solely from its amino acid sequence has long been a greatchallenge in modern molecular biology1, 2. Although remarkable progress was maderecently by Baker and his collaborators3, 4, the problem remains unsolved. Recent studieson peptides5, 6 by molecular dynamics simulation tools using explicit water model, aidedby advanced sampling strategies, such as replica exchange molecular dynamics (REMD)method7, show a prospective way towards solving the folding puzzle. However whenREMD is applied to larger protein system a huge computation facility will be neededbecause the number of replicas needed increases simultaneously with the increase of thedegrees of freedom of the system8. Endeavors to overcome the shortage of REMD havebeen undertook by several groups8–12.The Hamiltonian of a protein solvated in an aqueous environment can be written asHsys = Hp + Hpw + Hw, where p, pw and w denote protein-protein, protein-waterand water-water interactions respectively. In standard REMD scheme, replicas arepropagating independently at certain temperatures with Monte Carlo (MC) exchange atcertain intervals. The MC exchange corresponds to the Metropolis criterion following thedetailed balance condition: pi↔j = min(1, e∆β∆E), where ∆β = βj − βi, β = 1kBT ,∆E = E (χj − χi) and χi, χi are the configurations of the neighboring replicas. The totalenergy,E (χ), should be the same as the system Hamiltonian and can also be decomposedinto three terms: E (χ) = Ep + Epw + Ew as well. In explicit water REMD simulationmany water molecules are employed to solvate well the unfolded protein. This makes thetotal energy, E, and its differenceds between neighbours, ∆E, huge. In order to get areasonable high exchange probability between the neighboring temperature ladders, pi↔j ,the temperature jumps, ∆β, should be small. Therefore a large number of replicas areneeded to cover a broad temperature range, which from 300K to 600K . Realizing that119the water-water interaction makes the largest contribution to the total energy a possibleworkaround could be that decouple the water-water interaction from the thermal bathsand make the ∆E small though the total energy E nearly unchanged. Recent works doneby Berne group9 and Simmerling group11 were pursued in this direction using differentstrategies. However such decoupling not only causes unphysical Hamiltonian on mostreplicas in MD simulation, the resulting hot protein and cold water effect would alsohinder the conformation sampling efficiency at high temperatures.This letter presents an improved solution to this problem. The standard REMD methodcan be thought as a hybrid method of MD and MC simulation. The Hamiltonians forboth MD and MC should be the same for the sake of consistency. However if one checksthe roles of both Hamiltonians in details, such consistency is found to be not a necessity.The full system Hamiltonian, Hsys, is necessary for MD simulations to correctly sampleconformational space. While for MC steps we find it is not necessary to use the sameHamiltonian, Hsys. In Hsys the water-water interaction is dominant but is not directlyrelated to the protein folding. The key of the MC steps in REMD is to help protein jumpout of its local minima. Evidently the REMD could be more efficient in folding proteinif the Hamiltonian for MC steps represents the energies of the protein more directly. Inlight of that we introduce an approximation to the Hamiltonian for MC steps. Insteadof using the Hsysas in the MD steps a second Hamiltonian with the same energies forprotein-protein interaction but protein-water and water-water interactions approximatedby a continuum solvent model, surface-generalized Born model (GBSA)13 is used in MCsteps. In this way the dominate water-water interaction is averaged out and therefore asmall ∆E is obtained. At current moment the implicit water model or continuum solventmodel, although very computationally efficient, was found not to be able to fold peptidecorrectly14. Our method is a complement for such deficiency in the GBSA model.We call this new version of REMD as REMD with hybrid Hamiltonians (REMDhH).Hybrid Hamiltonians costs the violation of detailed balance of the whole system to somedegree. Considering the interested temperature range for REMD is from 300K to 600Kthe water behaves well as a liquid (in constant volume REMD) and should be quicklyrelaxed to its equilibrium state when the system jumps to a new temperature.Here REMDhH is applied to fold the C terminus (residue 41-56) of protein G, a 16amino acids peptide which was found to be able to fold into a beta-hairpin in vitro15.The GROMACS program suite16 and the full atomic OPLS-AA force field17 are used.The peptide is capped with the normal ACE and NME groups with 256 atoms in total.It is solvated by 5469 water molecules plus three K+ ions to neutralize the chargedmolecular system. The whole system consists of 16666 atoms in a cubic simulation boxof 5.3 nm length. Sixteen replicas are used whose temperatures are 300.0, 317.2, 335.4,354.6, 374.9, 396.4, 419.1, 443.1, 468.5, 495.4, 523.8, 553.8, 585.5, 619.1, 654.5 and692.1 Kelvin. While 64 replicas had been used for the same peptide with less watermolecules (1361) studied by standard REMD5. The GBSA energies are calculated usingTinker program18. The interface between GROMACS and Tinker is built by modifyingGROMACS source code. A twin-range cutoff of 0.9/1.4 nm is used for the non-bondedinteractions and a reaction-field correction with permittivity is employed. The integration1200 20 40 60 80 100time [ns]0481204812RMSD [A] 048120 20 40 60 80 100time [ns]200400600800400600800temprature[K]400600800(a)(b)(c)(d)(e)(f)Figure 1. (a), (b) and (c) are the three folded replica trajectories of RMSD from NMR structure, (d), (e) and (f)are the trajectories of temperature jumps of the related replicas, respectively.step in all simulations is 0.002 ps. Non-bonded pair lists are updated every 10 integrationsteps. The system is coupled to an external heat bath with a relaxation time of 0.1 ps.All bonds including hydrogen atoms are constrained in length. Each replica is run for 98ns with replica exchange attempted every 2 ps. The final acceptance ration for replicaexchanges is found to be 35%-43%.To test the validity of the method the initial configurations of peptide are extendedstructures, which means only the information of amino acids sequence is supplied. Threefolding events are found during the simulation. Figure 1a, 1b and 1c show the trajectoriesof root mean squared deviations (RMSD) from the NMR structure19 calculated withall atoms. The fastest folding happens in replica 9 (Figure 1b) around simulation timeof 18 ns. The following folding events happen in replica 4 and 15 (Figure 1a and 1c)around simulation time of 60 ns. All the folded structures maintain stable until the endof simulations. Here we denote a structure as folded when its RMSD drops below 4A.A more detailed analysis is given below. The folding processes revealed by REMD are aprocess of annealing and relaxation. The temperature trajectories (Figure 1d, 1e and 1f) ofthe three folded replicas indicate this view clearly. Accompanying with the decrease ofRMSDs the temperatures are cooling down.In order to find the quantitative reasons why the employed REMDhH method issuperior to the standard REMD, the correlation between RMSDs and the total energyof the system with explicit water molecules, and the correlation between RMSDs andthe GBSA energies are calculated and displayed in Figure 2a and 2b respectively. It isclearly shown that the correlation between the RMSDs and the total energies used byregular REMD is poor. Its correlation coefficient is 0.01. While the correlation between1210 2 4 6 8 10 12 14RMSD[A]−3500−3400−3300−3200−3100−3000E gbsa[kJ/mol]−2320−2310−2300−2290−2280E tot[100kJ/mol](a)(b)Figure 2. (a) Correlations of the total energies, including protein-protein, protein-water, and water-water inter-actions, with RMSD from NMR structure. The correlation coefficient is 0.01. (b) Correlations of the GBSAenergies with RMSD. The correlation coefficient is 0.3. The white lines are the linear regression results.the RMSDs and the energies from the GBSA model is much stronger whose correlationcoefficient is 0.3. As mentioned earlier the folding process in REMD method couldbe regarged as an annealing process. To perform such annealing efficiently two factorsfor the MC steps must be considered. First the lower energy conformations could beidentified in the high temperature conformational ensemble. Secondly these lower energyconformations have to be quickly annealed to the lower temperature simulations. The firstfactor requires that the energies for MC steps should interrogate the protein conformationdirectly. Otherwise the lower energy conformations would not easily be identified.The second factor claims that the fewer temperature ladders the quicker the annealingprocess which implies that the number of replicas in REMD should be minimized. Theimplementation of the current REMDhH could make both requirements satisfactory. If theREMD efficiency is approximated to be linearly correlated with the correlation coefficientof RSMD/energy for MC steps the speedup of folding by REMDhH compared with thatof standard REMD is a factor of 30 for this system.In protein structure prediction usually the folded structure is unknown. The RMSDsfrom the native structure and the fraction of native contacts which are commonly usedas folding reaction coordinates are not available in such situations. In order to identifythe global minimum new suitable reaction coordinates are desirable. We find that therecently invented dihedral principal component analysis (dPCA) method works20 finehere. Figure 3 reveals the folding free energy curve for this beta-hairpin peptide obtainedby projecting conformational ensemble of T = 300K onto the first eigenvector of dPCA122−4 −2P10246810Free Energy [kJ/mol]1 3 5 7 90123eigenvalue0 5 10 1500.10.2probabilityrmsdlabel(a) (b)0 2Figure 3. Folding free energy curve of the GB1 hairpin peptide obtained by projection onto the eigenvector oflargest eigenvalue using dihedral PCA method. Inset (a) shows the first 10 largest eigenvalues of dPCA result.Inset (b) shows the distribution of RMSD of ensemble at T = 300K (solid line) and the distribution of RMSDof the ensemble within the global minimum using such criteria P1 > 0.7 (dashed line).whose eigenvalue is shown to be much larger than other eigenvalues in inset (a). The freeenergy curve is funnel-like and it is easy to identify the global minimum which is locatedaround P1 = 1.2. Here P1 is the projection value on the first eigenvector. To checkwhether this global minimum is related to the native state the distribution of RMSDsfrom this ensemble is plotted in Figure 3 inset (b) by the dashed line. For comparisonthe distribution of RMSDs for all the structures obtained at 300K is shown by the solidline. Evidently the structures from the global minimum are native-like with averageRMSD = 3.5A. And P1 is a better reaction coordinate than RMSD in the sense that theformer can distinguish folded states from unfolded ones more clearly.In summary a new version of REMD method is suggested. By decomposing REMDinto two steps, one is MD step and the other is MC step, different Hamiltonians can beapplied to them. In MD step the full Hamiltonian with all protein-protein, water-water,water-protein interactions included guarantees that the conformations of protein aresampled correctly. In MC step, however, only the protein-protein interaction is explicitlyconsidered, the water-water and water-protein interactions are approximated by a contin-uum GBSA model. Such a hybrid Hamiltonian scheme reduces the number of replicagreatly and meanwhile increases the folding efficiency by more than a factor of 10. Onthe other hand, dPCA analysis makes the global minimum easily identified if protein isfolded in simulation. Combining REMDhH with dPCA a genenal framework for ab initioprotein structure prediction is emerging. Together with the available of accurate forcefields and powerful computational facility the era of solving protein folding by brute forceMD simulations will be coming soon21.123AcknowledgmentsThe support of a Lee Kuan Yew Research Fellowship to Y.G.M. is acknowledged. Thesupport of PhD study in Nanyang Technological Univerisity to Yang Ye is appreciated. Thesimulations were performed on the supercomputer of Bioinformation Research Center andthe Frankfurt Center of Scientific Computing, which are acknowledged by the generousallocation of CPU time.References1. D. J. Wales and H. A. Scheraga, Science 285, 1368 (1999).2. J. N. Onuchic and P. G. Wolynes, Curr. Opin. Struct. Biol. 14, 70 (2004).3. K. T. Simons, C. Kooperberg, E. Huang, and D. Baker, J. Mol. Biol. 268, 209 (1997).4. C. A.Rohl, C. E.Strauss, K. M.Misura, and D. Baker, Methods Enzymol 383, 66(2004).5. R. Zhou, B. J. Berne, and R. Germain, Proc. Natl. Acad. Sci. USA 98, 14931 (2001).6. A. E. Garcia and J. N. Onuchic, Proc. Natl. Acad. Sci. USA 100, 13898 (2003).7. Y. Sugita and Y. Okamoto, Chem. Phys. Lett. 314, 141 (1999).8. H. Fukunishi, O. Watanabe, and S. Takada, J. Chem. Phys. 116, 9058 (2002).9. P. Liu, B. Kim, R. A. Friesner, and B. J. Berne, Proc. Natl. Acad. Sci. USA 102,13749 (2005).10. T. Z. Lwin and R. Luo, J. Chem. Phys. 123, 194904 (2005).11. X. Cheng, G. Cui, V. Hornak, and C. Simmerling, J. Phys. Chem. B 109, 8220 (2005).12. E. Lyman, F. M. Ytreberg, and D. M. Zuckerman, Phys. Rev. Lett. 96, 028105 (2006).13. Q. Liu, P. S. Shenkin, F. P. Hollinger, and W. C. Still, J. Phys. Chem. A 101, 3005(1997).14. R. Zhou and B. J. Berne, Proc. Natl. Acad. Sci. USA 99, 12777 (2002).15. F. J. Blanco, G. Rivas, and L. Serrano, Nat. Struct. Biol. 1, 584 (1994).16. H. J. C. Berendsen, D. van der Spoel, and R. van Drunen, Comp. Phys. Comm. 91,43 (1995).17. W. L. Jorgensen, D. S. Maxwell, and J. Tirado-Rives, J. Am. Chem. Soc. 118, 11225(1996).18. R. V. Pappu, R. K. Hart, and J. W. Ponder, J. Phys. Chem. B 102, 9725 (1998).19. A. M. Gronenborn, D. R. Filpula, N. Z. Essig, A. Achari, M. Whitlow, P. T. Wingfield,and G. M. Clore, Science 253, 657 (1991).20. Y. G. Mu, P. H. Nguyen, and G. Stock, Proteins 58, 45 (2005).21. M. Karplus and J. A. McCammon, Nat. Struct. Biol. 9, 646 (2002).124

Modeling of Water Explicitly in the Replica-Exchange Simulation Method for Protein Folding

http://juser.fz-juelich.de/record/152786/files/FZJ-2014-02422.pdf

Modeling of Water Explicitly in the Replica-Exchange Simulation Method for Protein Folding

Abstract

Similar works

Full text

Available Versions

Juelich Shared Electronic Resources