Reliability analysis and optimization of the consensus docking approach for the development of virtual screening studies

Abstract

Ligand-protein docking is one of the most common techniques used in virtual screening campaigns. Despite the large number of docking software available, there is still the need of improving the efficacy of docking-based virtual screenings. To date, only very few studies evaluated the possibility of combining the results of different docking methods to achieve higher success rates in virtual screening studies (consensus docking). In order to better understand the range of applicability of this approach, we carried out an extensive enriched database analysis using the DUD dataset. The consensus docking protocol was then refined by applying modifications concerning the calculation of pose consensus and the combination of docking methods included in the procedure. The results obtained suggest that this approach performs as well as the best available methods found in literature, confirming the idea that this procedure can be profitably used for the identification of new hit compounds.

Keywords:

• Consensus docking
• hit identification
• virtual screening

Introduction

Molecular docking has a pivotal role in the computer-aided drug design field. Besides being the prototype of the receptor-based approaches, it is probably the principal in silico technique applied to drug discovery. The widespread use of docking is due to the different purposes it can be employed for, which target different steps and issues of a drug discovery campaign. Ligand-protein docking represents the key tool to predict the most energetically favored disposition of a ligand bound to its corresponding protein target, and the predicted binding mode of hit and lead compounds can represent a valuable starting point for SAR analyzes and lead optimization studies aimed at improving ligand activity^1–3. Docking is also largely used for hit identification, representing the main example of receptor-based virtual screening (VS) procedure: a docking-based approach is commonly reported in several successful VS studies, both alone and in combination with other computational strategies^4–8. For these reasons, the exponential attention that docking is gaining over time, as proved by the number of novel docking procedures and scoring functions reported in literature only in the last few years, cannot be blamed^9–20. However, the remarkable interest for docking and for the development of new and optimized approaches is mainly due to the fact that the prediction of the actual binding mode of a ligand into a protein target is still far from being an easy task. The main problems connected to ligand pose prediction are represented by the flexibility of the ligand (and protein) as well as by the estimation of the ligand–protein interaction energy relative to the various ligand binding conformations evaluated. These two issues are treated through different approaches in the various docking sofware currently available. The sampling of the ligand conformational space within a defined receptor binding pocket is mainly performed by either stochastic (e.g. evolutionary and Monte Carlo algorithms) or systematic search methods (e.g. incremental construction, exhaustive search)^21,22. The binding energy associated to the different docking poses is calculated through the application of scoring functions (SF) that can be classified into knowledge-based, force field-based and empirical SFs^23,24. Each of these and other methods implemented in the most popular docking programs has its own strenghts and limitations, therefore there is not a particular algorithm/SF combination outperforming the other ones in terms of reliablity and accuracy. As a consequence, it is not possible to know a priori if a docking procedure is more suitable than another one for a particular task. Several attempts to increase docking accuracy and reliability, thus improving its VS performance, reported in literature concerned the development of optimized SFs, rescoring procedures and pose filtering methods^13,25–31. Among these strategies, the combination of different scoring functions (consenus scoring) to improve VS hit rates is a well known practice introduced by Charifson et al. in 1999³¹ that is still commonly applied^33–38. A very recent approach concerns the introduction of a ligand-based component within the docking protocol that consists in exploiting the ligand information relative to crystallographic ligand-protein complexes for guiding pose selection and ranking, in a sort of combined ligand/receptor-based approach^12,39–41. Surprisingly, only very few studies evaluated the possibility of combining the results of different docking methods (consensus docking) to obtain reliable ligand binding poses and to improve the success rate in VS studies, as recently reported by Houston and Walkinshaw⁴². Inspired by their analysis, we very recently assessed the reliability of a consensus docking (CD) protocol employing ten different docking procedures at the same time to improve the performance of docking from both a qualitative and a quantitative point of view⁴³. Cross-docking studies carried out on three different protein targets showed that the CD approach was able to predict the ligand binding poses better than the single docking methods and that the higher the consensus level of a docking pose (i.e. the number of docking methods yielding the same pose) the higher was its reliability. Moreover, an enriched database analysis performed using three different data sets belonging from the Directory of Useful Decoys (DUD) suggested that the CD approach could be a valuable tool for the identification of new hit compounds through VS studies. In fact, it was found that filtering the database compounds based on their consensus level and thus selecting only those for which a common binding pose was predicted by most of the docking procedures, resulted in high enrichment factors and showed to perform as well as the best available methods reported in literature. The reliability of the CD approach was also experimentally assessed through its applications in a VS campaign aimed at identifying new non-covalent fatty acid amide hydrolase (FAAH) inhibitors, where a pre-filtered library of commercial compounds was subjected to the CD protocol refined by molecular dynamic (MD) simulations, which were performed for the top-ranked compounds. The final 80% hit rate achieved and the identification of two low micromolar non-covalent FAAH inhibitors confirmed the efficacy of the CD applied to VS studied⁴⁴. Encouraged by the results obtained with both the in silico and experimental studies we decided to further test the CD approach in order to better understand the range of applicability of this procedure in VS studies and with the aim of refining the protocol in the attempt to maximize its performance and efficiency. For these purposes, we extended the enriched database analysis based on DUD data set and evaluated the results in terms of enrichment factor. Then, we applied modifications to the protocol concerning the calculation of pose consensus as well as the type and number of docking methods employed, studying the impact of these changes in the overall performance of the CD approach.

