Absolute Binding Free Energy Research Articles

Molecular dynamics and absolute binding free energy studies of piperine derivatives as potential inhibitors of SARS-CoV-2 main protease

The work presents a library of piperine derivatives as potential inhibitors of the main protease protein (Mpro) functionality using Docking Studies, Molecular Dynamics (MD) Simulations and Absolute Binding Free-Energy calculations. 342 ligands were selected for this work and docked with Mpro protein. Among all the ligands studied, PIPC270, PIPC299, PIPC252, PIPC63, PIPC311 were the top five docked conformations having significant hydrogen bonding and hydrophobic interactions inside the active pocket of Mpro. These top five ligands were subjected to MD simulations for 100 ns using GROMACS. Root Mean Square Deviation (RMSD), Root Mean Square Fluctuation (RMSF), Radius of Gyration (Rg), Solvent Accessible Surface Area (SASA) and hydrogen bond analysis revealed that the ligands bounded to protein remain stable without significant deviations during the course of MD simulations. Absolute binding free energy (ΔGb) was calculated for theses complexes and found that the ligand PIPC299 shows the prevalent binding affinity with binding free-energy of about −113.05 Kcal/mol. Thus, these molecules can be further tested by in vitro and in vivo studies on Mpro. This study lays a path to explore the new functionality of piperine derivatives as novel drug like molecules.Communicated by Ramaswamy H. Sarma

Development and test of highly accurate endpoint free energy methods. 1: Evaluation of ABCG2 charge model on solvation free energy prediction and optimization of atom radii suitable for more accurate solvationfreeenergy prediction by the PBSA method.

Accurate estimation of solvation free energy (SFE) lays the foundation for accurate prediction of binding free energy. The Poisson-Boltzmann (PB) or generalized Born (GB) combined with surface area (SA) continuum solvation method (PBSA and GBSA) have been widely used in SFE calculations because they can achieve good balance between accuracy and efficiency. However, the accuracy of these methods can be affected by several factors such as the charge models, polar and nonpolar SFE calculation methods and the atom radii used in the calculation. In this work, the performance of the ABCG2 (AM1-BCC-GAFF2) charge model as well as other two charge models, that is, RESP (Restrained Electrostatic Potential) and AM1-BCC (Austin Model 1-bond charge corrections), on the SFE prediction of 544 small molecules in water by PBSA/GBSA was evaluated. In order to improve the performance of the PBSA prediction based on the ABCG2 charge, we further explored the influence of atom radii on the prediction accuracy and yielded a set of atom radius parameters for more accurate SFE prediction using PBSA based on the ABCG2/GAFF2 by reproducing the thermodynamic integration (TI) calculation results. The PB radius parameters of carbon, oxygen, sulfur, phosphorus, chloride, bromide and iodine, were adjusted. New atom types, on, oi, hn1, hn2, hn3, were introduced to further improve the fitting performance. Then, we tuned the parameters in the nonpolar SFE model using the experimental SFE data and the PB calculation results. By adopting the new radius parameters and new nonpolar SFE model, the root mean square error (RMSE) of the SFE calculation for the 544 molecules decreased from 2.38 to 1.05 kcal/mol. Finally, the new radius parameters were applied in the prediction of protein-ligand binding free energies using the MM-PBSA method. For the eight systems tested, we could observe higher correlation between the experiment data and calculation results and smaller prediction errors for the absolute binding free energies, demonstrating that our new radius parameters can improve the free energy calculation using the MM-PBSA method.

Recent Advances in Alchemical Binding Free Energy Calculations for Drug Discovery.

Rigorous physics-based methods to calculate binding free energies of protein-ligand complexes have become a valued component of structure-based drug design. Relative and absolute binding free energy calculations have been deployed prospectively in support of solving diverse drug discovery challenges. Here we review recent applications of binding free energy calculations to fragment growing and linking, scaffold hopping, binding pose validation, virtual screening, covalent enzyme inhibition, and positional analogue scanning. Furthermore, we discuss the merits of using protein models and highlight recent efforts to replace costly binding free energy calculations with predictions from machine learning models trained on a limited number of free energy perturbation or thermodynamic integration calculations thereby allowing for extended chemical space exploration.

Absolute binding free energy calculations for molecules binding to the GPCR GCGR at the protein/bilayer/water interface.

