Gaussian accelerated molecular dynamics (GaMD): principles and applications

Abstract

Gaussian accelerated molecular dynamics (GaMD) is a robust computational method for simultaneous unconstrained enhanced sampling and free energy calculations of biomolecules. It works by adding a harmonic boost potential to smooth biomolecular potential energy surface and reduce energy barriers. GaMD greatly accelerates biomolecular simulations by orders of magnitude. Without the need to set predefined reaction coordinates or collective variables, GaMD provides unconstrained enhanced sampling and is advantageous for simulating complex biological processes. The GaMD boost potential exhibits a Gaussian distribution, thereby allowing for energetic reweighting via cumulant expansion to the second order (i.e., "Gaussian approximation"). This leads to accurate reconstruction of free energy landscapes of biomolecules. Hybrid schemes with other enhanced sampling methods, such as the replica exchange GaMD (rex-GaMD) and replica exchange umbrella sampling GaMD (GaREUS), have also been introduced, further improving sampling and free energy calculations. Recently, new "selective GaMD" algorithms including the ligand GaMD (LiGaMD) and peptide GaMD (Pep-GaMD) enabled microsecond simulations to capture repetitive dissociation and binding of small-molecule ligands and highly flexible peptides. The simulations then allowed highly efficient quantitative characterization of the ligand/peptide binding thermodynamics and kinetics. Taken together, GaMD and its innovative variants are applicable to simulate a wide variety of biomolecular dynamics, including protein folding, conformational changes and allostery, ligand binding, peptide binding, protein-protein/nucleic acid/carbohydrate interactions, and carbohydrate/nucleic acid interactions. In this review, we present principles of the GaMD algorithms and recent applications in biomolecular simulations and drug design.

Figure 1.

Scheme illustration of Gaussian accelerated molecular dynamics (GaMD). When the threshold energy is set to the maximum potential (E=Vmax), the system potential energy surface is smoothened by adding a harmonic boost potential that follows Gaussian distribution. The coefficient k₀ in the range of 0−1 determines the magnitude of the applied boost potential. With greater k₀, higher boost potential is added to the original energy surface in conventional molecular dynamics (cMD), which provides enhanced sampling of biomolecules across decreased energy barriers. Adapted with permission from Miao et al. (2015). Copyright 2015 American Chemical Society. https://pubs.acs.org/doi/abs/10.1021/acs.jctc.5b00436. Further permissions related to the material excerpted should be directed to the American Chemical Society.

Figure 2.

(A) Overview of the Streptococcus pyogenes CRISPR-Cas9 system. The Cas9 protein is represented in molecular surface, showing individual domains in different colors. The RNA (yellow), target DNA (TS, violet), and non-target DNA (NTS, cyan) are also shown. (B) Energetic landscape associated to the conformational transition of the Cas9 protein from the apo form to the RNA-bound state, computed using GaMD. The potential of mean force (PMF), which describes the free energy landscape, was computed along the E945-D435 FRET distance and the root mean square deviation (RMSD) with respect to the apo state. The simulations identified three minima: M1 corresponding to the crystallographic apo, M2 that is the RNA-bound structure and M3, which is an intermediate characterized by the solvent exposure of an arginine-rich helix. For selected states, the ensemble averaged electrostatic potential has been computed, revealing the formation of a positively charged cavity (blue) in the intermediate states. Adapted with permission from Palermo et al. (2017). https://www.pnas.org/content/114/28/7260. Copyright 2017 National Academy of Sciences.

Figure 3.

