How to Dock a Ligand to a GPCR or Membrane Protein with AutoDock Vina (Step by Step) - StemSkills Lab
Skip to content

How to Dock a Ligand to a GPCR or Membrane Protein with AutoDock Vina (Step by Step)

How to Dock a Ligand to a GPCR or Membrane Protein with AutoDock Vina (Step by Step)

Docking into a GPCR follows the same AutoDock Vina workflow as any protein, with four changes: choose a structure whose activation state matches your ligand class, delete the BRIL or T4 lysozyme fusion and the crystallisation lipids, centre the grid box on the co-crystallised ligand inside the seven-transmembrane bundle, and validate poses against known anchor contacts.

G protein-coupled receptors are the target class most MSc students end up picking for a thesis, and for a defensible reason: 475 approved drugs, roughly 34% of everything the FDA has cleared, act at 108 different GPCRs, according to Hauser and colleagues in Nature Reviews Drug Discovery. The problem is that almost every beginner docking tutorial, including the official one, uses a soluble globular enzyme. A kinase sits in water, its pocket opens onto the surface, and you can see the cleft in any viewer. A GPCR does none of that. Its pocket is buried inside a bundle of seven helices, the file you download contains a protein that does not exist in nature, and the ligand you are chasing may not even bind the state you downloaded.

This guide covers the parts that are actually different. It was written by the StemSkills Lab team, from 10+ years in sequence and structural bioinformatics, drug discovery and design, and multiscale molecular modeling. If you have not run a docking job at all yet, start with our molecular docking pillar guide and come back.

What actually changes when the target is a GPCR instead of a soluble protein?

Four things, and only four. The Vina binary, the scoring function and the PDBQT format are identical.

  1. The structure is a chimera. Most crystal structures of GPCRs were only obtainable by fusing a small soluble protein into an intracellular loop or the N-terminus. That fusion is in the coordinate file and must go.
  2. The construct is mutated. Thermostabilising point mutations were often introduced to get a diffracting crystal. Some of them sit close enough to the pocket to matter.
  3. The pocket is buried and the file is full of lipid. Cholesterol, monoolein, detergents and polyethylene glycol fragments are packed against the outside of the helices, and a blind search box will happily dock your ligand into one of the greasy surface grooves they leave behind.
  4. The receptor has states. An antagonist docked into an agonist-bound active structure, or the reverse, produces a number that means nothing.

Which GPCR structure should you dock into, active or inactive?

Match the state to the pharmacology of the compound you are testing. An inactive (antagonist or inverse agonist bound) structure has a slightly contracted orthosteric pocket and a closed intracellular face. An active structure, usually captured with an agonist plus a G protein or a nanobody, has a rearranged transmembrane helix 6 and a subtly different pocket shape. Docking an agonist series into an inactive structure systematically under-ranks the agonists.

Structure typeWorked exampleWhat is fused or boundDock this when
Inactive, antagonist bound2RH1, human beta-2 adrenergic receptor, 2.40 AT4 lysozyme chimera, carazolol (ligand code CAU), 3 cholesterols, palmitic acid, maltoseYour compounds are antagonists, inverse agonists or unknown blockers
Inactive, antagonist bound3EML, human adenosine A2A receptor, 2.60 AT4 lysozyme chimera, ZM241385You want a well-characterised, heavily benchmarked test case
Active, agonist and transducer bound3SN6, beta-2 adrenergic receptor and Gs complex, 3.20 AGs heterotrimer (3 chains) plus a camelid nanobody, agonist BI-167107Your compounds are agonists or partial agonists

Two practical notes on that table. First, resolution is not a tie-breaker you can ignore: 3SN6 is a 3.20 A structure, so its side chain rotamers are less reliable than the 2.40 A of 2RH1, and you should be more forgiving about small pose deviations. Second, an active-state entry brings extra baggage. In 3SN6 you are not just deleting a fusion, you are deleting three G protein chains and a nanobody before the receptor is dockable.

