Module 6: Molecular Docking

Understand how computer algorithms predict ligand binding configurations inside a receptor pocket. Explore conformational searches, scoring functions, and validation strategies.


Learning outcomes

  • Compare systematic, stochastic, and incremental conformational search strategies.
  • Distinguish force-field, empirical, and knowledge-based scoring functions.
  • Prepare a receptor and ligand correctly, including protonation and charge assignment.
  • Validate a docking protocol by redocking and interpreting RMSD.
  • Explain why covalent and PROTAC ternary docking need different machinery.

1. What is Molecular Docking?

Molecular docking is a computational simulation technique that models the interaction between a small molecule (ligand) and a macromolecular target (typically a protein).

The docking process has two primary components:

  1. Search Algorithm: Explores the conformational space of the ligand inside the active site.
  2. Scoring Function: Evaluates the binding energy of each generated pose to identify the most stable orientation.

2. Conformational Search Methods

Systematic Search

Samples a finite torsion grid (e.g. 60° increments). Exhausting that grid does not guarantee the continuous global minimum, the best receptor conformation, or the experimental bound pose; the grid also grows combinatorially.

Stochastic Search

Performs random conformational moves. Prominent algorithms include Monte Carlo and Genetic Algorithms, where coordinates mutate and cross over.

Incremental Construction

Splits the ligand into fragments. The anchor fragment is docked first, and the remaining segments are attached step-by-step.

Interactive Playground: Pose Evidence Challenge

Choose which synthetic pose you would advance. Make the decision from the interaction pattern and geometry—not the score alone.

Guided experiment

Which pose deserves priority for expert inspection and experimental follow-up?

Explore further: move a pose yourself

Use the coordinate sliders or numeric inputs to move the ligand. The dashed outline is the fixed reference pose; the colored ligand is your current pose.

-+
Synthetic teaching scoreApproaching pocket...

-8.28score units

Teaching boundary: the pose geometry, interaction labels, and score in this playground are synthetic. The score is not a calibrated binding free energy, and a favorable value does not establish binding or identify the experimental pose.

3. Target Coordinate Selection: Navigating the PDB

Before running a docking simulation, researchers must select a protein target coordinate set. The Worldwide Protein Data Bank (PDB) houses hundreds of thousands of structures, but they are not created equal. Poor target coordinates will invalidate all downstream scoring.

Checklist for Systematic PDB Selection:

  1. Sequence Fidelity via UniProt: Search by UniProt accession ID, not common protein names, to filter out homologous species variants and ensure exact sequence identity.
  2. Experimental Method & Resolution: X-ray crystal structures are the default choice for docking as they represent frozen energy snapshots with defined coordinate densities. NMR provides structural ensembles that require complex ensemble docking. Resolution should ideally be below 2.5 Å.
  3. Refinement Metrics:Inspect the structure's R-free (relative to R-work) and Ramachandran outlier percentages.
  4. Engineered Mutations: Crystallographers frequently introduce mutations to prompt crystallization. Ensure the active site loop residues are wild-type.
  5. Loop Completeness: Weak electron densities in flexible loops lead crystallographers to omit loops entirely. If gaps occur in the active site, loops must be modeled in using loop refinement tools before docking.
  6. Biological Assembly vs. Asymmetric Unit: Download the biological assembly (the functional form under physiological conditions), not the asymmetric unit (a mathematical crystallographic convenience). Extracting a monomer from a functional dimer exposes hydrophobic interfaces to solvent, causing artificial coordinate collapses.
The mmCIF Format Transition:

The legacy .pdb file format (designed in the 1970s for punch cards) suffers from a hard limit of 99,999 atoms and 62 chains. In the era of cryo-EM macromolecular assemblies, the format breaks down. The wwPDB's standard format is PDBx/mmCIF (.cif). From 21 July 2027, new depositions will receive only 12-character extended IDs (e.g. pdb_1000axyz; existing entries keep their ID with a pdb_0000 prefix, so 2v5z becomes pdb_00002v5z), and legacy .pdb files will not be generated for those new entries. Modern pipelines should migrate to .cif inputs exclusively.

When Experimental Structures Are Unavailable

