How to Analyze Protein-Ligand Interactions Over an MD Trajectory (ProLIF and MDAnalysis, Step by Step)

To analyze protein-ligand interactions over an MD trajectory, load a PBC-corrected GROMACS trajectory into MDAnalysis, build a ProLIF Fingerprint, and run it over the frames. Converting the result to a DataFrame and taking the column mean gives contact occupancy: the percentage of frames each binding-site residue held each interaction type.

Almost every docking-plus-MD thesis reaches the same moment. You have your RMSD plot, your RMSF plot, maybe an MM-PBSA number, and then your supervisor asks the question that none of those answer: which residues actually held the ligand, and for how long? A docking pose gives you a snapshot of contacts that may not survive the first nanosecond. This guide covers the analysis that answers the question properly, written by the StemSkills Lab team from 10+ years in sequence and structural bioinformatics, drug discovery and design, and multiscale molecular modeling.

What is an interaction fingerprint, and why not just draw a LigPlot+ diagram?

An interaction fingerprint encodes, per frame, whether each ligand-residue pair satisfies the geometric definition of each interaction type. The ProLIF paper defines it directly: “Interaction fingerprints are vector representations that summarize the three-dimensional nature of interactions in molecular complexes, typically formed between a protein and a ligand” (Bouysset and Fiorucci, Journal of Cheminformatics 13:72, 2021).

The word doing the work is per frame. LigPlot+ (Laskowski and Swindells, JCIM 51:2778-2786, 2011) and PLIP (Adasme et al., Nucleic Acids Research 49:W530-W534, 2021) both produce excellent diagrams, but of one structure. Run either on your docked pose and you get the contacts the docking program thought existed. Run it on the last frame instead and you get a different set, with no way to tell which one is representative. The fingerprint replaces that guess with a count over the whole run.

ProLIF is the library built for exactly this. As of version 2.2.1 it ships nine interaction types in the default fingerprint, and the 2021 paper has been cited more than 400 times (Crossref, August 2026), so it is a defensible citation in a methods section.

What do you need before you start?

Three things, and skipping any of them produces silently wrong numbers rather than an error.

1. A PBC-corrected trajectory. If your ligand hops across the periodic box, ProLIF sees it kilometres from the protein and records zero contacts. Fix it first with the standard three-step chain covered in our gmx trjconv guide: -pbc whole, then -pbc nojump, then -center -pbc mol -ur compact.
2. Explicit hydrogens. ProLIF infers bond orders and formal charges from connectivity, and the documentation is blunt about the requirement: the inference “requires all atoms including hydrogens to be present in the input file”. A heavy-atom-only trajectory will not give you correct hydrogen bond donors.
3. An elements attribute. The RDKit converter underneath ProLIF needs element symbols, not force-field atom types. A GROMACS .tpr is a good topology choice here because the TPR parser reads bonds directly from the force field, and it exposes elements when at least one atom carries a valid symbol.

How do you install ProLIF?

Conda-forge is the recommended route in the official installation docs, because it resolves the RDKit dependency for you:

conda install -c conda-forge prolif

The pip route works too, but you must ask for RDKit explicitly:

pip install rdkit prolif

ProLIF 2.2.1 requires Python 3.10 or newer. Check what you got before going further:

import prolif as plf
print(plf.__version__)

How do you load the trajectory and select the ligand?

ProLIF sits on top of MDAnalysis (Michaud-Agrawal et al., Journal of Computational Chemistry 32:2319-2327, 2011), so the loading step is the same one you already know from our MDAnalysis trajectory analysis tutorial:

import MDAnalysis as mda
import prolif as plf

u = mda.Universe("md.tpr", "md_center.xtc")

ligand_selection = u.select_atoms("resname LIG")
protein_selection = u.select_atoms("protein")

print(ligand_selection.n_atoms, protein_selection.n_atoms)

Replace LIG with your own residue name. If you are not sure what it is, print the set of residue names in the system rather than guessing:

print(set(u.atoms.resnames))

On a large system you can restrict the protein side to the neighbourhood of the ligand and save a lot of time. The ProLIF documentation gives the selection and the warning that goes with it, that the shell is chosen “on the first frame” only:

protein_selection = u.select_atoms(
    "protein and byres around 20.0 group ligand", ligand=ligand_selection
)

A 20 angstrom shell around a ligand that stays in its pocket is safe. If your ligand is expected to unbind or crawl along a surface, keep the full protein selection.

How do you compute the interaction fingerprint?

Two lines. Build the fingerprint object, then run it over the trajectory slice you want:

fp = plf.Fingerprint()
fp.run(u.trajectory[::10], ligand_selection, protein_selection)

The [::10] takes every tenth frame. Start there while you are debugging, then drop to [::1] for the run that goes in the thesis. A default plf.Fingerprint() tracks nine interaction types: Hydrophobic, HBDonor, HBAcceptor, PiStacking, Anionic, Cationic, CationPi, PiCation and VdWContact. The halogen-bond and metal-donor classes exist but are not in that default set, and the water-bridge class is hidden unless you ask for it. To list everything the library can detect:

print(plf.Fingerprint.list_available())
print(plf.Fingerprint.list_available(show_bridged=True))

You can pass your own list, which is the sensible move when van der Waals contacts flood the output:

fp = plf.Fingerprint(
    ["Hydrophobic", "HBDonor", "HBAcceptor", "PiStacking", "Anionic", "Cationic"]
)
fp.run(u.trajectory, ligand_selection, protein_selection, n_jobs=4)

The n_jobs argument parallelises across frames and is the single easiest speed-up on a multi-core laptop.

How do you turn the fingerprint into a contact-occupancy table?

This is the figure your thesis actually needs. Convert to a DataFrame, then take the mean down each column. Because the default fingerprint stores booleans, the mean of a column is the fraction of frames in which that interaction was present, and multiplying by 100 gives occupancy as a percentage:

df = fp.to_dataframe()
(df.mean().sort_values(ascending=False).to_frame(name="%").T * 100)

Read the output as a ranked list. A hydrogen bond to ASP129 present in 78% of frames is a real, reportable finding. The same hydrogen bond at 4% is noise that your docking pose happened to catch. This one table separates the two, and it is the sentence you can defend in a viva.

Three follow-up views are worth knowing. Drop the interaction type that dominates by sheer count, keep a single residue, or keep a single interaction type:

df.drop("Hydrophobic", level="interaction", axis=1)
df.xs("ASP129.A", level="protein", axis=1)
df.xs("PiStacking", level="interaction", axis=1)

For the visual version, the barcode plot puts frames on one axis and residue-interaction pairs on the other, so a contact that dies at 20 ns is immediately visible as a band that stops:

fp.plot_barcode()

How does this differ from the gmx hbond count?

Both are legitimate, and they answer different questions. gmx hbond applies a geometric definition with a default donor-acceptor cutoff of 0.35 nm and a default acceptor-donor-hydrogen angle cutoff of 30 degrees, and its -num output writes hbnum.xvg, described in the GROMACS manual as the “Number of hydrogen bonds as a function of time”.

That is a total, not an attribution. If the count sits at three for the whole run, you still do not know whether it was the same three hydrogen bonds throughout or a rotating cast. ProLIF answers the attribution question and covers hydrophobic, pi-stacking, ionic and van der Waals contacts that gmx hbond never looks at. The honest thesis reports both: the global count as a stability trace, alongside the occupancy table as the per-residue breakdown. Our guide to radius of gyration, SASA and hydrogen bonds in GROMACS covers the global side.

ProLIF vs PLIP vs LigPlot+ vs gmx hbond: which should you use?

The pragmatic combination for a docking-plus-MD chapter is ProLIF for the occupancy table, then PLIP or LigPlot+ on the representative structure from your clustering or MM-PBSA step for the pretty 2D diagram. The occupancy table tells the reader which contacts to trust; the diagram shows them what those contacts look like.

