DOX: A new computational protocol for accurate prediction of the protein–ligand binding structures

Molecular docking techniques have now been widely used to predict the protein–ligand binding modes, especially when the structures of crystal complexes are not available. Most docking algorithms are able to effectively generate and rank a large number of probable binding poses. However, it is hard for them to accurately evaluate these poses and identify the most accurate binding structure. In this study, we first examined the performance of some docking programs, based on a testing set made of 15 crystal complexes with drug statins for the human 3‐hydroxy‐3‐methylglutaryl coenzyme A reductase (HMGR). We found that most of the top ranking HMGR–statin binding poses, predicted by the docking programs, were energetically unstable as revealed by the high theoretical‐level calculations, which were usually accompanied by the large deviations from the geometric parameters of the corresponding crystal binding structures. Subsequently, we proposed a new computational protocol, DOX, based on the joint use of molecular Docking, ONIOM, and eXtended ONIOM (XO) methods to predict the accurate binding structures for the protein–ligand complexes of interest. Our testing results demonstrate that the DOX protocol can efficiently predict accurate geometries for all 15 HMGR‐statin crystal complexes without exception. This study suggests a promising computational route, as an effective alternative to the experimental one, toward predicting the accurate binding structures, which is the prerequisite for all the deep understandings of the properties, functions, and mechanisms of the protein–ligand complexes. © 2015 Wiley Periodicals, Inc.     

Pose prediction accuracy in docking studies and enrichment of actives in the active site of GSK-3beta

We present molecular docking studies on the inhibitors of GSK-3beta kinase in the enzyme binding sites of the X-ray complexes (1H8F, 1PYX, 1O9U, 1Q4L, 1Q5K, and 1UV5) using the Schr?dinger docking tool Glide. Cognate and cross-docking studies using standard precision (SP) and extraprecision (XP) algorithms have been carried out. Cognate docking studies demonstrate that docked poses similar to X-ray poses (root-mean-square deviations of less than 2 A) are found within the top four ranks of the GlideScore and E-model scores. However, cross-docking studies typically produce poses that are significantly deviated from X-ray poses in all but a couple of cases, implying potential for induced fit effects in ligand binding. In this light, we have also carried out induced fit docking studies in the active sites of 1O9U, 1Q4L, and 1Q5K. Specifically, conformational changes have been effected in the active sites of these three protein structures to dock noncognate ligands. Thus, for example, the active site of 1O9U has been induced to fit the ligands of 1Q4L, 1Q5K, and 1UV5. These studies produce ligand docked poses which have significantly lower root-mean-square deviations relative to their X-ray crystallographic poses, when compared to the corresponding values from the cross-docking studies. Furthermore, we have used an ensemble of the induced fit models and X-ray structures to enhance the retrieval of active GSK-3beta inhibitors seeded in a decoy database, normally used in Glide validation studies. Thus, our studies provide valuable insights into computational strategies useful for the identification of potential GSK-3beta inhibitors.     

Comprehensive 3D-RISM analysis of the hydration of small molecule binding sites in ligand-free protein structures

Hydration is a critical factor in the ligand binding process. Herein, to examine the hydration states of ligand binding sites, the three-dimensional distribution function for the water oxygen site, g_O( r ) , is computed for 3,706 ligand-free protein structures based on the corresponding small molecule–protein complexes using the 3D-RISM theory. For crystallographic waters (CWs) close to the ligand, g_O( r ) reveals that several CWs are stabilized by interaction networks formed between the ligand, CW, and protein. Based on the g_O( r ) for the crystallographic binding pose of the ligand, hydrogen bond interactions are dominant in the highly hydrated regions while weak interactions such as CH-O are dominant in the moderately hydrated regions. The polar heteroatoms of the ligand occupy the highly hydrated and moderately hydrated regions in the crystallographic (correct) and wrongly docked (incorrect) poses, respectively. Thus, the g_O( r ) of polar heteroatoms may be used to distinguish the correct binding poses.     