Docking may begin from a predicted or comparative receptor model, but generating and validating that model is a structural-bioinformatics task. This module therefore focuses on the narrower question: does the resulting pocket support the proposed docking experiment?

Docking-specific acceptance gate

  • Inspect confidence and completeness locally at the pocket, not only across the global fold.
  • Challenge protonation, side-chain rotamers, cofactors, metals, waters, and alternative pocket states.
  • Use an ensemble when plausible receptor conformations materially change the site.
  • Where reference ligands exist, test pose recovery and screening enrichment before prospective use.

A well-formed model is not automatically a useful docking receptor. Task-specific retrospective performance is the final acceptance test.

4. Ligand Coordinate Preparation: Electrostatic Charge Selection

Force-field and AutoDock 4 scores use ligand partial charges in a Coulomb term, so the charge model changes the number those programs report. AutoDock Vina does not. Its score has no electrostatic term, and the charge column in a PDBQT file is kept for format compatibility. Atom type still matters in Vina, because hydrogen-bond donors and acceptors are typed atoms.

While full ab initio quantum mechanics methods (like RESP) are highly accurate, they are computationally prohibitive for virtual screening libraries. Semi-empirical charge models approximate this quality at a fraction of the cost:

Charge ModelElectrostatic QualityCalculation SpeedRecommended Use Case
AM1-BCCVery High (errors < 0.1 e)FastCharge-using scores and explicit solvent, where a Coulomb term is evaluated. A Vina score does not use it.
PM6High (improves H-bonding)MediumChallenging systems requiring precise hydrogen-bonding parameters.
RESPExcellent (Gold Standard)Very SlowLead optimization and final pose verification. Too slow for screening.
GasteigerPoorExtremely FastQuick initial pre-filtering or triage of massive compound libraries.

The PDBQT Format: Charges and Atom Types

A PDB file stores element symbols and Cartesian coordinates. AutoDock 4 and Vina read PDBQT, which adds a charge column, an AutoDock atom type, and a torsion tree. AutoDock 4 uses the charge. Vina uses the atom types and the torsion tree.

Partial Charge (Q)

Stores a partial charge, often Gasteiger or AM1-BCC. AutoDock 4 uses it in the Coulomb term. Vina does not read it into the score.

AutoDock Atom Type (T)

Maps atoms to specific docking categories (e.g., HD for donor hydrogen, A for aromatic carbon, OA for acceptor oxygen) to look up force field parameters.

Torsional Topology

Defines the molecule as a tree with ROOT and BRANCH statements, specifying which bonds are frozen and which can rotate.

5. Types of Scoring Functions

A

Force-Field-Based Functions

Calculate binding energy from electrostatics (Coulomb) and van der Waals interactions, using mechanical parameters. Examples: DOCK, GOLDScore.

B

Empirical Scoring Functions

Sum of weighted terms (hydrogen bonds, contacts, rotatable bonds) parameterized using experimental binding affinities. Examples: ChemScore, GlideScore.

C

Knowledge-Based Functions

Derived from statistical distributions of atom-atom contact pairs in structural databases. Favorable shapes yield lower potentials. Examples: PMF, DrugScore.

AutoDock 4 and AutoDock Vina use different scores

AutoDock 4 (Huey et al.) estimates affinity as a weighted sum of physical terms, including a Coulomb term and a desolvation term:

ΔG_AD4 = w_vdW × Σ vdW + w_hbond × Σ hbond + w_elec × Σ elec + w_desolv × Σ desolv + w_rot × N_rot
  • vdW: a Lennard-Jones contact between atoms i and j.
  • hbond: a directional donor–acceptor term.
  • elec: Coulomb interaction of the partial charges, with a distance-dependent dielectric.
  • desolv: a penalty for burying polar groups, estimated from the volume of atoms around each pair.
  • N_rot: a penalty proportional to rotatable bonds frozen on binding.

Vina (Trott and Olson, 2010) is a different function. It sums two attractive steric Gaussians, a repulsion, a hydrophobic contact, a hydrogen-bond term, and the same kind of torsion count:

score_Vina = gauss₁ + gauss₂ + repulsion + hydrophobic + hbond + w_rot × N_rot

