Ulmschneider results for POPC without ions

[ohsOllila]

I have calculated the order parameters from two different trajectories ran with Ulmscneider model and the results are essentially different.

The results from this trajectory: http://dx.doi.org/10.5281/zenodo.13392 are here: https://github.com/NMRLipids/lipid_ionINTERACTION/blob/master/Data/POPC/NaCl/ULM/0na/OrderParametersMATTI.dat and copypasted here: -0.0676466 0.00295541 -0.0644059 0.00294201 0.0287083 0.00415642 0.0585144 0.00541744 -0.316206 0.00424214 -0.301216 0.00679239 -0.29812 0.00717626 -0.282629 0.00473567 0.18337 0.011728 These are close to what was reported in the NMRlipids I project.

The results from this trajectory: http://dx.doi.org/10.5281/zenodo.30904 are here: https://github.com/NMRLipids/lipid_ionINTERACTION/blob/master/Data/POPC/NaCl/ULM/0na/OrderParameters.dat and copypasted here: -0.0666512 0.00282814 -0.0569859 0.00328485 0.0395524 0.00475013 0.0489784 0.00486469 -0.226971 0.0134899 -0.21175 0.0104782 -0.224464 0.00942215 -0.249064 0.00435139 -0.0408404 0.00670562 These are closer to the experimental values.

The differences in alpha and beta are in reasonable range but difference in glycerol backbone are significant.

I have added these results in the ion manuscript since the alpha and beta are in focus here. However, we should find out the origin in the difference.

Potential reasons I have in mind:

• Mistake in *itp files
• Gromacs version differences. However, in the latter data set the difference between verlet and group scheme was tested and no difference found, see https://github.com/NMRLipids/lipid_ionINTERACTION/tree/master/Data/POPC/NaCl/ULM/0na
• Difference in initial structure, in the first one is from Lipidbook, the second from CHARMM GUI.

[MykhailoGirych]

I've checked all *itp files, they seem to be absolutely identical except that Matti's gbsa.itp has two additional strings at the end: MCH3A 0 0 0 0 0
MCH3B 0 0 0 0 0

There are some differences in *mdp files: Matti: nstlist = 10 fourierspacing = 0.1 pme_order = 6 Tcoupl = berendsen tc-grps = system Pcoupl = berendsen constraints = h-bonds

Mykhailo/Samuli: nstlist = 5 fourierspacing = 0.16 pme_order = 4 tcoupl = Nose-Hoover tc-grps = POPC SOL pcoupl = Parrinello-Rahman constraints = all-bonds refcoord_scaling = com nstcomm = 100 comm-mode = Linear comm-grps = POPC SOL

The starting *gro files seem to be identical (except of the coordinates and the number of molecules).

I have a question to Matti: how long did you equllibrate the system?

If the difference in mdp doesn't affect the order parameters than the most probable reasons are GROMACS version or initial structure. We should check the angles in the starting gro files.

[ohsOllila]

According to NMRlipid I paper Matti's run is equilibrated 50ns and analyzed 50ns.

I would first check the initial structure. One option is also start a run with Mykhailo's files from the same starting structure which Matti used.

[mattijavanainen]

How did you choose the run parameters you use? The original paper uses Berendsen barostat and thermostat. It also says that it uses constraints only in bonds involving hydrogen but DPPC/POPC do not have any such bonds in the UA description.

Our simulation temperatures are also different (310 K vs 298 K) and I don't couple the membrane and solvent separately since this was not mentioned in the original paper. Still, I find it hard to believe that these differences in the itp would lead to substantial differences in the order parameters.

[ohsOllila]

I calculated some dihedral angle distribution from both simulation trajectories. There is significant difference in the O-g3-g2-g1 dihedral:

This is most likely the reason for the order parameter difference. If there is no difference in *itp this probably originates from the initial structures. I did check the distribution in the initial structure by Matti and the distribution was similar. If Mykhailo would push the initial structure I could check that one as well.

