Understanding the impact of binding free energy and kinetics calculations in modern drug discovery.

INTRODUCTION: For rational drug design, it is crucial to understand the receptor-drug binding processes and mechanisms. A new era for the use of computer simulations in predicting drug-receptor interactions at an atomic level has begun with remarkable advances in supercomputing and methodological breakthroughs. AREAS COVERED: End-point free energy calculation methods such as Molecular Mechanics/Poisson Boltzmann Surface Area (MM/PBSA) or Molecular-Mechanics/Generalized Born Surface Area (MM/GBSA), free energy perturbation (FEP), and thermodynamic integration (TI) are commonly used for binding free energy calculations in drug discovery. In addition, kinetic dissociation and association rate constants (koff and kon) play critical roles in the function of drugs. Nowadays, Molecular Dynamics (MD) and enhanced sampling simulations are increasingly being used in drug discovery. Here, the authors provide a review of the computational techniques used in drug binding free energy and kinetics calculations. EXPERT OPINION: The applications of computational methods in drug discovery and design are expanding, thanks to improved predictions of the binding free energy and kinetic rates of drug molecules. Recent microsecond-timescale enhanced sampling simulations have made it possible to accurately capture repetitive ligand binding and dissociation, facilitating more efficient and accurate calculations of ligand binding free energy and kinetics.

Structure and function of an intermediate GPCR-Gαßγ complex.

Unraveling the signaling roles of intermediate complexes is pivotal for G protein-coupled receptor (GPCR) drug development. Despite hundreds of GPCR-Gαßγ structures, these snapshots primarily capture the fully activated end-state complex. Consequently, a comprehensive understanding of the conformational transitions during GPCR activation and the roles of intermediate GPCR-G protein complexes in signaling remain elusive. Guided by a conformational landscape profiled by 19 F quantitative NMR ( 19 F-qNMR) and Molecular Dynamics (MD) simulations, we resolved the structure of an unliganded GPCR-G protein intermediate complex by blocking its transition to the fully activated end-state complex. More importantly, we presented direct evidence that the intermediate GPCR-Gαsßγ complex initiates a rate-limited nucleotide exchange without progressing to the fully activated end-state complex, thereby bridging a significant gap in our understanding the complexity of GPCR signaling. Understanding the roles of individual conformational states and their complexes in signaling efficacy and bias will help us to design drugs that discriminately target a disease-related conformation.

Characterization of natural product inhibitors of quorum sensing in Pseudomonas aeruginosa reveals competitive inhibition of RhlR by ortho-vanillin.

Quorum sensing (QS) is a cell-cell signaling system that enables bacteria to coordinate population density-dependent changes in behavior. This chemical communication pathway is mediated by diffusible N-acyl L-homoserine lactone signals and cytoplasmic signal-responsive LuxR-type receptors in Gram-negative bacteria. As many common pathogenic bacteria use QS to regulate virulence, there is significant interest in disrupting QS as a potential therapeutic strategy. Prior studies have implicated the natural products salicylic acid, cinnamaldehyde and other related benzaldehyde derivatives as inhibitors of QS in the opportunistic pathogen Pseudomonas aeruginosa, yet we lack an understanding of the mechanisms by which these compounds function. Herein, we evaluate the activity of a set of benzaldehyde derivatives using heterologous reporters of the P. aeruginosa LasR and RhlR QS signal receptors. We find that most tested benzaldehyde derivatives can antagonize LasR or RhlR reporter activation at micromolar concentrations, although certain molecules also caused mild growth defects and nonspecific reporter antagonism. Notably, several compounds showed promising RhlR or LasR specific inhibitory activities over a range of concentrations below that causing toxicity. Ortho-Vanillin, a previously untested compound, was the most promising within this set. Competition experiments against the native ligands for LasR and RhlR revealed that ortho-vanillin can interact competitively with RhlR but not with LasR. Overall, these studies expand our understanding of benzaldehyde activities in the LasR and RhlR receptors and reveal potentially promising effects of ortho-vanillin as a small molecule QS modulator against RhlR.

Mechanistic insights into ligand dissociation from the SARS-CoV-2 spike glycoprotein.