Fast and accurate predictions of binding free energies using MM‐PBSA and MM‐GBSA

In the drug discovery process, accurate methods of computing the affinity of small molecules with a biological target are strongly needed. This is particularly true for molecular docking and virtual screening methods, which use approximated scoring functions and struggle in estimating binding energies in correlation with experimental values. Among the various methods, MM‐PBSA and MM‐GBSA are emerging as useful and effective approaches. Although these methods are typically applied to large collections of equilibrated structures of protein‐ligand complexes sampled during molecular dynamics in water, the possibility to reliably estimate ligand affinity using a single energy‐minimized structure and implicit solvation models has not been explored in sufficient detail. Herein, we thoroughly investigate this hypothesis by comparing different methods for the generation of protein‐ligand complexes and diverse methods for free energy prediction for their ability to correlate with experimental values. The methods were tested on a series of structurally diverse inhibitors of Plasmodium falciparum DHFR with known binding mode and measured affinities. The results showed that correlations between MM‐PBSA or MM‐GBSA binding free energies with experimental affinities were in most cases excellent. Importantly, we found that correlations obtained with the use of a single protein‐ligand minimized structure and with implicit solvation models were similar to those obtained after averaging over multiple MD snapshots with explicit water molecules, with consequent save of computing time without loss of accuracy. When applied to a virtual screening experiment, such an approach proved to discriminate between true binders and decoy molecules and yielded significantly better enrichment curves. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2010     

Exponential repulsion improves structural predictability of molecular docking

Molecular docking is a powerful tool for theoretical prediction of the preferred conformation and orientation of small molecules within protein active sites. The obtained poses can be used for estimation of binding energies, which indicate the inhibition effect of designed inhibitors, and therefore might be used for in silico drug design. However, the evaluation of ligand binding affinity critically depends on successful prediction of the native binding mode. Contemporary docking methods are often based on scoring functions derived from molecular mechanical potentials. In such potentials, nonbonded interactions are typically represented by electrostatic interactions between atom‐centered partial charges and standard 6–12 Lennard–Jones potential. Here, we present implementation and testing of a scoring function based on more physically justified exponential repulsion instead of the standard Lennard–Jones potential. We found that this scoring function significantly improved prediction of the native binding modes in proteins bearing narrow active sites such as serine proteases and kinases. © 2016 Wiley Periodicals, Inc.     

Blind docking method combining search of low‐resolution binding sites with ligand pose refinement by molecular dynamics‐based global optimization

This study describes the development of a new blind hierarchical docking method, bhDock, its implementation, and accuracy assessment. The bhDock method uses two‐step algorithm. First, a comprehensive set of low‐resolution binding sites is determined by analyzing entire protein surface and ranked by a simple score function. Second, ligand position is determined via a molecular dynamics‐based method of global optimization starting from a small set of high ranked low‐resolution binding sites. The refinement of the ligand binding pose starts from uniformly distributed multiple initial ligand orientations and uses simulated annealing molecular dynamics coupled with guided force‐field deformation of protein–ligand interactions to find the global minimum. Assessment of the bhDock method on the set of 37 protein–ligand complexes has shown the success rate of predictions of 78%, which is better than the rate reported for the most cited docking methods, such as AutoDock, DOCK, GOLD, and FlexX, on the same set of complexes. © 2009 Wiley Periodicals, Inc. J Comput Chem 2010     

Coarse‐grained molecular dynamics simulations of protein–ligand binding

