Vecura Biotech Insiders Spotlight #12: Building and Optimizing a Virtual Screening Pipeline for HIV-1 Protease: Lessons from Real-World Benchmarking with Darunavir

Standard drug-likeness filters threw Darunavir out of the screen. Darunavir is an FDA-approved HIV-1 protease inhibitor, and it was sitting in the very crystal structure the campaign was docking into. Relaxing the filters to keep it repaired the benchmark and compromised it in the same move: a pipeline whose thresholds were tuned around one approved drug cannot then use that same drug as proof it works. So two more approved inhibitors went in, ones the pipeline had never seen, and were allowed to rank wherever they fell. One came fourth of eight. The other came last.

Why HIV-1 protease, and why an approved drug is the right thing to check against

HIV-1 protease cuts newly made viral proteins into the pieces the virus needs to assemble new copies of itself. Block the enzyme and the virus can still infect a cell but cannot finish building anything infectious, which is why protease inhibitors have been core HIV treatment for decades.

The enzyme is a homodimer, and the structure used here, PDB 2IEN at 1.30 Å resolution, already has Darunavir solved into its active site as ligand 017. That makes it useful twice over. It is exactly the kind of structure a screening campaign would dock new candidates into, and the right answer is already in the file, so every stage of the pipeline can be measured against a molecule known to work instead of trusted on faith.

Viet Nguyen Duc set the screen up around that structure: catalytic dyad at Asp25 on chain A and Asp25 on chain B, chain B residues mapped across 101 to 199, docking grid centered on the co-crystallized ligand at (15.733, 22.437, 16.919) with a 19.0 x 18.8 x 15.5 Å box. The library was 4,719 active compounds pulled from ChEMBL.

Figure 1. The HIV-1 protease homodimer (PDB 2IEN) with the predicted binding pose of the calibration control, Darunavir, in the active site.

Three controls, and why the count matters

The campaign carried three approved inhibitors, in two different roles.

• Darunavir (ChEMBL1323) is the in-run calibration control. It is the molecule the filters were tested against and adjusted around.

• Atazanavir (ChEMBL1163) and Lopinavir (ChEMBL729) are held-out controls. They were kept out of every calibration decision and then screened like any unknown in the library.

That split is the structural fix for a problem most pipelines never name. If a filter is loosened because it rejected Darunavir, and the pipeline is then declared validated because it recovers Darunavir, the argument is circular and the recovery was guaranteed. A compound the thresholds were never tuned on cannot be recovered by construction, so whatever happens to it is information.

Stage 1. The filter that rejected the answer

The first pass applied conventional rule-based filters, Lipinski's Rule of 5 and Veber's Rule, the kind almost every screening pipeline opens with. They eliminated Darunavir itself: 11 rotatable bonds and a polar surface area of 140.4 Ų, both outside Veber's usual cutoffs.

That is not a quirk of one molecule. Later-generation protease inhibitors tend to be larger and more polar than the rule of thumb assumes. A filter that discards the reference compound does not just lose one row, it invalidates the benchmark, because there is no longer a known-good answer in the pool to check anything else against.

The fix was to demote the Veber criteria from a hard cutoff to an advisory flag. Compound survival more than doubled, from 14.6% to 30.7%, leaving 1,449 of the original 4,719 ChEMBL hits in play, Darunavir included.

Stage 2. Toxicity cutoffs measured against drugs that already ship

ADMET prediction ran into a sharper version of the same problem. Applying ADMET-AI's raw thresholds directly, failing anything at or above 0.5 on predicted hERG or DILI risk, eliminated nearly all 1,449 surviving compounds.

Benchmarking against FDA-approved drugs in DrugBank showed why. Approved HIV protease inhibitors as a class sit at the 61st percentile for predicted liver injury, with 39% of all approved drugs scoring above 0.5, and at the 94th percentile for CYP3A4 inhibition. The second of those is not a flaw at all but a known property of the class, which is why protease inhibitors are routinely co-dosed with a CYP3A4-inhibiting booster such as ritonavir.

If the drugs already on the market score worse than your cutoff on a given risk, the cutoff is measuring the wrong thing for this target class.

