Binding Analysis Using Accelerated Molecular Dynamics Simulations and Future Perspectives

Introduction

Cellular functions such as gene translation, protein folding, ligand binding and signaling often require conformational transitions of biological macromolecules such as proteins, nucleic acids and lipids into various low-energy states.^1–4 Free energy landscapes of biomolecules essentially determine their structure, dynamics and functions.^5,6 It is important to comprehend the dynamic nature of biomolecules. Experimental techniques including nuclear magnetic resonance (NMR), X-ray crystallography and cryo-electron microscopy (cryo-EM) have greatly advanced structural characterizations of the biomolecular binding processes, but these efforts are time-consuming and very expensive. Moreover, these techniques often provide a rather static picture of the biomolecular interactions and conformations, whereas biomolecules are dynamic in vivo to perform the necessary cellular functions. It remains challenging and critical for us to understand the mechanisms, thermodynamics and kinetics of biomolecular recognition.

Molecular dynamics (MD) is a powerful “computational microscope” that helps us to visualize atomistic time evolution of biomolecules.⁷ Remarkable advances have been made in computing hardware (eg, Anton supercomputers and GPUs) and algorithm developments over the past several decades for running longer and cheaper MD simulations.^8–12 With microsecond-timescale conventional MD (cMD) simulations, spontaneous ligand/drug binding was successfully captured to the Src protein kinase,¹³ the β₁- and β₂-adrenergic receptors (β₁AR and β₁AR)¹⁴ and the M₂ and M₃ muscarinic receptors¹⁵ by the DE Shaw Research Group using the Anton supercomputer, for which the simulation predicted ligand binding poses closely matched those determined in the X-ray crystal structures. Moreover, hundreds-of-nanosecond cMD simulations were combined in Markov State Models to characterize the benzamidine ligand binding to the trypsin model protein by the De Fabritiis group.¹⁶

In addition, microsecond-timescale cMD simulations have been successful to capture peptide binding in several peptide-protein complexes, including proline-rich motif (PRM)-SH3,¹⁷ BAD-BclX_L and p53-MDM2.¹⁸ The “folding-upon-binding” mechanism of peptides was investigated using a measles virus nucleoprotein fragment and X domain of its phosphoprotein complex.¹⁹ Anton supercomputer was used to run ~200 μs MD simulations at 400 K elevated temperature. The simulations could capture 70 binding and unbinding events of the peptide fragment to the protein domain. Therefore, cMD has been successful to capture binding of ligands and peptides on the microsecond timescales. However, cMD with limited simulation timescales still lacks the power to study slower biological processes that would occur over longer timescales (eg, milliseconds and beyond).¹⁵

Metadynamics,^20,21 umbrella sampling,^22,23 adaptive biasing force (ABF) calculations^24,25 and conformational flooding^26,27 represent more advanced simulation methods for enhanced sampling of biomolecules developed over the past few decades. Metadynamics, in particular, has enabled investigating ligand binding to receptors in terms of the kinetic rates of binding and unbinding, residence time^28,29 and binding free energies.^30–33 The binding rate constant of peptide p53 to MDM2 protein predicted by ~120 μs weighted ensemble simulations was in good agreement with the experimentally calculated values.³⁴ Enhanced sampling implemented with a combination of long microsecond cMD and Hamiltonian replica exchange method could simulate rare events such as binding and unbinding of the PMI peptide inhibitor-MDM2 oncoprotein fragment complex³⁵ and the resulting kinetic rates agreed well with experiments.³⁵ Moreover, infrequent metadynamics could predict peptide association and dissociation rates of the p53-MDM2 peptide-protein complex, whereas, bias-exchange metadynamics could calculate binding free energy of the same system with 27 μs simulations.³⁶ Although these methods have shown remarkable improvements in sampling and capturing rare events that happen over exceedingly long timescales, they often pose a limitation for the requirement of predefined reaction coordinates and collective variables (CVs).^27,37,38

To address the above limitations, more enhanced sampling methods that do not require CVs have been developed,³⁹ including replica exchange MD,^40,41 self-guided MD,^42–45 and accelerated MD (aMD),^46–55 etc. In particular, aMD adds a boost potential ΔV to smooth the system potential energy surface to sample biomolecular transitions across reduced energy barriers.^55,56 The aMD simulations have been successful in capturing spontaneous ligand binding to a muscarinic G-protein-coupled receptor (GPCR).⁵⁷ However, aMD poses a drawback of large energetic noise during reweighting because of very high boost potential that is normally on the order of tens to hundreds of kcal/mol.⁵⁸ Recovering the original free energy landscapes represents a major limitation of aMD simulations.^57,59,60 To address this issue, Gaussian aMD (GaMD)^61,62 has been developed using a harmonic function to construct the boost potential that exhibits a Gaussian distribution. Through cumulant expansion to the second order, accurate reweighting can be achieved from GaMD simulations.⁶³ Without the need of predefined CVs, GaMD provides “unconstrained” enhanced sampling and can be easier to use for simulations of complex biomolecules than CV-based enhanced sampling methods. GaMD has been also shown to clearly identify distinct low-energy intermediate states during large protein conformational changes, ligand binding and protein–protein interactions. These intermediate states may not be sufficiently sampled in the CV-biasing simulation methods. On the other hand, relatively longer GaMD simulations could be needed in order to achieve sufficient sampling. Nevertheless, hundreds-of-nanosecond to microsecond timescale GaMD simulations have been successfully applied to sample large conformational changes of various biomolecules, including fast-folding proteins,^61,62,64 the CRISPR-Cas9 gene-editing system,^65,66 GPCRs,^62,67–71 protein kinases,^72,73 ligand binding^{61,62,67–71,74,75} and protein-protein/membrane/nucleic acid interactions.^66,76–80