Numerous structures have been solved of small molecules binding at the protein-lipid interface of membrane proteins, such as GPCRs and ion channels. Allosteric modulators that bind at these lipid-exposed sites provide a novel way to target these therapeutically important proteins with greater specificity. Computational tools play a major role in modern drug discovery pipelines, but these tools must be validated, extended, and developed before being applied prospectively to lipid-exposed sites. The creation of a small molecule with a high affinity for its target is a key step in the drug discovery process and the accurate prediction of affinity remains a central challenge in computational chemistry. Alchemical absolute binding free energy (ABFE) calculations have been shown to be highly accurate in predicting the binding affinity of small molecules to soluble proteins. However, ABFE calculations remain challenging, with little work assessing their application to membrane proteins. Our work has focused on extending ABFE calculations to lipid exposed sites. The binding free energy calculations for these sites are complicated by the fact the binding site is ‘solvated’ by lipids. We present a study on the GPCR GCGR in complex with the small molecule MK-0893. MK-0893 sits at the protein/bilayer/water interface. This, coupled with the size and charge of the molecule, makes ABFE calculations for this system particularly challenging. We discuss the adaptations we made to our ABFE workflow to overcome these challenges. Particularly, we make suggestions as to how to monitor the generation of a ‘correct’ apo state for the non-interacting ligand, with regards to the occupancy of the binding site and side-chain rotamers. These adaptations significantly improved our results and show a good agreement with experimental values from a publicly available SAR series.

Binding free energies from the polarizable AMOEBA force field.

A protocol for computing absolute binding free energies with the AMOEBA polarizable force field is presented. Examples discussed include various host-guest complexes from the SAMPL community exercises, binding of ions to an RNA G-quadruplex and a series of protein-ligand interactions. In many cases where calculations are not limited by sampling issues, the AMOEBA model is able to provide binding values within “chemical accuracy,” i.e., roughly 1 kcal/mol from experiment.

Computational fragment screening of dimeric SARS-CoV-2 main protease.

During the coronavirus pandemic the lack of antiviral therapeutics drove intense interest in rapid drug discovery by the scientific, and medical communities. The Voelz lab took part in the COVID Moonshot initiative, which focused on virtually screening small drug-like fragments as leads for inhibitors of the main protease (Mpro) of SARS-CoV-2, identified using high throughput crystallographic techniques. The Voelz lab performed absolute binding free energy calculations on the [email protected] distributed computing platform to massively parallelize the molecular dynamic (MD) simulations. Crystallographic fragment-based screening techniques can be time consuming, costly, and dependent on the availability of protein crystals. Could fragment based drug discovery be performed completely in silico? To address this question, we aim to perform simulations to virtually screen the original fragment library used in the COVID Moonshot, and assess the accuracy of the predictions. The GROMACS simulation package was used to construct a simulation topology of the apo form of dimeric Mpro from its crystal structure. The correct protonation states of titratable residues were taken from published experimental results. As a first step, we prepared and performed a set of simulations of Mpro in solution, each with one of 11 drug fragments selected from the crystallographic fragment screen, all of which are known to co-crystallize with Mpro experimentally. The simulations are analyzed to determine where the drug fragments like to bind to the Mpro receptor. These results are compared results to the experimental co-crystal structures to determine how accurately MD simulations can predict the known binding pose of small drug like fragments.

Expanded ensemble predictions of absolute binding free energies in the SAMPL9 host-guest challenge.

As part of the SAMPL9 community-wide blind host-guest challenge, we implemented an expanded ensemble workflow to predict absolute binding free energies for 13 small molecules against pillar[6]arene. Notable features of our protocol include consideration of a variety of protonation and enantiomeric states for both host and guests, optimization of alchemical intermediates, and analysis of free energy estimates and their uncertainty using large numbers of simulation replicates performed using distributed computing. Our predictions of absolute binding free energies resulted in a mean absolute error of 2.29 kcal mol-1 and an R2 of 0.54. Overall, results show that expanded ensemble calculations using all-atom molecular dynamics simulations are a valuable and efficient computational tool in predicting absolute binding free energies.

Calculations of the binding free energies of the Comprehensive in vitro Proarrhythmia Assay (CiPA) reference drugs to cardiac ion channels.

