Classical molecular dynamics solves Newton's equations of motion for a system composed of many particles/atoms that interact through a given potential energy. The equations of motion are integrated using time-reversible algorithms that employ discrete time-steps. Based on the atoms positions and momenta at time "t", new positions at time "t+dt" are computed according to Newton's law and the corresponding force field (see [1]; pdf ). The use of a discrete time step introduces an intrinsic numerical error that can be monitored through the total energy of the system. This quantity should be conserved in a classical system described by Newton's equations. However, a molecular dynamics calculation may not conserve energy due to intrinsic numerical errors. In fact, variations in the total energy of the system can be used to test the accuracy of the numerical algorithm implemented and the time-step dt chosen [7,8]. Algorithms that employ multiple-time-step schemes are known to introduce a further uniform energy drift [5]. The numerical drift in energy, in general, does not distort the dynamic behavior of simulated systems, except for increasing temperature.

Figure 4 Figure 1 End-to-end distance for four constant-force SMD simulations performed on the solvated and equilibrated (sim5a) model of all 24 repeats of human ankyrin-R [2]. The N terminus of the protein was fixed while the C terminus was stretched with forces of 25 pN (blue, sim5b), 50 pN (red, sim5c) and 100 pN (black, sim5d) [2]. The end-to-end distance of the first control simulation is shown in green (50 pN, sim-c1; thicker for visualization purposes only).

Here we consider various numerical schemes for the simulation of the stretching on ankyrin reported in [2] (see our website on ankyrin here). An increase in the temperature of simulations sim5b and sim5c (up to 28K) in [2], caused by an energy drift related to the integration time steps chosen [3,4], raised concerns about the effect of the multiple-time step algorithm utilized on the elastic response of ankyrin (Joe Howard, private communication). The program NAMD and the algorithms it implements have been widely tested before (see for instance [4-6] and multiple simulations by independent groups). Nevertheless, we explicitly verified the elastic response of ankyrin in control simulations performed using more conservative sets of time steps. As expected, the control simulations yielded the same ankyrin elastic properties as reported previously with a choice of longer integration time steps [2], while the increase in temperature was drastically reduced (less than one degree Celsius during six nanoseconds) and energy was better conserved (less than 0.5% drift over six nanoseconds). The work done on ankyrin during stretching agreed well with the expected theoretical value. Naturally, the numerically more precise simulations decreased significantly the wall clock performance of the code; this was the reason why longer integration time steps had been chosen in [2].

The new simulations presented in detail below, as expected, justify the previously [2] used integration scheme.

Simulation Setup

Ankyrin
Figure 2 Front and side views of the initially equilibrated 24-repeat structure (sim5a) and the same structure after 3 ns and 6 ns of a constant force (50 pN) simulation using time steps of 1 fs, 1 fs, and 2 fs (sim-c1, see Simulation Setup).

Control simulations presented here were performed as in [2] but using a different set of integration time steps as described below. We used the program NAMD 2.5 and the CHARMM27 force field for proteins along with the TIP3P model for water molecules. A cutoff of 12 Å (switching function starting at 10 Å) for van der Waals interactions was assumed. Periodic boundary conditions (281 x 90 x 128 Å3) were used and distances between periodic images of the protein were greater than 30 Å over the trajectories. The particle mesh Ewald method (PME) was used to compute long-range electrostatic forces without cutoff (density of grid points for PME was at least 1/Å3 in all cases).

In the first two control simulations (sim-c1, sim-c2) the integration time steps utilized were 1 fs for interactions involving covalent bonds, 1 fs for short-range nonbonded interactions, and 2 fs for long-range electrostatic forces (as opposed to 1 fs, 2 fs, and 4 fs used previously [2]). An external force of 50 pN was applied to the C terminus while the N terminus was held fixed in simulation sim-c1. No external force was applied during sim-c2.

Three additional control simulations (sim-c3, sim-c4, sim-c5) employed a uniform integration time step of 1 fs for all interactions. An external force of 50 pN was applied to the C terminus while the N terminus was held fixed in simulation sim-c3. After 2.5 ns of dynamics in sim-c3, the external force was decreased from 50 pN to 25 pN in simulation sim-c4. No external force was applied during sim-c5.

Results

Temperature Figure 3 End-to-end distance for five constant-force SMD simulations performed on the solvated and equilibrated (sim5a) model of all 24 repeats of human ankyrin-R [2]. The N terminus of the protein was fixed while the C terminus was stretched with forces of 25 pN (blue, sim5b), 50 pN (red, sim5c) and 100 pN (black, sim5d) [2]. The end-to-end distance for sim-c3 (50 pN) is shown in light green while the end-to-end distance for sim-c4 (25 pN) is shown in turquoise (thicker for visualization purposes only).

