High-Throughput, High-Quality: Unlocking GNINA for Precision Virtual Screening

1. Introduction

Drug discovery is a complex and expensive process aimed at identifying and developing new drugs, which has seen an increasing integration of computational methodologies in recent decades. These methodologies have become crucial components in many phases of drug discovery programs, from identifying “hits” to optimizing “leads” [1]. Among these, molecular docking has emerged as a fundamental and widely used computational tool [2,3,4].

Currently, molecular docking is an established and popular in-silico tool for identifying new compounds with therapeutic potential. It is primarily used in the early stages of drug discovery due to its cost-effectiveness and potential to identify novel chemotypes, as well as provide mechanistic insights into ligand-protein interactions [2,3]. With the rapid improvement of computational platforms and the substantial increase in structural, chemical, and biological data available for an expanding number of therapeutic targets, in-silico approaches, such as docking, have undergone significant growth. Molecular docking has become a crucial component of computational procedures employed in modern drug discovery. Structure-based virtual screening techniques, which often employ docking, enable the screening of digital libraries comprising millions of compounds within affordable timescales, thereby reducing the initial costs of hit identification and increasing the likelihood of discovering promising drug candidates. Automated workflows integrating docking have been developed for screening large libraries of compounds and targets. The advancement of high-performance computing and Graphics Processing Units (GPUs) has enabled large-scale screening to become possible [2,3,5]. Although initially used as a stand-alone method, docking is now predominantly integrated into workflows that combine different computational methodologies, such as ligand-based approaches, molecular dynamics (MD), binding free energy calculations, and artificial intelligence (AI). This integration aims to overcome some of the inherent limitations of docking and better leverage information from various sources, typically resulting in improved predictive performance in terms of hit rates [2].

Despite its widespread use and considerable progress, molecular docking has inherent limitations that limit its predictive accuracy. The two most relevant and widely recognized limitations of molecular docking are the sampling of conformations and the accuracy of scoring functions. Docking involves finding the most favorable reciprocal configurations, known as “poses,” between a ligand and a protein target. However, sampling all possible conformations of the ligand and receptor is often restricted. Additionally, the target protein is frequently treated as a rigid body during docking for computational reasons, which is a significant limitation. Proteins are dynamic systems with some flexibility, which can be affected by ligand binding, known as “induced fit” effects. Approaches to address protein flexibility, such as utilizing multiple protein structures or integrating with molecular dynamics simulations, are more complex and require greater computational resources [2,3]. The scoring function is a component of docking that quantifies the interaction between the ligand and the target protein, evaluating the goodness-of-fit of a given pose or ranking different ligands. To achieve the computational efficiency required for screening many compounds, scoring functions are often simplified and based on approximate models of protein-ligand interactions. This simplified modeling is a substantial limitation, as approximate scoring functions very rarely correlate well with experimental binding affinities. Imperfections in the scoring functions remain a significant limiting factor in docking. They lead to inaccuracies in the ranking of predicted poses and poor performance in predicting binding free energy. The need to improve predictive accuracy has led to the development of several strategies to mitigate the limitations of individual scoring functions, such as consensus scoring, which combines the results of several scoring functions, and rescoring methods based on different techniques [1,2,3,4].

In the context of rational drug development, molecular docking represents one of the benchmark computational strategies, particularly in the early stages of virtual screening (VS). The use of in-silico screening techniques on libraries containing thousands of compounds enables the identification of potential ligands active toward a protein target, resulting in significant time and cost savings compared to traditional in vitro biological assays, which would be impractical on such a large scale. Optimizing computational workflows to ensure a high hit rate as early as the preliminary stage of VS is a crucial area of importance. Improving the quality of results in these early stages can, in fact, significantly increase the probability of success in the later stages of pharmaceutical development.

