
==== Front
Sci Rep
Sci Rep
Scientific Reports
2045-2322
Nature Publishing Group UK London

39237737
71699
10.1038/s41598-024-71699-3
Article
Teaching old docks new tricks with machine learning enhanced ensemble docking
Bhatt Roshni
Wang Ann
Durrant Jacob D. durrantj@pitt.edu

https://ror.org/01an3r305 grid.21925.3d 0000 0004 1936 9000 Department of Biological Sciences, University of Pittsburgh, Pittsburgh, PA 15260 USA
5 9 2024
5 9 2024
2024
14 2072225 5 2024
30 8 2024
© The Author(s) 2024
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.
We here introduce Ensemble Optimizer (EnOpt), a machine-learning tool to improve the accuracy and interpretability of ensemble virtual screening (VS). Ensemble VS is an established method for predicting protein/small-molecule (ligand) binding. Unlike traditional VS, which focuses on a single protein conformation, ensemble VS better accounts for protein flexibility by predicting binding to multiple protein conformations. Each compound is thus associated with a spectrum of scores (one score per protein conformation) rather than a single score. To effectively rank and prioritize the molecules for further evaluation (including experimental testing), researchers must select which protein conformations to consider and how best to map each compound’s spectrum of scores to a single value, decisions that are system-specific. EnOpt uses machine learning to address these challenges. We perform benchmark VS to show that for many systems, EnOpt ranking distinguishes active compounds from inactive or decoy molecules more effectively than traditional ensemble VS methods. To encourage broad adoption, we release EnOpt free of charge under the terms of the MIT license.

Keywords

Ensemble virtual screening
Computer-aided drug design
Molecular docking
User-friendly software
Machine learning
Decision trees
Subject terms

Drug discovery
Computational biology and bioinformatics
Machine learning
Software
Virtual drug screening
National Institute of General Medical Sciences1R01GM132353 issue-copyright-statement© Springer Nature Limited 2024
==== Body
pmcIntroduction

Many biological processes, from signaling to cellular development, rely on interactions between proteins and small-molecule ligands. The inherent flexibility of many proteins is key to understanding these interactions, the processes they enable, and strategies for designing molecules that bind to them and alter their activity. If a protein is naturally rigid, a single (e.g., crystallographic) structure may be sufficient for effective computer-aided ligand design. However, many flexible proteins sample multiple biologically relevant conformations in vivo. Such conformational changes are often critical to biological function, and different conformations can accommodate very different small-molecule ligands.

To account for these conformational changes, ensemble virtual screening (VS) uses molecular docking to predict ligand binding to multiple distinct but representative binding-pocket conformations, collectively known as a conformational ensemble1–3. Such ensembles can include structures identified through X-ray crystallography, NMR, molecular dynamics simulations, and deep-learning structure-prediction methods such as AlphaFold1,4. Docking into an ensemble generates a spectrum of docking scores for each compound, one score per protein conformation. To prioritize compounds for subsequent experimental testing, researchers use these scores to rank compounds from the best-predicted binder to the worst.

To use ensemble VS effectively, researchers must address two challenges. First, they must determine which pocket conformations to include in the ensemble. Some previous work suggests that accuracy strictly increases with the number of conformations included5. In contrast, in some instances, identifying an optimal subensemble may be best6. Optimal subensembles should arguably include conformations that plausibly represent potential holo states compatible with favorable docked poses. For example, one might reasonably exclude pocket conformations that are entirely collapsed or so open that the protein is nearly unfolded. Similarly, one might reasonably favor conformations that are similar to those known to bind small-molecule ligands (e.g., per crystallography)7.

Second, researchers must determine how to map each compound’s spectrum of docking scores to a single composite score suitable for ranking. Common approaches include calculating the average score across all conformations (ensemble average) or identifying the single best score among all conformations (ensemble best). These ranking approaches are unsupervised and so do not require users to consider experimental activity labels for training8. However, they sometimes fail to accurately capture different conformation’s contributions to protein-ligand binding. For example, the ensemble-average score assumes that all conformations contribute equally to binding; in reality, compounds often bind preferentially to a limited subset of conformational states1. On the other hand, the ensemble-best score disregards most ensemble members and assumes that only a single conformation is compatible with binding. Few ligands bind their protein targets via such a rigid “lock and key” mechanism1. The ensemble-average and ensemble-best schemes thus represent two extremes, but the reality of binding lies somewhere in between. Optimal methods must identify relevant conformational subensembles and weigh their members appropriately.

Various machine-learning (ML) approaches aim to intelligently use ensemble docking scores to rank candidate ligands, with promising accuracy in many cases5,9,10. However, these methods have several limitations. For one, ML-based ensemble docking methods are often validated using benchmarking datasets that contain far more protein conformations than are typically available for practical drug discovery. Methods that are more aligned with practical research applications are often tailored only for specific targets or protein families. Additionally, many users otherwise interested in using ensemble docking cannot implement ML methods due to time and resource limitations, and may not be familiar with best practices for interpreting the resulting data.

To address these challenges, we introduce Ensemble Optimizer (EnOpt), an accessible tool that streamlines ensemble-docking analyses. EnOpt provides an easy-to-use, target-agnostic platform for training dataset-specific models, based on ML methods known to perform well in ensemble virtual screening. Given the docking scores associated with a completed ensemble VS, EnOpt uses an ML method called gradient boosted trees to identify appropriate subensembles and ensemble composite scores. We use multiple target proteins docked with compound sets including known ligands (positive controls) to assess how well EnOpt (1) ranks active above inactive/decoy compounds and (2) identifies which conformations most contribute to accurate binding predictions. We also show that EnOpt compares favorably to more traditional ensemble VS approaches.

Methods

Methodological details

Input format and options