The principles of GaMD and notably applications in understanding drug binding pathways were reviewed recently.^74,81 Here, we will review recent method developments of GaMD, Ligand GaMD (LiGaMD)⁸² and Peptide GaMD (Pep-GaMD),⁸³ which have greatly expanded our capabilities to simulate biomolecular binding processes. More recently, spontaneous binding of various biomolecules to their receptors such as drug binding to membrane proteins like chemokine GPCRs,⁷¹ human angiotensin-converting enzyme 2 (ACE2),⁸⁴ and RNA-protein binding^85,86 has been demonstrated by GaMD simulations. Furthermore, microsecond LiGaMD and Pep-GaMD simulations were able to simulate repetitive binding and dissociation of small-molecule ligands⁸² and highly flexible peptides,⁸³ and thus enabled ligand/peptide binding thermodynamics and kinetics calculations. We will also present relevant application studies in simulations of important drug targets and future perspectives for rational computer-aided drug design.

Materials and Methods

Gaussian Accelerated Molecular Dynamics (GaMD)

GaMD is an unconstrained enhanced sampling technique that reduces the system energy barriers by adding a harmonic boost potential to the system potential energy terms.⁶¹ The added boost potential smooths the system potential energy surface and allows biomolecules to sample large conformational space. The method is briefly illustrated here. For a system with N atoms at positions , when the system potential is lower than a threshold energy E, a boost potential is added as:

(1)

where k is the harmonic force constant. The two adjustable parameters E and k are automatically determined based on three enhanced sampling principles.⁶¹ The reference energy needs to be set in the following range:

(2)

where V_max and V_min are the system minimum and maximum potential energies. To ensure that Eqn. (2) is valid, k has to satisfy: Let us define , then . The standard deviation of needs to be small enough (ie, narrow distribution) to ensure proper energetic reweighting:⁶³ where and are the average and standard deviation of the system potential energies, is the standard deviation of with as a user-specified upper limit (eg, 10k_BT) for proper reweighting. When E is set to the lower bound E=V_max, can be calculated as:

(3)

Alternatively, when the threshold energy E is set to its upper bound , is set to:

(4)

if is found to be between 0 and 1. Otherwise, is calculated using Eqn. (3).

Like aMD, GaMD also provides options to add only the total potential boost, only dihedral potential boost, or the dual potential boost. The dual-boost GaMD (GaMD_Dual) simulation generally provides higher acceleration than the other two types of simulations for enhanced sampling.⁵⁶ The simulation parameters comprise of the threshold energy values and the effective harmonic force constants for the total and dihedral potential boost, respectively.

Energetic Reweighting of GaMD for Free Energy Calculations

For energetic reweighting to calculate the potential of mean force (PMF), the probability distribution along a reaction coordinate is written as . Given the boost potential of each frame, can be reweighted to recover the canonical ensemble distribution, , as:

(5)

where M is the number of bins, and is the ensemble-averaged Boltzmann factor of for simulation frames found in the jth bin. The ensemble-averaged reweighting factor can be approximated using cumulant expansion:

(6)

where the first two cumulants are given by:

(7)

The boost potential obtained from GaMD simulations usually follows near-Gaussian distribution.⁷⁴ Cumulant expansion to the second order thus provides a good approximation for computing the reweighting factor.^61,63 The reweighted free energy is calculated as:

(8)

where is the modified free energy obtained from GaMD simulation and is a constant. A free Python scripts toolkit “PyReweighting”⁶³ is developed and distributed at http://miao.compbio.ku.edu/PyReweighting/.

Ligand Gaussian Accelerated Molecular Dynamics (LiGaMD)

LiGaMD is developed based on GaMD⁶¹ and efficiently sample ligand binding and unbinding by selectively boosting non-bonded interaction potential energy between the ligand-receptor and ligand/environment. A detailed description of the method has been provided in the previous LiGaMD report.⁸² A brief summary is provided below. Let us assume a system containing protein P with ligand L and a biological environment E. The system has N atoms with their coordinates and momenta . The system Hamiltonian can be expressed as:

(9)

where and are the system kinetic and total potential energies, respectively. Now, the potential energy can be expressed as:

(10)

where , and represent the bonded potential energies in protein P, ligand L and environment E, respectively. , and are the self non-bonded potential energies in protein P, ligand L and environment E, respectively. , and represent the non-bonded interaction energies between P-L, P-E and L-E, respectively. As we know, in the classical molecular mechanics force fields,^87–90 the non-bonded potential energies are usually calculated as:

(11)

Where and represent the system electrostatic and van der Waals potential energies. Let us assume, while the ligand binds to the receptor, it mainly involves the non-bonded interaction energies, . In accordance with the previous GaMD algorithm, we will add boost potential to the essential ligand non-bonded potential energy in the LiGaMD method:

(12)

where k_L,nb represents the harmonic constant and E_L,nb represents the threshold energy for applying boost potential. These parameters can be based on three enhanced sampling principles as described in GaMD.⁶¹

Next, one can add multiple ligand molecules in the solvent to facilitate ligand binding to proteins in MD simulations.^14,70 This is based on the fact that the ligand-binding rate constant k_on is inversely proportional to the ligand concentration. The higher the ligand concentration, the faster the ligand binds, provided that the ligand concentration is still within its solubility limit. In addition to selectively boosting the bound ligand, another boost potential could thus be applied on the unbound ligand molecules, protein and solvent to facilitate both ligand dissociation and rebinding. The second boost potential is calculated using the total system potential energy other than the non-bonded potential energy of the bound ligand in Eqn. (10) as:

(13)

where E_D and k_D are the corresponding threshold energy for applying the second boost potential and the harmonic constant, respectively. This leads to dual-boost LiGaMD (LiGaMD_Dual) with the total boost potential .

Peptide Gaussian Accelerated Molecular Dynamics (Pep-GaMD)

Peptides often undergo large conformational changes during binding to target proteins, being distinct from small-molecule ligand binding or protein–protein interactions (PPIs). Another algorithm called peptide GaMD or “Pep-GaMD” has been developed to enhance sampling of peptide binding.⁸³ In Pep-GaMD, we consider a system of peptide L binding to a protein P in a biological environment E. We decompose the potential energy into similar terms as in Eqn. (10). Presumably, peptide binding mainly involves in both the bonded and non-bonded interaction energies of the peptide since peptides often undergo large conformational changes during binding to the target proteins. Thus, the essential peptide potential energy is In Pep-GaMD, we add boost potential selectively to the essential peptide potential energy according to the GaMD algorithm:

(15)

where E_L is the threshold energy for applying boost potential and k_L is the harmonic constant. The Pep-GaMD simulation parameters are derived similarly as in GaMD.

A second boost is added to the interaction between the protein and the environment/solvent to accelerate the peptide binding, along with the selective boost added to the peptide. This boost represents the total system potential energy without the essential peptide potential energy included:

(16)

Where V_D represents the total system potential energy without the essential peptide potential energy included, E_D represents the second boost potential threshold energy and k_D represents the harmonic constant. Hence, this contributes to the dual-boost Pep-GaMD as the total boost potential .

Applications

Mechanism of Drug Binding to Chemokine Receptor

Chemokine receptors are key GPCRs with involvement in many biological functions including immune cell response, growth and development of bodily structure like cardiac tissues and central nervous system, and several cancers.^91,92 In particular, CXCR4 facilitates the entry of HIV into host immune cells.^93–95 Plerixafor (PLX) was discovered as the first selective small-molecule inhibitor of CXCR4 (Mozobil^®, AMD3100), however, is suggested to have multiple binding modes to CXCR4.^96–98 Here, GaMD simulations were performed to study the pathway and binding mechanism of the PLX drug to the CXCR4 receptor.⁷¹

The X-ray structure of CXCR4 (PDB: 3ODU) with the highest resolution (2.5 Å) was used to build the simulation system. The conformation of PLX with positively charged N1, N3, N4 and N7 of the bicyclam rings at physiological pH was also added to the system. GaMD simulations were performed for ligand in both unbound and bound state with ten ligand molecules >15 Å away from CXCR4 (Figure 1A), and with the lowest energy bound conformation of PLX obtained from rigid-body docking using Autodock.⁹⁹ GaMD_Dual_NB boost scheme was used in which system nonbonded and dihedral energy terms were boosted. GaMD_Dual_NB was applied to five independent 800–1000 ns GaMD production runs to capture PLX drug binding to CXCR4 from an unbound state.

Figure 1 (A) Computational representation of the CXCR4 receptor (green ribbons) with 10 PLX drug molecules (magenta spheres) placed away in the solvent used in the simulation. The receptor was placed in a POPC lipid bilayer (sticks) and solvated in 0.15 M NaCl aqueous solution. (B) GaMD predicted binding pose of PLX (magenta sticks) at the orthosteric site of CXCR4 (green ribbons) with their interacting residues labeled and highlighted in green ball and sticks. Antagonist IT1t (yellow) and docking conformation of PLX (orange) are shown for reference. The seven transmembrane helices I–VII and three extracellular loops (ECL) 1–3 are labeled in the CXCR4 receptor; (C) Intermediate I1 state with interacting residues highlighted in green sticks, including D187ECL2 and D2626.58 that formed ionic interactions with PLX atoms N4 and N3, respectively; (D) Intermediate I2 state with interacting residues highlighted in green sticks, including D187ECL2, D1935.32 and D2626.58 that formed salt bridges with positively charged N7, N4 and N1 of PLX, respectively.

The PLX bound simulations showed stable binding of PLX in the CXCR4 orthosteric site. The simulations were used for hierarchical agglomerative clustering to obtain five clusters; the top-ranked cluster was selected as a reference, termed the “bound state”. Spontaneous binding of PLX was captured to the orthosteric site of the receptor during one of the five replicas. Complete binding of PLX was observed at ~480 ns with minimum RMSD of 2.76 Å relative to the X-ray structure of CXCR4 (PDB: 3ODU). The positively charged PLX formed stable salt bridge interactions with residues Asp97^2.63, Asp262^6.58 and Glu288^7.39 occupying both major and minor sub-pockets of the receptor (Figure 1B). Other simulations revealed important drug–receptor interactions significant in the recognition and binding of PLX to the CXCR4 receptor.

