Advanced techniques for constrained internal coordinate molecular dynamics

Abstract

Internal coordinate molecular dynamics (ICMD) methods provide a more natural description of a protein by using bond, angle, and torsional coordinates instead of a Cartesian coordinate representation. Freezing high-frequency bonds and angles in the ICMD model gives rise to constrained ICMD (CICMD) models. There are several theoretical aspects that need to be developed to make the CICMD method robust and widely usable. In this article, we have designed a new framework for (1) initializing velocities for nonindependent CICMD coordinates, (2) efficient computation of center of mass velocity during CICMD simulations, (3) using advanced integrators such as Runge-Kutta, Lobatto, and adaptive CVODE for CICMD simulations, and (4) cancelling out the "flying ice cube effect" that sometimes arises in Nosé-Hoover dynamics. The Generalized Newton-Euler Inverse Mass Operator (GNEIMO) method is an implementation of a CICMD method that we have developed to study protein dynamics. GNEIMO allows for a hierarchy of coarse-grained simulation models based on the ability to rigidly constrain any group of atoms. In this article, we perform tests on the Lobatto and Runge-Kutta integrators to determine optimal simulation parameters. We also implement an adaptive coarse-graining tool using the GNEIMO Python interface. This tool enables the secondary structure-guided "freezing and thawing" of degrees of freedom in the molecule on the fly during molecular dynamics simulations and is shown to fold four proteins to their native topologies. With these advancements, we envision the use of the GNEIMO method in protein structure prediction, structure refinement, and in studying domain motion.

Figure 1

Plot of center of mass(CM) kinetic energy as a function of time for a) cluster model(Cl) using integration time-step 20fs and b) all atom(AA) Cartesian simulations using integration time-step 1fs. The kinetic energy of CM from simulations with GB/SA solvation (red line) and vacuum simulation (blue dashed line) for both cluster and all-atom models have been shown here. The formula c exp(−2 ln s_t) for a fitted value of c is shown as squares and triangles for cases with and without GB/SA solvation respectively. All runs shown here use RK4 integrator.

Figure 2

Plot of KE_CM at the end of 90ps for different integration time-steps with GB/SA solvation (red full lines) and in vacuum (dashed blue curve). Nosé-Hoover thermostat at a bath temperature of 300K and bath relaxation constant 500fs was used with RK4 integrator.

Figure 3

Colored representation of the automated GNEIMO clustering scheme. Each group of same-colored atoms represents a rigid “cluster” which is connected to its neighbors via torsional hinges.

Figure 4

(a): Backbone CRMSD histogram of a “Dynamic Clustering” replica-exchange simulation of 1BDD, beginning from an extended conformation containing only predicted secondary structure elements. Helices that are treated as rigid bodies are shown as broad ribbons. (b), (c), and (d): Backbone CRMSD histograms of dynamic clustering simulations of proteins from predicted helical structure.

Figure 5

(a-c) The crystal structures of the three proteins of various sizes used for equilibrium dynamics simulations and (d-g) the crystal structures of the four proteins (captioned with the experimentally-resolved residue subrange used in simulations) used for ab-initio structure prediction.

Figure 6

Standard deviation of temperature vs. timestep size over a 5ns simulation for various molecule-integrator combinations.

Figure 7

Average CRMSD in coordinates vs. timestep size over a 5ns simulation for various molecule-integrator combinations.

Figure 8

Simulation RMSF - B factor derived RMSF for various proteins, integrators, and timesteps. RMSFs are calculated from the entire trajectory of each 5ns simulation. A value of 0 indicates that the RMSF of the residue on the x-axis observed in a simulation using the timestep on the y-axis matches the RMSF derived from the crystallographic B-factor. All distances in Å.