How to Run an MD Simulation of a Vaccine Construct-TLR Complex in GROMACS: Setup, Length and What to Report

Run molecular dynamics on a docked vaccine construct and TLR receptor by carrying the top cluster representative forward, rebuilding the docked PDB with pdb2gmx, solvating and neutralizing to physiological ionic strength, then running several independent replicas rather than one long trajectory. Report RMSD after discarding the equilibration window, per-residue RMSF, interface hydrogen bonds and buried surface area.

You finished docking your multi-epitope construct against a Toll-like receptor and you have a ranked list of complexes. A docked pose is a static snapshot from a search that treated most of both partners as rigid. Molecular dynamics is how you find out whether that interface survives when both proteins are free to move in water.

This post covers the decisions that belong to the vaccine project: which pose to simulate, what the docked file does to your topology, how long to run, and which numbers a reviewer will expect. It does not re-teach GROMACS. If you have never run a simulation before, work through the molecular dynamics with GROMACS pillar on a single small protein first, then come back with a complex.

Which docked pose should you carry into MD?

Carry forward the representative structure of the top-ranked cluster, not the single lowest-energy model in the results list. Cluster rank and score rank are different columns and they frequently disagree.

ClusPro’s own help page is blunt about this. It states that “the best way to rank models is by cluster size, which is how the models are ranked coming out of Cluspro”, and separately warns users “not judge models based on these scores because that is not what the scoring was designed for”. The energy column exists to drive the search, not to rank the output. Students routinely report the lowest energy value in a paper and it is the wrong column.

HADDOCK users face the same choice, and the published precedent from our own group is explicit. In the monkeypox multi-epitope study (Akhtar et al., Viruses 2022, PMID 36366472), “the top ranked cluster containing the lowest HADDOCK score was selected as a final structure”. Cluster first, score as the tie-breaker inside that cluster.

Two practical checks before you commit a pose to a week of compute:

• Are the epitopes at the interface? If the receptor is contacting your adjuvant and linkers while every predicted epitope points into solvent, the pose is not the one your paper is about, whatever its rank.
• Does the top cluster stand apart from the second? Compare the two by member count and by interface composition. There is no published cutoff that turns a cluster size into a verdict, so do not invent one. A top cluster that barely beats the runner-up is a result worth reporting honestly.

What breaks when you feed a docked complex to pdb2gmx?

A docked output file is not a crystal structure, and the failures are predictable. The GROMACS manual’s common errors page names most of them, and a two-chain docked complex hits a specific subset.

Chain identifiers. Docking servers write the receptor and ligand as separate chains, and merged output files often reuse a letter or split a chain across the file. GROMACS reports “Chain identifier ‘X’ was used in two non-sequential blocks”. The manual’s fix is exactly what it sounds like: move the inserted molecule so it is not splitting another, or rename the second chain to a unique identifier.

Non-standard residues and heteroatoms. If your receptor came from a co-crystal structure, it arrived with ligands, glycans and ions. The human TLR4 ectodomain entry 3FXI carries six non-polymer entities alongside its protein chains. GROMACS answers with “Residue ‘XXX’ not found in residue topology database”. Strip anything you are not deliberately parameterizing before you build the topology.

Hydrogens. Docked models carry whatever hydrogen naming their upstream tool used, and it rarely matches the force field’s expectations. You get “WARNING: atom X is missing in residue XXX Y in the pdb file” or “Atom X in residue YYY not found in rtp entry”. The manual’s recommendation is to use -ignh so that pdb2gmx strips and rebuilds the hydrogens itself.

Missing residues. If your receptor is a crystal structure with unresolved loops, those residues are absent, and GROMACS reports long bonds or missing atoms. The manual is unambiguous that there is no GROMACS tool to fix this: “look for REMARK 465 and REMARK 470 entries in the pdb file. These atoms will have to be modeled in using external software.” It is equally clear that reaching for -missing “is almost always inappropriate” for a protein, because the resulting topology “is likely physically unrealistic”. Model the gaps before you simulate, or pick a more complete structure. If your construct itself came from a structure predictor rather than a template, check its per-residue confidence before you trust the terminal regions.

Which force field and water model should you choose?

Choose a force field that has published use for protein-protein complexes and state the exact variant in your methods. “AMBER” is not a force field name; AMBER99SB-ILDN is.

Both published multi-epitope studies from our team used Amber ff99SB-ILDN (Lindorff-Larsen et al. 2010, PMID 20408171) with TIP3P water, which is a defensible pairing for this class of system and gives you a citable precedent when a reviewer asks why.

GROMACS ships several AMBER variants natively and the force fields page lists them, along with two constraints worth knowing before you choose. If you use AMBER19SB, the manual instructs you to pair it with the OPC or OPC3 water model rather than TIP3P. And CHARMM36 is not bundled: the files come from the MacKerell lab, and the manual prescribes an exact .mdp block to go with them, including vdw-modifier = force-switch, rvdw = 1.2, rvdw-switch = 1.0 and rlist = 1.2. Copying CHARMM36 parameter files without those settings is a silent error, not a loud one.

