Chapter 11: Molecular Dynamics Simulations with GROMACS
Chapter 11: Molecular Dynamics Simulations with GROMACS
While docking provides a static snapshot of binding, biological macromolecules are dynamic. Molecular Dynamics (MD) simulations numerically integrate Newton's equations of motion (\(\vec{F} = m \vec{a}\)) across femtosecond timesteps, validating whether a docked drug pose remains stable in an aqueous physiological environment over nanoseconds to microseconds.
The 5 Sequential Stages of an MD Simulation
- Topology & Force Field Assignment: Assign atomic masses, charges, and bond parameters (e.g. CHARMM36m or AMBER14SB).
- Solvation & Ion Neutralization: Place protein in a TIP3P water box and add \(\text{Na}^+ / \text{Cl}^-\) ions to reach physiological 0.15 M salinity.
- Energy Minimization (EM): Relax steric clashes using steepest descent.
- Equilibration: NVT (constant Volume and Temperature, 300 K) followed by NPT (constant Pressure, 1 bar).
- Production MD: Unrestrained simulation run (e.g. 100 ns) saving trajectory coordinates.
1. End-to-End GROMACS Simulation Pipeline
# 1. Generate GROMACS topology from PDB
gmx pdb2gmx -f protein.pdb -o processed.gro -p topol.top -ff amber99sb-ildn -water tip3p
# 2. Define rhombic dodecahedron solvent box with 1.0 nm clearance
gmx editconf -f processed.gro -o boxed.gro -c -d 1.0 -bt dodecahedron
# 3. Solvate with water
gmx solvate -cp boxed.gro -cs spc216.gro -o solvated.gro -p topol.top
# 4. Add ions to neutralize charge at 0.15 M NaCl
gmx grompp -f ions.mdp -c solvated.gro -p topol.top -o ions.tpr
echo "SOL" | gmx genion -s ions.tpr -o neutralized.gro -p topol.top -pname NA -nname CL -neutral -conc 0.15
# 5. Energy Minimization (EM)
gmx grompp -f minim.mdp -c neutralized.gro -p topol.top -o em.tpr
gmx mdrun -v -deffnm em
# 6. NVT and NPT Equilibration (100 ps each)
gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr
gmx mdrun -deffnm nvt
gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -p topol.top -o npt.tpr
gmx mdrun -deffnm npt
# 7. Production MD (100 ns run on GPU)
gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -o md_100ns.tpr
gmx mdrun -deffnm md_100ns -nb gpu -pme gpu2. Analyzing Trajectory Metrics (RMSD & RMSF)
- Root Mean Square Deviation (RMSD): Measures overall structural stability compared to the starting frame. A plateau indicates an equilibrated system.
- Root Mean Square Fluctuation (RMSF): Measures per-residue flexibility, highlighting rigid core \(\alpha\)-helices vs flexible outer binding loops.
# Calculate backbone RMSD over time
gmx rms -s md_100ns.tpr -f md_100ns.xtc -o rmsd.xvg -tu ns