Step Two: Modify the Topology
The lipid bilayer we will be simulating is DPPC (dipalmitoylphosphatidylcholine), so we need parameters for it as well. But if pdb2gmx is designed to handle only proteins, nucleic acids, and a finite amount of cofactors, where do we get parameters for this molecule, and how do we apply them to our system?
For this tutorial, we will use a united-atom force field to describe the lipids derived by Berger, Edholm, and Jähnig (ref). These parameters can be combined with a GROMOS representation of the protein.
The so-called "Berger lipids" are somewhat of a hybrid between GROMOS atomtypes and OPLS partial charges. Since the long alkane chains are poorly represented by GROMOS bonded parameters, a Ryckaert-Bellemans dihedral potential is used, and a scaling factor of 0.125 is applied to Lennard-Jones 1-4 interactions. These lipid parameters are distributed by D. Peter Tieleman, through his website. While you are there, download the following files:
So what exactly is lipid.itp, and how do you use it? Think of this analogy: gromos53a6.ff/forcefield.itp is to your protein what lipid.itp is to the lipids. Essentially, lipid.itp contains all the atom types, nonbonded parameters, and bonded parameters for a large class of lipids, much like how forcefield.itp, ffnonbonded.itp, and ffbonded.itp function in relation to a protein. That said, we cannot simply
To use the parameters in lipid.itp, we will have to make some changes to our pre-packaged GROMOS96 53A6 force field files (in $GMXLIB/gromos53a6.ff). Make a copy of this directory in your working directory called "gromos53a6_lipid.ff" (assuming you have GROMACS installed in /usr/local/gromacs):
cp -r /usr/local/gromacs/share/gromacs/top/gromos53a6.ff/ ./gromos53a6_lipid.ff
Inside this directory, you will find the following files:
aminoacids.c.tdb aminoacids.hdb aminoacids.n.tdb aminoacids.r2b aminoacids.rtp aminoacids.vsd atomtypes.atp ff_dum.itp ffbonded.itp ffnonbonded.itp forcefield.doc forcefield.itp ions.itp spc.itp spce.itp tip3p.itp tip4p.itp watermodels.dat
GROMOS96 53A6 force field, extended to include Berger lipid parameters
The files contained in the gromos53a6_lipid.ff directory constitute the complete description of the force field. The files serve the following purposes:
Now, to add the lipid parameters into the parent force field, we will need to copy and paste the entries in the
LO 8 15.9994 0.000 A 2.36400e-03 1.59000e-06 ;carbonyl O, OPLS LOM 8 15.9994 0.000 A 2.36400e-03 1.59000e-06 ;carboxyl O, OPLS LNL 7 14.0067 0.000 A 3.35300e-03 3.95100e-06 ;Nitrogen, OPLS LC 6 12.0110 0.000 A 4.88800e-03 1.35900e-05 ;Carbonyl C, OPLS LH1 6 13.0190 0.000 A 4.03100e-03 1.21400e-05 ;CH1, OPLS LH2 6 14.0270 0.000 A 7.00200e-03 2.48300e-05 ;CH2, OPLS LP 15 30.9738 0.000 A 9.16000e-03 2.50700e-05 ;phosphor, OPLS LOS 8 15.9994 0.000 A 2.56300e-03 1.86800e-06 ;ester oxygen, OPLS LP2 6 14.0270 0.000 A 5.87400e-03 2.26500e-05 ;RB CH2, Bergers LJ LP3 6 15.0350 0.000 A 8.77700e-03 3.38500e-05 ;RB CH3, Bergers LJ LC3 6 15.0350 0.000 A 9.35700e-03 3.60900e-05 ;CH3, OPLS LC2 6 14.0270 0.000 A 5.94700e-03 1.79000e-05 ;CH2, OPLS
Append the contents of the
To recap, you must make all of the following changes for the force field to be functional:
Finally, change the
Lastly, we need to include the specific parameters for our DPPC molecules. The process for doing so is quite simple. Just add the line
; Include Position restraint file #ifdef POSRES #include "posre.itp" #endif ; Include DPPC chain topology #include "dppc.itp" ; Include water topology #include "gromos53a6_lipid.ff/spc.itp"
Notes on Force Fields
Since the Berger parameters were drawn from aspects of both the GROMOS and OPLS-UA force fields, it is also possible to use the OPLS-AA force field to represent your protein and have a compatible model. You will have to make some modifications to lipid.itp so that it is consistent with OPLS conventions. For information on how to do this, read the following links:
Regardless of which setup you choose, you must demonstrate that your model is reasonable. The particular parameter combinations described here have been shown to work in-house and according to reports by other users. It is up to you to convince your audience (i.e., reviewers) that you know what you are doing and that your model is valid. For instance, if you intend to do long simulations, it is also known that GROMOS96 53A6 under-stabilizes α-helices and you may end up with a random coil rather than a correct helix. In this case, you may want to use the newer GROMOS96 54A7 parameter set, or another force field entirely.
If you have correctly followed all of the steps above, you will have a fully functional force field that can be used to process other membrane proteins with pdb2gmx. Doing so removes the need to manually hack topol.top after the fact. Placing the new gromos53a6_lipid.ff directory in $GMXLIB will allow you to use this force field system-wide.
Site design and content copyright Justin Lemkul
Problems with the site? Send them to the Webmaster