Gaussian Accelerated Molecular Dynamics in NAMD

Overcoming Free-Energy Barriers with a Seamless Combination of a Biasing Force and a Collective Variable-Independent Boost Potential

Amid collective-variable (CV)-based importance-sampling algorithms, a hybrid of the extended adaptive biasing force and the well-tempered metadynamics algorithms (WTM-eABF) has proven particularly cost-effective for exploring the rugged free-energy landscapes that underlie biological processes. However, as an inherently CV-based algorithm, this hybrid scheme does not explicitly accelerate sampling in the space orthogonal to the chosen CVs, thereby limiting its efficiency and accuracy, most notably in those cases where the slow degrees of freedom of the process at hand are not accounted for in the model transition coordinate. Here, inspired by Gaussian-accelerated molecular dynamics (GaMD), we introduce the same CV-independent harmonic boost potential into WTM-eABF, yielding a hybrid algorithm coined GaWTM-eABF. This algorithm leans on WTM-eABF to explore the transition coordinate with a GaMD-mollified potential and recovers the unbiased free-energy landscape through thermodynamic integration followed by proper reweighting. As illustrated in our numerical tests, GaWTM-eABF effectively overcomes the free-energy barriers in orthogonal space and correctly recovers the unbiased potential of mean force (PMF). Furthermore, applying both GaWTM-eABF and WTM-eABF to two biologically relevant processes, namely, the reversible folding of (i) deca-alanine and (ii) chignolin, our results indicate that GaWTM-eABF reduces the uncertainty in the PMF calculation and converges appreciably faster than WTM-eABF. Obviating the need of multiple-copy strategies, GaWTM-eABF is a robust, computationally efficient algorithm to surmount the free-energy barriers in orthogonal space and maps with utmost fidelity the free-energy landscape along selections of CVs. Moreover, our strategy that combines WTM-eABF with GaMD can be easily extended to other biasing-force algorithms.

Structural basis for chemokine recognition and receptor activation of chemokine receptor CCR5

The chemokine receptor CCR5 plays a vital role in immune surveillance and inflammation. However, molecular details that govern its endogenous chemokine recognition and receptor activation remain elusive. Here we report three cryo-electron microscopy structures of G_i1 protein-coupled CCR5 in a ligand-free state and in complex with the chemokine MIP-1α or RANTES, as well as the crystal structure of MIP-1α-bound CCR5. These structures reveal distinct binding modes of the two chemokines and a specific accommodate pattern of the chemokine for the distal N terminus of CCR5. Together with functional data, the structures demonstrate that chemokine-induced rearrangement of toggle switch and plasticity of the receptor extracellular region are critical for receptor activation, while a conserved tryptophan residue in helix II acts as a trigger of receptor constitutive activation.

Identification of pyrogallol as a warhead in design of covalent inhibitors for the SARS-CoV-2 3CL protease

The ongoing pandemic of coronavirus disease 2019 (COVID-19) caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) urgently needs an effective cure. 3CL protease (3CL^pro) is a highly conserved cysteine proteinase that is indispensable for coronavirus replication, providing an attractive target for developing broad-spectrum antiviral drugs. Here we describe the discovery of myricetin, a flavonoid found in many food sources, as a non-peptidomimetic and covalent inhibitor of the SARS-CoV-2 3CL^pro. Crystal structures of the protease bound with myricetin and its derivatives unexpectedly revealed that the pyrogallol group worked as an electrophile to covalently modify the catalytic cysteine. Kinetic and selectivity characterization together with theoretical calculations comprehensively illustrated the covalent binding mechanism of myricetin with the protease and demonstrated that the pyrogallol can serve as an electrophile warhead. Structure-based optimization of myricetin led to the discovery of derivatives with good antiviral activity and the potential of oral administration. These results provide detailed mechanistic insights into the covalent mode of action by pyrogallol-containing natural products and a template for design of non-peptidomimetic covalent inhibitors against 3CL^pros, highlighting the potential of pyrogallol as an alternative warhead in design of targeted covalent ligands.

SARS-CoV-2 3CL protease (3CL^pro) is essential for coronavirus replication and of great interest as an antiviral drug target. Here, the authors show that the naturally occurring flavonoid myricetin is a non-peptidomimetic and covalent inhibitor of 3CL^pro, and they solve crystal structures of 3CL^pro with myricetin and derivatives, which reveal that the pyrogallol group covalently modifies the catalytic cysteine.

Peptide Gaussian accelerated molecular dynamics (Pep-GaMD): Enhanced sampling and free energy and kinetics calculations of peptide binding

Peptides mediate up to 40% of known protein-protein interactions in higher eukaryotes and play an important role in cellular signaling. However, it is challenging to simulate both binding and unbinding of peptides and calculate peptide binding free energies through conventional molecular dynamics, due to long biological timescales and extremely high flexibility of the peptides. Based on the Gaussian accelerated molecular dynamics (GaMD) enhanced sampling technique, we have developed a new computational method "Pep-GaMD," which selectively boosts essential potential energy of the peptide in order to effectively model its high flexibility. In addition, another boost potential is applied to the remaining potential energy of the entire system in a dual-boost algorithm. Pep-GaMD has been demonstrated on binding of three model peptides to the SH3 domains. Independent 1 µs dual-boost Pep-GaMD simulations have captured repetitive peptide dissociation and binding events, which enable us to calculate peptide binding thermodynamics and kinetics. The calculated binding free energies and kinetic rate constants agreed very well with available experimental data. Furthermore, the all-atom Pep-GaMD simulations have provided important insights into the mechanism of peptide binding to proteins that involves long-range electrostatic interactions and mainly conformational selection. In summary, Pep-GaMD provides a highly efficient, easy-to-use approach for unconstrained enhanced sampling and calculations of peptide binding free energies and kinetics.

Pathways and Mechanism of Caffeine Binding to Human Adenosine A2A Receptor

Caffeine (CFF) is a common antagonist to the four subtypes of adenosine G-protein-coupled receptors (GPCRs), which are critical drug targets for treating heart failure, cancer, and neurological diseases. However, the pathways and mechanism of CFF binding to the target receptors remain unclear. In this study, we have performed all-atom-enhanced sampling simulations using a robust Gaussian-accelerated molecular dynamics (GaMD) method to elucidate the binding mechanism of CFF to human adenosine A_2A receptor (A_2AAR). Multiple 500–1,000 ns GaMD simulations captured both binding and dissociation of CFF in the A_2AAR. The GaMD-predicted binding poses of CFF were highly consistent with the x-ray crystal conformations with a characteristic hydrogen bond formed between CFF and residue N6.55 in the receptor. In addition, a low-energy intermediate binding conformation was revealed for CFF at the receptor extracellular mouth between ECL2 and TM1. While the ligand-binding pathways of the A_2AAR were found similar to those of other class A GPCRs identified from previous studies, the ECL2 with high sequence divergence serves as an attractive target site for designing allosteric modulators as selective drugs of the A_2AAR.