What goes wrong, and how do you fix it?

The ligand selection returns 0 atoms

Your residue name is not LIG. Print set(u.atoms.resnames) and use the real one. Ligands parameterised through CGenFF or ACPYPE often keep a three-letter code you did not choose.

An error complaining that the elements attribute is missing

The topology gave atom types instead of element symbols. The ProLIF documentation fixes this with MDAnalysis guessers. On MDAnalysis 2.8 and newer the current API is:

u.guess_TopologyAttrs(to_guess=["elements"])

Every occupancy is 0%, or absurdly low

Nine times out of ten this is periodic boundary conditions. Load the same trajectory, jump to a middle frame, and check the distance between the ligand centre of mass and the protein centre of mass. If it is tens of angstroms, go back and rerun gmx trjconv. The other cause is a heavy-atom-only trajectory, where hydrogen bond detection quietly fails.

The DataFrame is enormous and dominated by VdWContact

Van der Waals contacts are in the default set from ProLIF 2.x and they fire constantly. Pass an explicit interaction list without it, or drop the level afterwards with df.drop("VdWContact", level="interaction", axis=1).

Bond orders look wrong for an unusual ligand

The MDAnalysis RDKit converter is documented as experimental, and the docs quantify it: the converter “accurately infers the structures of approximately 99% of the ChEMBL27 dataset”. Metal complexes, unusual tautomers and some charged heterocycles are where the remaining 1% lives. Check the inferred SMILES against your input ligand before trusting the fingerprint.

The run takes hours

Slice the trajectory with u.trajectory[::10], restrict the protein selection to a 20 angstrom shell, and pass n_jobs. A 100 ns trajectory saved every 10 ps is 10,000 frames, which is far more resolution than an occupancy percentage needs.

How do you report this in your methods section?

Name the library, the version, the interaction set, the frames analysed, and the definition of occupancy. One honest sentence is enough: interaction fingerprints were computed with ProLIF 2.2.1 on the PBC-corrected trajectory using MDAnalysis, over every frame of the production run, and contact occupancy was defined as the percentage of analysed frames in which a given residue-ligand interaction was detected. Cite the ProLIF and MDAnalysis papers. Our guide to writing the methods section for an MD study covers the surrounding paragraphs, and the full sequence from setup to analysis sits in the molecular dynamics with GROMACS pillar guide. If you are still deciding which skills to build first, the computational biology skills roadmap shows where this analysis sits in the wider workflow.

Frequently asked questions

What occupancy percentage counts as a “stable” interaction?

There is no universal threshold, so state yours. Many published protein-ligand MD studies treat contacts above roughly 50% of frames as persistent and report everything above about 20% for completeness. Whatever you pick, apply it consistently and put the number in the methods, because the cutoff is a choice and reviewers will ask.

Can I run ProLIF on docking poses instead of a trajectory?

Yes. The ProLIF paper describes the library as generating fingerprints “for molecular complexes extracted from molecular dynamics trajectories, experimental structures, and docking simulations”. The docking route uses an SDF supplier rather than an MDAnalysis Universe, which makes it useful for comparing a ranked pose list against the contacts your MD run actually sustained.

Do I need to strip water before running ProLIF?

Not for the standard analysis, since a protein selection excludes solvent anyway. You do need water if you want water-bridged interactions, which ProLIF provides through a separate WaterBridge class that is not part of the default fingerprint.

Does ProLIF handle periodic boundary conditions for me?

No. Correct the trajectory with gmx trjconv before you load it. This is the most common reason a beginner’s occupancy table comes back empty.

Can I use a PDB file as the topology instead of a .tpr?

You can, but a PDB carries no bond orders or formal charges, so ProLIF has to infer them from connectivity and needs every hydrogen present. A .tpr gives the force-field connectivity directly, which is why it is the safer choice for a GROMACS run.