Large difference in PC alpha order parameter between Lipid14 POPC and POPC:POPS (5:1) simulations

[ohsOllila]

Recent commit https://github.com/NMRLipids/NMRlipidsIVotherHGs/commit/141a1b4fbbae8d657ea87590de864ed02c58c0a2 by J. Melcr gives order parameters for POPC at 298K: POPC_a1 error POPC_a2 error POPC_b1 error POPC_b2 error 0.076372 0.002768 0.077614 0.002430 -0.000137 0.001962 0.002520 0.002198

These values are close to the ones reported in the NMRlipids I for Lipid14 POPC: https://github.com/NMRLipids/nmrlipids.blogspot.fi/blob/master/DATAreportediINblog/POPC/AMBERLIPID14-303K_blogged-26-03-14.dat

However, when I was comparing these to the POPC:POPS (5:1) simulations from Lipid14/17 with ff99 ions, I noticed that there is a difference of about 0.1 in alpha carbon order parameters between pure POPC and the mixture. Results from POPC:POPS (5:1) mixture with sodium counterions (https://github.com/NMRLipids/MATCH/blob/master/Data/Lipid_Bilayers/POPS%2B83%25popc/T298K/MODEL_LIPID17/OrdParsPOPC.dat):

POPC_a1 error POPC_a2 error POPC_b1 error POPC_b2 error
0.160485 0.005235 0.17849 0.00437 -0.00214 0.00293 0.0018 0.002815

In conclusion, these results would suggest that the addition of 17% POPS would cause a major change in alpha carbon order parameter when ff99 ions are used, but not when Dang ions are used.

I can see two possible explanations. 1) The result is real and for some reason PC headgroup has very different order parameter in the mixture with ff99 ions, which is not seen with Dang ions: https://github.com/NMRLipids/NMRlipidsIVotherHGs/blob/141a1b4fbbae8d657ea87590de864ed02c58c0a2/Data/POPCvsPOPSlipid17.dat 2) Mixtures with ff99 ions are ran with Amber program, but pure system and mixtures with dang ions are ran with Gromacs. It could be that the programs give different results or our analysis code fails somehow for Amber trajectories (I tried to look at this but I did not found any obvious reason). Pure POPC system ran with Amber could be helpful in resolving this issue.

[markussmiettinen]

Pure POPC system ran with Amber could be helpful in resolving this issue.

@hsantila has pure POPC ran with Amber, I believe.

[batukav]

I have re-run the order parameters on ff99 trajectories to check if I made a mistake, but no luck. All the order parameters for ff99 trajectories are OK.

I have POPC at 303K. I also started to run them at 298K. This week I will be able to calculate the order parameters again.

A quick update: I calculated the OP for pure POPC at 303K (run with Amber engine, Lipid14 forcefield). Here are the results

alpha2 0.0766307166357 0.0214733725622 0.00190545369222 alpha1 0.073706762966 0.0244717412655 0.00217151589087

which are similar to the original values posted. Seems like there's something else going on in POPC+POPS mixtures with ff99 ions.

I will re-check all the input structures and parameters to see if there's any error. If any of you can suggest further tests, I'll be happy to run them.

I will update this comment when I finish my checks.

[hsantila]

Hi all,

I do have amber trajs but only using c36, ECC, and slipids force fields. I guess those wont be helpfull in resolving this.

I can say that I reproduced the same POPC order parameters using amber and gromacs for c36.

[jmelcr]

Hi! There is an emerging possibility of running the very same PC:PS 5:1 mix with both Dang ions (done) and also ff99 ions - those by Johan Aqvist (to be done if desired).

This is very technical and I have a feeling that there will be something primitive in the end behind the scenes, but it has to be solved -- small step by a small step, one change at a time.

Let me know if you deem this helpful. Josef

[ohsOllila]

Indeed, these results support the explanation 1. above. Also, the density profiles from simulations with ff99 ions look somewhat strange and the ion binding affinity is low: https://github.com/NMRLipids/NMRlipidsIVotherHGs/blob/master/Figs/CIdensPSOCmixt-eps-converted-to.pdf and http://nmrlipids.blogspot.com/2018/09/nmrlipids-iv-challenges-in-evaluating.html?showComment=1538508682203#c6772152190626178698

Independent simulation with ff99 ions suggested by @jmelcr would be a quite straightforward test to exclude simple mistakes. So if you can do this with little effort, at least simulations with only counterions and one larger concentration would be useful.

[jmelcr]

The simulations are in progress now (only counterions, and added 1M NaCl).