The COVID-19 pandemic, driven by the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), has spurred an urgent need for effective therapeutic interventions. The spike glycoprotein of the SARS-CoV-2 is crucial for infiltrating host cells, rendering it a key candidate for drug development. By interacting with the human angiotensin-converting enzyme 2 (ACE2) receptor, the spike initiates the infection of SARS-CoV-2. Linoleate is known to bind the spike glycoprotein, subsequently reducing its interaction with ACE2. However, the detailed mechanisms underlying the protein-ligand interaction remain unclear. In this study, we characterized the pathways of ligand dissociation and the conformational changes associated with the spike glycoprotein by using ligand Gaussian accelerated molecular dynamics (LiGaMD). Our simulations resulted in eight complete ligand dissociation trajectories, unveiling two distinct ligand unbinding pathways. The preference between these two pathways depends on the gate distance between two α-helices in the receptor binding domain (RBD) and the position of the N-linked glycan at N343. Our study also highlights the essential contributions of K417, N121 glycan, and N165 glycan in ligand unbinding, which are equally crucial in enhancing spike-ACE2 binding. We suggest that the presence of the ligand influences the motions of these residues and glycans, consequently reducing accessibility for spike-ACE2 binding. These findings enhance our understanding of ligand dissociation from the spike glycoprotein and offer significant implications for drug design strategies in the battle against COVID-19.

Familial Alzheimer mutations stabilize synaptotoxic γ-secretase-substrate complexes.

Mutations that cause familial Alzheimer's disease (FAD) are found in amyloid precursor protein (APP) and presenilin, the catalytic component of γ-secretase, that together produce amyloid ß-peptide (Aß). Nevertheless, whether Aß is the primary disease driver remains controversial. We report here that FAD mutations disrupt initial proteolytic events in the multistep processing of APP substrate C99 by γ-secretase. Cryoelectron microscopy reveals that a substrate mimetic traps γ-secretase during the transition state, and this structure aligns with activated enzyme-substrate complex captured by molecular dynamics simulations. In silico simulations and in cellulo fluorescence microscopy support stabilization of enzyme-substrate complexes by FAD mutations. Neuronal expression of C99 and/or presenilin-1 in Caenorhabditis elegans leads to synaptic loss only with FAD-mutant transgenes. Designed mutations that stabilize the enzyme-substrate complex and block Aß production likewise led to synaptic loss. Collectively, these findings implicate the stalled process-not the products-of γ-secretase cleavage of substrates in FAD pathogenesis.

Mechanism of Ligand Binding to Theophylline RNA Aptamer.

Studying RNA-ligand interactions and quantifying their binding thermodynamics and kinetics are of particular relevance in the field of drug discovery. Here, we combined biochemical binding assays and accelerated molecular simulations to investigate ligand binding and dissociation in RNA using the theophylline-binding RNA as a model system. All-atom simulations using a Ligand Gaussian accelerated Molecular Dynamics method (LiGaMD) have captured repetitive binding and dissociation of theophylline and caffeine to RNA. Theophylline's binding free energy and kinetic rate constants align with our experimental data, while caffeine's binding affinity is over 10,000 times weaker, and its kinetics could not be determined. LiGaMD simulations allowed us to identify distinct low-energy conformations and multiple ligand binding pathways to RNA. Simulations revealed a "conformational selection" mechanism for ligand binding to the flexible RNA aptamer, which provides important mechanistic insights into ligand binding to the theophylline-binding model. Our findings suggest that compound docking using a structural ensemble of representative RNA conformations would be necessary for structure-based drug design of flexible RNA.

Activation of Polycystin-1 Signaling by Binding of Stalk-derived Peptide Agonists.

