Tuesday, January 19, 2016

Does the glycerol backbone structure depend on initial structure?

For model by Ulmschneiders it does. This was observed while finalizing publication from NMRlipids II project.

First it was observed that two independent simulations ran with Ulmschneiders model (http://dx.doi.org/10.5281/zenodo.13392 and http://dx.doi.org/10.5281/zenodo.30904) gave different order parameters for glycerol backbone. It turned out that the essential difference between simulations was the initial structure, other one was dowloaded from Lipidbook and other one was generated with CHARMM GUI. More specifically, the difference in O-g3-g2-g1 dihedral in the initial structure. The results from different initial configurations are shown in Fig. 1 and the whole discussion is here.

Figure 1: Order parameter results for POPC headgroup and glycerol backbone from model by Ulmscneiders simulated with two different initial configurations (Lipidbook and CHARMM GUI). The segments labeled 4-8 belong to glycerol backbone and segments 0-3 belong to headgroup.

Importantly, superficial comparison to the experimental results in NMRlipids I publication (Fig. 2 in publication) indicates that the result with CHARMM GUI initial structure gives order parameters closer to experiments. However, in NMRlipids I publication the results with Lipidbook initial structure was used (the other one was not yet ran) and the force field quality (Fig. 4 in publication) was ranked based on these results. In conclusion, using more realistic initial structure, especially for O-g3-g2-g1 dihedral, gives better structure with Ulmschneiders model and increases its ranking (Fig. 4 in NMRlipids I publication).

Now, the question arises if the glycerol backbone (or other) structures depend on initial configurations also in other models? This has been somewhat discussed in the project but not systemically tested, see About glycerol conformations post and comments in Accuracy of order parameter measurements post. My argument has been that based on comparison to the NMR relaxation data, the typical simulation times are sufficient to sample full structural phase space. If this does not happen, the simulation is too rigid. Also, independent simulations have been ran with several models tested during the project (at least with CHARMM, Berger, Slipids, GAFFlipid) and there has never been significant differences in order parameter results between different runs. However, it is not clear if in all of these simulations the initial structures have been really independent.

In conclusion, it seems likely to me that in most models the glycerol backbone structure do not depend on the initial structure, but in Ulmschneiders model it does. Further tests would be needed to confirm this argument. I think that this is relevant topic for users and developers of atomistic resolution MD simulations.

Also another question arises; should we update the results in NMRlipids I publication? At least we could update the version in GitHub. However, science always goes forward and this will not be the only issue to be updated. It may be reasonable to update more changes at once. Traditionally, we could submit a correction to a journal or make a new article. However, due to Internet we have the possibility to update outdated information directly into the article without losing history; but I am not sure how general readers would realize that updated version is available. While waiting for journals publishing modifiable articles, we have to use some temporary approaches. Any ideas?