Inhibition of Staphylococcus aureus biofilm-forming functional amyloid by molecular tweezers

Biofilms are rigid and largely impenetrable three-dimensional matrices constituting virulence determinants of various pathogenic bacteria. Here, we demonstrate that molecular tweezers, unique supramolecular artificial receptors, modulate biofilm formation of Staphylococcus aureus. In particular, the tweezers affect the structural and assembly properties of phenol-soluble modulin α1 (PSMα1), a biofilm-scaffolding functional amyloid peptide secreted by S. aureus. The data reveal that CLR01, a diphosphate tweezer, exhibits significant S. aureus biofilm inhibition and disrupts PSMα1 self-assembly and fibrillation, likely through inclusion of lysine side chains of the peptide. In comparison, different peptide binding occurs in the case of CLR05, a tweezer containing methylenecarboxylate units, which exhibits lower affinity for the lysine residues yet disrupts S. aureus biofilm more strongly than CLR01. Our study points to a possible role for molecular tweezers as potent biofilm inhibitors and antibacterial agents, particularly against untreatable biofilm-forming and PSM-producing bacteria, such as methicillin-resistant S. aureus.

Perspectives on High-Throughput Ligand/Protein Docking With Martini MD Simulations

Molecular docking is central to rational drug design. Current docking techniques suffer, however, from limitations in protein flexibility and solvation models and by the use of simplified scoring functions. All-atom molecular dynamics simulations, on the other hand, feature a realistic representation of protein flexibility and solvent, but require knowledge of the binding site. Recently we showed that coarse-grained molecular dynamics simulations, based on the most recent version of the Martini force field, can be used to predict protein/ligand binding sites and pathways, without requiring any a priori information, and offer a level of accuracy approaching all-atom simulations. Given the excellent computational efficiency of Martini, this opens the way to high-throughput drug screening based on dynamic docking pipelines. In this opinion article, we sketch the roadmap to achieve this goal.

Specific inhibition of the Survivin-CRM1 interaction by peptide-modified molecular tweezers

Survivin's dual function as apoptosis inhibitor and regulator of cell proliferation is mediated via its interaction with the export receptor CRM1. This protein-protein interaction represents an attractive target in cancer research and therapy. Here, we report a sophisticated strategy addressing Survivin's nuclear export signal (NES), the binding site of CRM1, with advanced supramolecular tweezers for lysine and arginine. These were covalently connected to small peptides resembling the natural, self-complementary dimer interface which largely overlaps with the NES. Several biochemical methods demonstrated sequence-selective NES recognition and interference with the critical receptor interaction. These data were strongly supported by molecular dynamics simulations and multiscale computational studies. Rational design of lysine tweezers equipped with a peptidic recognition element thus allowed to address a previously unapproachable protein surface area. As an experimental proof-of-principle for specific transport signal interference, this concept should be transferable to any protein epitope with a flanking well-accessible lysine.

Disease-related mutations in PI3Kγ disrupt regulatory C-terminal dynamics and reveal a path to selective inhibitors

Class I Phosphoinositide 3-kinases (PI3Ks) are master regulators of cellular functions, with the class IB PI3K catalytic subunit (p110γ) playing key roles in immune signalling. p110γ is a key factor in inflammatory diseases and has been identified as a therapeutic target for cancers due to its immunomodulatory role. Using a combined biochemical/biophysical approach, we have revealed insight into regulation of kinase activity, specifically defining how immunodeficiency and oncogenic mutations of R1021 in the C-terminus can inactivate or activate enzyme activity. Screening of inhibitors using HDX-MS revealed that activation loop-binding inhibitors induce allosteric conformational changes that mimic those in the R1021C mutant. Structural analysis of advanced PI3K inhibitors in clinical development revealed novel binding pockets that can be exploited for further therapeutic development. Overall, this work provides unique insights into regulatory mechanisms that control PI3Kγ kinase activity and shows a framework for the design of PI3K isoform and mutant selective inhibitors.

Elucidation of interactions regulating conformational stability and dynamics of SARS-CoV-2 S-protein

The ongoing COVID-19 pandemic caused by the new coronavirus, SARS-CoV-2, calls for urgent developments of vaccines and antiviral drugs. The spike protein of SARS-CoV-2 (S-protein), which consists of trimeric polypeptide chains with glycosylated residues on the surface, triggers the virus entry into a host cell. Extensive structural and functional studies on this protein have rapidly advanced our understanding of the S-protein structure at atomic resolutions, although most of these structural studies overlook the effect of glycans attached to the S-protein on the conformational stability and functional motions between the inactive down and active up forms. Here, we performed all-atom molecular dynamics simulations of both down and up forms of a fully glycosylated S-protein in solution as well as targeted molecular dynamics simulations between them to elucidate key interdomain interactions for stabilizing each form and inducing the large-scale conformational transitions. The residue-level interaction analysis of the simulation trajectories detects distinct amino acid residues and N-glycans as determinants on conformational stability of each form. During the conformational transitions between them, interdomain interactions mediated by glycosylated residues are switched to play key roles on the stabilization of another form. Electrostatic interactions, as well as hydrogen bonds between the three receptor binding domains, work as driving forces to initiate the conformational transitions toward the active form. This study sheds light on the mechanisms underlying conformational stability and functional motions of the S-protein, which are relevant for vaccine and antiviral drug developments.

Membrane Permeability in Cyclic Peptides is Modulated by Core Conformations

Cyclic peptides have the potential to bind to challenging targets, which are undruggable with small molecules, but their application is limited by low membrane permeability. Here, using a series of cyclic pentapeptides, we showed that established physicochemical criteria of permeable peptides are heavily violated. We revealed that a dominant core conformation, stabilized by amides' shielding pattern, could guide the design of novel compounds. As a result, counter-intuitive strategies, such as incorporation of polar residues, can be beneficial for permeability. We further find that core globularity is a promising descriptor, which can extend the capability of standard predictive models.

Biphasic Force-Regulated Phosphorylation Site Exposure and Unligation of ERM Bound with PSGL-1: A Novel Insight into PSGL-1 Signaling via Steered Molecular Dynamics Simulations

