Ligand Binding Free Energy Research Articles

Efficient calculation of protein-ligand binding free energy using GFN methods: the power of the cluster model.

Protein-ligand interactions are crucial in many biochemical processes and biomedical applications, yet accurately calculating the binding free energy of the interactions still remains challenging. In this work, we systematically investigate the performance of a generic force field GFN-FF and some semi-empirical quantum mechanical (SQM) methods (GFNn, n = 0, 1, 2) in terms of the accuracy of the calculated binding free energy. It is found that the performance of the GFN-FF method is quite good in a neutral-ligand system since the Pearson correlation coefficient (rp) is 0.70 and the mean absolute error (MAE) is 5.49 kcal mol-1. However, it may fail in a charged-ligand system (the MAE is 18.98 kcal mol-1). Moreover, we also propose a cluster model (i.e., truncating the protein at a given cutoff) along with the SQM method in the GFN family. Importantly, the GFN2-xTB shows the best performance among the SQM methods (the MAE is 4.91 kcal mol-1 and 10.25 kcal mol-1 in the neutral-ligand and charged-ligand systems, respectively), much better than GFN-FF in the charged-ligand system. Notably, the computing cost of the GFN2-xTB in the appropriate cluster model is even lower than that of the GFN-FF (in the entire complex). The present study sheds some light on the potential power of the GFN family in the efficient calculation of the binding free energy in bio-systems.

An open-source molecular builder and free energy preparation workflow

Automated free energy calculations for the prediction of binding free energies of congeneric series of ligands to a protein target are growing in popularity, but building reliable initial binding poses for the ligands is challenging. Here, we introduce the open-source FEgrow workflow for building user-defined congeneric series of ligands in protein binding pockets for input to free energy calculations. For a given ligand core and receptor structure, FEgrow enumerates and optimises the bioactive conformations of the grown functional group(s), making use of hybrid machine learning/molecular mechanics potential energy functions where possible. Low energy structures are optionally scored using the gnina convolutional neural network scoring function, and output for more rigorous protein–ligand binding free energy predictions. We illustrate use of the workflow by building and scoring binding poses for ten congeneric series of ligands bound to targets from a standard, high quality dataset of protein–ligand complexes. Furthermore, we build a set of 13 inhibitors of the SARS-CoV-2 main protease from the literature, and use free energy calculations to retrospectively compute their relative binding free energies. FEgrow is freely available at https://github.com/cole-group/FEgrow, along with a tutorial.

Ligand-Induced Conformational and Dynamical Changes in a GT-B Glycosyltransferase: Molecular Dynamics Simulations of Heptosyltransferase I Complexes.

Understanding the dynamical motions and ligand recognition motifs of heptosyltransferase I (HepI) can be critical to discerning the behavior of other glycosyltransferase (GT) enzymes. Prior studies in our lab have demonstrated that GTs in the GT-B structural class, which are characterized by their connection of two Rossman-like domains by a linker region, have conserved structural fold and dynamical motions, despite low sequence homology, therefore making discoveries found in HepI transferable to other GT-B enzymes. Through molecular dynamics simulations and ligand binding free energy analysis of HepI in the apo and bound complexes (for all kinetically relevant combinations of the native substrates/products), we have determined the energetically favored enzymatic pathway for ligand binding and release. Our principal component, dynamic cross correlation, and network analyses of the simulations have revealed correlated motions involving residues within the N-terminal domain communicating with C-terminal domain residues via both proximal amino acid residues and also functional groups of the bound substrates. Analyses of the structural changes, energetics of substrate/product binding, and changes in pKa have elucidated a variety of inter and intradomain interactions that are critical for enzyme catalysis. These data corroborate our experimental observations of protein conformational changes observed in both presteady state kinetic and circular dichroism analyses of HepI. These simulations provided invaluable structural insights into the regions involved in HepI conformational rearrangement upon ligand binding. Understanding the specific interactions governing conformational changes is likely to enhance our efforts to develop novel dynamics disrupting inhibitors against GT-B structural enzymes in the future.

Probing Allosteric Modulation of Membrane Receptor in the Native State by Data Mining-Integrated Tracking Microscopy

