Computational prediction of monosaccharide binding free energies to lectins with linear interaction energy models

The linear interaction energy (LIE) method to compute binding free energies is applied to lectin‐monosaccharide complexes. Here, we calculate the binding free energies of monosaccharides to the Ralstonia solanacearum lectin (RSL) and the Pseudomonas aeruginosa lectin‐II (PA‐IIL). The standard LIE model performs very well for RSL, whereas the PA‐IIL system, where ligand binding involves two calcium ions, presents a major challenge. To overcome this, we explore a new variant of the LIE model, where ligand–metal ion interactions are scaled separately. This model also predicts the saccharide binding preference of PA‐IIL on mutation of the receptor, which may be useful for protein engineering of lectins. © 2012 Wiley Periodicals, Inc.     

Solvated interaction energy (SIE) for scoring protein-ligand binding affinities. 2. Benchmark in the CSAR-2010 scoring exercise

Solvated interaction energy (SIE) is an end-point physics-based scoring function for predicting binding affinities from force-field nonbonded interaction terms, continuum solvation, and configurational entropy linear compensation. We tested the SIE function in the Community Structure-Activity Resource (CSAR) scoring challenge consisting of high-resolution cocrystal structures for 343 protein-ligand complexes with high-quality binding affinity data and high diversity with respect to protein targets. Particular emphasis was placed on the sensitivity of SIE predictions to the assignment of protonation and tautomeric states in the complex and the treatment of metal ions near the protein-ligand interface. These were manually curated from an originally distributed CSAR-HiQ data set version, leading to the currently distributed CSAR-NRC-HiQ version. We found that this manual curation was a critical step for accurately testing the performance of the SIE function. The standard SIE parametrization, previously calibrated on an independent data set, predicted absolute binding affinities with a mean-unsigned-error (MUE) of 2.41 kcal/mol for the CSAR-HiQ version, which improved to 1.98 kcal/mol for the upgraded CSAR-NRC-HiQ version. Half-half retraining-testing of SIE parameters on two predefined subsets of CSAR-NRC-HiQ led to only marginal further improvements to an MUE of 1.83 kcal/mol. Hence, we do not recommend altering the current default parameters of SIE at this time. For a sample of SIE outliers, additional calculations by molecular dynamics-based SIE averaging with or without incorporation of ligand strain, by MM-PB(GB)/SA methods with or without entropic estimates, or even by the linear interaction energy (LIE) formalism with an explicit solvent model, did not further improve predictions.     

Evaluation of several two-step scoring functions based on linear interaction energy, effective ligand size, and empirical pair potentials for prediction of protein-ligand binding geometry and free energy

The performances of several two-step scoring approaches for molecular docking were assessed for their ability to predict binding geometries and free energies. Two new scoring functions designed for "step 2 discrimination" were proposed and compared to our CHARMM implementation of the linear interaction energy (LIE) approach using the Generalized-Born with Molecular Volume (GBMV) implicit solvation model. A scoring function S1 was proposed by considering only "interacting" ligand atoms as the "effective size" of the ligand and extended to an empirical regression-based pair potential S2. The S1 and S2 scoring schemes were trained and 5-fold cross-validated on a diverse set of 259 protein-ligand complexes from the Ligand Protein Database (LPDB). The regression-based parameters for S1 and S2 also demonstrated reasonable transferability in the CSARdock 2010 benchmark using a new data set (NRC HiQ) of diverse protein-ligand complexes. The ability of the scoring functions to accurately predict ligand geometry was evaluated by calculating the discriminative power (DP) of the scoring functions to identify native poses. The parameters for the LIE scoring function with the optimal discriminative power (DP) for geometry (step 1 discrimination) were found to be very similar to the best-fit parameters for binding free energy over a large number of protein-ligand complexes (step 2 discrimination). Reasonable performance of the scoring functions in enrichment of active compounds in four different protein target classes established that the parameters for S1 and S2 provided reasonable accuracy and transferability. Additional analysis was performed to definitively separate scoring function performance from molecular weight effects. This analysis included the prediction of ligand binding efficiencies for a subset of the CSARdock NRC HiQ data set where the number of ligand heavy atoms ranged from 17 to 35. This range of ligand heavy atoms is where improved accuracy of predicted ligand efficiencies is most relevant to real-world drug design efforts.     

