Suggestions for creating a helix in octane simulation
system in Gromacs.
1. Get a helix
First make a helix pdb-file. Either
take it from a membrane protein, make one with a molecular graphics
program, or use X-plor scripts kicking around in the
group to build an ideal helix. I have used init.in, ann.in, dyn.in and noe.intra.in
for this. Edit init.in
and noe.intra.in
and run Xplor on init.in, then on ann.in, and
finally on dyn.in.
2. Make an .itp file for your helix
Run pdb2gmx on your pdb file as usual, e.g.
pdb2gmx -f helix.pdb
-ter
When using dummies (see below), use instead:
pdb2gmx -f helix.pdb
-ter -dummy hydrogens -heavyh
When you also have aromatic groups, use:
pdb2gmx -f helix.pdb
-ter -dummy aromatics -heavyh
Take the resulting .top file
and delete the appropriate lines: basically the first couple of lines and
everything after the bit that includes position restraints. Do have a look at
existing .itp
files if you're not sure, or ask someone. Change the name of the molecule in
the .itp
file from Protein to something descriptive, say BR_helix1 for the first helix
of bacteriorhodopsin. You’ll need to run pdb2gmx, and
have a separate .itp file for every chain of your
protein. Make a topology file something.top
that looks like:
; topology for alm in octane plus water, G96
#include "ffgmx.itp"
#include "almH.itp"
#include "octane.itp"
#ifdef FLEX_SPC
#include "flexspc.itp"
#else
#include "spc.itp"
#endif
[ system ]
; name
Alm octane, starting topology
[ molecules ]
; name number
AlmH 1
OCTA 182
SOL 1618
AlmH is the name of the helix, which would be
BR_helix1 for the example above. The number of OCTA and SOL molecules might
need to be changed later on, depending on the procedure you're following. where you adjust the number of OCTA and SOL molecules later
on. If you use the old Gromacs forcefield, ffG43a2.itp would be ffgmx.itp, and the g96_*
files would be replaced by their non-g96 counterparts.
3. Making your system
This involves a number of steps, which can be done by a combination of editconf, genbox and a
bit of editing. Summarized, the steps are:
- Decide on the size of the box
you want. 4.5 x 4.5 x 6.5 nm is a nice starting value for a single helix.
4.0 x 4.0 x 6.5 nm is probably the minimal size that is still useful.
- Decide on the thickness of
your octane slab. This depends on the hydrocarbon thickness of the
biological membrane you are interested in. 2.5 - 3.0 nm are reasonable
values for most purposes.
- Make an octane slab with the
required thickness and x and y size, set the boxsize to the size you want for the whole system, and
translate it to the center of your box.
- Put your helix in a box the
size you want, aligned along z, centered in the middle of the box.
- Put the helix in the octane
slab.
- Add water.
- If necessary, add ions.
- Energy minimize.
Let's take an example (Alamethicin – PDB-id: 1AMT).
I'll assume you have a helix alm.pdb, which is running along some random diagonal
in a rectangular box. You want an octane slab of size 4 x 4 x 3 nm, in a box of
4 x 4 x 6.4 nm. The hydrophobic length of Alm is
roughly 3 nm. I already have .itp files for Alm and
for octane, and by chance the example topology file given above works for Alm.
Align the helix along the z-axis. This is easily done with
editconf -princ -f alm.pdb -o aligned_helix.pdb.
Set the box to 4 x 4 x 6.4 nm and move the helix to the center of the box.
This can be done with:
editconf -f aligned_helix.pdb -o boxed_helix.pdb -box 4 4 6.4 -center 2 2 3.2
Now make the octane slab. Because octane molecules are rather large and
we're going to throw away all octane molecules with any overlap with the helix,
we need a few spare ones. Instead of being smart, I just take a slightly
thicker slab, say 3.2 nm thickness, and see how many octanes I'm left with. (As a guide, the slab will shrink by
~10% during equilibration
- 10% in z if only the dimensions of the box are allowed to
change.) This depends also on the x,y size of the system. If in the
end there are too many octanes (hydrophobic layer
gets too thick), just delete a few more molecules. Adding a few molecules is
much harder. To make the slab, we need an equilibrated octane box, like
spc216.gro for water. You can use octane.g96.
In theory, the following should work:
genbox -cs octane.g96 -o octane_slab.gro -box 4 4 3.2
editconf -f octane_slab.gro -o octane_slab.pdb
The extra step from .gro
to .pdb
is to work around a bug in genbox
(fixed in newer versions of Gromacs).
Note that files ending in .g96
are GROMOS96 coordinate files, look at them with an
editor or pager for more information.
Alternatively, use a premade octane slab (just add
water!). Here are some sizes that might be useful:
Now use editconf:
editconf -f octane_slab.pdb -o octane.pdb -box 4 4 6.4 -center 2 2 3.2
Next use genbox
to put the helix inside the octane, using the correctly oriented helix as
solute and the correctly oriented octane slab as solvent:
genbox -cs octane.pdb
-cp boxed_helix -o helix+octane.pdb
Use genbox
again to add water as usual:
genbox -cs spc216 -cp helix+octane.pdb -o system.pdb
Adjust your something.top
file to the right number of water molecules and octane molecules, and make an
intelligent guess whether or not too many octane molecules have been thrown
away.
If desired, run grompp
and genion
to add ions. Adjust something.top
accordingly.
Energy minimize the system and look at it
intelligently with rasmol
4. Running the helix system
Position restraints. As usual, creating the
starting structure involves unphysical situations. It might be a good idea to
do a position restraint run first, with the helix position restrained so that
the octane and water molecules can properly rearrange. 100 ps
should be long enough given the fast dynamics of both octane and water. On the
other hand, if in simulations of many nanoseconds the details of how the
starting structure was created really matter in these systems we might as well
give up. This is a difference with phospholipid
simulations, where we currently don't have an alternative for careful creation
of starting structures.
Force field. I suggest using Gromos96 unless
you have a really good excuse. There are a number of standard input parameters.
Have a look at this example grompp.mdp file. Most
interesting: the neighbourlist update frequency is
higher, and a twin range cutoff of 0.8/1.4 nm for both LJ and Coulomb is
standard.
Long range electrostatics. You have a choice
of twin range cutoff, reaction field, generalized reaction field, PPPM or PME
for electrostatics. Reasonable options are twin range cutoff only, or twin
range plus PME. I don't know what is better at this point.
Dummies. An interesting option for rough
simulations like these is to use dummy hydrogens and
aromatics. This allows a 5 fs timestep (possibly
larger, if you feel adventurous) by using a special treatment of bonded
interactions involving hydrogens and aromatic rings.
You need to select the -dummy hydrogens or -dummy
aromatics and -heavyh options in pdb2gmx for making your helix .itp file to use this, and adjust
the timestep and neighbourlist update frequency in
your grompp.mdp
file (0.005 for timestep, 3 for nstlist).
5. Problems
Can't best fit helices because the number of C-alphas is different. Make 2 indexfiles,
one for old.pdb,
one for helix.pdb. Delete enough C-alpha atoms from the
C-alpha index group for the helix with more C-alpha atoms to make the size of
the index groups you are fitting with the same. Then use the -n1 and -n2
options for g_confrms
to fit on a subset of C-alphas. If the length of the old and the new helix is
very different just start from scratch rather than trying to use this fitting
stuff.
System does not minimize, or crashes with a LINCS or SHAKE error. If
the system does not minimize after step 3a, best fitting, or minimizes but
crashes with a LINCS or SHAKE error right at the beginning of a normal run, you
probably have an octane molecule sticking through an aromatic ring or some fun artefact like that. Find out which molecule it is (check
the log file which gives the deviations from normal geometry and look for large
deviations), delete the octane or water molecule from the starting
configuration, adjust the number of octane and water
molecules accordingly in the something.top file, and try minimizing and running
again.
grompp
complains about missing parameters. This still happens with some choices of
termini, notably ACE and NH2. There are updated forcefield files floating
around that include the necessary information but I might have missed updating
the ones you happen to be using. Complain loudly.
The box deforms terribly. Octane does not have an intrinsic area per
molecule like phospholipids. Make sure you have in your .mdp file the
compressibility in x and y set to 0, anisotropic pressure coupling on, and a
reasonable value for the compressibility in z. I just use 4.6E-5, although
that's not likely to be correct, strictly speaking.