Docking of small molecules to farnesoid X receptors using AutoDock Vina with the Convex-PL potential: lessons learned from D3R Grand Challenge 2

The 2016 D3R Grand Challenge 2 provided an opportunity to test multiple protein–ligand docking protocols on a set of ligands bound to farnesoid X receptor that has many available experimental structures. We participated in the Stage 1 of the Challenge devoted to the docking pose predictions, with the mean RMSD value of our submission poses of 2.9 Å. Here we present a thorough analysis of our docking predictions made with AutoDock Vina and the Convex-PL rescoring potential by reproducing our submission protocol and running a series of additional molecular docking experiments. We conclude that a correct receptor structure, or more precisely, the structure of the binding pocket, plays the crucial role in the success of our docking studies. We have also noticed the important role of a local ligand geometry, which seems to be not well discussed in literature. We succeed to improve our results up to the mean RMSD value of 2.15–2.33 Å  dependent on the models of the ligands, if docking these to all available homologous receptors. Overall, for docking of ligands of diverse chemical series we suggest to perform docking of each of the ligands to a set of multiple receptors that are homologous to the target.     

Computational fragment-based screening using RosettaLigand: the SAMPL3 challenge

SAMPL3 fragment based virtual screening challenge provides a valuable opportunity for researchers to test their programs, methods and screening protocols in a blind testing environment. We participated in SAMPL3 challenge and evaluated our virtual fragment screening protocol, which involves RosettaLigand as the core component by screening a 500 fragments Maybridge library against bovine pancreatic trypsin. Our study reaffirmed that the real test for any virtual screening approach would be in a blind testing environment. The analyses presented in this paper also showed that virtual screening performance can be improved, if a set of known active compounds is available and parameters and methods that yield better enrichment are selected. Our study also highlighted that to achieve accurate orientation and conformation of ligands within a binding site, selecting an appropriate method to calculate partial charges is important. Another finding is that using multiple receptor ensembles in docking does not always yield better enrichment than individual receptors. On the basis of our results and retrospective analyses from SAMPL3 fragment screening challenge we anticipate that chances of success in a fragment screening process could be increased significantly with careful selection of receptor structures, protein flexibility, sufficient conformational sampling within binding pocket and accurate assignment of ligand and protein partial charges.     

Surflex-Dock: Docking benchmarks and real-world application

Benchmarks for molecular docking have historically focused on re-docking the cognate ligand of a well-determined protein-ligand complex to measure geometric pose prediction accuracy, and measurement of virtual screening performance has been focused on increasingly large and diverse sets of target protein structures, cognate ligands, and various types of decoy sets. Here, pose prediction is reported on the Astex Diverse set of 85 protein ligand complexes, and virtual screening performance is reported on the DUD set of 40 protein targets. In both cases, prepared structures of targets and ligands were provided by symposium organizers. The re-prepared data sets yielded results not significantly different than previous reports of Surflex-Dock on the two benchmarks. Minor changes to protein coordinates resulting from complex pre-optimization had large effects on observed performance, highlighting the limitations of cognate ligand re-docking for pose prediction assessment. Docking protocols developed for cross-docking, which address protein flexibility and produce discrete families of predicted poses, produced substantially better performance for pose prediction. Performance on virtual screening performance was shown to benefit by employing and combining multiple screening methods: docking, 2D molecular similarity, and 3D molecular similarity. In addition, use of multiple protein conformations significantly improved screening enrichment.     

Docking pose selection by interaction pattern graph similarity: application to the D3R grand challenge 2015

