[gmx-users] CHARMM36 - Smaller Area per lipid for POPE - Why?

[Sebastien Cote]

Dear Gromacs users,
 
My simulations on a POPE membrane using the CHARMM36 parameters are giving ''area per lipid'' values well below the experimental value (59.75-60.75 Angstroms2). Is their someone else experiencing a similar problem? If yes, how did you solved it? 

I did the following :

I used the CHARMM36 parameters kindly provided by Thomas J. Piggot on the Users contribution section on Gromacs website.
My starting configuration was taken from : http://terpconnect.umd.edu/~jbklauda/research/download.html
It is a POPE membrane of 80 lipids equilibrated in NPT at T=310K and P=1atm for 40 ns. It is taken from the article Klauda, J. B. et al. 2010 J. Phys. Chem. B, 114, 7830-7843.

At first, I tested normal TIP3P vs. CHARMM TIP3P and saw that normal TIP3P gives smaller Area per lipid of about 2-3 Angstroms. This was also observed by T.J. Piggot (personnal communication) and Tieleman (Sapay, N. et al. 2010 J. Comp. Chem. 32, 1400-1410). So, I will present only the simulations using CHARMM TIP3P. As in Klauda's paper, my simulations are at 310K and 1 atm. As them, I used a switch cutoff for vdw, and I used normal cutoff for PME. The simulations are 20 ns. I can send my .mdp file for more details. I varied the switch condition on vdw : 
 
1- For a switch from 0.8 to 1.2 (as in Klauda's paper), I got Area per lipid of about 56.5 Angstroms2; whereas they got 59.2 in their paper, matching the experimental value of 59.75-60.75. 
2- For a switch from 1.0 to 1.2, I got Area per lipid of about 53.5 Angstroms2, which is smaller than the previous cutoff. This is surprising since a previous thread on gromacs-users mailing lists said that increasing the lower cutoff, increased the Area per lipid or had not impact on POPC of DPPC. 
3- For a switch from 1.1 to 1.2, I got Area per lipid of about 55 Angstroms2. 
4- For a hard cutoff at 1.4, I got Area per lipid of about 52 Angstroms2.  

I also tried to re-equilibrate the membrane in the NPAT ensemble for 10 ns at 310K and 1 atm. Then, when I launched the simulation in NPT, I ended up with different results :

1- Switch from 0.8 to 1.2 gave a smaller area per lipid of 54 Angstroms2.
2- Switch from 1.0 to 1.2 gave a larger area per lipid of 55 Angstroms2. 
4- Hard cutoff at 1.4 gave a similar area per lipid of 52.5 Angstroms2.

I looked at the POPE paramaters for CHARMM36 in Gromacs, and they agree with the published parameters.

Am I doing anything wrong? Is their someone else experiencing a similar problem for POPE? If yes, how did you solved it?

Should I instead use CHARMM27 parameters in the NPAT ensemble? I want to study the interaction between a peptide and the POPE membrane. I am troubled that the NPAT ensemble might influence my results in a bad way. Also, I can not use OPLS AA nor GROMOS for the protein interactions because these force fields are not giving the correct structural ensemble for my peptide in solution. 

I am willing to send more information if you need. 

Thanks a lot, 
Sincerely,

Sébastien