(A) Conformations of the HNH domain (green) in its inactive (#1, #2), pre-active (#3) and active (#4) states, as experimentally determined through single-molecule FRET and structural approaches (top panel). The free energy landscape (i.e., Potential of Mean Force, PMF) associated to the conformational changes of the HNH domain from its inactive to active states is shown in the bottom panel. The minima correspond to the four states experimentally found (top). The PMF was computed along the S867-S355 and N1054-S867 FRET distances. Adapted with permission from Palermo et al. (2017). https://www.pnas.org/content/114/28/7260. Copyright 2017 National Academy of Sciences. (B) Conformational change of the HNH domain from its pre-active conformation (captured in the PDB ID: 5F9R, left) to the active state identified through GaMD (right). The active state displays the catalytic residue H983 close to the DNA target strand (TS). The distance between the catalytic H840 and the scissile phosphate (H840–PDNA) has been computed along ~400 ns GaMD (central panel) and ~16 μs of continuous MD using the specialized supercomputer Anton-2 (bottom panel). The black dashed line indicates the pre-active conformation (PDB ID: 5F9R) used as a starting point for MD simulations, while the magenta dashed line indicates the active conformation more recently captured through cryo-EM (PDB ID: 6O0Y). Reprinted with permission from Palermo et al. (2018). Copyright 2018 Cambridge University Press. https://doi.org/10.1017/S0033583518000070.

Figure 4.

(A) Extended opening of the RNA:DNA hybrid and newly formed interactions with the L2 loop (magenta) of the HNH domain (green), observed during GaMD simulations of CRISPR-Cas9 in the presence of four base pair mismatches at the RNA:DNA hybrid ends. (B) RNA:DNA minor groove width computed along MD simulations of CRISPR-Cas9 bound to an on-target DNA (black) and in the presence of one to four mismatches at the hybrid ends. A vertical bar indicates the experimental minor groove width (i.e., 11 Å from X-ray crystallography). The minor groove width has been measured at the level of base pair 17 (shown on the right). Adapted with permission from Ricci et al. (2019). Copyright 2019 American Chemical Society. https://pubs.acs.org/doi/full/10.1021/acscentsci.9b00020. Further permissions related to the material excerpted should be directed to the American Chemical Society.

Figure 5.

Pep-GaMD simulations have captured repetitive dissociation and binding of three model peptides to the SH3 domains: (A-C) X-ray structures of the SH3 domains bound by peptides (A) “PAMPAR” (PDB: 1SSH), (B) “PPPALPPKK” (PDB: 1CKA) and (C) “PPPVPPRR” (PDB: 1CKB). The SH3 domains and peptides are shown in green and magenta cartoon, respectively. Key protein residues Asp19 and Trp40 in the 1SSH structure and Asp150 and Trp169 in the 1CKA and 1CKB structures, and peptide residues Arg10 in the 1SSH structure, Lys8 in the 1CKA structure and Arg7 in the 1CKB structure are highlighted in sticks. The “N” and “C” labels denote the N-terminus and C-terminus of the peptides. (D–F) time courses of peptide backbone RMSDs relative to X-ray structures with the protein aligned calculated from three independent 1 μs Pep-GaMD simulations of the (D) 1SSH, (E) 1CKA and (F) 1CKB structures. (G–I) The corresponding PMF profiles of the peptide backbone RMSDs averaged over three Pep-GaMD simulations of the (G) 1SSH, (H) 1CKA and (I) 1CKB structures. Error bars are standard deviations of the free energy values calculated from three Pep-GaMD simulations. Reprinted from “Jinan Wang, Yinglong Miao, J Chem Phys 2020, 153:154109”, with the permission of AIP Publishing.

Figure 6.

Binding of agonist IXO and G_i protein mimetic nanobody Nb9–8 to the M2 muscarinic GPCR was captured in one of five GaMD simulations: (A) Trajectories of a nitrogen atom in IXO (beads) and the β8 strand of Nb9–8 (ribbons) colored by simulation time in a blue (0 ns)–white (2250 ns)–red (4500 ns) scale. (B) RMSDs of the IXO and Nb9–8 relative to the X-ray structure, Tyr104^3.33-Tyr403^6.51-Tyr426^7.39 triangle perimeter and Arg121^3.50-Thr386^6.34 distance calculated from the simulation. Dashed lines indicate X-ray structural values of the M2 receptor (3UON: green and 4MQS: red). (C) Binding pose of IXO (spheres) in the receptor extracellular vestibule with 13.84 Å RMSD relative to the X-ray conformation (yellow spheres). Residues found within 5 Å of IXO are highlighted in sticks. (D) Binding of Nb9–8 (cyan), which exhibits only 2.48 Å RMSD in the protein core (the β2, β3, β6, β7 and β8 strands). X-ray conformations of the M2 receptor and nanobody are shown in orange and purple ribbons, respectively. Adapted with permission from Miao et al. (2018). https://www.pnas.org/content/115/12/3036. Copyright 2018 National Academy of Sciences.

Figure 7.

2D potential of mean force (PMF) profiles of the (A) ADO-bound A₁AR-G_i, (B) NECA bound A_2AAR-G_s, (C) NECA bound A_2AAR-G_i and (D) ADO bound A₁AR-G_s complex systems regarding the agonist RMSD relative to the cryo-EM conformation and AR:NPxxY-G:α5 distance. The white triangles indicate the cryo-EM or simulation starting structures. Summary of specific AR-G protein interactions: (E) the ADO-bound A₁AR prefers to bind the G_i protein to the G_s. The latter could not stabilize binding of agonist ADO in the A₁AR and tended to dissociate from the receptor. (F) The A_2AAR could bind both the G_s and G_i proteins, which adopted distinct conformations in the complexes. Adapted with permission from Wang et al. (2019). Copyright 2019 American Chemical Society. https://pubs.acs.org/doi/10.1021/acs.jpcb.9b04867. Further permissions related to the material excerpted should be directed to the American Chemical Society.

Figure 8.

GaMD simulations revealed the activation and its ε cleavage mechanisms of γ-secretase in the wildtype and mutant APP substrates. Summary of the (A) inactive cryo-EM, (B) active (wildtype), and (C) shifted active (M51F) conformational states of the APP-bound γ-secretase. Distinct AICD products were generated from the wildtype and M51F mutant APP. GaMD free energy profiles of (D) wildtype and (E) M51F APP-bound γ-secretase regarding the Asp257:Cγ–Asp385:Cγ and Asp257:protonated O–Leu49:O distances. Adapted with permission from Bhattarai et al. (2020). Copyright 2020 American Chemical Society. https://pubs.acs.org/doi/abs/10.1021/acscentsci.0c00296. Further permissions related to the material excerpted should be directed to the American Chemical Society.