To give a flavour of the possible artifacts with ff99 ions, see the snapshot from the above mentioned ongoing simulations - the overestimated correlations of Na+ and Cl- motions are highlighted with bonds (PS pink, PC cyan, Na+ green, Cl- orange). This will be visible for example in radial distribution functions of the ions, which could be a possible check for @batukav - what do you think?

[jmelcr]

The simulations (Gromacs 2018 was used) with the ff99 ions (Na+ counterions, and also NaCl 1M added solution) for PC:PS 5:1 bilayer at 298K have reached lengths over 400 ns, which is sufficient to reach some conclustions.

The order parameters of POPC are...

• simulation with only counterions:

  POPC_a1 error POPC_a2 error POPC_b1 error POPC_b2 error
  0.084512  0.002483   0.083016  0.002072    0.007816  0.001747   0.003882  0.001989     

• simulation with added 1M NaCl concentration (buffer concentration)

  POPC_a1 error POPC_a2 error POPC_b1 error POPC_b2 error
  0.067843  0.002318   0.070836  0.001840    -0.000353  0.001725  -0.003821  0.001675     

In contrast to the Amber-generated data in this Fig (yellow), these numbers suggest a decrease of the order parameters, as was seen in NMRlipids2 paper in this Fig with pure POPC but for other force fields. The magnitude of the decrease with ff99 ions is weaker than the decrease observed with Dang ions, and it is similar to the decrease observed with other force field combinations (again see this fig )

This is obviously very troublesome - the new data do not support the rationale no.1, yet I do not see a reason for such a dramatic difference (delta~0.1!) between the MD engines as noted in no.2.

It may be, however, that the analysis code for the order parameters, which uses external libraries for reading-in the data, may be the source of the discrepancy - but I do not see, how this could happen either, e.g. there are internal checks for bond lengths, which shall catch problems with PBC or OP-definitions....

Any thoughts on possible Amber-related specifics, which may have been overlooked?

[hsantila]

Hi,

see below:)

to 11. lokak. 2018 klo 12.15 Josef Melcr (notifications@github.com) kirjoitti:

It may be, however, that the analysis code for the order parameters, which uses external libraries for reading-in the data, may be the source of the discrepancy - but I do not see, how this could happen either, e.g. there are internal checks for bond lengths, which shall catch problems with PBC or OP-definitions....

This could be checked by converting amber data to gromacs equivalents. I can provide the code if you are interested.

Any thoughts on possible Amber-related specifics, which may have been overlooked?

Which barostat are you using?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/NMRLipids/NMRlipidsIVotherHGs/issues/21#issuecomment-428901007, or mute the thread https://github.com/notifications/unsubscribe-auth/AMFkilXsg2SOZwTMb36g9j7_GnS-cTehks5ujxpFgaJpZM4W-1xL .

[jmelcr]

Hi =]

Which barostat are you using?

Parinello-Rahman. This and other details are in the mdp file. Unless the membrane undergoes a major transition with the switch of a barostat, this shall otherwise not yield such dramatic changes.

This could be checked by converting amber data to gromacs equivalents. I can provide the code if you are interested.

Interesting, I tried to do just that in the past, but ran into problems with "magic numbers" - not kidding! =] Is the code available somewhere? I think I still have the Amber trajectories downloaded somewhere, so I could just try this option straight away...

[ohsOllila]

@jmelcr could you share also the ion density profiles from these simulations? Then we can see if those are similar to the results from Amber.

[jmelcr]

The density profiles from both simulations (follow the link to see) are different from what is reported in the comparison figure in the manuscript for Lipid17, however, they are much more similar to what can be judged as "common" from the other force fields.

Can this suggest that there is some problem with the trajectory generated by Amber? To give an example, when I open the trajectory n VMD, I obtain a the following picture (blue box represents the simulation box dimensions): Does an output from Amber look like this in common? Is it a correct representation of a simulation box in Amber?

[ohsOllila]

I was also looking the trajectories with vmd previously and was a little bit puzzled. However, I managed to move the atoms inside the periodic box and the trajectory seemed reasonable after this. However, the density profiles indeed seem more reasonable from the Gromacs trajectories indicating that there may be something weird going on with the Amber data and that it would be good to take a little bit closer look to the trajectories.

@jmelcr the density profiles in data format would be useful as well.

[jmelcr]

See the new pull request #44

[hsantila]

It seems like my amber conversion code might not be terribly useful after all.

I've used it to create amber inputs from gromacs topology and .gro, and then convert the amber trajectory back using mdconvert. Mdconvert is used instead of amber cpptraj because I ran into problems with preriodicity not being interpreted correctly with cpptraj.

You can take a look here: https://github.com/hsantila/ambertools