High affinity ligands for a given target tend to share key molecular interactions with important anchoring amino acids and therefore often present quite conserved interaction patterns. This simple concept was formalized in a topological knowledge-based scoring function (GRIM) for selecting the most appropriate docking poses from previously X-rayed interaction patterns. GRIM first converts protein–ligand atomic coordinates (docking poses) into a simple 3D graph describing the corresponding interaction pattern. In a second step, proposed graphs are compared to that found from template structures in the Protein Data Bank. Last, all docking poses are rescored according to an empirical score (GRIMscore) accounting for overlap of maximum common subgraphs. Taking the opportunity of the public D3R Grand Challenge 2015, GRIM was used to rescore docking poses for 36 ligands (6 HSP90α inhibitors, 30 MAP4K4 inhibitors) prior to the release of the corresponding protein–ligand X-ray structures. When applied to the HSP90α dataset, for which many protein–ligand X-ray structures are already available, GRIM provided very high quality solutions (mean rmsd = 1.06 Å, n = 6) as top-ranked poses, and significantly outperformed a state-of-the-art scoring function. In the case of MAP4K4 inhibitors, for which preexisting 3D knowledge is scarce and chemical diversity is much larger, the accuracy of GRIM poses decays (mean rmsd = 3.18 Å, n = 30) although GRIM still outperforms an energy-based scoring function. GRIM rescoring appears to be quite robust with comparison to the other approaches competing for the same challenge (42 submissions for the HSP90 dataset, 27 for the MAP4K4 dataset) as it ranked 3rd and 2nd respectively, for the two investigated datasets. The rescoring method is quite simple to implement, independent on a docking engine, and applicable to any target for which at least one holo X-ray structure is available.     

Protein-ligand docking against non-native protein conformers

In the validation of protein-ligand docking protocols, performance is mostly measured against native protein conformers, i.e. each ligand is docked into the protein conformation from the structure that contained that ligand. In real-life applications, however, ligands are docked against non-native conformations of the protein, i.e. the apo structure or a structure of a different protein-ligand complex. Here, we have constructed an extensive test set for assessing docking performance against non-native protein conformations. This new test set is built on the Astex Diverse Set (which we recently constructed for assessing native docking performance) and contains 1112 non-native structures for 65 drug targets. Using the protein-ligand docking program GOLD, the Astex Diverse Set and the new Astex Non-native Set, we established that, whereas docking performance (top-ranked solution within 2 A rmsd of the experimental binding mode) is approximately 80% for native docking, this drops to 61% for non-native docking. A similar drop-off is observed for sampling performance (any solution within 2 A): 91% for native docking vs 72% for non-native docking. No significant differences were observed between docking performance against apo and nonapo structures. We found that, whereas small variations in protein conformation are generally tolerated by our rigid docking protocol, larger protein movements result in a catastrophic drop-off in performance. Some docking performance and nearly all sampling performance can be recovered by considering dockings produced against a small number of non-native structures simultaneously. Docking against non-native structures of complexes containing ligands that are similar to the docked ligand also significantly improves both docking performance and sampling performance.     

Potential and limitations of ensemble docking

A major problem in structure-based virtual screening applications is the appropriate selection of a single or even multiple protein structures to be used in the virtual screening process. A priori it is unknown which protein structure(s) will perform best in a virtual screening experiment. We investigated the performance of ensemble docking, as a function of ensemble size, for eight targets of pharmaceutical interest. Starting from single protein structure docking results, for each ensemble size up to 500,000 combinations of protein structures were generated, and, for each ensemble, pose prediction and virtual screening results were derived. Comparison of single to multiple protein structure results suggests improvements when looking at the performance of the worst and the average over all single protein structures to the performance of the worst and average over all protein ensembles of size two or greater, respectively. We identified several key factors affecting ensemble docking performance, including the sampling accuracy of the docking algorithm, the choice of the scoring function, and the similarity of database ligands to the cocrystallized ligands of ligand-bound protein structures in an ensemble. Due to these factors, the prospective selection of optimum ensembles is a challenging task, shown by a reassessment of published ensemble selection protocols.     

Considerations in compound database preparation--"hidden" impact on virtual screening results

