Counterscreening Protocol

Step 1: Target Identification and Validation

The input identifier for the desired target is validated against three databases in priority order: RCSB PDB, UniProt, and HGNC. This normalization step determines whether the identifier is a PDB ID, a UniProt accession, a UniProt entry name, or a gene symbol, and confirms that it resolves to a known, unambiguous entry. Validation before any download prevents silent failures where a mistyped or deprecated identifier would propagate through the pipeline and produce meaningless results.

Step 2: Protein Metadata Download

Protein metadata is fetched from UniProt. This step resolves cross-references between databases and records gene names, organism, functional annotations, and known isoforms. Having this metadata before structure download allows the pipeline to select the biologically relevant chain and to flag cases where multiple orthologues or isoforms exist — information that is particularly relevant in counterscreening, where isoform context is often central to the selectivity question.

Step 3: Protein Structure Acquisition

The 3D structure is obtained in one of two ways depending on the input identifier.

Experimental structures (RCSB PDB): The pipeline fetches the PDB record and chain information, downloads the mmCIF file alongside the PDB file, converts mmCIF to PDB when no PDB file is available using the Gemmi library, and splits multi-model structures (e.g., NMR ensembles) into individual conformations for independent evaluation. It also queries the OPM (Orientations of Proteins in Membranes) database to flag membrane proteins, and fetches binding site annotations from NCBI for each chain.

Starting from an experimental structure is strongly preferred when one is available. X-ray crystal structures and cryo-EM structures capture the true three-dimensional arrangement of atoms in the binding site, including water molecules, metal ions, and co-crystallized ligands. In a selectivity context, even small structural differences between an isoform pair — a single residue substitution, a shifted loop, a different bound cofactor — can be decisive, so structure quality has a disproportionate effect on the reliability of selectivity predictions.

Predicted structures (AlphaFold DB): For targets without an experimentally determined structure, the pipeline fetches the predicted structure from the AlphaFold Protein Structure Database.^[1] Per-residue confidence scores (pLDDT) are available and can flag regions where the predicted geometry is uncertain — particularly disordered loops that may border the binding site.

Alternatively, users may upload their own structure file directly in PDB or mmCIF format. Uploaded structures bypass the download steps and enter the pipeline at the preparation stage.

Step 4: Receptor Preparation — Fixing, Reducing, and Protonating

Raw PDB files from public databases are almost always unsuitable for docking in their downloaded form. The problems are systematic and, if left unaddressed, would produce artifactual selectivity scores:

• Missing hydrogens. X-ray crystallography does not resolve hydrogen atoms at typical resolutions. Hydrogen atoms are essential for docking: they participate in hydrogen bonds, alter the shape and electrostatic surface of the binding site, and are required by AutoDock-family scoring functions to compute partial charges.
• Missing atoms and residues. Flexible loops and termini are often disordered in crystal structures. If a missing loop borders the binding site, leaving it absent creates an artificial opening that distorts the pocket geometry.
• Non-standard residues. Modified amino acids and crystallization artifacts must be converted to their standard equivalents or removed, because force fields do not have parameters for them.
• Incorrect protonation states. The protonation state of titratable residues (histidine, aspartate, glutamate, lysine, arginine, cysteine, tyrosine) depends on the local pH and electrostatic environment. The wrong protonation state at the active site can shift docking scores by several kcal/mol. In selectivity analysis, a protonation error on one protein but not the other would produce a spurious apparent selectivity, making this step particularly important.

The preparation pipeline addresses all of these with a four-stage protocol:

PDBFixer is applied first. It adds missing heavy atoms using template-based modeling, resolves missing terminal residues, replaces non-standard residues with their standard equivalents, removes heterogens, and performs an initial coarse hydrogen placement.

Reduce (TRIM mode) is then applied to optimize hydrogen placement for key titratable and ambiguous groups.^[2] This is especially important for histidine (which can be protonated at Nδ, Nε, or both) and for asparagine and glutamine (whose amide groups can be flipped 180° without changing heavy-atom positions but dramatically changing hydrogen-bond geometry). The TRIM variant removes the hydrogens it placed while preserving their optimized orientations, avoiding double-protonation in the following step.

