Towards a new version of the manuscript

In the opening post of the nmrlipids project, published simultaneously with the accompanying manuscript, it was written that:  

"The ultimate goal of this blog is to find an atomistic (preferably united-atom) force field that reproduces the experimental properties discussed in the manuscript. Naturally the optimal situation would be that some of the already available force fields would fulfill this goal. If this, however, turns out not to be the case, the goal will be to find the appropriate modifications."

As discussed in The lipid forcefield comparison post, none of the tested united-atom force fields fulfills the goal\(^1\). The same post concluded that three all-atom force fields (CHARMM36, GAFFlipids and MacRog) gave a reasonable agreement with experiments for a fully hydrated pure bilayer, when the absolute values of headgroup and glycerol order parameters were compared. Since then, the comparison has been refined to include the sign of the order parameters as well as the R/S labeling of hydrogens. Fig. 1 shows an up to date comparison of the top three all-atom force fields, the Berger force field, and the Lamberg parameters (currently in development in this blog) to experiments.

Figure 1. A set of simulation results for the POPC headgroup and glycerol order parameters reported in the blog is compared to experiments.

With these refinements (sign and R/S labeling) it becomes clear that the order parameters for the g\(_1\) group in the GAFFlipid force field are not in a reasonable  agreement with experiments. Thus there are only two all-atom models left to consider: CHARMM36 and MacRog. Interestingly, as discussed in the post Response of headgroup and glycerol order parameters to changing conditions: Results, the responses of these two force fields to changing conditions (dehydration, ion concentration and cholesterol content) are in qualitative agreement with experiments. However, the Na\(^+\) association with the bilayer is clearly too strong for MacRog; the association is less for CHARMM36, but more simulations are needed to complete the comparison to experiments.

In conclusion, CHARMM36 is a force field that relatively well satisfies the original goal of the nmrlipids project. In this respect, we could now just write up the manuscript, submit and finish the project.

Having said that, I still think that it would be very useful to have also a united-atom model that fulfills the original goal. Regarding the studies of pure lipid bilayers, CHARMM36 appears significantly slower than united-atom models and also slower than other all-atom models. Based on my personal experience (I have not investigated it further) I think the slowness is due to the Lennard-Jones interactions of the water hydrogens in the CHARMM water model\(^2\). In addition to speed, there will probably be interest in using a lipid model with correct molecular structure in combination with proteins (or other molecules) that lack a model combatible with the CHARMM36 lipids.

If we manage to get correct structures for one united atom model, the same procedure could probably be used for other models as well. A promising method to tune the potentials to get the correct order parameters was already reported in this blog by Antti Lamberg in March. However, at the time, we had not considered the signs and the R/S labeling and thus the order parameters were matched to wrong values.

On the basis of these perpectives, I suggest that we make one further attempt to correct the order parameters for the Berger model in the following way:

1. We calculate the dihedral distributions from the CHARMM36 (and MacRog) for the glycerol and headgroup region.
2. We set up dihedral potentials for the Berger model such that they resemble the ones in the CHARMM36 (and/or MacRog).
3. We use these potentials as a starting point for a Lamberg-style optimization.

I think that this procedure has potential to produce a united-atom model that will have the correct glycerol and headgroup structures. If this procedure does not succeed, I suggest we finish this project and publish the results; further attempts to fix united-atom model can be continued in another project in any case. I will already start to prepare the next version of the manuscript.

First results for step 1 of the procedure. As preliminary work for the above procedure I plotted 40 consecutive lipid structures observed in Berger and CHARMM36 simulations such that the glycerol groups were aligned. The results are shown in Figs. 2 and 3.

Figure 2. 40 consecutive (\(\Delta\)t=1.25 ns) structures observed in a Berger simulations (POPC, T=298K). Hydrogens have been added to the visualization by the Gromacs protonate program in the same way as for the order parameter calculation.

Figure 3. 40 consecutive (\(\Delta\)t=1.25 ns) structures observed in a CHARMM36 simulation (POPC, T=303K).

Comparing the Figs. 2 and 3 one can clearly see that the conformations sampled by the g\(_1\) and g\(_3\) groups significantly differ between the models - which is no surprise as the models predict different order parameters. In all configurations of the CHARMM36 model the hydrogens of the g\(_1\) group are pointing to the same direction, indicating a rather rigid dihedral in the g\(_1\)-g\(_2\) bond. In the Berger model, the g\(_1\)-g\(_2\) bond seems to be significantly more flexible. The situation is opposite for the g\(_3\) group: the hydrogens are pointing in the same direction in almost all configurations for the Berger model, while CHARMM36 seems to be significantly more flexible. This difference seems to have significant effect to the conformations of the whole PC headgroup. (When looking at the figures, it should be noted that the molecules are shown from different angles to make the above points more clear.)

I am currently processing some ideas how to proceed to the step 2 in the procedure. I will report the progress in the near future. It would be also interesting to see structural figures with aligned glycerol from the MacRog model.




\(^1\)The recent united atom CHARMM has not been tested yet. The authors of the model state that "The headgroup is still an AA potential in C36-UA, so only the chain S_CDs were compared." Thus similar results than in CHARMM36 are expected. However, the united atom CHARMM might still be quite slow due to the water model. Some results with the model would be interesting.

\(^2\)Using regular TIP3P model does not seem to be a solution since significant changes in phase behaviour was reported by Piggot et al.