Probing Allosteric Modulation of Membrane Receptor in the Native State by Data Mining-Integrated Tracking Microscopy

Investigation of phytoconstituents of Enicostemma littorale as potential glucokinase activators through molecular docking for the treatment of type 2 diabetes mellitus.

Glucokinase (GK) is an enzyme involved in synthesising glucose into glucose-6 phosphate and serves a crucial function in glucose sensing. Therefore, agents thatinduce GK activation could be used to treat T2DM. The present work has been carried out to investigate the GK activation potential of phytoconstituents of Enicostemma littorale through molecular docking. All the phytoconstituents have been screened through the Lipinski rule of 5, Veber's rule, and ADMET properties. From these initial screening, only Apigenin, Ferulic acid, Genkwanin, p-coumaric acid, Protocatechuic acid, Syringic acid, and Vanillic acid have been selected to perform molecular docking studies. The binding free energy and binding mode of the native ligand in the allosteric site of the enzyme have been considered the reference for the other molecules' validation. The native ligand has exhibited - 7.2kcal/mol binding free energy, whereas; it has formed four hydrogen bonds with THR-228, LYS-169, ASP-78, and GLY-81. Based on these findings, the interactions of phytoconstituents have been justified. Apigenin, genkwanin, and swertiamarin exhibited - 8.7, - 7.5, and - 8.3kcal/mol binding free energy, respectively, which indicates better enzyme activation than the native ligand. Swertiamarin has formed 08 hydrogen bonds with allosteric amino acid residues, which confirms the excellent enzyme activation by these phytoconstituents. We concluded that if we can isolate and consume the exact active phytoconstituents (GK activators) from this plant, we can use them effectively to treat T2DM. More GK activators can be developed by considering them as a natural lead moiety.

Prediction of Binding Free Energy of Protein-Ligand Complexes with a Hybrid Molecular Mechanics/Generalized Born Surface Area and Machine Learning Method.

Accurate prediction of protein–ligand binding free energies is important in enzyme engineering and drug discovery. The molecular mechanics/generalized Born surface area (MM/GBSA) approach is widely used to estimate ligand-binding affinities, but its performance heavily relies on the accuracy of its energy components. A hybrid strategy combining MM/GBSA and machine learning (ML) has been developed to predict the binding free energies of protein–ligand systems. Based on the MM/GBSA energy terms and several features associated with protein–ligand interactions, our ML-based scoring function, GXLE, shows much better performance than MM/GBSA without entropy. In particular, the good transferability of the GXLE model is highlighted by its good performance in ranking power for prediction of the binding affinity of different ligands for either the docked structures or crystal structures. The GXLE scoring function and its code are freely available and can be used to correct the binding free energies computed by MM/GBSA.

CHARMM-GUI Free Energy Calculator for Practical Ligand Binding Free Energy Simulations with AMBER.

Alchemical free energy methods, such as free energy perturbation (FEP) and thermodynamic integration (TI), become increasingly popular and crucial for drug design and discovery. However, the system preparation of alchemical free energy simulation is an error-prone, time-consuming, and tedious process for a large number of ligands. To address this issue, we have recently presented CHARMM-GUI Free Energy Calculator that can provide input and postprocessing scripts for NAMD and GENESIS FEP molecular dynamics systems. In this work, we extended three submodules of Free Energy Calculator to work with the full suite of GPU-accelerated alchemical free energy methods and tools in AMBER, including input and postprocessing scripts. The BACE1 (β-secretase 1) benchmark set was used to validate the AMBER-TI simulation systems and scripts generated by Free Energy Calculator. The overall results of relatively large and diverse systems are almost equivalent with different protocols (unified and split) and with different timesteps (1, 2, and 4 fs), with R2 > 0.9. More importantly, the average free energy differences between two protocols are small and reliable with four independent runs, with a mean unsigned error (MUE) below 0.4 kcal/mol. Running at least four independent runs for each pair with AMBER20 (and FF19SB/GAFF2.1/OPC force fields), we obtained a MUE of 0.99 kcal/mol and root-mean-square error of 1.31 kcal/mol for 58 alchemical transformations in comparison with experimental data. In addition, a set of ligands for T4-lysozyme was used to further validate our free energy calculation protocol whose results are close to experimental data (within 1 kcal/mol). In summary, Free Energy Calculator provides a user-friendly web-based tool to generate the AMBER-TI system and input files for high-throughput binding free energy calculations with access to the full set of GPU-accelerated alchemical free energy, enhanced sampling, and analysis methods in AMBER.