The evaluation of the inhibitory activities of drugs on multiple cardiac ion channels is required for the accurate assessment of proarrhythmic risks. Moreover, the in silico prediction of such inhibitory activities of drugs on cardiac channels can improve the efficiency of the drug-development process. Here, we performed molecular docking simulations to predict the complex structures of 25 reference drugs that were proposed by the Comprehensive in vitro Proarrhythmia Assay consortium using two cardiac ion channels, the human ether-a-go-go-related gene (hERG) potassium channel and human NaV1.5 (hNaV1.5) sodium channel, with experimentally available structures. The absolute binding free energy (ΔGbind) values of the predicted structures were calculated by a molecular dynamics-based method and compared with the experimental half-maximal inhibitory concentration (IC50) data. Furthermore, the regression analysis between the calculated values and negative of the common logarithm of the experimental IC50 values (pIC50) revealed that the calculated values of four and ten drugs deviated significantly from the regression lines of the hERG and hNaV1.5 channels, respectively. We reconsidered the docking poses and protonation states of the drugs based on the experimental data and recalculated their ΔGbind values. Finally, the calculated ΔGbind values of 24 and 19 drugs correlated with their experimental pIC50 values (coefficients of determination=0.791 and 0.613 for the hERG and hNaV1.5 channels, respectively). Thus, the regression analysis between the calculated ΔGbind and experimental IC50 data ensured the realization of an increased number of reliable complex structures.

The role of loop dynamics in the prediction of ligand-protein binding enthalpy.

The enthalpic and entropic components of ligand-protein binding free energy reflect the interactions and dynamics between ligand and protein. Despite decades of study, our understanding and hence our ability to predict these individual components remains poor. In recent years, there has been substantial effort and success in the prediction of relative and absolute binding free energies, but the prediction of the enthalpic (and entropic) contributions in biomolecular systems remains challenging. Indeed, it is not even clear what kind of performance in terms of accuracy could currently be obtained for such systems. It is, however, relatively straight-forward to compute the enthalpy of binding. We thus evaluated the performance of absolute enthalpy of binding calculations using molecular dynamics simulation for ten inhibitors against a member of the bromodomain family, BRD4-1, against isothermal titration calorimetry data. Initial calculations, with the AMBER force-field showed good agreement with experiment (R2 = 0.60) and surprisingly good accuracy with an average of root-mean-square error (RMSE) = 2.49 kcal mol-1. Of the ten predictions, three were obvious outliers that were all over-predicted compared to experiment. Analysis of various simulation factors, including parameterization, buffer concentration and conformational dynamics, revealed that the behaviour of a loop (the ZA loop on the periphery of the binding site) strongly dictates the enthalpic prediction. Consistent with previous observations, the loop exists in two distinct conformational states and by considering one or the other or both states, the prediction for the three outliers can be improved dramatically to the point where the R2 = 0.95 and the accuracy in terms of RMSE improves to 0.90 kcal mol-1. However, performance across force-fields is not consistent: if OPLS and CHARMM are used, different outliers are observed and the correlation with the ZA loop behaviour is not recapitulated, likely reflecting parameterization as a confounding problem. The results provide a benchmark standard for future study and comparison.

Taming multiple binding poses in alchemical binding free energy prediction: the β-cyclodextrin host-guest SAMPL9 blinded challenge.

We apply the Alchemical Transfer Method (ATM) and a bespoke fixed partial charge force field to the SAMPL9 bCD host-guest binding free energy prediction challenge that comprises a combination of complexes formed between five phenothiazine guests and two cyclodextrin hosts. Multiple chemical forms, competing binding poses, and computational modeling challenges pose significant obstacles to obtaining reliable computational predictions for these systems. The phenothiazine guests exist in solution as racemic mixtures of enantiomers related by nitrogen inversions that bind the hosts in various binding poses, each requiring an individual free energy analysis. Due to the large size of the guests and the conformational reorganization of the hosts, which prevent a direct absolute binding free energy route, binding free energies are obtained by a series of absolute and relative binding alchemical steps for each chemical species in each binding pose. Metadynamics-accelerated conformational sampling was found to be necessary to address the poor convergence of some numerical estimates affected by conformational trapping. Despite these challenges, our blinded predictions quantitatively reproduced the experimental affinities for the β-cyclodextrin host and, to a lesser extent, those with a methylated derivative. The work illustrates the challenges of obtaining reliable free energy data in in silico drug design for even seemingly simple systems and introduces some of the technologies available to tackle them.