The PSGL-1-actin cytoskeleton linker proteins ezrin/radixin/moesin (ERM), an adaptor between P-selectin glycoprotein ligand-1 (PSGL-1) and spleen tyrosine kinase (Syk), is a key player in PSGL-1 signal, which mediates the adhesion and recruitment of leukocytes to the activated endothelial cells in flow. Binding of PSGL-1 to ERM initials intracellular signaling through inducing phosphorylation of Syk, but effects of tensile force on unligation and phosphorylation site exposure of ERM bound with PSGL-1 remains unclear. To answer this question, we performed a series of so-called “ramp-clamp” steered molecular dynamics (SMD) simulations on the radixin protein FERM domain of ERM bound with intracellular juxtamembrane PSGL-1 peptide. The results showed that, the rupture force of complex pulled with constant velocity was over 250 pN, which prevented the complex from breaking in front of pull-induced exposure of phosphorylation site on immunoreceptor tyrosine activation motif (ITAM)-like motif of ERM; the stretched complex structure under constant tensile forces <100 pN maintained on a stable quasi-equilibrium state, showing a high mechano-stabilization of the clamped complex; and, in consistent with the force-induced allostery at clamped stage, increasing tensile force (<50 pN) would decrease the complex dissociation probability but facilitate the phosphorylation site exposure, suggesting a force-enhanced biophysical connectivity of PSGL-1 signaling. These force-enhanced characters in both phosphorylation and unligation of ERM bound with PSGL-1 should be mediated by a catch-slip bond transition mechanism, in which four residue interactions on binding site were involved. This study might provide a novel insight into the transmembrane PSGL-1 signal, its biophysical connectivity and molecular structural basis for cellular immune responses in mechano-microenvironment, and showed a rational SMD-based computer strategy for predicting structure-function relation of protein under loads.

Agonist Binding and G Protein Coupling in Histamine H2 Receptor: A Molecular Dynamics Study

The histamine H2 receptor (H₂R) plays an important role in the regulation of gastric acid secretion. Therefore, it is a main drug target for the treatment of gastroesophageal reflux or peptic ulcer disease. However, there is as of yet no 3D-structural information available hampering a mechanistic understanding of H₂R. Therefore, we created a model of the histamine-H₂R-G_s complex based on the structure of the ternary complex of the β2-adrenoceptor and investigated the conformational stability of this active GPCR conformation. Since the physiologically relevant motions with respect to ligand binding and conformational changes of GPCRs can only partly be assessed on the timescale of conventional MD (cMD) simulations, we also applied metadynamics and Gaussian accelerated molecular dynamics (GaMD) simulations. A multiple walker metadynamics simulation in combination with cMD was applied for the determination of the histamine binding mode. The preferential binding pose detected is in good agreement with previous data from site directed mutagenesis and provides a basis for rational ligand design. Inspection of the H₂R-G_s interface reveals a network of polar interactions that may contribute to H₂R coupling selectivity. The cMD and GaMD simulations demonstrate that the active conformation is retained on a μs-timescale in the ternary histamine-H₂R-G_s complex and in a truncated complex that contains only G_s helix α5 instead of the entire G protein. In contrast, histamine alone is unable to stabilize the active conformation, which is in line with previous studies of other GPCRs.

Advances in Molecular Dynamics Simulations and Enhanced Sampling Methods for the Study of Protein Systems

Molecular dynamics (MD) simulation is a rigorous theoretical tool that when used efficiently could provide reliable answers to questions pertaining to the structure-function relationship of proteins. Data collated from protein dynamics can be translated into useful statistics that can be exploited to sieve thermodynamics and kinetics crucial for the elucidation of mechanisms responsible for the modulation of biological processes such as protein-ligand binding and protein-protein association. Continuous modernization of simulation tools enables accurate prediction and characterization of the aforementioned mechanisms and these qualities are highly beneficial for the expedition of drug development when effectively applied to structure-based drug design (SBDD). In this review, current all-atom MD simulation methods, with focus on enhanced sampling techniques, utilized to examine protein structure, dynamics, and functions are discussed. This review will pivot around computer calculations of protein-ligand and protein-protein systems with applications to SBDD. In addition, we will also be highlighting limitations faced by current simulation tools as well as the improvements that have been made to ameliorate their efficiency.

Ligand Gaussian Accelerated Molecular Dynamics (LiGaMD): Characterization of Ligand Binding Thermodynamics and Kinetics

Calculations of ligand binding free energies and kinetic rates are important for drug design. However, such tasks have proven challenging in computational chemistry and biophysics. To address this challenge, we have developed a new computational method, ligand Gaussian accelerated molecular dynamics (LiGaMD), which selectively boosts the ligand nonbonded interaction potential energy based on the Gaussian accelerated molecular dynamics (GaMD) enhanced sampling technique. Another boost potential could be applied to the remaining potential energy of the entire system in a dual-boost algorithm (LiGaMD_Dual) to facilitate ligand binding. LiGaMD has been demonstrated on host-guest and protein-ligand binding model systems. Repetitive guest binding and unbinding in the β-cyclodextrin host were observed in hundreds-of-nanosecond LiGaMD_Dual simulations. The calculated guest binding free energies agreed excellently with experimental data with <1.0 kcal/mol errors. Compared with converged microsecond-time scale conventional molecular dynamics simulations, the sampling errors of LiGaMD_Dual simulations were also <1.0 kcal/mol. Accelerations of ligand kinetic rate constants in LiGaMD simulations were properly estimated using Kramers' rate theory. Furthermore, LiGaMD allowed us to capture repetitive dissociation and binding of the benzamidine inhibitor in trypsin within 1 μs simulations. The calculated ligand binding free energy and kinetic rate constants compared well with the experimental data. In summary, LiGaMD provides a powerful enhanced sampling approach for characterizing ligand binding thermodynamics and kinetics simultaneously, which is expected to facilitate computer-aided drug design.

Scalable molecular dynamics on CPU and GPU architectures with NAMD

NAMDis a molecular dynamics program designed for high-performance simulations of very large biological objects on CPU- and GPU-based architectures. NAMD offers scalable performance on petascale parallel supercomputers consisting of hundreds of thousands of cores, as well as on inexpensive commodity clusters commonly found in academic environments. It is written in C++ and leans on Charm++ parallel objects for optimal performance on low-latency architectures. NAMD is a versatile, multipurpose code that gathers state-of-the-art algorithms to carry out simulations in apt thermodynamic ensembles, using the widely popular CHARMM, AMBER, OPLS, and GROMOS biomolecular force fields. Here, we review the main features of NAMD that allow both equilibrium and enhanced-sampling molecular dynamics simulations with numerical efficiency. We describe the underlying concepts utilized by NAMD and their implementation, most notably for handling long-range electrostatics; controlling the temperature, pressure, and pH; applying external potentials on tailored grids; leveraging massively parallel resources in multiple-copy simulations; and hybrid quantum-mechanical/molecular-mechanical descriptions. We detail the variety of options offered by NAMD for enhanced-sampling simulations aimed at determining free-energy differences of either alchemical or geometrical transformations and outline their applicability to specific problems. Last, we discuss the roadmap for the development of NAMD and our current efforts toward achieving optimal performance on GPU-based architectures, for pushing back the limitations that have prevented biologically realistic billion-atom objects to be fruitfully simulated, and for making large-scale simulations less expensive and easier to set up, run, and analyze. NAMD is distributed free of charge with its source code at www.ks.uiuc.edu.