Continuum solvation models in the linear interaction energy method

The linear interaction energy (LIE) method in combination with two different continuum solvent models has been applied to calculate protein-ligand binding free energies for a set of inhibitors against the malarial aspartic protease plasmepsin II. Ligand-water interaction energies are calculated from both Poisson-Boltzmann (PB) and Generalized Born (GB) continuum models using snapshots from explicit solvent simulations of the ligand and protein-ligand complex. These are compared to explicit solvent calculations, and we find close agreement between the explicit water and PB solvation models. The GB model overestimates the change in solvation energy, and this is caused by consistent underestimation of the effective Born radii in the protein-ligand complex. The explicit solvent LIE calculations and LIE-PB, with our standard parametrization, reproduce absolute experimental binding free energies with an average unsigned error of 0.5 and 0.7 kcal/mol, respectively. The LIE-GB method, however, requires a constant offset to approach the same level of accuracy.     

Accessing the applicability of polarized protein‐specific charge in linear interaction energy analysis

The reliability of the linear interaction energy (LIE) depends on the atomic charge model used to delineate the Coulomb interaction between the ligand and its environment. In this work, the polarized protein‐specific charge (PPC) implementing a recently proposed fitting scheme has been examined in the LIE calculations of the binding affinities for avidin and β‐secretase binding complexes. This charge fitting scheme, termed delta restrained electrostatic potential, bypasses the prevalent numerical difficulty of rank deficiency in electrostatic‐potential‐based charge fitting methods via a dual‐step fitting strategy. A remarkable consistency between the predicted binding affinities and the experimental measurement has been observed. This work serves as a direct evidence of PPC's applicability in rational drug design. © 2014 Wiley Periodicals, Inc.     

Efficient free energy calculations on small molecule host-guest systems - a combined linear interaction energy/one-step perturbation approach

Two efficient methods to calculate binding affinities of ligands with proteins have been critically evaluated by using sixteen small ligand host-guest complexes. It is shown that both the one-step (OS) perturbation method and the linear interaction energy (LIE) method have complementing strengths and weaknesses and can be optimally combined in a new manner. The OS method has a sound theoretical basis to address the free energy of cavity formation, whereas the LIE approach is more versatile and efficient to calculate the free energy of adding charges to such cavities. The off-term, which is neglected in the original LIE equation, can be calculated without additional costs from the OS, offering a powerful synergy between the two methods. The LIE/OS approach presented here combines the best of two worlds and for the model systems studied here, is more accurate than and as efficient as the original methods. It has a sound theoretical background and no longer requires any empirical parameters. The method appears very well suited for application in lead-optimization programmes in drug research, where the structure and dynamics of a series of molecules is of interest, together with an accurate calculation of the binding free energy.     

Rapid prediction of solvation free energy. 3. Application to the SAMPL2 challenge

The SAMPL2 hydration free energy blind prediction challenge consisted of a data set of 41 molecules divided into three subsets: explanatory, obscure and investigatory, where experimental hydration free energies were given for the explanatory, withheld for the obscure, and not known for the investigatory molecules. We employed two solvation models for this challenge, a linear interaction energy (LIE) model based on explicit-water molecular dynamics simulations, and the first-shell hydration (FiSH) continuum model previously calibrated to mimic LIE data. On the 23 compounds from the obscure (blind) dataset, the prospectively submitted LIE and FiSH models provided predictions highly correlated with experimental hydration free energy data, with mean-unsigned-errors of 1.69 and 1.71 kcal/mol, respectively. We investigated several parameters that may affect the performance of these models, namely, the solute flexibility for the LIE explicit-solvent model, the solute partial charging method, and the incorporation of the difference in intramolecular energy between gas and solution phases for both models. We extended this analysis to the various chemical classes that can be formed within the SAMPL2 dataset. Our results strengthen previous findings on the excellent accuracy and transferability of the LIE explicit-solvent approach to predict transfer free energies across a wide spectrum of functional classes. Further, the current results on the SAMPL2 test dataset provide additional support for the FiSH continuum model as a fast yet accurate alternative to the LIE explicit-solvent model. Overall, both the LIE explicit-solvent model and the FiSH continuum solvation model show considerable improvement on the SAMPL2 data set over our previous continuum electrostatics-dispersion solvation model used in the SAMPL1 blind challenge.     

