GROMACS Tutorial

Step Three: Defining the Unit Cell & Adding Solvent

Defining a unit cell for a membrane protein is considerably more complicated than for a protein in water. There are several key steps in building the unit cell:

  1. Orient the protein and membrane in the same coordinate frame
  2. Pack the lipids around the protein
  3. Solvate with water

1. Orient the protein and membrane

We have already aligned the KALP peptide using editconf. The bilayer lies in the x-y plane, with the normal along the z-axis. Convert the dppc128.pdb to .gro format with editconf and remove the initial periodicity. The latter step is easily accomplished through the use of trjconv, following these steps:

(1) Generate a .tpr file for a DPPC-only system using grompp. You can use any valid .mdp file and a topology corresponding to pure DPPC. An example .mdp file can be found here and such a topology can be found here. Note how simple the topology is. It includes dppc.itp and spc.itp to read parameters for DPPC and water; it's that simple! Run grompp:

gmx grompp -f minim.mdp -c dppc128.gro -p topol_dppc.top -o em.tpr

You may receive a fatal error like this one, but in this case it is safe to use -maxwarn 1 with the above command. The reason that doing so is acceptable is explained here. Please note that this step is the only one in which topol_dppc.top will be used. It serves no other purpose and you should not use it in any remaining step in this tutorial.

It is not necessary to run a complete energy minimization procedure on the bilayer, although you can if you want. The .tpr file contains information about bonding and periodicity, so it can, in a sense, be used to reconstruct "broken" molecules.

(2) Use trjconv to remove periodicity:

gmx trjconv -s em.tpr -f dppc128.gro -o dppc128_whole.gro -pbc mol -ur compact

Now, take a look at the last line of this .gro file; it corresponds to the x/y/z box vectors of the DPPC unit cell. We need to orient the KALP peptide within this same coordinate frame, and place the center of mass of the peptide at the center of this box:

gmx editconf -f KALP-15_processed.gro -o KALP_newbox.gro -c -box 6.41840 6.44350 6.59650

The center of our system now lies at (3.20920, 3.22175, 3.29825), half of each box vector. This is a GROMACS convention. Note that other systems you wish to simulate may not be symmetrical with respect to the membrane, and thus the above command must be modified to something like the following:

gmx editconf -f protein.gro -o protein_newbox.gro -box (membrane box vectors) -center x y z

In the above command, x y z represents the center of mass such that the protein is properly placed. Placement should be based on experimental knowledge of membrane positioning, or intuition based on the chemical composition of your particular protein.


2. Pack the lipids around the protein

The easiest method I have found so far for packing lipids around an embedded protein is the InflateGRO methodology (ref), available here. Please note that I am distributing my own copy of the original version of InflateGRO, not InflateGRO2 available elsewhere from the authors. Download the linked file and rename it "inflategro.pl" to continue. First, concatenate the protein and bilayer structure files:

cat KALP_newbox.gro dppc128_whole.gro > system.gro

Remove unnecessary lines (the box vectors from the KALP structure, the header information from the DPPC structure) and update the second line of the coordinate file (total number of atoms) accordingly.

The authors of the InflateGRO script recommend using a very strong position-restraining force on protein heavy atoms to ensure that the position of the protein does not change during EM. Add a new #ifdef statement to your topology, one that will call a special set of position restraints, such that your topology now contains a section like:

; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif

; Strong position restraints for InflateGRO
#ifdef STRONG_POSRES
#include "strong_posre.itp"
#endif

; Include DPPC chain topology
#include "dppc.itp"

; Include water topology
#include "gromos53a6_lipid.ff/spc.itp"

Now, generate this new position restraint file using genrestr:

genrestr -f KALP_newbox.gro -o strong_posre.itp -fc 100000 100000 100000

In the .mdp file used for the minimizations, add a line "define = -DSTRONG_POSRES" to make use of these new position restraints. Then, simply follow the InflateGRO instructions (contained in the script itself), a procedure that is easily scripted. Scale the lipid positions by a factor of 4:

perl inflategro.pl system.gro 4 DPPC 14 system_inflated.gro 5 area.dat

Since the InflateGRO script requires a very specific order to its command-line arguments, a brief explanation is warranted here. They are:

  1. system.gro - the input coordinate file name to which scaling will be applied
  2. 4 - the scaling factor to apply. A value > 1 indicates inflation, a value < 1 indicates shrinking/compression
  3. DPPC - the residue name of the lipid to which scaling is applied
  4. 14 - the cutoff radius (in Å) for searching for lipids to delete
  5. system_inflated.gro - the output file name
  6. 5 - grid spacing (also in Å) for area per lipid calculation
  7. area.dat - output file with area per lipid information, useful in assessing the suitability of the structure

Note how many lipids were deleted and update the [ molecules ] directive of your topology accordingly. Run energy minimization. Then, scale down the lipids by a factor of 0.95 (assuming you have used default names, the result of the minimization is called "confout.gro"):

perl inflategro.pl confout.gro 0.95 DPPC 0 system_shrink1.gro 5 area_shrink1.dat

Follow this up by another round of EM. During the "shrinking" steps, be sure to change the cutoff value to 0, or else you will continue to delete lipids! After 26 iterations of scaling down by 0.95, I reached an area per lipid of ~71 Å2, above the experimental value of ~62 Å2. Since the script tends to overestimate the area per lipid, this value is good enough to continue to equilibration.


3. Solvate with water

Solvating a protein is a trivial task with gmx solvate. Solvating a membrane protein system is not so simple, since solvate has a tendency to fill gaps in the lipid acyl chains with random water molecules. The quickest way I have found to solvate membrane systems is to make a local copy of vdwradii.dat and change the value of C from 0.15 to 0.375. Doing so will cause solvate to assign carbon atoms an artificially large van der Waals radius, thus making addition of water within the lipids less likely. After solvating, check your structure to be sure no water molecules are present in the hydrophobic core of the bilayer. If there are a few stray molecules, you can either delete them manually, continue to adjust the van der Waals radius of carbon, or use one of the scripts found at the GROMACS site.

Be aware that using a large radius for carbon may also create artificial voids around proteins that protrude outside the membrane (since they also contain carbon) as well as the lipid headgroups (if they contain carbon, i.e. PC, PE, etc). These voids may require a bit of finesse during equilibration. Delete the local copy of vdwradii.dat before continuing.

Back: The Topology Next: Adding Ions

Bevan Lab Homepage

Virginia Tech Homepage

Virginia Tech Biochemistry


Site design and content copyright 2008-2015 by Justin Lemkul
Problems with the site? Send them to the Webmaster