AutoDock Vina (often abbreviated just as Vina) is the most widely used and most integrated molecular docking algorithm within molecular docking software. Trott and Olson developed Vina at the Scripps Research Institute in California as an alternative to its predecessor, AutoDock 4. Key features that have contributed to Vina’s widespread adoption include its speed, improved accuracy over AutoDock 4, and free availability under an open-source license. Its high computational efficiency and ease of use led to very rapid uptake in the docking community, well evidenced by the high number of citations in the original article. AutoDock Vina is considered one of the fastest and most widely used open-source programs for molecular docking [6,7]. Despite its undoubted qualities and broad adoption, Vina has limitations. It has difficulties with increasingly flexible ligands, and its empirical scoring function may have a size-related bias [8,9]. GNINA represents a significant evolution in the field of molecular docking, building on Vina and its fork, Smina. The main innovation of GNINA is the integration of convolutional neural networks (CNNs) for scoring protein-ligand poses. In its typical workflow, GNINA samples ligand poses using Markov Chain Monte Carlo (MCMC) sampling, initially driven by AutoDock Vina’s empirical scoring function. Subsequently, GNINA uses the CNN-based scoring function to “rescore” and rank the poses obtained from MCMC sampling. Unlike empirical or knowledge-based scoring functions, which often assume a linear relationship between structural features and binding affinity, CNNs can model nonlinear relationships and potentially interpret molecular interactions in a more sophisticated way [7,10,11,12].

2. Results and Discussion

2.1. Energy Values and Conversions

Evaluations of the predicted energy from molecular docking were treated differently for the two algorithms analyzed because they yield different outputs. As for Vina, this is reflected in the output values of predicted energy, expressed in terms of the free energy of binding (ΔG), which is calculated in kcal/mol. The ΔG value was then converted to the corresponding binding affinity constant (K) value using the thermodynamic relationship reported in Equation (1).

K = e^(ΔG/RT)

(1)

Equation (1) was solved at a temperature of 310 K, with a gas constant R of 0.0019872036 kcal/(K × mol). Subsequently, the K value was transformed into its negative logarithm, pK.

GNINA returns several output metrics related to CNN assessment, namely: CNN score, CNN affinity, and CNN_VS. The CNN score is an assessment of the goodness of the generated pose, ranging from 0 to 1, where a score of 1 indicates a pose of higher reliability. The CNN affinity is an expected bond affinity and is expressed in pK. Finally, the CNN_VS is the product of the CNN score and CNN affinity and is used to rank compounds in virtual screening, having been found effective in retrospective evaluations [7,10,11,12].

2.2. Target Validation

The protein target validation phase involved identifying and selecting the most suitable protein model for use in the subsequent VS steps. Only protein models that met the following criteria were considered: (i) presence of a co-crystallized ligand; (ii) availability of an experimental binding affinity value (K_i or K_d); and (iii) crystallographic resolution of less than 3 Å to ensure an adequate level of structural detail. For each selected model, the docking algorithm’s ability to faithfully reproduce the pose of the co-crystallized ligand was evaluated by calculating the Root Mean Square Deviation (RMSD), as well as its accuracy in estimating the binding affinity by comparing the predicted inhibition constant with the experimentally determined value. In cases where the model had multiple co-crystallized domains, the analysis was conducted on each domain individually, and the one that showed the best performance in terms of structural and energetic accuracy was finally selected. Preliminary studies were performed using GNINA to select the PDB file and ultimately the domain to be used. The choice of GNINA at the PDB file selection stage was based on the fact that the CNN score calculated by the algorithm allows a distinction between receptors with high-quality binding sites and those with poor-quality binding sites. Therefore, receptors with the highest possible CNN score, but never below the threshold of 0.90 recommended in the literature, were selected. In the case of multiple domains present, the protein domain with the highest CNN score was selected [13]. The domains chosen and the CNN scores calculated for each target are shown in Table 1.

All selected protein targets were subjected to docking simulations using both Vina and GNINA to compare their performance directly. The results, presented in Table 2, demonstrate a significant performance improvement achieved with GNINA compared to Vina, in both structural and predictive accuracy.

GNINA produced RMSD values consistently below 2 Å for the predicted pose, indicating an excellent ability to reproduce the conformation of the co-crystallized ligand [14]. In contrast, Vina showed greater variability: in only 4 out of 10 targets, the RMSD value was less than 2 Å, while in the remaining cases it exceeded 3 Å, showing limited reliability in reproducing the co-crystallized pose. The variability in the redocking of Vina and the stability of GNINA can be observed in the two representative examples shown in Figure 1; the remaining comparisons are presented in Figures S1 and S2.

Similarly, prediction of binding affinity (expressed as pK) was also more accurate with GNINA in all targets examined, except in the case of Cyclin-dependent kinase 2, for which both algorithms returned overlapping values.

2.3. Virtual Screening Analysis