New parameterization approaches of the LIE method to improve free energy calculations of PlmII‐Inhibitors complexes

The standard parameterization of the Linear Interaction Energy (LIE) method has been applied with quite good results to reproduce the experimental absolute binding free energies for several protein–ligand systems. However, we found that this parameterization failed to reproduce the experimental binding free energy of Plasmepsin II (PlmII) in complexes with inhibitors belonging to four dissimilar scaffolds. To overcome this fact, we developed three approaches of LIE, which combine systematic approaches to predict the inhibitor‐specific values of α, β, and γ parameters, to gauge their ability to calculate the absolute binding free energies for these PlmII‐Inhibitor complexes. Specifically: (i) we modified the linear relationship between the weighted nonpolar desolvation ratio (WNDR) and the α parameter, by introducing two models of the β parameter determined by the free energy perturbation (FEP) method in the absence of the constant term γ, and (ii) we developed a new parameterization model to investigate the linear correlation between WNDR and the correction term γ. Using these parameterizations, we were able to reproduce the experimental binding free energy from these systems with mean absolute errors lower than 1.5 kcal/mol. © 2010 Wiley Periodicals, Inc. J Comput Chem 2010     

The performance of ANI-ML potentials for ligand-n(H2O) interaction energies and estimation of hydration free energies from end-point MD simulations

Here, we investigate the performance of “Accurate NeurAl networK engINe for Molecular Energies” (ANI), trained on small organic compounds, on bulk systems including non-covalent interactions and applicability to estimate solvation (hydration) free energies using the interaction between the ligand and explicit solvent (water) from single-step MD simulations. The method is adopted from ANI using the Atomic Simulation Environment (ASE) and predicts the non-covalent interaction energies at the accuracy of wb97x/6-31G(d) level by a simple linear scaling for the conformations sampled by molecular dynamics (MD) simulations of ligand-n(H₂O) systems. For the first time, we test ANI potentials' abilities to reproduce solvation free energies using linear interaction energy (LIE) formulism by modifying the original LIE equation. Our results on ~250 different complexes show that the method can be accurate and have a correlation of R² = 0.88–0.89 (MAE <1.0 kcal/mol) to the experimental solvation free energies, outperforming current end-state methods. Moreover, it is competitive to other conventional free energy methods such as FEP and BAR with 15-20 × fold reduced computational cost.     

Improving binding mode and binding affinity predictions of docking by ligand-based search of protein conformations: evaluation in D3R grand challenge 2015

The growing number of protein–ligand complex structures, particularly the structures of proteins co-bound with different ligands, in the Protein Data Bank helps us tackle two major challenges in molecular docking studies: the protein flexibility and the scoring function. Here, we introduced a systematic strategy by using the information embedded in the known protein–ligand complex structures to improve both binding mode and binding affinity predictions. Specifically, a ligand similarity calculation method was employed to search a receptor structure with a bound ligand sharing high similarity with the query ligand for the docking use. The strategy was applied to the two datasets (HSP90 and MAP4K4) in recent D3R Grand Challenge 2015. In addition, for the HSP90 dataset, a system-specific scoring function (ITScore2_hsp90) was generated by recalibrating our statistical potential-based scoring function (ITScore2) using the known protein–ligand complex structures and the statistical mechanics-based iterative method. For the HSP90 dataset, better performances were achieved for both binding mode and binding affinity predictions comparing with the original ITScore2 and with ensemble docking. For the MAP4K4 dataset, although there were only eight known protein–ligand complex structures, our docking strategy achieved a comparable performance with ensemble docking. Our method for receptor conformational selection and iterative method for the development of system-specific statistical potential-based scoring functions can be easily applied to other protein targets that have a number of protein–ligand complex structures available to improve predictions on binding.     

Solvent interaction energy calculations on molecular dynamics trajectories: increasing the efficiency using systematic frame selection

