Performance of HADDOCK and a simple contact-based protein–ligand binding affinity predictor in the D3R Grand Challenge 2

We present the performance of HADDOCK, our information-driven docking software, in the second edition of the D3R Grand Challenge. In this blind experiment, participants were requested to predict the structures and binding affinities of complexes between the Farnesoid X nuclear receptor and 102 different ligands. The models obtained in Stage1 with HADDOCK and ligand-specific protocol show an average ligand RMSD of 5.1 Å from the crystal structure. Only 6/35 targets were within 2.5 Å RMSD from the reference, which prompted us to investigate the limiting factors and revise our protocol for Stage2. The choice of the receptor conformation appeared to have the strongest influence on the results. Our Stage2 models were of higher quality (13 out of 35 were within 2.5 Å), with an average RMSD of 4.1 Å. The docking protocol was applied to all 102 ligands to generate poses for binding affinity prediction. We developed a modified version of our contact-based binding affinity predictor PRODIGY, using the number of interatomic contacts classified by their type and the intermolecular electrostatic energy. This simple structure-based binding affinity predictor shows a Kendall’s Tau correlation of 0.37 in ranking the ligands (7th best out of 77 methods, 5th/25 groups). Those results were obtained from the average prediction over the top10 poses, irrespective of their similarity/correctness, underscoring the robustness of our simple predictor. This results in an enrichment factor of 2.5 compared to a random predictor for ranking ligands within the top 25%, making it a promising approach to identify lead compounds in virtual screening.     

LEAP: Highly accurate prediction of protein loop conformations by integrating coarse‐grained sampling and optimized energy scores with all‐atom refinement of backbone and side chains

Prediction of protein loop conformations without any prior knowledge (ab initio prediction) is an unsolved problem. Its solution will significantly impact protein homology and template‐based modeling as well as ab initio protein‐structure prediction. Here, we developed a coarse‐grained, optimized scoring function for initial sampling and ranking of loop decoys. The resulting decoys are then further optimized in backbone and side‐chain conformations and ranked by all‐atom energy scoring functions. The final integrated technique called loop prediction by energy‐assisted protocol achieved a median value of 2.1 Å root mean square deviation (RMSD) for 325 12‐residue test loops and 2.0 Å RMSD for 45 12‐residue loops from critical assessment of structure‐prediction techniques (CASP) 10 target proteins with native core structures (backbone and side chains). If all side‐chain conformations in protein cores were predicted in the absence of the target loop, loop‐prediction accuracy only reduces slightly (0.2 Å difference in RMSD for 12‐residue loops in the CASP target proteins). The accuracy obtained is about 1 Å RMSD or more improvement over other methods we tested. The executable file for a Linux system is freely available for academic users at http://sparks‐lab.org . © 2013 Wiley Periodicals, Inc.     

logD7.4 modeling using Bayesian Regularized Neural Networks. Assessment and correction of the errors of prediction

Bayesian Regularized Neural Networks (BRNNs) employing Automatic Relevance Determination (ARD) are used to construct a predictive model for the distribution coefficient logD7.4 from an in-house data set of 5000 compounds with experimental endpoints. A method for assessing the accuracy of prediction is established based upon a query compound's distance to the training set. logD7.4 predictions are also dynamically corrected with an associated library of compounds of continuously updated, experimentally measured logD7.4 values. A comparison of local models and associated libraries comprising separate ionization class subsets of compounds to compounds of a homogeneous ionization class reveals in this case that local models and libraries have no advantage over global models and libraries.     

Measuring experimental cyclohexane-water distribution coefficients for the SAMPL5 challenge