VS analyses were conducted using sets of ligands for which experimental affinity values (K_i or K_d) were available. The criterion adopted for selecting candidate molecules (“hits”) was based on the pK value estimated by the docking algorithm. Specifically, a selective filter was applied, considering only molecules with a pK ≥6.3, a threshold that corresponds to a binding affinity of approximately 500 nM, a value commonly accepted in literature as indicative of good biological activity. At the end of the filtering process, the predictive validity of the system was assessed by comparing the docking results with available experimental data. All molecules with experimental K_i/K_d values less than 1000 nM were considered effectively active. Based on this comparison, the percentage success rate of virtual screening was calculated according to the formula reported in Equation (2).

$$\mathrm{s}\mathrm{u}\mathrm{c}\mathrm{c}\mathrm{e}\mathrm{s}\mathrm{s} \mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e}(\mathrm{\%})={\displaystyle \frac{\mathrm{n}\mathrm{u}\mathrm{m}\mathrm{b}\mathrm{e}\mathrm{r} \mathrm{o}\mathrm{f} \mathrm{e}\mathrm{x}\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{i}\mathrm{m}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{a}\mathrm{l}\mathrm{l}\mathrm{y} \mathrm{a}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e} \mathrm{h}\mathrm{i}\mathrm{t}\mathrm{s}}{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l} \mathrm{n}\mathrm{u}\mathrm{m}\mathrm{b}\mathrm{e}\mathrm{r} \mathrm{o}\mathrm{f} \mathrm{h}\mathrm{i}\mathrm{t}\mathrm{s} \mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h} \mathrm{p}K \ge 6.3}}\times 100$$

(2)

This approach enables a quantitative evaluation of the algorithm’s ability to identify truly active molecules, providing an objective measure of the effectiveness of the screening protocol employed. As in the validation phase, the expected activity value considered in the VS analyses was CNN_VS for the GNINA algorithm and the pK value derived from the Vina ΔG transformation.

The data shown in Table 3 indicate that, in most cases, the use of GNINA resulted in a higher success rate percentage (SR%) compared with Vina. For 8 out of 10 targets analyzed, GNINA showed a greater ability to identify experimentally active ligands from the screening sets. A particularly significant example is Carbonic anhydrase II, for which GNINA identified 149 hits, 146 (98%) of which were active, compared with only two hits obtained from Vina, of which only one was active.

However, the difference was not marked in all cases. For the SYK kinase target, for example, Vina achieved a hit rate of 85%, slightly higher than that of GNINA (83%). Furthermore, for the Dopamine D3 receptor, both algorithms achieved essentially equivalent results, with a high predictive accuracy of 97%.

2.4. ROC Curves and Enrichment Factors Analysis

To statistically evaluate the ability of the two docking methods compared and obtain quantitative estimates of their sensitivity in distinguishing between active and inactive compounds, Receiver Operating Characteristic (ROC) curves with their Area Under the Curve (AUC) and Enrichment Factors (EF) were calculated.

ROC curves are an established graphical and analytical tool. In the context of virtual screening, an ROC curve is constructed by representing the true positive fraction (TPF) on the y-axis versus the false positive fraction (FPF) on the x-axis. Each point on the curve represents a TPF/FPF pair corresponding to a specific fraction of the molecular dataset classified according to the scoring function scores. A scoring function that can perfectly discriminate between active and inactive compounds (without overlapping score distributions) would have an ROC curve that passes through the upper left corner of the graph (TPF = 1, FPF = 0), indicating perfect sensitivity and specificity. In contrast, a scoring function that has no discriminatory power would produce a curve coincident with the 45° diagonal, corresponding to random selection (AUC = 0.5). The Area Under the ROC Curve (AUC) is a scalar value that summarizes the overall performance of a virtual screening. An AUC value closer to 1 indicates a high discrimination ability of the scoring function in distinguishing between active and inactive compounds over the entire data set. AUC > 0.5 suggests that the prediction ability is better than a random distribution model. Thus, AUC corresponds to the probability of correctly classifying a random pair of active ligand and decoy [15,16,17,18,19].