There was only small difference in g3-g2-g1-O(sn-1) dihedral, see *xvg files in commit: https://github.com/NMRLipids/lipid_ionINTERACTION/commit/743e440292bb56390b4f1dde293cd78b787979b7

The scripts using the mapping files to calculate dihedrals are *sh files in this commit: https://github.com/NMRLipids/lipid_ionINTERACTION/commit/743e440292bb56390b4f1dde293cd78b787979b7

[ohsOllila]

Mykhailo uploaded the intial structure and I did plot the dihedrals from those as well:

The same difference seen in simulation results is present already in the starting structure. I think that this indicates that also some other results are potentially affected by the initial structure. I think that this needs some more attention.

Just to be sure it may be good if Mykhailo could run a new simulation using the starting structure of Matti.

[ohsOllila]

Mykhailo ran the simulation by using the same initial structure as Matti (from lipidbook) and I calculated the order parameters. Numerical result are here: -0.0706268 0.00305048 -0.0804867 0.00296083 0.034743 0.00576739 0.0663432 0.0059467 -0.318142 0.00706363 -0.301088 0.00860436 -0.313162 0.00812947 -0.293217 0.00506629 0.197927 0.0115895

and plot:

It seem clear that the results depends significantly on initial configuration. Based on NMR relaxation data I have argued in the blog and in NMRlipid papers I and V that realistic model should sample all configurations in these timescales. This would indicate that in realistic model this should not happen. However, there are also other opinions. This issue needs more careful study from many respects. If we would use the results from CHARMM GUI initial structure in the NMRlipids I paper, the ranking of Ulmschneider force field would increase. Further, question arises if the results in other models depends on initial structure as well. For several models we have similar results inpendently reported by different authors (At least CHARMM, Slipids and Berger based models. Different initial configurations are also ran for GAFFlipid due to stereoisomer issue in Lipidbook) and these results have given systemically the same results independetly on author. However, more careful investigation would be needed to check if initial configurations were really independent in the different sources.

I am not sure how much time I will have to invest on this, but people who will do these simulations also in the future should take this issue into account. I will write a blog post about this anyway as soon as I have time.

[ohsOllila]

There is a blog post written about this. It is yet open what should be done with this.

[TomPiggot]

Hi,

Firstly I should apologise for joining in late to all of the great work you have been doing in these projects. I have kept meaning to try and make the time to contribute to the work but one thing after another has kept me too busy to find the time.