Small molecule distribution coefficients between immiscible nonaqueuous and aqueous phases—such as cyclohexane and water—measure the degree to which small molecules prefer one phase over another at a given pH. As distribution coefficients capture both thermodynamic effects (the free energy of transfer between phases) and chemical effects (protonation state and tautomer effects in aqueous solution), they provide an exacting test of the thermodynamic and chemical accuracy of physical models without the long correlation times inherent to the prediction of more complex properties of relevance to drug discovery, such as protein-ligand binding affinities. For the SAMPL5 challenge, we carried out a blind prediction exercise in which participants were tasked with the prediction of distribution coefficients to assess its potential as a new route for the evaluation and systematic improvement of predictive physical models. These measurements are typically performed for octanol-water, but we opted to utilize cyclohexane for the nonpolar phase. Cyclohexane was suggested to avoid issues with the high water content and persistent heterogeneous structure of water-saturated octanol phases, since it has greatly reduced water content and a homogeneous liquid structure. Using a modified shake-flask LC-MS/MS protocol, we collected cyclohexane/water distribution coefficients for a set of 53 druglike compounds at pH 7.4. These measurements were used as the basis for the SAMPL5 Distribution Coefficient Challenge, where 18 research groups predicted these measurements before the experimental values reported here were released. In this work, we describe the experimental protocol we utilized for measurement of cyclohexane-water distribution coefficients, report the measured data, propose a new bootstrap-based data analysis procedure to incorporate multiple sources of experimental error, and provide insights to help guide future iterations of this valuable exercise in predictive modeling.     

A new approach to predict the biological activity of molecules based on similarity of their interaction fields and the logP and logD values: application to auxins

The activity of a biological compound is dependent both on specific binding to a target receptor and its ADME (Absorption, Distribution, Metabolism, Excretion) properties. A challenge to predict biological activity is to consider both contributions simultaneously in deriving quantitative models. We present a novel approach to derive QSAR models combining similarity analysis of molecular interaction fields (MIFs) with prediction of logP and/or logD. This new classification method is applied to a set of about 100 compounds related to the auxin plant hormone. The classification based on similarity of their interaction fields is more successful for the indole than the phenoxy compounds. The classification of the phenoxy compounds is however improved by taking into account the influence of the logP and/or the logD values on biological activity. With the new combined method, the majority (8 out of 10) of the previously misclassified derivatives of phenoxy acetic acid are classified in accord with their bioassays. The recently determined crystal structure of the auxin-binding protein 1 (ABP1) enabled validation of our approach. The results of docking a few auxin related compounds with different biological activity to ABP1 correlate well with the classification based on similarity of MIFs only. Biological activity is, however, better predicted by a combined similarity of MIFs + logP/logD approach.     

Statistically based reduced representation of amino acid side chains

Preferred conformations of amino acid side chains have been well established through statistically obtained rotamer libraries. Typically, these provide bond torsion angles allowing a side chain to be traced atom by atom. In cases where it is desirable to reduce the complexity of a protein representation or prediction, fixing all side-chain atoms may prove unwieldy. Therefore, we introduce a general parametrization to allow positions of representative atoms (in the present study, these are terminal atoms) to be predicted directly given backbone atom coordinates. Using a large, culled data set of amino acid residues from high-resolution protein crystal structures, anywhere from 1 to 7 preferred conformations were observed for each terminal atom of the non-glycine residues. Side-chain length from the backbone C(alpha) is one of the parameters determined for each conformation, which should itself be useful. Prediction of terminal atoms was then carried out for a second, nonredundant set of protein structures to validate the data set. Using four simple probabilistic approaches, the Monte Carlo style prediction of terminal atom locations given only backbone coordinates produced an average root mean-square deviation (RMSD) of approximately 3 A from the experimentally determined terminal atom positions. With prediction using conditional probabilities based on the side-chain chi(1) rotamer, this average RMSD was improved to 1.74 A. The observed terminal atom conformations therefore provide reasonable and potentially highly accurate representations of side-chain conformation, offering a viable alternative to existing all-atom rotamers for any case where reduction in protein model complexity, or in the amount of data to be handled, is desired. One application of this representation with strong potential is the prediction of charge density in proteins. This would likely be especially valuable on protein surfaces, where side chains are much less likely to be fixed in single rotamers. Prediction of ensembles of structures provides a method to determine the probability density of charge and atom location; such a prediction is demonstrated graphically.     

A pose prediction approach based on ligand 3D shape similarity

