Wild-Type vs Mutant MD Simulation in GROMACS

A valid wild type versus mutant MD comparison changes one thing only: the mutated residue. Run both systems through an identical protocol, use five or more replicas each, compare per-residue RMSF rather than global RMSD, and report means with error bars across replicas. Overlapping error bars mean no claim.

This guide is part of our molecular dynamics with GROMACS series, written by the StemSkills Lab team, who have spent more than 10 years in sequence and structural bioinformatics, drug discovery and design, and multiscale molecular modeling. It assumes you already know how to mutate a residue in PyMOL and how to run a single simulation. The hard part is not either simulation. The hard part is the design that lets you say the mutation did something.

This is the most common MSc and PhD project shape in Indian departments: an nsSNP from ClinVar, a drug resistance mutation, a thermostability variant. The site already teaches every individual analysis. What follows is the part nobody writes down, which is how to assemble those analyses into a comparison that an examiner or a reviewer cannot dismiss.

What makes a wild-type versus mutant comparison valid?

One rule governs everything: the only difference between your two systems is the residue you changed. Every other difference is a confound, and a confound is indistinguishable from the effect you are claiming.

That means both systems share, identically:

• the same force field and the same water model
• the same box type and the same padding, set the same way with editconf
• the same salt concentration and the same ion-placement procedure when you add ions
• the same minimisation and the same NVT and NPT equilibration protocol, from the same .mdp settings
• the same production length, the same number of replicas per system, and the same seed-generation procedure
• the same analysis groups and, critically, the same fitting selection

Build the mutant from the same starting coordinates as the wild type, not from a separate PDB entry or a separate model. If the wild type needed missing loops rebuilt, rebuild them once, then mutate the finished structure and send both through pdb2gmx. Two independently prepared structures differ in dozens of ways before the mutation is even considered.

The confound almost everyone hides: net charge

A charge-changing mutation, say glutamate to lysine, shifts the net charge of the solute by two units. GROMACS neutralises the box for you, so the two systems end up with different counterion counts. This is unavoidable and it is not a mistake. Hiding it is the mistake.

Keep the salt concentration identical, let the neutralising counterion count differ as the chemistry demands, and state both ion counts in your methods. An examiner who finds an unreported ion-count difference will reasonably ask whether your effect is electrostatic bookkeeping rather than biology. An examiner who sees it reported and discussed has nothing to attack.

How many replicas does each system need?

This decides whether you have a result or an anecdote, and it is where most student projects fail.

The reason is physical, not statistical pedantry. Two runs of the same system, differing only in initial velocities, diverge. In Avoiding False Positive Conclusions in Molecular Simulation: The Importance of Replicas (Journal of Chemical Theory and Computation, 2018, volume 14, pages 6127 to 6138), Knapp, Ospina and Deane examined this directly across 310,000 ns of simulation time. They ran 100 identically parametrised replicas of 3,000 ns each for a 10 amino acid system, and 100 identically parametrised replicas of 100 ns each for an 827 residue T-cell receptor and MHC complex, then sampled random subgroups to see what a given number of replicas actually buys you.

Their findings are blunt: conclusions from single simulations are often not reproducible, and several shorter replicas are more reliable than one longer run. Their stated rule of thumb is a minimum of five to 10 replicas.

Read that against your own project. If the spread between two replicas of your wild type is as large as the wild type versus mutant difference you are about to report, you have measured noise. One run each is not a comparison, it is two samples from two distributions with no estimate of either. Treat the per-system replica count as an explicit design parameter, decided before you launch anything, and see our guide on how many replicate MD simulations to run in GROMACS for how to generate and manage them.

A practical floor for a thesis: five replicas per system, same length, with the replica count and the total aggregate time stated in your methods. Two systems at five replicas is 10 simulations, which is why the design matters before the compute does.

Which observables should you actually compare?

Not all of them equally. Each observable answers a different question, and each carries a standard misreading that reviewers look for.

Global RMSD is the weakest evidence in the set, and it is the one students lead with. The Lysozyme in Water GROMACS tutorial by Justin Lemkul warns on its analysis page that RMSD cannot be used to assess convergence or judge stability, calling it “a degenerate metric of structural change” that is also extrinsic. Two different conformational changes can give the same RMSD, and a larger protein accumulates a larger number from smaller per-atom deviations. A difference in mean global RMSD between your systems tells you something changed somewhere, which is not a finding.

Per-residue RMSF is the strongest and most interpretable signal, because it localises the effect. If the mutation matters, the fluctuation profile should change near the mutation site, or in a region you can argue is coupled to it. gmx rmsf computes per-residue averages with the -res option, which is off by default, and can write the fluctuations as B-factors into a PDB file with -oq so you can colour the structure by mobility.