Structure of the Human Respiratory Syncytial Virus M2-1 Protein in Complex with a Short Positive-Sense Gene-End RNA

The M2-1 protein of human respiratory syncytial virus (HRSV) is a transcription anti-terminator that regulates the processivity of the HRSV RNA-dependent RNA polymerase (RdRP). Here, we report a crystal structure of HRSV M2-1 bound to a short positive-sense gene-end RNA (SH7) at 2.7 Å resolution. We identified multiple critical residues of M2-1 involved in RNA interaction and examined their roles using mutagenesis and MicroScale Thermophoresis (MST) assay. We found that hydrophobic residue Phe23 is indispensable for M2-1 to recognize the base of RNA. We also captured spontaneous binding of RNA (SH7) to M2-1 in all-atom simulations using a robust Gaussian accelerated molecular dynamics (GaMD) method. Both experiments and simulations revealed that the interactions of RNA with two separate domains of M2-1, the zinc-binding domain (ZBD) and the core domain (CD), are independent of each other. Collectively, our results provided a structural basis for RNA recognition by HRSV M2-1.

Mechanisms of γ-Secretase Activation and Substrate Processing

Amyloid β-peptide, the principal component of characteristic cerebral plaques of Alzheimer's disease (AD), is produced through intramembrane proteolysis of the amyloid precursor protein (APP) by γ-secretase. Despite the importance in the pathogenesis of AD, the mechanisms of intramembrane proteolysis and substrate processing by γ-secretase remain poorly understood. Here, complementary all-atom simulations using a robust Gaussian accelerated molecular dynamics (GaMD) method and biochemical experiments were combined to investigate substrate processing of wildtype and mutant APP by γ-secretase. The GaMD simulations captured spontaneous activation of γ-secretase, with hydrogen bonded catalytic aspartates and water poised for proteolysis of APP at the ε cleavage site. Furthermore, GaMD simulations revealed that familial AD mutations I45F and T48P enhanced the initial ε cleavage between residues Leu49-Val50, while M51F mutation shifted the ε cleavage site to the amide bond between Thr48-Leu49. Detailed analysis of the GaMD simulations allowed us to identify distinct low-energy conformational states of γ-secretase, different secondary structures of the wildtype and mutant APP substrate, and important active-site subpockets for catalytic function of the enzyme. The simulation findings were highly consistent with experimental analyses of APP proteolytic products using mass spectrometry and Western blotting. Taken together, the GaMD simulations and biochemical experiments have enabled us to elucidate the mechanisms of γ-secretase activation and substrate processing, which should facilitate rational computer-aided drug design targeting this functionally important enzyme.

Probing the Molecular Mechanism of Rifampin Resistance Caused by the Point Mutations S456L and D441V on Mycobacterium tuberculosis RNA Polymerase through Gaussian Accelerated Molecular Dynamics Simulation

Rifampin is the first-line antituberculosis drug, with Mycobacterium tuberculosis RNA polymerase as the molecular target. Unfortunately, M. tuberculosis strains that are resistant to rifampin have been identified in clinical settings, which limits its therapeutic effects. In clinical isolates, S531L and D516V (in Escherichia coli) are two common mutated codons in the gene rpoB, corresponding to S456L and D441V in M. tuberculosis However, the resistance mechanism at the molecular level is still elusive. In this work, Gaussian accelerated molecular dynamics simulations were performed to uncover the resistance mechanism of rifampin due to S456L and D441V mutations at the atomic level. The binding free energy analysis revealed that the reduction in the ability of two mutants to bind rifampin is mainly due to a decrease in electrostatic interaction, specifically, a decrease in the energy contribution of the R454 residue. R454 acts as an anchor and forms stable hydrogen bond interaction with rifampin, allowing rifampin to be stably incorporated in the center of the binding pocket. However, the disappearance of the hydrogen bond between R454 and the mutated residues increases the flexibility of the side chain of R454. The conformation of R454 changes, and the hydrogen bond interaction between it and rifampin is disrupted. As result, the rifampin molecule moves to the outside of the pocket, and the binding affinity decreases. Overall, these findings can provide useful information for understanding the drug resistance mechanism of rifampin and also can give theoretical guidance for further design of novel inhibitors to overcome the drug resistance.

Pathway and mechanism of drug binding to chemokine receptors revealed by accelerated molecular simulations

Background: Chemokine GPCRs play key roles in biology and medicine. Particularly, CXCR4 promotes cancer metastasis and facilitate HIV entry into host cells. Plerixafor (PLX) is a CXCR4 drug, but the pathway and binding site of PLX in CXCR4 remain unknown. Results & methodology: We have performed molecular docking and all-atom simulations using Gaussian accelerated molecular dynamics (GaMD), which are consistent with previous mutation experiments, suggesting that PLX binds to the orthosteric site of CXCR4 as an antagonist. The GaMD simulations further revealed an intermediate allosteric binding site at the extracellular mouth of CXCR4. Conclusion: The newly identified allosteric site can be targeted for novel drug design targeting CXCR4 and other chemokine receptors.

Retrospective ensemble docking of allosteric modulators in an adenosine G-protein-coupled receptor

BACKGROUND

Ensemble docking has proven useful in drug discovery and development. It increases the hit rate by incorporating receptor flexibility into molecular docking as demonstrated on important drug targets including G-protein-coupled receptors (GPCRs). Adenosine A₁ receptor (A₁AR) is a key GPCR that has been targeted for treating cardiac ischemia-reperfusion injuries, neuropathic pain and renal diseases. Development of allosteric modulators, compounds binding to distinct and less conserved GPCR target sites compared with agonists and antagonists, has attracted increasing interest for designing selective drugs of the A₁AR. Despite significant advances, more effective approaches are needed to discover potent and selective allosteric modulators of the A₁AR.

METHODS

Ensemble docking that integrates Gaussian accelerated molecular dynamic (GaMD) simulations and molecular docking using Autodock has been implemented for retrospective docking of known positive allosteric modulators (PAMs) in the A₁AR.

RESULTS

Ensemble docking outperforms docking of the receptor cryo-EM structure. The calculated docking enrichment factors (EFs) and the area under the receiver operating characteristic curves (AUC) are significantly increased.

CONCLUSIONS

Receptor ensembles generated from GaMD simulations are able to increase the success rate of discovering PAMs of A₁AR. It is important to account for receptor flexibility through GaMD simulations and flexible docking.

GENERAL SIGNIFICANCE