Molecular docking predicts the best pose of a ligand in the target protein binding site by sampling and scoring numerous conformations and orientations of the ligand. Failures in pose prediction are often due to either insufficient sampling or scoring function errors. To improve the accuracy of pose prediction by tackling the sampling problem, we have developed a method of pose prediction using shape similarity. It first places a ligand conformation of the highest 3D shape similarity with known crystal structure ligands into protein binding site and then refines the pose by repacking the side-chains and performing energy minimization with a Monte Carlo algorithm. We have assessed our method utilizing CSARdock 2012 and 2014 benchmark exercise datasets consisting of co-crystal structures from eight proteins. Our results revealed that ligand 3D shape similarity could substitute conformational and orientational sampling if at least one suitable co-crystal structure is available. Our method identified poses within 2 Å RMSD as the top-ranking pose for 85.7 % of the test cases. The median RMSD for our pose prediction method was found to be 0.81 Å and was better than methods performing extensive conformational and orientational sampling within target protein binding sites. Furthermore, our method was better than similar methods utilizing ligand 3D shape similarity for pose prediction.     

Correlation and prediction of liquid-liquid distribution coefficients in aqueous systems using a modified Wilson model

A modified Wilson model is tested for its ability to correlate and predict distribution coefficients in two representative systems: 1-butanol-water and cyclohexanewater. The model is fitted to ternary equilibrium data for various solutes in these systems using a procedure involving minimization of the least-squares distance between calculated and experimental logarithmic distribution ratios. In addition, benzene-water, hexane-water, and cyclohexane-water distribution coefficients for infinitely diluted liquid solutes are predicted using only binary system information. All computations involve using both van der Waals and molar volumes as structural parameters to account for the geometry of the molecules studied. Satisfactory representations of experimental distribution ratios and fairly accurate distribution coefficients at infinite dilution are obtained for both systems. However, in a number of cyclohexane-water systems, miscibilities of constituent binary mixtures are poorly predicted from ternary system information when van der Waals volumes are used. Replacement of van der Waals volumes by molar volumes has little influence on the fit, but significant improvement is observed for the prediction of both binary miscibility properties and for distribution coefficients at infinite dilution in all the solvent-water systems considered.Presentation to First International Symposium on Solubility Phenomena, University of Western Ontario, London, Ontario, August 21–23, 1984.     

Replica exchange Monte-Carlo simulations of helix bundle membrane proteins: rotational parameters of helices

We propose a united-residue model of membrane proteins to investigate the structures of helix bundle membrane proteins (HBMPs) using coarse-grained (CG) replica exchange Monte-Carlo (REMC) simulations. To demonstrate the method, it is used to identify the ground state of HBMPs in a CG model, including bacteriorhodopsin (BR), halorhodopsin (HR), and their subdomains. The rotational parameters of transmembrane helices (TMHs) are extracted directly from the simulations, which can be compared with their experimental measurements from site-directed dichroism. In particular, the effects of amphiphilic interaction among the surfaces of TMHs on the rotational angles of helices are discussed. The proposed CG model gives a reasonably good structure prediction of HBMPs, as well as a clear physical picture for the packing, tilting, orientation, and rotation of TMHs. The root mean square deviation (RMSD) in coordinates of C_α atoms of the ground state CG structure from the X-ray structure is 5.03 Å for BR and 6.70 Å for HR. The final structure of HBMPs is obtained from the all-atom molecular dynamics simulations by refining the predicted CG structure, whose RMSD is 4.38 Å for BR and 5.70 Å for HR.     

Protein side chain modeling with orientation-dependent atomic force fields derived by series expansions

We describe the development of new force fields for protein side chain modeling called optimized side chain atomic energy (OSCAR). The distance‐dependent energy functions (OSCAR‐d) and side‐chain dihedral angle potential energy functions were represented as power and Fourier series, respectively. The resulting 802 adjustable parameters were optimized by discriminating the native side chain conformations from non‐native conformations, using a training set of 12,000 side chains for each residue type. In the course of optimization, for every residue, its side chain was replaced by varying rotamers, whereas conformations for all other residues were kept as they appeared in the crystal structure. Then, the OSCAR‐d were multiplied by an orientation‐dependent function to yield OSCAR‐o. A total of 1087 parameters of the orientation‐dependent energy functions (OSCAR‐o) were optimized by maximizing the energy gap between the native conformation and subrotamers calculated as low energy by OSCAR‐d. When OSCAR‐o with optimized parameters were used to model side chain conformations simultaneously for 218 recently released protein structures, the prediction accuracies were 88.8% for χ₁, 79.7% for χ_{1 + 2}, 1.24 Å overall root mean square deviation (RMSD), and 0.62 Å RMSD for core residues, respectively, compared with the next‐best performing side‐chain modeling program which achieved 86.6% for χ₁, 75.7% for χ_{1 + 2}, 1.40 Å overall RMSD, and 0.86 Å RMSD for core residues, respectively. The continuous energy functions obtained in this study are suitable for gradient‐based optimization techniques for protein structure refinement. A program with built‐in OSCAR for protein side chain prediction is available for download at http://sysimm.ifrec.osaka‐u.ac.jp/OSCAR/ . © 2011 Wiley Periodicals, Inc. J Comput Chem 2011     