Free energy profiles of PLX binding to the CXCR4 receptor were calculated using the PLX RMSD and the distance between the PLX and the receptor as reaction coordinates of all five GaMD_Dual_NB simulations combined. The 2D PMF profile helped identify four low-energy conformational states of PLX including “Bound (B)”, “Intermediate 1 (I1)”, “Intermediate 2 (I2)” and “Unbound (U)”. In the I1 and I2 state, polar residues in the ECL2-TM5-TM6 region of the CXCR4, namely, D187^ECL2, D193^5.32 and D262^6.58 formed favourable interactions with the nitrogen atoms of PLX which are positively charged (Figure 1C and D). A single novel intermediate drug binding site was identified in the CXCR4 as the two intermediate states shared the common key residues including D187^ECL2, D193^5.32 and D262^6.58 because the symmetric ligand flips its conformation at the same site of the receptor in the ECL2-TM5-TM6. GaMD simulations also identified PLX drug binding pathway to the CXCR4. Starting from diffusion in the solvent, PLX binds to an intermediate site between ECL2-TM5-TM6 and then the final orthosteric site. This study reveals a novel allosteric site that is expected to facilitate drug design of CXCR4.

Mechanism of Ligand Recognition by Human ACE2 Receptor

Human ACE2 receptor has been of critical importance in biology and medicine. It is highly expressed in different organ epithelial cells such as heart, kidney, pancreas and lungs. Human ACE2 receptor is an excellent therapeutic target to severe acute respiratory syndrome coronavirus (SARS-Cov) as it serves as receptor for both SARS-CoV and SARS-CoV-2. It is critically important to study drug binding mechanism to this receptor, which remains elusive till today, for the development of effective treatment against SARS-CoV-2.

In order to study the molecular and structural basic of ligand binding to the human ACE2, LiGaMD simulations were performed using the X-ray structure of the inhibitor MLN-4760-bound ACE2 as model (PDB: 1R4L).⁸⁴ The AMBER ff19SB force field and GAFF-2 were used to model the protein and ligand, respectively. In addition to the MLN-4760 in the X-ray structure, 9 MLN-4760 molecules were added to the solvent to build the initial simulation systems. 64 ns LiGaMD equilibration followed by ten independent 700-2000 ns LiGaMD productions were performed. Dissociation of MLN-4760 from the ACE2 active site was captured during LiGaMD equilibration. Protein conformational change was observed upon ligand dissociation where the receptor subdomain I transited from “Closed” to “Open” state. Ligand binding and unbinding to the ACE2 active site was captured in three of the ten LiGaMD simulations (“Sim1”, “Sim2” and “Sim3”) with minimum ligand RMSD relative to the X-ray structure of ~0.99 Å (Figure 2A). Large-scale conformational changes of the ACE2 were sampled in LiGaMD. In addition, LiGaMD identified four low-energy conformations of the ACE2 receptor, viz. “Open”, “Partially Open”, “Closed”, and “Fully Closed”. The “Fully Closed” conformation was newly discovered through LiGaMD simulations, whereas, the three “Open”, “Partially Open” and “Closed” conformations were highly consistent with the experimentally resolved structures.

Particularly, in Sim2, within ~100 ns, one MLN-4760 initially interacted to the receptor interface between the α5 and 3₁₀ H4 helices. In ~100–160 ns, the ligand further entered the ACE2 active site through the cavity between the subdomain I and II, bound to the active site up to ~500 ns and then dissociated (Figure 2B). Ligand dissociation pathway was completely different than that of the ligand binding (Figure 2C and D). When the ligand dissociated, the conformation of receptor subdomain I was changed from “Closed” to “Open” through opening between the receptor helices α2 and α4.

Figure 2 (A) “Bound” state conformation of the ACE2 receptor and MLN-4760 inhibitor (green) in comparison to the X-ray conformation (Orange, PDB: 1R4L) (B) RMSD of MLN-4760 ligands relative to the bound X-ray conformation (PDB: 1R4L) calculated from “Sim 2” LiGaMD simulations. (C) Rotated views of ligand binding pathway with MLN-4760 central ring (lines) colored in a blue-white-red (BWR) color scale along simulation time. (D) Rotated views of ligand dissociation pathway with MLN-4760 central ring (lines) colored in a blue-white-red (BWR) color scale along simulation time. (E) Free energy profile of the ligand RMSD and interdomain distance highlighting nine low-energy conformational states: one Bound (B), four Intermediate (I-1, I-2, I-3, I-4), and four Unbound (U-1, U-2, U-3, U-4). Adapted with permission from Bhattarai A, Pawnikar S, Miao Y. Mechanism of ligand recognition by human ACE2 receptor. J Phys Chem Lett. 12:4814-4822. Copyright 2021 American Chemical Society.84

Using the ligand RMSD relative to the X-ray structure and interdomain distance as reaction coordinates, free energy landscape was calculated by combining all LiGaMD simulation trajectories. The free energy landscape identified nine low-energy conformations of the receptor-ligand complex, namely, “Bound (B)”, “Intermediate-1 (I-1)”, “Intermediate-2 (I-2)”, “Intermediate-3 (I-3)”, “Intermediate (I-4)”, “Unbound-1 (U-1)”, “Unbound-2 (U-2)”, “Unbound-3 (U-3)” and “Unbound-4 (U-4)” (Figure 2E). Ligand bound at the active site with minimum ligand RMSD of 0.99 Å was termed as the “Bound” conformation, in which the receptor interdomain distance was ~12-14 Å (Figure 2A). In the I-1 state, ligand interacted with residues in the α5 helix, α11 helix, α14 helix, α18 helix and 3₁₀H4 helix loop in the two protein subdomains. The ligand RMSD and interdomain distance in the I-1 state were ~9.8 Å and ~18–22 Å, respectively. The receptor adopted a “Partially-open” conformation in this state. In the I-2 state, ligand interacted with receptor subdomain II involving the α17 helix, α18 helix, and α19 helix. The ligand RMSD and interdomain distance was located at (~25.6 Å, ~5-7 Å). In the I-3 state, ligand interactions were observed with the α2 helix and α4 helix in the protein subdomain I and II, respectively, in which the ligand RMSD and interdomain distance was located at (~30.4 Å, ~13-20 Å). In the I-4 state, ligand interacted with the protein subdomain II involving the α8 helix, α14 helix, and 3₁₀H4 with ligand RMSD and the interdomain distance located at (~31.7 Å, ~25-26 Å). In the unbound states U-1, U-2, U-3 and U-4 exhibited interdomain distances of ~5-7, ~10-12, ~20-21, and ~25 Å, respectively, at ligand RMSD ~80 Å.

