How to Do Ensemble Docking with MD Snapshots (Receptor Flexibility Beyond Flexible Side Chains)

Answer: Ensemble docking replaces one rigid receptor with several receptor conformations taken from a molecular dynamics trajectory. You simulate the apo protein, cluster the trajectory to pick representative structures, dock your whole ligand library into each one separately, then combine the per-conformation scores into a single ranking.

Why isn’t flexible side-chain docking enough?

Because it moves the wrong things. When you turn on flexible residues in AutoDock Vina, you pick a handful of side chains and let their torsions rotate inside a pocket whose backbone stays frozen. That handles a rotamer flip. It does not handle a loop swinging shut, a helix shifting, or a subpocket that only exists for part of the time.

The AutoDock developers say this plainly in their own flexible docking documentation: “the lack of receptor flexibility is arguably the greatest limitation in these types of docking methods.” Flexible side chains narrow that gap. They do not close it.

There is a measurable cost to getting the receptor conformation wrong. A kinase cross-docking benchmark by Schaller, Christ, Chodera and Volkamer tested what happens when you dock a ligand into a protein structure that was not solved with that ligand, which is the honest version of a real prospective project. Docking each ligand into a single randomly chosen kinase structure reproduced the experimental binding mode within 2 A RMSD only 23.8% of the time with a physics-based docking program. Allowing the same program to try multiple structures and keeping the best pose raised that to 84.3%. The authors’ summary of the trend is direct: “Increasing the number of protein structures used for docking steadily increased the chance to generate a low RMSD docking pose.”

That is the entire argument for ensemble docking in one sentence. If you have only ever run single-structure docking, start with our molecular docking guide and the flexible docking walkthrough before you attempt this.

What is ensemble docking, and where did it come from?

The method has a name and a paper. Lin, Perryman, Schames and McCammon called it the relaxed complex scheme in Journal of the American Chemical Society (2002, volume 124, pages 5632 to 5633). The logic behind it is that a ligand does not always bind the conformation you happen to have crystallised. It may bind a shape the protein only visits occasionally, and molecular dynamics is how you generate those shapes.

The most cited demonstration came two years later. Schames and co-workers docked the inhibitor 5CITEP into snapshots from a 2 ns MD trajectory of HIV-1 integrase and found a trench next to the active site that opens intermittently and was not apparent in the crystal structure. That paper, “Discovery of a Novel Binding Trench in HIV Integrase” (Journal of Medicinal Chemistry, 2004, volume 47, pages 1879 to 1881), is the standard example of MD revealing a pocket that static structures miss.

Amaro, Baron and McCammon later published an improved version of the scheme (Journal of Computer-Aided Molecular Design, 2008, volume 22, pages 693 to 705) that focuses on the practical problem you are about to hit: a trajectory has thousands of frames, and you cannot dock a library into all of them. It compares ways of reducing the trajectory to a small representative set, including evenly spaced snapshots, RMSD-based clustering and QR factorization.

What do you simulate, and for how long?

Simulate the apo receptor, meaning the protein with the ligand removed. If you simulate the holo complex, the pocket stays moulded around the ligand you already have, and every snapshot you extract will be biased toward it. The point of the exercise is to see the shapes the pocket adopts on its own.

Set the system up exactly as you would for any production run: solvate, neutralise, minimise, equilibrate under NVT then NPT, and only then collect the production trajectory. Our GROMACS protein simulation tutorial covers that pipeline, and protein preparation for MD covers the setup decisions that break runs later.

On length, be honest about what you can defend. Sampling is the limiting factor in this whole method, and a short trajectory samples only the motions that are fast enough to appear in it. The HIV integrase result above came from 2 ns, which was substantial in 2004 and is small today. A practical student-scale target is a production run long enough that your backbone RMSD has plateaued and your pocket-lining residues have visibly explored more than one arrangement. Check both before you go further, using the approach in RMSD and RMSF analysis. If the pocket never changes shape across your trajectory, ensemble docking will not help you, and that is a real finding worth reporting rather than a failure to hide.

One caution that saves people weeks: discard the equilibration portion before you cluster. Frames from a system still relaxing toward its equilibrium are not representative conformations. They are artefacts of your starting structure.

How do you turn a trajectory into a handful of receptor conformations?

You cluster it. Clustering groups structurally similar frames together and hands you one representative structure per group, which is exactly the reduction you need.

The decision that matters most is what you fit and measure on. Do not cluster on the whole protein. A trajectory clustered on all backbone atoms is dominated by global motion, and you will get representatives that differ in a distant loop while the binding site stays identical. Fit on the protein backbone to remove overall rotation and translation, then compute the clustering RMSD over the binding-site residues only. That produces representatives that differ where you care.

In GROMACS, this is gmx cluster. Its official documentation lists five algorithms: linkage, jarvis-patrick, monte-carlo, diagonalization and gromos. The gromos method is the usual default choice, and the manual describes it as: “Count number of neighbors using cut-off, take structure with largest number of neighbors with all its neighbors as cluster and eliminate it from the pool of clusters. Repeat for remaining structures in pool.” It cites Daura et al., Angewandte Chemie International Edition, 1999, volume 38, pages 236 to 240. Select the algorithm with -method and the RMSD cutoff with -cutoff, which GROMACS expects in nanometres, not angstroms. The -cl option writes out the representative structures.