Coarse‐grained molecular dynamics (CGMD) simulations with the MARTINI force field were performed to reproduce the protein–ligand binding processes. We chose two protein–ligand systems, the levansucrase–sugar (glucose or sucrose), and LinB–1,2‐dichloroethane systems, as target systems that differ in terms of the size and shape of the ligand‐binding pocket and the physicochemical properties of the pocket and the ligand. Spatial distributions of the Coarse‐grained (CG) ligand molecules revealed potential ligand‐binding sites on the protein surfaces other than the real ligand‐binding sites. The ligands bound most strongly to the real ligand‐binding sites. The binding and unbinding rate constants obtained from the CGMD simulation of the levansucrase–sucrose system were approximately 10 times greater than the experimental values; this is mainly due to faster diffusion of the CG ligand in the CG water model. We could obtain dissociation constants close to the experimental values for both systems. Analysis of the ligand fluxes demonstrated that the CG ligand molecules entered the ligand‐binding pockets through specific pathways. The ligands tended to move through grooves on the protein surface. Thus, the CGMD simulations produced reasonable results for the two different systems overall and are useful for studying the protein–ligand binding processes. © 2014 Wiley Periodicals, Inc.     

In Silico Prediction and Validation of CB2 Allosteric Binding Sites to Aid the Design of Allosteric Modulators

Although the 3D structures of active and inactive cannabinoid receptors type 2 (CB2) are available, neither the X-ray crystal nor the cryo-EM structure of CB2-orthosteric ligand-modulator has been resolved, prohibiting the drug discovery and development of CB2 allosteric modulators (AMs). In the present work, we mainly focused on investigating the potential allosteric binding site(s) of CB2. We applied different algorithms or tools to predict the potential allosteric binding sites of CB2 with the existing agonists. Seven potential allosteric sites can be observed for either CB2-CP55940 or CB2-WIN 55,212-2 complex, among which sites B, C, G and K are supported by the reported 3D structures of Class A GPCRs coupled with AMs. Applying our novel algorithm toolset-MCCS, we docked three known AMs of CB2 including Ec2la (C-2), trans-β-caryophyllene (TBC) and cannabidiol (CBD) to each site for further comparisons and quantified the potential binding residues in each allosteric binding site. Sequentially, we selected the most promising binding pose of C-2 in five allosteric sites to conduct the molecular dynamics (MD) simulations. Based on the results of docking studies and MD simulations, we suggest that site H is the most promising allosteric binding site. We plan to conduct bio-assay validations in the future.     

Interaction identification of Zif268 and TATAZF proteins with GC‐/AT‐rich DNA sequence: A theoretical study

Molecular dynamics (MD) simulations for Zif268 (a zinc‐finger‐protein binding specifically to the GC‐rich DNA)‐d(A₁G₂C₃G₄T₅G₆G₇G₈C₉A₁₀C₁₁)₂ and TATA_ZF (a zinc‐finger‐protein recognizing the AT‐rich DNA)‐d(A₁C₂G₃C₄T₅A₆T₇A₈A₉A₁₀A₁₁G₁₂G₁₃)₂ complexes have been performed for investigating the DNA binding affinities and specific recognitions of zinc fingers to GC‐rich and AT‐rich DNA sequences. The binding free energies for the two systems have been further analyzed by using the molecular mechanics Poisson‐Boltzmann surface area (MM‐PBSA) method. The calculations of the binding free energies reveal that the affinity energy of Zif268‐DNA complex is larger than that of TATA_ZF‐DNA one. The affinity between the zinc‐finger‐protein and DNA is mainly driven by more favorable van‐der‐Waals and nonpolar/solvation interactions in both complexes. However, the affinity energy difference of the two binding systems is mainly caused by the difference of van‐der‐Waals interactions and entropy components. The decomposition analysis of MM‐PBSA free energies on each residue of the proteins predicts that the interactions between the residues with the positive charges and DNA favor the binding process; while the interactions between the residues with the negative charges and DNA behave in the opposite way. The interhydrogen‐bonds at the protein‐DNA interface and the induced intrafinger hydrogen bonds between the residues of protein for the Zif268‐DNA complex have been identified at some key contact sites. However, only the interhydrogen‐bonds between the residues of protein and DNA for TATA_ZF‐DNA complex have been found. The interactions of hydrogen‐bonds, electrostatistics and van‐der‐Waals type at some new contact sites have been identified. Moreover, the recognition characteristics of the two studied zinc‐finger‐proteins have also been discussed. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2011     