Blind prediction of distribution in the SAMPL5 challenge with QM based protomer and p<Emphasis Type="Italic">K</Emphasis><Subscript>a</Subscript> corrections

The computation of distribution coefficients between polar and apolar phases requires both an accurate characterization of transfer free energies between phases and proper accounting of ionization and protomerization. We present a protocol for accurately predicting partition coefficients between two immiscible phases, and then apply it to 53 drug-like molecules in the SAMPL5 blind prediction challenge. Our results combine implicit solvent QM calculations with classical MD simulations using the non-Boltzmann Bennett free energy estimator. The OLYP/DZP/SMD method yields predictions that have a small deviation from experiment (RMSD = 2.3 \(\log\) D units), relative to other participants in the challenge. Our free energy corrections based on QM protomer and \({\text{p}}K_{\text{a}}\) calculations increase the correlation between predicted and experimental distribution coefficients, for all methods used. Unfortunately, these corrections are overly hydrophilic, and fail to account for additional effects such as aggregation, water dragging and the presence of polar impurities in the apolar phase. We show that, although expensive, QM-NBB free energy calculations offer an accurate and robust method that is superior to standard MM and QM techniques alone.     

Blind prediction of host–guest binding affinities: a new SAMPL3 challenge

The computational prediction of protein-ligand binding affinities is of central interest in early-stage drug-discovery, and there is a widely recognized need for improved methods. Low molecular weight receptors and their ligands--i.e., host-guest systems--represent valuable test-beds for such affinity prediction methods, because their small size makes for fast calculations and relatively facile numerical convergence. The SAMPL3 community exercise included the first ever blind prediction challenge for host-guest binding affinities, through the incorporation of 11 new host-guest complexes. Ten participating research groups addressed this challenge with a variety of approaches. Statistical assessment indicates that, although most methods performed well at predicting some general trends in binding affinity, overall accuracy was not high, as all the methods suffered from either poor correlation or high RMS errors or both. There was no clear advantage in using explicit versus implicit solvent models, any particular force field, or any particular approach to conformational sampling. In a few cases, predictions using very similar energy models but different sampling and/or free-energy methods resulted in significantly different results. The protonation states of one host and some guest molecules emerged as key uncertainties beyond the choice of computational approach. The present results have implications for methods development and future blind prediction exercises.     

Prediction of the interaction of metallic moieties with proteins: An update for protein‐ligand docking techniques

In this article, we present a new approach to expand the range of application of protein‐ligand docking methods in the prediction of the interaction of coordination complexes (i.e., metallodrugs, natural and artificial cofactors, etc.) with proteins. To do so, we assume that, from a pure computational point of view, hydrogen bond functions could be an adequate model for the coordination bonds as both share directionality and polarity aspects. In this model, docking of metalloligands can be performed without using any geometrical constraints or energy restraints. The hard work consists in generating the convenient atom types and scoring functions. To test this approach, we applied our model to 39 high‐quality X‐ray structures with transition and main group metal complexes bound via a unique coordination bond to a protein. This concept was implemented in the protein‐ligand docking program GOLD. The results are in very good agreement with the experimental structures: the percentage for which the RMSD of the simulated pose is smaller than the X‐ray spectra resolution is 92.3% and the mean value of RMSD is < 1.0 Å. Such results also show the viability of the method to predict metal complexes–proteins interactions when the X‐ray structure is not available. This work could be the first step for novel applicability of docking techniques in medicinal and bioinorganic chemistry and appears generalizable enough to be implemented in most protein‐ligand docking programs nowadays available. © 2017 Wiley Periodicals, Inc.     

D3R Grand Challenge 4: ligand similarity and MM-GBSA-based pose prediction and affinity ranking for BACE-1 inhibitors