To find the right entry in the first place, use GPCRdb rather than a plain PDB keyword search. It classifies structures by receptor, activation state and coupled transducer, which is precisely the filter you need and precisely the one RCSB does not give you.

What do you have to delete from a GPCR PDB file before docking?

Work through the raw file in this order. Every step is a deletion, and every deletion is checkable.

1. The fusion partner

Two fusions dominate the GPCR structures you will meet. T4 lysozyme, spliced into the third intracellular loop, is the one in 2RH1 and 3EML. BRIL, a thermostabilised apocytochrome b562RIL carrying the M7W, H102I and R106L mutations, was introduced as a smaller and more rigid alternative in the fusion partner survey published in Structure in 2012, and now appears in a large share of newer entries.

Whichever one you have, it is a separate stretch of residues in the same chain, usually with a residue numbering jump. Open the entry page, read the macromolecule description (2RH1 is annotated as a “beta-2-adrenergic receptor/T4-lysozyme chimera”), find the residue range of the fusion, and delete it. In PyMOL:

fetch 2rh1, async=0
remove solvent
remove not polymer and not resn CAU
# inspect the numbering break, then remove the lysozyme range
remove chain A and resi 1002-1161
save receptor_clean.pdb

Do not skip the visual check afterwards. A fusion left in place gives you a receptor with a large soluble domain hanging off the intracellular side, and if your box is even slightly misplaced, Vina will find a perfectly plausible pocket in it.

2. The thermostabilising mutations

Mutations introduced for crystallisation are recorded in the SEQADV records of the PDB file, which flag every difference between the deposited sequence and the reference database entry:

grep SEQADV 2rh1.pdb

Read the output against the UniProt sequence for your receptor. Most thermostabilising mutations sit in the transmembrane core or the intracellular face and are far from the orthosteric site, in which case leave them alone; reverting them by hand introduces a modelled side chain that is less trustworthy than the crystallographic one. Revert only if a mutated position is inside your binding site, and say so in your methods either way. This is a sentence your examiner will look for.

3. Lipids, detergents and crystallisation additives

List every non-water heteroatom before deciding what to keep:

grep '^HETATM' 2rh1.pdb | cut -c18-20 | sort | uniq -c | sort -rn

For 2RH1 that returns cholesterol, palmitic acid, maltose, a dodecaethylene glycol fragment, sulfate, 1,4-butanediol and acetamide alongside the carazolol. None of them belong in a docking receptor. They are annotated in the entry as crystallisation components, and cholesterol in particular is described there as necessary for crystallisation rather than as a functional part of the receptor. Delete all of them, and delete the waters too unless you have a specific, defensible reason to keep a named bridging water.

4. Everything that is not the receptor chain

For an active-state complex, this is the big one. In 3SN6, chains A, B and C are the Gs heterotrimer and chain E is the nanobody. Only chain D is the receptor. Keep chain D and discard the rest.

Want the guided, hands-on version?

Our live Molecular Modeling & MD Simulations cohort bootcamp takes you from zero to running real docking and MD workflows, with a portfolio project for your grad-school applications.

Join the waitlist (free) →

Should you keep the co-crystallised ligand?

Keep it in a separate file, never in the receptor. It has two jobs after that.

It defines your box centre, because its centroid is the single most reliable description of where the orthosteric site is. And it is your positive control: re-docking that ligand back into the cleaned receptor is the experiment that tells you whether your whole preparation is sound. If the redocked pose does not come back close to the crystallographic one, the field’s usual pass mark being a heavy-atom RMSD under 2.0 A, your preparation is wrong and no result on a novel compound is worth reporting.

How do you place the grid box on a buried orthosteric pocket?

You do not place it by eye. Compute the centroid of the extracted reference ligand from the fixed-width PDB coordinate columns:

awk '/^HETATM/ && substr($0,18,3)=="CAU" {x+=substr($0,31,8); y+=substr($0,39,8); \
     z+=substr($0,47,8); n++} END {printf "%.2f %.2f %.2f\\n", x/n, y/n, z/n}' 2rh1.pdb

