We evaluated eight docking programs available on a free license, so our work can be reproduced without monetary expense: ADFR, v. 1.1 (Ravindranath et al., 2015), AutoDock Vina, v1.1.2 (Trott and Olson, 2010), UCSF DOCK, v6.9 (Allen et al., 2015), GNINA, v1.0.3 (McNutt et al., 2021), jdock, v2.2.3c (Li et al., 2012), PLANTS, v1.2 (Korb et al., 2009), RxDock, v0.1.0 (Ruiz-Carmona et al., 2014), and smina, v1.1.2 (Koes et al., 2013). The selection includes the two most widely used packages: UCSF DOCK and AutoDock Vina. UCSF DOCK and RxDock were downloaded as source code and compiled locally. The other programs are available as precompiled executables. ADFR, AutoDock Vina, jdock, smina, and RxDock return an affinity in kcal/mol; UCSF DOCK gives a gridscore, while GNINA returns three values: GNINA affinity in kcal/mol, an adimensional CNN score (0–1), and a CNN affinity [−log (K)].

The binding target was the human cGMP-specific 3′,5′-cyclic phosphodiesterase (PDE5A, UNIPROT ID: O76074) (Zhang et al., 2004), a metalloenzyme containing zinc and magnesium. There are 38 X-ray diffraction (XRD) PDE5A structures deposited in RCSB, four crystalized with sildenafil (PDBID: 2H42, 1UDT, 3JWQ, and 1TBF). An AlphaFold model was obtained using the PDE5A sequence from UNIPROT through UCSF ChimeraX, v1.7. All 38 PDE5A structures and the AlphaFold model were evaluated, and structures 1T9R, 1T9S, and the AlphaFold model were tested in a screen against a ligand database.

Because ADFR, Vina, GNINA, jdock, and smina use pdbqt receptors and ligands, both sets were prepared using the scripts included with ADFR (prepare_receptor.py and prepare_ligand.py). UCSF DOCK and PLANTS require mol2 files; receptors for UCSF DOCK were prepared with UCSF Chimera (Pettersen et al., 2004). For PLANTS, receptors were processed with SPORES, v1.3 (Brink and Exner, 2010).

Ligands were obtained from UCSF ZINC 12: a total of 2115 FDA-approved drugs, including stereoisomers (Irwin et al., 2012), were screened, that is, they were all docked to each of our three receptors (1T9R, 1T9S, and the AlphaFold model). Other than conversion to mol2 or pdbqt, no other preprocessing was performed. PLANTS failed with four ligands; jdock failed with 41, and ADFR failed with one. Failure reasons were not identified or corrected. A total of 513 experimentally determined binders for PDE5A were obtained from UCSF ZINC by selecting the results with affinities better than 10 nM. These were used as true positive results for ROC calculation.

Ligands were converted into mol2 using Racoon, v1.0f (Forli et al., 2016). RxDock requires sd files for ligands; these were prepared with Open Babel, v3.1.0 (O’Boyle et al., 2011), while the receptor is a mol2 file prepared with UCSF Chimera, v alpha 1.19. ADFR requires an additional step to create a .trg file.

File formats for the ligands and receptors, as well as instructions to prepare them, are summarized in Supplementary Table S1.

The coordinates and size of the search were defined using the largest cavity in the PDE5A models, found by AGFGRUI and used for all docking programs except UCSF DOCK and RxDock. For Dock and RxDock, the search box was generated using GMP (called 5 GP in the PDB) as a reference. Search coordinates were the same for all systems, and grid sizes are specified in Supplementary Table S2. All programs were run with default values for the number of results and other parameters (see supplemental material Command Line for specific details). Running time for docking ranged from 4 h to 18 h for each program. Only the enhanced searches (using ADFR, Vina, and GNINA) required longer running times. All docking was conducted without Zn because only PLANTS can handle Zn atoms. All ligands and water molecules in XRD structures were removed before docking experiments.

ADFR, Vina, and PLANTS were executed in the “Laboratorio Nacional de Supercómputo del Sureste de Mexico” facility; GNINA, jdock, RxDock, smina, and UCSF DOCK were run on a workstation under Debian 11. An RTX 3060 GPU was used for GNINA. The most recent GNINA version, v1.3, can be run on CPUs.

Most docking programs return several results per input ligand. For analysis, we used only the one with the highest binding energy.

ROC calculations were performed using R (R Core Team, 2021) and RStudio (RStudio, 2020) using libraries pROC (Robin et al., 2011) and cutpointr (Thiele and Hirschfeld, 2021). All relevant files (initial files and results) are available at figshare with doi: 10.6084/m9.figshare.24217872. Ligand RMSD was calculated using DockRMSD (Bell and Zhang, 2019).