In AmberTools, the equivalent is the cluster command in cpptraj, which supports hierarchical agglomerative clustering via hieragglo with epsilon and clusters keywords, as well as dbscan and kmeans. Representative structures come out through repout and repfmt, and cluster statistics through summary and info. The official AMBER clustering tutorial shows the full syntax, including the sieve keyword that subsamples frames so the pairwise RMSD matrix fits in memory.

Choose your cutoff so that you end up with a workable number of clusters that together account for most of your frames. If three clusters cover 90% of the trajectory and the remaining forty cover the rest, dock the three and say so. Our dedicated walkthrough on extracting representative structures with gmx cluster goes through the cutoff scan in detail.

How do you dock the library into every snapshot?

Mechanically, this is ordinary docking repeated in a loop. Conceptually, there is one rule that people break constantly.

Keep the grid box identical across every receptor conformation. Same centre, same size. If you recentre the box on each snapshot’s pocket, the conformations are no longer scored on a common footing, and comparing their scores becomes meaningless. Define the box once on a reference structure using the method in our grid box guide, then reuse those numbers everywhere. Make the box slightly generous, because a pocket that opens in one snapshot needs room.

The rest is bookkeeping. Prepare each representative structure as a receptor file the way you normally would, keep one output directory per conformation, and record which snapshot produced which score. Raise exhaustiveness above the default, because you are going to compare scores between runs and you cannot afford search noise to look like a conformational effect. The loop structure itself is the same one in our batch docking tutorial, with the outer loop over receptors instead of a single fixed one.

Budget the cost before you start. A library of 500 ligands against 8 conformations is 4,000 docking runs, not 500. Our note on hardware requirements will tell you whether your laptop can carry that.

How do you combine the per-conformation scores?

Each ligand now has one score per receptor conformation. There are three defensible ways to collapse that into a ranking, and they answer different questions.

• Best score across the ensemble. Take each ligand’s most favourable score from any conformation. This asks “can this ligand bind any accessible shape of the pocket?” It is the most sensitive and the most permissive, and it is what the kinase benchmark above was measuring when the success rate climbed to 84.3%.
• Mean score across the ensemble. Average over conformations. This asks “does this ligand bind consistently, not just to one lucky snapshot?” It suppresses single-conformation artefacts and is generally the more conservative ranking.
• Weighted average. Weight each conformation by the fraction of trajectory frames its cluster contained, so rare conformations count less than common ones. This is closer in spirit to the original relaxed complex scheme, which treats the ensemble as a population rather than a list.

Report the ranking you used and show the spread. A ligand that scores well on every conformation is a different proposition from one that scores well on a single rare snapshot, and the second is not automatically wrong. The HIV integrase trench was exactly that kind of transient opportunity. What you must not do is compute all three, pick whichever promotes your favourite compound, and report only that. Decide the combination rule before you look at the results and state it in your methods section.

What usually goes wrong?

• Every cluster representative looks the same. Your cutoff is too large, or you clustered on the whole protein instead of the binding site. Rerun the clustering with a tighter cutoff and a binding-site atom selection.
• Hundreds of clusters, each with a few frames. The opposite problem, and sometimes a sign the trajectory has not converged. Loosen the cutoff, and check whether your production run is long enough.
• Scores swing wildly between conformations for every ligand. Check that the grid box is genuinely identical across runs, and that exhaustiveness is high enough to be reproducible. See exhaustiveness and reproducibility.
• The GROMACS cutoff behaves strangely. Almost always the nanometre versus angstrom mistake. A cutoff of 2 in gmx cluster means 2 nm, which is 20 A, and will merge your entire trajectory into one cluster.
• No way to tell whether the ensemble helped. Redock a known ligand into every conformation and check which ones recover the crystallographic pose, following our redocking validation protocol. If one representative recovers a pose that the crystal structure alone did not, you have your answer.

Frequently asked questions

How many receptor conformations should I dock into?

There is no fixed number. Use however many clusters are needed to account for the bulk of your trajectory, typically a handful rather than dozens for a student project. Report the cluster populations alongside the count so a reader can judge the coverage themselves.

Can I do ensemble docking with AutoDock Vina?

Yes. Vina has no built-in ensemble mode, but ensemble docking is simply the same docking repeated across several receptor files with an identical grid box, then combined afterwards. Any docking program you can script in a loop will do it.

Should I still use flexible side chains inside each snapshot?

You can, and the two are complementary: the MD snapshots supply backbone and loop movement, flexible residues add local rotamer adjustment. It multiplies your runtime, so add it only for the specific residues you have a reason to suspect, not by default.

Is ensemble docking the same as the relaxed complex scheme?

They overlap heavily. The relaxed complex scheme is the named protocol from Lin and McCammon that formalised docking into MD-derived conformations. “Ensemble docking” is the broader term, and the ensemble can also come from multiple crystal structures or NMR models rather than a simulation.

Do I need a GPU for this?

For the MD stage, effectively yes for anything beyond a very short run. The docking stage is CPU-bound and parallelises trivially across snapshots, so a multi-core machine handles it. Free GPU time through GROMACS on Google Colab is enough to get a student project started.

Where this fits in your skill set

Ensemble docking sits at the junction of two skills you build separately: running a stable MD simulation and running a reliable docking protocol. It is worth attempting once you are comfortable with both, and it makes a strong MSc project because the workflow is defensible, the failure modes are diagnosable, and the result is interpretable either way. Our computational biology skills roadmap shows the order to pick up the prerequisites, and the molecular docking pillar collects the docking side of it.

This guide was written by the StemSkills Lab team, whose members have more than ten years of combined research experience in sequence and structural bioinformatics, drug discovery and design, and multiscale molecular modeling.