Step Two: Modify the Topology
The lipid bilayer we will be simulating is DPPC (dipalmitoyl phosphatidylcholine), 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?
There is no universal force field for simulating all proteins, nucleic acids, lipids, carbohydrates, and arbitrary small molecules. Since we are simulating a system that contains both protein and lipids, how do we solve the force field problem? Under GROMACS, the most widely-used parameters for the lipid component of membrane-protein simulations are commonly called "Berger 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 gromos53a6.ff/forcefield.itp. Make a new directory in your working directory called "gromos53a6_lipid.ff" and copy the following files from gromos53a6.ff into it:
aminoacids.rtp aminoacids.hdb aminoacids.c.tdb aminoacids.n.tdb aminoacids.r2b aminoacids.vsd ff_dum.itp ffnonbonded.itp ffbonded.itp forcefield.itp ions.itp spc.itp watermodels.dat
Next, create a
GROMOS96 53A6 force field, extended to include Berger lipid parameters
Next, 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
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
It is also possible to use the OPLS-AA force field to represent your protein. 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:
For a bit more theory on force field combination rules and theory, consult this document. Thanks to Chris Neale for providing this link.
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.
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.
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