"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.|
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:
- We calculate the dihedral distributions from the CHARMM36 (and MacRog) for the glycerol and headgroup region.
- We set up dihedral potentials for the Berger model such that they resemble the ones in the CHARMM36 (and/or MacRog).
- We use these potentials as a starting point for a Lamberg-style optimization.
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).|
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 SCDs 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.