There is no Coulomb term and no desolvation term. Polar recognition enters through the typed hydrogen-bond term. A strong Vina score is therefore not evidence that the electrostatics were evaluated.

Docking Software Comparison

Different docking programs employ distinct algorithms and scoring paradigms. The choice of software depends on the target system, licensing constraints, and screening scale.

SoftwareAlgorithmScoringLicenseBest For
AutoDock VinaIterated Local SearchEmpirical hybridOpen-sourceAcademic VS, large screens
GlideHierarchical HTVS/SP/XPGlideScoreCommercialPharma lead optimization
GOLDGenetic AlgorithmChemScore, GoldScoreCommercialMetalloenzymes, flexible
FlexXIncremental ConstructionBohm empiricalCommercialFragment-based
rDockCavity-based + GAEmpirical + vdWOpen-sourceRNA-ligand docking

Practical Docking Protocol

A systematic end-to-end docking workflow. Each step records assumptions and controls to support reproducible pose hypotheses, not guaranteed binding predictions.

1

Target Preparation

Download the PDB structure, remove crystallographic waters (except conserved ones mediating key interactions), assign protonation states using PropKa or H++ at physiological pH, and add missing hydrogens.

2

Binding Site Definition

Define a grid box centered on the co-crystallized ligand with 0.375 Å spacing and a 10–15 Å buffer around the binding pocket to allow full ligand exploration.

3

Ligand Preparation

Generate 3D coordinates from SMILES, enumerate relevant tautomers and protonation states at pH 7.4, and generate diverse low-energy conformers. Assign charges that match the scoring function: Meeko writes Gasteiger charges into the PDBQT, AutoDock 4 uses them in its Coulomb term, and the Vina score ignores them.

4

Docking Run

Set the exhaustiveness parameter to control thoroughness: 8 for fast screening, 32 for thorough sampling, and 128 for exhaustive poses. Generate 9–20 poses per ligand for adequate conformational coverage.

5

Post-Processing

Rank poses by docking score, filter candidates by interaction fingerprint similarity, visually inspect the top 20 hits for chemically sensible binding modes, and apply consensus scoring across multiple scoring functions to improve reliability.

6. Covalent Docking: Modeling Chemical Warheads