Structure-based virtual screening (SBVS) utilizing docking algorithms has become an essential tool in the drug discovery process, and significant progress has been made in successfully applying the technique to a wide range of receptor targets. In silico validation of virtual screening protocols before application to a receptor target using a corporate or commercially available compound collection is key to establishing a successful process. Ultimately, retrieval of a set of active compounds from a database of inactives is required, and the metric of enrichment (E) is habitually used to discern the quality of separation of the two. Numerous reports have addressed the performance of docking algorithms with regard to the quality of binding mode prediction and the issue of postprocessing "hit lists" of docked ligands. However, the impact of ligand database preprocessing has yet to be examined in the context of virtual screening and prioritization of compounds for biological evaluation. We provide an insight into the implications of cheminformatic preprocessing of a validation database of compounds where multiple protonated, tautomeric, stereochemical, and conformational states have been enumerated. Several commonly used methods for the generation of ligand conformations and conformational ensembles are examined, paired with an exhaustive rigid-body algorithm for the docking of different "multimeric" compound representations to the ligand binding site of the human estrogen receptor alpha. Chemgauss, a shapegaussian scoring function with intrinsic chemical knowledge, was combined with PLP as a consensus-scoring scheme to rank output from the docking protocol and enrichment rates calculated for each screen. The overheads of CPU consumption and the effect on relative database size (disk requirement) for each of the protocols employed are considered. Assessment of these parameters indicates that SBVS enrichments are highly dependent on the initial cheminformatic treatment(s) used in database construction. The interplay of SMILES representations, stereochemical information, protonation state enumeration, and ligand conformation ensembles are critical in achieving optimum enrichment rates in such screening.     

Comparing ligand interactions with multiple receptors via serial docking

Standard uses of ligand-receptor docking typically focus on the association of candidate ligands with a single targeted receptor, but actual applications increasingly require comparisons across multiple receptors. This study demonstrates that comparative docking to multiple receptors can help to select homology models for virtual compound screening and to discover ligands that bind to one set of receptors but not to another, potentially similar, set. A serial docking algorithm is furthermore described that reduces the computational costs of such calculations by testing compounds against a series of receptor structures and discarding a compound as soon as it fails to satisfy specified bind/no bind criteria for each receptor. The algorithm also realizes substantial efficiencies by taking advantage of the fact that a ligand typically binds in similar conformations to similar receptors. Thus, once detailed docking has been used to fit a ligand into the first of a series of similar receptors, much less extensive calculations can be used for the remaining structures.     

Multiple ligand docking by Glide: implications for virtual second-site screening

Performance of Glide was evaluated in a sequential multiple ligand docking paradigm predicting the binding modes of 129 protein-ligand complexes crystallized with clusters of 2-6 cooperative ligands. Three sampling protocols (single precision-SP, extra precision-XP, and SP without scaling ligand atom radii-SP hard) combined with three different scoring functions (GlideScore, Emodel and Glide Energy) were tested. The effects of ligand number, docking order and druglikeness of ligands and closeness of the binding site were investigated. On average 36?% of all structures were reproduced with RMSDs lower than 2??. Correctly docked structures reached 50?% when docking druglike ligands into closed binding sites by the SP hard protocol. Cooperative binding to metabolic and transport proteins can dramatically alter pharmacokinetic parameters of drugs. Analyzing the cytochrome P450 subset the SP hard protocol with Emodel ranking reproduced two-thirds of the structures well. Multiple ligand binding is also exploited by the fragment linking approach in lead discovery settings. The HSP90 subset from real life fragment optimization programs revealed that Glide is able to reproduce the positions of multiple bound fragments if conserved water molecules are considered. These case studies assess the utility of Glide in sequential multiple docking applications.     

Distilling the essential features of a protein surface for improving protein-ligand docking,scoring, and virtual screening

