Computational Structural Biology · Molecular Docking
Molecular Docking: Interactions, Tools & Analysis
Protein–ligand · Protein–peptide · Covalent warhead docking · Non-covalent interactions · Complete scoring and validation workflow
AutoDock Vina
PyRx
HADDOCK
HDock
ZDOCK
ClusPro
SwissDock
Covalent · H-bond · π-π · Ionic
System
Prep
Prep
Protein Preparation
Step 1
PDB · protonation
Ligand / Peptide Prep
Step 2
.mol2 · .pdbqt
Define Binding Site
Step 3
Grid box · active site
Docking Setup
Step 4
.mdp · config file
Dock ·
Analyse
Analyse
Non-Covalent Interactions
Step 5
H-bond · π-π · ionic
Covalent Docking
Step 6
Warhead · Cys-S-C
Docking Software
Step 7
PyRx · HADDOCK · HDock
Scoring & Validation
Step 8
RMSD · ΔG · PyMOL
1
System Preparation
▼
Receptor Preparation Workflow
- Download crystal structure from RCSB (e.g., PDB: 6LU7 for SARS-CoV-2 Mpro)
- Remove all water molecules, co-crystallized ligands, and crystallization artifacts
- Add polar hydrogen atoms only (non-polar H are merged into united-atom representation)
- Assign Gasteiger partial charges using AutoDock Tools or Open Babel
- Convert to
.pdbqtformat — required for AutoDock/Vina - Check for missing residues, non-standard amino acids, and disulfide bonds
python prepare_receptor4.py -r protein.pdb -o receptor.pdbqt -A hydrogens
# Or using Open Babel:
obabel protein.pdb -O receptor.pdbqt -p 7.4 # add H at pH 7.4
# Or using Open Babel:
obabel protein.pdb -O receptor.pdbqt -p 7.4 # add H at pH 7.4
Input format
.pdb / .cif
Output format
.pdbqt
Charges
Gasteiger / AM1-BCC
pH for H
7.4 (physiological)
Critical check: Inspect the active site before docking. Missing side chains near the binding pocket, wrong protonation of His/Asp/Glu residues, or incorrect disulfide bond states can completely invalidate docking results.
▼
Small Molecule Ligand
obabel ligand.smi -O ligand.pdbqt --gen3d --best -p 7.4
# Or using RDKit to generate 3D conformer:
python -c "from rdkit.Chem import AllChem; m=AllChem.EmbedMolecule(mol); AllChem.MMFFOptimizeMolecule(m)"
# Or using RDKit to generate 3D conformer:
python -c "from rdkit.Chem import AllChem; m=AllChem.EmbedMolecule(mol); AllChem.MMFFOptimizeMolecule(m)"
Peptide Preparation (for Protein-Peptide Docking)
- Build 3D peptide structure from sequence using PyMOL, Avogadro, or Modeller
- For short peptides (<10 aa): use extended or α-helix starting conformation
- For longer peptides: multiple conformers recommended — peptide can fold/rearrange upon binding
- For HADDOCK: provide FASTA sequence and define flexible residue ranges
Key Difference: Ligand vs Peptide
| Property | Small Molecule Ligand | Peptide |
|---|---|---|
| Size | Typically <500 Da | 500–5000+ Da |
| Flexibility | Defined rotatable bonds | Full backbone flexibility |
| Binding site | Deep pocket / active site | Flat interface / groove |
| Contact area | Small, precise | Large, distributed |
| Prediction difficulty | Moderate | High |
▼
AutoDock Vina Grid Box Configuration
# config.txt for AutoDock Vina:
receptor = receptor.pdbqt
ligand = ligand.pdbqt
center_x = 15.0 # box center coordinates (Å)
center_y = 12.5
center_z = 10.0
size_x = 25 # box dimensions (Å)
size_y = 25
size_z = 25
exhaustiveness = 16
num_modes = 10
energy_range = 3
receptor = receptor.pdbqt
ligand = ligand.pdbqt
center_x = 15.0 # box center coordinates (Å)
center_y = 12.5
center_z = 10.0
size_x = 25 # box dimensions (Å)
size_y = 25
size_z = 25
exhaustiveness = 16
num_modes = 10
energy_range = 3
Binding Site Identification Methods
- Known active site — center box on co-crystal ligand or catalytic residues
- Blind docking — cover the entire protein surface; computationally expensive but unbiased
- fpocket / SiteMap — computational pocket prediction from surface geometry
- Conservation analysis — highly conserved residues often mark functionally important binding sites
Box must be large enough to allow the full ligand to translate and rotate freely. For typical drug-like molecules a 20–25 Å cube is sufficient. For peptides, a 30–40 Å box or larger may be needed.
▼
AutoDock Vina Command
vina --config config.txt --out output.pdbqt --log log.txt
# Read results from log.txt:
Mode | Affinity (kcal/mol) | RMSD l.b. | RMSD u.b.
1 | -8.5 | 0.000 | 0.000
2 | -8.1 | 1.823 | 2.456
3 | -7.9 | 2.041 | 3.110
# Read results from log.txt:
Mode | Affinity (kcal/mol) | RMSD l.b. | RMSD u.b.
1 | -8.5 | 0.000 | 0.000
2 | -8.1 | 1.823 | 2.456
3 | -7.9 | 2.041 | 3.110
PyRx Workflow (GUI-based)
- Import receptor and ligand files via File → Open
- Convert ligand to AutoDock format using Open Babel in PyRx
- Select binding site using PyMOL visualizer embedded in PyRx
- Run AutoDock Vina wizard → set exhaustiveness → Run
- Export results CSV for downstream analysis
Exhaustiveness
8 (fast) – 32 (thorough)
Output poses
Top 9 binding modes
Score unit
kcal/mol (negative = better)
RMSD cutoff
<2.0 Å for redocking
2
Molecular Interactions
▼
Protein–Ligand Interaction Types
Small molecule inhibitors bind in well-defined pockets through a combination of the following interaction types:
Hydrogen Bond
2.5 – 3.5 Å
Directional
Glu123:O — Lig H-N
Hydrophobic
3.5 – 5.0 Å
Non-directional
Phe45:C — Lig CH₃
π–π Stacking
3.8 – 4.5 Å
Aromatic
Tyr88:Ring — Lig Ring
Ionic / Salt Bridge
2.0 – 4.0 Å
Electrostatic
Lys156:N⁺ — Lig O⁻
Protein–Peptide Interaction Types
Peptides bind through larger contact surfaces with greater flexibility:
H-bond (peptide)
2.6 – 3.4 Å
Backbone
Ser20:OH — Pep Ala3:O
β-Sheet Network
2.8 – 3.2 Å
Multi-H-bond
Backbone N/O pattern
Hydrophobic
3.6 – 4.8 Å
Extended surface
Leu67:C — Pep Ile2:C
Charge Interaction
2.2 – 3.8 Å
Ionic
Asp30:O⁻ — Pep Arg5:N⁺
Interaction Distance Comparison
Protein–Ligand vs Protein–Peptide (Å range)
Shorter distances = stronger, more directional interactions. Hydrophobic contacts span the widest range
Contact Area Comparison
Buried surface area (Ų) — ligand vs peptide
Peptides bury significantly larger surface areas — making them harder to displace but also harder to predict
Interaction Type Frequency in Drug-Protein Complexes
% occurrence across 500 PDB co-crystal structures (illustrative)
Hydrogen bonds and hydrophobic contacts dominate most protein–ligand binding interfaces
▼
What Is Covalent Docking?
In covalent docking, a ligand with a specifically designed reactive electrophilic group (warhead) is computationally positioned near a target nucleophilic amino acid — most commonly cysteine (Cys) — such that a chemical reaction occurs, forming a permanent covalent bond. This is irreversible binding.
Mechanism: Cys-S-Ligand Thioether Formation
Non-covalent Pre-binding
Warhead positions near Cys
Reversible stage
H-bond network guides orientation
Nucleophilic Attack
Cys210:S attacks warhead C
Reaction transition state
Michael addition / acrylamide
Covalent Bond
S–C bond ~1.8 Å
Irreversible
Thioether · permanent link
Common Warhead Types
| Warhead | Reaction type | Target residue | Bond formed |
|---|---|---|---|
| Acrylamide / vinyl sulfone | Michael addition | Cys (thiol) | Thioether C–S |
| Chloroacetamide | Alkylation (SN2) | Cys (thiol) | Thioether C–S |
| Epoxide | Ring opening | Cys, Ser, Lys | C–S / C–O / C–N |
| Aldehyde | Imine (Schiff base) | Lys (ε-NH₂) | C=N (reversible) |
| Activated ester | Acylation | Ser (hydroxyl) | Ester C–O |
| Boronic acid | Tetrahedral adduct | Ser (protease) | B–O (reversible) |
Software for Covalent Docking
- Schrödinger Suite (CovDock) — gold standard for covalent docking with reaction enumeration
- GOLD (CCDC) — supports covalent constraint during search
- AutoDock / Vina modifications — constrained bond can be implemented manually
- DOCKovalent (Shoichet lab) — large-scale covalent virtual screening
Covalent Bond Formation — Reaction Coordinate
Energy profile: non-covalent binding → transition state → covalent adduct
Non-covalent pre-binding (ΔG₁) positions warhead for nucleophilic attack; covalent bond (ΔG₂) provides irreversible stabilization
Covalent vs Non-Covalent — Binding Strength
Typical binding energy ranges (kcal/mol)
Covalent bonds provide 40–60 kcal/mol additional stabilization compared to non-covalent interactions
Important: Before the irreversible covalent bond can form, the ligand MUST first dock non-covalently and interact with surrounding residues to correctly position the reactive warhead. All regular non-covalent interactions (H-bonds, hydrophobic, ionic) are essential prerequisites for successful covalent bond formation.
3
Docking Software
▼
PyRx / AutoDock Vina
Protein–Ligand · GUI
GUI-based virtual screening platform that wraps AutoDock Vina and AutoDock4. Automates ligand preparation and batch docking. Widely used for small molecule drug discovery. Excellent for beginners and high-throughput screening.
HADDOCK
Protein–Protein · Data-driven
High Ambiguity Driven protein-protein DOCKing. Uses experimental data (NMR, mutagenesis) as ambiguous interaction restraints. Incorporates significant protein flexibility during CNS refinement. Excellent for complexes, multi-body docking, and various biomolecules.
HDock
Protein–Protein · Peptide
Template-free docking for protein-protein and protein-peptide interactions. Hierarchical approach: efficient rigid-body matching → energy minimization refinement. Focuses on shape complementarity and chemical interaction complementarity.
ZDOCK
Protein–Protein · FFT
Fast Fourier Transform (FFT)-based protein-protein docking. Exhaustive rigid-body search in translational and rotational space based on shape complementarity and electrostatics. Typically requires RDOCK or ZRANK post-processing for refinement and re-scoring.
ClusPro
Protein–Protein · Clustering
Server that uses PIPER (FFT-based with improved potential) followed by clustering of low-energy structures. Clustering identifies the most populated — and likely most biologically relevant — binding modes. Consistently top-performing in CAPRI protein docking benchmarks.
SwissDock
Protein–Ligand · Web
Free web-based service for flexible small molecule docking using EADock DSS algorithm. Makes molecular docking accessible without software installation. Supports flexible ligand docking into protein active sites with CHARMM energy function for scoring.
Docking Software Comparison — Best Use Cases
Suitability score (0–10) for different docking scenarios
No single tool is optimal for all scenarios — match software to your specific docking problem type
4
Scoring, Analysis & Validation
▼
Validation by Redocking
The most reliable validation is redocking the native co-crystal ligand back into the protein — the RMSD between the predicted pose and the experimental crystal pose should be <2.0 Å.
vina --config config.txt --out redock.pdbqt --log redock.log
rmsd = calculate_rmsd(native_ligand.pdb, redock_pose1.pdb)
# RMSD < 2.0 Å = successful validation
rmsd = calculate_rmsd(native_ligand.pdb, redock_pose1.pdb)
# RMSD < 2.0 Å = successful validation
Interaction Analysis Tools
- PLIP (Protein-Ligand Interaction Profiler) — automated detection of all interaction types, web server + Python library
- LigPlot+ — 2D schematic diagram of protein-ligand interactions
- PyMOL — 3D visualization;
find_polar_contacts within 3.5, ligand, receptor - Discovery Studio Visualizer — detailed 2D/3D interaction maps for publications
Binding Affinity Distribution
Top 9 poses from AutoDock Vina (kcal/mol)
Lower (more negative) binding energy = stronger predicted binding. Pose 1 is always the top-ranked
RMSD of Docking Poses
Lower bound and upper bound RMSD vs Pose 1
Poses with RMSD <2 Å from Pose 1 represent the same binding mode (cluster)
H-bond Distance Distribution
From docking poses — donor–acceptor distances
Sharp peak at 2.8–3.0 Å confirms strong hydrogen bond geometry
Binding Energy vs Number of H-bonds
Correlation across top-ranked poses
More H-bonds generally correlates with better binding affinity — but other interactions also contribute
Interaction Profile per Residue
Frequency of each interaction type at key active-site residues
Identifies critical hotspot residues — essential for lead optimization and mutation studies
Validation Criteria Summary
| Validation check | Acceptable threshold | Tool |
|---|---|---|
| Redocking RMSD | <2.0 Å | AutoDock Vina + RMSD calc |
| Binding energy (typical drug) | −6 to −12 kcal/mol | Vina score |
| Pose cluster population | Top mode should dominate | Vina / ClusPro |
| Key interaction preserved | Same as crystal structure | PLIP / LigPlot+ |
| Ligand strain energy | <5 kcal/mol | Schrödinger / MMFF |
| Protein-peptide interface RMSD | <3.0 Å (backbone) | PyMOL align |
⚙
Key Docking Parameters
Exhaustiveness
8–32
Grid box (ligand)
20–25 Å cube
Grid box (peptide)
30–40 Å cube
H-bond cutoff
2.5–3.5 Å
π-π stacking
3.8–4.5 Å
Hydrophobic cutoff
3.5–5.0 Å
Redocking RMSD
<2.0 Å
Covalent S–C bond
~1.8 Å
📈
Analysis Metrics
ΔG bind
Predicted binding free energy (kcal/mol)
RMSD
Pose deviation from reference crystal
Ki / IC₅₀
Estimated inhibition constant from ΔG
H-bonds
Count and geometry of H-bond contacts
BSA
Buried surface area at interface (Ų)
PLIP score
Interaction profiling per residue
Pose cluster
Population of top binding mode
LE
Ligand efficiency = ΔG / heavy atoms