LiGaMD has successfully captured the dissociation and rebinding pathways of MLN-4760 to the human ACE2 receptor. With detailed analysis of the receptor and protein conformations, LiGaMD simulations suggested that the binding of MLN-4760 to human ACE2 mainly followed the conformational selection mechanism as the ligand binding shifted the receptor conformation to the “Closed” state. These simulation findings provide important insights for design and development of drugs for therapeutic treatment of COVID-19 and other diseases associated with the ACE2 receptor.

In case of SARS-CoV-2 infection, RBD of the viral spike protein attaches to the top region of subdomain I in the ACE2 receptor, which is distant from the receptor peptidase active site. Thus, the MLN-4760 inhibitor does not compete directly with the RBD of SARS-CoV-2. However, there have been evidences of allosteric modulation in the ACE2 receptor, which could be utilized for antiviral drug design. Hernandez Prada et al¹⁰⁰ identified two small-molecule compounds from in silico screening that bound to the hinge region of ACE2 and could activate the receptor ~2 fold in an in vitro ACE2 activity assay, with EC₅₀ values of approximately 20 μM. Huentelman et al¹⁰¹ also discovered a novel ACE2 inhibitor from virtual screening, which bound the receptor active site but also allosterically blocked SARS-CoV spike protein-mediated receptor attachment. Here, we have characterized the mechanism of ligand binding/dissociation and further highlighted the dynamic nature of the ACE2 receptor in terms of the large-scale movement of subdomain I upon ligand binding. LiGaMD simulations also showed that MLN-4760 could form hydrophobic interactions with the RBD binding site in subdomain I of the receptor in the closed state. This information could be helpful in designing drugs that could block RBD binding to the receptor. Furthermore, since the human ACE2 receptor shows conformational selection for ligand binding as revealed from the LiGaMD simulations, virtual screening using ensemble docking^102–104 with receptor structural ensembles generated from the LiGaMD simulations will be a promising approach to designing potent drug molecules of the ACE2 receptor against coronavirus infection.

Mechanism of RNA Binding to M2-1 Protein

As one of non-segmented negative-sense (NNS) RNA viruses, human respiratory syncytial virus (HRSV) including rabies and Ebola virus could cause severe lower track respiratory diseases.^105–108 The function of RNA-dependent RNA polymerase (RdRP) is regulated by a zinc-binding transcription terminator, M2-1 protein, in HRSV. Thus, M2-1 protein is an important therapeutic target for drug design. The structure of a short positive-sense gene-end RNA (SH7) binding to the HRSV M2-1 was resolved, suggesting multiple M2-1 residues were interacted with the RNA.⁸⁵ Five independent 300 ns dual-boost GaMD simulations were performed to investigate the binding mechanism between RNA and the M2-1 protein.⁸⁵ The initial simulation model was prepared by placing the seven-nucleotide SH7 RNA (AGUUAAU) at a distance >20 Å from the M2-1 protein (PDB: 6PZQ). Spontaneous binding of the SH7 RNA to the M2-1 protein was captured. Structural clustering of the SH7 RNA was performed to obtain the top five clusters using a hierarchical agglomerative algorithm. The first and second top-ranked cluster showed RNA binding to the ZBD and the CD of M2-1, respectively. In the ZBD of M2-1, protein residue R4 formed polar contact with the nucleotide U-4 of SH7 RNA (Figure 3B and D). In the CD of M2-1, protein residue K92 formed ionic interaction with the phosphodiester backbone of U-4 RNA nucleotide (Figure 3A and C), being consistent with the crystal structure that a direct interaction of K92 and the phosphodiester backbone of RNA1ʹ was observed.

Figure 3 (A and B) RNA (orange) bound to the Core domain (CD) and RNA (red) bound to the Zn-binding domain (ZBD) of the M2-1 protein (blue ribbons) with highlighted M2-1 interacting residues (magenta sticks) (C) Combined time course distance plots of CD M2-1 residue K92:NZ and RNA base U-4:P (D) Combined time course distance plots of ZBD M2-1 residue R4 and RNA base U-4:P (E) Free energy profiles with radius of gyration (Rg) of RNA and the distance between M2-1 residue K92 and RNA base U4 revealed low-energy Unbound and Bound states of the RNA on M2-1 CD surface (F) Free energy profiles with Rg of RNA and the distance between protein residue R4 and RNA base U4 revealed low-energy Unbound and Bound states of the RNA on M2-1 ZBD. Reprinted from Structure, 28(9), Gao Y, Cao D, Pawnikar S, et al. Structure of the Human Respiratory Syncytial Virus M2-1 Protein in Complex with a Short Positive-Sense Gene-End RNA. 979-990. e974. Copyright © 2020 Elsevier Ltd., with permission from Elsevier.85