PDB2PQR with the SWANSON force field is the core protonation engine.^[3] It assigns protonation states at pH 7.4 using empirical pK_a models and adds all hydrogen atoms accordingly. The SWANSON force field assigns partial charges compatible with the AutoDock/Vina scoring function; using charges parameterized for a different force field would create a systematic mismatch that degrades scoring accuracy. PDB2PQR outputs a PQR file — a variant of PDB format that includes per-atom partial charges and radii required by the docking engine.

OpenBabel converts the PQR file back to standard PDB format and then to PDBQT format, the native input format for AutoDock-family docking engines, encoding torsional flexibility information for ligands through BRANCH/ENDBRANCH annotation blocks.

Step 5: Receptor Energy Minimization

Even after fixing and protonation, the protein geometry carries artifacts from template-based atom additions and protonation state changes. These manifest as local steric clashes and bond strain that would bias docking results. The receptor_minimization pipeline applies energy minimization using OpenBabel with the MMFF94 (Merck Molecular Force Field).^[4] The minimization runs a two-stage algorithm — steepest descent first to escape large energy barriers quickly, followed by conjugate gradient for efficient convergence near the local minimum. This relaxes strain and resolves clashes without driving the structure away from the experimental geometry.

Step 6: Structural Quality Control with MolProbity

Before committing to the computationally expensive docking stage, the pipeline evaluates the quality of the prepared receptor using MolProbity,^[5] the standard tool in structural biology for assessing model quality (used by the wwPDB validation pipeline for all deposited structures). The pipeline calls a MolProbity server with four sequential analyses:

• Oneline analysis — Computes an overall MolProbity score expressed as a percentile relative to structures of comparable resolution in the PDB.
• Clashscore — Counts serious steric overlaps (> 0.4 Å) per 1000 atoms. A high clashscore in the binding site is particularly damaging for selectivity analysis, as it can create artificial voids or occlusions that skew scores relative to the other protein.
• Ramachandran analysis — Evaluates backbone φ/ψ dihedral angles against distributions derived from ultra-high-resolution crystal structures.
• Rotamer quality — Assesses whether side-chain conformations match preferred rotamer distributions. Rotamer outliers in or near the binding site directly affect the shape and electrostatic character of the pocket.

When multiple conformations are available (e.g., from NMR ensembles), the pipeline automatically selects the conformation with the best (lowest) MolProbity score as the representative structure for docking.

Step 7: Binding Pocket Prediction with P2Rank

With a validated, high-quality receptor structure, the pipeline identifies druggable binding sites using P2Rank (version 2.5).^[6] P2Rank is a machine learning-based pocket predictor that uses a random forest model trained on the structural and physicochemical properties of protein surfaces. Unlike purely geometric methods that identify concave surface regions, P2Rank learns which surface features are associated with observed ligand binding in the PDB training set.

P2Rank operates on the protein surface as sampled by solvent-accessible surface (SAS) points. For each point, it computes a feature vector encoding local geometry (curvature, burial depth), physicochemical properties (hydrophobicity, charge, hydrogen-bond donor/acceptor capacity), and evolutionary conservation. The random forest outputs a per-point probability of belonging to a binding pocket; nearby high-probability points are clustered into pocket predictions, each ranked by a combined score.

The pipeline then: (1) extracts center coordinates of each predicted pocket; (2) calculates a bounding box using MDAnalysis, sized to enclose all SAS points plus a 2 Šmargin; (3) enforces a maximum box volume of 27,000 ų (the limit imposed by UniDock), scaling oversized boxes down while preserving their center and aspect ratio; and (4) checks new predictions against previously stored sites (within a 5 Šcenter tolerance) to avoid duplicate pockets.

In counterscreening, P2Rank is run with a larger max_number_pockets setting (up to 20 pockets per receptor by default) compared to standard single-target docking. This is intentional: when looking for a pocket to avoid, rare or cryptic pockets that might not appear in a low-budget search need to be surfaced so the user can make an informed selection.

Step 8: Binding Site Selection

After pocket detection, the agent presents all identified pockets to the user through a 3D visualization overlaid on the protein structure and a summary table of scores and positions.