As input, EnOpt accepts a CSV-formatted n×m matrix of docking scores, organized as a set of n compounds (rows) docked into m target-protein conformations (columns). All compounds should be labeled uniquely, and a separate CSV file should list the known ligands that will serve as positive controls. This format allows users to easily import compounds and scores from different sources with different naming standards. The git repository (https://github.com/durrantlab/EnOpt) provides example input CSV files to illustrate the necessary format.

Model selection and hyperparameters

EnOpt trains a machine-learning model to predict the probability that a given compound is active. By default, it uses gradient-boosted trees, as implemented in XGBoost11. Users can also select random forest, as implemented in scikit-learn12. We incorporated these two methods because others have applied them successfully in the context of ensemble docking5. EnOpt’s default hyperparameters for model fitting mostly align with these two packages’ default values. Specifically, the learning rate is set to 0.3, and the maximum tree depth is set to 6. Both alpha (the L1 regularization term) and gamma (the minimum loss reduction required to split a leaf) are set to 0. By default, all features are used in the generation of each new tree node (per the default values of the colsample_by_node and max_features parameters in XGBoost and scikit-learn, respectively).

EnOpt’s default number of estimators is set to 15, the lowest value that maintained satisfactory performance in our examples. We chose this low value to avoid overfitting and ensure generalizability. We acknowledge that random forest models are generally robust against overfitting because they rely on large ensembles of trees13, which mitigates the impact of individual tree variability by aggregating predictions. Nevertheless, we thought it best to exercise caution when making predictions within a dataset that has so few positive controls. We expect that models trained with a larger number of estimators may perform better in some cases, as shown by Ricci et al.5. However, given that compound libraries often include few positive controls (known ligands), we worry that using a much larger number of estimators would risk overfitting, biasing the model towards the included known ligands.

Users may specify different hyperparameters via a JSON input file. If users are uncertain which specific hyperparameters are best suited to their system, they can also use EnOpt’s built-in hyperparameter-optimization feature. See the documentation for more information, and Supplementary Table S1 for the optimized XGBoost hyperparameters used in the present study.

Cross-validation training

EnOpt generates 3-fold cross-validation training and testing sets from the user’s input data. Specifically, the user’s data is divided into three non-overlapping parts with equivalent known/unknown ratios. Each of these three parts (folds) is used once as a testing set, while the remaining two are combined to form the corresponding training set. This cross-validation approach allows the user to verify that EnOpt performance does not depend on the specific examples in a given training set (i.e., the three EnOpt models trained on each fold’s training set should have similar accuracies when evaluated on the respective testing sets).

During training, EnOpt assumes only known ligands are positive examples (actives) and that all other compounds do not bind. This assumption is not necessarily valid; indeed, the goal is often to identify novel ligands among the other compounds. But it is valid for the great majority of uncharacterized compounds and is a common assumption made when training novel scoring functions or generating benchmarking datasets5,14.

Prediction strategy

After training, EnOpt assesses the probability that each compound is active using a leave-one-out approach within the cross-validation framework. This approach seeks to address two challenges. First, given that many VS include few known actives (positive controls), withholding data from training (for the final evaluation) is often impractical. At the same time, EnOpt must also prevent data leakage, i.e., it must ensure that each compound is evaluated using a model that has not been trained on that same compound.

To address these challenges, EnOpt separately considers each compound in the VS dataset. For each, it selects the one model of the three generated during cross-validation whose training set did not include the compound being evaluated. EnOpt then uses that model to assess the probability that the compound is active (i.e., to assign the final EnOpt score). This approach effectively balances the need to use all available data with the requirement for valid, unbiased predictions. Conceptually, EnOpt can be seen as training a regressor on the input data itself, rather than a generalized predictor applicable to outside (withheld) data.Figure 1 EnOpt’s interactive visual summary. (A) Top N compounds by predicted activity probability (EnOpt score), including positive controls. (B) Docking-score distributions of these compounds. (C) How often compound docking into the specified conformation gave the best docking score of the ensemble (y-axis: number of compounds; x-axis: conformation IDs; green bars: top M most predictive conformations, per the tree models). N and M are user-defined parameters.

Readable and pipeline-friendly output

Activity probabilities and performance metrics

As output, EnOpt reports the predicted probability that each compound is active per its trained machine-learning models (i.e., the EnOpt score). It ranks the compounds by these predicted probabilities and uses that ranking to calculate performance metrics such as AUROC, PRAUC, BEDROC, and enrichment factors (see the Evaluating predictive performance section below for details). For comparison, EnOpt also assesses the predictive performance of each individual conformation in the ensemble by ranking the compounds by their associated single-conformation docking scores and calculating the same performance metrics. Both the EnOpt machine-learning predictions and the single-conformation analyses are saved to CSV files, enabling easy integration into existing computational pipelines and allowing for straightforward comparison of different scoring methods.

Feature importance metrics

Aside from calculating activity probabilities and performance metrics, EnOpt also uses feature importance methods to assign numerical scores to each conformation based on the contribution of each to the predictions. For gradient-boosted tree models, the importance metric is information gain15. For random forest models, the importance metric is Gini impurity16. The CSV output reports the per-conformation feature importance values associated with each of the three cross-validation models, allowing users to assess the consistency of importance rankings. Identifying the most important conformations is useful because these conformations can then be prioritized in subsequent analyses (e.g., re-docking into a smaller ensemble with more computationally intensive programs or settings, examining the protein-ligand interactions of the most predictive conformations in more detail, etc.).

Interactive summary

Finally, EnOpt provides an interactive summary (Fig. 1) using Plotly17. This summary focuses on the top-ranked docked compounds so the user can easily assess docking results and select compounds for further study. Specifically, EnOpt provides information about all known active compounds as well as the top decoy/inactive molecules (50 by default). It also provides the AUROC metric obtained when the compounds are ranked by their EnOpt-predicted probabilities of being active, allowing the user to easily evaluate the model. Further, the interactive summary includes the average feature importance value associated with each conformation so users can quickly identify which protein conformations contribute most to the models’ predictivity.

EnOpt testing

Target and compound selection

AChE benchmark

To fairly test EnOpt, we sought to replicate the challenges typical of a real-world VS. We generated an initial docking dataset using a protocol designed to emulate a target-specific docking inquiry focused on acetylcholinesterase (AChE). We chose AChE because it is a drug target implicated in several human diseases, has multiple high-quality crystal structures available in the Protein Data Bank18 for ensemble VS, and has a large amount of bioactivity data available in PubChem19. We compiled a set of small molecules that included known ligands from PubChem as well as chemically similar decoy compounds identified using DUD-E’s “Generate DUD-E Decoys” feature14. The ratio of decoys to known ligands was approximately 50:1.

DUD-E benchmark

As a second assessment, we generated a larger test set of seven crystallographic protein-receptor ensembles taken from the DUD-E database. We first considered the 102 receptors cataloged in the DUD-E database. We discarded those receptors with two or fewer wild-type structures deposited in the PDB because they allowed only for trivially small ensemble sizes. We similarly discarded those receptors with 100 or more wild-type PDB structures because most real-world docking projects will not have access to so large an ensemble. We further focused on human proteins implicated in human diseases, per UniProt20, with five or more high-resolution PDB structures (< 2.5 Å). To ensure the test set covered diverse protein receptors, we eliminated redundant proteins belonging to the same SCOP2 protein family21,22. Another six proteins were removed from the list because, given the sizes of their ensembles and associated DUD-E compound sets, the ensemble VS would have required >200,000 dockings. Further information on selected target PDB entries can be found in the supporting information (see Supplementary Table S2). Lists of all docked compounds are available on zenodo (see Availability of data and materials section, below).

LIT-PCBA benchmark

As a third assessment, we created an additional test set of 15 receptor/ligand systems derived from the LIT-PCBA dataset (Supplementary Table S2). This dataset is specifically designed to provide a realistic and unbiased evaluation of VS methods. Unlike DUD-E, which has been shown to contain hidden biases that can lead to over-optimistic performance estimates, LIT-PCBA is based on carefully curated experimental data from PubChem bioassays23.

We first downloaded the files for the 15 protein/ligand systems from the LIT-PCBA server. To ensure EnOpt had suitable ensembles for analysis, we expanded each system by downloading all structures from the Protein Data Bank18 that shared the same UniProt ID, thus representing the same protein. We then discarded those structures with mutations, resolutions greater than 3 Å (including all NMR structures, which had no specified resolutions), etc. The remaining structures were categorized as either apo or holo, depending on the presence of a non-polymer ligand.

To assemble a manageable benchmark set from these files, we next identified at most 20 active compounds for each system from among the actives that the LIT-PCBA database provides. Two of the LIT-PCBA systems have fewer than 20 actives (ADRB2 and ESR1_ago); for these, we retained all actives. For the remaining 13 systems, we selected 20 representative actives by clustering the LIT-PCBA active compounds using the Butina algorithm24 and selecting one compound from each of 20 distinct clusters. After selecting active compounds for each system, we randomly selected 40 LIT-PCBA inactive compounds for each active.

To create ensembles of protein structures for each of the 15 LIT-PCBA systems, we considered the identified holo (ligand-bound) and apo (unbound) protein structures separately. If there were more than 10 structures of either type, we clustered them based on structural similarity (RMSD) and selected representative conformations; otherwise, we used all available structures.

Protein and ligand preparation

For the AChE and DUD-E protein targets, we removed (1) duplicate domains or multimers, (2) all co-crystallized small molecules and peptides, and (3) all water molecules. We aligned (superimposed) all structures of each ensemble using PyMol25 or UCSF Chimera26. We used OpenBabel27 to prepare target-protein and small-molecule structures for docking, protonating at pH 7.4 (e.g., “obabel input.pdb -O output.pdbqt -p 7.4 -xr” for proteins, and “obabel input.mol2 -O output.pdbqt -p 7.4” for ligands).

For the LIT-PCBA protein targets, we removed all co-crystallized small molecules and peptides. We also removed crystallographic water molecules but retained metal cations, which are often important for small-molecule binding. We aligned the structures using Biopython28 and used PDB2PQR29 to prepare the target-protein structures for docking, protonating at pH 7. Similarly, we used OpenBabel27 to prepare small-molecule structures for docking, also protonating at pH 7.

Docking

For the initial AChE and subsequent DUD-E screens, we used smina30 for docking. Each smina run predicts several docked poses; we considered only the top-scoring pose for further analysis. We specified the docking box using smina’s autoselection feature with a padding of 4 Å on all six sides, centered on a crystallographic ligand taken from an appropriate holo complex. For the initial docking dataset using AChE, the ligand structure and position were taken from a PDB entry for acetylcholinesterase bound to Donepezil (ID 4EY7)18,31. For the additional test cases, the ligand structure was provided by the DUD-E database14. To improve run times, we set smina’s exhaustiveness parameter to 4 for all docking runs, reduced from the default (8).

We used a somewhat different docking protocol for the LIT-PCBA screens, with the goal of demonstrating EnOpt’s broad effectiveness in the context of various pipelines. Specifically, we used AutoDock Vina 1.1.2 (Vina)32 for the LIT-PCBA docking. Each Vina run also predicts several docked poses; we again considered only the top-scoring pose. We determined docking box centers by calculating the geometric center of selected crystallographic ligands. We then used MolModa33 to determine the docking box sizes required to entirely encompass these reference ligands. Finally, we used Vina to dock the prepared ligands (both actives and inactives) into each of the corresponding protein structures. For the LIT-PCBA runs, we used Vina’s default exhaustiveness setting (8).

Evaluating predictive performance

For each benchmark VS, we used in-house scripts to reformat the docking outputs into an EnOpt-compatible n×m matrix (CSV). We then trained EnOpt models to predict the probability that a given compound is active and ranked all compounds by those probabilities (i.e., EnOpt scores), from most to least probable. To assess EnOpt’s ability to identify true ligands, we calculated multiple evaluation metrics using these ranked lists.

The primary metric was the area under the receiver operating characteristic curve (AUROC). ROC curves relate true positive rates (TPR) to false positive rates (FPR) across various classification thresholds. We considered all possible thresholds, each partitioning the ordered list into presumed actives and inactives. Because we knew which compounds were true actives and inactives (decoys), we could calculate the true and false positive rates for each threshold. Plotting each (FPR, TPR) pair defines the ROC curve. The area under this curve indicates how well the model classifies the candidate ligands. More precisely, it describes the probability that a randomly chosen true active will rank better than a randomly chosen inactive/decoy.

As a second metric, we considered the area under the precision-recall curve (PRAUC) as implemented in scikit-learn12,34. Precision quantifies the accuracy of positive predictions as the ratio of true positives to all predicted positives. Recall measures the completeness of positive predictions as the ratio of true positives to all actual positives. The PRAUC provides a single value that summarizes the trade-off between precision and recall across various classification thresholds, offering a comprehensive view of the model’s performance, particularly in cases with unbalanced datasets.

As a third metric, we considered the Boltzmann-enhanced discrimination of receiver operating curve (BEDROC), as implemented in RDKit35. This metric modifies the traditional ROC by emphasizing the early detection of active compounds. BEDROC applies an exponential weight to compounds ranked higher in the list, reflecting the practical importance of identifying active compounds early in VS campaigns. The BEDROC calculation uses a parameter α to control the emphasis on early recognition. A higher α value focuses the metric on a smaller percentage of top-ranked compounds36. We selected α=5, which effectively focuses on the top ∼ 20% of ranked results37.

As a fourth metric, we considered enrichment factors. The enrichment factor quantifies a model’s ability to preferentially rank active compounds higher than would be expected by random selection. Specifically, it is the ratio of the proportion of actives in the top N percent of ranked compounds to the proportion of actives in the entire dataset. We set N to 20% to maintain consistency with our BEDROC α parameter. An enrichment factor greater than 1 indicates that the model performs better than random selection in identifying active compounds within the top-ranked results.

To compare our use of predicted probabilities to more traditional approaches for ranking by ensemble scores, we separately re-ranked the compounds by their ensemble-average and ensemble-best docking scores. We similarly calculated AUROC, PRAUC, BEDROC, and enrichment factor values for these compound rankings. While most of our discussed results rely on the AUROC metric, we have included results for all four metrics in Supplementary Table S3.

Scientific contribution

EnOpt is a user-friendly virtual-screening (VS) tool that trains machine learning (ML) models to predict the probability that a given molecule is active. It learns to map the docking scores associated with an ensemble VS (i.e., multiple compounds, each docked into multiple protein-receptor conformations) to probability scores for compound ranking. EnOpt builds on previous research that uses ML to process ensemble-VS scores, but it wraps its methods in an easy-to-use tool that outputs predictions in an intuitive, visual format.

Results

EnOpt ranking compares favorably to traditional ranking schemes

We calculated AUROC values to compare EnOpt compound rankings to ensemble-average and ensemble-best rankings, which are commonly used to map ensemble docking scores associated with a given ligand to a single value for ranking. Ensemble-average and ensemble-best rankings tend to perform similarly within examples, but EnOpt outperforms them in most cases (Table 1 and Fig. 2). Interestingly, some protein receptors are much better suited to ensemble docking using EnOpt’s methods than others, e.g., ACHE, XIAP, OPRK1, and VDR (Table 1).

Interestingly, EnOpt’s performance improvement over the ensemble-average and ensemble-best methods is more pronounced when applied to the DUD-E dataset than the LIT-PCBA dataset.

This disparity suggests that EnOpt tree models may have detected artifactual patterns that distinguish DUD-E decoys from actives, which may not generalize to real-world scenarios. The LIT-PCBA compounds were selected through rigorous curation, making that dataset a more challenging and realistic test. Specifically, LIT-PCBA is derived from dose-response PubChem bioassays and has undergone additional processing to remove compounds that are likely false positives, associated with potential assay artifacts, outliers in terms of molecular properties, etc.23 Further, the LIT-PCBA inactive compounds are experimentally confirmed inactives rather than presumed-inactive decoys.

Regardless, EnOpt still performs well overall, with an average AUROC of 0.745 across all targets, compared to 0.588 for ensemble-average and 0.599 for ensemble-best methods. The XGBoost results are reported in Table 1 and Fig. 2, but we found no substantial difference in performance between XGBoost and Random Forest models (see Supplementary Fig. S1).

To further contextualize EnOpt’s performance, we also evaluated the predictivity of each individual conformation within the ensembles. For both the DUD-E and LIT-PCBA datasets, we calculated AUROC values after ranking the compounds by the docking scores associated with each single conformation (Supplementary Fig. S2). Notably, EnOpt’s performance was almost always at least comparable with the best single-conformation screen, and in many cases (e.g., ACHE, GLCM, HXK4, XIAP, ADRB2, GBA, OPKR2, VDR), it was dramatically better (Supplementary Fig. S2).

These results demonstrate EnOpt’s ability to effectively synthesize information from multiple conformations, often achieving better performance than could be obtained from any single conformation alone.Table 1 Descriptions of the benchmark ensemble VS. For each target protein, we show the ensemble size, active-to-decoy ratio, and range of conformational importance values across each ensemble (“Importance Range”; see Supplementary Fig. S3 for distributions). We also report the AUROC when compounds are ranked by the EnOpt score (“EnOpt AUROC”, gradient-boosted trees), as well as the AUROC values associated with ensemble-average (“EA AUROC”) and ensemble-best (“EB AUROC”) rankings for comparison. The best AUROC is highlighted in bold.

Target protein	Ensemble size	Active-to-decoy ratio	Importance values (range)	EnOpt AUROC	EA AUROC	EB AUROC	
ACHE	7	0.0340	0.272	0.859	0.312	0.290	
COMT	9	0.0105	1.000	0.511	0.585	0.566	
GLCM	12	0.0140	0.080	0.873	0.535	0.528	
HXK4	21	0.0192	0.060	0.941	0.769	0.745	
KIF11	23	0.0166	0.118	0.988	0.876	0.898	
TGFR1	17	0.0154	0.108	0.983	0.850	0.865	
TRYB1	9	0.0189	0.683	0.731	0.723	0.736	
XIAP	22	0.0190	0.094	0.959	0.495	0.557	
ADRB2	10	0.0245	0.526	0.616	0.332	0.330	
ALDH1	12	0.0245	0.988	0.600	0.617	0.634	
ESR1_ago	10	0.0247	0.564	0.807	0.624	0.683	
ESR1_ant	10	0.0248	0.826	0.590	0.625	0.682	
FEN1	5	0.0244	0.481	0.534	0.581	0.567	
GBA	17	0.0244	0.132	0.875	0.650	0.641	
IDH1	11	0.0245	0.413	0.671	0.682	0.670	
KAT2A	8	0.0245	0.648	0.622	0.508	0.506	
MAPK1	12	0.0246	0.753	0.621	0.474	0.491	
MTORC1	12	0.0244	0.928	0.678	0.408	0.403	
OPRK1	8	0.0244	0.146	0.864	0.402	0.574	
PKM2	12	0.0244	0.995	0.598	0.690	0.636	
PPARG1	19	0.0248	0.649	0.717	0.735	0.719	
TP53	16	0.0248	0.645	0.651	0.644	0.644	
VDR	10	0.0244	0.139	0.850	0.405	0.421	

Figure 2 AUROC values for AChE/DUD-E (top) and LIT-PCBA (bottom) benchmark ensemble VS. EnOpt’s gradient-boosted-trees models (yellow) generally outperform ensemble average/best (light/dark purple). Red and blue lines indicate AUROC values of 0.5 and 0.75, respectively. See Supplemental Fig. S1 for a similar comparison of EnOpt’s random forest models.

EnOpt suggests that binding sometimes depends on only a few conformations

Our EnOpt tree models suggest that certain standout conformations are often far more important for identifying true ligands than other conformations.“Importance” in this context refers to feature importance as computed by the trained models. This importance metric provides helpful pharmacological insights because conformations with high importance are likely the most useful for discerning between actives and other compounds (i.e., the active-compound docking scores associated with these conformations tend to differ meaningfully from the non-active/decoy scores).

Most distributions of conformation importance values in our benchmark VS had asymmetrical distributions (Supplementary Fig. S3), suggesting a small subset of conformations (sometimes comprised of a single conformation) often contributed disproportionately to EnOpt’s ability to distinguish between active and inactive compounds. While this pattern appears to be receptor-dependent, it indicates that identifying a smaller subensemble of conformations could be an effective strategy for interpreting VS results in many cases.

Impact of ensemble size and active-to-decoy ratio on EnOpt performance

Although increasing the number of features and data used for training often improves machine-learning accuracy, many ensemble VS must draw on only limited numbers of available receptor conformations, positive controls (known ligands), and/or negative controls (compounds known to bind poorly). To mimic these challenges, we performed benchmark VS using ensembles of fewer than 25 protein conformations. Additionally, our benchmark datasets include differing active-to-decoy ratios (Table 1).

We first investigated whether EnOpt’s performance correlates with ensemble size. The 23 benchmark VS have ensemble sizes ranging from 5 to 23 conformations (Table 1). We performed a linear regression analysis to characterize the relationship between ensemble size and EnOpt AUROC (Table 1). For the seven-member DUD-E set (COMT, GLCM, HXK4, KIF11, TGFR1, TRYB1, and XIAP), we found only a weak but statistically significant correlation (R2=0.675, p = 0.023).

For the more challenging LIT-PCBA dataset, which better mimics real-world VS scenarios, we found no such correlation (R2=0.063, p = 0.368).

These results demonstrate the robustness of EnOpt’s predictions across varying ensemble sizes. The lack of strong correlation, particularly in regards to the more realistic LIT-PCBA dataset, suggests that EnOpt can effectively handle ensembles of different sizes without compromising performance.

We also examined the impact of active-to-decoy ratio on EnOpt’s performance. Interestingly, we found that performance is relatively consistent across varying ratios. The benchmark VS have ratios ranging from 1.05–3.40% actives (Table 1). A linear regression analysis between this ratio and the DUD-E AUROC values yielded an R2 of 0.375, with a p-value of 0.144 (i.e., not significant). A similar analysis applied to the LIT-PCBA dataset yielded an R2 = 0.027, with a p-value of 0.558 (not significant). It is worth noting that the LIT-PCBA screens had similar active-to-decoy ratios across all systems. Nevertheless, these results suggest that EnOpt’s performance is generally robust to variations in the proportion of active compounds as well.

Discussion

Model interpretability and ease of use

EnOpt’s core unique features are model simplicity, generality, and ease of use. We show that simple and interpretable tree methods accurately predict actives without overfitting. These methods are transparent, allowing users to easily understand which data features are integral to the model’s decisions. The importance values the models assign to each conformation are key to their interpretability, as they can suggest smaller but effective sub-ensembles for follow-up ensemble screens using more computationally intensive scoring functions.

EnOpt’s tree models are nonlinear, consistent with how ligand binding occurs in vivo. When a ligand interacts with a dynamic protein, multiple pocket conformations contribute to binding, to varying degrees. Neither a single “best” score nor an unweighted linear combination (mean) properly represents such a binding event. By identifying the importance of different pocket conformations, tree models can more accurately represent binding to a dynamic protein.

EnOpt output accommodates a diverse user base

EnOpt is also designed to accommodate a diverse user base. Expert computationalists can parse its text-based CSV output using command-line tools that can be easily incorporated into automated pipelines. Students and scientists less familiar with the command line will benefit from EnOpt’s visual, interactive plot summary (Fig. 1).

The summary provides information about all known actives as well as the top decoy compounds (which may, in fact, be uncharacterized actives), ranked by their EnOpt-predicted probabilities of being active (i.e., EnOpt scores). It also includes the range of docking scores, predicted probabilities for each compound, and the ensemble-average and ensemble-best values.

Hovering over any compound label reveals detailed information as a pop-up, allowing users to visually compare likely actives among the decoys and identify conformations with unusually poor or favorable scores. Users can then remove poorly scoring outlier conformations when generating subensembles or investigate favorably scoring outliers further as particularly suitable conformations for binding.

The summary also provides the AUROC values of the three models constructed during 3-fold cross-validation. Comparing these values allows users to verify that none of the models is overfitting to its subset of the docking matrix. Because EnOpt randomly splits the training data into three folds, the distribution of compounds across these folds can sometimes lead to imbalances that affect model performance. If any of the three models performs substantially better or worse than the others, users may opt to re-run EnOpt, which will generate a new random data split. This simple action can often result in more balanced training sets and, consequently, better-trained models with more consistent performance across folds.

The output also allows users to assess whether EnOpt’s conformation importance rankings are reliable. Both random forest and gradient-boosted trees incorporate randomness into their training, so EnOpt is nondeterministic. This randomness has only a small impact on the predicted activity probabilities. However, it may more substantially impact EnOpt’s conformational rankings if many (or all) ensemble conformations have roughly equal importance. In such cases, small random differences in the assigned importance values can reorder the conformations.

Recommendations for use

Ensemble generation

We evaluated EnOpt using ensembles of experimentally derived structures. However, in cases where limited or no experimental structures are available, researchers may consider predicted structures. Molecular dynamics (MD) simulations are rich sources of diverse computationally derived ensembles that capture a target protein’s atomic motions. Such simulations can effectively complement conformations available through experiments1. Methods for predicting protein structures, such as homology-modeling and machine-learning approaches (e.g., AlphaFold24), can also provide conformations for ensemble docking, though the suitability of the latter for docking remains a subject of debate38–40. Regardless, we recommend including at least one experimental structure to ensure that the conformational ensemble encompasses physically real pocket conformations.

Compound selection

Although there was little correlation between the active-to-decoy ratio and EnOpt predictivity in the range of ratios we tested (1.05–3.40%, i.e., ∼ 98% decoy or unknown compounds; see Table 1), we generally recommend including as many positive controls (known ligands) in the VS as possible. EnOpt uses 3-fold cross-validation to evaluate for overfitting, assigning an equal proportion of the provided true ligands to each fold. Random splits ensure a balanced distribution of knowns across the three sets, but if each model is trained on only a few positive controls, the AUROC calculations used to evaluate model performance may not be meaningful.

Conclusion

This study introduces Ensemble Optimizer (EnOpt), a tool to rapidly train protein-specific machine-learning models that can be easily applied to ensemble virtual screens. We show that EnOpt model performance is not heavily dependent on the number of features (conformations in the ensemble) or the proportion of positive controls in the dataset, illustrating its applicability across diverse docking datasets.

EnOpt’s modular codebase can be easily expanded with additional model types beyond gradient-boosted trees (the default). That said, EnOpt’s tree models already provide a flexible platform for users to analyze their ensemble-VS results. EnOpt is available under the MIT license. Users can download a copy free of charge from http://git.durrantlab.com/jdurrant/EnOpt.

Supplementary Information

Supplementary Information.

Abbreviations

EnOpt Ensemble optimizer

VS Virtual screening

AUROC Area under the receiver operating characteristic

PRAUC Area under the precision-recall curve

BEDROC Boltzmann-enhanced discrimination of receiver operating characteristic

EF Enrichment factor

PDB Protein Data Bank

AChE Acetylcholinesterase

COMT Catechol O-methyltransferase

GLCM Acid-beta-glucosidase

HXK4 Hexokinase 4

KIF11 Kinesin Eg5

TGFR1 TGF-beta receptor type 1

TRYB1 Beta-tryptase

XIAP X-linked inhibitor of apoptosis protein

FPR False positive rate

TPR True positive rate

MD Molecular dynamics

Supplementary Information

The online version contains supplementary material available at 10.1038/s41598-024-71699-3.

Acknowledgements

We would like to thank Yogindra Raghav for his contributions in generating the initial proof-of-concept code. We also thank Darian Yang for assistance in collating and pruning ideas.

Author contributions

R.B. and J.D.D. designed the tool (EnOpt) and contributed to the codebase. R.B., A.W., and J.D.D. ran docking calculations and evaluated EnOpt analyses. R.B. and J.D.D. composed the paper. All authors read and approved the manuscript.

Funding

This work was supported by the National Institute of Health (1R01GM132353) and the University of Pittsburgh’s Center for Research Computing, RRID:SCR_022735 (supported by NSF OAC-2117681). The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health or the National Science Foundation. The funders had no role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Data availability

The EnOpt codebase, including testing scripts and demo outputs, is available at http://git.durrantlab.com/durrantlab/EnOpt. All PDB structures described in the text can be downloaded from the RCSB PDB. The active and decoy/inactive molecules used in all screens are available via zenodo.

Competing interests

The authors declare no competing interests.

Use of assistive writing technologies

We used assistive writing tools such as Grammarly, OpenAI’s ChatGPT, and Anthropic’s Claude during manuscript preparation to edit and refine our text. In all cases, the authors critically reviewed, revised, and selectively implemented any recommended edits to ensure accuracy and clarity. The authors alone are responsible for the paper’s quality and content.

Publisher's note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
==== Refs
References

1. Amaro RE Li WW Emerging methods for ensemble-based virtual screening Curr. Top. Med. Chem. 2010 10 3 13 10.2174/156802610790232279 19929833
Amaro, R. E. & Li, W. W. Emerging methods for ensemble-based virtual screening. Curr. Top. Med. Chem. 10, 3–13 (2010).19929833 10.2174/156802610790232279
2. Lionta E Spyrou G Vassilatis KD Cournia Z Structure-based virtual screening for drug discovery: Principles, applications and recent advances Curr. Top. Med. Chem. 2014 14 1923 1938 10.2174/1568026614666140929124445 25262799
Lionta, E., Spyrou, G., Vassilatis, K. D. & Cournia, Z. Structure-based virtual screening for drug discovery: Principles, applications and recent advances. Curr. Top. Med. Chem. 14, 1923–1938 (2014).25262799 10.2174/1568026614666140929124445
3. Amaro RE Ensemble docking in drug discovery Biophys. J. 2018 114 2271 2278 10.1016/j.bpj.2018.02.038 29606412
Amaro, R. E. et al. Ensemble docking in drug discovery. Biophys. J. 114, 2271–2278 (2018).29606412 10.1016/j.bpj.2018.02.038
4. Sala D Engelberger F Mchaourab H Meiler J Modeling conformational states of proteins with alphafold Curr. Opin. Struct. Biol. 2023 81 102645 10.1016/j.sbi.2023.102645 37392556
Sala, D., Engelberger, F., Mchaourab, H. & Meiler, J. Modeling conformational states of proteins with alphafold. Curr. Opin. Struct. Biol. 81, 102645. 10.1016/j.sbi.2023.102645 (2023).37392556 10.1016/j.sbi.2023.102645
5. Ricci-Lopez J Aguila SA Gilson MK Brizuela CA Improving structure-based virtual screening with ensemble docking and machine learning J. Chem. Inf. Model. 2021 61 5362 5376 10.1021/acs.jcim.1c00511 34652141
Ricci-Lopez, J., Aguila, S. A., Gilson, M. K. & Brizuela, C. A. Improving structure-based virtual screening with ensemble docking and machine learning. J. Chem. Inf. Model. 61, 5362–5376 (2021).34652141 10.1021/acs.jcim.1c00511
6. Rao S Improving database enrichment through ensemble docking J. Comput. Aided Mol. Des. 2008 22 621 627 10.1007/s10822-008-9182-y 18253700
Rao, S. et al. Improving database enrichment through ensemble docking. J. Comput. Aided Mol. Des. 22, 621–627 (2008).18253700 10.1007/s10822-008-9182-y
7. Kumar A Zhang KY A cross docking pipeline for improving pose prediction and virtual screening performance J. Comput. Aided Mol. Des. 2018 32 163 173 10.1007/s10822-017-0048-z 28836076
Kumar, A. & Zhang, K. Y. A cross docking pipeline for improving pose prediction and virtual screening performance. J. Comput. Aided Mol. Des. 32, 163–173 (2018).28836076 10.1007/s10822-017-0048-z
8. Willett P Combination of similarity rankings using data fusion J. Chem. Inf. Model. 2013 53 1 10 10.1021/ci300547g 23297768
Willett, P. Combination of similarity rankings using data fusion. J. Chem. Inf. Model. 53, 1–10 (2013).23297768 10.1021/ci300547g
9. Morris CJ Stern JA Stark B Christopherson M Della Corte D Milcdock: Machine learning enhanced consensus docking for virtual screening in drug discovery J. Chem. Inf. Model. 2022 1 1 10.1021/acs.jcim.2c00705
Morris, C. J., Stern, J. A., Stark, B., Christopherson, M. & Della Corte, D. Milcdock: Machine learning enhanced consensus docking for virtual screening in drug discovery. J. Chem. Inf. Model. 1, 1. 10.1021/acs.jcim.2c00705 (2022).10.1021/acs.jcim.2c00705
10. Tian S Assessing an ensemble docking-based virtual screening strategy for kinase targets by considering protein flexibility J. Chem. Inf. Model. 2014 54 2664 2679 10.1021/ci500414b 25233367
Tian, S. et al. Assessing an ensemble docking-based virtual screening strategy for kinase targets by considering protein flexibility. J. Chem. Inf. Model. 54, 2664–2679 (2014).25233367 10.1021/ci500414b
11. Chen, T. & Guestrin, C. XGBoost: A scalable tree boosting system. In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, 785–794 (2016).
12. Pedregosa F Scikit-learn: Machine learning in python J. Mach. Learn. Res. 2011 12 2825 2830
Pedregosa, F. et al. Scikit-learn: Machine learning in python. J. Mach. Learn. Res. 12, 2825–2830 (2011).
13. Breiman L Random forests Mach. Learn. 2001 45 5 32 10.1023/A:1010933404324
Breiman, L. Random forests. Mach. Learn. 45, 5–32 (2001).10.1023/A:1010933404324
14. Huang N Shoichet BK Irwin JJ Benchmarking sets for molecular docking J. Med. Chem. 2006 49 6789 6801 10.1021/jm0608356 17154509
Huang, N., Shoichet, B. K. & Irwin, J. J. Benchmarking sets for molecular docking. J. Med. Chem. 49, 6789–6801. 10.1021/jm0608356 (2006).17154509 10.1021/jm0608356
15. XGBoost developers. XGBoost python API documentation (2022).
16. Scikit-learn developers. Tree model mathematical formulation: Classification criteria (2023).
17. plotly technologies Inc. Collaborative data science (2015).
18. Berman HM The protein data bank Nucleic Acids Res. 2000 28 235 242 10.1093/nar/28.1.235 10592235
Berman, H. M. et al. The protein data bank. Nucleic Acids Res. 28, 235–242 (2000).10592235 10.1093/nar/28.1.235
19. Kim S PubChem in 2021: New data content and improved web interfaces Nucleic Acids Res. 2021 49 D1388 D1395 10.1093/nar/gkaa971 33151290
Kim, S. et al. PubChem in 2021: New data content and improved web interfaces. Nucleic Acids Res. 49, D1388–D1395 (2021).33151290 10.1093/nar/gkaa971
20. Consortium T. U. Uniprot: The universal protein knowledgebase in 2021. Nucleic Acids Res. 49, D480–D489 (2021).
21. Andreeva A Howorth D Chothia C Kulesha E Murzin AG SCOP2 prototype: A new approach to protein structure mining Nucleic Acids Res. 2014 42 D310 D314 10.1093/nar/gkt1242 24293656
Andreeva, A., Howorth, D., Chothia, C., Kulesha, E. & Murzin, A. G. SCOP2 prototype: A new approach to protein structure mining. Nucleic Acids Res. 42, D310–D314 (2014).24293656 10.1093/nar/gkt1242
22. Andreeva A Kulesha E Gough J Murzin AG The SCOP database in 2020: Expanded classification of representative family and superfamily domains of known protein structures Nucleic Acids Res. 2020 48 D376 D382 10.1093/nar/gkz1064 31724711
Andreeva, A., Kulesha, E., Gough, J. & Murzin, A. G. The SCOP database in 2020: Expanded classification of representative family and superfamily domains of known protein structures. Nucleic Acids Res. 48, D376–D382 (2020).31724711 10.1093/nar/gkz1064
23. Tran-Nguyen V-K Jacquemard C Rognan D LIT-PCBA: An unbiased data set for machine learning and virtual screening J. Chem. Inf. Model. 2020 60 4263 4273 10.1021/acs.jcim.0c00155 32282202
Tran-Nguyen, V.-K., Jacquemard, C. & Rognan, D. LIT-PCBA: An unbiased data set for machine learning and virtual screening. J. Chem. Inf. Model. 60, 4263–4273 (2020).32282202 10.1021/acs.jcim.0c00155
24. Butina D Unsupervised data base clustering based on daylight’s fingerprint and Tanimoto similarity: A fast and automated way to cluster small and large data sets J. Chem. Inf. Comput. Sci. 1999 39 747 750 10.1021/ci9803381
Butina, D. Unsupervised data base clustering based on daylight’s fingerprint and Tanimoto similarity: A fast and automated way to cluster small and large data sets. J. Chem. Inf. Comput. Sci. 39, 747–750 (1999).10.1021/ci9803381
25. Schrodinger, L. The PyMOL molecular graphics system, version 1.8 (2015).
26. Pettersen EF UCSF ChimeraX: Structure visualization for researchers, educators, and developers Protein Sci. 2021 30 70 82 10.1002/pro.3943 32881101
Pettersen, E. F. et al. UCSF ChimeraX: Structure visualization for researchers, educators, and developers. Protein Sci. 30, 70–82 (2021).32881101 10.1002/pro.3943
27. O’Boyle NM Open babel: An open chemical toolbox J. Cheminform. 2011 3 1 14 21214931
O’Boyle, N. M. et al. Open babel: An open chemical toolbox. J. Cheminform. 3, 1–14 (2011).21214931
28. Cock PJ Biopython: Freely available python tools for computational molecular biology and bioinformatics Bioinformatics 2009 25 1422 10.1093/bioinformatics/btp163 19304878
Cock, P. J. et al. Biopython: Freely available python tools for computational molecular biology and bioinformatics. Bioinformatics 25, 1422 (2009).19304878 10.1093/bioinformatics/btp163
29. Dolinsky TJ PDB2PQR: Expanding and upgrading automated preparation of biomolecular structures for molecular simulations Nucleic Acids Res. 2007 35 W522 W525 10.1093/nar/gkm276 17488841
Dolinsky, T. J. et al. PDB2PQR: Expanding and upgrading automated preparation of biomolecular structures for molecular simulations. Nucleic Acids Res. 35, W522–W525 (2007).17488841 10.1093/nar/gkm276
30. Koes DR Baumgartner MP Camacho CJ Lessons learned in empirical scoring with smina from the CSAR 2011 benchmarking exercise J. Chem. Inf. Model. 2013 53 1893 1904 10.1021/ci300604z 23379370
Koes, D. R., Baumgartner, M. P. & Camacho, C. J. Lessons learned in empirical scoring with smina from the CSAR 2011 benchmarking exercise. J. Chem. Inf. Model. 53, 1893–1904 (2013).23379370 10.1021/ci300604z
31. Cheung, J. et al. Crystal structure of recombinant human acetylcholinesterase in complex with donepezil (2012).
32. Trott O Olson AJ AutoDock vina: Improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading J. Comput. Chem. 2010 31 455 461 10.1002/jcc.21334 19499576
Trott, O. & Olson, A. J. AutoDock vina: Improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading. J. Comput. Chem. 31, 455–461 (2010).19499576 10.1002/jcc.21334
33. Kochnev Y Ahmed M Maldonado AM Durrant JD MolModa: Accessible and secure molecular docking in a web browser Nucleic Acids Res. 2024 52 gkae406 10.1093/nar/gkae406
Kochnev, Y., Ahmed, M., Maldonado, A. M. & Durrant, J. D. MolModa: Accessible and secure molecular docking in a web browser. Nucleic Acids Res. 52, gkae406 (2024).10.1093/nar/gkae406
34. scikit-learn developers. Metrics and scoring: quantifying the quality of predictions (2023).
35. developers, R. Rdkit: Open-source cheminformatics.
36. Truchon J-F Bayly CI Evaluating virtual screening methods: Good and bad metrics for the “early recognition” problem J. Chem. Inf. Model. 2007 47 488 508 10.1021/ci600426e 17288412
Truchon, J.-F. & Bayly, C. I. Evaluating virtual screening methods: Good and bad metrics for the “early recognition’’ problem. J. Chem. Inf. Model. 47, 488–508 (2007).17288412 10.1021/ci600426e
37. Chaput L Martinez-Sanz J Saettel N Mouawad L Benchmark of four popular virtual screening programs: Construction of the active/decoy dataset remains a major determinant of measured performance J. Cheminform. 2016 8 1 17 10.1186/s13321-016-0112-z 26807156
Chaput, L., Martinez-Sanz, J., Saettel, N. & Mouawad, L. Benchmark of four popular virtual screening programs: Construction of the active/decoy dataset remains a major determinant of measured performance. J. Cheminform. 8, 1–17 (2016).26807156 10.1186/s13321-016-0112-z
38. Zhang Y Benchmarking refined and unrefined alphafold2 structures for hit discovery J. Chem. Inf. Model. 2023 63 1656 1667 10.1021/acs.jcim.2c01219 36897766
Zhang, Y. et al. Benchmarking refined and unrefined alphafold2 structures for hit discovery. J. Chem. Inf. Model. 63, 1656–1667 (2023).36897766 10.1021/acs.jcim.2c01219
39. Holcomb M Chang Y-T Goodsell DS Forli S Evaluation of alphafold2 structures as docking targets Protein Sci. 2023 32 e4530 10.1002/pro.4530 36479776
Holcomb, M., Chang, Y.-T., Goodsell, D. S. & Forli, S. Evaluation of alphafold2 structures as docking targets. Protein Sci. 32, e4530 (2023).36479776 10.1002/pro.4530
40. Scardino, V., Di Filippo, J. I. & Cavasotto, C. N. How good are alphafold models for docking-based virtual screening? iScience 26 (2023).