I vaguely remember @Batuhan Kav batuhan.kav@mpikg.mpg.de had some experience with testing different barostats in amber and seeing deviating results?

pe 12. lokak. 2018 klo 10.11 Josef Melcr (notifications@github.com) kirjoitti:

See the new pull request #44

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/NMRLipids/NMRlipidsIVotherHGs/issues/21#issuecomment-429242998, or mute the thread https://github.com/notifications/unsubscribe-auth/AMFkipTWOHUXdIdAa_1TKcvSiBp7JpFAks5ukE63gaJpZM4W-1xL .

[ohsOllila]

@jmelcr thanks for the densities. I would also like to know the exact number of ions and water molecules in the 1M NaCl simulation to double check the concentration here and other simulations.

[jmelcr]

• 3684 water molecules,
• 39 Na added ions, i.e., 51 total (with counterions).
• buffer concentration of added NaCl = 586.8 mM (apparently I somehow miscounted before, sorry!)

[batukav]

Hello,

1)

I vaguely remember @batuhan Kav batuhan.kav@mpikg.mpg.de had some experience with testing >different barostats in amber and seeing deviating results?

Indeed I observed different behavior with the membranes by simply changing the barostat. For all my simulations I am using Berendsen barostat (other option is Monte Carlo barostat which does not give satisfactory outputs). Plus, in every simulation I apply semi-isotropic pressure coupling and set the surface tension to be exactly zero. This combination of barostat-pressure coupling was necessary when I started to simulate membranes two years ago (with Amber16 engine), and so far I think this is the best pair to use with Amber.

Still, change in barostat does not explain what's happening. My simulations with pure POPC give the correct OP. Results start to get weird when the ions are included.

2)

Can this suggest that there is some problem with the trajectory generated by Amber? To give an example, when I open the trajectory n VMD, I obtain a the following picture (blue box >represents the simulation box dimensions): Does an output from Amber look like this in common? Is it a correct representation of a simulation box in Amber?

I don't see any problem here, as this is common. Amber has different options for wrapping the trajectory and depending on how you wrap the atoms the box calculated with VMD can look weird.

3)

To give a flavour of the possible artifacts with ff99 ions, see the snapshot from the above mentioned >ongoing simulations - the overestimated correlations of Na+ and Cl- motions are highlighted with bonds >(PS pink, PC cyan, Na+ green, Cl- orange). This will be visible for example in radial distribution >functions of the ions, which could be a possible check for @batukav - what do you think?

Please look at the snapshot I obtained from PS:PC mixture with 1000mM NaCl using ff99 ions. This is the last frame of the simulation, I don't think radial distribution functions will show any better proof than there's something wrong. With 500mM NaCl, though, I don't observe any clustering.

I am looking deeper into this problem now. It is ridiculous to me that the ions cluster like this. I hope that this is a particular issue with my simulations but not a generic one. This also explains the weir behior we observed with the ion density profiles. I'll reply again when I make some progress.

[ohsOllila]

Thanks, now the density profiles are better in line with the results from other force fields. I have now updated these results in the figures: https://github.com/NMRLipids/NMRlipidsIVotherHGs/blob/master/Figs/CIdensPSOCmixt-eps-converted-to.pdf https://github.com/NMRLipids/NMRlipidsIVotherHGs/blob/master/Figs/CHANGESwithMONVALENTwithPS-eps-converted-to.pdf https://github.com/NMRLipids/NMRlipidsIVotherHGs/blob/master/Figs/HGorderparametersPCvsPS-eps-converted-to.pdf

The results for PS lipids in these figures are still from the Amber simulations.

I am now trying to figure out if we can determine from these which model has the most reasonable counterion binding affinity, but this seems to be very difficult. I will write this discussion somewhere else.

[batukav]

As a follow up for my previous comment:

1) clustering of NaCl is visible in 1000mM, 2000mM, 3000mM, and 4000mM simulations. In all of these simualations, the initial structure does not have such clusters.

Clustering of Åqvist ions is, as far as I know, is a known issue. What puzzles me is that such clustering is not visible in the same simulations using Gromacs, am I right @jmelcr ?

2) about the script for OP calculations: I've checked the results again with cpptraj (amber analysis tool) and calculate the OP for 500mM NaCl. The results I get from cpptraj and our script are identical. As a result, I don't think that there is a problem with the original OP script for amber.

[jmelcr]

ad no.2

-- I am very glad to hear that!

ad no.1

I think that the issue is even imprinted in the literature, eg. paper by Cheatham for the case of KCl.