Discovery of allosteric SHP2 inhibitors through ensemble-based consensus molecular docking, endpoint and absolute binding free energy calculations

SHP2 (Src homology-2 domain-containing protein tyrosine phosphatase-2) is a cytoplasmic protein -tyrosine phosphatase encoded by the gene PTPN11. It plays a crucial role in regulating cell growth and differentiation. Specifically, SHP2 is an oncoprotein associated with developmental pathologies and several different cancer types, including gastric, leukemia and breast cancer and is of great therapeutic interest. Given these roles, current research efforts have focused on developing SHP2 inhibitors. Allosteric SHP2 inhibitors have been shown to be more selective and pharmacologically appealing compared to competitive catalytic inhibitors targeting SHP2. Nevertheless, there remains a need for novel allosteric inhibitor scaffolds targeting SHP2 to develop compounds with improved selectivity, cell permeability, and bioavailability. Towards this goal, this study applied various computational tools to screen over 6 million compounds against the allosteric site within SHP2. The top-ranked hits from our in-silico screening were validated using protein thermal shift and biolayer interferometry assays, revealing three potent compounds. Kinetic binding assays were employed to measure the binding affinities of the top-ranked compounds and demonstrated that they all bind to SHP2 with a nanomolar affinity. Hence the compounds and the computational workflow described herein provide an effective approach for identifying and designing a generation of improved allosteric inhibitors of SHP2.

Absolute binding free energies of the antiviral peptide ATN-161 with protein targets of SARS-CoV-2

The interactions of the antiviral pentapeptide ATN-161 with the closed and open conformations of the α5β1 integrin, the SARS-CoV-2 major protease, and the omicron variant spike protein complexed with hACE2 were studied using molecular docking and molecular dynamics simulation. Molecular docking was performed to obtain ATN-161 binding poses with these studied protein targets. Subsequently, molecular dynamics simulations were performed to verify the ligand stability at the binding site of each protein target. Pulling simulations, umbrella sampling, and weighted histogram analysis method were used to obtain the potential of mean force of each system and calculate the Gibbs free energy of binding for the ATN-161 peptide in each binding site of these protein targets. The results showed that ATN-161 binds to α5β1 integrin in its active and inactive form, binds weakly to the omicron variant spike protein complexed with hACE2, and strongly binds to the main protease target. Communicated by Ramaswamy H. Sarma

Interaction preferences between protein side chains and key epigenetic modifications 5-methylcytosine, 5-hydroxymethycytosine and N6-methyladenine

Covalent modifications of standard DNA/RNA nucleobases affect epigenetic regulation of gene expression by modulating interactions between nucleic acids and protein readers. We derive here the absolute binding free energies and analyze the binding modalities between key modified nucleobases 5-methylcytosine (5mC), 5-hydroxymethylcytosine (5hmC) and N6-methyladenine (m6A) and all non-prolyl/non-glycyl protein side chains using molecular dynamics simulations and umbrella sampling in both water and methanol, the latter mimicking the low dielectric environment at the dehydrated nucleic-acid/protein interfaces. We verify the derived affinities by comparing against a comprehensive set of high-resolution structures of nucleic-protein complexes involving 5mC. Our analysis identifies protein side chains that are highly tuned for detecting cytosine methylation as a function of the environment and can thus serve as microscopic readers of epigenetic marks. Conversely, we show that the relative ordering of sidechain affinities for 5hmC and m6A does not differ significantly from those for their precursor bases, cytosine and adenine, respectively, especially in the low dielectric environment. For those two modified bases, the effect is more nuanced and manifests itself primarily at the level of absolute changes in the binding free energy. Our results contribute towards establishing a quantitative foundation for understanding, predicting and modulating the interactions between modified nucleic acids and proteins at the atomistic level.

Free Energy Calculations Using the Movable Type Method with Molecular Dynamics Driven Protein-Ligand Sampling.