In Mode 1 (two pockets on the same protein): The user selects two distinct pockets — the desired binding site and the off-target site to avoid — from the same receptor structure. The agent enforces that the two selected pocket IDs are different, preventing the trivial and uninformative case where the same pocket is designated as both target and counter-target.

In Mode 2 (two different proteins): The user selects one pocket from the desired protein's pocket list. This completes Phase 1. Phase 2 then begins independently for the off-target protein.

Step 9: Virtual Screening with UniDock Against Both Sites

With both binding sites prepared and identified, the compound library is docked against each site using UniDock,^[8] a GPU-accelerated molecular docking engine derived from AutoDock Vina.^[9] UniDock reimplements Vina's Monte Carlo–based conformational search on GPU hardware, enabling batch processing of thousands of ligands simultaneously.

Before docking, the compound library is filtered at the database level:

• Drug-likeness (QED ≥ 0.2): The Quantitative Estimate of Drug-likeness^[7] combines multiple Lipinski-style descriptors into a single score from 0 to 1. A threshold of 0.2 removes obvious non-drug-like compounds while retaining fragments and lead-like molecules.
• Nanoparticles: Carbon-only structures (fullerenes, nanotubes, etc.) are excluded.
• Metabolites: Endogenous metabolites are excluded to focus the screen on exogenous drug-like molecules.

A two-stage docking strategy balances throughput against accuracy:

Stage 1 — Fast screening across the full library: UniDock runs in fast search mode. Each ligand's best-scoring pose is evaluated by Vina's hybrid scoring function. Ligands scoring worse than −3.0 kcal/mol against the desired site are discarded.

Stage 2 — Detail re-docking of top candidates: The top candidates from Stage 1 are re-docked in detail search mode, which allocates significantly more Monte Carlo steps and evaluations per ligand. This stage uses PDBQT files with explicit hydrogen atoms rather than the implicit-hydrogen representation used in Stage 1, improving scoring accuracy for hydrogen-bond interactions.

Critically, in counterscreening, each compound is docked against both the desired site and the off-target site. The scoring against the off-target uses the same two-stage protocol and the same compound representations, ensuring that the scores from the two sites are directly comparable on the same energy scale.

Step 10: Selectivity Ranking

After docking is complete for both sites, compounds are ranked by a selectivity score, defined as the difference between the Vina score at the off-target site and the Vina score at the desired site:

Because Vina scores are negative (more negative = stronger predicted binding), a large positive selectivity score means the compound binds the desired site much more strongly than the off-target site. Compounds are sorted by this score in descending order, with the most selective compounds appearing first.

This ranking directly reflects the drug discovery objective: the ideal hit not only binds the desired target with high predicted affinity (a strongly negative desired-site score), but also fails to engage the off-target to any significant degree (a near-zero or weakly negative off-target score). The selectivity score encapsulates both criteria in a single quantity.

The top candidates are stored in the database with both docking scores, the selectivity score, SDF coordinates of both docked poses, and cross-references to the receptor structures and binding sites.

Step 11: Docking Quality Control with PoseBusters

High selectivity scores are a necessary but not sufficient condition for a compound to be a viable hit. The final step runs PoseBusters^[10] on the docked poses for all top-ranked compounds. PoseBusters applies a battery of deterministic, physics-based checks that include:

• Molecular sanitization: Correct valence, no charged atoms, no radicals.
• Ligand bond geometry: Bond lengths, bond angles, and torsion angles within expected chemical ranges.
• Internal clashes: No steric overlaps between atoms within the ligand.
• Receptor clashes: No steric overlaps between the docked ligand and the protein heavy atoms.
• Binding site proximity: The ligand center of mass is within the declared binding site, not docked to a remote surface patch.
• Flat ring systems: Aromatic rings are planar, as required by their electronic structure.

In a selectivity context, PoseBusters checks are applied to the poses from both binding sites. A compound that passes PoseBusters for its desired-site pose but fails for the off-target pose suggests that the off-target score reflects a spurious or unphysical interaction — potentially inflating the apparent selectivity. Conversely, a compound that passes both checks is a physically plausible binder at both sites, and its selectivity score reflects a genuine geometric and energetic preference for the desired site over the off-target.

The final results table displays the top-ranked selective compounds with their target scores, off-target scores, selectivity scores, and PoseBusters validation status for each pose.