For the successful identification and docking of new ligands to a protein target by virtual screening, the essential features of the protein and ligand surfaces must be captured and distilled in an efficient representation. Since the running time for docking increases exponentially with the number of points representing the protein and each ligand candidate, it is important to place these points where the best interactions can be made between the protein and the ligand. This definition of favorable points of interaction can also guide protein structure-based ligand design, which typically focuses on which chemical groups provide the most energetically favorable contacts. In this paper, we present an alternative method of protein template and ligand interaction point design that identifies the most favorable points for making hydrophobic and hydrogen–bond interactions by using a knowledge base. The knowledge-based protein and ligand representations have been incorporated in version 2.0 of SLIDE and resulted in dockings closer to the crystal structure orientations when screening a set of 57 known thrombin and glutathione S–transferase (GST) ligands against the apo structures of these proteins. There was also improved scoring enrichment of the dockings, meaning better differentiation between the chemically diverse known ligands and a 15,000-molecule dataset of randomly-chosen small organic molecules. This approach for identifying the most important points of interaction between proteins and their ligands can equally well be used in other docking and design techniques. While much recent effort has focused on improving scoring functions for protein-ligand docking, our results indicate that improving the representation of the chemistry of proteins and their ligands is another avenue that can lead to significant improvements in the identification, docking, and scoring of ligands.(These authors contributed equally to this work)     

Scoring functions for protein-ligand docking

Virtual screening by molecular docking has become established as a method for drug lead discovery and optimization. All docking algorithms make use of a scoring function in combination with a method of search. Two theoretical aspects of scoring function performance dominate operational performance. The first is the degree to which a scoring function has a global extremum within the ligand pose landscape at the proper location. The second is the degree to which the magnitude of the function at the extremum is accurate. Presuming adequate search strategies, a scoring function's location performance will dominate behavior with respect to docking accuracy: the degree to which a predicted pose of a ligand matches experimental observation. A scoring function's magnitude performance will dominate behavior with respect to screening utility: enrichment of true ligands over non-ligands. Magnitude estimation also controls pure scoring accuracy: the degree to which bona fide ligands of a particular protein may be correctly ranked. Approaches to the development of scoring functions have varied widely, with a number of functions yielding similarly high levels of performance relating to the location issue. However, even among functions performing equally well on location, widely varying performance is observed on the question of magnitude. In many cases, performance is good enough to yield high enrichments of true ligands versus non-ligands in screening across a wide variety of protein types. Generally, performance is not good enough to correctly rank among true ligands. Strategies for improvement are discussed.     

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.     

Virtual screening of the SAMPL4 blinded HIV integrase inhibitors dataset

Several combinations of docking software and scoring functions were evaluated for their ability to predict the binding of a dataset of potential HIV integrase inhibitors. We found that different docking software were appropriate for each one of the three binding sites considered (LEDGF, Y3 and fragment sites), and the most suitable two docking protocols, involving Glide SP and Gold ChemScore, were selected using a training set of compounds identified from the structural data available. These protocols could successfully predict respectively 20.0 and 23.6 % of the HIV integrase binders, all of them being present in the LEDGF site. When a different analysis of the results was carried out by removing all alternate isomers of binders from the set, our predictions were dramatically improved, with an overall ROC AUC of 0.73 and enrichment factor at 10 % of 2.89 for the prediction obtained using Gold ChemScore. This study highlighted the ability of the selected docking protocols to correctly position in most cases the ortho-alkoxy-carboxylate core functional group of the ligands in the corresponding binding site, but also their difficulties to correctly rank the docking poses.     

Rescoring of docking poses under Occam’s Razor: are there simpler solutions?