Fast and accurate biomolecular free energy estimation has been a significant interest for decades, and with recent advances in computer hardware, interest in new method development in this field has even grown. Thorough configurational state sampling using molecular dynamics (MD) simulations has long been applied to the estimation of the free energy change corresponding to the receptor-ligand complexing process. However, performing large-scale simulation is still a computational burden for the high-throughput hit screening. Among molecular modeling tools, docking and scoring methods are widely used during the early stages of the drug discovery process in that they can rapidly generate discrete receptor-ligand binding modes and their individual binding affinities. Unfortunately, the lack of thorough conformational sampling in docking and scoring protocols leads to difficulty discovering global minimum binding modes on a complicated energy landscape. The Movable Type (MT) method is a novel absolute binding free energy approach which has demonstrated itself to be robust across a wide range of targets and ligands. Traditionally, the MT method is used with protein-ligand binding modes generated with rigid-receptor or flexible-receptor (induced fit) docking protocols; however, these protocols are by their nature less likely to be effective with more highly flexible targets or with those situations in which binding involves multiple step pathways. In these situations, more thorough samplings are required to better explain the free energy of binding. Therefore, to explore the prediction capability and computational efficiency of the MT method when using more thorough protein-ligand conformational sampling protocols, in the present work, we introduced a series of binding mode modeling protocols ranging from conventional docking routines to single-trajectory conventional molecular dynamics (cMD) and parallel Monte Carlo molecular dynamics (MCMD). Through validation against several structurally and mechanistically diverse protein-ligand test sets, we explore the performance of the MT method as a virtual screening tool to work with the docking protocols and as an MD simulation-based binding free energy tool.

Electrostatic Contributions to the Binding Free Energy of Nicotine to the Acetylcholine Binding Protein.

Biomolecular binding relies on specific attractive interactions between two partner molecules, including electrostatics, dispersion, hydrophobicity, and solvation. Assessing the contributions of electrostatic interactions to binding is key to the understanding of ligand binding mechanisms and the design of improved biomolecular binders. For example, nicotine is a well-known agonist of nicotinic acetylcholine receptors (nAChRs), but the molecular mechanisms for the differential action of nicotine on brain and muscle nAChRs remain elusive. In this work, we have chosen the acetylcholine binding protein (AChBP) in complex with nicotine as a model system to interrogate the electrostatic contributions to nicotine binding. Our absolute binding free energy simulations confirm that nicotine binds AChBP predominantly in its protonated (charged) form. By comparing energetic contributions from decomposed interactions for either neutral or charged nicotine, our calculations shed light on the nature of the binding of nicotine to the AChBP. The preferred binding of charged nicotine over neutral nicotine originates from its stronger electrostatic interactions with AChBP, a cation-π interaction to a tryptophan residue and a hydrogen bond between nicotine and the backbone carbonyl of the tryptophan, whereas the major force driving the binding process appears to be van der Waals interactions. The various nonelectrostatic terms can also indirectly modulate the electrostatic interactions through fine-tuning the binding pose of the ligand in the binding site, providing an explanation of why the binding specificity of nicotine to the brain versus muscle nAChRs is driven by electrostatic interaction, given that the immediate binding site residues, including the key tryptophan residue, are identical in the two receptors.

Absolute Binding Free Energy Calculations for Buried Water Molecules.

Water often plays a key role in mediating protein-ligand interactions. Understanding contributions from active-site water molecules to binding thermodynamics of a ligand is important in predicting binding free energies for ligand optimization. In this work, we tested a non-equilibrium switching method for absolute binding free energy calculations on water molecules in binding sites of 13 systems. We discuss the lessons we learned about identified issues that affected our calculations and ways to address them. This work fits with our larger focus on how to do accurate ligand binding free energy calculations when water rearrangements are very slow, such as rearrangements due to ligand modification (as in relative free energy calculations) or ligand binding (as in absolute free energy calculations). The method studied in this work can potentially be used to account for limited water sampling via providing endpoint corrections to free energy calculations using our calculated binding free energy of water.

Meta-Analysis Reveals That Absolute Binding Free-Energy Calculations Approach Chemical Accuracy.