So the rule was rewritten in relative terms: a candidate fails on toxicity only if its predicted risk exceeds the 90th percentile of approved drugs within the same target class. Those percentiles were computed from eight approved HIV-1 protease inhibitors, Darunavir, Saquinavir, Ritonavir, Indinavir, Nelfinavir, Amprenavir, Tipranavir and Fosamprenavir. Atazanavir and Lopinavir were deliberately left out, so neither held-out control had any hand in setting a cutoff it would later be judged by. Of the 1,449 Veber-compliant compounds, 68 cleared every calibrated ADMET criterion.

Stage 3. AutoDock Vina disagrees with itself

With a workable pool in hand, Viet Nguyen Duc ran a check that is easy to skip: dock the same molecules twice, under identical settings, and see whether the scores come back the same.

They do not. Across duplicate runs the mean run-to-run difference was 1.06 kcal/mol, with a peak of 2.30.

Table 1. Run-to-run noise on repeated Vina docking of the same molecules under identical settings.

The row worth sitting with is the last one. Darunavir, the co-crystallized reference, swung 1.1 kcal/mol between runs, which places the calibration standard among the least reproducible compounds in the set. A molecule that is known to bind, whose pose is in the crystal structure, moving that much under stochastic sampling is the clearest possible argument that a single Vina run cannot be asked to rank a top tier.

The operational consequence was applied across the whole campaign: duplicate runs for every final candidate, scores reported as mean and spread, and any two compounds separated by less than 1.06 kcal/mol treated as statistically tied rather than ranked against each other.

Stage 4. An orthogonal check, and the bar it had to clear first

The answer to Vina's noise was not to trust Vina harder. It was to bring in a method built on different principles and see where the two agreed. GNINA's CNN rescoring ran across the 68 ADMET-safe hits plus Darunavir as an in-run reference, using the default crossdock_default2018 model, exhaustiveness 16, 9 modes per ligand, a fixed seed of 42, rigid receptor docking and the same grid box, in roughly 17.4 minutes of compute for the full set.

Final ranking was done on GNINA's own built-in composite, CNN_VS, which is CNNscore multiplied by CNNaffinity: a pose probability between 0 and 1 against a predicted affinity in pKd units. Pose probability and predicted strength have to agree before a compound rises.

Before that metric was pointed at unknowns it had to clear its own bar. Redocking Darunavir reproduced the crystal pose at 0.45 Å centroid RMSD, comfortably inside the 1.0 Å pass threshold. Centroid shift rather than atom-by-atom RMSD is the right measure for a symmetric homodimer like 2IEN, where equivalent atoms in the two chains can be numbered either way round and a strict atom mapping invents error that is not there. Separately, a QSAR model (Morgan fingerprints with Random Forest and XGBoost, 5-fold cross-validated R² = 0.600, MAE = 0.721) predicted Darunavir's potency at pIC50 8.95 against a measured 9.25, an error of 0.30, with Darunavir held out of training through out-of-fold predictions so the model was blind to its scaffold.

Pose strain was filtered as well. Poses with intramolecular energy above +1.0 kcal/mol were deprioritized as too distorted to be physical, which is how several otherwise high-affinity poses for Lopinavir, poses 7 through 9, came to be rejected on internal energy alone.

Figure 2. GNINA CNN rescoring in Vecura, showing the top binding pose of the control, Darunavir. Darunavir appears as "Lig.1" here because it was analyzed in an isolated verification run, separate from the multi-ligand batch in which Atazanavir occupied the first index. The interface reports an internal raw Vina energy of -9.00 kcal/mol for this pose, matching the -9.00 ± 0.55 kcal/mol mean in Table 1.

Where the pipeline landed

Table 2. Final composite ranking across the two held-out approved inhibitors, the in-run control and the five surviving library leads. Tier-1 means scoring above Lopinavir, the highest-ranking held-out approved drug; Tier-2 falls below it.

Darunavir tops the list at 10.58, which is what a functioning reference is supposed to look like once the filters stop fighting it. Two compounds out of the 68-compound library, ChEMBL292043 and ChEMBL301213, place second and third, above Lopinavir, and both are ADMET clean under the calibrated thresholds. ChEMBL301213 carries the highest model confidence in the whole set at 0.98 along with comfortably negative strain. Three more library compounds land in Tier-2, below Lopinavir but still ahead of the other held-out control.