Results and discussion

The initial purpose of our study was to assess if the CD approach could be considered as a general valuable strategy to be applied in VS campaigns, as suggested by the promising results obtained through both the enriched database analysis and the experimental evaluation of CD that we recently reported. On these bases, we tested the reliability of the consensus docking virtual screening protocol by applying it on the Directory of Useful Decoys (DUD)⁴⁵. For each enriched data set, actives and decoy compounds were docked into the crystal structure of the corresponding protein target employing the ten different docking methods included in the already published CD approach⁴³, for a total of more than one million docking calculations. For each analyzed molecule, 10 different docking poses were thus obtained and their reciprocal distances in terms of root-mean-square deviation (RMSD) were then calculated in order to perform pose clustering and identify common binding modes (see Methods for details). The consensus level, i.e. the number of docking methods predicting the same binding pose, was calculated for each docked compound and it was used as a parameter to filter the enriched databases. Initially, the RMSD threshold for pose consensus was set to 2.0 Å as in our previous analysis⁴³; therefore, the consensus level corresponded to the number of ligand poses with reciprocal RMSD < 2.0 Å, which were considered as common binding modes. The enrichment factor (EF) was used as a metric for the performance assessment of the CD approach. For each enriched database tested, the EF was calculated at the different levels of pose consensus shown by the docked compounds and it was compared to the maximum EF achievable for the specific data set, which ranged from 35.0 to 37.0. In agreement with the results obtained in our previous CD enriched database analysis, by increasing the consensus level considered for the selection of compounds we obtained an increment of the corresponding EF. The performance of the CD approach was evaluated based on the highest EF achieved for each data set (EF_MAX). Overall, the CD filter showed satisfying results, since the average value of EF_MAX (aEF_MAX) obtained for the different databases was 16.7, corresponding to about 45% of the maximum EF achievable, and thus to an average 45% hit rate for the subsets of compounds selected from the whole databases. In particular, 17 out of 35 data sets showed an EF_MAX corresponding to a hit rate ranging from 50% to 100% (Figure 1).

Figure 1. Direct comparison between the EF_MAX results obtained using the CD virtual screening (dark grey) and the DDFA approach (grey). The average EF_MAX are reported as dotted lines.

In 2014, Arciniega and Lange reported the evaluation of the docking data feature analysis (DDFA) as a new intriguing approach for carrying out virtual screening experiments⁴⁶. This approach, based on artificial neural network, was tested by using DUD data sets and the analysis highlighted that it achieved results similar to those shown by the best performing methods available in literature. Since the authors reported in their analysis the EF_MAX obtained for each data set, it was possible to carry out a direct comparison between the two methods. As shown in Figure 1, although the two procedures showed a similar performance, the CD approach proved to be more effective, since the average value of EF_MAX obtained for CD (EF_MAX = 16.7) was higher than that reported for DDFA, which corresponded to 13.8 (Figure 1).