Polycystin-1 (PC1) is the membrane protein product of the PKD1 gene whose mutation is responsible for 85% of the cases of autosomal dominant polycystic kidney disease (ADPKD). ADPKD is primarily characterized by the formation of renal cysts and potential kidney failure. PC1 is an atypical G protein-coupled receptor (GPCR) consisting of 11 transmembrane helices and an autocatalytic GAIN domain that cleaves PC1 into extracellular N-terminal (NTF) and membrane-embedded C-terminal (CTF) fragments. Recently, signaling activation of the PC1 CTF was shown to be regulated by a stalk tethered agonist (TA), a distinct mechanism observed in the adhesion GPCR family. A novel allosteric activation pathway was elucidated for the PC1 CTF through a combination of Gaussian accelerated molecular dynamics (GaMD), mutagenesis and cellular signaling experiments. Here, we show that synthetic, soluble peptides with 7 to 21 residues derived from the stalk TA, in particular, peptides including the first 9 residues (p9), 17 residues (p17) and 21 residues (p21) exhibited the ability to re-activate signaling by a stalkless PC1 CTF mutant in cellular assays. To reveal molecular mechanisms of stalk peptide-mediated signaling activation, we have applied a novel Peptide GaMD (Pep-GaMD) algorithm to elucidate binding conformations of selected stalk peptide agonists p9, p17 and p21 to the stalkless PC1 CTF. The simulations revealed multiple specific binding regions of the stalk peptide agonists to the PC1 protein including an "intermediate" bound yet inactive state. Our Pep-GaMD simulation findings were consistent with the cellular assay experimental data. Binding of peptide agonists to the TOP domain of PC1 induced close TOP-putative pore loop interactions, a characteristic feature of the PC1 CTF signaling activation mechanism. Using sequence covariation analysis of PC1 homologs, we further showed that the peptide binding regions were consistent with covarying residue pairs identified between the TOP domain and the stalk TA. Therefore, structural dynamic insights into the mechanisms of PC1 activation by stalk-derived peptide agonists have enabled an in-depth understanding of PC1 signaling. They will form a foundation for development of PC1 as a therapeutic target for the treatment of ADPKD.

Deep Learning Dynamic Allostery of G-Protein-Coupled Receptors.

G-protein-coupled receptors (GPCRs) make up the largest superfamily of human membrane proteins and represent primary targets of ∼1/3 of currently marketed drugs. Allosteric modulators have emerged as more selective drug candidates compared with orthosteric agonists and antagonists. However, many X-ray and cryo-EM structures of GPCRs resolved so far exhibit negligible differences upon the binding of positive and negative allosteric modulators (PAMs and NAMs). The mechanism of dynamic allosteric modulation in GPCRs remains unclear. In this work, we have systematically mapped dynamic changes in free energy landscapes of GPCRs upon binding of allosteric modulators using the Gaussian accelerated molecular dynamics (GaMD), deep learning (DL), and free energy prOfiling Workflow (GLOW). GaMD simulations were performed for a total of 66 µs on 44 GPCR systems in the presence and absence of the modulator. DL and free energy calculations revealed significantly reduced dynamic fluctuations and conformational space of GPCRs upon modulator binding. While the modulator-free GPCRs often sampled multiple low-energy conformational states, the NAMs and PAMs confined the inactive and active agonist-G-protein-bound GPCRs, respectively, to mostly only one specific conformation for signaling. Such cooperative effects were significantly reduced for binding of the selective modulators to "non-cognate" receptor subtypes. Therefore, GPCR allostery exhibits a dynamic "conformational selection" mechanism. In the absence of available modulator-bound structures as for most current GPCRs, it is critical to use a structural ensemble of representative GPCR conformations rather than a single structure for compound docking ("ensemble docking"), which will potentially improve structure-based design of novel allosteric drugs of GPCRs.

Molecular Dynamics Activation of γ-Secretase for Cleavage of the Notch1 Substrate.

γ-Secretase is an intramembrane aspartyl protease complex that cleaves the transmembrane domain of over 150 peptide substrates, including amyloid precursor protein (APP) and the Notch family of receptors, via two conserved aspartates D257 and D385 in the presenilin-1 (PS1) catalytic subunit. However, while the activation of γ-secretase for cleavage of APP has been widely studied, the cleavage of Notch by γ-secretase remains poorly explored. Here, we combined Gaussian accelerated molecular dynamics (GaMD) simulations and mass spectrometry (MS) analysis of proteolytic products to present the first dynamic models for cleavage of Notch by γ-secretase. MS showed that γ-secretase cleaved the WT Notch at Notch residue G34, while cleavage of the L36F mutant Notch occurred at Notch residue C33. Initially, we prepared our simulation systems starting from the cryoEM structure of Notch-bound γ-secretase (PDB: 6IDF) and failed to capture the proper cleavages of WT and L36F Notch by γ-secretase. We then discovered an incorrect registry of the Notch substrate in the PS1 active site through alignment of the experimental structure of Notch-bound (PDB: 6IDF) and APP-bound γ-secretase (PDB: 6IYC). Every residue of the APP substrate was systematically mutated to the corresponding Notch residue to prepare a resolved model of Notch-bound γ-secretase complexes. GaMD simulations of the resolved model successfully captured γ-secretase activation for proper cleavages of both WT and L36F mutant Notch. Our findings presented here provided mechanistic insights into the structural dynamics and enzyme-substrate interactions required for γ-secretase activation for cleavage of Notch and other substrates.

