Protecting High Energy Barriers: A New Equation to Regulate Boost Energy in Accelerated Molecular Dynamics Simulations

Abstract

Molecular dynamics (MD) is one of the most common tools in computational chemistry. Recently, our group has employed accelerated molecular dynamics (aMD) to improve the conformational sampling over conventional molecular dynamics techniques. In the original aMD implementation, sampling is greatly improved by raising energy wells below a predefined energy level. Recently, our group presented an alternative aMD implementation where simulations are accelerated by lowering energy barriers of the potential energy surface. When coupled with thermodynamic integration simulations, this implementation showed very promising results. However, when applied to large systems, such as proteins, the simulation tends to be biased to high energy regions of the potential landscape. The reason for this behavior lies in the boost equation used since the highest energy barriers are dramatically more affected than the lower ones. To address this issue, in this work, we present a new boost equation that prevents oversampling of unfavorable high energy conformational states. The new boost potential provides not only better recovery of statistics throughout the simulation but also enhanced sampling of statistically relevant regions in explicit solvent MD simulations.

Figure 1

Hypothetical one-dimensional potential representing the effect of ΔV^c. In all charts, α₁ = 200 and E₁ = −250. The upper and lower solid black lines represent the original potential and the modified potential generated with ΔV^b, respectively. This color scheme is used throughout. (A) Effects of different parameters E₂ (dashed colored lines) on the modified potential generated with ΔV^c (solid colored lines). (B) Boost levels ΔV^b (solid black line) and ΔV^c (colored lines) as V(r) moves away from E₁. For both A and B, α₂=15 and E₂ = −100 (red), 0 (blue), and 150 (green). (C) Effect of varying α₂ parameter on ΔV^c: α₂= 3 (red), 15 (blue), and 75 (green) with E₂ = 0.

Figure 2

Alanine dipeptide simulation results. (A) Ψ angle values obtained from cMD and five different aMD simulations. From top to bottom, aMD parameters were set to E₂ = E₁ + 15 and α₁ = 5, E₂ = E₁ + 20 and α₁ = 5, E₂ = E₁ + 25 and α₁ = 5, E₂ = E₁ + 25 and α₁ = 2.5, and E₂ = E₁ + 25 and α₁ = 1.25. In all simulations, E₁ and α₂ were set to 10 and 5, respectively. Weighted free energy density plots obtained from cMD (B, C, D) and aMD with ΔV^c (E). All values are in kcal/mol.

Figure 3

ψ and Ω angle values obtained from aMD simulations with boost potentials ΔV^b and ΔV^c. In all simulations, E₁ = 10.0 and α₁ were set as shown on the far right. Additional parameters for aMD with ΔV^c were set to E₂ = E₁ + 15 and α₂ = 5. All values are in kcal/mol.

Figure 4

Decalanine Φ–Ψ angles distribution obtained from cMD and aMD simulations. For the aMD simulations with ΔV^c parameters were set to E₁ = 74, E₂ = E₁ + 25, α₂ = 5, α₁ = 30 (C), and α₁ = 15 (D). All values are in kcal/mol.

Figure 5

Principle component analysis obtained from decalanine MD simulations. (A) 50 ns of aMD simulation with ΔV^c. Parameters were set to E₁ = 74, E₂ = E₁ + 25, α₂ = 5, and α₁ = 15, same as in Figure 4D. (B) 50 ns of cMD simulation and (C) 350 ns of cMD simulation. Structures 1, 2, and 3 shown in yellow represent relevant populated states in PC subspace sampled by aMD and cMD.

Figure 6

(A) Relative free energy of binding between Ac_2-L-Lys_-D-Ala_–-D-Ala and Ac_2-L-Lys_-D-Ala_-D-Lac to vancomycin calculated from cMD (solid black line) and aMD with ΔV^c (solid red line). A dashed line displays the experimental value, 4.1 kcal/mol. (B) Cumulative free energy curves calculated from simulations of 600 ps (left) and 1000 ps (right) per λ point. The ∗ shows points where there is no overlapping between error bars.