First, we compared the 38 available PDE5A XRD structures and the AlphaFold model through ensemble docking. Four structures could not be used (1TBF, 3HC8, 3HDZ, and 6L6E) because they were not successfully converted to pdbqt or otherwise processed through agfrgui; this was due to sidechains with multiple occupancies. No corrections were attempted due to a lack of a clear heuristic, resulting in 35 target structures. Ensemble docking with sildenafil was performed for these structures with the docking packages that allowed it (Figure 1A). UCSF DOCK, PLANTS, and RxDock were excluded due to command line limitations. ADFR produced the highest affinities, followed by Vina, smina, and GNINA, with similar performance, which was expected because the latter two are based on Vina. jdock returned the same affinity for sildenafil on every receptor, suggesting that it is not suitable for ensemble docking (Figure 1A). The ensemble experiment results are shown in Figure 1B, graphed by the receptor used. The figure shows that receptor selection impacts docking results. Structures displayed different binding affinities for sildenafil, with a small standard deviation of 0.58–0.75 kcal/mol (Figure 1C). The best affinity was for receptor 2H44 (* in Figure 1B), a structure determined in the presence of icarisid II, that displayed high affinity when using smina, Vina, and ADFR (Figures 1B,C). The three structures crystallized with sildenafil (vertical lines) also showed high affinities for sildenafil (Figures 1B,C).

To perform a virtual screen against all the FDA-approved drugs (screen described in the next section), we selected three PDE5A structures that represent the conformational diversity of the binding site: 1T9S, a chimeric structure crystalized with the enzyme’s catalytic product GMP that could bind inhibitors like sildenafil and similar molecules; 1T9R crystallized with an occluded, empty binding site; and the AlphaFold model that we created for PDE5A to simulate a docking study where targets have no previous crystallographic reference. 1T9R was the most different, with an RMSD of ∼3.5 Å compared to the other two (Supplemental Table 3) due to its “H-loop” (residues L672-H685), located over the GMP binding site (Zhang et al., 2004). We did not select a structure crystallized with sildenafil because most docking campaigns for repurposing would not have such a reference available. In the ensemble experiment, GNINA, smina, Vina, and ADFR docked sildenafil to 1T9S (†) and to the AlphaFold model (‡) with affinities close to the ensemble average (Figure 1C). 1T9R, with a closed H-loop over the binding site, displayed less affinity on all the programs evaluated (Figures 1B,C).

GNINA calculates not only a “GNINA affinity” but also a CNN score, which grades the quality of docking results through network-trained image recognition algorithms. GNINA also produces a second affinity, based on the CNN model (“CNN affinity”). In Figure 1D, we graphed each receptor against its CNN score for sildenafil from the ensemble docking in GNINA. The three structures crystalized with sildenafil (vertical lines) had higher scores than an arbitrary CNN score cutoff of 0.9 that we hereby defined as “high-quality docking.” 1T9S also scored as high quality, but the CNN score was lower than this threshold for the AlphaFold model and 11 other structures. Visual examination of receptors below this 0.9 CNN score cutoff reveals receptors with occluded binding sites or missing residues. The worst-scoring was 1T9R, likely due to its occluded binding site. These results reaffirm our selection of 1T9S, 1T9R, and the AlphaFold model to represent the conformational diversity of PDE5A’s binding site.

Next, we docked a dataset of 2,115 approved drugs on the three selected PDE5A structures (1T9S, AlphaFold model, and 1T9R) using each of the eight docking packages. We ranked results by affinity from highest to lowest, then looked for the known PDE5A inhibitors: sildenafil and vardenafil, which are chemically similar, sharing a guanosine-like core with GMP, and tadalafil, avanafil, and dipyridamole, which diverge structurally (Figure 2). The ordinal position of these inhibitors in the ranking was compared (Figure 2). For 1T9S, UCSF DOCK, and RxDock detected sildenafil in the 44th and 46th place (Figure 2, 1T9S). CNN affinity picked sildenafil at rank 101, better than rank 143 by GNINA’s affinity. In addition, vardenafil and tadalafil were ranked in the top 100 best results by several packages. Tadalafil was 16th by GNINA, 51st by jdock, 69th by Vina, and 88th by smina. Vardenafil was 51st by UCSF DOCK and 64th by CNN affinity. Cartesian minimization after docking with UCSF DOCK improved sildenafil’s rank to 39 and vardenafil’s to 48.

When the AlphaFold model was used as a receptor, UCSF DOCK identified sildenafil as the 58th best ligand. No other program ranked it better than 115th. Tadalafil ranked remarkably well, first by JDock, second by GNINA and Vina, and seventh by smina, while vardenafil was 64th with UCSF DOCK (Figure 2, AlphaFold). CNN affinity did not improve the ranking of known ligands in the AlphaFold model.

When docking to 1T9R (closed binding site), the best ranking for sildenafil was 106, with ADFR (Figure 2, 1T9R). ADFR picked vardenafil in place 37. Avanafil was 66th by CNN Affinity and 35th by PLANTS. As expected, 1T9R did not find sildenafil in the top ligands, and known PDE5A ligands consistently ranked worse than in the other two target structures.

