A molecular dynamics (MD) simulation produces a trajectory: thousands of snapshots of every atom in your system, saved over nanoseconds or microseconds of simulated time. The simulation itself is the easy part — modern engines automate most of it. The hard part is answering the question that made you run it: did anything meaningful happen?
This guide covers the five analyses that answer that question for almost every protein system, what each one actually measures, and how to tell a real signal from noise.
First, check that your simulation is even valid
Before interpreting anything, confirm the system equilibrated. Two quick checks catch most failures:
- Energy and temperature stability. Potential energy should fall during minimisation and then plateau. Temperature and pressure should fluctuate around their target values, not drift.
- Density. For a solvated system, water density should settle near 1.0 g/cm³ during NPT equilibration. If it does not, your box or ion count is wrong.
Analysing a trajectory that never equilibrated is the single most common way to reach a confident, wrong conclusion. Discard the equilibration phase before you measure anything — typically the first 5–20 ns, judged from the plots below rather than a fixed rule.
1. RMSD — has the structure settled?
Root mean square deviation measures how far the backbone has moved from a reference structure (usually the starting or minimised structure) after removing overall translation and rotation.
What you are looking for is a plateau. A trace that rises quickly and then flattens means the protein relaxed away from the crystal structure and found a stable basin — normal and expected. A trace still climbing at the end of the run means you have not simulated long enough to conclude anything about stability.
Typical values: a well-behaved globular protein settles somewhere between 1 and 3 Å. Above roughly 5 Å, look closely — you may be seeing genuine large-scale motion such as domain opening, or you may be seeing partial unfolding or a bad starting model.
A frequent mistake is comparing RMSD values between different proteins as if they were interchangeable. RMSD scales with system size and flexibility; it is meaningful within a system across conditions, not as an absolute quality score.
2. RMSF — which residues are actually moving?
RMSD collapses the whole protein into one number, which hides where the motion is. Root mean square fluctuation gives you a per-residue value instead, so you can see which parts move and which stay rigid.
Expect high RMSF at termini and in long surface loops — that is normal and not interesting on its own. What is interesting:
- A binding-site residue that becomes rigid when a ligand is bound but flexible in the apo form — direct evidence of induced fit.
- A loop whose flexibility changes between wild type and mutant, pointing to an allosteric effect.
- A secondary structure element with unexpectedly high fluctuation, which often signals local unfolding.
RMSF is most powerful as a comparison. A single RMSF plot tells you the protein has flexible loops; two plots side by side tell you what your mutation or ligand actually did.
3. Radius of gyration — is the protein staying folded?
Radius of gyration (Rg) measures how compact the structure is: the mass-weighted spread of atoms about the centre of mass. It is the cleanest single indicator of global folding state.
A stable, folded protein shows a flat Rg trace with small fluctuations. A steady rise means the structure is expanding — unfolding, or a domain swinging open. A fall can indicate collapse, sometimes an artefact of a poorly parameterised force field or an over-strong protein–water interaction.
Rg is particularly useful in thermal stability work: run the same construct at several temperatures and compare where the trace stops being flat.
4. SASA — what is exposed to solvent?
Solvent accessible surface area measures how much of the protein the solvent can reach. Split into hydrophobic and hydrophilic contributions, it becomes considerably more informative.
Rising hydrophobic SASA is a strong unfolding signal: residues that should be buried in the core are becoming exposed. In protein–ligand work, the buried surface area at the interface is a useful proxy for how extensive the contact is, and it correlates loosely with binding strength.
5. Hydrogen bonds — is the interaction actually holding?
For a protein–ligand or protein–protein complex, hydrogen bond analysis is often the most directly interpretable result you have. Count them over time and you learn whether the interaction you docked is durable or whether it dissolves within a few nanoseconds.
Two things matter more than the raw count:
- Occupancy — the percentage of frames in which a specific bond exists. A bond present in 80% of frames is a real feature of the complex. One present in 10% is a transient contact.
- Which residues — a persistent bond to a catalytic residue means something quite different from a persistent bond to a surface lysine.
This is the analysis that most often overturns a docking result. Docking scores a single static pose; MD asks whether that pose survives contact with physics. Poses that look excellent by docking score and disappear in the first nanosecond are common, which is exactly why the two methods belong together.
Going further: free energy surfaces
The five analyses above describe what happened. A free energy surface asks which conformational states the protein actually occupies and how they relate energetically.
By projecting the trajectory onto two collective variables — commonly the first two principal components from PCA, or a pair of physically meaningful distances — and converting the population of each region into an energy, you get a map of the accessible conformational landscape. Deep, well-separated basins indicate distinct stable states; a single broad basin indicates one flexible state.
This is how you distinguish “the protein is floppy” from “the protein switches between two defined conformations”, which is often the biologically important question.
How long does a simulation need to be?
Long enough for the process you care about, which is not a satisfying answer but is the honest one. As rough guidance:
- Side-chain rearrangement and local relaxation: tens of nanoseconds.
- Loop motions and ligand pose stability: 100–500 ns.
- Domain motions and larger conformational change: microseconds, or enhanced sampling methods such as metadynamics or replica exchange.
Run replicas. Three independent 200 ns simulations from different starting velocities are far more trustworthy than one 600 ns run, because they reveal whether your result is reproducible or an accident of one trajectory. Single-replica conclusions are among the most common weaknesses in published MD work.
A practical analysis checklist
- Confirm energy, temperature, pressure and density are stable; discard the equilibration phase.
- Plot RMSD and identify where it plateaus.
- Plot per-residue RMSF, comparing conditions rather than reading one plot alone.
- Check Rg for global compactness.
- Split SASA into hydrophobic and hydrophilic components.
- For complexes, measure hydrogen bond occupancy per residue pair.
- Run at least three replicas and confirm your conclusion appears in all of them.
- Build a free energy surface if the question involves distinct conformational states.
Interpretation is where MD projects succeed or fail. The simulation generates the data; these analyses are what turn it into a result you can defend in review.