One entry needs a caveat attached to it permanently. ChEMBL292043 is ranked on its fourth pose rather than its first, because its 9.98 pKd affinity shows up in that pose alone and nowhere else in its ensemble, paired with a slightly positive intramolecular energy of +0.26 kcal/mol. It was carried forward precisely because that pattern is ambiguous: it is either a genuinely novel binding mode or a CNN artifact, and the campaign's own next step is designed to tell the two apart.

The control that finished last, and why it is in the write-up

Atazanavir is an approved HIV-1 protease inhibitor. Held out of every calibration step and screened like an unknown, it placed eighth of eight.

The two scoring methods split sharply on it. Raw Vina liked it: a mean of -9.93 ± 0.53 kcal/mol, with pose 7 reaching -10.71, strong by any standard in this campaign. The CNN did not, assigning scores of only 0.48 to 0.75 across its poses, 0.75 at pose 1. Multiplied out, that gives a composite of 7.19, below all five library leads and well below the other held-out control.

Two explanations are on the table, and neither is flattering to the pipeline. Atazanavir has a bulky, flexible peptidomimetic scaffold, and docking it into a rigid 2IEN receptor may have forced it into strained geometries that the CNN then read as physically unrealistic even while Vina was reporting favorable interaction energies. The other possibility is simpler: CNN_VS carries noise and scaffold bias of its own, and does not recover every known active.

What the result rules out is using CNN_VS as an absolute exclusion criterion. What it requires is molecular dynamics with a flexible receptor before any candidate goes near a bench. An approved drug coming last is not a validation failure that needed to be smoothed over, it is the one outcome a held-out control exists to produce, and a pipeline that prints it is worth more than one that only ever recovers what it was tuned on.

What these numbers do not support

• No GNINA CNN_VS replicates were run for the final scoring phase, so the composite metric's own noise floor is unknown. The gaps that set the tiers are narrow: rank 2 and rank 3 differ by 0.03 composite units, and the Tier-1 to Tier-2 boundary rests on margins as small as 0.14. The tier assignments should be read as indicative rather than discriminating.

• The 90th-percentile toxicity cutoff was computed from a class of eight approved inhibitors, which means in practice the threshold is effectively set by a single molecule.

• Redocking Darunavir into its own cognate crystal structure at 0.45 Å is the easiest validation scenario that exists. It is a baseline sanity check passed, not evidence of predictive performance on novel ligands.

What runs next

• Use Darunavir's composite score of 10.58 and CNN affinity of 10.91 pKd as the contextual reference point against which any new hit is judged.

• Pull the top poses for ChEMBL292043, ChEMBL301213, ChEMBL305125 and ChEMBL294331, alongside Darunavir, into PyMOL or Chimera and verify real hydrogen bonding with the catalytic dyad (Asp25 and Asp25'), hydrophobic contact with the flap residues (Ile50 and Ile50'), and occupancy across the S1, S1', S2 and S2' sub-pockets.

• Re-dock ChEMBL292043 at exhaustiveness 32 with 20 modes, to settle whether its standout single-pose affinity is a real global minimum or a sampling artifact.

• Analyze the rank shifts between raw Vina energy and the CNN output across the leads, and identify which structural features, aromatic stacking or directional hydrogen bonding, drove compounds up.

• Run 50 to 100 ns explicit-solvent molecular dynamics on the leads and on Darunavir, assess trajectory stability and flap dynamics, and quantify binding free energies with MM-GBSA or MM-PBSA.

• Filter the surviving scaffolds for PAINS, reactive groups and aggregation risk before anything is procured or synthesized.

Scope

• This campaign is computational end to end. Nothing was synthesized or assayed, and no HIV-1 protease inhibition has been measured in a biological system.

• All rankings are built from predicted values: docking energies, CNN pose scores, QSAR potency and calibrated ADMET flags. They describe a prioritization, not measured potency or safety.

• The follow-up work listed above, pose inspection, high-exhaustiveness redocking and molecular dynamics, had not been run at the time of writing.