Evaluation of the searching abilities of HBOP and HBSITE for binding pocket detection

We evaluated the pocket‐searching abilities of the computer programs HBOP and HBSITE, in which hydrophobic potentials calculated on grid points generated around a protein function as an indicator of the pocket‐like property‐using a test set of 458 experimentally observed structures of protein–ligand complexes. The results were compared with those obtained using the alternative algorithms PASS and SiteID, which only consider geometric properties, and Q‐SiteFinder, which considers not only geometry but also physicochemical properties. The comparison revealed that HBOP and HBSITE could detect experimentally observed ligand‐binding pockets for more test complexes than PASS and SiteID. In addition, the success rate of HBOP for detecting binding pockets was higher than that of Q‐SiteFinder, and that of HBSITE was comparable with that of Q‐SiteFinder. Results of tests for ligand‐unbound state proteins indicated that HBSITE was more appropriate than Q‐SiteFinder for pocket searches of ligand‐unbound proteins, and HBSITE was more robust than Q‐SiteFinder against structural differences between ligand‐bound and ‐unbound state proteins. © 2009 Wiley Periodicals, Inc. J Comput Chem 2009     

TargetATPsite: A template‐free method for ATP‐binding sites prediction with residue evolution image sparse representation and classifier ensemble

Understanding the interactions between proteins and ligands is critical for protein function annotations and drug discovery. We report a new sequence‐based template‐free predictor (TargetATPsite) to identify the Adenosine‐5′‐triphosphate (ATP) binding sites with machine‐learning approaches. Two steps are implemented in TargetATPsite: binding residues and pockets predictions, respectively. To predict the binding residues, a novel image sparse representation technique is proposed to encode residue evolution information treated as the input features. An ensemble classifier constructed based on support vector machines (SVM) from multiple random under‐samplings is used as the prediction model, which is effective for dealing with imbalance phenomenon between the positive and negative training samples. Compared with the existing ATP‐specific sequence‐based predictors, TargetATPsite is featured by the second step of possessing the capability of further identifying the binding pockets from the predicted binding residues through a spatial clustering algorithm. Experimental results on three benchmark datasets demonstrate the efficacy of TargetATPsite. © 2013 Wiley Periodicals, Inc.     

A scoring function for docking ligands to low-resolution protein structures

We present a docking method that uses a scoring function for protein-ligand docking that is designed to maximize the docking success rate for low-resolution protein structures. We find that the resulting scoring function parameters are very different depending on whether they were optimized for high- or low-resolution protein structures. We show that this docking method can be successfully applied to predict the ligand-binding site of low-resolution structures. For a set of 25 protein-ligand complexes, in 76% of the cases, more than 50% of ligand-contacting residues are correctly predicted (using receptor crystal structures where the binding site is unspecified). Using decoys of the receptor structures having a 4 A RMSD from the native structure, for the same set of complexes, in 72% of the cases, we obtain at least one correctly predicted ligand-contacting residue. Furthermore, using an 81-protein-ligand set described by Jain, in 76 (93.8%) cases, the algorithm correctly predicts more than 50% of the ligand-contacting residues when native protein structures are used. Using 3 A RMSD from native decoys, in all but two cases (97.5%), the algorithm predicts at least one ligand-binding residue correctly. Finally, compared to the previously published Dolores method, for 298 protein-ligand pairs, the number of cases in which at least half of the specific contacts are correctly predicted is more than four times greater.     

Divalent Sialylated Precision Glycooligomers Binding to Polyomaviruses and the Effect of Different Linkers