Structural dynamics of chemokine receptors.

Membrane proteins such as G protein-coupled receptors (GPCRs) are involved in awide range of physiological and pathological cellular processes. Binding of extracellular signals to GPCRs, including hormones, neurotransmitters, peptides and proteins, can activate intracellular signaling cascades via G protein interaction. Chemokine receptors are key GPCRs implicated in cancers, immune responses, cell migration and inflammation. Specifically, the CCR5 and CXCR4 chemokine receptors serve as important therapeutic targets against Human Immunodeficiency virus (HIV) entry into human cells. Maraviroc and Vicriviroc, two clinically used HIV entry inhibitors, are antagonists of the CCR5 receptor. These drugs block HIV entry, but ultimately resistance develops, due to emergence of viruses that can utilize the CXCR4 co-receptor. Unfortunately, development of chemokine receptor antagonists as selective drugs of HIV infection has been greatly hindered as their target orthosteric site is conserved among different receptor subtypes. Accordingly, it is important to understand the structural dynamics of these receptors to develop more effective therapeutics. In this chapter, we describe the latest advances in studies of these two key chemokine receptors with respect to their structures, dynamics and function.

Structural basis of agonist specificity of α1A-adrenergic receptor.

α1-adrenergic receptors (α1-ARs) play critical roles in the cardiovascular and nervous systems where they regulate blood pressure, cognition, and metabolism. However, the lack of specific agonists for all α1 subtypes has limited our understanding of the physiological roles of different α1-AR subtypes, and led to the stagnancy in agonist-based drug development for these receptors. Here we report cryo-EM structures of α1A-AR in complex with heterotrimeric G-proteins and either the endogenous common agonist epinephrine or the α1A-AR-specific synthetic agonist A61603. These structures provide molecular insights into the mechanisms underlying the discrimination between α1A-AR and α1B-AR by A61603. Guided by the structures and corresponding molecular dynamics simulations, we engineer α1A-AR mutants that are not responsive to A61603, and α1B-AR mutants that can be potently activated by A61603. Together, these findings advance our understanding of the agonist specificity for α1-ARs at the molecular level, opening the possibility of rational design of subtype-specific agonists.

Pharmacological hallmarks of allostery at the M4 muscarinic receptor elucidated through structure and dynamics.

Allosteric modulation of G protein-coupled receptors (GPCRs) is a major paradigm in drug discovery. Despite decades of research, a molecular-level understanding of the general principles that govern the myriad pharmacological effects exerted by GPCR allosteric modulators remains limited. The M4 muscarinic acetylcholine receptor (M4 mAChR) is a validated and clinically relevant allosteric drug target for several major psychiatric and cognitive disorders. In this study, we rigorously quantified the affinity, efficacy, and magnitude of modulation of two different positive allosteric modulators, LY2033298 (LY298) and VU0467154 (VU154), combined with the endogenous agonist acetylcholine (ACh) or the high-affinity agonist iperoxo (Ipx), at the human M4 mAChR. By determining the cryo-electron microscopy structures of the M4 mAChR, bound to a cognate Gi1 protein and in complex with ACh, Ipx, LY298-Ipx, and VU154-Ipx, and applying molecular dynamics simulations, we determine key molecular mechanisms underlying allosteric pharmacology. In addition to delineating the contribution of spatially distinct binding sites on observed pharmacology, our findings also revealed a vital role for orthosteric and allosteric ligand-receptor-transducer complex stability, mediated by conformational dynamics between these sites, in the ultimate determination of affinity, efficacy, cooperativity, probe dependence, and species variability. There results provide a holistic framework for further GPCR mechanistic studies and can aid in the discovery and design of future allosteric drugs.