In order to investigate whether the consensus docking results were influenced by the molecular properties of the docked compounds, the EF_MAX value obtained for each of the analyzed DUD data sets was compared with the average molecular weight, number of heavy atoms and number of rotatable bonds of the corresponding database compounds. Furthermore, the same analysis was carried out taking into account the volume of the binding site cavities, thus studying the possible correlation with the binding site characteristics of the target proteins. As shown in Figure S1 (Supplementary Material), this analysis highlighted a complete lack of correlation between the EF_MAX and any of the properties considered, suggesting that the performance of CD in the filtering of enriched databases was not biased by the size and the conformational freedom of the screened compounds, or by the size of the target binding sites. This was in agreement with the results obtained in our previous consensus cross-docking analyzes, which evidenced no correlation between consensus pose level and ligand properties⁴³. Overall, these results confirmed that the CD approach represents a valuable tool that can be useful in VS studies to identify ligand hits for various protein targets and it showed to perform better than one of the best available methods reported in literature.

The second part of our study was aimed at improving the performance and efficiency of the CD docking approach. As a first step, we questioned about the possibility of achieving better VS results by using different RMSD thresholds to cluster the ten docking poses obtained for each ligand and thus to calculate their pose consensus level. A RMSD value of 2.0 Å is a well-accepted cutoff used to evaluate the success rate of docking and conformational sampling methods in reproducing experimental ligand poses, since two small-molecule conformations showing RMSD < 2.0 Å are commonly considered as equivalent. In order to study the effect of the RMSD tolerance used for pose clustering on the VS performance of CD, the docking results obtained for all the different database compounds were re-clustered using different RMSD thresholds ranging from 0.5 to 4.0 Å, with an increment of 0.1 Å, and the new corresponding pose consensus levels were calculated. The VS results relative to each new CD analysis were then evaluated and compared in terms of aEF_MAX and Figure 2 shows the distribution of the aEF_MAX values obtained as a function of the RMSD cutoff used for pose clustering. As shown in this figure, the analysis revealed that the optimal RMSD threshold corresponded to 1.5 Å, since it produced an aEF_MAX of 19.6, which was the highest aEF_MAX obtained among all the different analyzes. An increase of the RMSD cutoff with respect to the optimal value generated a proportional decrease of the aEF_MAX, which was probably due to a corresponding reduced efficacy of the clustering protocol in discriminating qualitatively different ligand binding poses. Conversely, the lower aEF_MAX values obtained with a further reduction of the RMSD threshold derived from an excessively strict clustering procedure identifying too many different binding modes for each docked compound.

Figure 2. Average EF_MAX values obtianed for the DUD dataset by modifying the RMSD cutoff used for pose clustering in the CD calculations.

Once optimized the pose clustering procedure, we investigated if it was possible to achieve an improvement in the efficiency and also a further gain in the performance of the CD virtual screening approach by changing the combination of docking methods included in the protocol. For this purpose, we first evaluated the contribution of each of the ten docking procedures to the final performance of the CD docking by analyzing the effect of their exclusion from the protocol. The docking results obtained for the different database compounds were re-clustered using the optimized RMSD threshold and considering only 9 out of the 10 original docking methods, thus removing a docking procedure from the calculation of pose clustering. Following this approach, ten different analyzes were carried out, so that in each of them a different docking method was excluded from the pose clustering and from the calculation of the consensus level. In this way, it was possible to understand the importance of each single docking procedure in the overall CD performance. The results relative to each new CD analysis were evaluated and compared in terms of aEF_MAX. As shown in Figure 3(A), the removal of GOLD-ChemPLP method from the consensus protocol was found to produce no deleterious effect on the VS performance of the CD approach. Indeed, the corresponding value of aEF_MAX observed (19.8) was slightly higher than that obtained with the full consensus protocol (19.6). Therefore, this analysis highlighted that the efficiency of the CD virtual screening could be improved by the exclusion of GOLD-ChemPLP procedure from the set of docking methods originally used in the protocol, since a considerable reduction of the computation time required to perform the CD calculations could be obtained without negatively affecting its filtering efficacy.