Ligand affinity prediction from docking simulations is usually performed by means of highly empirical and diverse protocols. These protocols often involve the re-scoring of poses generated by a force field (FF) based Hamiltonian to provide either estimated binding affinities—or alternatively, some empirical goodness score. Re-scoring is performed by so-called scoring functions—typically, a reweighted sum of FF terms augmented by additional terms (e.g., desolvation/entropic penalty, hydrophobicity, aromatic interactions etc.). Sometimes, the scoring function actually drives ligand positioning, but often it only operates on the best scoring poses ranked top by the initial ligand positioning tool. In either of these rather intricate scenarios, scoring functions are docking-specific models, and most require machine-learning-based calibration. Therefore, docking simulations are less straightforward when compared to “standard” molecular simulations in which the FF Hamiltonian defines the energy, and affinity emerges as an ensemble average property over pools of representative conformers (i.e., the trajectory). Paraphrasing on Occam’s Razor principle, additional model complexity is only acceptable if demonstrated to bring a significant improvement of prediction quality. In this work we therefore examined whether the complexity inherent to scoring functions is indeed justified. For this purpose we compared sampler for multiple protein–ligand entities, a general purpose conformation sampler based on the AMBER/GAFF FF, complemented with continuum solvation terms, with several state of the art docking tools that rely on calibrated scoring functions (Glide, Gold, Autodock-Vina) in terms of its ability to top-rank the actives from large and diverse ligand series associated with various proteins. There is no clear winner of this study, where each program performed well on most of the targets, but also failed with respect to at least one of them. Therefore, a well-parameterized force field with a simple, energy-based ligand ranking protocol appears to be an as effective docking protocol as intricate rescoring strategies based on scoring functions. A tool that can sample the conformational space of the free ligand, the bound ligand and the protein binding site using the same force field may avoid many of the approximations common to contemporary docking protocols and allow e.g., for docking into highly flexible active sites, when current scoring functions are not well suited to estimate receptor strain energies.     

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.     

Fully automated molecular mechanics based induced fit protein-ligand docking method

We describe a method for docking a ligand into a protein receptor while allowing flexibility of the protein binding site. The method employs a multistep procedure that begins with the generation of protein and ligand conformations. An initial placement of the ligand is then performed by computing binding site hotspots. This initial placement is followed by a protein side-chain refinement stage that models protein flexibility. The final step of the process is an energy minimization of the ligand pose in the presence of the rigid receptor. Thus the algorithm models flexibility of the protein at two stages, before and after ligand placement. We validated this method by performing docking and cross docking studies of eight protein systems for which crystal structures were available for at least two bound ligands. The resulting rmsd values of the 21 docked protein-ligand complexes showed values of 2 A or less for all but one of the systems examined. The method has two critical benefits for high throughput virtual screening studies. First, no user intervention is required in the docking once the initial binding site selection has been made in the protein. Second, the initial protein conformation generation needs to be performed only once for a given binding region. Also, the method may be customized in various ways depending on the particular scenario in which dockings are being performed. Each of the individual steps of the method is fully independent making it straightforward to explore different variants of the high level workflow to further improve accuracy and performance.     

Lead Finder docking and virtual screening evaluation with Astex and DUD test sets

Lead Finder is a molecular docking software. Sampling uses an original implementation of the genetic algorithm that involves a number of additional optimization procedures. Lead Finder's scoring functions employ a set of semi-empiric molecular mechanics functionals that have been parameterized independently for docking, binding energy predictions and rank-ordering for virtual screening. Sampling and scoring both utilize a staged approach, moving from fast but less accurate algorithm versions to computationally more intensive but more accurate versions. Lead Finder includes tools for the preparation of full atom protein and ligand models. In this exercise, Lead Finder achieved 72.9% docking success rate on the Astex test set when the original author-prepared full atom models were used, and 74.1% success rate when the structures were prepared by Lead Finder. The major cause of docking failures were scoring errors resulting from the use of imperfect solvation models. In many cases, docking errors could be corrected by the proper protonation and the use of correct cyclic conformations of ligands. In virtual screening experiments on the DUD test set the early enrichment factor of several tens was achieved on average. However, the area under the ROC curve ("AUC ROC") ranged from 0.70 to 0.74 depending on the screening protocol used, and the separation from the null model was not perfect-0.12-0.15 units of AUC ROC. We assume that effective virtual screening in the whole range of enrichment curve and not just at the early enrichment stages requires more accurate solvation modeling and accounting for the protein backbone flexibility.     