Calculation of CYP450 protein-ligand binding and dissociation free energy paths.

The function of an enzyme depends on its dynamic structure, and the catalytic mechanism has long been an active focus of research. The principle for interpreting protein selectivity and fidelity stems from optimization of the active site upon protein-substrate complexation, i.e., a lock-and-key configuration, on which most protein-substrate molecule binding recognition, and hence drug discovery, relies. Yet another thought has been to incorporate the protein folding interior tunnels for stereo- and regio-selectivity along the protein-substrate or protein-ligand/inhibitor binding process. Free energy calculations provide valuable information for molecular recognition and protein-ligand binding dynamics and kinetics. In this study, we focused on the kinetics of cytochrome P450 proteins (CYP450s) and the protein interior tunnel structure-dynamics relationship in terms of the substrate binding and leaving mechanism. A case in point is given by the prostaglandin H2 (PGH2) homologous isomerase of prostacyclin synthase. To calculate the reactant and product traversing the tunnels to and from the heme site, the free energy paths and tunnel potentials of mean force are constructed from steered molecular dynamics simulations and adaptive basing force umbrella sampling simulations. We explore the binding tunnels and critical residue lining characteristics for the ligand traverse and the underlying mechanism of CYP450 activity. Our theoretical analysis provides insights into the decisive role of the substrate tunnel binding process of the CYP450 mechanism and may be useful in drug design and protein engineering contexts.

VAD-MM/GBSA: A Variable Atomic Dielectric MM/GBSA Model for Improved Accuracy in Protein-Ligand Binding Free Energy Calculations.

The molecular mechanics/generalized Born surface area (MM/GBSA) has been widely used in end-point binding free energy prediction in structure-based drug design (SBDD). However, in practice, it is usually being treated as a disputed method mostly because of its system dependence. Here, combining with machine-learning optimization, we developed a novel version of MM/GBSA, named variable atomic dielectric MM/GBSA (VAD-MM/GBSA), by assigning variable dielectric constants directly to the protein/ligand atoms. The new strategy exhibits markedly improved accuracy in binding affinity calculations for various protein-ligand systems and is promising to be used in the postprocessing of structure-based virtual screening. Moreover, VAD-MM/GBSA outperformed prime MM/GBSA in Schrödinger software and showed remarkable predictive performance for specific protein targets, such as POL polyprotein, human immunodeficiency virus type 1 (HIV-1) protease, etc. Our study showed that the VAD-MM/GBSA method with little extra computational overhead provides a potential replacement of the MM/GBSA in AMBER software. An online web server of VAD-MMGBSA has been developed and is now available at http://cadd.zju.edu.cn/vdgb.

Benchmark of Popular Free Energy Approaches Revealing the Inhibitors Binding to SARS-CoV-2 Mpro.

The COVID-19 pandemic has killed millions of people worldwide since its outbreak in December 2019. The pandemic is caused by the SARS-CoV-2 virus whose main protease (Mpro) is a promising drug target since it plays a key role in viral proliferation and replication. Currently, developing an effective therapy is an urgent task, which requires accurately estimating the ligand-binding free energy to SARS-CoV-2 Mpro. However, it should be noted that the accuracy of a free energy method probably depends on the protein target. A highly accurate approach for some targets may fail to produce a reasonable correlation with the experiment when a novel enzyme is considered as a drug target. Therefore, in this context, the ligand-binding affinity to SARS-CoV-2 Mpro was calculated via various approaches. The molecular docking approach was manipulated using Autodock Vina (Vina) and Autodock4 (AD4) protocols to preliminarily investigate the ligand-binding affinity and pose to SARS-CoV-2 Mpro. The binding free energy was then refined using the fast pulling of ligand (FPL), linear interaction energy (LIE), molecular mechanics-Poisson–Boltzmann surface area (MM-PBSA), and free energy perturbation (FEP) methods. The benchmark results indicated that for docking calculations, Vina is more accurate than AD4, and for free energy methods, FEP is the most accurate method, followed by LIE, FPL, and MM-PBSA (FEP > LIE ≈ FPL > MM-PBSA). Moreover, atomistic simulations revealed that the van der Waals interaction is the dominant factor. The residues Thr26, His41, Ser46, Asn142, Gly143, Cys145, His164, Glu166, and Gln189 are essential elements affecting the binding process. Our benchmark provides guidelines for further investigations using computational approaches.