Figure 3. Effect on the aEF_MAX obtained by removing a docking method (A) from the initial set of ten docking procedures, (B) in addition to the exclusion of the GOLD-ChemPLP method, (C) in addition to the exclusion of the GOLD-ChemPLP and GLIDE-XP methods. The reference aEF_MAX obtained by using the ten docking procedures is reported as a dotted line.

This analysis was then recursively applied to evaluate if a further gain in efficiency could be obtained by removing other docking methods from the CD protocol. The results showed that the exclusion of GLIDE-XP method in addition to GOLD-ChemPLP was tolerated, since the corresponding aEF_MAX achieved (19.8) was identical to that obtained using the CD protocol where only GOLD-ChemPLP method was removed (Figure 3B). On the contrary, it was found that no other docking procedure could be excluded without significantly decreasing the performance of the CD protocol; in fact, all the EF analyzes carried out considering only 7 out of the 10 initial docking methods (i.e. removing GLIDE-XP, GOLD-ChemPLP and another of the 8 remaining methods) showed lower aEF_MAX values (Figure 3C). Overall, this analysis allowed to substantially improve the efficiency of the consensus docking approach through the identification of the best combination of docking methods to be employed among the ten procedures originally included in the CD approach. In fact, such CD protocol comprising only 8 docking methods (CD8) achieved the same VS performance than that obtained with the full protocol in much less computation time.

The last step of our study concerned the introduction in the consensus docking protocol of new procedures that were not considered within the original set of selected docking methods, in order to evaluate the possibility of a further improvement of the VS performance. Three different docking software were selected for this analysis: GLAMDOCK⁴⁷, PLANTS (Protein-Ligand ANT System)⁴⁸ and rDOCK⁴⁹. Three different EF evaluations were thus carried out, each employing 9 docking methods for pose clustering and consensus level calculations, i.e. the 8 procedures selected through the previous analysis and one of the three new docking methods taken into account. Table 1 shows the values of aEF_MAX achieved by including each of the new docking procedures compared with the aEF_MAX obtained by using the CD8 protocol. While the use of GLAMDOCK and rDOCK led to a reduced efficacy in the filtering of the enriched data sets, the introduction of PLANTS produced a positive effect, since it was associated with a slight increase of the aEF_MAX value, which corresponded to 20.3.

Table 1. Consensus docking results obtained by adding the reported docking procedure to the optimized CD8 protocol (i.e. the combination of the following docking procedures: AUTODOCK, DOCK, FRED, GOLD-ASP, GOLD-CSCORE, GOLD-GSCORE, GLIDE-SP, VINA).

Docking software	aEFMAX
GLAMDOCK	17.6
PLANTS	20.3
rDOCK	19.0
CD8	19.8

The results of this analysis highlighted that the consensus docking approach could be further optimized by broadening the set of docking procedures to be considered as potential members of the CD protocol and testing different combinations of docking methods in terms of both number and type of procedures employed. Overall, the study herein reported already allowed a substantial optimization of the consensus docking protocol. Figure S2 (Supplementary Material) shows the values of EF_MAX obtained by using the optimized CD9 protocol (including PLANTS) and the original one for the 35 DUD data sets herein analyzed. For nine targets, the modifications applied to the initial protocol produced an increase in the EF_MAX corresponding to more than 50% of the initial value. In particular, five targets showed a 2-fold or higher improvement in EF_MAX and the number of data sets showing an EF_MAX corresponding to a hit rate ranging from 50% to 100% passed from 17 (as obtained using the original CD protocol) to 21. On the contrary, only for very few targets a lower EF_MAX was obtained using the optimized procedure. Overall, the final aEF_MAX obtained with the CD9 protocol (20.3) corresponded to about 55% of the maximum EF achievable for the different databases, and thus to an average 55% hit rate for the subsets of selected compounds.

Conclusions