Free energy profiles were calculated from GaMD simulations by using the polar interactions and the radius of gyration (Rg) of the SH7 RNA as reaction coordinates. With the RNA molecule being very small, the end-to-end distance can also be used to characterize the molecular conformational changes although it often undergoes high fluctuations. Here, using Rg as a reaction coordinate to calculate free energy profiles allowed us to capture different conformational states of the system. For SH7 binding to the CD of M2-1 protein, the distance between M2-1 residue K92 and RNA base U4 was ~3 Å in the free energy minimum (Figure 3E); While for RNA binding to the ZBD of M2-1, free energy minimum was located at M2-1 residue R4 and RNA base U4 distance ~5 Å (Figure 3F).

In agreement with the experimental results, GaMD simulations revealed that binding of RNA to the CD or ZBD of the M2-1 protein are independent events. Together, the results provided important insights into the recognition mechanism between RNA and HRSV M2-1 protein, which can be valuable for effective drug design and development for its related diseases.

Mechanism of RNA Binding to a Musashi RNA-Binding Proteins

The Musashi 1 (MSI1) protein serves as a therapeutic drug target for treating several cancers such as colorectal, ovarian, bladder and myeloid leukemia. It is known to bind and suppresses translation of the 3ʹ-UTR of Numb mRNA,¹⁰⁹ however, the molecular mechanism of this interaction remains elusive which is important for effective drug design.

We performed all-atom GaMD simulations to study the binding mechanism between Numb RNA and MSI1.⁸⁶ For system setup, the Numb RNA was placed ~30 Å away from the MSI1. The AMBER force fields were used with ff14SBonlysc for protein, RNA.LJbb^110–112 for RNA and TIP3P model for water molecules. Spontaneous binding of Numb RNA to the MSI1 protein was successfully captured in 6 out of the 19 independent 1200 ns of GaMD simulations. In Sim1, RNA binding was observed at ~100 ns where the RMSD of RNA relative to the NMR structure reduced to ~2.50 Å (Figure 4A). In Sim2, RNA binding was observed at ~1010–1130 ns followed by RNA dissociation into the bulk solvent (Figure 4A). The RNA bound to MSI1 after ~800 ns in Sim3, Sim4 and Sim5 (Figure 4A). In Sim6, spontaneous binding of RNA wasobserved after ~1000 ns (Figure 4A).

Figure 4 GaMD simulation of Numb RNA-MSI1 protein interaction. (A) Time courses plots of backbone RMSD of core RNA (central three nucleotides GUA of the Numb RNA) relative to the NMR structure (PDB: 2RS2); (B) Free energy profiles of the backbone RMSD of core RNA relative to the NMR structure and number of native contacts between MSI1 and Numb RNA; (C) Free energy profiles of Rg of the Numb RNA and backbone RMSD of core RNA relative to the NMR structure; (D–F) Low-energy intermediate conformational states I1 (D), I2 (E) and I3 (F) as identified from the free energy profile of the Numb RNA-MSI1 simulation system. Reprinted from Wang J, Lan L, Wu X, Xu L, Miao Y. Mechanism of RNA recognition by a Musashi RNA-binding protein. Curr Res Struct Biol. 2022.4:10-20. Creative Commons license and disclaimer available from: http://creativecommons.org/licenses/by/4.0/legalcode.86

Free energy profiles were calculated by using reaction coordinates including backbone RMSD of RNA core relative to its NMR bound conformation, radius of gyration (Rg) and the number of native contacts (N_contact). Five low-energy minima were characterized - “Bound”, “Intermediate I1”, “Intermediate I2”, “Intermediate I3” and “Unbound”. These states were identified at the backbone RMSD of the RNA core and N_contacts being (2.0 Å, 1500), (5.2 Å, 480), (9.5 Å, 200), (25.0 Å, 10) and (40 Å, 0), respectively (Figure 4B). The “Unbound”, “Intermediate I1”, “Intermediate I2” and “Intermediate I3” states were identified at the backbone RMSD of the RNA core and R_g of Numb being (40.0 Å, 6.2 Å), (5.0 Å, 7.2 Å), (6.9 Å, 6.2 Å), respectively (Figure 4C). From the GaMD simulations, the R_g of the Numb RNA in the “Bound” state was observed to have a wider range as compared to the “Unbound” and “Intermediate” conformations (Figure 4B and C) which suggested an induced fit mechanism of Numb RNA binding to MSI1. The I1 state showed interactions of Numb RNA with the β2-β3 loop and C terminus of MSI1 (Figure 4D). Three hydrogen bonds were formed between MS1 C terminus residue Arg99 and the RNA nucleotide A106 (Figure 4D). The I2 state showed flipping of the MS1 residue Arg61 sidechain towards the solvent leading to the formation of hydrogen bond and salt-bridge interactions with the phosphate oxygen of the sidechain and backbone of RNA nucleotide A106, respectively (Figure 4E). The I3 state showed large conformational changes in the MS1 C terminus where hydrogen bond and salt bridge interactions were observed between the C terminus residue Arg99 and the sidechain and backbone of RNA nucleotide A106, respectively (Figure 4F). These important understandings of the RNA binding mechanism to MSI1 provided by GaMD simulations can aid rational structure-based drug design against MSI1 and other related diseases.

Ligand Binding Thermodynamics and Kinetics Characterized by LiGaMD