The Drug Design Data Resource (D3R) Grand Challenges present an opportunity to assess, in the context of a blind predictive challenge, the accuracy and the limits of tools and methodologies designed to help guide pharmaceutical drug discovery projects. Here, we report the results of our participation in the D3R Grand Challenge 4 (GC4), which focused on predicting the binding poses and affinity ranking for compounds targeting the $$\beta$$-amyloid precursor protein (BACE-1). Our ligand similarity-based protocol using HYBRID (OpenEye Scientific Software) successfully identified poses close to the native binding mode for most of the ligands with less than 2 Å RMSD accuracy. Furthermore, we compared the performance of our HYBRID-based approach to that of AutoDock Vina and DOCK 6 and found that using a reference ligand to guide the docking process is a better strategy for pose prediction and helped HYBRID to perform better here. We also conducted end-point free energy estimates on molecules dynamics based ensembles of protein-ligand complexes using molecular mechanics combined with generalized Born surface area method (MM-GBSA). We found that the binding affinity ranking based on MM-GBSA scores have poor correlation with the experimental values. Finally, the main lessons from our participation in D3R GC4 are: (i) the generation of the macrocyclic conformers is a key step for successful pose prediction, (ii) the protonation states of the BACE-1 binding site should be treated carefully, (iii) the MM-GBSA method could not discriminate well between different predicted binding poses, and (iv) the MM-GBSA method does not perform well at predicting protein–ligand binding affinities here.     

Comparison of the umbrella sampling and the double decoupling method in binding free energy predictions for SAMPL6 octa-acid host–guest challenges

We calculate the absolute binding free energies of tetra-methylated octa-acids host–guest systems as a part of the SAMPL6 blind challenge (receipt ID vq30p). We employed two different free energy simulation methods, i.e., the umbrella sampling (US) and double decoupling method (DDM). The US method was used with the weighted histogram analysis method (WHAM) (US-WHAM scheme). In the DDM scheme, Hamiltonian replica-exchange method (HREM) was combined with the Bennett acceptance ratio (BAR) (HREM-BAR scheme). We obtained initial binding poses via molecular docking using GalaxyDock-HG program, which is developed for the SAMPL challenge. The root mean square deviation (RMSD) and the mean absolute deviations (MAD) using US-WHAM scheme were 1.33 and 1.02 kcal/mol, respectively. The MAD was the top among all submissions, however the correlation with respect to experiment was unexceptional. While the RMSD and MAD via HREM-BAR scheme were greater than US-WHAM scheme, (i.e., 2.09 and 1.76 kcal/mol), their correlations were slightly better than US-WHAM. The correlation between the two methods was high. Further discussion on the DDM method can be found in a companion paper by Han et al. (receipt ID 3z83m) in the same issue.     

Efficient computation of root mean square deviations under rigid transformations

The computation of root mean square deviations (RMSD) is an important step in many bioinformatics applications. If approached naively, each RMSD computation takes time linear in the number of atoms. In addition, a careful implementation is required to achieve numerical stability, which further increases runtimes. In practice, the structural variations under consideration are often induced by rigid transformations of the protein, or are at least dominated by a rigid component. In this work, we show how RMSD values resulting from rigid transformations can be computed in constant time from the protein's covariance matrix, which can be precomputed in linear time. As a typical application scenario is protein clustering, we will also show how the Ward‐distance which is popular in this field can be reduced to RMSD evaluations, yielding a constant time approach for their computation. © 2014 Wiley Periodicals, Inc.     

Fast prediction of hydration free energies for SAMPL4 blind test from a classical density functional theory

We report the performance of a classical density functional theory (CDFT) in the competition for the solvation free-energy category of the SAMPL4 blind prediction event. The theoretical calculations were carried out with the TIP3P water model and different combinations of solute configurations and molecular force fields. In comparison with the experimental data, the blind test yields an average unsigned error of 2.38 kcal/mol and the root mean square deviation of 2.99 kcal/mol. Whereas these numbers are significantly larger than the best results from explicit-solvent MD simulations, we find that the theoretical performance is sensitive to both the molecular force fields and solute configurations and that a comparable level of accuracy can be achieved by a judicious selection of the solute configurations and the force-field parameters. Most importantly, CDFT reduces the computational cost of MD simulation by almost 3 orders of magnitude, making it very attractive for large-scale hydration free-energy calculations (e.g., screening the aqueous solubility of drug-like molecules).