The study herein reported strongly suggests that the CD approach can represent a powerfool tool for the development of virtual screening studies, consistently with the results obtained in our already published computational and experimental validations. The extensive enriched database analysis performed using DUD data set highlighted the general reliability of the CD approach and its efficacy in identifying active ligands for a wide panel of different protein targets, showing even better results than one of the best available methods reported in literature. Moreover, by tuning the pose clustering procedure and the combination of docking methods included in the consensus docking protocol, we obtained a substantial improvement of the virtual screening performance accompanied by a reduction of the total CPU time required for calculations, which was the major limitation of this approach. Such new optimized consensus docking procedure showed very promising results in terms of enrichment factors and hit rates achieved, encouraging us to apply it in further experimental studies aimed at the identification of new ligand hits through virtual screening campaigns. However, it is worth considering that for a small number of DUD targets it was not possible to achieve satisfying results despite the optimization of the protocol. Since the performance of the consensus docking showed a complete lack of correlation with both the molecular properties of the screened compounds and the size of the target binding sites, further investigations are still necessary to identify possible elements that can limit the efficacy of this approach. On the other hand, the results of our study suggest that by applying other modifications to the protocol (i.e. evaluating new combinations of docking methods) it could be possible to obtain further improvements in the performance of consensus docking. Therefore, additional analysis employing different enriched databases, like the Maximum Unbiased Validation (MUV) data sets⁵⁰, will be carried out to address these issues and to further refine the protocol. Overall, our study suggests that the CD approach represents a powerful virtual screening strategy for the identification of new ligand hits for a broad range of different targets. Moreover, like all computational procedures, a validation test using a small enriched data set of compounds would be recommended before starting a virtual screening campaign, in order to confirm the efficacy of the protocol for the target of interest.

Methods

Docking procedures

For all docking analyzes, only the best scored pose was taken into account.

Autodock 4.2.3

AUTODOCK Tools utilities⁵¹ were used in order to identify the torsion angles in the ligands, to add the solvent model and assign the Gasteiger atomic charges to proteins and ligands. The regions of interest used by AUTODOCK⁵² were defined by considering the reference ligand as the central group of a grid box of 10 Å in the x, y and z directions. A grid spacing of 0.375 Å and a distance dependent function of the dielectric constant were used for the energetic map calculations. By using the Lamarckian genetic algorithm, the docked compounds were subjected to 20 runs of the AUTODOCK search using 2500000 steps of energy evaluation and the default values of the other parameters.

Dock 6.7

The molecular surface of the binding site was calculated by means of the MS program⁵³, generating the Connolly surface with a 1.4 Å radius probe. The points of the surface and the vectors normal to it were used by the Sphgen program in order to build a set of spheres with radii varying from 1.4 to 4 Å, which describe from a stereoelectronic point of view the negative image of the site. Spheres within a radius of 10 Å from the reference ligand were used to represent the site. For each ligand, DOCK 6.7 calculated 1000 orientations; of these, the best grid scored was taken into consideration. The grid-based score is based on the nonbonded terms of the molecular mechanic force field. The ligand charge was calculated using the AM1-BCC method, as implemented in the MOLCHARGE program⁵⁴.

Fred 3.0

FRED⁵⁵ requires a set of input conformers for each ligand. The conformers were generated by OMEGA2^56–58. We applied the following modifications to the default settings of OMEGA2: the energy window was set at 50.0, the maximum number of output conformers was set at 10000, the time limit was set at 1200, and the root mean square deviation (RMSD) value below which two conformations were considered to be similar was set at 0.3 Å. The region of interest for the docking studies was defined in such a manner that it contained all residues staying within 10 Å from the ligand in the X-ray structures. FRED default parameters were used setting the high dock_resolution.

Glide 5.0

The binding site was defined by a rectangular box of 10 Å in the x, y and z directions centered on the ligand. The option allowing only the docking of ligands containing a defined range of atoms were deactivated, so that all the ligands were docked independently from the number of their atoms, whereas GLIDE⁵⁹ defaults were used for all other parameters. Docking calculations were carried out using the standard precision (SP) and extra precision (XP) methods.

Gold 5.1