The Enrichment Factor (EF) is used to evaluate the ability of a scoring function or virtual screening method to enrich a subset of selected molecules (typically those with the best scores) in active compounds compared to a random selection from the entire data set. In practical terms, EF at a given percentage X% indicates how many times more active compounds are recovered in the top X% of the library than would be expected from a random selection of the same percentage of compounds. Performance evaluation often focuses on EF at low percentages of the ranked database, since the goal of virtual screening is to quickly identify a small fraction of promising compounds for experimental testing. For example, an EF at 1% (EF1%) equal to 10 means that 10 times more active compounds are found in the top 1% of compounds ranked by the scoring function than in a random selection of 1% of compounds from the entire data set [16,17,18,19,20].

Analysis of the ROC curves (Figure 2) clearly shows that GNINA outperforms Vina for all targets considered, except for the dopamine D3 receptor, for which the two methods exhibit overlapping performance. It is particularly relevant to note that GNINA achieved AUC values above 0.70 in almost all cases, never approaching the random classification threshold (AUC ≈ 0.50). In contrast, Vina exceeded the threshold of AUC > 0.5 in only six protein targets.

Analysis of EFs (Table 4 and Figure S3), assessed at the thresholds of 1%, 5%, and 10%, also confirmed the superiority of GNINA over Vina. Out of a total of 30 comparisons, GNINA showed better performance in 25 cases, and the two algorithms showed the same performance in 4 cases. Vina performed better than GNINA in only one analysis. Only in three circumstances – SYK kinase at EF1%, HSP90α at EF1% and dopamine D3 receptor at EF10% – does the performance of the two tools appear equivalent.

3. Materials and Methods

3.1. Protein Preparation

Protein structure selection was done by evaluating the RCSB PDB Databank (https://www.rcsb.org/) and the PDB-REDO Databank (https://pdb-redo.eu/) structure quality, and the best one was chosen considering crystallographic refinement and model quality parameters shown in the PDB-REDO Databank. The organism chosen was Homo sapiens, and the resolution was always less than 3 Å. Each protein model selected (Table 1) was prepared in YASARA software (v. 25.1.13, YASARA Biosciences GmbH, Vienna, Austria).

Hydrogens were added through ‘Clean > All’ YASARA’s option, waters and other non-necessary molecules for docking, and the co-crystallized ligand was removed using YASARA software. The protein was saved in PDB file format.

3.2. Dataset and 3D Structures Generation

Ligand datasets were downloaded from Binding Database (https://www.bindingdb.org/) using the protein name as a research query. Data were filtered by selecting only Homo sapiens as the target source organism, and ligands with known K_i or K_d values were chosen. Duplicate molecules were deleted, merging equivalent rows in DataWarrior (powered by openchemlib, v. 06.00.00) [21] using SMILES code as reference. Salts were also deleted using the largest fragment selection option in DataWarrior. Ligands’ 3D structures were generated and minimized in Open Babel (v. 3.1.1) [22], reflecting physiological pH (7.4) states. Then, 3D structures were optimized at a semi-empirical level using the xTB (extended tight-binding, v. 6.7.1) program package [23]. Geometry optimization was performed in xTB using the Analytical Linearized Poisson-Boltzmann (ALPB) model for water, with charge states specified for each molecule according to physiological pH conditions.

3.3. Molecular Docking Analysis

Molecular docking analysis was performed using both the molecular docking program GNINA (v. 1.3) [11] and the AutoDock Vina (v. 1.2.5) [24] docking algorithm. Regarding GNINA, we utilized the “rescore” docking mode setting, which involved 15 ligand pose rotations, and the poses were sorted by CNN score. The protein input format was in PDB file format, and the ligand input format was in SDF file format. AutoDock Vina docking analysis was performed by setting the exhaustiveness parameters to 32 and converting the input files to the PDBQT file format required by Vina.

In both cases, the simulation boxes were built using AutoDock Tools (v. 1.5.7), with the co-crystallized ligand pose serving as the reference. Grid parameters – shown in SI – were chosen wide enough not to force the ligand-receptor interaction.

3.4. RMSD Calculation

RMSD calculation was performed using the tool DockRMSD (v. 1.0.0) [25]. The co-crystallized ligand was converted into a MOL2 file format through Open Babel, as required by the tool. The docking output poses were merged into a single pose, and the first one was converted into a MOL2 file format. The co-crystallized ligand pose and the first docked one were used as input to calculate RMSD.

3.5. Decoys Generation

Decoys were prepared via the LUDe web server [26] from the 20 most active ligands for each target, which constituted the true positives. Of the decoys generated by LUDe, 400 were selected for each target, which constituted the false positives. Decoys are decoy molecules constructed from ligands with known experimental activity toward the protein target, having similar physicochemical properties but different 2D topology. Decoys have been prepared as reported in paragraph 3.2.

3.6. ROC Curves and Enrichment Factors

ROC curves and related AUC values were calculated using the roc_curve() and auc() functions of the sklearn.metrics module (scikit-learn v. 1.3.0). A set of 420 compounds was prepared for each protein target, including the 20 most active ligands according to experimental affinity values and 400 decoy molecules. The active ligands represented True Positives, while the 400 decoys represented False Positives. False Positive Rate (FPR) and True Positive Rate (TPR) values were obtained by comparing the binary labels with the predictive scores of each docking method. The AUC was subsequently calculated by trapezoidal numerical integration of the ROC curve. ROC curves were generated using matplotlib (v. 3.7.1), plotting TPR vs FPR for each method with indication of the respective AUC values in the legends. The diagonal reference line (random classifier) was included for visual comparison of performance.

EFs were calculated at the thresholds of 1%, 5%, and 10% of the total dataset. For each predictive method, the dataset was sorted in descending order of score. The number of active compounds in the top-ranked subset was divided by the expected number of active compounds in a random selection of the same size, according to the formula reported in Equation (3).

$$\mathrm{E}\mathrm{F}={\displaystyle \frac{(\mathrm{N}\_\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e}\_\mathrm{t}\mathrm{o}\mathrm{p}/\mathrm{N}\_\mathrm{t}\mathrm{o}\mathrm{p})}{(\mathrm{N}\_\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e}\_\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}\mathrm{s}/\mathrm{N}\_\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}\mathrm{s})}}$$