Traditional small-molecule drugs bind via reversible, non-covalent interactions. However, a major class of modern therapeutics consists of covalent inhibitors(e.g. Sotorasib targeting KRAS G12C, or Ibrutinib targeting Bruton's Tyrosine Kinase). These molecules contain a mildly electrophilic "warhead" that forms a covalent bond (irreversible or reversible, depending on chemistry) with a nucleophilic amino acid residue in the binding pocket.

Because standard docking scoring functions are parameterized only for non-covalent forces (electrostatics, van der Waals), simulating covalent binding requires specialized covalent docking protocols:

1. Non-Covalent Fit

The ligand first docks non-covalently. The scoring function evaluates whether the scaffolds fit the pocket topography and positions the reactive warhead (e.g. acrylamide, sulfonyl fluoride) within reaction distance (< 3.5 Å) of the target residue.

2. Chemical Linkage

The algorithm forms a virtual covalent bond between the nucleophile (typically Cys, Ser, or Lys) and the electrophilic carbon. The system conformational search resolves the dihedral angles of the newly formed linkage.

3. Complex Minimization

The bound ligand-protein complex is energy-minimized. Scoring functions evaluate the local strain of the covalent linker alongside the remaining non-covalent interactions of the scaffold to rank candidate covalent hits.

7. Ternary Complex Docking: PROTACs & Molecular Glues

Covalent docking bent the rules by adding a bond. Targeted Protein Degradation (TPD) breaks a deeper assumption: that docking means fitting one ligand into one pocket. Here you must dock a three-body complex.

As Module 2 showed, most disease-driving proteins have no druggable pocket at all. TPD changes the mechanism, but not every target becomes tractable: rather than inhibiting the target, it recruits the cell's own Ubiquitin-Proteasome System to destroy it. A deep inhibitory pocket may not be needed, but target engagement, productive ligase geometry, cellular delivery, ubiquitination, and proteasome access still matter.

PROTACs vs. Molecular Glues

PROTACs (Proteolysis Targeting Chimeras)Bifunctional molecules: one head binds the target, a flexible linker spans the gap, and the other head recruits an E3 ubiquitin ligase (Cereblon, VHL). The ligase tags the target with ubiquitin, marking it for proteasomal destruction.
Molecular GluesMonovalent and compact. Rather than tethering two proteins like a leash, a glue sits at the interface, remodelling the E3 surface to create a composite pocket that captures the target.

Cooperativity (α) — why ternary docking is not just docking twice

The two binding events are not independent. Cooperativity relates the ternary dissociation constant (binding of the second protein to the binary complex) to the binary one (binding to the isolated protein):

α = K_D(binary) / K_D(ternary)

α > 1 is positive cooperativity: favourable protein-protein contacts induced by the PROTAC make the ternary complex more stable than either binary interaction predicts. α < 1 is negative cooperativity; unfavorable geometry or other energetic effects may contribute, not necessarily a direct steric clash. This is why a PROTAC built from two excellent binders can still fail completely: the complex, not the compound, is the drug.

Why the whole toolkit strains here

Ternary docking breaks assumptions from three different modules at once. The search problem explodes: you sample two rigid-body placements plus a floppy linker with many rotatable bonds (Module 4). The scoring problem becomes protein-protein, not protein-ligand. And PROTACs sit far outside Lipinski space — routinely 800–1,100 Da — so the Rule of Five you meet in Module 12 simply does not apply. TPD is the clearest demonstration in this course that the standard methods are approximations with edges.

Interactive Playground: PROTAC Ternary Assembly & the Hook Effect

Tune the linker, cooperativity, and concentration in an imposed bell-shaped teaching curve. It illustrates a possible hook effect, not a fitted equilibrium model or predicted degradation response. Real systems may not show a hook within the tested range.

PROTAC concentration100 nM
Linker length10 atoms
Cooperativity α2.0

α > 1 positive · α < 1 negative (surfaces clash)

Illustrative assembly score

99.4 / 100

Illustratively favorable region: The illustrative geometry and concentration terms are favorable. This toy model does not calculate binding, ubiquitination, or degradation.

2550751001e-101e-91e-81e-71e-61e-51e-4[PROTAC] (M, log scale)illustrative score / 100hook effect →
Why the curve turns over

At low dose there is too little PROTAC to bridge both proteins. At the optimum, each molecule links a target to an E3. But push the concentration high enough and every protein gets its own PROTAC molecule — target-PROTAC binaries and E3-PROTAC binaries — and none are bridged. The productive ternary complex is competed away by the drug itself. A hook can arise when binary saturation competes with ternary assembly. Its presence and location depend on affinities, abundance, kinetics, and assay context; not every degradation curve is bell-shaped.

8. Redocking Validation (RMSD)

Before trust is placed in docking scores, protocols must be validated. The standard validation method is self-docking (redocking):

Root Mean Square Deviation (RMSD): The measure of spatial distance between atoms of the predicted docking pose and the actual co-crystallized native structure in the crystal. A heavy-atom RMSD ≤ 2.0 Å is a common pose-recovery heuristic, not general protocol validation. Define atom correspondence, symmetry handling, receptor alignment, and reference coordinates; test alternate receptor states, cross-docking, and held-out activity ranking where relevant.

Knowledge check

Checkpoint

Research decision checkpoint

Helpful background: Module 3: Fundamentals of Ligand-Receptor Interactions · Module 4: Molecular Mechanics & Force Fields · Module 5: Cheminformatics & Molecular Representations · Module 13: Structural Bioinformatics & Molecular Visualization

Redocking gives 1.5 Å RMSD. Is the protocol ready for prospective ranking?

Your output: Write a pose-validation report including atom mapping, receptor state, settings, alternative poses, and remaining tests.

Compare your reasoning

This supports one pose-recovery test only. Check symmetry, preparation sensitivity, cross-docking, known actives/decoys, and experimental follow-up appropriate to the decision.

Evidence anchor: AutoDock Vina: accuracy and scoring assumptions. A first reading and this prompt do not replace a validated practical exercise.

Knowledge checkQuestion 1 of 3

When a score actually evaluates a Coulomb term, why is AM1-BCC preferred over Gasteiger charges?