The region of interest for the docking studies was defined in such a manner that it contained all residues that stayed within 10 Å from the ligand in the X-ray structures; the “allow early termination” command was deactivated, while the possibility for the ligand to flip ring corners was activated. For all other parameters, GOLD⁶⁰ defaults were used and the ligands were subjected to 30 genetic algorithm runs. Four docking analyzes were carried out. The four fitness functions implemented in GOLD, i.e., GoldScore (GS), ChemScore (CS), Astex Statistical Potential (ASP) and ChemPLP (PLP), were used.

Autodock vina 1.1

The input files for the protein and ligands originated from the AUTODOCK Tools utilities for the AUTODOCK calculations were also used for the AUTODOCK VINA⁶¹ calculations, including the grid box dimensions. The exhaustiveness parameter was set to 10 and the Energy_range to 1, whereas for all other parameters, AUTODOCK VINA defaults were used. The VINA scoring function combines certain advantages of knowledge-based potentials and empirical scoring functions, extracting information from both the conformational preferences of the receptor-ligand complexes and the experimental affinity measurements.

Glamdock 1.0

The GLAMDOCK docking protocol consisted of 5 docking runs each comprising 650 Monte Carlo minimization (MCM) steps, with 15 steps of Levenberg–Marquardt minimization in torsion space at each MCM step. A maximum of 40 poses were finally post-minimized by 150 steps of Levenberg–Marquardt⁴⁷.

Plants

This docking software uses Ant Colony Optimization, a state-of-the-art global optimization algorithm to find minima of a scoring function representing favorable complex structures⁶². ChemPLP scoring function was employed to score protein-ligand interactions as well as intra-ligand clash terms. Standard settings for all parameters were used for the scoring function as well as the optimization algorithm (search speed setting: “speed1”). The regions of interest used by PLANTS⁴⁸ were defined by considering the reference ligand as the central group of a grid box of 10 Å in the x, y and z directions.

rDOCK 1.0

This docking software uses a combination of stochastic and deterministic search techniques to generate low energy ligand poses⁴⁹. The docking protocol generates a single ligand pose using 3 stages of Genetic Algorithm search (GA1, GA2, GA3), followed by low temperature Monte Carlo (MC) and Simplex minimization (MIN) stages. The GA stages are interdependent and are designed to be used sequentially. The cavity within a radius of 10 Å from the reference ligand was used to represent the binding site; for all the other parameters rDOCK defaults were used.

Consensus docking evaluation

To test the CD approach, we analyzed the broadly used Directory of Useful Decoys (DUD)⁴⁵. The four databases relative to targets presenting metal ion prosthetic groups (i.e. ACE, ADA, COMT, PDE5) were not taken into account in the analysis, since the presence of ions within the receptor binding site considered for docking could negatively affect the performance of some of the docking methods employed in the CD approach. Moreover, the PDGFrb target was not considered in the study because the protein structure included in the DUD database and supposed to be used for docking was not an X-ray crystal structure but only a homology model. For each data set considered in the study, the compounds were docked into their corresponding target structure by using the different docking procedures described above. The RMSD of each docking pose against the remaining docking results obtained for the same ligand was evaluated by using the rms_analysis software of the GOLD suite. On this basis, for each database ligand an N × N matrix was generated reporting the RMSD results (N = number of docking procedures). By using an in-house program, these results were clustered so that, among the N results, all the similar docking poses were clustered together according to the chosen RMSD threshold. As clustering algorithm we used the complete-linkage method which is an agglomerative type of hierarchical clustering. This method starts considering each element in a cluster of its own. The clusters are then sequentially combined into larger ones, until all elements are in the same cluster. At each step, the two clusters separated by the shortest distance are combined. The consensus docking virtual screening results were evaluated by using the enrichment factor (EF) value where t_p is the number of high affinity compounds retrieved (true positives); f_n is the number of high affinity compounds rejected during the VS filtering (false negatives); NC_tot is the total number of molecules of the database; NC is the total number of compounds obtained by the VS protocol.

Declaration of interest

The authors report no declarations of interest. We are grateful to the University of Pisa (Progetti di Ricerca di Ateneo, PRA 2016) for funding.

Supplementary materials available online only

Acknowledgements

Many thanks are due to Prof. Maurizio Botta for the use of the GLIDE program in his computational laboratory (University of Siena, Italy).