Ensemble docking is a promising approach for drug discovery targeting flexible receptors.

Revealing the Positive Binding Cooperativity Mechanism between the Orthosteric and the Allosteric Antagonists of CCR2 by Metadynamics and Gaussian Accelerated Molecular Dynamics Simulations

CC chemokine receptor 2 (CCR2) and its endogenous CC chemokine ligands are associated with numerous inflammatory, neurodegenerative diseases, and cancer. CCR2 is becoming an attractive target in the treatment of autoimmune disease and neurodegenerative diseases. The orthosteric antagonist BMS-681 and allosteric antagonist CCR2-RA-[R] of CCR2 show positive binding cooperativity. We performed well-tempered metadynamics simulations and Gaussian accelerated MD simulations to reveal the influence of the orthosteric antagonist on the unbinding of allosteric antagonist of CCR2. We revealed different unbinding pathways of CCR2-RA-[R] in binary complex CCR2-VT5 and ternary complex CCR2-73R-VT5. The different unbinding pathways of CCR2-RA-[R] are due to the conformational dynamics of TM6. We obtained the significant conformational differences of the intracellular side of TM6 upon CCR2 binding to different ligands by GaMD simulation. The conformational dynamics of TM6 are consistent with the unbinding pathway analysis. GaMD simulations indicate that BMS-681 binding restricts the bend of intracellular side of TM6 by stabilizing the extracellular sides of TM6 and TM7. The charged residues Arg206^5.43 of TM5 and Glu291^7.39 of TM7 play key roles in stabling TM7 and TM6. TM6 and TM7 are crucial components in the orthosteric and allosteric binding sites. Our results illustrate the conformational details about the effect of the orthosteric antagonist on the allosteric antagonist of CCR2. The conformational dynamics of CCR2 upon binding to different ligands can provide a rational basis for development of allosteric ligands of CCR2.

G-Protein-Coupled Receptor-Membrane Interactions Depend on the Receptor Activation State

G-protein-coupled receptors (GPCRs) are the largest family of human membrane proteins and serve as primary targets of approximately one-third of currently marketed drugs. In particular, adenosine A1 receptor (A1 AR) is an important therapeutic target for treating cardiac ischemia-reperfusion injuries, neuropathic pain, and renal diseases. As a prototypical GPCR, the A1 AR is located within a phospholipid membrane bilayer and transmits cellular signals by changing between different conformational states. It is important to elucidate the lipid-protein interactions in order to understand the functional mechanism of GPCRs. Here, all-atom simulations using a robust Gaussian accelerated molecular dynamics (GaMD) method were performed on both the inactive (antagonist bound) and active (agonist and G-protein bound) A1 AR, which was embedded in a 1-palmitoyl-2-oleoyl-glycero-3-phosphocholine (POPC) lipid bilayer. In the GaMD simulations, the membrane lipids played a key role in stabilizing different conformational states of the A1 AR. Our simulations further identified important regions of the receptor that interacted distinctly with the lipids in highly correlated manner. Activation of the A1 AR led to differential dynamics in the upper and lower leaflets of the lipid bilayer. In summary, GaMD enhanced simulations have revealed strongly coupled dynamics of the GPCR and lipids that depend on the receptor activation state. © 2019 Wiley Periodicals, Inc.

Konjac glucomannan/polyvinyl alcohol nanofibers with enhanced skin healing properties by improving fibrinogen adsorption

Skin tissue engineering aims to develop the effective healing strategy to repair the wound by optimizing skin scaffold materials. During the skin wound healing process, fibrin plays an important role due to the specific blood coagulation effect. In this study, the outstanding fibrin capability of konjac glucomannan (KGM) is demonstrated by the molecular dynamics simulation and confirmed by the protein adsorption experiments. A series of konjac glucomannan/polyvinyl alcohol (KGM/PVA) composites with different ratio are fabricated and their role in enhancing the skin repair is tested by in vitro cell culture and in vivo study. The E_ads (adsorption energy) between fibrin and KGM is about 30% larger than that between fibrin and PVA. The fibrinogen adsorption rates of PVA and KGM/PVA (5:5) composites can reach about 20% and 60%, respectively. The results show the blood adsorption capacity of KGM/PVA (5:5) composite can reach about 13 g/g. After 7 days of cell culture, the optical density values of 3T3 fibroblasts on KGM/PVA (5:5) composite could reach 0.8. The mechanical properties of the composites are also verified to meet the practical needs. Thus, we propose a potential wound dressing material strategy based on the materials design and the intrinsic properties of KGM.

Dynamics and Location of the Allosteric Midazolam Site in Cytochrome P4503A4 in Lipid Nanodiscs

Promiscuous and allosteric drug interactions with cytochrome P450 3A4 (CYP3A4) are ubiquitous but incompletely understood at the molecular level. A classic allosteric CYP3A4 drug interaction includes the benzodiazepine midazolam (MDZ). MDZ exhibits homotropic and heterotropic allostery when metabolized to 1'-hydroxy and 4-hydroxy metabolites in varying ratios. The combination of hydrogen-deuterium exchange mass spectrometry (HDX-MS) and Gaussian accelerated molecular dynamics (GaMD) simulations of CYP3A4 in lipid nanodiscs and in a lipid bilayer, respectively, reveals MDZ-dependent changes in dynamics in a membrane environment. The F-, G-, and intervening helices, as well as the loop preceding the β1-sheets, display the largest observed changes in HDX. The GaMD suggests a potential allosteric binding site for MDZ in the F'- and G'-regions, which undergo significant increases in HDX at near-saturating MDZ concentrations. The HDX-MS and GaMD results confirm that changes in dynamics are most significant near the developing consensus allosteric site, and these changes are distinct from those observed previously with the nonallosteric inhibitor ketoconazole. The results suggest that the allosteric MDZ remains mobile in its binding site at the Phe-cluster. The results further suggest that this binding site remains dynamic or changes the depth of insertion in the membrane.

Mechanism of allosteric activation of human mRNA cap methyltransferase (RNMT) by RAM: insights from accelerated molecular dynamics simulations

The RNA guanine-N7 methyltransferase (RNMT) in complex with RNMT-activating miniprotein (RAM) catalyses the formation of a N7-methylated guanosine cap structure on the 5' end of nascent RNA polymerase II transcripts. The mRNA cap protects the primary transcript from exonucleases and recruits cap-binding complexes that mediate RNA processing, export and translation. By using microsecond standard and accelerated molecular dynamics simulations, we provide for the first time a detailed molecular mechanism of allosteric regulation of RNMT by RAM. We show that RAM selects the RNMT active site conformations that are optimal for binding of substrates (AdoMet and the cap), thus enhancing their affinity. Furthermore, our results strongly suggest the likely scenario in which the cap binding promotes the subsequent AdoMet binding, consistent with the previously suggested cooperative binding model. By employing the network community analyses, we revealed the underlying long-range allosteric networks and paths that are crucial for allosteric regulation by RAM. Our findings complement and explain previous experimental data on RNMT activity. Moreover, this study provides the most complete description of the cap and AdoMet binding poses and interactions within the enzyme's active site. This information is critical for the drug discovery efforts that consider RNMT as a promising anti-cancer target.