(3)

where N_active_top represents the number of active compounds in the selected subset, N_top the size of the subset, N_active_totals the total number of active compounds in the dataset, and N_totals the total size of the dataset. Figure S3 was generated with the online Screening Explorer Tool (http://stats.drugdesign.fr/, accessed 05/31/2025) [27].

4. Conclusions

The present study compared the performance of Vina with that of GNINA, which integrates a CNN-based system used for rescoring poses obtained from molecular docking. Vina was chosen as the comparison algorithm because it is widely used and integrated into standard molecular docking software. Therefore, the primary application we aimed to evaluate is a retrospective analysis of VS. We sought to determine whether the ability of CNNs to interpret molecular interactions provides advantages in the search for active compounds when performing VS of large libraries of molecules.

GNINA showed higher accuracy in reproducing the co-crystallized ligand pose and predicting protein binding affinity compared to Vina. RMSD values consistently less than 2 Å were recorded in all analyses performed with GNINA, while Vina showed greater variability with RMSD values often exceeding 3 Å. In predicting binding affinity, GNINA also returned more accurate values than the experimental values reported in the literature.

In VS analyses, GNINA demonstrated a greater ability to identify experimentally active ligands than Vina, achieving excellent success rates in most cases that were superior to those of Vina. A noteworthy result is that of Carbonic anhydrase II, where GNINA recorded several hits almost exclusively populated by experimentally active ligands. To reinforce these observations, more accurate statistical investigations were conducted by generating a decoy library for each protein target. ROC curves were generated, and respective AUC values – an indicator of the algorithm’s ability to discriminate between active and inactive compounds – were calculated. For GNINA, consistently high AUC values were recorded (AUC >0.70) and far from the threshold of AUC ≈ 0.50, which would indicate random classification. Vina results were less homogeneous, with values often below the AUC threshold. EF at different rates of 1%, 5%, and 10% was also evaluated to assess the effectiveness of the methods in recovering active compounds from the library fractions. Again, there were results in favor of GNINA, which consistently outperformed Vina.

In summary, we can state that the use of GNINA CNNs in the rescoring phase of the poses leads to benefits in performance in all phases of VS: from pose reproducibility to prediction of binding affinity and discriminatory ability against active compounds. Therefore, in future computational drug searches, we recommend adopting the workflow for VS discussed in this study, which includes a rigorous target validation phase and the use of the GNINA algorithm for VS of compounds, to maximize its effectiveness and significantly increase the chance of success in identifying promising compounds. The integration of AI-based algorithms allows overcoming the limitations of traditionally employed tools, offering new perspectives for in-silico approaches used in drug discovery.