Vulnerability in Popular Molecular Dynamics Packages Concerning Langevin and Andersen Dynamics.

We report a serious problem associated with a number of current implementations of Andersen and Langevin dynamics algorithms. When long simulations are run in many segments, it is sometimes possible to have a repeating sequence of pseudorandom numbers enter the calcuation. We show that, if the sequence repeats rapidly, the resulting artifacts can quickly denature biomolecules and are then easily detectable. However, if the sequence repeats less frequently, the artifacts become subtle and easily overlooked. We derive a formula for the underlying cause of artifacts in the case of the Langevin thermostat, and find it vanishes slowly as the inverse square root of the number of time steps per simulation segment. Numerous examples of simulation artifacts are presented, including dissociation of a tetrameric protein after 110 ns of dynamics, reductions in atomic fluctuations for a small protein in implicit solvent, altered thermodynamic properties of a box of water molecules, and changes in the transition free energies between dihedral angle conformations. Finally, in the case of strong thermocoupling, we link the observed artifacts to previous work in nonlinear dynamics and show that it is possible to drive a 20-residue, implicitly solvated protein into periodic trajectories if the thermostat is not used properly. Our findings should help other investigators re-evaluate simulations that may have been corrupted and obtain more accurate results.

[1]  William Feller,et al.  The fundamental limit theorems in probability , 1945 .

[2]  G. Ciccotti,et al.  Numerical Integration of the Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of n-Alkanes , 1977 .

[3]  H. C. Andersen.  Molecular dynamics simulations at constant pressure and/or temperature , 1980 .

[4]  J. Mccammon,et al.  Molecular dynamics with stochastic boundary conditions , 1982 .

[5]  M. Karplus,et al.  CHARMM: A program for macromolecular energy, minimization, and dynamics calculations , 1983 .

[6]  M. Karplus,et al.  Stochastic boundary conditions for molecular dynamics simulations of ST2 water , 1984 .

[7]  H. Berendsen,et al.  Molecular dynamics with coupling to an external bath , 1984 .

[8]  Hoover,et al.  Canonical dynamics: Equilibrium phase-space distributions. , 1985, Physical review. A, General physics.

[9]  T. Straatsma,et al.  THE MISSING TERM IN EFFECTIVE PAIR POTENTIALS , 1987 .

[10]  C R Cantor,et al.  Cooperative biotin binding by streptavidin. Electrophoretic behavior and subunit association of streptavidin in the presence of 6 M urea. , 1990, The Journal of biological chemistry.

[11]  M. Wilchek,et al.  The quaternary structure of streptavidin in urea. , 1991, The Journal of biological chemistry.

[12]  P. Kollman,et al.  Settle: An analytical version of the SHAKE and RATTLE algorithm for rigid water models , 1992 .

[13]  Hamann,et al.  Transition from chaotic to nonchaotic behavior in randomly driven systems. , 1992, Physical review letters.

[14]  Alan M. Ferrenberg,et al.  Monte Carlo simulations: Hidden errors from "good" random number generators. , 1992, Physical review letters.

[15]  Whitlock,et al.  Pseudorandom number generator for massively parallel molecular-dynamics simulations. , 1994, Physical review. E, Statistical physics, plasmas, fluids, and related interdisciplinary topics.

[16]  D. van der Spoel,et al.  GROMACS: A message-passing parallel molecular dynamics implementation , 1995 .

[17]  T. Darden,et al.  A smooth particle mesh Ewald method , 1995 .

[18]  P. Kollman,et al.  A Second Generation Force Field for the Simulation of Proteins, Nucleic Acids, and Organic Molecules , 1995 .

[19]  C. Cantor,et al.  Streptavidins with intersubunit crosslinks have enhanced stability , 1996, Nature Biotechnology.

[20]  O. Tapia,et al.  On the sensitivity of MD trajectories to changes in water‐protein interaction parameters: The potato carboxypeptidase inhibitor in water as a test case for the GROMOS force field , 1996, Proteins.

[21]  D. Peter Tieleman,et al.  Molecular dynamics simulations of peptides from BPTI: A closer look at amide—aromatic interactions , 1996, Journal of biomolecular NMR.

[22]  B. Katz.  Binding of biotin to streptavidin stabilizes intersubunit salt bridges between Asp61 and His87 at low pH. , 1997, Journal of molecular biology.

[23]  W. L. Jorgensen,et al.  Temperature dependence of TIP3P, SPC, and TIP4P water from NPT Monte Carlo simulations: Seeking temperatures of maximum density , 1998 .

[24]  W. C. Still,et al.  Approximate atomic surfaces from linear combinations of pairwise overlaps (LCPO) , 1999 .

[25]  I Vattulainen,et al.  Framework for testing random numbers in parallel calculations. , 1999, Physical review. E, Statistical physics, plasmas, fluids, and related interdisciplinary topics.

[26]  P. Kollman,et al.  How well does a restrained electrostatic potential (RESP) model perform in calculating conformational energies of organic and biological molecules? , 2000 .

[27]  T. N. Bhat,et al.  The Protein Data Bank , 2000, Nucleic Acids Res..

[28]  D. Case,et al.  Modification of the Generalized Born Model Suitable for Macromolecules , 2000 .

[29]  Synchronization induced by Langevin dynamics. , 2001, Physical review. E, Statistical, nonlinear, and soft matter physics.

[30]  R. Skeel,et al.  Langevin stabilization of molecular dynamics , 2001 .

[31]  Berk Hess,et al.  GROMACS 3.0: a package for molecular simulation and trajectory analysis , 2001 .

[32]  A. Roitberg,et al.  All-atom structure prediction and folding simulations of a stable protein. , 2002, Journal of the American Chemical Society.

[33]  P. M. Rodger,et al.  DL_POLY: Application to molecular simulation , 2002 .

[34]  J. W. Neidigh,et al.  Designing a 20-residue protein , 2002, Nature Structural Biology.

[35]  M. Wilchek,et al.  Dimer-Tetramer Transition between Solution and Crystalline States of Streptavidin and Avidin Mutants , 2003, Journal of Bacteriology.

[36]  B. Brooks,et al.  Self-guided Langevin dynamics simulation method , 2003 .

[37]  Laxmikant V. Kalé,et al.  Scalable molecular dynamics with NAMD , 2005, J. Comput. Chem..

[38]  Barry Honig,et al.  Comparative study of generalized Born models: protein dynamics. , 2005, Proceedings of the National Academy of Sciences of the United States of America.

[39]  Holger Gohlke,et al.  The Amber biomolecular simulation programs , 2005, J. Comput. Chem..

[40]  E. A. Koopman,et al.  Advantages of a Lowe-Andersen thermostat in molecular dynamics simulations. , 2006, The Journal of chemical physics.

[41]  Y. Duan,et al.  Two-stage folding of HP-35 from ab initio simulations. , 2007, Journal of molecular biology.