Improved Modeling of Peptide-Protein Binding Through Global Docking and Accelerated Molecular Dynamics Simulations

Peptides mediate up to 40% of known protein-protein interactions in higher eukaryotes and play a key role in cellular signaling, protein trafficking, immunology, and oncology. However, it is challenging to predict peptide-protein binding with conventional computational modeling approaches, due to slow dynamics and high peptide flexibility. Here, we present a prototype of the approach which combines global peptide docking using ClusPro PeptiDock and all-atom enhanced simulations using Gaussian accelerated molecular dynamics (GaMD). For three distinct model peptides, the lowest backbone root-mean-square deviations (RMSDs) of their bound conformations relative to X-ray structures obtained from PeptiDock were 3.3–4.8 Å, being medium quality predictions according to the Critical Assessment of PRediction of Interactions (CAPRI) criteria. GaMD simulations refined the peptide-protein complex structures with significantly reduced peptide backbone RMSDs of 0.6–2.7 Å, yielding two high quality (sub-angstrom) and one medium quality models. Furthermore, the GaMD simulations identified important low-energy conformational states and revealed the mechanism of peptide binding to the target proteins. Therefore, PeptiDock+GaMD is a promising approach for exploring peptide-protein interactions.

Evolutionary, genetic, structural characterization and its functional implications for the influenza A (H1N1) infection outbreak in India from 2009 to 2017

Influenza A (H1N1) continues to be a major public health threat due to possible emergence of a more virulent H1N1 strain resulting from dynamic changes in virus adaptability consequent to functional mutations and antigenic drift in the hemagglutinin (HA) and neuraminidase (NA) surface proteins. In this study, we describe the genetic and evolutionary characteristics of H1N1 strains that circulated in India over a period of nine years from 2009 to 2017 in relation to global strains. The finding is important from a global perspective since previous phylogenetic studies have suggested that the tropics contributed substantially to the global circulation of influenza viruses. Bayesian phylogenic analysis of HA sequences along with global strains indicated that there is a temporal pattern of H1N1 evolution and clustering of Indian isolates with globally circulating strains. Interestingly, we observed four new amino acid substitutions (S179N, I233T, S181T and I312V) in the HA sequence of H1N1 strains isolated during 2017 and two (S181T and I312V) were found to be unique in Indian isolates. Structurally these two unique mutations could lead to altered glycan specificity of the HA gene. Similarly, sequence and structural analysis of NA domain revealed that the presence of K432E mutation in H1N1 strains isolated after 2015 from India and in global strains found to induce a major loop shift in the vicinity of the catalytic site. The findings presented here offer an insight as to how these acquired mutations could be associated to an improved adaptability of the virus for efficient human transmissibility.

Replica-Exchange Umbrella Sampling Combined with Gaussian Accelerated Molecular Dynamics for Free-Energy Calculation of Biomolecules

We have developed an enhanced conformational sampling method combining replica-exchange umbrella sampling (REUS) with Gaussian accelerated molecular dynamics (GaMD). REUS enhances the sampling along predefined reaction coordinates, while GaMD accelerates the conformational dynamics by adding a boost potential to the system energy. The method, which we call GaREUS (Gaussian accelerated replica-exchange umbrella sampling), enhances the sampling more efficiently than REUS or GaMD, while the computational resource for GaREUS is the same as that required for REUS. The two-step reweighting procedure using the multistate Bennett acceptance ratio method and the cumulant expansion for the exponential average is applied to the simulation trajectories for obtaining the unbiased free-energy landscapes. We apply GaREUS to the calculations of free-energy landscapes for three different cases: conformational equilibria of N-glycan, folding of chignolin, and conformational change of adenyl kinase. We show that GaREUS speeds up the convergences of free-energy calculations using the same amount of computational resources as REUS. The free-energy landscapes reweighted from the trajectories of GaREUS agree with previously reported ones. GaREUS is applicable to free-energy calculations of various biomolecular dynamics and functions with reasonable computational costs.

Broken force dispersal network in tip-links by the mutations at the Ca2+-binding residues induces hearing-loss

Tip-link as force-sensor in hearing conveys the mechanical force originating from sound to ion-channels while maintaining the integrity of the entire sensory assembly in the inner ear. This delicate balance between structure and function of tip-links is regulated by Ca²⁺-ions present in endolymph. Mutations at the Ca²⁺-binding sites of tip-links often lead to congenital deafness, sometimes syndromic defects impairing vision along with hearing. Although such mutations are already identified, it is still not clear how the mutants alter the structure-function properties of the force-sensors associated with diseases. With an aim to decipher the differences in force-conveying properties of the force-sensors in molecular details, we identified the conformational variability of mutant and wild-type tip-links at the single-molecule level using FRET at the endolymphatic Ca²⁺ concentrations and subsequently measured the force-responsive behavior using single-molecule force spectroscopy with an Atomic Force Microscope (AFM). AFM allowed us to mimic the high and wide range of force ramps (10³-10⁶ pN s^-1) as experienced in the inner ear. We performed in silico network analysis to learn that alterations in the conformations of the mutants interrupt the natural force-propagation paths through the sensors and make the mutant tip-links vulnerable to input forces from sound stimuli. We also demonstrated that a Ca²⁺ rich environment can restore the force-response of the mutant tip-links which may eventually facilitate the designing of better therapeutic strategies to the hearing loss.

Mechanistic Insights into Specific G Protein Interactions with Adenosine Receptors