Structure-Based Virtual Screening and De Novo Design of PIM1 Inhibitors with Anticancer Activity from Natural Products.

Background: the proviral insertion site of Moloney murine leukemia (PIM) 1 kinase has served as a therapeutic target for various human cancers due to the enhancement of cell proliferation and the inhibition of apoptosis. Methods: to identify effective PIM1 kinase inhibitors, structure-based virtual screening of natural products of plant origin and de novo design were carried out using the protein–ligand binding free energy function improved by introducing an adequate dehydration energy term. Results: as a consequence of subsequent enzyme inhibition assays, four classes of PIM1 kinase inhibitors were discovered, with the biochemical potency ranging from low-micromolar to sub-micromolar levels. The results of extensive docking simulations showed that the inhibitory activity stemmed from the formation of multiple hydrogen bonds in combination with hydrophobic interactions in the ATP-binding site. Optimization of the biochemical potency by chemical modifications of the 2-benzylidenebenzofuran-3(2H)-one scaffold led to the discovery of several nanomolar inhibitors with antiproliferative activities against human breast cancer cell lines. Conclusions: these new PIM1 kinase inhibitors are anticipated to serve as a new starting point for the development of anticancer medicine.

Variational Method for Networkwide Analysis of Relative Ligand Binding Free Energies with Loop Closure and Experimental Constraints.

We describe an efficient method for the simultaneous solution of all free energies within a relative binding free-energy (RBFE) network with cycle closure and experimental/reference constraint conditions using Bennett Acceptance Ratio (BAR) and Multistate BAR (MBAR) analysis. Rather than solving the BAR or MBAR equations for each transformation independently, the simultaneous solution of all transformations are obtained by performing a constrained minimization of a global objective function. The nonlinear optimization of the objective function is subjected to affine linear constraints that couple the free energies between the network edges. The constraints are used to enforce the closure of thermodynamic cycles within the RBFE network, and to enforce an additional set of linear constraint conditions demonstrated here to be subsets of (1 or 2) experimental values. We describe details of the practical implementation of the network BAR/MBAR procedure, including use of generalized coordinates in the minimization of the free-energy objective function, propagation of bootstrap errors from those coordinates, and performance and memory optimization. In some cases it is found that use of restraints in the optimization is more practical than use of generalized coordinates for enforcing constraint conditions. The fast BARnet and MBARnet methods are used to analyze the RBFEs of six prototypical protein-ligand systems, and it is shown that enforcement of cycle closure conditions reduces the error in the predictions only modestly, and further reduction in errors can be achieved when one or two experimental RBFEs are included in the optimization procedure. These methods have been implemented into FE-ToolKit, a new free-energy analysis toolkit. The BARnet/MBARnet framework presented here opens the door to new, more efficient and robust free-energy analysis with enhanced predictive capability for drug discovery applications.

Solvation Contribution to the Free Energy of Ligand Binding is Determined by the Coupling Between Surface Hydration Structure and Side Chain Motion

Surface hydration significantly influences binding of a ligand to a protein. Yet its contribution to the binding free energy has been difficult to quantify. We introduce a computational method to calculate the desolvation penalty upon ligand binding based on the analysis of the surface hydration structure in all-atom molecular dynamics simulation. As a test system, we use the Src homology 3 (SH3) domain and two intrinsically disordered proline-rich motifs as ligands, one from the cellular protein cAbl and the other from NS1 of the 1918 Spanish influenza a virus.

Ligand Gaussian Accelerated Molecular Dynamics (LiGaMD) for 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 “LiGaMD”, which selectively boosts the ligand non-bonded 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.

