[gmx-users] How to use Höltje's cholesterol parameters ?

[Sim gmx]

Hello everyone,

I am currently working on interactions between small biomolecules and
bilayers. Phospholipids parameters come from Peter Tieleman's website
(Berger lipids forcefield) and the small compounds are parametrized for
gromos53a6. I would like to add cholesterol to my bilayers, which is very
frequently done with Höltje cholesterol parameters (see e.g.
http://pubs.acs.org/doi/full/10.1021/ja211929h ). However, these Höltje
parameters were originally designed for the ffgmx forcefield, which is
quite old.

The way we should include these parameters in a Berger lipids - gromos53a6
mixed forcefield has already raised some questions on the mailing list:
http://comments.gmane.org/gmane.science.biology.gromacs.user/68478

However, it remains unclear to me. I guess we can use the bonded parameters
as they are written in the original topology from Höltje, but it becomes
more complicated when talking about the non bonded interactions. Every
atomtypes but two (CB and CR61) from ffgmx also exist in gromos53a6 (and
are thus included in the gromos53a6 forcefield). Hence, I think there are 3
possible ways to make a simulation run without crashing when including this
cholesterol to a berger lipid - gromos53a6 forcefield:

1) Keeping the cholesterol topology file unchanged, and adding the
atomtypes 'CB' and 'CR61' to the forcefield file 'ffnonbonded.itp', with
their parameters coming from the ffgmx forcefield. It means that each
cholesterol will be seen as a "hybrid object", with most of the atomtypes
being gromos53a6 ones, and CB and CR61 being ffgmx atomtypes. Non bonded
interaction involving CB or CR61 atomtypes will be computed with the
standard combination rule.

2) Keeping the Berger lipid - gromos53a6 forcefield unchanged, and changing
the atomtype 'CB' to 'C' and the atomtype 'CR61' to 'CR1' into the
cholesterol topology file, C and CR1 being the corresponding atomtypes
found in gromos53a6. It means that each cholesterol will be seen as a
"gromos53a6" object for the non bonded interactions.

3) Changing every atomtype from the cholesterol topology file to make them
different from gromos53a6 atomtypes (for instance: CH2 becomes CH2F (for
ffgmx)) and adding all these ffgmx specific atomtypes in the forcefield
file 'ffnonbonded.itp' with the proper parameters. It means that each
cholesterol will be seen as a "ffgmx object", each non bonded interaction
being computed with the standard combination rule.

Solution 1 seems to be quite obvious, but it sounds a bit weird to me,
because it mixes up ffgmx and 53a6 interactions. Solution 2 seems more
consistent, but if we use "pure" gromos53a6 non bonded interactions, maybe
we should also "translate" bonded parameters from ffgmx to gromos53a6 ?
Solution 3 is maybe the most logical solution, but it seems that ffgmx is
now considered to be deprecated...

The question becomes even more complicated when considering the fact that
CH2 and CH3 atomtypes could lead to overcondensed bilayers, according to
some authors who advise to switch them to 'LP2' and 'LP3' Berger lipid
atomtypes.

It is possible that those solutions give quite similar results (if ffgmx
and gromos53a6 forcefields are "similar enough"), but I am very curious to
know the 'usual' protocol that is followed when people only write "Höltje
parameters were used".

Another solution would be to use the cholesterol topology file found on ATB
(manual validation), frequently used by Pr. Alan E. Mark, but never used
with Berger lipids phospholipids...

Thank you in advance for your help!