Avoid the GROMOS family for new work. The same page carries a warning that these force fields “have been parametrized with a physically incorrect multiple-time-stepping scheme for a twin-range cut-off”, and that physical properties such as density may differ from intended values under modern integrators.

How big does the box need to be, and what goes in it?

Big enough that the complex never sees its own periodic image. If the box is too small, GROMACS stops with “The cut-off length is longer than half the shortest box vector or longer than the smallest box diagonal element”, and the manual explains why this matters: the box must be “at an absolute minimum twice the cut-off length in all three dimensions”. Under pressure coupling the box shrinks and grows, so leave margin rather than sitting on the limit.

A TLR ectodomain plus a multi-epitope construct is a large system, and the published numbers make that concrete. The monkeypox study solvated its complex in a cubic box with a 10 Angstrom solvent layer and ended up with 228,291 water molecules, 615 potassium ions and 612 chloride ions, for a total of 701,665 atoms. The dengue study (Kaushik et al., Frontiers in Immunology 2022, PMID 35799781) used a 12 Angstrom layer and 224,932 water molecules. Plan your compute allocation around a system of that size, not around the 30,000-atom lysozyme box from an introductory tutorial.

Neutralize the charge, then add salt on top of it. Both studies added potassium ions to neutralize and then extra potassium and chloride using the Joung-Cheatham ion model to reach a bulk ionic strength of 0.15 M. Adding only enough ions to zero the charge leaves you simulating a construct in pure water, which is not the environment you are claiming to model. If the charge that GROMACS reports is not close to an integer, that is a topology problem rather than an ion problem, and the manual tells you where to look: the cumulative charge column in the pdb2gmx atom listing should hit an integer after every residue.

How long should the simulation run, and how many replicas?

Run several independent replicas rather than one long trajectory, and start each from different initial velocities. Both published studies from our team ran three independent 100 ns simulations per complex, which is 300 ns of total sampling.

The reason is visible in the results rather than the methods. In the monkeypox study the three runs “produced very similar results” overall, but the construct in trial 2 behaved differently: the paper reports that its RMSD values “are still not stable during the 100 ns simulation run, which may indicate a longer simulation run is required for attaining stability”. A single trajectory would have given you either the reassuring answer or the worrying one, with no way to tell which you had drawn. Replicas are what convert “my complex was stable” into a claim with evidence behind it.

Equilibrate in stages before you start counting production time. The published protocol restrains the protein heavy atoms through minimization and both equilibration phases, then releases them:

Every value in that table comes from the methods sections of the two papers cited above. Do not copy a simulation length or an RMSD threshold out of a blog post, including this one, without a source you can cite. There is no universal “run it for 100 ns” rule; there is a precedent for systems of this kind, and precedent is what you cite.

Which analyses actually belong in your results section?

Four analyses carry the argument for a receptor complex, and a fifth is worth adding when a reviewer asks what holds the interface together. The GROMACS documentation for each tool is the citable source for what it computes.

Two details decide whether these numbers mean anything.

Discard the equilibration window before you quote an average. The monkeypox complex “stabilizes after 30ns simulation time”, and the paper reports its RMSD for the post-30 ns portion: 1.08 plus or minus 0.1 nm. Averaging from time zero would have folded the settling transient into the number. Identify where your own plateau begins, say so in the text, and use -b to start the analysis there.

Split the complex into groups first. The interesting quantities are per chain and across the interface, which means building index groups with gmx make_ndx before running any analysis. Reporting one RMSD for the merged system hides the pattern that matters. In the monkeypox study the receptor alone averaged 0.52 plus or minus 0.04 nm while the construct sat at 1.17 plus or minus 0.09 nm, and the per-chain RMSF told the same story: the TLR5 receptor averaged 0.22 plus or minus 0.12 nm, against 0.41 plus or minus 0.20 nm for the construct’s N-terminal region. A flexible construct on a rigid receptor is the expected shape of the result, and you can only show it if you separate the chains.

Also note which receptor the precedent used. Both of these studies docked and simulated against TLR5, not TLR4. If your own work targets TLR4, cite these papers for the workflow and say plainly which receptor they used. Quietly transplanting a TLR5 result onto a TLR4 claim is the kind of error a reviewer in the field catches immediately.

For a deeper argument about what holds the interface together, gmx energy gives you the short-range Coulomb and Lennard-Jones terms between two groups. The monkeypox study reported an electrostatic contribution of -714.65 plus or minus 90.9 kcal/mol against a Lennard-Jones contribution of -665.89 plus or minus 61.5 kcal/mol, and concluded that both components contribute almost equally. That is a claim about relative contributions within one system, not a binding affinity, and it should be written that way.

What goes wrong during the run, and how do you fix it?

Docked complexes fail in a narrower set of ways than a general MD system, because the starting structure is a prediction rather than an experiment.