Potential effective inhibitory compounds against Prostate Specific Membrane Antigen (PSMA): A molecular docking and molecular dynamics study

One of the most prevalent cancers in men is prostate cancer and could be managed with immunotoxins or antibody treatment. Because of the substantial rise of the Prostate-Specific Antigen and the Prostate-Specific Membrane Antigen (PSMA), cancer vaccination should be rendered with these antigens. Through pharmacodynamic experiments in a library of natural compounds from ZINC database, the current research sought to identify compounds that could suppress PSMA protein. To test the most productive compounds for further research, the Library has been scanned with Pharmacophore and ADMET analysis followed by molecular docking methods in the first phase. After selecting 15 ligands with the best pose related to docking results, to evaluate the stability of the ligand-protein bounds of the compounds, a molecular dynamics simulation considering the effect of the presence of zinc ions on the protein structure was performed. The measurement of ligand binding modes and free energy has shown that four compounds, including Z10, Z06, Z01, and Z03, have formed critical interactions with the active site's residues. Besides, multiple approaches were employed to determine their inhibition rating and describe the variables that facilitate the attachment of ligands to the protein active site. The results are obtained from the MMPBSA/GBSA analysis of four selected small molecules (Z10, Z06, Z01, and Z03), which are very close to the IC50 value of reference ligand (DCIBzl); they are −13.85 kcal/mol, −12.58 kcal/mol, −10.71 kcal/mol and −9.39 kcal/mol respectively. Finally, we evaluate the results obtained from selected ligands using hydrogen bond and decomposition analyzes. We have examined the effective interactions between ligands and S1/S1ˊpockets in protein. Our computational results illustrate the design of more efficient inhibitors of PSMA.

The role of water in host-guest interaction

One of the main applications of atomistic computer simulations is the calculation of ligand binding free energies. The accuracy of these calculations depends on the force field quality and on the thoroughness of configuration sampling. Sampling is an obstacle in simulations due to the frequent appearance of kinetic bottlenecks in the free energy landscape. Very often this difficulty is circumvented by enhanced sampling techniques. Typically, these techniques depend on the introduction of appropriate collective variables that are meant to capture the system’s degrees of freedom. In ligand binding, water has long been known to play a key role, but its complex behaviour has proven difficult to fully capture. In this paper we combine machine learning with physical intuition to build a non-local and highly efficient water-describing collective variable. We use it to study a set of host-guest systems from the SAMPL5 challenge. We obtain highly accurate binding free energies and good agreement with experiments. The role of water during the binding process is then analysed in some detail.

Alchemical absolute protein-ligand binding free energies for drug design.

The recent advances in relative protein–ligand binding free energy calculations have shown the value of alchemical methods in drug discovery. Accurately assessing absolute binding free energies, although highly desired, remains a challenging endeavour, mostly limited to small model cases. Here, we demonstrate accurate first principles based absolute binding free energy estimates for 128 pharmaceutically relevant targets. We use a novel rigorous method to generate protein–ligand ensembles for the ligand in its decoupled state. Not only do the calculations deliver accurate protein–ligand binding affinity estimates, but they also provide detailed physical insight into the structural determinants of binding. We identify subtle rotamer rearrangements between apo and holo states of a protein that are crucial for binding. When compared to relative binding free energy calculations, obtaining absolute binding free energies is considerably more challenging in large part due to the need to explicitly account for the protein in its apo state. In this work we present several approaches to obtain apo state ensembles for accurate absolute ΔG calculations, thus outlining protocols for prospective application of the methods for drug discovery.

Rational Computational Design of Fourth-Generation EGFR Inhibitors to Combat Drug-Resistant Non-Small Cell Lung Cancer.

