How to Prepare a Protein Structure for MD Simulation (Missing Residues, Hydrogens, Protonation)

To prepare a protein for molecular dynamics, clean the raw PDB before you build a topology: rebuild missing residues and loops, remove crystallographic waters unless they are functional, resolve alternate conformations, then add hydrogens and assign protonation states for histidine, aspartate, glutamate and cysteine at your target pH. Skipping any of these steps can crash the run or produce meaningless results.

You have downloaded a structure from the RCSB Protein Data Bank, or you have just built a fresh model with AlphaFold or homology modelling. It looks complete in the viewer, so it is tempting to feed it straight into pdb2gmx or a docking box. Do not. A raw coordinate file is almost never ready for simulation. This guide walks an MSc student through every cleaning step, in order, with the real tools and the specific residue cases that quietly ruin a run. It sits inside the wider computational biology skills roadmap, one step after you validate the structure and one step before you build a topology and solvate.

Why does a raw PDB file need preparation before MD?

A crystal structure is a model fitted to experimental density, and that density is incomplete. Flexible loops often have no visible electron density, so the crystallographer leaves those residues out. Hydrogen atoms are almost always absent because X-rays scatter weakly off a single electron, so hydrogens only become visible at ultra-high resolution, better than about 1.2 Angstrom, which very few structures reach. The file also carries crystallisation artefacts: ordered waters, buffer ions, cryoprotectants, and sometimes two modelled positions for the same side chain.

A molecular dynamics force field, on the other hand, needs a chemically complete and unambiguous molecule: every heavy atom present, every hydrogen placed, one conformation per atom, and a defined charge on every ionisable group. Hand it a structure with a broken backbone or a wrong histidine charge and one of two things happens. The topology builder fails with an error about a missing atom, or worse, it runs, and you spend a week simulating a protein that carries the wrong electrostatics. The second failure is the dangerous one because it looks like success.

What is the difference between missing atoms and missing residues?

Missing atoms means a residue is present but incomplete, usually a side chain whose tip was disordered, so a few atoms past the C-beta are absent. Missing residues means the residue is gone entirely, backbone and all, leaving a gap in the chain. The Protein Data Bank records both. In a PDB file, REMARK 465 lists whole residues that were not modelled, and REMARK 470 lists individual atoms that are missing from residues that are otherwise present. The full intended sequence lives in the SEQRES records, so comparing SEQRES against the coordinates tells you exactly what is absent.

How do you find and rebuild missing residues and loops?

Start by reading REMARK 465 and REMARK 470 in the header, or open the mmCIF version and inspect the pdbx_unobs_or_zero_occ_residues category, which is the modern equivalent. Short internal gaps of a few residues usually need rebuilding because a broken backbone cannot be simulated as one chain. Long, highly flexible terminal tails are a judgement call: rebuilding a 30-residue disordered tail adds a floppy region that contributes noise and can dominate your trajectory, so many workflows simply cap the chain at the last resolved residue rather than model coordinates the experiment never saw.

Two tools do most of this work. PDBFixer, part of the OpenMM project, adds missing residues and heavy atoms automatically and is ideal for short gaps and a quick, reproducible clean-up. For longer loops where geometry matters, MODELLER builds and optimises multiple loop conformations against spatial restraints. MODELLER is described by its authors as a program for “comparative modeling of protein three-dimensional structures”, and its loop-refinement routine is the standard way to rebuild a missing segment when you care about the conformation. The original method is set out by Sali and Blundell in their 1993 paper in the Journal of Molecular Biology.

After rebuilding a loop, always re-check the geometry of the new segment. A modelled loop can contain steric clashes or backbone outliers, so run it back through the same stereochemical checks you use to validate any protein structure before trusting it.

Should you keep or remove crystallographic waters, ions, and hetero atoms?

Everything that is not a standard amino acid sits in HETATM records: waters (residue name HOH), ions, ligands, cofactors, and crystallisation additives. The default is to strip most of them, but do it deliberately, one class at a time, rather than deleting all hetero atoms in one blind command.

• Crystallographic waters: remove them by default. Your MD box will be filled with a fresh solvent model, so keeping crystal waters is usually redundant. The exception is a water that mediates a specific, functional interaction, for example bridging a ligand to the active site. Keep those few, delete the rest.
• Structural and catalytic ions: keep them. A zinc in a zinc-finger, a magnesium coordinating a phosphate, or a calcium in an EF-hand is part of the folded protein. Removing it changes the structure and the chemistry. Deleting a sodium or chloride that only sat in the crystallisation buffer is fine.
• Ligands and cofactors: decide by intent. An apo simulation strips the ligand; a holo simulation keeps it, but the ligand then needs its own force-field parameters, which standard protein topology builders do not provide.
• Cryoprotectants and buffer molecules: glycerol, PEG, sulfate, and similar additives are almost always artefacts. Remove them.

How do you handle alternate conformations?

When a side chain was modelled in two positions, the PDB stores both with an alternate-location indicator in column 17, usually A and B, each with a partial occupancy that sums to roughly 1.0. A force field cannot place one atom in two spots, so you must choose one conformation, normally the higher-occupancy A, and remove the alternate. Many preparation tools do this automatically, but confirm it happened. A leftover alternate location is a common cause of a pdb2gmx error about a duplicate atom, and it is easy to miss because the structure still looks fine on screen.

How do you add hydrogens, and why are they missing?

Hydrogens are missing because standard-resolution crystallography cannot see them, yet they carry much of the hydrogen bonding and the entire notion of protonation, so the force field needs every one. There are two routes to add them.