The PCA trap in that table deserves its own warning. If you run PCA separately on each system, the first two eigenvectors describe different motions, so the two free energy landscapes are plotted on axes that are not the same axes. The comparable version is to build the covariance matrix on a combined or common reference and project both trajectories onto that single set of eigenvectors. The same logic applies to DCCM and to gmx cluster: a shared basis, or the comparison is decorative.

For the remaining mechanics, including the exact commands and flags, the existing guides cover them: RMSD and RMSF, radius of gyration, SASA and hydrogen bonds, secondary structure with gmx dssp, and MM-PBSA binding free energy. One current detail worth knowing before you script anything: the GROMACS manual notes that gmx hbond is a new implementation added in GROMACS 2024, driven by selections through -r and -t, and that the previous interactive version is still available as gmx hbond-legacy. If a tutorial from 2019 does not match your prompts, that is why.

How do you decide whether the difference is real?

This is the section that separates a defensible project from a long description of two trajectories, and it is the part almost no page on this topic writes.

The procedure is simple to state. For every observable you report, compute the per-replica value, then the mean across replicas for each system, then an error estimate on that mean. Compare the two means with their error bars. If the error bars overlap, you have no claim, and you say so.

GROMACS can do the averaging for you. gmx analyze reads the .xvg files your analysis tools wrote, always reports the average and standard deviation of each set, and offers two things you want here. The -av option produces the average over sets with error bars added through -errbar, where the error bars can represent the standard deviation, the error assuming the points are independent, or the interval containing 90% of the points. The -ee option produces error estimates by block averaging, dividing the set into blocks and computing the error on the total average from the variance between block averages. Block averaging matters because consecutive frames in a trajectory are correlated, so a naive standard error computed over every frame is far too small and will manufacture significance that is not there.

On why this is not optional, Grossfield and Zuckerman set out the statistical tools for exactly this problem in Quantifying uncertainty and sampling quality in biomolecular simulations (Annual Reports in Computational Chemistry, 2009, volume 5, pages 23 to 48), which asks how one determines the statistical significance of observed results and concludes that such analyses are of paramount importance in establishing the reliability of simulation data in any given study.

There is a related habit to break. Judging equilibration by looking at an RMSD plot is not reliable either. Knapp, Frantal, Cibena, Schreiner and Bauer tested this with a survey of practising simulation scientists, shown randomised RMSD plots, and reported in Is an intuitive convergence definition of molecular dynamics simulations solely based on the root mean square deviation possible? (Journal of Computational Biology, 2011, volume 18, pages 997 to 1005) that there was no mutual consent about the point of equilibrium and that the judgements were severely biased by plot parameters. Their conclusion is that equilibration should not be discussed on the basis of an RMSD plot at all. Decide your equilibration and discard windows by a stated rule applied identically to both systems, before you look at the result.

The negative result is a result

If your five wild-type replicas and your five mutant replicas give overlapping RMSF profiles everywhere including the mutation site, the correct sentence is that the mutation did not measurably alter the backbone dynamics of the fold on this timescale with this sampling. That is a finding. It is examinable, it is honest, and it is far stronger than a claimed destabilisation that lives inside the replica spread. Students fabricate effects because they believe a null outcome fails the thesis. Overclaiming is what fails the thesis.

How do you report a wild-type versus mutant comparison?

Report in three layers, and build them in this order.

The figure set. Per-residue RMSF for both systems on one axis, mean across replicas, shaded error band, with the mutation site marked. That is your lead figure because it carries the localisation argument. Then the site-specific evidence: hydrogen bond or salt bridge occupancy around the mutated residue. Then the collective picture: both systems projected onto a shared PCA basis. Global RMSD belongs in the supplementary material, as evidence that both systems behaved, not as evidence of the effect. Our guide to plotting GROMACS .xvg output in Python covers the mechanics of overlaying two systems with error bands.

The table. One row per observable, one column per system, each cell a mean plus an error estimate across replicas, and a final column saying whether the difference is resolved by your sampling. Say “not resolved” where it is not resolved.

The methods sentences. Name the force field and water model, the box type and padding, the salt concentration and both neutralising ion counts, the minimisation and equilibration protocol, the production length, the number of replicas per system, the equilibration window you discarded and the rule you used, the fitting selection used for every RMSD and RMSF calculation, and the error estimator. Our MD methods section guide has the sentence patterns, and the GROMACS manual is the citation for every tool default you relied on, including the current release documentation at version 2026.3.