The stack of 24 ankyrin repeats during sim-c1 showed the same elastic behavior as that reported for simulation sim5b in [2]. The end-to-end distance increased from ~ 112 Å to values over ~ 210 Å after six nanoseconds (see Fig. 1). No significant changes in the secondary structure of individual repeats were observed, as reported in [2] (see Fig. 2). A similar behavior was observed for sim-c3 (Fig. 3). When the external force applied during sim-c3 was reduced from 50 pN to 25 pN in sim-c4, the end-to-end distance decreased corroborating the reversible elastic behavior of ankyrin reported in [2] (Fig. 3).

Energy Drift, Temperature Increase and Work

Temperature Figure 4 Temperature as a function of time during sim-c1 (50 pN) and sim-c2 (0 pN). A small increase (less than one degree Celsius) is observed.

As in every experimental setup, different artifacts and sometimes unavoidable sources of error need to be considered. In order to quantify the magnitude of the error introduced by the numerical integration performed by NAMD, the temperature and total energy of the system were monitored during our control simulations (see Figs. 4, 5 and 6). An increase in the temperature of the system of less than one degree Celsius was observed for sim-c1 and sim-c2. The variation of the total energy of the system during the simulations was found to be less than 0.5 % in both cases. Moreover, the difference between the total energy for sim-c1 (50 pN) and sim-c2 (0 pN) agrees well with the work done by the external force (~ 38 kcal/mol observed in Fig. 5, 50 pN x 58 Å ~ 2.9 x 10-19 J).

When a uniform time step of 1 fs is employed (sim-c3, sim-c4 and sim-c5) the uniform energy drift and respective increase in temperature previously observed are practically negligible and therefore variations in the total energy of the system are caused mainly by the work done through external forces. This is clearly illustrated in Fig. 6. When no external forces were present, the energy remained at a steady value. When an external force of 50 pN was applied the end-to-end distance of ankyrin increased from 112 to 182 Å (sim-c3) and the energy of the system increased by about 60 kcal/mol (Fig. 6). When the external force was reduced from 50 pN to 25 pN, the end-to-end distance of ankyrin decreased from 183 to 164 Å and the total energy observed decreased by about 10 kcal/mol.

Total Energy Figure 5 Total energy of the simulated system during sim-c1 (50 pN, black) and sim-c2 (0 pN, red). An energy drift of less than 0.05%/nanosecond is observed. The inset shows a zoom of the plot at the end of the sim-c2 trajectory. The difference between these quantities accounts for the work done by the external force on ankyrin.
Total Energy Figure 6 Total energy (average TOTAL3, see NAMD user's guide) of the simulated system during simulations sim-c3 (50 pN, black), sim-c4 (25 pN, blue) and sim-c5 (0 pN, red). The variations in energy are consistent with the work done by the external forces.

Configuration and log files

  • Click here (12 Kb) to download the configuration files utilized for simulation sim-c1 (see the NAMD tutorial and NAMD user's guide for details).
  • Click here (15 Mb) to download the log files generated during simulation sim-c1.
  • Click here (3.2 Mb) to download the configuration and log files generated during simulation sim-c2.
  • Click here (6.4 Mb) to download the configuration and log files generated during simulation sim-c3.
  • Click here (3.5 Mb) to download the configuration and log files generated during simulation sim-c4.
  • Click here (1.8 Mb) to download the configuration and log files generated during simulation sim-c5.
  • Trajectory files containing the position of all atoms during the simulation are available upon request (~ 50 Gigabytes).

References

[1] Scalable Molecular Dynamics with NAMD. James C. Phillips, Rosemary Braun, Wei Wang, James Gumbart, Emad Tajkhorshid, Elizabeth Villa, Robert Skeel, Laxmikant Kale, and Klaus Schulten. Journal of Computational Chemistry, 2005. In press. pdf
[2] In search of the hair-cell gating spring: Elastic properties of ankyrin and cadherin repeats. Marcos Sotomayor, David P. Corey, and Klaus Schulten. Structure, 13:669-682, 2005.
[3] Nonlinear Resonance Artifacts in Molecular Dynamics Simulations. T. Schlick, M. Mandziuk, R. D. Skeel and K. Srinivas. J. Comput. Phys., 140:1-29, 1998
[4] Verlet-I/r-RESPA Is Limited by Nonlinear Instability. Q. Ma, J. Izaguirre, and R. D. Skeel. SIAM J. Sci. Comput., 24(6), 1951-1973, 2003 (pdf)
[5] Difficulties with Multiple Timestepping and the Fast Multipole Algorithm in Molecular Dynamics. T. Bishop and Robert D. Skeel and K. Schulten. J. Comp. Chem., 18(14):1785-1791, 1997.
[6] Langevin Stabilization of Molecular Dynamics J. A. Izaguirre, D. P. Catarello, J. M. Wozniak, and R. D. Skeel. J. Chem. Phys., 114(4), 2090-2098, 2000.
[7] Understanding Molecular Simulation, from algorithms to Applications. D. Frenkel and B. Smith. Academic Press, London 2002.
[8] Computer Simulation of Liquids. M. P. Allen and D. J. Tildesley. Oxford Science Publications, New York 1987.

Page created and maintained by Marcos Sotomayor.