As for this issue, I have some old trajectories of DPPC/POPC with some united-atom force fields (the ones in http://pubs.acs.org/doi/abs/10.1021/ct3003157, so Berger and a few other GROMOS based force fields) where I there are simulations that definitely have different initial starting structures of the membranes. Happy to either share the trajectories (2 x 200 ns for each of the two starting structures IIRC) or repeat the above type of analysis of the order parameters. One of these two structures was converted from a CHARMM bilayer and so should provide a good initial starting position for the carbons. I wouldn't be at all surprised if most of the united-atom force field simulations you have will have been started from the same structures, as typically people would just download the bilayers from a few limited locations rather than either build them from scratch or convert structures between force fields (even though both of these are easier nowadays).

Cheers

Tom

[ohsOllila]

Welcome to the project and sorry for the slow reply, I was in conference.

I think that comparing Berger simulations from different initial structures as you describe is very good idea. Currently the most convenient way to do the calculations would be probably by using this script: https://github.com/NMRLipids/lipid_ionINTERACTION/blob/master/scripts/calcORDP_BERGER_0na.sh

If you clone the repo, only thing you need to do (hopefully) is to specify: trajname= tprname402= starttime= endtime= numberOFlipids= outFILE=

If you do not have time to this, you can also share the trajectory and hope that someone else runs this.

It is true that the starting structures may not be fully independent due to the reasons you mention. However, we have quite a lot of data from different sources and I think that there are several sources of initial structures. On the other hand, I am almost sure that some simulations with different models have been started from the structures which have the same origin. Anyway, from the analysis you suggest we learn more.

[jmelcr]

I wouldn't be at all surprised if most of the united-atom force field simulations you have will have been started from the same structures, as typically people would just download the bilayers from a few limited locations rather than either build them from scratch or convert structures between force fields (even though both of these are easier nowadays).

In an ergodic simulation, this should not be a problem, of course.

Is it a sure thing that the discrepancy between the simulations is really caused by a different initial configuration? Couldn't it be that there is a mis-setup of some kind? I'd suggest re-running the same thing with the same software and setup to be really sure... In addition -- if the results really depend on the initial configuration, then the simulation is not ergodic and the drawn statistics are biased in some unknown way.

[ohsOllila]

Please read the above discussion and post http://nmrlipids.blogspot.fi/2016/01/does-glycerol-backbone-structure-depend.html. It seems pretty clear to me that in the case of Ulmschneiders model the difference becomes from the initial structure.

[TomPiggot]

I've finally managed to find and dig out all the old simulation files from archives, etc. I will do the analysis and get back to you.

[TomPiggot]

With Berger, the good news is that the calculated order parameters are the essentially same using two completely different starting structures (I have only done the calculation for POPC):

Converted from a CHARMM bilayer:

0.0376519 0.00314741 0.0478337 0.00315695 0.093418 0.00486981 0.157072 0.00609077 -0.181986 0.00602361 -0.244166 0.00378593 -0.130981 0.00630915 0.141791 0.00981575 0.105315 0.00691997

And with a starting structure converted from a GROMOS43A1-S3 bilayer:

0.0451934 0.00443239 0.0526291 0.00422757 0.0928523 0.00631501 0.157778 0.00878606 -0.18851 0.00716714 -0.250237 0.00537912 -0.138085 0.009007 0.148824 0.0124284 0.0953234 0.00913057

For this latter structure, a range of simulation settings were also used. None of the changes in cut-off change these order parameters (with either longer vdW or with a dispersion correction). I can provide these numbers too, if needed.

Slightly off topic (related more to the original work), but I thought it was worth noting that the only thing that does change these order parameters by any real amount is whether or not the additional dihedrals around the glycerol region, found in the downloadable 'Berger' parameters from the Tieleman website, are included or not. The simulations above have these 'extra' dihedrals included (and by extra I mean in terms of correctly applying the original GROMOS based bonded parameters, as was reported to be done by Berger et al, with just one dihedral applied across the same two central atoms). The results below show the differences when you remove these extra dihedrals (results provided for three different cut-off schemes to highlight that these differences are consistent)

0.0480216 0.00311658 0.0549797 0.00338997 0.112526 0.00526322 0.182798 0.00630141 -0.271834 0.00362744 -0.192003 0.00550344 -0.195213 0.00547021 0.00865883 0.0090997 0.24195 0.00666299

0.0549084 0.00264333 0.0623501 0.00287378 0.119894 0.00514792 0.192775 0.00515024 -0.279034 0.00299184 -0.193927 0.00485098 -0.215712 0.00406648 0.00765309 0.00764178 0.258365 0.00561748

0.0469954 0.00383024 0.0577935 0.00415169 0.120138 0.00643224 0.188602 0.00639054 -0.275375 0.00405178 -0.190575 0.00591569 -0.198502 0.00571969 -0.0142326 0.008787 0.241964 0.0068721

And (apologies for all the numbers) a plot to highlight these differences regarding the dihedrals:

Finally, I'm also checking the simulations I have done for the POPC with Kukol/GROMOS-CKP parameters as I also have these to hand and there are simulations with both different starting structures and there are simulations with and without extraneous dihedrals in the glycerol region (which are found in the original Kukol parameters). Removing these and also changing the dihedrals around the double-bond which results in bad order parameters gives the GROMOS-CKP POPC lipid parameters. I imagine the CKP order parameters should look more like the Poger ones than the Kukol ones you reported in your original work, but I'll let you know once I've done this.