End-point methods such as linear interaction energy (LIE) analysis, molecular mechanics generalized Born solvent-accessible surface (MM/GBSA), and solvent interaction energy (SIE) analysis have become popular techniques to calculate the free energy associated with protein-ligand binding. Such methods typically use molecular dynamics (MD) simulations to generate an ensemble of protein structures that encompasses the bound and unbound states. The energy evaluation method (LIE, MM/GBSA, or SIE) is subsequently used to calculate the energy of each member of the ensemble, thus providing an estimate of the average free energy difference between the bound and unbound states. The workflow requiring both MD simulation and energy calculation for each frame and each trajectory proves to be computationally expensive. In an attempt to reduce the high computational cost associated with end-point methods, we study several methods by which frames may be intelligently selected from the MD simulation including clustering and address the question of how the number of selected frames influences the accuracy of the SIE calculations.     

Free energy simulations of uncatalyzed DNA replication fidelity: structure and stability of T.G and dTTP.G terminal DNA mismatches flanked by a single dangling nucleotide

A reference system for DNA replication fidelity was studied by free energy perturbation (FEP) and linear interaction energy (LIE) methods. The studied system included a hydrated duplex DNA with the 5'-CG dangling end of the templating strand, and dCTP4-.Mg2+ or dTTP4-.Mg2+ inserted opposite the dangling G to form a correct (i.e., Watson-Crick) or incorrect (i.e., wobble) base pair, respectively. The average distance between the 3'-terminal oxygen of the primer strand and the alpha-phosphorus of dNTP was found to be 0.2 A shorter for the correct base pair than for the incorrect base pair. Binding of the incorrect dNTP was found to be disfavored by 0.4 kcal/mol relative to the correct dNTP. We estimated that improved binding and more near-attack configurations sampled by the correct base pair should translate in aqueous solution and in the absence of DNA polymerase into a six times faster rate for the incorporation of the correct dNTP into DNA. The accuracy of the calculated binding free energy difference was verified by examining the relative free energy for melting duplex DNA containing GC and GT terminal base pairs flanked by a 5' dangling C. The calculated LIE and FEP free energies of 1.7 and 1.1 kcal/mol, respectively, compared favorably with the experimental estimate of 1.4 kcal/mol obtained using the nearest neighbor parameters. To decompose the calculated free energies into additive electrostatic and van der Waals contributions and to provide a set of rigorous theoretical data for the parametrization of the LIE method, we suggested a variant of the FEP approach, for which we coined a binding-relevant free energy (BRFE) acronym. BRFE approach is characterized by its unique perturbation pathway and by its exclusion of the intramolecular energy of a rigid part of the ligand from the total potential energy.     

Binding affinity prediction with different force fields: examination of the linear interaction energy method

A systematic study of the linear interaction energy (LIE) method and the possible dependence of its parameterization on the force field and system (receptor binding site) is reported. We have calculated the binding free energy for nine different ligands in complex with P450cam using three different force fields (Amber95, Gromos87, and OPLS-AA). The results from these LIE calculations using our earlier parameterization give relative free energies of binding that agree remarkably well with the experimental data. However, the absolute energies are too positive for all three force fields, and it is clear that an additional constant term (gamma) is required in this case. Out of five examined LIE models, the same one emerges as the best for all three force fields, and this, in fact, corresponds to our earlier one apart from the addition of the constant gamma, which is almost identical for the three force fields. Thus, the present free energy calculations clearly indicate that the coefficients of the LIE method are independent of the force field used. Their relation to solvation free energies is also demonstrated. The only free parameter of the best model is gamma, which is found to depend on the hydrophobicity of the binding site. We also attempt to quantify the binding site hydrophobicity of four different proteins which shows that the ordering of gamma's for these sites reflects the fraction of hydrophobic surface area.     

A challenging system: free energy prediction for factor Xa