Deep Boosted Molecular Dynamics: Accelerating Molecular Simulations with Gaussian Boost Potentials Generated Using Probabilistic Bayesian Deep Neural Network.

We have developed a new deep boosted molecular dynamics (DBMD) method. Probabilistic Bayesian neural network models were implemented to construct boost potentials that exhibit Gaussian distribution with minimized anharmonicity, thereby allowing for accurate energetic reweighting and enhanced sampling of molecular simulations. DBMD was demonstrated on model systems of alanine dipeptide and the fast-folding protein and RNA structures. For alanine dipeptide, 30 ns DBMD simulations captured up to 83-125 times more backbone dihedral transitions than 1 µs conventional molecular dynamics (cMD) simulations and were able to accurately reproduce the original free energy profiles. Moreover, DBMD sampled multiple folding and unfolding events within 300 ns simulations of the chignolin model protein and identified low-energy conformational states comparable to previous simulation findings. Finally, DBMD captured a general folding pathway of three hairpin RNAs with the GCAA, GAAA, and UUCG tetraloops. Based on a deep learning neural network, DBMD provides a powerful and generally applicable approach to boosting biomolecular simulations. DBMD is available with open source in OpenMM at https://github.com/MiaoLab20/DBMD/.

Deep Boosted Molecular Dynamics (DBMD): Accelerating molecular simulations with Gaussian boost potentials generated using probabilistic Bayesian deep neural network.

We have developed a new Deep Boosted Molecular Dynamics (DBMD) method. Probabilistic Bayesian neural network models were implemented to construct boost potentials that exhibit Gaussian distribution with minimized anharmonicity, thereby allowing for accurate energetic reweighting and enhanced sampling of molecular simulations. DBMD was demonstrated on model systems of alanine dipeptide and the fast-folding protein and RNA structures. For alanine dipeptide, 30ns DBMD simulations captured up to 83-125 times more backbone dihedral transitions than 1µs conventional molecular dynamics (cMD) simulations and were able to accurately reproduce the original free energy profiles. Moreover, DBMD sampled multiple folding and unfolding events within 300ns simulations of the chignolin model protein and identified low-energy conformational states comparable to previous simulation findings. Finally, DBMD captured a general folding pathway of three hairpin RNAs with the GCAA, GAAA, and UUCG tetraloops. Based on Deep Learning neural network, DBMD provides a powerful and generally applicable approach to boosting biomolecular simulations. DBMD is available with open source in OpenMM at https://github.com/MiaoLab20/DBMD/.

Deep Learning Dynamic Allostery of G-Protein-Coupled Receptors.

G-protein-coupled receptors (GPCRs) are the largest superfamily of human membrane proteins and represent primary targets of ~ 1/3 of currently marketed drugs. Allosteric modulators have emerged as more selective drug candidates compared with orthosteric agonists and antagonists. However, many X-ray and cryo-EM structures of GPCRs resolved so far exhibit negligible differences upon binding of positive and negative allosteric modulators (PAMs and NAMs). Mechanism of dynamic allosteric modulation in GPCRs remains unclear. In this work, we have systematically mapped dynamic changes in free energy landscapes of GPCRs upon binding of allosteric modulators using the Gaussian accelerated molecular dynamics (GaMD), Deep Learning (DL) and free energy prOfiling Workflow (GLOW). A total of 18 available high-resolution experimental structures of allosteric modulator-bound class A and B GPCRs were collected for simulations. A number of 8 computational models were generated to examine selectivity of the modulators by changing their target receptors to different subtypes. All-atom GaMD simulations were performed for a total of 66 µs on 44 GPCR systems in the presence/absence of the modulator. DL and free energy calculations revealed significantly reduced conformational space of GPCRs upon modulator binding. While the modulator-free GPCRs often sampled multiple low-energy conformational states, the NAMs and PAMs confined the inactive and active agonist-G protein-bound GPCRs, respectively, to mostly only one specific conformation for signaling. Such cooperative effects were significantly reduced for binding of the selective modulators to "non-cognate" receptor subtypes in the computational models. Therefore, comprehensive DL of extensive GaMD simulations has revealed a general dynamic mechanism of GPCR allostery, which will greatly facilitate rational design of selective allosteric drugs of GPCRs.