Coupling between G-protein-coupled receptors (GPCRs) and the G proteins is a key step in cellular signaling. Despite extensive experimental and computational studies, the mechanism of specific GPCR-G protein coupling remains poorly understood. This has greatly hindered effective drug design of GPCRs that are primary targets of ∼¹/₃ of currently marketed drugs. Here, we have employed all-atom simulations using a robust Gaussian accelerated molecular dynamics (GaMD) method to decipher the mechanism of the GPCR-G protein interactions. Adenosine receptors (ARs) were used as model systems based on very recently determined cryo-EM structures of the A₁AR and A_2AAR coupled with the G_i and G_s proteins, respectively. Changing the G_i protein to the G_s led to increased fluctuations in the A₁AR and agonist adenosine (ADO), while agonist 5'-N-ethylcarboxamidoadenosine (NECA) binding in the A_2AAR could be still stabilized upon changing the G_s protein to the G_i. Free energy calculations identified one stable low-energy conformation for each of the A₁AR-G_i and A_2AAR-G_s complexes as in the cryo-EM structures, similarly for the A_2AAR-G_i complex. In contrast, the ADO agonist and G_s protein sampled multiple conformations in the A₁AR-G_s system. GaMD simulations thus indicated that the A₁AR preferred to couple with the G_i protein to the G_s, while the A_2AAR could couple with both the G_s and G_i proteins, being highly consistent with experimental findings of the ARs. More importantly, detailed analysis of the atomic simulations showed that the specific AR-G protein coupling resulted from remarkably complementary residue interactions at the protein interface, involving mainly the receptor transmembrane 6 helix and the G_α α5 helix and α4-β6 loop. In summary, the GaMD simulations have provided unprecedented insights into the dynamic mechanism of specific GPCR-G protein interactions at an atomistic level.

How Does Agonist and Antagonist Binding Lead to Different Conformational Ensemble Equilibria of the κ-Opioid Receptor: Insight from Long-Time Gaussian Accelerated Molecular Dynamics Simulation

The opioid receptors belong to the class A seven transmembrane-spanning (7TM) G protein-coupled receptors (GPCRs). The κ-opioid receptor (KOR) is a subfamily of four opioid receptors. The endogenous peptide and a variety of selective agonists and antagonists of KOR have been developed. The structurally similar ligands at the same site cause completely opposite biological functions and induce different conformational changes. To shed light on the conformation ensembles and conformational dynamics in activation and deactivation processes of KOR, we performed all-atom, long-time Gaussian accelerated molecular dynamics simulation (GaMD) on KOR binding with agonist epoxymorphinan MP1104 and antagonist JDTic, respectively. Our results revealed different conformation ensembles of KOR binding with agonist and with antagonist. Agonist binding stabilizes the active state of key motifs including DYYNM motif and CWxP motif, and biases the conformation equilibria toward the active state. Antagonist binding will not destroy inactive conformation equilibria, by keeping the stable inactive state of these crucial motifs. We found that the inactive apo form of KOR is the most stable state, while the active apo form relaxes readily to inactive state. Our results also revealed a stable intermediate (I), which is attributed to the hydrophobic interactions between Tyr246^5.58 and TM6, as well as the steric hindrance of them. Our results not only show the conformation equilibria bias of KOR by binding with agonist and antagonist, but also provide the structural information for the design and discovery of potential ligands with different functions.

Recent advances in computational studies of GPCR-G protein interactions

Protein-protein interactions are key in cellular signaling. G protein-coupled receptors (GPCRs), the largest superfamily of human membrane proteins, are able to transduce extracellular signals (e.g., hormones and neurotransmitters) to intracellular proteins, in particular the G proteins. Since GPCRs serve as primary targets of ~1/3 of currently marketed drugs, it is important to understand mechanisms of GPCR signaling in order to design selective and potent drug molecules. This chapter focuses on recent advances in computational studies of the GPCR-G protein interactions using bioinformatics, protein-protein docking and molecular dynamics simulation approaches.

Structural Basis for Binding of Allosteric Drug Leads in the Adenosine A1 Receptor

Despite intense interest in designing positive allosteric modulators (PAMs) as selective drugs of the adenosine A₁ receptor (A₁AR), structural binding modes of the receptor PAMs remain unknown. Using the first X-ray structure of the A₁AR, we have performed all-atom simulations using a robust Gaussian accelerated molecular dynamics (GaMD) technique to determine binding modes of the A₁AR allosteric drug leads. Two prototypical PAMs, PD81723 and VCP171, were selected. Each PAM was initially placed at least 20 Å away from the receptor. Extensive GaMD simulations using the AMBER and NAMD simulation packages at different acceleration levels captured spontaneous binding of PAMs to the A₁AR. The simulations allowed us to identify low-energy binding modes of the PAMs at an allosteric site formed by the receptor extracellular loop 2 (ECL2), which are highly consistent with mutagenesis experimental data. Furthermore, the PAMs stabilized agonist binding in the receptor. In the absence of PAMs at the ECL2 allosteric site, the agonist sampled a significantly larger conformational space and even dissociated from the A₁AR alone. In summary, the GaMD simulations elucidated structural binding modes of the PAMs and provided important insights into allostery in the A₁AR, which will greatly facilitate the receptor structure-based drug design.

Gaussian accelerated molecular dynamics for elucidation of drug pathways

INTRODUCTION

Understanding pathways and mechanisms of drug binding to receptors is important for rational drug design. Remarkable advances in supercomputing and methodological developments have opened a new era for application of computer simulations in predicting drug-receptor interactions at an atomistic level. Gaussian accelerated molecular dynamics (GaMD) is a computational enhanced sampling technique that works by adding a harmonic boost potential to reduce energy barriers. GaMD enables free energy calculations without the requirement of predefined collective variables. GaMD has proven useful in biomolecular simulations, in particular, the prediction of drug-receptor interactions. Areas covered: Herein, the authors review recent GaMD simulation studies that elucidated pathways of drug binding to proteins including the G-protein-coupled receptors and HIV protease. Expert opinion: GaMD is advantageous for enhanced simulations of, amongst many biological processes, drug binding to target receptors. Compared with conventional molecular dynamics, GaMD speeds up biomolecular simulations by orders of magnitude. GaMD enables routine drug binding simulations using personal computers with GPUs or common computing clusters. GaMD and, more broadly, enhanced sampling simulations are expected to dramatically increase our capabilities to determine the mechanisms of drug binding to a wide range of receptors in the near future. This will greatly facilitate computer-aided drug design.

Identification of SLAC1 anion channel residues required for CO2/bicarbonate sensing and regulation of stomatal movements

Increases in CO2 concentration in plant leaves due to respiration in the dark and the continuing atmospheric [CO2] rise cause closing of stomatal pores, thus affecting plant-water relations globally. However, the underlying CO2/bicarbonate (CO2/HCO3-) sensing mechanisms remain unknown. [CO2] elevation in leaves triggers stomatal closure by anion efflux mediated via the SLAC1 anion channel localized in the plasma membrane of guard cells. Previous reconstitution analysis has suggested that intracellular bicarbonate ions might directly up-regulate SLAC1 channel activity. However, whether such a CO2/HCO3- regulation of SLAC1 is relevant for CO2 control of stomatal movements in planta remains unknown. Here, we computationally probe for candidate bicarbonate-interacting sites within the SLAC1 anion channel via long-timescale Gaussian accelerated molecular dynamics (GaMD) simulations. Mutations of two putative bicarbonate-interacting residues, R256 and R321, impaired the enhancement of the SLAC1 anion channel activity by CO2/HCO3- in Xenopus oocytes. Mutations of the neighboring charged amino acid K255 and residue R432 and the predicted gate residue F450 did not affect HCO3- regulation of SLAC1. Notably, gas-exchange experiments with slac1-transformed plants expressing mutated SLAC1 proteins revealed that the SLAC1 residue R256 is required for CO2 regulation of stomatal movements in planta, but not for abscisic acid (ABA)-induced stomatal closing. Patch clamp analyses of guard cells show that activation of S-type anion channels by CO2/HCO3-, but not by ABA, was impaired, indicating the relevance of R256 for CO2 signal transduction. Together, these analyses suggest that the SLAC1 anion channel is one of the physiologically relevant CO2/HCO3- sensors in guard cells.

