Overcoming potential energy distortions in constrained internal coordinate molecular dynamics simulations

Abstract

The Internal Coordinate Molecular Dynamics (ICMD) method is an attractive molecular dynamics (MD) method for studying the dynamics of bonded systems such as proteins and polymers. It offers a simple venue for coarsening the dynamics model of a system at multiple hierarchical levels. For example, large scale protein dynamics can be studied using torsional dynamics, where large domains or helical structures can be treated as rigid bodies and the loops connecting them as flexible torsions. ICMD with such a dynamic model of the protein, combined with enhanced conformational sampling method such as temperature replica exchange, allows the sampling of large scale domain motion involving high energy barrier transitions. Once these large scale conformational transitions are sampled, all-torsion, or even all-atom, MD simulations can be carried out for the low energy conformations sampled via coarse grained ICMD to calculate the energetics of distinct conformations. Such hierarchical MD simulations can be carried out with standard all-atom forcefields without the need for compromising on the accuracy of the forces. Using constraints to treat bond lengths and bond angles as rigid can, however, distort the potential energy landscape of the system and reduce the number of dihedral transitions as well as conformational sampling. We present here a two-part solution to overcome such distortions of the potential energy landscape with ICMD models. To alleviate the intrinsic distortion that stems from the reduced phase space in torsional MD, we use the Fixman compensating potential. To additionally alleviate the extrinsic distortion that arises from the coupling between the dihedral angles and bond angles within a force field, we propose a hybrid ICMD method that allows the selective relaxing of bond angles. This hybrid ICMD method bridges the gap between all-atom MD and torsional MD. We demonstrate with examples that these methods together offer a solution to eliminate the potential energy distortions encountered in constrained ICMD simulations of peptide molecules.

FIG. 1.

Harmonic dihedral potential with k_α = 0.30 kcal and α₀ = 0.

FIG. 2.

(a) Torsional probability density function for C4 for a harmonic torsional potential (Eq. (15)) with the barrier center at α₀ = 90° and T = 800 K. (b) RMS difference from the FLEXIBLE pdf for the TMD and FIXMAN simulations for various temperatures. (c) Transition state barrier crossing rates for C4 at α₀ = 90° and α₀ = 180°. (d) Ratio of the FLEXIBLE barrier crossing rates to the TMD and FIXMAN barrier crossing rates as a function of the location of the barrier peak.

FIG. 3.

(a) Torsional pdf’s for C4, with charges 0.2e and −0.2e at the terminal atoms, for k_θ = 95 kcal. (b) RMS deviation from the FLEXIBLE pdf as a function of k_θ. The FIXMAN and FLEXIBLE simulations agree only at very high values of k_θ. The TMD simulation remains divergent throughout.

FIG. 4.

(a) Structure of alanine dipeptide. The torsions ϕ = C–N–C_α–C and ψ = N–C_α–C–N effectively represent the conformational space explored. (b) Schema for opening a bond angle θ₁—we open the two terminal angles (dashed) to effectively open the non-terminal (solid) angle.

FIG. 5.

Free energy surface calculated for the backbone dihedral angles of alanine dipeptide: FLEXIBLE and TMD simulations at 300 K (top) and 800 K (bottom).

FIG. 6.

FES for FIXMAN, hybrid ICMD, and FLEXIBLE simulations of alanine dipeptide. The hybrid ICMD simulation has the bond angles C–N–C_α and N–C_α–C open.

FIG. 7.

Structures of (a) valine (VAL) dipeptide, (b) leucine (LEU) dipeptide, (c) isoleucine (ILE) dipeptide, (d) methionine (MET) dipeptide, (e) phenylalanine (PHE) dipeptide, (f) tryptophan (TRP dipeptide), (g) proline (PRO) dipeptide, and (h) tyrosine (TYR) dipeptide.

FIG. 8.

(a) Free energy surfaces for FIXMAN simulations(top), hybrid ICMD simulations (mid), and FLEXIBLE simulations (bottom) for dipeptide molecules. (b) Hellinger distances from the FLEXIBLE probability density functions for the FIXMAN and hybrid ICMD simulations. For VAL, LEU, MET, PHE, TRP, and TYR, the FIXMAN simulations are unable to sample the conformations in the first quadrant, the sampling of which requires that the C–N–C_α and N–C_α–C angles be kept open in the hybrid ICMD simulations. For ILE, opening the C–N–C_α and N–C_α–C angles does not improve the free energy surface.

FIG. 9.

Free energy surfaces for hybrid ICMD simulations with the open backbone angles C–N–C_α, N–C_α–C, and C_α–C–N for valine (VAL), methionine (MET), tyrosine (TYR), and isoleucine (ILE) dipeptides. With these angles open, the Hellinger distances to the FLEXIBLE probability density function are now 0.14, 0.13, 0.17, and 0.16, respectively.

FIG. 10.

Free energy surfaces for FIXMAN, hybrid ICMD, and FLEXIBLE simulations for proline dipeptide (PRO). The FIXMAN probability density function has a Hellinger distance of 0.45 from the FLEXIBLE density function. The hybrid ICMD simulation with the backbone angle N–C_α–C and the sidechain angle C_α–C_β–C_γ relaxed yields a Hellinger distance of 0.26 from the FLEXIBLE simulation.

FIG. 11.

Starting structure for the simulations of the 10-residue peptide CLN025. The PDB ID for the NMR structure of CLN025 is 2RVD.

FIG. 12.

Hellinger distances from the FLEXIBLE (ϕ, ψ) probability density functions for the FIXMAN and hybrid ICMD simulations for each of the residues in CLN025.

FIG. 13.

Free energy surfaces for hybrid ICMD simulations with open backbone angles C–N–C_α and N–C_α–C for alanine, leucine, phenylalanine, and tryptophan, with the additional angle C_α–C–N open for valine, isoleucine, methionine, and tyrosine dipeptides. For the proline dipeptide, the N–C_α–C and C_α–C_β–C_γ are kept open. The simulations use time steps of 5 fs each, with each simulation 50 ns long.