LiGaMD simulations performed on model systems such as host-guest and protein-ligand enabled ligand binding thermodynamics and kinetics calculations through repetitive events of binding and dissociation observed of small-molecule ligands.⁸² Calculations were done using LiGaMD and LiGaMD_Dual schemes where a boost potential is applied to the non-bonded and bonded interactions between the ligand-protein/environment, respectively, along with GAFF and q4MD force field parameters used to simulate binding of aspirin and 1-butanol (two guest molecules) to cyclodextrin (CD) (host). PMF profiles were calculated regarding the distance between the guest and host molecules using multiple independent 300 ns LiGaMD simulations. For most of the systems, the bound state was located at the global energy minima. The PMF plot showed similar results as compared to 6500ns cMD simulations in terms of low energy states including “Bound”, “Intermediate” and “Unbound”. Further, to characterize the conformational changes in the host, PMF profiles were calculated using the radius of gyration (R_g) as the reaction coordinate. The apo system and bound guest 1-butanol biased the system to adopt “Compact” (lower R_g) low energy conformational state whereas aspirin bound system predominately adopted “Open” (higher R_g) state. LiGaMD PMF profiles showed that the q4MD force field preferred open conformation where the host predominantly adopted compact conformation using GAFF forcefield. Similarly, the binding free energies of the host-guest systems were calculated using LiGaMD 3D PMF profiles. The binding free energies calculated from hundreds of nanoseconds of LiGaMD simulations were in excellent agreement with that of the experiments and previous “converged” microsecond-timescale cMD simulations. In addition to cMD, LiGaMD simulations’ errors in these values were less than 1.0 kcal/mol. Moreover, LiGaMD simulations were also used to calculate the kinetic binding and dissociation constants. Reweighting of kinetics of the LiGaMD simulations was done using the Transition state theory and Kramer’s rate theory. The kinetic binding rate calculated from the LiGaMD simulations decreased to some degree as compared to the previous cMD simulations while guest dissociation rates increased significantly in LiGaMD_Dual simulations.

Similarly, repetitive binding and dissociation of benzamidine inhibitor in trypsin were observed during microsecond LiGaMD simulations (Figure 5A, B and F). 2D reweighted PMF profile was plotted regarding the atom distances benzamidine:C-Asp189:CG and Trp215:NE-Asp189:CG to characterize conformational changes of protein during ligand binding/unbinding (Figure 5C). Five low-energy wells were characterized from the free energy profiles- Bound, Intermediate 1, Intermediate 2, Unbound 1, and Unbound 2. When the inhibitor was bound to the protein, benzamidine:C-Asp189:CG distance was ~4 Å whereas Trp215:NE-Asp189:CG distance was ~12.5 Å. Here, the inhibitor was bound to the Asp189 residue of the protein S1 pocket as seen in the X-ray crystal conformation. The Intermediate 1 showed inhibitor bound near Trp215 “gate” close to His57-Asp102-Ser214 catalytic whereas the Intermediate 2 showed ligand bound to the open S1* pocket (Figure 5D and E). Rearrangement of Trp215 loop and sidechain flipping of Trp215 residue facilitated the opening of S1* pocket in the Intermediate 2 state. Additionally, the ligand binding free energy of trypsin-benzamidine complex calculated from the 3D PMF profile was −6.13 ± 0.35 kcal/mol which agreed excellently with the experimental value of −6.2 kcal/mol.¹¹¹ Similarly, with the recorded time periods of ligand binding and unbinding, the reweighted k_on and ko_ff values were calculated as 1.15 ± 0.79 × 10⁷ M⁻¹·s⁻¹ and 3.53 ± 1.41 s⁻¹, respectively. Similar protocol was applied, as in the host-guest binding, to calculate these values using Kramer’s rate theory. These data were comparable to the values calculated from experiments.¹¹¹

Figure 5 LiGaMD simulations of benzamidine binding/dissociation to trypsin. (A) Time course plot of distance between benzamidine amidine:C atom and trypsin Asp189:CG atom calculated from representative LiGaMD simulation. (B) Reweighted (LiGaMD_Dual) and non-reweighted (LiGaMD_Dual_Modified) PMF profiles regarding distance between benzamidine amidine:C atom and trypsin Asp189:CG atom. (C) Reweighted 2D PMF profile regarding benzamidine amidine:C – trypsin Asp189:CG and trypsin Trp215:NE – Asp189:CG atom distances. (D and E) Intermediate 1 (I1) and Intermediate 2 (I2) low energy conformational states identified from 2D PMF profiles. (F) Benzamidine (sticks) binding pathway to trypsin (ribbon) obtained from LiGaMD simulations. Catalytic triad (His57, Asp102 and Ser214), Asp189 and Trp215 residues important for trypsin–benzamidine interaction are shown in sticks. Reprinted with permission from Miao Y, Bhattarai A, Wang J. Ligand Gaussian accelerated molecular dynamics (LiGaMD): characterization of ligand binding thermodynamics and kinetics. J Chem Theory Comput. 2020;16(9):5526–5547. Copyright © 2020, American Chemical Society.82

Peptide Binding Thermodynamics and Kinetics Characterized by Pep-GaMD

Pep-GaMD simulations have allowed us to capture repetitive binding and dissociation of highly flexible peptides,⁸³ and thus enabled peptide binding thermodynamics and kinetics calculations. In Pep-GaMD, potential energy of the peptide is selectivity boosted allowing it to sample higher flexibility during the simulation. In the dual boost scheme, the entire system potential energy is boosted. Three model peptide systems including “PAMPAR” (PDB: 1SSH), “PPPALPPKK” (PDB: 1CKA) and “PPPVPPRR” (PDB: 1CKB) were used to capture the binding/dissociation to the SH3 domains.^113,114 1 μs Pep-GaMD allowed us to capture repetitive binding and dissociation of the three model peptides during three independent production runs. In the 1 μs Pep-GaMD simulations, the 1SSH peptide bound 3–4 times and dissociated 3 times to the SH3 domains. The 1CKA peptide was observed having 3–5 binding events whereas 2–5 dissociation events. Similarly, in the 1CKB simulation, the peptide bound 3–4 times to the receptor whereas it dissociated 2–3 times (Figure 6).