Run against the deposited 2RH1 file, that returns -29.52 9.23 6.94 over carazolol’s 22 heavy atoms, and those three numbers become center_x, center_y and center_z. Substitute your own ligand code for CAU. PyMOL’s centerofmass command on the same selection gives you the equivalent value.

For the size, the official AutoDock Vina FAQ gives the rule in five words: “As small as possible, but not smaller.” It adds a hard ceiling, that you “should probably avoid search spaces bigger than 30 x 30 x 30 Angstrom, unless you also increase –exhaustiveness”. For a GPCR orthosteric site this advice is unusually easy to follow, because the pocket really is small and enclosed. A cube of 20 to 25 A on a side, centred on the reference ligand, covers the orthosteric site and the extracellular vestibule above it without spilling out past the helices into the lipid-facing surface.

The failure mode to watch for is the opposite of the soluble-protein one. Beginners docking a kinase make the box too small and clip the pocket. Beginners docking a GPCR make the box too big, it escapes the bundle, and top poses come back sitting on the outside of helices 3 and 4 where the cholesterol used to be. If your best pose is on the outer surface of the receptor, the box is the bug.

If you also want to know which part of the receptor is inside the membrane, submit the cleaned structure to the PPM server from the OPM database (Lomize et al., Nucleic Acids Research 2012). It returns the calculated position of the protein relative to the lipid bilayer, which makes it obvious at a glance whether your box is inside the bundle or out in the lipid.

What does the full command sequence look like?

These are the commands from the official AutoDock Vina basic docking tutorial, with your own filenames and the centre you just computed. Prepare the receptor and write the box definition in one step:

mk_prepare_receptor.py -i receptor_cleanH.pdb -o receptor -p -v \
  --box_size 22 22 22 --box_center -29.52 9.23 6.94

Prepare the ligand from a 3D SDF with correct protonation:

mk_prepare_ligand.py -i ligand.sdf -o ligand.pdbqt

Then run the docking. The default exhaustiveness is 8; the official tutorial raises it to 32 for a difficult ligand to get a more consistent result, and a flexible drug-like compound in a GPCR pocket is exactly that case:

vina --receptor receptor.pdbqt --ligand ligand.pdbqt \
     --config receptor.box.txt \
     --exhaustiveness=32 --out ligand_vina_out.pdbqt

Run the reference ligand through the identical pipeline first. Its redocking RMSD is the number that licenses everything after it.

How do you know the pose is right if you cannot trust the affinity number?

You check contacts, not kcal/mol. The Vina developers are explicit about this in the FAQ: “The predictive accuracy varies a lot depending on the target, so it makes sense to evaluate AutoDock Vina against your particular target first.” A single affinity value from a single receptor conformation is a ranking hint, not a measurement.

What a GPCR gives you instead is a set of published anchor contacts. Aminergic receptors are the clearest case. In the beta-2 adrenergic receptor, carazolol makes close contacts under 3 A with Asp113 (position 3.32 in Ballesteros-Weinstein numbering) and Asn312 (7.39). Asp113 is not a decorative contact: mutating it to asparagine abolishes detectable antagonist affinity and drops agonist potency by more than four orders of magnitude. If your docked amine-containing ligand does not engage Asp3.32, the pose is wrong regardless of what the score says.

So build a three-part check for every pose you keep:

  • Anchor contact. Does the pose reproduce the known salt bridge or hydrogen bond for that receptor family? Look it up in the paper for your reference structure.
  • Depth. Is the ligand inside the helical bundle, between the extracellular ends of helices 3, 5, 6 and 7, or is it perched on the lipid-facing outside?
  • Consistency. Do the top few modes converge on one binding mode, or does every mode look different? Scattered modes usually mean the box is too large or exhaustiveness is too low.