Acceleration of biomolecular kinetics in Gaussian accelerated molecular dynamics

Recent studies demonstrated that Gaussian accelerated molecular dynamics (GaMD) is a robust computational technique, which provides simultaneous unconstrained enhanced sampling and free energy calculations of biomolecules. However, the exact acceleration of biomolecular dynamics or speedup of kinetic rates in GaMD simulations and, more broadly, in enhanced sampling methods, remains a challenging task to be determined. Here, the GaMD acceleration is examined using alanine dipeptide in explicit solvent as a biomolecular model system. Relative to long conventional molecular dynamics simulation, GaMD simulations exhibited ∼36-67 times speedup for sampling of the backbone dihedral transitions. The acceleration depended on level of the GaMD boost potential. Furthermore, Kramers' rate theory was applied to estimate GaMD acceleration using simulation-derived diffusion coefficients, curvatures and barriers of free energy profiles. In most cases, the calculations also showed significant speedup of dihedral transitions in GaMD, although the GaMD acceleration factors tended to be underestimated by ∼3-96 fold. Because greater boost potential can be applied in GaMD simulations of systems with increased sizes, which potentially leads to higher acceleration, it is subject to future studies on accelerating the dynamics and recovering kinetic rates of larger biomolecules such as proteins and protein-protein/nucleic acid complexes.

Isolation of a Structural Mechanism for Uncoupling T Cell Receptor Signaling from Peptide-MHC Binding

TCR-signaling strength generally correlates with peptide-MHC binding affinity; however, exceptions exist. We find high-affinity, yet non-stimulatory, interactions occur with high frequency in the human T cell repertoire. Here, we studied human TCRs that are refractory to activation by pMHC ligands despite robust binding. Analysis of 3D affinity, 2D dwell time, and crystal structures of stimulatory versus non-stimulatory TCR-pMHC interactions failed to account for their different signaling outcomes. Using yeast pMHC display, we identified peptide agonists of a formerly non-responsive TCR. Single-molecule force measurements demonstrated the emergence of catch bonds in the activating TCR-pMHC interactions, correlating with exclusion of CD45 from the TCR-APC contact site. Molecular dynamics simulations of TCR-pMHC disengagement distinguished agonist from non-agonist ligands based on the acquisition of catch bonds within the TCR-pMHC interface. The isolation of catch bonds as a parameter mediating the coupling of TCR binding and signaling has important implications for TCR and antigen engineering for immunotherapy.

Ligand Binding Pathways and Conformational Transitions of the HIV Protease

It is important to determine the binding pathways and mechanisms of ligand molecules to target proteins to effectively design therapeutic drugs. Molecular dynamics (MD) is a promising computational tool that allows us to simulate protein-drug binding at an atomistic level. However, the gap between the time scales of current simulations and those of many drug binding processes has limited the usage of conventional MD, which has been reflected in studies of the HIV protease. Here, we have applied a robust enhanced simulation method, Gaussian accelerated molecular dynamics (GaMD), to sample binding pathways of the XK263 ligand and associated protein conformational changes in the HIV protease. During two of 10 independent GaMD simulations performed over 500-2500 ns, the ligand was observed to successfully bind to the protein active site. Although GaMD-derived free energy profiles were not fully converged because of insufficient sampling of the complex system, the simulations still allowed us to identify relatively low-energy intermediate conformational states during binding of the ligand to the HIV protease. Relative to the X-ray crystal structure, the XK263 ligand reached a minimum root-mean-square deviation (RMSD) of 2.26 Å during 2.5 μs of GaMD simulation. In comparison, the ligand RMSD reached a minimum of only ∼5.73 Å during an earlier 14 μs conventional MD simulation. This work highlights the enhanced sampling power of the GaMD approach and demonstrates its wide applicability to studies of drug-receptor interactions for the HIV protease and by extension many other target proteins.

Mechanism of the G-protein mimetic nanobody binding to a muscarinic G-protein-coupled receptor

Protein-protein binding is key in cellular signaling processes. Molecular dynamics (MD) simulations of protein-protein binding, however, are challenging due to limited timescales. In particular, binding of the medically important G-protein-coupled receptors (GPCRs) with intracellular signaling proteins has not been simulated with MD to date. Here, we report a successful simulation of the binding of a G-protein mimetic nanobody to the M₂ muscarinic GPCR using the robust Gaussian accelerated MD (GaMD) method. Through long-timescale GaMD simulations over 4,500 ns, the nanobody was observed to bind the receptor intracellular G-protein-coupling site, with a minimum rmsd of 2.48 Å in the nanobody core domain compared with the X-ray structure. Binding of the nanobody allosterically closed the orthosteric ligand-binding pocket, being consistent with the recent experimental finding. In the absence of nanobody binding, the receptor orthosteric pocket sampled open and fully open conformations. The GaMD simulations revealed two low-energy intermediate states during nanobody binding to the M₂ receptor. The flexible receptor intracellular loops contribute remarkable electrostatic, polar, and hydrophobic residue interactions in recognition and binding of the nanobody. These simulations provided important insights into the mechanism of GPCR-nanobody binding and demonstrated the applicability of GaMD in modeling dynamic protein-protein interactions.

Gaussian Accelerated Molecular Dynamics: Theory, Implementation, and Applications

A novel Gaussian Accelerated Molecular Dynamics (GaMD) method has been developed for simultaneous unconstrained enhanced sampling and free energy calculation of biomolecules. Without the need to set predefined reaction coordinates, GaMD enables unconstrained enhanced sampling of the biomolecules. Furthermore, by constructing a boost potential that follows a Gaussian distribution, accurate reweighting of GaMD simulations is achieved via cumulant expansion to the second order. The free energy profiles obtained from GaMD simulations allow us to identify distinct low energy states of the biomolecules and characterize biomolecular structural dynamics quantitatively. In this chapter, we present the theory of GaMD, its implementation in the widely used molecular dynamics software packages (AMBER and NAMD), and applications to the alanine dipeptide biomolecular model system, protein folding, biomolecular large-scale conformational transitions and biomolecular recognition.