Figure 6 Pep-GaMD simulations of binding/dissociation of three model peptides to the SH3 domains: (A–C) X-ray crystal structures of the SH3 domains (green cartoon) bound by three different peptides (magenta cartoon) including (A) “PAMPAR” (PDB: 1SSH), (B) “PPPALPPKK” (PDB: 1CKA) and (C) “PPPVPPRR” (PDB: 1CKB). Key interaction residues are highlighted in sticks. Time courses of (D) 1SSH, (E) 1CKA and (F) 1CKB peptide RMSDs relative to X-ray structures calculated from three independent Pep-GaMD simulations. (G–I) 1D PMF profiles regarding the peptide backbone RMSDs of the (G) 1SSH, (H) 1CKA and (I) 1CKB structures calculated from the Pep-GaMD simulations. Reprinted from Wang J, Miao Y. Peptide Gaussian accelerated molecular dynamics (Pep-GaMD): enhanced sampling and free energy and kinetics calculations of peptide binding. J Chem Phys, 2020;153(15):154109, with the permission of AIP Publishing. © 2020 Author(s). Published under license by AIP Publishing.83

3D Pyreweighting scripts were used to calculate the peptide binding free energies. These values agreed well with the experimental findings. Pep-GaMD predicted the binding free energy value to be −7.72 ± 0.54 kcal/mol for the 1CKA peptide which agreed well with the experimental value of −7.84 kcal/mol.¹¹² The 1CKB peptide had −6.84 ± 0.14 kcal/mol binding free energy which is close to −7.24 kcal/mol¹¹² as calculated from the experiments. The differences of binding free energies between the simulation predictions and experiments were less than 1.0 kcal/mol. In addition, with efficient sampling of peptide binding/dissociation Pep-GaMD helped us to predict the peptide binding kinetics. For the 1CKA system, the k_on and k_off were calculated as 4.06 ± 2.26×10¹⁰ M⁻¹⋅s⁻¹ and 1.45 ± 1.07×10³ s⁻¹, respectively. This was in excellent agreement with the experimental values¹¹⁵ of = 1.5 × 10⁹ M⁻¹⋅s⁻¹ and = 8.9 × 10³ s⁻¹. Hence, the Pep-GaMD simulation has enabled us to study the mechanism of peptide binding/dissociation to proteins by accelerating the interactions between them. This is a significant contribution to the field as it is difficult to attain correct sampling of biomolecules that involves long-range electrostatic interactions and conformational selection.

Conclusions

GaMD presents to be an unconstrained enhanced sampling technique, which can sample the conformational space of biomolecules and their interactions without the need of predefined CVs. Using cumulant expansion to the second order, GaMD can achieve accurate reweighting of the simulations because the added boost potential exhibits a Gaussian distribution. GaMD can be implemented to a broad range of biological systems. In particular, recent GaMD simulations have revealed important binding mechanisms of small-molecule ligands,^71,84 peptides,^116,117 proteins⁷⁸ and nucleic acids^85,86 to the target receptors. 300 ns of GaMD simulations has successfully captured complete folding of chignolin (10 residues) into its native conformation⁶¹ with a minimum RMSD of 0.2 Å relative to the NMR native conformation (PDB: 1UAO). The reweighted free energy profiles calculated showed distinct low-energy conformational states of chignolin as “Folded”, “Intermediate” and “Unfolded”. With the same simulation time, however, aMD failed to distinguish between the unfolded and intermediate states.¹¹⁸ GaMD has been further developed for more powerful techniques of improved sampling such as LiGaMD⁸² and Pep-GaMD.⁸³ Altogether, novel GaMD methods and advancement in computer hardware will enable us to approach and simulate complex biomolecular interactions at atomic level.

Future Perspectives

GaMD has been a successful enhanced sampling method to characterize biomolecular interactions and binding processes. The novel selective GaMD algorithms including LiGaMD and Pep-GaMD can efficiently calculate the binding free energy of ligand/peptide and their kinetic rate constants by adding a selective boost to the essential ligand/peptide potential energy terms. Important mechanistic insights into drug binding to the ACE2 receptor have provided a significant perspective for rational drug development against the COVID-19 pandemic. RNA binding to the M2-1 and MSI1 model proteins depicted by GaMD simulations was in good agreement with the experimentally resolved structures. However, only very few binding events were captured in GaMD simulations of the CXCR4-drug (Figure 1) and protein-RNA (Figures 3 and 4) systems and LiGaMD simulations of the ACE2-inhibitor system (Figure 2). The simulations did not sample enough binding and unbinding events (eg, at least 3–5 or more) in order to estimate the ligand/RNA binding free energies and kinetics for comparison with experimental data. Therefore, these simulations still suffered from insufficient sampling. Further method developments and advances in supercomputing are still needed to obtain converged binding simulations and thus calculate accurate biomolecular binding thermodynamics and kinetics. In addition to binding of the small-molecule ligands and flexible peptides/RNA, new platforms also need to be explored to study biomolecular interactions of larger complexes such as protein–protein interactions. The binding analysis using GaMD and further improved simulation methods shall provide significant contributions in rational computer-aided drug design.