Optimization of high throughput virtual screening by combining shape-matching and docking methods

Receptor flexibility is a critical issue in structure-based virtual screening methods. Although a multiple-receptor conformation docking is an efficient way to account for receptor flexibility, it is still too slow for large molecular libraries. It was reported that a fast ligand-centric, shape-based virtual screening was more consistent for hit enrichment than a typical single-receptor conformation docking. Thus, we designed a "distributed docking" method that improves virtual high throughput screening by combining a shape-matching method with a multiple-receptor conformation docking. Database compounds are classified in advance based on shape similarities to one of the crystal ligands complexed with the target protein. This classification enables us to pick the appropriate receptor conformation for a single-receptor conformation docking of a given compound, thereby avoiding time-consuming multiple docking. In particular, this approach utilizes cross-docking scores of known ligands to all available receptor structures in order to optimize the algorithm. The present virtual screening method was tested for reidentification of known PPARgamma and p38 MAP kinase active compounds. We demonstrate that this method improves the enrichment while maintaining the computation speed of a typical single-receptor conformation docking.     

Large scale free energy calculations for blind predictions of protein–ligand binding: the D3R Grand Challenge 2015

We describe binding free energy calculations in the D3R Grand Challenge 2015 for blind prediction of the binding affinities of 180 ligands to Hsp90. The present D3R challenge was built around experimental datasets involving Heat shock protein (Hsp) 90, an ATP-dependent molecular chaperone which is an important anticancer drug target. The Hsp90 ATP binding site is known to be a challenging target for accurate calculations of ligand binding affinities because of the ligand-dependent conformational changes in the binding site, the presence of ordered waters and the broad chemical diversity of ligands that can bind at this site. Our primary focus here is to distinguish binders from nonbinders. Large scale absolute binding free energy calculations that cover over 3000 protein–ligand complexes were performed using the BEDAM method starting from docked structures generated by Glide docking. Although the ligand dataset in this study resembles an intermediate to late stage lead optimization project while the BEDAM method is mainly developed for early stage virtual screening of hit molecules, the BEDAM binding free energy scoring has resulted in a moderate enrichment of ligand screening against this challenging drug target. Results show that, using a statistical mechanics based free energy method like BEDAM starting from docked poses offers better enrichment than classical docking scoring functions and rescoring methods like Prime MM-GBSA for the Hsp90 data set in this blind challenge. Importantly, among the three methods tested here, only the mean value of the BEDAM binding free energy scores is able to separate the large group of binders from the small group of nonbinders with a gap of 2.4 kcal/mol. None of the three methods that we have tested provided accurate ranking of the affinities of the 147 active compounds. We discuss the possible sources of errors in the binding free energy calculations. The study suggests that BEDAM can be used strategically to discriminate binders from nonbinders in virtual screening and to more accurately predict the ligand binding modes prior to the more computationally expensive FEP calculations of binding affinity.     

Docking to heme proteins

In silico screening has become a valuable tool in drug design, but some drug targets represent real challenges for docking algorithms. This is especially true for metalloproteins, whose interactions with ligands are difficult to parametrize. Our docking algorithm, EADock, is based on the CHARMM force field, which assures a physically sound scoring function and a good transferability to a wide range of systems, but also exhibits difficulties in case of some metalloproteins. Here, we consider the therapeutically important case of heme proteins featuring an iron core at the active site. Using a standard docking protocol, where the iron–ligand interaction is underestimated, we obtained a success rate of 28% for a test set of 50 heme‐containing complexes with iron‐ligand contact. By introducing Morse‐like metal binding potentials (MMBP), which are fitted to reproduce density functional theory calculations, we are able to increase the success rate to 62%. The remaining failures are mainly due to specific ligand–water interactions in the X‐ray structures. Testing of the MMBP on a second data set of non iron binders (14 cases) demonstrates that they do not introduce a spurious bias towards metal binding, which suggests that they may reliably be used also for cross‐docking studies. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2009