The first is to let your topology builder do it. In GROMACS, gmx pdb2gmx adds hydrogens as it builds the topology, and by default it ignores any hydrogens already in the file so it can place them consistently with the chosen force field. The second route is to add and optimise hydrogens with a dedicated tool such as PDB2PQR before you build the topology, which is the better choice when you also want careful protonation-state assignment, covered next. Whichever route you take, the placement of a hydrogen depends entirely on the protonation state of its residue, so hydrogen addition and protonation are really one decision, not two.

How do you assign correct protonation states at physiological pH?

Every ionisable side chain has a pKa, and whether it is protonated at your simulation pH depends on how that pKa shifts inside the folded protein. Tools predict these shifted pKa values from the structure and assign a protonation state for each residue. Two widely used options are PDB2PQR, which wraps the empirical predictor PROPKA, and the H++ server, which uses a continuum-electrostatics Poisson-Boltzmann model. The PROPKA method behind PDB2PQR is described by Olsson and colleagues in their 2011 paper in the Journal of Chemical Theory and Computation, and the H++ approach is documented by Anandakrishnan and colleagues in Nucleic Acids Research (2012).

The reason this matters is that a residue’s charge is set by comparing its pKa to the pH. The imidazole side chain of histidine has a model pKa near 6.0, which sits close to physiological pH 7.4, so histidine is the single most ambiguous residue in any protein: small environmental shifts can flip it between neutral and positively charged, and the neutral form itself has two tautomers. Get it wrong and you place a proton on the wrong ring nitrogen or add a full positive charge that should not be there, right next to a binding site.

Which residues cause the most trouble?

Four cases account for most silent errors:

• Histidine: choose between HID (proton on the delta nitrogen), HIE (proton on the epsilon nitrogen), and HIP (both, positively charged). Base the choice on the local hydrogen-bonding network or metal coordination, not on a coin flip. A histidine coordinating a catalytic metal is usually the neutral tautomer that points its lone pair at the metal.
• Cysteine in a disulfide: two cysteines in an S-S bond must be declared as bonded and are not protonated. In the AMBER naming scheme they become CYX. If you leave them as ordinary reduced cysteines, the topology builder adds two thiol hydrogens that should not exist and the disulfide is lost.
• Buried acidic residues: an Asp or Glu inside a hydrophobic pocket or forming a strong salt bridge can carry a proton at neutral pH (ASH or GLH). Assuming every acidic residue is negative is a common shortcut that changes the electrostatics of the active site.
• Termini: the free N-terminus and C-terminus are ionisable too. Standard preparation caps them as charged (NH3+ and COO-), which is correct for a full-length chain but wrong if you are simulating a fragment that should be capped with a neutral acetyl or amide group.

In GROMACS you can set these interactively. Running gmx pdb2gmx with the flags -his, -asp, -glu, -lys, and -ter makes the tool ask you for each ambiguous residue rather than guessing, which is the safest way to apply the states a predictor like PROPKA recommended.

What is a clean order of operations?

Do these steps in order. Changing the order, for example adding hydrogens before you have rebuilt a gap, forces you to redo work.

1. Inspect REMARK 465 and REMARK 470, and compare coordinates against SEQRES, to list every missing residue and atom.
2. Rebuild internal gaps and missing heavy atoms with PDBFixer or MODELLER; cap or model long disordered termini deliberately.
3. Remove alternate conformations, keeping the higher-occupancy set.
4. Strip crystallographic waters, cryoprotectants, and buffer ions; keep structural or catalytic ions and any functional water.
5. Predict pKa values with PROPKA/PDB2PQR or H++ and decide the state of every His, Asp, Glu, Cys, and terminus.
6. Add hydrogens consistent with those states, then build the topology and continue to solvation and energy minimisation.

A model from AlphaFold or homology modelling skips steps 1 to 4, since it has no missing residues, no alternate conformations, and no crystal waters, but it still needs protonation assignment and hydrogen placement, and its low-confidence regions deserve the same geometric scrutiny as a rebuilt loop.

This whole clean-up is the difference between a trajectory a reviewer trusts and one they reject. The StemSkills Lab team, with more than a decade in sequence and structural bioinformatics, drug discovery and design, and multiscale molecular modeling, has seen the same handful of preparation mistakes sink otherwise good projects, and almost every one traces back to a skipped step above.

Frequently asked questions

Can I just run pdb2gmx on a raw PDB and skip preparation?

Sometimes it runs, which is exactly the trap. pdb2gmx will fail on a missing atom or a leftover alternate conformation, but it will happily assign default protonation states and ignore an unrebuilt loop by treating the chain as broken. You get a topology, but it can encode the wrong charges or a discontinuous backbone, so prepare the structure first.

Do I need to rebuild every missing loop?

No. Rebuild internal gaps, because a broken backbone cannot be simulated as one continuous chain. For long, disordered terminal tails that the experiment never resolved, capping the chain at the last ordered residue is often more honest than modelling coordinates you cannot justify.

How do I decide the histidine tautomer?

Look at the local environment. Use a predictor such as PROPKA or H++ for a first pass, then inspect the hydrogen-bonding partners and any nearby metal. Put the proton on the nitrogen that donates a hydrogen bond, or use the charged HIP form only when the residue clearly accepts a proton, for example near two acidic groups.

What pH should I prepare the structure at?

Use the pH your system is meant to represent. Physiological pH 7.4 is the common default for a cytoplasmic or extracellular protein, but a lysosomal enzyme works near pH 4.5 to 5, and that lower pH changes which acidic residues are protonated. Set the pH in your pKa tool to match the biology.

Does an AlphaFold model still need preparation?

Yes, a different subset. It has no missing residues, alternate conformations, or crystal waters to fix, but it has no hydrogens and no assigned protonation states, so you still run the protonation and hydrogen steps. Treat any low-confidence region with the same care as a rebuilt loop before you simulate it.