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?
The 2015 paper by Ferreira et al., I think, showed that PC lipids sample the relevant parts of their conformational ensemble in 50 ns. And, therefore, a correct MD simulation starting from a typical conformational distribution of PC lipids should also sample the relevant conformational ensemble in 50 ns, and thus produce the correct C-H order parameters.
ReplyDeleteHowever, if the simulation is started from a pathological conformational distribution, I think we are looking at a process that is (or at least could be) very different from the sampling process studied in the NMR experiments of Ferreira et al.
That is, when we start from the pathological case, we have an equilibration process from an unlikely part of the conformational space to the typically populated part. This process could (at least in principle) have considerable different time scale than the (50 ns) equilibrium sampling process.
Now, it might be of interest to test if (and how fast) the better force fields (such as C36) are able to find the realistic conformational ensemble, when started from a pathological state. However, as we do not know the time scale expected for such a process in real life, these studies would be mostly of interest in terms of knowing how sensible the conformational ensembles observed in MD simulations are on the initial conformations. (But they would not tell us if this returning process happens on a physically correct time scale.)
However, what seems to clearly be the case is that we should test the various models starting from initial conformations that are as close to the real life conformations as possible. If the models are able to keep the lipids in this part of the conformational space and sample it properly (give correct OPs), then they should be judged as good force fields.
Therefore, I would suggest running simulations of all the badly ranking force fields starting from the final conformations of the C36 simulations.
I am no sure anyone would like to include a "systematic study of conformational sampling of lipids" into the current publication.
ReplyDeleteMy personal opinion on this is that it would be the best if we use good initial conditions (meaning close to the observed mean conformation) and add a note of awareness to the discussion -- that we see for example Ulmschneiders model not to equilibrate within xx ns into the proper average conformation and hence the model might easily provide non-ergodic systems.
I wouldn't do much about the previous paper until we have at least some numbers to add and make the argument stronger.
Can you check, what the area per lipid is and how it fluctuates in either of the two "inconsistent" simulations? This would be an easy source of the observed discrepancy.
ReplyDeleteIf interested, you can do this yourself. The links to datasets in Zenodo are in the beginning of the post:
Deletehttp://dx.doi.org/10.5281/zenodo.13392 http://dx.doi.org/10.5281/zenodo.30904
Anyway, by looking at the differences in dihedral angle distributions reported in https://github.com/NMRLipids/lipid_ionINTERACTION/issues/8,
it seems quite clear to me that the difference comes from the dihedral.