Factor Xa (fXa) is a promising target for antithrombotic drugs. Recently, we presented a molecular dynamics study on fXa, which highlighted the need for a careful system setup to obtain stable simulations. Here, we show that these simulations can be used to predict the free energy of binding of several fXa inhibitors. We tested molecular mechanics/Poisson-Boltzmann surface area, molecular mechanics/Generalized Born surface area, and linear interaction energy (LIE) on a small data set of fXa ligands. The continuum solvent approaches only yield satisfying correlations to the experimental results if some of the water molecules are explicitly included in the free energy calculations. LIE gave reasonable results if a sufficiently large data set is used. In general, our procedure of setting up the fXa simulation system enabled MD simulations, which produce adequate ensembles for free energy calculations.     

In silico binding free energy predictability by using the linear interaction energy (LIE) method: bromobenzimidazole CK2 inhibitors as a case study

Protein kinase CK2 is essential for cell viability, and its control regards a broad series of cellular events such as gene expression, RNA, and protein synthesis. Evidence of its involvement in tumor development and viral replication indicates CK2 as a potential target of antineoplastic and antiviral drugs. In this study the Linear Interaction Energy (LIE) Method with the Surface Generalized Born (SGB) continuum solvation model was used to study several bromobenzimidazole CK2 inhibitors. This methodology, developed by Aqvist, finds a plausible compromise between accuracy and computational speed in evaluating binding free energy (DeltaGbind) values. In this study, two different free binding energy models, named "CK2scoreA" and "CK2scoreB", were developed using 22 inhibitors as the training set in a stepwise approach useful to appropriately select both the tautomeric form and the starting binding position of each inhibitor. Both models are statistically acceptable. Indeed, the better one is characterized by a correlation coefficient (r2) of 0.81, and the predictive accuracy was 0.65 kcal/mol. The corresponding validation, using an external test set of 16 analogs, showed a correlation coefficient (q2) of 0.68 and a prediction root-mean-square error of 0.78 kcal/mol. In this case, the LIE approach has been proved to be an efficient methodology to rationalize the difference of activity, the key interactions, and the different possible binding modes of this specific class of potent CK2 inhibitors.     

Are automated molecular dynamics simulations and binding free energy calculations realistic tools in lead optimization? An evaluation of the linear interaction energy (LIE) method

An extensive evaluation of the linear interaction energy (LIE) method for the prediction of binding affinity of docked compounds has been performed, with an emphasis on its applicability in lead optimization. An automated setup is presented, which allows for the use of the method in an industrial setting. Calculations are performed for four realistic examples, retinoic acid receptor gamma, matrix metalloprotease 3, estrogen receptor alpha, and dihydrofolate reductase, focusing on different aspects of the procedure. The obtained LIE models are evaluated in terms of the root-mean-square (RMS) errors from experimental binding free energies and the ability to rank compounds appropriately. The results are compared to the best empirical scoring function, selected from a set of 10 scoring functions. In all cases, good LIE models can be obtained in terms of free-energy RMS errors, although reasonable ranking of the ligands of dihydrofolate reductase proves difficult for both the LIE method and scoring functions. For the other proteins, the LIE model results in better predictions than the best performing scoring function. These results indicate that the LIE approach, as a tool to evaluate docking results, can be a valuable asset in computational lead optimization programs.     

Assessment and acceleration of binding energy calculations for protein–ligand complexes by the fragment molecular orbital method

In the field of drug discovery, it is important to accurately predict the binding affinities between target proteins and drug applicant molecules. Many of the computational methods available for evaluating binding affinities have adopted molecular mechanics‐based force fields, although they cannot fully describe protein–ligand interactions. A noteworthy computational method in development involves large‐scale electronic structure calculations. Fragment molecular orbital (FMO) method, which is one of such large‐scale calculation techniques, is applied in this study for calculating the binding energies between proteins and ligands. By testing the effects of specific FMO calculation conditions (including fragmentation size, basis sets, electron correlation, exchange‐correlation functionals, and solvation effects) on the binding energies of the FK506‐binding protein and 10 ligand complex molecule, we have found that the standard FMO calculation condition, FMO2‐MP2/6‐31G(d), is suitable for evaluating the protein–ligand interactions. The correlation coefficient between the binding energies calculated with this FMO calculation condition and experimental values is determined to be R = 0.77. Based on these results, we also propose a practical scheme for predicting binding affinities by combining the FMO method with the quantitative structure–activity relationship (QSAR) model. The results of this combined method can be directly compared with experimental binding affinities. The FMO and QSAR combined scheme shows a higher correlation with experimental data (R = 0.91). Furthermore, we propose an acceleration scheme for the binding energy calculations using a multilayer FMO method focusing on the protein–ligand interaction distance. Our acceleration scheme, which uses FMO2‐HF/STO‐3G:MP2/6‐31G(d) at R_int = 7.0 Å, reduces computational costs, while maintaining accuracy in the evaluation of binding energy. © 2015 Wiley Periodicals, Inc.     