Systematic and quantitative analysis of the reliability of formally exact methods that in silico calculate absolute protein-ligand binding free energies remains lacking. Here, we provide, for the first time, evidence-based information on the reliability of these methods by statistically studying 853 cases from 34 different research groups through meta-analysis. The results show that formally exact methods approach chemical accuracy (error = 1.58 kcal/mol), even if people are challenging difficult tasks such as blind drug screening in recent years. The geometrical-pathway-based methods prove to possess a better convergence ability than the alchemical ones, while the latter have a larger application range. We also reveal the importance of always using the latest force fields to guarantee reliability and discuss the pros and cons of turning to an implicit solvent model in absolute binding free-energy calculations. Moreover, based on the meta-analysis, an evidence-based guideline for in silico binding free-energy calculations is provided.

Identification of Novel Ribonucleotide Reductase Inhibitors for Therapeutic Application in Bile Tract Cancer: An Advanced Pharmacoinformatics Study

Biliary tract cancer (BTC) is constituted by a heterogeneous group of malignant tumors that may develop in the biliary tract, and it is the second most common liver cancer. Human ribonucleotide reductase M1 (hRRM1) has already been proven to be a potential BTC target. In the current study, a de novo design approach was used to generate novel and effective chemical therapeutics for BTC. A set of comprehensive pharmacoinformatics approaches was implemented and, finally, seventeen potential molecules were found to be effective for the modulation of hRRM1 activity. Molecular docking, negative image-based ShaEP scoring, absolute binding free energy, in silico pharmacokinetics, and toxicity assessments corroborated the potentiality of the selected molecules. Almost all molecules showed higher affinity in comparison to gemcitabine and naphthyl salicylic acyl hydrazone (NSAH). On binding interaction analysis, a number of critical amino acids was found to hold the molecules at the active site cavity. The molecular dynamics (MD) simulation study also indicated the stability between protein and ligands. High negative MM-GBSA (molecular mechanics generalized Born and surface area) binding free energy indicated the potentiality of the molecules. Therefore, the proposed molecules might have the potential to be effective therapeutics for the management of BTC.

Evaluating the use of absolute binding free energy in the fragment optimisation process

Key to the fragment optimisation process within drug design is the need to accurately capture the changes in affinity that are associated with a given set of chemical modifications. Due to the weakly binding nature of fragments, this has proven to be a challenging task, despite recent advancements in leveraging experimental and computational methods. In this work, we evaluate the use of Absolute Binding Free Energy (ABFE) calculations in guiding fragment optimisation decisions, retrospectively calculating binding free energies for 59 ligands across 4 fragment elaboration campaigns. We first demonstrate that ABFEs can be used to accurately rank fragment-sized binders with an overall Spearman’s r of 0.89 and a Kendall τ of 0.67, although often deviating from experiment in absolute free energy values with an RMSE of 2.75 kcal/mol. We then also show that in several cases, retrospective fragment optimisation decisions can be supported by the ABFE calculations. Comparing against cheaper endpoint methods, namely Nwat-MM/GBSA, we find that ABFEs offer better ranking power and correlation metrics. Our results indicate that ABFE calculations can usefully guide fragment elaborations to maximise affinity.

Structure Based in Silico Screening Revealed a Potent Acinetobacter Baumannii Ftsz Inhibitor From Asinex Antibacterial Library.

The superbug Acinetobacter baumannii is an increasingly prevalent pathogen of the intensive care units where its treatment is challenging. The identification of newer drug targets and the development of propitious therapeutics against this pathogen is of utmost importance. A drug target, cell division enzyme (FtsZ), involved in A. baumannii cytokinesis is a promising avenue for antibacterial therapy. Structure based virtual screening illustrated a lead-like molecule from Asinex antibacterial library to have the best binding affinity for the FtsZ active pocket. Computational pharmacokinetics predicted the compound to have the safest pharmacokinetics profile, thus maximizing the chances of the molecule reaching the market with enhanced efficacy and lesser toxicity. Molecular dynamics simulations in an aqueous environment revealed the flexibility of protein loop regions, and upward extension followed by the backward movement of the inhibitor N, N-dimethylpyridazin-3-amine ring on its axis. The active pocket residue Thr310 demonstrated to play significant role in inhibitor binding. The binding free energy predicted by MM/GBSA and MM/PBSA reflected system stability with a total value of -62.15 kcal/mol and -10.60 kcal/mol, respectively. The absolute binding free energy estimated by WaterSwap was -16 kcal/mol that validates and affirms complex stability. The inhibitor represents a promising scaffold as a lead optimization for the FtsZ enzyme.