Divalent precision glycooligomers terminating in N‐acetylneuraminic acid (Neu5Ac) or 3′‐sialyllactose (3′‐SL) with varying linkers between scaffold and the glycan portions are synthesized via solid phase synthesis for co‐crystallization studies with the sialic acid‐binding major capsid protein VP1 of human Trichodysplasia spinulosa‐associated Polyomavirus. High‐resolution crystal structures of complexes demonstrate that the compounds bind to VP1 depending on the favorable combination of carbohydrate ligand and linker. It is found that artificial linkers can replace portions of natural carbohydrate linkers as long as they meet certain requirements such as size or flexibility to optimize contact area between ligand and receptor binding sites. The obtained results will influence the design of future high affinity ligands based on the structures presented here, and they can serve as a blueprint to develop multivalent glycooligomers as inhibitors of viral adhesion.     

Evaluation of the performance of four molecular docking programs on a diverse set of protein‐ligand complexes

Many molecular docking programs are available nowadays, and thus it is of great practical value to evaluate and compare their performance. We have conducted an extensive evaluation of four popular commercial molecular docking programs, including Glide, GOLD, LigandFit, and Surflex. Our test set consists of 195 protein‐ligand complexes with high‐resolution crystal structures (resolution ≤2.5 Å) and reliable binding data [dissociation constant (K_d) or inhibition constant (K_i)], which are selected from the PDBbind database with an emphasis on diversity. The top‐ranked solutions produced by these programs are compared to the native ligand binding poses observed in crystal structures. Glide and GOLD demonstrate better accuracy than the other two on the entire test set. Their results are also less sensitive to the starting structures for docking. Comparison of the results produced by these programs at three different computation levels reveal that their accuracy are not always proportional to CPU cost as one may expect. The binding scores of the top‐ranked solutions produced by these programs are in low to moderate correlations with experimentally measured binding data. Further analyses on the outcomes of these programs on three suites of subsets of protein‐ligand complexes indicate that these programs are less capable to handle really flexible ligands and relatively flat binding sites, and they have different preferences to hydrophilic/hydrophobic binding sites. Our evaluation can help other researchers to make reasonable choices among available molecular docking programs. It is also valuable for program developers to improve their methods further. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2010     

An NMR Method To Pinpoint Supramolecular Ligand Binding to Basic Residues on Proteins

Targeting protein surfaces involved in protein–protein interactions by using supramolecular chemistry is a rapidly growing field. NMR spectroscopy is the method of choice to map ligand‐binding sites with single‐residue resolution by amide chemical shift perturbation and line broadening. However, large aromatic ligands affect NMR signals over a greater distance, and the binding site cannot be determined unambiguously by relying on backbone signals only. We herein employed Lys‐ and Arg‐specific H2(C)N NMR experiments to directly observe the side‐chain atoms in close contact with the ligand, for which the largest changes in the NMR signals are expected. The binding of Lys‐ and Arg‐specific supramolecular tweezers and a calixarene to two model proteins was studied. The H2(C)N spectra track the terminal CH₂ groups of all Lys and Arg residues, revealing significant differences in their binding kinetics and chemical shift perturbation, and can be used to clearly pinpoint the order of ligand binding.     

Inhibitor docking screened by the modified SAFE_p scoring function: application to cyclic urea HIV-1 PR inhibitors

Our laboratory has in the past developed a method for the prediction of ligand binding free energies to proteins, referred to as SAFE_p (Solvent free energy predictor). Previously, we have applied this protocol for the prediction of the binding free energy of peptidic and cyclic urea HIV-1 PR inhibitors, whose X-ray structures bound to enzyme are known. In this work, we present the first account of a docking simulation, where the ligand conformations were screened and inhibitor ranking was predicted on the basis of a modified SAFE_p approach, for a set of cyclic urea-HIV-1 PR complexes whose structures are not known. We show that the optimal dielectric constant for docking is rather high, in line with the values needed to reproduce some protein residue properties, like pKa's. Our protocol is able to reproduce most of the observed binding ranking, even in the case that the components of the equation are not fitted to experimental data. Partition of the binding free energy into pocket and residue contributions sheds light into the importance of the inhibitor's fragments and on the prediction of "hot spots" for resistance mutations.     