Up to this point, all programs had been executed using moderate screening settings for exhaustiveness and number of runs. For 1T9S, an additional set of dockings was produced with a ten-fold increase in runs at high exhaustiveness, with the three programs that supported it, ADFR, GNINA, and Vina (Figure 2, 1T9S x10). Generally, rankings for the known PDE5A inhibitors were worse with more exhaustiveness, while CNN affinity remained unaffected. Thus, more computational power and time did not improve the identification of known binders.

We used the CNN score produced by GNINA to extract the best poses according to the arbitrary ≥0.90 threshold for high-quality results presented in Figure 1D. Only the poses that surpassed this quality threshold were selected and then ranked by regular GNINA affinity or by CNN affinity. This heuristic, incorporating the CNN quality score, did not find sildenafil or any other known PDE5A inhibitors as a hit for the AlphaFold model or for 1T9R (Figure 2, CNN score). In contrast, for 1T9S, it ranked sildenafil and vardenafil better than when using docking affinity: sildenafil was 26th and 101st, while vardenafil was 52nd and 64th using GNINA’s affinity and CNN affinity, respectively (Figure 2, see row: “CNN score”). Figure 1D shows that the CNN score reveals a useful distinction between receptors with poor (1T9R) and high-quality (1T9S) binding sites. Our AlphaFold model is closer to high quality, yet critically below the 0.9 CNN score threshold. Thus, when a ligand that is a known true positive (like sildenafil) is available to evaluate a receptor, the CNN score is a useful metric to select a receptor to find other true positives. Once a good receptor is selected, the CNN score can be used to further filter docking screening results, improving their ranking.

Receiver operating characteristic (ROC) curves have been used to evaluate docking performance, so we used them to compare the three target structures and eight docking programs analyzed here. To create the ROC analyses, we docked 513 experimentally determined binders (all curves are shown in Supplementary Figure S1) and the set of 2,115 FDA-approved drugs into each structure (1T9S, AlphaFold model and 1T9R) with each program (eight suites plus the CNN affinity from GNINA) resulting in 27 initial ROC curves (Table 1). Overall, the AUC values and Youden’s index were better for 1T9S (Table 1). The best AUC values were found with CNN affinity, PLANTS, and ADFR for all targets. Accuracy, precision, and sensitivity displayed acceptable performance across all programs. However, specificity, that is, the portion of true negatives correctly identified, was low (under 0.5) except when using CNN affinity. Cutoff points represent the best method of distinguishing true positives from true negatives, and they were higher for 1T9S, but values are dependent on the scale used by each docking package and were thus not directly comparable between them. CNN affinity cutoff values are positive because they correspond to −log(K).

The heuristic described in the previous section, which incorporated a CNN score threshold of ≥0.9 before affinity ranking, was used to produce another four ROC analyses, two for 1T9S, and two for the AlphaFold model, using only the results that passed the quality threshold. The use of this threshold produced a large increase in specificity (0.15–0.97 and 0.53–0.92 for 1T9S; 0.15–0.95 and 0.46–0.96 for the AlphaFold model) along with a small decrease in sensitivity (0.95–0.92 and 0.95–0.92 for 1T9S; 0.94–0.81 and 0.94–0.81 for the AlphaFold model) with GNINA affinity and CNN affinity. Youden’s index and the AUC also increased. This heuristic, based on the CNN score, could not be used for 1T9R because no docking results passed the quality threshold. These results highlight that docking affinities alone are poor classifiers that can be supplemented with information about quality, such as the CNN score.

Superposition of the structures evaluated here (1T9S, 1T9R, and the AlphaFold model) with a counterpart crystallized with sildenafil (2H42) shows conformational changes in PDE5A’s H-loop (Figure 3, black ovals) that explain ligand accessibility. The H-loop is extended in the wild-type sildenafil-bound 2H42 (orange) (Figures 3E,J). In 1T9S (green), this region forms two small helices away from the catalytic site (Figures 3B,G), allowing ligand access; while in 1T9R (purple), histidine 678 on the extended H-loop seems to block ligand access (Figures 3C,H). With GNINA, sildenafil docked into 1T9S in a pose remarkably like the XRD determined (Figure 3G vs. J), with a resulting RMSD of 1.04 Å when excluding the ethylpiperazine motif (Figure 2, red square) and 3.07 Å when including it. This motif is stabilized by the H-loop in 2H24 but is solvent exposed in our docking to 1T9S. PDE5A contains numerous water molecules, magnesium, and zinc trapped in the binding site (Figure 3F). Their removal was compatible with successful docking with GNINA, aided by the CNN score, and produced a sildenafil pose very close to the crystallographic. The AlphaFold model has an accessible binding site (Figure 3D) like 1T9S and docks sildenafil in a similar site (Figure 3I), yet docking results were better for 1T9S. In turn, docking into 1T9R happens at a different site (Figure 3H). By comparing the parameters displayed in Figures 3G–J, we propose that a high CNN score identifies correct poses for potential binders that have high affinity and low RMSD (Figures 3G–J).