Troubleshooting: real failures and their fixes

  • Top pose sits on the outside of the helices. The box escaped the bundle. Recompute the centre from the reference ligand centroid and shrink the size to 22 A or less.
  • Top pose sits in a compact globular domain with sensible-looking contacts. You docked into the T4 lysozyme or BRIL fusion. Go back and delete it.
  • Redocking the reference ligand gives an RMSD well above 2.0 A. Usually protonation. Check that the receptor was prepared with hydrogens added and that the ligand SDF has the right charge state at physiological pH; a neutral amine will not form the Asp3.32 salt bridge.
  • Every run gives a different top pose. Raise exhaustiveness and shrink the box. Note that Vina is stochastic, so set a seed with --seed when you need a reproducible run for your thesis.
  • Suspiciously strong affinities across your whole compound set. Check for leftover heteroatoms. A cholesterol or monoolein you forgot to delete creates an artificial hydrophobic cavity that scores well and means nothing.
  • Agonists rank below antagonists in a series you know is agonist-driven. Wrong activation state. Redo the run against an active-state structure.

Where this fits in your training

Docking is one checkpoint on a longer path. A GPCR docking pose is a hypothesis, and the standard next move is a molecular dynamics simulation in an explicit lipid bilayer to see whether the pose survives. Our computational biology skills roadmap lays out the order to learn these in, and the molecular docking guide covers the fundamentals this article assumes.

Want the guided, hands-on version?

Our live Molecular Modeling & MD Simulations cohort bootcamp takes you from zero to running real docking and MD workflows, with a portfolio project for your grad-school applications.

Join the waitlist (free) →

Frequently asked questions

Do I need to include the lipid bilayer when docking to a GPCR?

No. AutoDock Vina docks a ligand into a rigid receptor and has no representation of a membrane. The bilayer matters for deciding where your box goes and for any molecular dynamics you run afterwards, not for the docking calculation itself. Delete the crystallographic lipids and detergents from the receptor file.

Can I dock into a homology model or an AlphaFold model of a GPCR?

You can, with two cautions. Predicted models usually come out in an inactive-like state, and side chain rotamers in the pocket are the least reliable part of the model. Always run the redocking control using a ligand from an experimental structure of a close homologue before you trust any result on a model.

What box size should I use for a GPCR orthosteric pocket?

Start at 22 x 22 x 22 A centred on the co-crystallised ligand centroid. That covers the orthosteric site and the vestibule above it. The Vina FAQ advises staying under 30 x 30 x 30 A unless you also raise exhaustiveness, and for a buried GPCR pocket you should be comfortably below that ceiling.

Should I remove the co-crystallised ligand before docking?

Yes, remove it from the receptor file, but save it separately. You need it to compute the box centre and to run the redocking control that validates your preparation.

Which AutoDock Vina version should I cite in my thesis?

Cite both the original method paper by Trott and Olson (Journal of Computational Chemistry, 2010) and, if you used version 1.2 or later, Eberhardt, Santos-Martins, Tillack and Forli, “AutoDock Vina 1.2.0: New Docking Methods, Expanded Force Field, and Python Bindings” (Journal of Chemical Information and Modeling, 2021). Record the exact version string your binary reports, along with your exhaustiveness and random seed.

Do thermostabilising mutations have to be reverted?

Only if a mutated residue lies inside your binding site. Most sit in the transmembrane core or the intracellular face and are irrelevant to orthosteric docking. Check the SEQADV records, state what you found in your methods, and justify whatever you decided.

Think you know Molecular Docking?
Take the free StemSkills assessment and earn a verifiable certificate you can download and add to your LinkedIn profile.
Start the free assessment

Keep going

AutoDock Vina Errors and How to Fix Them: The 10 That Stop Beginners Decode the 10 AutoDock Vina errors that stop beginners, from PDBQT parsing failures to grid box mistakes, with… How to Predict a Protein-Protein Complex with AlphaFold-Multimer in ColabFold (Step by Step) Learn how to model a protein-protein complex free in ColabFold, set the chain break correctly, and judge ipTM,… Discovery Studio Visualizer Tutorial: How to Analyse Docking Results and Make 2D Interaction Diagrams Load your docked pose, list every contact with its distance, and export the 2D interaction diagram your thesis…
See live workshops