GROMACS Tutorial

Step Six: Equilibration

Equilibration will be conducted much like in the case of a solvated protein. Generally a short NVT equilibration phase is followed by a longer NPT phase. The reason for this procedure is that we are now dealing with a heterogenous system, with both water and DPPC acting as solvent. Such heterogeneity requires a longer equilibration process. Water has to re-orient around the lipid headgroups and any exposed parts of the protein, and the lipids have to orient themselves around the protein as well. Such processes take some time, and lipid equilibration may take several ns of simulation time.

For membrane protein simulations, we will need to create a special index group consisting of protein and lipids (explained below). To do this, use make_ndx:

gmx make_ndx -f em.gro -o index.ndx
...
 > 1 | 13
 > q

Entering "1 | 13" at the make_ndx prompt merges the Protein (1) and DPPC (13) groups. This new group will be used for center-of-mass motion removal and temperature coupling (more on this shortly).

We will once again start with NVT (with the script found here), calling grompp and mdrun just as we did at the EM step:

gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -n index.ndx -o nvt.tpr

gmx mdrun -deffnm nvt

Most of the parameters we are using are comparable to those in the lysozyme tutorial, with a few changes:

  • rcoulomb, rvdw = 1.2: We are using a 1.2-nm short-range cutoff for electrostatics, and van der Waals interactions, which differs from the original parametrization. The Berger lipids were parametrized with a 1.0-nm cutoff and GROMOS96 with 0.9 nm for electrostatics and 1.4 nm for van der Waals interactions. The value used here has been employed by the author and verified to produce reasonable physical behavior. The value of rcoulomb is somewhat flexible when using PME, but rvdw is very important for balance in force fields. In any case, it is always an exercise for the user to demonstrate validity of any setup, and it is important to consult the relevant literature regarding force field comparisons and systematic evaluations of cutoff values.
  • ref_t, gen_temp = 323: We must use a temperature that is above the phase transition temperature of the lipid. For DPPC, 323 K is commonly used.
  • tc-grps = Protein_DPPC Water_and_ions: The two groups are coupled separately for due to the different rates of diffusion of the two phases. With an aqueous protein, we specify "Protein Non-Protein," with Non-Protein containing solvent and ions. For a membrane protein, Non-Protein would contain the lipids as well, so we must explicitly couple the lipids and aqueous solvent separately. Do not couple the ions separately from the solvent; there are insufficient degrees of freedom in an ion-only group to justify a separate temperature coupling group.
  • A new section pertaining to center-of-mass (COM) motion removal. Since interfacial systems (i.e., membrane-water systems) have a tendency to move laterally, the motion of the bilayer COM and solvent COM must be reset separately. Otherwise, the phases may drift in opposite directions, such that the overall COM for the system is unchanged, but artificial motion is still present. Note that our comm-grps include Protein_DPPC; since the protein is embedded in the membrane, it is essentially part of the membrane. Removing its COM motion separately may lead to spurious collisions. It is for this reason that we do not remove COM motion in separate groups for proteins in water; diffusion occurs in 3 dimensions in these systems. In a bilayer, motion is largely restricted to 2-dimensional motion.

Again, use gmx energy to confirm that the temperature of the system has stabilized at 323 K before continuing. The choice of temperature should be based on the physical properties of the lipid, most notably the phase transition temperature. Some useful data that has been mined from the literature regarding area per lipid, phase transition temperatures, and/or parameter derivation for various lipids is presented here (see below). Please read the references and understand their implications. This list is not comprehensive; the references point to examples of literature wherein these particular lipids were used either for simulation or experimental work. The user should investigate further citations within these works, or subsequent papers that have cited these. The list should also not be viewed as comprehensive, as many other lipids have been successfully simulated. Presented here are some of the more common ones.

Lipid Name Area Per Lipid (Å2) Phase Transition (K) Reference

DPPC 62.9 - 64 315 J. Nagle (1993) Biophys. J. 64: 1476

DMPC 60.6 297 Wohlert and Edholm (2006) J. Chem. Phys. 125: 204703

POPG 53 269 Dickey and Faller (2008) Biophys. J. 95: 2636

POPA 51-52 301 Dickey and Faller (2008) Biophys. J. 95: 2636

POPC 65.8 271 Tieleman, et al. (1998) Biochem. 37: 17554

POPE 56 298 Tieleman, et al. (1998) Biochem. 37: 17554

DMTAP 71 310 Gurtovenko, et al. (2004) Biophys. J. 86: 3461

POPS 55 300 Mukhopadhyay, et al. (2004) Biophys. J. 86: 1601

Having trouble with your system? Is it crashing during equilibration? Please see the Advanced Troubleshooting page.


Back: Energy Minimization

Next: Equilibration, Phase 2


Site design and content copyright Justin Lemkul

Problems with the site? Send them to the Webmaster