Auffinger, P., Cheatham, T. E., & Vaiana, A. C. (2007). Spontaneous Formation of KCl Aggregates in Biomolecular Simulations: A Force Field Issue? Journal of Chemical Theory and Computation, 3(5), 1851–1859. https://doi.org/10.1021/ct700143s

and it could just be the lucky case, that the simulation in Gromacs did not experience any nucleation event during its course, but...

... the parameters for my simulations were taken from amber99sb-ildn ffnonbonded.itp (as shipped with Gromacs distribution):

;name  at.num    mass      charge   ptype   sigma     epsilon
  Cl_ff99   17      35.45    0.0000  A   4.40104e-01  4.18400e-01
  Na_ff99   11      22.99    0.0000  A   3.32840e-01  1.15897e-02

which I remember as "shall be" originally by Aqvist, but which I do not find to agree with the numbers from the original paper by J. Aqvist according to the formulae for converting of A-B formalism to sigma-epsilon.

In any case, the above mentioned parameter set I used for these simulations is what is called as "ff99 ions" in the distribution of Gromacs, its origin is only unclear to me now. This and also its sub-optimal performance is the reason, why I decided to abandon these parameters and use a better set by Dang in any future investigation. (as also mentioned in the paper by Cheatham )

;name  at.num    mass      charge   ptype   sigma     epsilon
  Na_d  11      22.98980    1.00      A     0.2350    0.5439   ; DOI: 10.1063/1.466363  ; Dang1994
  Cl_2d 17      35.45300   -1.00      A     0.43387   0.4184   ; DOI: 10.1021/jp064661f ; Dang2006

So, we may have had just different ion parameters in Gromacs .vs. Amber simulations in the end - but this will be sure and clear only once the actual parameters used in Amber simulations are revealed....

[ohsOllila]

As far as I understand @batukav does not see aggregation at 500mM NaCl in Amber simulations and @jmelcr has a gromacs simulation at 586.8 mM (https://github.com/NMRLipids/NMRlipidsIVotherHGs/issues/21#issuecomment-429287333). Therefore, we do not yet have a case where aggregation is seen in Amber but not in Gromacs.

However, there is the difference in the alpha carbon order parameter at this concentration. Do have any idea if this difference is related to the ion aggregation or something else?

The goal of these simulations was to evaluate which force field has the most realistic counterion binding affinity. This seems to be quite difficult for many reasons and this aggregation business makes it even more difficult. I think that we should omit the systems with aggregates from the main discussion and mention these systems aggregate. Maybe we could also show a figure in the SI.

[jmelcr]

I have also compared my simulations of pure POPS bilayer with ff99 ions to those, which are available from @batukav (Zenodo link). The plots of the order parameters from the different simulation engines (and force field implementations, too) are in the following figures

• Gromacs
• Amber (Note that the form factors do not belong to the same simulation in this figure)

The figures show closely matching order parameters supporting that the engines provide same results.

It looks now, that the possible sources of the large warp of the order parameters in the mixed simulations (mentioned in the first comment) may come from any possible differences during the simulation run and/or during the analysis assuming that the system and its topology was prepared correctly.

[mattijavanainen]

This information could be useful in comparison: I simulated a PC/PS membrane in Gromacs with Åqvist potassium I experienced crystallization at 4 M concentration, whereas no crystals formed at 3 M. Details and the files are available at https://doi.org/10.5281/zenodo.1404040.

[ohsOllila]

@jmelcr could you also share the sodium density profiles from the Gromacs simulations of pure POPS mentioned above: https://github.com/NMRLipids/NMRlipidsIVotherHGs/issues/21#issuecomment-431829430

[ohsOllila]

I have now excluded the Lipid17 data with additional counterions from the figures and discussion about the counterion binding affinity due the cluster formation. I also added the figures about clusters into the SI. The commit is this: https://github.com/NMRLipids/NMRlipidsIVotherHGs/commit/6af3f4b3689a794adb8e7cecefefe9a248539782

[jmelcr]

@ohsOllila The data from simulations with pure POPS are now in the pull request in the MATCH repository.

[jmelcr]

And here is the profile simply plotted: Lipid17 + ff99 Na counterions simulated with Gromacs (1000 ns)

[ohsOllila]

I plotted the sodium density profiles from pure POPS simulations with ff99 ions from Gromacs and Amber. The results look sufficiently similar to me:

[jmelcr]

Compares well to my eyes, too! At least we know for sure that the accuracy of the MD engines does not play a major role in this matter.

Are the differences at the edges coming from just a different size of the simulation box?

[ohsOllila]

All the issues discussed here may not be fully clear, but I believe that these are not relevant for the NMRlipids IV manuscripts. Therefore, I will close this issue now.