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
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
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
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:
Note how many lipids were deleted and update the
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.
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