• LINCS warnings, then the run stops. The manual describes this pattern: constraint algorithms “are usually the first thing to fail” when a system is blowing up, and it “does not necessarily mean you need to troubleshoot the constraint algorithm”. Look at the starting structure. Bad contacts at a docked interface are the usual cause, and another round of minimization is the usual fix.
• “1-4 interaction not within cut-off”. The manual’s guidance here is worth repeating exactly, because the instinct is wrong: “This is BAD. Most importantly, do not increase your cut-off!” It means atoms have acquired very large velocities. If it happens early, the system was not equilibrated well enough.
• “Number of coordinates in coordinate file does not match topology”. The usual cause is a [ molecules ] section that was not updated after solvating or adding ions. The counts and the order must match the coordinate file exactly.
• Position restraint files out of order. With two chains you have two topology includes and two restraint files, and each #include "posre_X.itp" must sit inside its own [ moleculetype ] block. Grouping all the restraint includes together at the end produces “Atom index n in position_restraints out of bounds”, and the manual shows the right and wrong orderings side by side.
• A complex that drifts apart across the periodic boundary in your trajectory viewer. This is almost always a visualization artifact rather than a physical event. Post-process with gmx trjconv to make molecules whole and re-center on the complex before you look at it or measure anything from it.
• Separate thermostats for every molecule type. The manual calls this “a particularly bad practice employed by many new users” and recommends one group, or at most Protein and Non-Protein. Do not give your receptor and your construct their own temperature coupling groups.

What can you claim from a 100 ns trajectory, and what can you not?

You can claim that the docked interface persisted under the conditions you simulated, for the length you simulated, across the replicas you ran. That is a real and useful result, and it is the standard evidence in this literature.

You cannot claim binding affinity, protection, or immunogenicity. A stable trajectory says the pose is not an artifact of a rigid-body search. It says nothing about whether the construct raises antibodies, which is a wet-lab question. Our dengue work went that extra step: one of the mapped epitope peptides was synthesized and tested in a rabbit model, and the in-vivo results are what support the immunogenicity claim in that paper. The simulation supported the structural claim only. You can see the shape of that evidence chain across the studies collected on our research page.

Write your methods so someone can repeat them. Name the GROMACS version (the current release series documented at manual.gromacs.org is 2026.3), the force field variant, the water model, the box shape and solvent layer, the ion model and ionic strength, the total atom count, the thermostat and barostat with their coupling times, the integrator and time step, the cut-offs, the replica count and the analysis windows. Every one of those appears in the two papers cited here, which is why their protocol was reusable enough to tabulate in this article.

This is the last computational step in the vaccine design chain. If you arrived here without the earlier ones, the sequence runs from epitope prediction through antigenicity, allergenicity and toxicity screening, then construct assembly with linkers and an adjuvant, then docking, then this. The immunoinformatics pillar holds that workflow together, the molecular docking pillar covers scoring and search in general terms, and if you are still assembling the underlying skills, the computational biology skills roadmap puts them in order.

Frequently asked questions

Can I run this simulation on a laptop?

Not comfortably. A TLR ectodomain complexed with a multi-epitope construct came to 701,665 atoms in the published monkeypox setup, and three 100 ns replicas of a system that size is a cluster or GPU-server workload. Test your entire pipeline end to end on a short run of a few hundred picoseconds first, so that a setup error surfaces in minutes rather than after a week of queue time.

Do I need to simulate the free receptor and the free construct as well?

It strengthens the argument. Both published studies plot RMSD and RMSF for the receptor and the construct separately from within the complex trajectory, which is cheaper than running additional simulations and answers most of the question. Separate simulations of the unbound partners are a reasonable addition if reviewers ask about conformational change on binding.

My RMSD is still climbing at 100 ns. What do I report?

Report exactly that. The monkeypox paper states plainly that one of its trials was “still not stable during the 100 ns simulation run, which may indicate a longer simulation run is required”. A trajectory that has not converged is a finding, and describing it accurately is better science than truncating the plot at the point where the line looks flat.

Which is the right thermostat and barostat for this system?

The published protocol here used velocity rescaling for temperature and Parrinello-Rahman for pressure during production. The GROMACS documentation notes that velocity rescaling is similar to Berendsen coupling but adds a stochastic term that “ensures that a proper canonical ensemble is generated”, which is why plain Berendsen coupling is not appropriate for a production run.

Can I skip docking and simulate the two proteins side by side until they bind?

No, not on this timescale. Binding events for a protein-protein complex are far outside what a few hundred nanoseconds of unbiased simulation will sample. Docking generates the hypothesis and MD tests whether it holds. They are different steps and neither substitutes for the other.

How do I count hydrogen bonds only at the interface?

Build index groups for the receptor and the construct with gmx make_ndx, then pass them to gmx hbond as the reference and target selections with -r and -t. The tool’s documentation gives 0.35 nm as the minimum and recommended donor-acceptor distance and 30 degrees as the recommended angle cut-off, which are the geometric criteria you should state in your methods.

The workflow described here comes from the StemSkills Lab team, which has more than ten years of combined work in sequence and structural bioinformatics, drug discovery and design, and multiscale molecular modeling, including the peer-reviewed multi-epitope vaccine studies cited throughout this article.