Accelerated Molecular Dynamics Simulations with the AMOEBA Polarizable Force Field on Graphics Processing Units.

The accelerated molecular dynamics (aMD) method has recently been shown to enhance the sampling of biomolecules in molecular dynamics (MD) simulations, often by several orders of magnitude. Here, we describe an implementation of the aMD method for the OpenMM application layer that takes full advantage of graphics processing units (GPUs) computing. The aMD method is shown to work in combination with the AMOEBA polarizable force field (AMOEBA-aMD), allowing the simulation of long time-scale events with a polarizable force field. Benchmarks are provided to show that the AMOEBA-aMD method is efficiently implemented and produces accurate results in its standard parametrization. For the BPTI protein, we demonstrate that the protein structure described with AMOEBA remains stable even on the extended time scales accessed at high levels of accelerations. For the DNA repair metalloenzyme endonuclease IV, we show that the use of the AMOEBA force field is a significant improvement over fixed charged models for describing the enzyme active-site. The new AMOEBA-aMD method is publicly available (http://wiki.simtk.org/openmm/VirtualRepository) and promises to be interesting for studying complex systems that can benefit from both the use of a polarizable force field and enhanced sampling.

Introduction

Molecular dynamics (MD) is one of the most prominent techniques used to study the dynamics and equilibrium properties of biomolecules. It solves Newton’s equations of motion for all atoms, using force fields that account for bonded and nonbonded atomic interactions. These force fields have been parametrized to agree with quantum mechanical simulations and experiments. Ever since the seminal first protein MD simulation of bovine pancreatic trypsin inhibitor (BPTI),[1] tremendous progress has been made both in terms of accuracy and sampling efficiency. Even larger macromolecular systems can now be simulated on the order of microseconds. However, many important biological processes occur on time scales far beyond this regime posing the need for millisecond or longer simulations. The development of the Anton supercomputer[2] certainly marked a large step in the right direction; however, that kind of computational fire power is not readily available to the average researcher. One of the most significant speed improvements of the last years was the utilization of graphics processing units (GPUs) for molecular dynamics. Efficient MD codes have been developed to aid with this endeavor.[3] Nonetheless, it will take many years before nonspecialized computer hardware can routinely simulate large molecular systems on the time scales of their slowest motions, which is needed to fully characterize their free energy landscapes.

In the meantime, enhanced sampling methods are one solution to this problem. A plethora of methods exist that perturb the underlying potential energy landscape to increase the chance of transitions between low energy states.[4] Three notable examples are the adaptive biasing force (ABF[5]), metadynamics[6] and driven adiabatic dynamics[7] methods. In addition to these methods, accelerated MD (aMD) is a promising technique that directly modifies regions of the potential energy landscape that are below a certain cutoff energy.[8] By energetically raising these regions and thus lowering barriers between energy wells, the landscape is perturbed to allow more frequent transitions between low energy states. This results in an enhanced sampling of the conformational space accessible to the simulation. aMD has been previously implemented for classical (nonpolarizable) simulations in AMBER[8a] and NAMD,[9] and for ab initio MD in CPMD.[10]

In addition to improvements in sampling (the precision problem), the force fields used to propagate atoms may also need to be improved (the accuracy problem). Moreover, sampling and force field accuracy are related challenges, since better sampling allows for more precise calculations of experimental quantities that in turn can lead to the detection of smaller force field deficiencies. Although the functional form of biomolecular force fields has remained largely unchanged since the first MD simulation of a protein conducted more than 35 years ago, recent studies highlight certain shortcomings present in all current force fields that are unlikely to be solved simply by tuning their parametrization. For example, none of the current methods seem to be able to accurately capture the temperature dependency of proteins secondary structure propensities.[11] In addition, polarizability may be essential to accurately simulate highly charged systems, such as nucleic acids in DNA and RNA.[12] For this reason, substantial efforts are under way to develop force fields with more sophisticated functional forms, including polarizable force fields that can capture the redistribution of electrons around each atom in response to changes in the environment. This approach may eventually replace conventional force fields in several key areas (e.g., the modeling of protein folding,[13] computer-aided drug design,[14] the calculation of ion channel properties,[15] or the description of allostery[16]).

The AMOEBA (Atomic Multipole Optimized Energetics for Biomolecular Applications) polarizable force field[17] was developed by Ponder and co-workers with the aim to move away from the well-established fixed point charge models and toward more expensive models that should allow a more accurate description of molecular properties. AMOEBA is one of the most widely used polarizable force fields to date, and it has been demonstrated to perform better than nonpolarizable force fields for describing solvation free energies of drug-like small molecules, and dynamical properties away from ambient condition.[18] However, despite the impressive acceleration provided by GPU computing,[19] the time performance of AMOEBA remains a significant limitation to its wider use to simulate large-scale biomolecular systems, since such simulations remain about 1–2 orders of magnitude more expensive than their nonpolarizable equivalents. For this reason, practical usage of AMOEBA has been limited, especially when target events are in micro- to millisecond time scales. AMOEBA has been previously combined with an enhanced-sampling free energy method—orthogonal space random walk (OSRW)—to achieve ab initio prediction of organic crystal structures and thermodynamics.[20] Similarly, combining AMOEBA with the enhanced sampling capabilities of aMD represents an attractive prospect. OpenMM is the perfect platform to implement the aMD method, as it already allows a wide user base to take advantage of GPU-accelerated MD simulations.[3c,19]

In this paper, we show that the aMD method and the AMOEBA force field can be synergistically combined in a way that conserves the accuracy of the polarizable force field, while significantly enhancing sampling. Using three examples, we show that (1) with larger energy boost levels, the AMEOBA-based aMD simulations can maintain structural stability, (2) the choice of aMD parameters is the same for AMOEBA and for nonpolarizable force fields, and (3) the AMOEBA-based aMD simulations allow metalloenzyme reactive sites to be better simulated.

Material and Methods

Implementation of aMD into OpenMM

Over the years, different variants of the aMD algorithm have been proposed.[8b] In the original implementation,[8a] aMD was used to only boost the dihedral potential since many protein conformational changes are governed by changes in the torsional degrees of freedom. Subsequent work introduced an aMD version that boosted the total potential.[21] In a more recent implementation, aMD uses a dual boosting approach that applies one boost potential only to the torsional terms and another separate boost potential to the total potential.[22] Here, we present the implementation of all three flavors of aMD into OpenMM using the CustomIntegrator, an integrator used to implement arbitrary, user defined integration algorithms. A detailed summary of the theory of aMD can be found in refs (8a and 22), so that we will only list the most important equations that were used for the implementation. In aMD, whenever the system’s potential energy, V(r), falls below a threshold energy, E, a boost potential, ΔV(r), is added to yield the biased potential V*(r).

Introducing the following abbreviationsThe total force on the system is calculated as a sum of constitutive forces, each obtained from the different components (“comp”) of the potential energy (e.g., dihedrals). A general expression can be derived for the aMD forces that is valid irrespectively of the choice made for the boosted components, as[8a,22]The overall force in the boosted system F* is then obtained aswhere F is the total force of the unaccelerated system and Fcomp is the unaccelerated force of the component that is accelerated (e.g., dihedrals). The forces for the three flavors of aMD then become

Note that, in the dual boost implementation, the dihedral force is only boosted once. This is in agreement with current aMD implementations in AMBER and NAMD. The Custom Integrator allows convenient implementation of aMD, by scaling the forces at every simulation step according to the boost. For this, the energy and forces associated with the component to be boosted have to be retrieved. Then, the total force is scaled depending on whether the system is above or below the threshold energy. The aMD CustomIntegrator is available online at http://wiki.simtk.org/openmm/VirtualRepository.

Systems and MD Simulations

Three different systems were studied: the alanine dipeptide (N-acetyl-N′-methyl-alaninamide), BPTI, and endonuclease IV. Alanine dipeptide was solvated with a water box containing a buffer region of 10 Å. To eliminate any steric clashes, 100 steps of conjugate gradient minimization were performed using SANDER.[23] A short 50 ps NPT simulation was used to heat up the system to a temperature of 300 K. BPTI and Endonuclease IV were prepared based on PDB entries 5PTI(24) and 1QTW,[25] respectively. Tleap[23] was used to neutralize the systems by adding 6 Cl– and 6 Na+ counterions respectively and solvate them with a water box. TIP3P water[26] was used for the nonpolarizable simulations. The fully solvated system contained 17 155 and 52 011 atoms, respectively. Minimization using SANDER[23] was carried out in two stages: an initial minimization of 1000 steps for the solvent and ions with the protein restrained by a force constant of 500 kcal/mol/Å2, followed by a 2500 step minimization of the entire system. A short initial 20 ps MD simulation with weak restraints (10 kcal/mol/Å2) on the protein residues was used to heat the system to a temperature of 300 K.

All production MD simulations for the alanine dipeptide, BPTI, and endonuclease IV were conducted in the NPT ensemble at 300 K, using the OpenMM Python-based application layer.[19] Periodic boundary conditions were used, along with a nonbonded interaction cutoff of 10 Å. Particle Mesh Ewald was used.[27] Nonpolarizable simulations were performed using the AMBER ff99SBildn force field[28] while polarizable simulations were performed using the AMOEBA force field.[18] For the AMBER simulations, bonds involving hydrogen atoms were constrained using the SHAKE algorithm,[29] allowing for a time step of 2 fs. A time step of 1 fs was used for the AMOEBA simulations, while no constraints were used. Mutual polarization was used along with an induced target ε of 0.01. Energy conservation was monitored and none of the simulations showed energy drift.

The acceleration level for the aMD method is defined in terms of E and α, where E is the threshold boost energy and α is a tuning parameter that determines the shape of the accelerated potential.[8a] To define the acceleration parameters, equilibrium MD simulations were first conducted, with the AMBER ff99SBildn and AMOEBA force fields, to calculate the average potential and dihedral energies. For the total potential/torsional energy term, the boost energy was defined as the average total/dihedral angle energy plus 3.5 times the number of residues in the solute, and α was defined as 20% of E. From previous experience, this empirical formula produces marked enhancement in the sampling, while allowing for the accurate reconstruction of the thermodynamic ensemble.[9,30] Despite the explicit inclusion of polarization and anharmonic corrections in the AMOEBA force field, we have found that similar parameters can be used to enhance the sampling in AMBER ff99SBildn and AMOEBA aMD simulations. Table 1 summarizes the boost parameters for our aMD simulations.

Table 1

Acceleration Parameters for the aMD Simulations

system	force field	αtot [kcal/mol]	Etot [kcal/mol]	αtorsion [kcal/mol]	Etorsion [kcal/mol]
alanine dipeptide	ff99SBildn	382.4	–5668.2	2.1	19.5
BPTI	ff99SBildn	3431.0	–50 335.4	40.6	821.3
	AMOEBA	3431.0	–56 857.9	40.6	310.1
endonuclease IV	ff99SBildn	10 402.2	–149 530.1	199.5	4,129.6
	AMOEBA	10 402.2	–169 944.8	199.5	1,971.8

Results and Discussion

Our results are presented as follows; first we use the workhorse of molecular simulations, the alanine dipeptide, to ensure that the aMD implementation does indeed significantly enhance the conformational space sampling. We then use BPTI to show that even at very high acceleration levels, AMOEBA-aMD simulations are stable, as indicated by low root-mean-square deviations (RMSD) with respect to the BPTI crystal structure. Finally, we study the enzyme endonuclease IV to showcase a large protein system, for which nonpolarizable simulations fail to accurately describe structural details.

aMD Samples All Conformations of the Alanine Dipeptide

To test our aMD implementation, four different types of simulations were performed for alanine dipeptide: 500 ns of conventional MD (cMD), 100 ns of dihedral boost aMD, 100 ns of total energy boost aMD, and 100 ns of dual boost aMD. The conformational space sampling is evaluated by deriving a free energy profile from occupancy of states in the φ–ψ dihedral backbone angle space. Pioneering theoretical work in the late 1980s and early 1990s established the four free energy minima accessible to alanine dipeptide in the φ–ψ conformational space:[31] the β region (φ < 0° and 120° < ψ < 180°, corresponding to a β-strand conformation), the αR region (φ < 0° and −60° < ψ < 60°, corresponding to a right-handed α-helical conformation), the αL region (φ ∼ 60° and ψ ∼ 60°, corresponding to a left-handed α-helical conformation) and the C7ax region (φ ∼60° and ψ ∼ −80°). For an excellent graphical representation of the free energy surface please refer to ref (32).

Figure 1 shows the φ–ψ free energy maps calculated from the cMD simulation and different simulation times. Despite the simplicity of the alanine dipeptide system, conventional MD cannot sample the entire φ–ψ dihedral backbone angle space. Within the first 50 ns, only the two deepest energy minima (β and αR) are sampled. Within 500 ns sampling is extended to cover roughly the αL region. No sampling of the C7ax region is seen within 500 ns of cMD simulation.

Figure 1

φ–ψ free energy maps (in kcal/mol) for alanine dipeptide calculated from cMD simulations of length (A) 50 ns, (B) 100 ns, (C) 300 ns, and (D) 500 ns.

In contrast, dual boost aMD samples the conformational space much more efficiently. Figure 2 shows the φ–ψ free energy maps for different simulation times of the dual boost aMD simulation. Even within 25 ns, all four regions (including the C7ax region) are sampled. After 100 ns of dual boost aMD the reweighted energy profile clearly matches previously reported profile generated by umbrella sampling.[32] This demonstrates the power of the aMD method to enhance sampling by orders of magnitude and its correct implementation in OpenMM. It is worth noting that an explicit boost of the dihedral energy is often beneficial for peptide and protein simulations. For instance, in this test, only boosting the total energy was not sufficient to sample the αL and C7ax regions within 100 ns (Supporting Information Figure 1), while the dihedral boost aMD samples the conformational space similar to the dual boost aMD (Supporting Information Figure 2).

Figure 2

Reweighted φ–ψ free energy maps (in kcal/mol) for alanine dipeptide calculated from dual boost aMD simulations of length (A) 25 ns, (B) 50 ns, (C) 75 ns, and (D) 100 ns.

AMOEBA-aMD Simulations of BPTI Do Not Exhibit Instabilities

The BPTI protein is often used to test new methods. In particular, it was used to test a previous implementation of aMD for the Amber code,[33] based on nonpolarizable AMBER force fields. It was shown that, at high acceleration levels, aMD simulations of ∼500 ns can replicate all the energy minima of BPTI that were observed by Shaw and co-workers using millisecond simulations,[34] thus leading to an acceleration in the sampling of 3–4 orders of magnitude.

Here, we were interested to show that BPTI is stable at similar high acceleration levels in our aMD OpenMM implementation, with both the AMBER ff99SBildn, and the AMOEBA force fields. For each force field, a total of four different simulations were performed: cMD, dihedral boost aMD, total energy boost aMD and dual boost aMD. The simulations were run for 100 ns, and 3 ns, for the AMBER ff99SBildn and AMOEBA force fields, respectively. Simulation speeds of ∼30 ns/day for AMBER and ∼300 ps/day for AMOEBA were observed. To allow for a fair comparison, we also performed a short 1 ns benchmark AMOEBA simulation using the RESPA integrator (allowing for a 2 fs time step for the multipole force evaluation, while maintaining a 1 fs time step for all other forces). The simulations ran at 632 ps/day. Thus, depending on the parameters, AMOEBA simulations are about 1–2 orders of magnitude slower than cMD; however, enhanced sampling with aMD has been previously shown for BPTI to enhance sampling by several orders of magnitude[33] and therefore should be able to bridge this gap. The BPTI protein stability was judged by its RMSD to the starting structure. Potential unfolding events or instabilities due to excessive boosting would show up as an increase in the RMSD. Figure 3 summarizes the RMSDs for the simulations and shows that none of the aMD simulations exceeds 3 Å in the RMSD, which is similar to the unaccelerated simulations.

Figure 3

RMSD vs time plots for simulations of BPTI. RMSD values over the course of the simulations are shown for eight different BPTI simulations: MD AMBER ff99SBildn, total energy boost aMD AMBER ff99SBildn, dual boost aMD AMBER ff99SBildn, dihedral boost aMD AMBER ff99SBildn, MD AMOEBA, total energy boost aMD AMOEBA, dual boost aMD AMOEBA, and dihedral boost aMD AMOEBA.

AMOEBA Improves the Description of Endonuclease IV

Next, we demonstrate that AMOEBA can achieve more accurate results than nonpolarizable force field simulations. For this, we chose to simulate the endonuclease IV metallo-enzyme. The enzyme has been determined at a very high resolution (1.02 Å) using X-ray crystallography.[25] Figure 4A,B shows the enzyme with its active site. It contains a stable three zinc cluster coordinated by histidine, aspartate, and glutamate residues. The zinc ions form a structure reminiscent of a right triangle with distinct relative distances. The trinuclear cluster has been previously studied with hybrid quantum mechanics/molecular mechanics (QM/MM) simulations and was found to be stable.[35] However, QM/MM simulations are typically time-consuming and can often probe only the picosecond time scale. Our initial hypothesis was that AMOEBA simulations should be able to more accurately describe the zinc cluster geometry than nonpolarizable AMBER simulations. To show this, we ran nonpolarizable and polarizable simulations on the protein and monitored the relative zinc distances. Again, depending on the simulation parameters, we observed about 1–2 orders of magnitude faster sampling for nonpolarizable force field simulations (∼9 ns/day, for AMBER, versus ∼90 ps/day, for AMOEBA without RESPA). However, the accuracy was improved with AMOEBA. For instance, Figure 5 shows the Zn–Zn distances over the course of the simulations. In all three independent AMBER ff99SBildn simulations, the relative positions of the zinc cluster were perturbed within 3 ns, demonstrating a failure of the nonpolarizable force field to correctly capture the coordination structure around these charged ions. In contrast, the AMOEBA force field maintained the zinc geometry over the course of the 3 ns simulation. Figure 4C,D summarizes the active site residues and zinc positions after 3 ns of simulation, superimposed with the crystal structure. The most pronounced perturbation in the nonpolarizable simulations is the displacement of Zn2 away from the active site, facilitated by the displacement of the coordinating side chain ASP 229, and transferring the coordination of Zn2 from HIS 182 to GLN 150.

Figure 4

Structure of endonuclease IV. (A) Ribbon representation of entire protein with its zinc cluster in the center. (B) Close-up view of the endonuclease IV active site in the high resolution crystal structure (1QTW, 1.02 Å resolution). The zinc ions are forming a right triangle geometry coordinated by histidine, aspartate, and glutamate residues. (C) Superposition of the crystal structure active site and the active site after 3 ns of AMBER ff99SBildn simulation. Crystal structure zinc is shown in gray, while the simulated zinc is shown in green. (D) Superposition of the crystal structure active site and the active site after 3 ns of AMOEBA simulation. Crystal structure zinc is shown in gray, while the simulated zinc is shown in red.

Figure 5

Zn–Zn distances in endonuclease IV over the course of the simulations. Endonuclease IV contains a three zinc cluster (Zn 286, 287, and 288). Distances between 286 and 287 (red), 286–288 (green), and 287–288 (blue) are shown. Horizontal dotted colored lines indicate the distance between the respective ion pair in the high-resolution crystal structure. Three independent MD simulations with the AMBER ff99SBildn force field are shown on top. A simulation with the AMOEBA force field is shown in the bottom.

Conclusions

In this paper, we presented the synergistic combination of the enhanced sampling method aMD with the polarizable force field AMOEBA. The implementation of the aMD sampling method into OpenMM should provide a fast and simple avenue for studying biomolecular systems that require extensive sampling (e.g., micro- to milliseconds), in addition to a polarizable force field. The tests presented here suggest that it is possible to maintain the accuracy of the AMOEBA polarizable force field while significantly improving the simulation sampling efficiency. At the acceleration levels used, the sampling efficiency for medium size protein simulations is typically increased by 2–3 orders in magnitude.[33,36] Finally, the superior performance of AMOEBA was demonstrated for describing the active-site of a zinc metalloenzyme, which is consistent with previous studies.[37] The AMOEBA-aMD implementation takes full advantage of GPU computing and is publicly available in OpenMM (https://simtk.org/home/openmm).