Some points concerning the POPC full-hydration results

Here is a summary of my current understanding on how the order parameters are reproduced by other force fields than the Berger:

In all-atom force fields CHARMM [Klauda et al. JPCB 114, 7830 (2010)] and GAFF [Dickson et al. Soft Matter 8, 9617 (2012)], the headgroup and glycerol group order parameters are relatively well reproduced. However, in the case of CHARMM it is not clearly reported if order parameter splitting is observed for the alpha and beta carbons.

For the Stockholm lipids [Jämbeck et al. JPCB 116, 3164 (2012)], I have not seen order parameters calculated for headgroup or glycerol. Similarly I have not seen comparison for the united atom models by Kukol [JCTC 5, 615 (2009)] or Chiu et al. [JPCB 113, 2748–2763 (2009)].

Poger et al. [J Comput Chem 31: 1117–1125, (2010)] write about their united atom model: "At 323K they found values of ν of approximately 5, 6, and 28 kHz for the methylenes α, β, and δ, respectively, yeilding the corresponding |SCD| values of 0.05, 0.04, and 0.22. The deuterium order parameters calculated from the simulations are consistent with the experiments with |SCD| values of 0.09 ± 0.01, 0.01 ± 0.01, and 0.16 ± 0.01 for the methylenes at positions α, β, and δ, respectively." They do not mention about splittings. To me it seems that these results are rather similar to those in the manuscript, but my interpretation was that the simulations and experiments are not consistent.

In my opinion, it would be useful to have a united atom model which would reproduce the correct structure for the headgroup and glycerol. It might be that either the Kukol or the Chiu force field already achieves this. This should, obviously, be tested first.

If these two do not work, there are many possible approaches that could be tried to fix the issue:

- Both CHARMM and GAFF developers relate the issue into the dihedral potentials. One way to approach the problem would be to take an all-atom dihedral potential that reproduces the experimental order parameters and transform it into the united-atom form. I do not have a clear idea, however, how this should be done in practise.

- It would be interesting to try to construct possible dihedral angle distributions from the experimental order parameters by some numerical fitting methods,  for example, with the methods presented by Thaning et al. [JPCB 111, 13638 (2007)].

- There are some suggestions in the literature how the headgroup should behave in order to generate the measured order parameters. Quoting from Ferreira et al. [PCCP, 15, 1976 (2013)]: " The first detailed model for the structure of the choline headgroup and glycerol backbone was built based on 2H and 31P NMR results from DPPC bilayers and crystallography studies [Seelig et al. BBA 467, 109 (1977), Pearson et al. Nature 281, 499 (1979)]. This model assumes rapid transitions between two enantiomeric choline states, a free rotation around the Cg1–Cg2–Cg3–O dihedral, and the assumption that the Cg2–Cg3 bond is on average perpendicular to the plane of the bilayer [Seelig et al. BBA 467, 109 (1977)]. Such description captures almost all NMR parameters; however, the last two assumptions are not compatible with the two distinct order parameters for Cg3, and thus other models were proposed e.g. in which the angle of the Cg1–Cg2–Cg3–O dihedral was completely fixed [Hong et al. Biochemistry 35, 8335 (1996)]." One option for us would be to directly include these ideas into the MD model.

- Finally, it might also be useful to remember that based on very recent NMR results, it seems that the headgroup and glycerol dynamics is too slow in the Berger force field [Ferreira PhD thesis].

Samuli Ollila
Aalto University