What goes wrong, and how do you fix it?

Six failures account for most of the wasted months. The first three are silent, which is what makes them dangerous.

The two systems have different atom counts, so the fits are not comparable

Error or symptom: gmx rms refuses to proceed with mismatched groups, or worse, accepts your selection and quietly fits different atom sets in the two systems. Remember that gmx rms asks you for two selections, one group for the least-squares fit and then a group for the RMSD calculation, and -fit defaults to rot+trans. Choosing “Protein” in each run selects different atoms if the mutation changed side-chain size.

Fix: build a matched index group once, with gmx make_ndx, containing only atoms present in both systems, usually backbone or C-alpha. Use that same group as the fit group and the calculation group in both systems, and pass it with -n.

A charge-changing mutation gave different ion counts, and it is not in the methods

Fix: keep the salt concentration identical between systems, accept the different neutralising counterion count, and report both. Do not quietly add ions to force the counts to match, because that changes the ionic strength instead, which is a worse confound.

The RMSD curves sit at different offsets because the references differ

Symptom: the mutant’s RMSD starts higher and stays higher by a constant amount. This is usually an artefact of comparing each trajectory to its own -s structure, where the two reference structures are not equivalent. Fix: decide one reference convention, either each system against its own equilibrated starting structure or both against a common experimental structure, apply it to both, and state which you used. The constant offset usually disappears.

Per-residue RMSF profiles will not align

Symptom: the two RMSF curves are shifted by one or more residues, or the mutation site lands at a different x value. The cause is residue renumbering, introduced by the mutation step or by a missing-loop rebuild. Fix: renumber both structures to one common scheme before production, not after analysis, and verify by aligning the two sequences rather than trusting the index.

The mutant looks less stable, but the difference is inside the replica spread

Fix: compute the between-replica spread for the wild type alone first. If that spread is comparable to the wild type versus mutant difference, the difference is not resolved by your sampling. Either add replicas or report the null outcome. Use gmx analyze -ee so the error estimate accounts for correlated frames.

A sudden RMSD jump appears in one system only

Symptom: a step change in RMSD at one time point, in one system, with no corresponding change in RMSF or secondary structure. This is usually a periodicity artefact rather than an event. Fix: reimage the trajectory with gmx trjconv using -pbc mol -center and redo the analysis. Then confirm with gmx mindist -pi, which plots the minimum distance of a group to its periodic image, considering one shift in each direction for a total of 26 shifts. If the molecule met its image, the box was too small and that system needs rerunning with the padding you used for the other one. Our list of common GROMACS errors covers the message-level failures.

Frequently asked questions

Can I compare a wild type and a mutant simulated for different lengths?

No, not as a primary claim. Different production lengths mean different sampling, so any difference you see is partly a difference in how much conformational space each system explored. Truncate both to the shorter length, apply the same equilibration discard, and treat the extra time as supplementary.

Do both systems have to use the same force field?

Yes, and the same water model, the same ion parameters, and the same cutoffs. Changing the force field changes the dynamics far more than most point mutations do, so a mixed comparison measures the force field, not the mutation.

How many replicas are enough for an MSc thesis?

Five per system is a reasonable floor, consistent with the minimum of five to 10 replicas recommended by Knapp, Ospina and Deane. Report the count and the aggregate simulated time. If compute limits you to fewer, state that explicitly as a limitation rather than presenting single runs as settled.

Is global RMSD useless then?

Not useless, just misused. It is a sanity check that each system stayed folded and reached a plateau. It is poor evidence that a mutation changed something, because it is extrinsic and degenerate, so two different structural changes can produce the same value.

Can MM-PBSA prove a mutation weakens ligand binding?

It can support that claim, not prove it. MM-PBSA values carry large variance, so compare the difference between systems with error bars across replicas rather than absolute numbers, and pair it with per-residue energy decomposition so the effect is localised to residues you can name.

What if the mutation site is disordered or in a rebuilt loop?

Then your localisation argument is weak, because the local flexibility you measure is partly a property of your model-building rather than of the protein. Say so, and lean on observables away from the rebuilt region, or choose a different variant for the project.

Where to go next

The comparison, not the simulation, is the skill being examined. Match the protocol, run replicas on both systems, lead with per-residue RMSF rather than global RMSD, and let error bars decide what you are allowed to claim. If you need a project that fits this shape, our MD project ideas for an MSc thesis list several variant comparisons, and the computational biology skills roadmap shows where this sits in the wider sequence of skills. If the mutation you want to study still needs choosing, start with predicting the effect of a mutation on protein stability and simulate the one with a prediction worth testing.