Computational Analysis of Triazole-Based Kojic Acid Analogs as Tyrosinase Inhibitors by Molecular Dynamics and Free Energy Calculations

Molecular docking, molecular dynamics (MD) simulations and the linear interaction energy (LIE) method were used here to predict binding modes and free energy for a set of 1,2,3-triazole-based KA analogs as potent inhibitors of Tyrosinase (TYR), a key metalloenzyme of the melanogenesis process. Initially, molecular docking calculations satisfactorily predicted the binding mode of evaluated KA analogs, where the KA part overlays the crystal conformation of the KA inhibitor into the catalytic site of TYR. The MD simulations were followed by the LIE method, which reproduced the experimental binding free energies for KA analogs with an r² equal to 0.97, suggesting the robustness of our theoretical model. Moreover, the van der Waals contributions performed by some residues such as Phe197, Pro201, Arg209, Met215 and Val218 are responsible for the binding recognition of 1,2,3-triazole-based KA analogs in TYR catalytic site. Finally, our calculations provide suitable validation of the combination of molecular docking, MD, and LIE approaches as a powerful tool in the structure-based drug design of new and potent TYR inhibitors.     

Estimation of relative binding free energy based on a free energy variational principle for the FKBP-ligand system

Predicting an accurate binding free energy between a target protein and a ligand can be one of the most important steps in a drug discovery process. Often, many molecules must be screened to find probable high potency ones. Thus, a computational technique with low cost is highly desirable for the estimation of binding free energies of many molecules. Several techniques have thus far been developed for estimating binding free energies. Some techniques provide accurate predictions of binding free energies but high large computational cost. Other methods give good predictions but require tuning of some parameters to predict them with high accuracy. In this study, we propose a method to predict relative binding free energies with accuracy comparable to the results of prior methods but with lower computational cost and with no parameter needing to be carefully tuned. Our technique is based on the free energy variational principle. FK506 binding protein (FKBP) with 18 ligands is taken as a test system. Our results are compared to those from other widely used techniques. Our method provides a correlation coefficient (r ² ) of 0.80 between experimental and calculated relative binding free energies and yields an average absolute error of 0.70 kcal/mol compared to experimental values. These results are comparable to or better than results from other techniques. We also discuss the possibility to improve our method further.     

Exhaustive search and solvated interaction energy (SIE) for virtual screening and affinity prediction

We carried out a prospective evaluation of the utility of the SIE (solvation interaction energy) scoring function for virtual screening and binding affinity prediction. Since experimental structures of the complexes were not provided, this was an exercise in virtual docking as well. We used our exhaustive docking program, Wilma, to provide high-quality poses that were rescored using SIE to provide binding affinity predictions. We also tested the combination of SIE with our latest solvation model, first shell of hydration (FiSH), which captures some of the discrete properties of water within a continuum model. We achieved good enrichment in virtual screening of fragments against trypsin, with an area under the curve of about 0.7 for the receiver operating characteristic curve. Moreover, the early enrichment performance was quite good with 50% of true actives recovered with a 15% false positive rate in a prospective calculation and with a 3% false positive rate in a retrospective application of SIE with FiSH. Binding affinity predictions for both trypsin and host-guest complexes were generally within 2 kcal/mol of the experimental values. However, the rank ordering of affinities differing by 2 kcal/mol or less was not well predicted. On the other hand, it was encouraging that the incorporation of a more sophisticated solvation model into SIE resulted in better discrimination of true binders from binders. This suggests that the inclusion of proper Physics in our models is a fruitful strategy for improving the reliability of our binding affinity predictions.