An estimation method of binding free energy in terms of ABEEMσπ/MM and continuum electrostatics fused into LIE method

A method is proposed for the estimation of absolute binding free energy of interaction between proteins and ligands. The linear interaction energy method is combined with atom‐bond electronegativity equalization method at σπ level Force field (fused into molecular mechanics) and generalized Born continuum model calculation of electrostatic solvation for the estimation of the absolute free energy of binding. The parameters of this method are calibrated by using a training set of 24 HIV‐1 protease–inhibitor complexes (PDB entry 1AAQ). A correlation coefficient of 0.93 was obtained with a root mean square deviation of 0.70 kcal mol^?1. This approach is further tested on seven inhibitor and protease complexes, and it provides small root mean square deviation between the calculated binding free energy and experimental binding free energy without reparametrization. By comparing the radii of gyration and the hydrogen bond distances between ligand and protein of three training model molecules, the consistent comparison result of binding free energy is obtained. It proves that this method of calculating the binding free energy with appropriate structural analysis can be applied to quickly assess new inhibitors of HIV‐1 proteases. To test whether the parameters of this method can apply to other drug targets, we have validated this method for the drug target cyclooxygenase‐2. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2011     

Scoring protein–protein docked structures based on the balance and tightness of binding

One main issue in protein-protein docking is to filter or score the putative docked structures. Unlike many popular scoring functions that are based on geometric and energetic complementarity, we present a set of scoring functions that are based on the consideration of local balance and tightness of binding of the docked structures. These scoring functions include the force and moment acting on one component (ligand) imposed by the other (receptor) and the second order spatial derivatives of protein-protein interaction potential. The scoring functions were applied to the docked structures of 19 test targets including enzyme/inhibitor, antibody/antigen and other classes of protein complexes. The results indicate that these scoring functions are also discriminative for the near-native conformation. For some cases, such as antibody/antigen, they show more discriminative efficiency than some other scoring functions, such as desolvation free energy (deltaG(des)) based on pairwise atom-atom contact energy (ACE). The correlation analyses between present scoring functions and the energetic functions also show that there is no clear correlation between them; therefore, the present scoring functions are not essentially the same as energy functions.     

FIPSDock: A new molecular docking technique driven by fully informed swarm optimization algorithm

The accurate prediction of protein–ligand binding is of great importance for rational drug design. We present herein a novel docking algorithm called as FIPSDock, which implements a variant of the Fully Informed Particle Swarm (FIPS) optimization method and adopts the newly developed energy function of AutoDock 4.20 suite for solving flexible protein–ligand docking problems. The search ability and docking accuracy of FIPSDock were first evaluated by multiple cognate docking experiments. In a benchmarking test for 77 protein/ligand complex structures derived from GOLD benchmark set, FIPSDock has obtained a successful predicting rate of 93.5% and outperformed a few docking programs including particle swarm optimization (PSO)@AutoDock, SODOCK, AutoDock, DOCK, Glide, GOLD, FlexX, Surflex, and MolDock. More importantly, FIPSDock was evaluated against PSO@AutoDock, SODOCK, and AutoDock 4.20 suite by cross‐docking experiments of 74 protein–ligand complexes among eight protein targets (CDK2, ESR1, F2, MAPK14, MMP8, MMP13, PDE4B, and PDE5A) derived from Sutherland‐crossdock‐set. Remarkably, FIPSDock is superior to PSO@AutoDock, SODOCK, and AutoDock in seven out of eight cross‐docking experiments. The results reveal that FIPS algorithm might be more suitable than the conventional genetic algorithm‐based algorithms in dealing with highly flexible docking problems. © 2012 Wiley Periodicals, Inc.