Although the inhibitors of singly mutated epidermal growth factor receptor (EGFR) kinase are effective for the treatment of non-small cell lung cancer (NSCLC), their clinical efficacy has been limited due to the emergence of various double and triple EGFR mutants with drug resistance. It has thus become urgent to identify potent and selective inhibitors of triple mutant EGFRs resistant to first-, second-, and third-generation EGFR inhibitors. Herein, we report the discovery of potent and highly selective inhibitors of EGFR exon 19 p.E746_A750del/EGFR exon 20 p.T790M/EGFR exon 20 p.C797S (d746-750/T790M/C797S) mutant, which were derived via two-track virtual screening and de novo design. This two-track approach was performed so as to maximize and minimize the inhibitory activity against the triple mutant and the wild type, respectively. Extensive chemical modifications of the initial hit compounds led to the identification of several low-nanomolar inhibitors of the d746-750/T790M/C797S mutant. Among them, two compounds exhibited more than 104-fold selectivity in the inhibition of EGFRd746-750/T790M/C797S over the wild type. The formations of a hydrogen bond with the mutated residue Ser797 and the van der Waals contact with the mutated residue Met790 were found to be a common feature in the interactions between EGFRd746-750/T790M/C797S and the fourth-generation inhibitors. Such an exceptionally high selectivity could also be attributed to the formation of the hydrophobic contact with a Gly loop residue or the hydrogen bond with Asp855 in the activation loop. The discovery of the potent and selective EGFRd746-750/T790M/C797S inhibitors were actually made possible by virtue of the modified protein–ligand binding free energy function involving a new hydration free energy term with enhanced accuracy. The fourth-generation EGFR inhibitors found in this work are anticipated to serve as a new starting point for the discovery of anti-NSCLC medicines to overcome the problematic drug resistance.

Druggability assessment of precursor membrane protein as a target for inhibiting the Zika virus

The Zika virus (ZIKV), a significant zoonotic flavivirus, was neglected as a human pathogen until the recent epidemic. The rapid geographic spread of the virus and association with neurological disorders has created a global public health concern pressing the need for anti-ZIKV drugs. Previous ZIKV drug discovery research has focused on three primary targets, RNA-dependent RNA polymerase, envelope protein, and viral proteases, and none has yet resulted in a commercially viable inhibitor. In the quest for finding effective inhibitors, it is important to expand the number of targets available for drug discovery research. To this end, the ZIKV precursor membrane protein (prM) comes to the forefront as a potential target due to its critical role in virus infectivity and pathogenicity. prM acts as a chaperone for envelope protein folding and prevents premature fusion of virions to the host membrane and has not been attempted as a drug target before. One critical requirement for a protein to be an effective target is the ability of the protein to be druggable, i.e. having active sites that can bind to specific ligands. In this work, the druggability of prM was assessed via molecular docking combined molecular dynamics simulations followed binding affinity kinetics studies. Compounds that had a high affinity to the prM protein were screened in silico and ligand-binding free energies were computed using molecular mechanics with generalized Born and surface area continuum solvation (MM-GBSA) method. In vitro binding kinetics via biolayer interferometry (BLI) and interaction analysis confirmed that prM could be targeted for drug discovery to combat ZIKV infection. Communicated by Ramaswamy H. Sarma

Kinetics and free energy of ligand dissociation using weighted ensemble milestoning.

We consider the recently developed weighted ensemble milestoning (WEM) scheme [D. Ray and I. Andricioaei, J. Chem. Phys. 152, 234114 (2020)] and test its capability of simulating ligand-receptor dissociation dynamics. We performed WEM simulations on the following host-guest systems: Na+/Cl- ion pair and 4-hydroxy-2-butanone ligand with FK506 binding protein. As a proof of principle, we show that the WEM formalism reproduces the Na+/Cl- ion pair dissociation timescale and the free energy profile obtained from long conventional MD simulation. To increase the accuracy of WEM calculations applied to kinetics and thermodynamics in protein-ligand binding, we introduced a modified WEM scheme called weighted ensemble milestoning with restraint release (WEM-RR), which can increase the number of starting points per milestone without adding additional computational cost. WEM-RR calculations obtained a ligand residence time and binding free energy in agreement with experimental and previous computational results. Moreover, using the milestoning framework, the binding time and rate constants, dissociation constants, and committor probabilities could also be calculated at a low computational cost. We also present an analytical approach for estimating the association rate constant (kon) when binding is primarily diffusion driven. We show that the WEM method can efficiently calculate multiple experimental observables describing ligand-receptor binding/unbinding and is a promising candidate for computer-aided inhibitor design.