Predicting Biomolecular Binding Kinetics: A Review.

Biomolecular binding kinetics including the association (kon) and dissociation (koff) rates are critical parameters for therapeutic design of small-molecule drugs, peptides, and antibodies. Notably, the drug molecule residence time or dissociation rate has been shown to correlate with their efficacies better than binding affinities. A wide range of modeling approaches including quantitative structure-kinetic relationship models, Molecular Dynamics simulations, enhanced sampling, and Machine Learning has been developed to explore biomolecular binding and dissociation mechanisms and predict binding kinetic rates. Here, we review recent advances in computational modeling of biomolecular binding kinetics, with an outlook for future improvements.

Critical Non-Covalent Binding Intermediate for an Allosteric Covalent Inhibitor of SUMO E1.

Post-translational modifications by small ubiquitin-like modifiers (SUMOs) are dysregulated in many types of cancers. The SUMO E1 enzyme has recently been suggested as a new immuno-oncology target. COH000 was recently identified as a highly specific allosteric covalent inhibitor of SUMO E1. However, a marked discrepancy was found between the X-ray structure of the covalent COH000-bound SUMO E1 complex and the available structure-activity relationship (SAR) data of inhibitor analogues due to unresolved noncovalent protein-ligand interactions. Here, we have investigated noncovalent interactions between COH000 and SUMO E1 during inhibitor dissociation through novel Ligand Gaussian accelerated molecular dynamics (LiGaMD) simulations. Our simulations have identified a critical low-energy non-covalent binding intermediate conformation of COH000 that agreed excellently with published and new SAR data of the COH000 analogues, which were otherwise inconsistent with the X-ray structure. Altogether, our biochemical experiments and LiGaMD simulations have uncovered a critical non-covalent binding intermediate during allosteric inhibition of the SUMO E1 complex.

Effects of presenilin-1 familial Alzheimer's disease mutations on γ-secretase activation for cleavage of amyloid precursor protein.

Presenilin-1 (PS1) is the catalytic subunit of γ-secretase which cleaves within the transmembrane domain of over 150 peptide substrates. Dominant missense mutations in PS1 cause early-onset familial Alzheimer's disease (FAD); however, the exact pathogenic mechanism remains unknown. Here we combined Gaussian accelerated molecular dynamics (GaMD) simulations and biochemical experiments to determine the effects of six representative PS1 FAD mutations (P117L, I143T, L166P, G384A, L435F, and L286V) on the enzyme-substrate interactions between γ-secretase and amyloid precursor protein (APP). Biochemical experiments showed that all six PS1 FAD mutations rendered γ-secretase less active for the endoproteolytic (ε) cleavage of APP. Distinct low-energy conformational states were identified from the free energy profiles of wildtype and PS1 FAD-mutant γ-secretase. The P117L and L286V FAD mutants could still sample the "Active" state for substrate cleavage, but with noticeably reduced conformational space compared with the wildtype. The other mutants hardly visited the "Active" state. The PS1 FAD mutants were found to reduce γ-secretase proteolytic activity by hindering APP residue L49 from proper orientation in the active site and/or disrupting the distance between the catalytic aspartates. Therefore, our findings provide mechanistic insights into how PS1 FAD mutations affect structural dynamics and enzyme-substrate interactions of γ-secretase and APP.

Small Molecules Targeting the RNA-Binding Protein HuR Inhibit Tumor Growth in Xenografts.

The RNA-binding protein Hu antigen R (HuR) is a post-transcriptional regulator critical in several types of diseases, including cancer, making it a promising therapeutic target. We have identified small-molecule inhibitors of HuR through a screening approach used in combination with fragment analysis. A total of 36 new compounds originating from fragment linking or structural optimization were studied to establish structure-activity relationships in the set. Two top inhibitors, 1c and 7c, were further validated by binding assays and cellular functional assays. Both block HuR function by directly binding to the RNA-binding pocket, inhibit cancer cell growth dependence of HuR, and suppress cancer cell invasion. Intraperitoneal administration of inhibitor 1c inhibits tumor growth as a single agent and shows a synergistic effect in combination with chemotherapy docetaxel in breast cancer xenograft models. Mechanistically, 1c interferes with the HuR-TGFB/THBS1 axis.