How to Pick the Best Compound After Docking and MD
Rank docked compounds on evidence, not on the raw Vina score. Convert each score to ligand efficiency so chemotypes of different sizes compare fairly, drop anything failing drug-likeness or PAINS filters, then carry three candidates into MD and decide on four criteria: ligand RMSD stability, contact occupancy, MM-PBSA binding energy, and per-residue agreement.
Every tutorial ends at the same place. You have a results table, twenty poses, and a column of numbers in kcal/mol. Nobody tells you what to do next, and the thesis chapter needs one compound with a reason attached. This guide is the decision step: how the StemSkills Lab team, with 10+ years in structural bioinformatics, drug design and multiscale molecular modeling, would build a shortlist that survives a viva.
Why can’t you just pick the compound with the best docking score?
Because a docking score is a fast approximation, and it was never designed to rank different chemical series against each other. The AutoDock Vina documentation is direct about this in its own FAQ, noting that “the predictive accuracy varies a lot depending on the target” and advising you to test Vina against your own system with known actives before trusting it.
Three specific problems break naive ranking:
- Size bias. Vina’s scoring function sums per-atom contributions, so a larger ligand collects a more negative total almost automatically. A 45 heavy-atom molecule beating a 20 heavy-atom molecule by 1.5 kcal/mol may simply be bigger, not better.
- A single rigid snapshot. Standard docking treats the receptor as rigid. A pose that fits one crystal conformation can be unreachable once side chains move.
- Score resolution. Differences of a few tenths of a kcal/mol between poses are inside the noise of the function. Treating 9.4 as meaningfully better than 9.2 is reading precision that is not there.
If you have not yet worked through what the numbers do and do not mean, read our guide on how to interpret molecular docking results before going further, and check that your protocol passed redocking validation. A protocol that cannot reproduce a known crystal pose within about 2 Å cannot be trusted to rank unknowns.
What is ligand efficiency, and how do you calculate it?
Ligand efficiency (LE) normalises binding energy by molecular size, which removes most of the size bias in one step. Hopkins, Groom and Alex introduced it in Drug Discovery Today in 2004 in Ligand efficiency: a useful metric for lead selection.
The calculation is arithmetic you can do in a spreadsheet:
LE = -ΔG / N
where ΔG is the predicted binding energy in kcal/mol (your Vina score) and N is the number of heavy atoms, meaning every non-hydrogen atom in the ligand. RDKit gives you N in one line with Descriptors.HeavyAtomCount(mol), and PyMOL or Open Babel will count them too.
A frequently used reference point in the medicinal chemistry literature is an LE near 0.3 kcal/mol per heavy atom, which is roughly what a 500 Da compound binding at 10 nM works out to. Treat it as a sanity line, not a pass mark. The useful comparison is between your own compounds, ranked by LE instead of by raw score.
One honest caveat to state in your methods: LE inherits every error in the docking score. It fixes the size bias, it does not turn an approximate score into a measured affinity.
Which filters should you apply before spending MD time?
Filters come before simulation, because a week of GPU time on a compound that can never be a drug is a week you do not get back. Two filters do most of the work.
Drug-likeness. Lipinski’s rule of five (Lipinski et al., Advanced Drug Delivery Reviews, 1997) flags poor oral absorption: molecular weight above 500, logP above 5, more than 5 hydrogen bond donors, more than 10 hydrogen bond acceptors. SwissADME from the Swiss Institute of Bioinformatics computes these free from a SMILES string, along with TPSA and synthetic accessibility. Our walkthrough of drug-likeness and ADMET for docked ligands covers the full workflow.
PAINS. Pan-assay interference compounds hit in many unrelated assays for reasons that have nothing to do with your target. Baell and Holloway catalogued the substructures in Journal of Medicinal Chemistry in 2010 (DOI: 10.1021/jm901137j). RDKit ships about 480 of these patterns in three families, usable through FilterCatalogParams.FilterCatalogs.PAINS, and SwissADME reports PAINS alerts in its output. A PAINS hit is not automatic disqualification, but if it is in your top three you must say so in the thesis and explain why you kept it.
Natural product screens are where this bites hardest. Curcumin, many flavonoids and several quinones carry PAINS substructures, which is exactly why they appear near the top of so many phytochemical docking studies. Examiners in this field know the pattern.
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.
How many compounds should you carry into MD?
Three, in most MSc projects. The reasoning is practical rather than theoretical.
One compound gives you nothing to compare against, so any result you report is unanchored. Two gives you a winner with no sense of spread. Three lets you show a ranking, and if the top two swap places between docking and MD, that inversion is itself a finding worth a paragraph. Beyond three, a 100 ns protein-ligand simulation each starts to dominate a one-year timeline, and the marginal insight is small.
Pick the three deliberately, rather than taking the top three rows of the LE column:
- The best LE compound that passes both filters.
- A compound from a different chemotype, even at slightly lower LE. If two scaffolds behave the same way in MD, your conclusion is much stronger.
- A reference: the co-crystallised ligand, or a known inhibitor of your target from ChEMBL. This is the single most useful control in the whole study, because it calibrates every number you report. Without it, an MM-PBSA value of -35 kcal/mol means nothing.
What four MD criteria actually decide the winner?
Run each candidate through the same protocol, same force field, same length, same water model, then judge on four things. Our protein-ligand MD tutorial in GROMACS covers the setup.
| Criterion | What it tells you | Tool | A strong result looks like | What it cannot tell you |
|---|---|---|---|---|
| Ligand RMSD stability | Whether the docked pose survives, or drifts and leaves | gmx rms with the ligand fitted on protein backbone | Rises early, then plateaus and stays flat for the production run | Nothing about affinity. A stable pose can still be a weak binder |
| Contact occupancy | Which interactions persist rather than appear once | gmx hbond, plus distance analysis or MDAnalysis | Key hydrogen bonds or salt bridges held for most of the trajectory | Strength of each contact, only its persistence |
| MM-PBSA binding energy | An end-state estimate of relative binding free energy | gmx_MMPBSA | Clear separation between candidates, with your reference in a sensible position | Absolute affinity. Entropy is usually omitted, so values are systematically too negative |
| Per-residue decomposition | Whether the energy comes from the residues you claim it does | gmx_MMPBSA decomposition output | Top contributing residues match the known catalytic or binding site | Causation. It describes the trajectory you ran, not the mechanism |
Criterion four is the one students skip and examiners ask about. If your MM-PBSA total is excellent but the energy is coming from a surface residue nowhere near the active site, the number is an artefact of the pose, not evidence of inhibition. We cover the workflow in per-residue energy decomposition with gmx_MMPBSA, and the totals in binding free energy with MM-PBSA in GROMACS.
On MM-PBSA itself, be careful how you word the claim. Genheden and Ryde’s review in Expert Opinion on Drug Discovery (2015) is the standard reference, and the consensus it reflects is that these end-state methods are useful for ranking related ligands and unreliable as absolute affinities. Write “ranked highest by MM-PBSA”, never “has a binding free energy of”.
How do you build a defensible final scoring table?
Put every criterion in one table, with the reference compound as a row. The values below are illustrative placeholders to show the shape of the table, not results from a real screen. Replace them with your own.
| Compound | Vina score (kcal/mol) | Heavy atoms | LE | Filters | Ligand RMSD | Key contacts held | MM-PBSA rank | Site residues match |
|---|---|---|---|---|---|---|---|---|
| Reference inhibitor | -8.9 | 32 | 0.28 | Pass | Plateaus | 3 of 3 | 2 | Yes |
| Candidate A | -9.6 | 44 | 0.22 | Pass | Drifts out | 1 of 3 | 3 | No |
| Candidate B | -9.1 | 29 | 0.31 | Pass | Plateaus | 3 of 3 | 1 | Yes |
| Candidate C | -8.7 | 26 | 0.33 | PAINS alert | Plateaus | 2 of 3 | 4 | Partial |
Read that table the way a reviewer would. Candidate A wins on raw score and loses everywhere else, which is the size-bias story in one row. Candidate C has the best LE and a PAINS alert, so it is reported and set aside with a reason. Candidate B is the pick: mid-pack on raw score, best LE among clean compounds, stable pose, contacts held, top MM-PBSA rank, and the energy is coming from the right residues. The reference sitting second confirms the pipeline behaves sensibly.
That last point matters more than the winner. A shortlist is credible when the control lands where a control should.
What do you write in the thesis chapter?
Write the criteria before you write the result, so the ranking reads as a method and not as a justification. Four short paragraphs are enough:
- Selection criteria, stated up front with thresholds and the filters applied.
- The table, with the reference compound included.
- The decision, naming the compound and citing the specific rows that support it.
- Limitations: scoring function approximation, MM-PBSA without entropy, single replica, no experimental validation.
Stating limitations plainly is not a weakness in a computational chapter, it is the thing that separates a careful study from an overclaimed one. Our guide on writing the methods section of a docking study has the phrasing. For where this sits in the wider skill set, see the molecular docking pillar guide and the computational biology skills roadmap.
Frequently asked questions
Is a more negative docking score always better?
Within one compound’s set of poses, yes, the more negative pose is the better-scoring one. Across different compounds it is unreliable, because the score scales with ligand size. Compare with ligand efficiency instead.
What ligand efficiency value counts as good?
There is no universal pass mark. Around 0.3 kcal/mol per heavy atom is the commonly cited reference from the medicinal chemistry literature, but the defensible use of LE is relative ranking within your own compound set against your own reference ligand.
Can I choose the final compound from docking alone, without MD?
You can, and plenty of published studies do, but the claim you are allowed to make shrinks accordingly. Without MD you can say a compound scored best in your docking protocol. You cannot say the interaction is stable. If MD is out of scope, say so in the limitations.
My MM-PBSA ranking disagrees with my docking ranking. Which do I report?
Report both, and treat the disagreement as a result rather than a problem. MM-PBSA samples a conformational ensemble while docking scores one rigid pose, so the MD-derived ranking is usually the better-supported one. Explain the difference instead of hiding it.
How long should the MD simulation be before I compare compounds?
Long enough that the ligand RMSD has plateaued and you can discard the pre-plateau region from the analysis. What matters for a comparison is that every candidate ran the identical protocol and length, so the numbers sit on the same footing.
Do I need a reference compound?
Yes. It is the cheapest credibility you can buy in the whole study. A co-crystallised ligand or a known inhibitor from ChEMBL turns an uncalibrated set of numbers into a ranking someone else can interpret.
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.