PN vector orientation not a good measure for evaluating phospholipid force field performance, use head group order parameters instead

The article below is a "guest-author" contribution by O. H. Samuli Ollila in reaction to a recent post on this blog:

"A property commonly calculated from molecular dynamics simulations of phosphatidycholine bilayers is the angle between the membrane normal and the vector from headgroup phosphorus to nitrogen (P-N angle). Sometimes the average angle is also compared to experimental estimates.

The most extensive comparison of this kind that I have seen is in the Martini group blog, where a color coded matrix visualizes the deviation between experimental and simulation values for several key properties of lipid bilayers. The simulation values in the blog were taken from extensive study by Piggot et al. (J. Chem. Theory Comput. 8, 4593 (2012)) and experimental values from various sources.

One of the shown properties is the P-N angle with the experimental value taken from Akutsu and Nakamori (Biochemistry 30, 4510 (1991)). When comparing this particular property between experiments and simulations, however, it is good to keep in mind how the experimental value is reached. The starting point of the Akutsu-Nakamori paper is that they have very accurate NMR spectroscopy and laser Raman measurements on the headgroup structures present in a phosphatidylcholine lipid bilayer.

More specifically, they know from NMR the quadrupole splittings of the deuterons in the choline group (i.e., the order parameters for the gamma, beta and alpha groups, see discussion at nmrlipids.blogspot.fi), phosphorus chemical shift anisotropy in the phosphate group, and from Raman experiments that the O-C-C-N backbone dihedral is in the gaughe conformation. In a nutshell, their approach is to construct different structural and dynamical models and then to check which ones reproduce the experimental observations. They conclude that in the most reasonable model the P-N angle is 72 degrees. Note that a similar approach, but with different experimental observables, has been used by Semchyschyn and Macdonald (Magn. Res. Chem. 42, 89 (2004)), and by Buldt and Wohlgemuth (J. Membr. Biol. 58, 81 (1981)). Semchyschyn and Macdonald ended up suggesting a model with a 60 degree P-N angle. Buldt and Wohlgemuth do not suggest a value for P-N angle but suggest several other conformational details characteristic for different headgroups.

Now how do MD simulations relate to this? Just as the models constructed by Akutsu and Nakamori, Semchyschyn and Macdonald or Buldt and Wohlgemuth, in the lipid bilayer MD simulations we have an atomistic resolution model. This model is now constructed by combining the QM calculation and various empirical data (details vary between the force fields but with the same aim that each lipid samples the correct conformations). And just as before, a model can be considered reasonable if it reproduces the experimental data used by Akutsu and Nakamori, Semchyshcyn and Macdonald, and Buld and Wohlgemuth. Once a model reproduces these data, then the P-N angle can be calculated from the model, as done by Akutsu and Nakamori or Semchyshcyn and Macdonald from their models.

In my opinion, the first check for the reasonability of an MD model would be to check if it reproduces the order parameters for the hydrocarbon segments in the headgroup and glycerol region. This kind of testing of MD force fields is in progress at nmrlipids.blogspot.fi. If some force field passes this test it should next be tested against the other experimental observables that can be directly measured. These include at least the phosphorus chemical shift anisotropy and the dipolar splittings between phosphorus and headgroup carbons. These observables are calculated from simulations and compared to experiments for example by Chowdhary et al. (J. Phys. Chem. B 117, 9142 (2013)), and Prakash and Sankararamakhrishnan (J. Comput. Chem. 31, 266 (2010)).

This line of thinking would have influence also on the color coded matrix in the MARTINI group blog: Maybe the P-N angle should be replaced by a measure that describes the deviation between calculated and measured order parameters? This is, of course, somewhat complicated since there is are large number of order parameters. One option could be to try and define a reasonable norm that would describe the average deviation over the hydrocarbon groups; alternatively one could simply put all the carbons into the color coded matrix. I think that this could change some conclusion drawn from the color coded matrix: For example, the P-N vector slot for the GAFF force field is red, even though it is the best published force field in terms of reproducing the headgroup and glycerol NMR order parameters (see nmrlipids.blogspot.fi). Above all, comparing the order parameters between simulations and experiments would have one significant advantage over the properties now used in the MARTINI group blog: The data are reliable and not sparse, and they do not depend on the technique used and the lab in which they were measured.

In conclusion, I think that to evaluate force field quality in the headgroup and glycerol regions one should primarily compare the measured and calculated NMR observables such as the order parameters. If these values match, then the P-N dipole angle calculated from simulations can be considered to be an "experimental value" just as much as the values estimated by Akutsu and Nakamori or Semchyschyn and Macdonald."