## Friday, February 6, 2015

### The first draft of the ion-lipid interaction manuscript

[update 14.12.2015] Different force field was used for CHARMM36 ions than I tought.

[update 16.4.2015] I have created an own GitHub Repository for this project. The links are updated.

I have now written the first draft of the manuscript considering the ion-lipid interactions discussed in this blog. All the manuscript files (written with LaTex) are in GitHub. The pdf version the manuscript with ToDo list is also available through link https://www.dropbox.com/s/ebiwdsoj2otg4hp/LIPIDionINTERACT.pdf?dl=0, in addition to GitHub.

I think that there are few simulations more which we definitely have to run for this manuscript:
• Orange model with higher concentration of NaCl (I will do this).
• CHARMM36 model with CaCl$$_2$$ (the data contribution from someone already having this data would be excellent. I if this does not happen, I will do it).
In addition, it would be very interesting to get some results with different NaCl and CaCl$$_2$$ concentrations from the modified CHARMM36 ions by Venable et al. According to the paper, the parameters improve the ion binding to the charged lipid bilayers. I think that the data with these parameters would be highly relevant, however I do not have time to do it now.  The force field by Venable et al. is already used in the current CHARMM36 simulations with NaCl, see comment on 14.12.2015.

All kinds of comments about the manuscript are welcomed from everybody (also from people who are not contributors yet). You can submit your comments by commenting this post, or you can directly fork the GitHub repository and make a pull request.

The authorship of the manuscript will be offered for everyone who has commented the blog, according to the "on the credits post". The final decision on authorship will be made by invited individuals themselves based on their self-assessment of their scientific contribution to the project. The separate authorship query is sent to to all contributors with email and authorlist will be updated according to the replies.

1. Here are the numbers to be inserted in Table I for the DPPC lipid bilayers simulated with the Slipids force field:

Nl Nw NNa NCl tsim (ns) tanal (ns)
600 18000 49 49 100 40

I have also recently posted a comment before with the details on the simulations, but it got deleted... let's see if this one makes it

1. Thanks,

Manuscript is updated. Let us know if you are planning to share some files in Zenodo. Your first message did go to the spam folder for some reason. This has happened previously also for Peter Heftberger for unknown reason.

2. you're welcome,

I thank you very much for explaining me the fate of my first comment. Ok, I will talk to my advisor about this and I will let you know if I can share some files in Zenodo as soon as possible.

2. and here are the details on the simulations of DPPC lipid bilayers performed with the Slipids force field:

The starting DPPC lipid bilayer, which was built with the online CHARMM-GUI
(http://www.charmm-gui.org/), contained 600 lipids, 30 water molecules/lipid, Na+ and Cl- ions (150 mM NaCl). The TIP3p water model was used to solvate the system. AA MD simulations of DPPC lipid bilayers were performed at ten different temperatures (283, 298, 303, 308, 312, 313, 314, 318, 323, and 333 K) using the GROMACS software package version 4.5.5 and the Stockholm lipids (Slipids) force field parameters for phospholipids. After energy minimization and a short equilibration run of 50 ps (time step 1fs), 100 ns production runs were performed using a time step of 2 fs with leap-frog integrator. All covalent bonds were constrained with the LINCS
algorithm. Coordinates were written every 100 ps. PME with real space cut-off at 1.0 nm was used for Coulomb interactions. Lennard-Jones (LJ) interactions were switched to zero between 1.0 nm and 1.4 nm. The neighbour lists were updated every 10th step with a cut-off of 1.6 nm. Temperature was coupled separately for upper and bottom leaflets of the lipid bilayer, and for water to one of the temperatures reported above with the Nosè-Hoover thermostat using a time constant of 0.5 ps. Pressure was semi-isotropically coupled to the atmospheric pressure with the Parrinello-Rahman barostat using a time constant of 10 ps.
The last 40 ns of each simulation was employed for the analysis of DPPC choline and glycerol backbone order parameters.

I hope this comment will stay longer than the previous one... anyway, please let me know if you need further details and/or analysis on these simulations.

Thank you very much for your consideration,

Andrea Catte

3. The MacRog ion data is available in Zenodo at http://dx.doi.org/10.5281/zenodo.14976

The values for the table are as follows:
0 mM - same as in the 1st manuscript for the 50 w/l system (i.e. same system)
100 mM - 103mM NaCl, 288 POPC, 14554 w, 27 Na, 27 Cl, 310 K, 90 ns, 50 ns.
200 mM - 207mM NaCl, 288 POPC, 14500 w, 54 Na, 54 Cl, 310 K, 90 ns, 50 ns.
300 mM - 311mM NaCl, 288 POPC, 14446 w, 81 Na, 81 Cl, 310 K, 90 ns, 50 ns.
400 mM - 416mM NaCl, 288 POPC, 14392 w, 108 Na, 108 Cl, 310 K, 90 ns, 50 ns.

And the values provided earlier are below with signs (for 0mM take from the 1st manuscript):

100 mM:
0.0236, -0.0180
0.0407, 0.0131
-0.1218, -0.2465
-0.2112
-0.1673, 0.0095

200 mM:
-0.0265, -0.0070
0.0106, 0.0046
-0.1121, -0.2532
-0.2279
-0.1814, 0.0239

300 mM:
-0.0355, -0.0108
0.0061, -0.0042
-0.1190, -0.2504
-0.2269
-0.1922, 0.0117

400 mM:
-0.0214, -0.0099
0.0099, 0.0071
-0.1240, -0.2535
-0.2291
-0.1821, 0.0253

And the description of these simulations in the method system could read something like:

The simulation parameters are identical to those employed in our earlier study [citation to 1st paper] for the full hydration and dehydration simulations. The initial structures with varying amounts of NaCl were constructed from an extensively hydrated bilayer by replacing water molecules with ions using the gromacs tool genion. Even at the highest considered salt concentration, the amount of water molecules per lipid after this replacement process was still greater than 50.

1. Thanks, the manuscript is updated. Two questions:

1) What do you mean when you write, for example, 100 mM - 103mM NaCl? What are the two different concentrations?

2) How do you calculate the ion concentrations? The MacRog system with 200mM NaCl has almost the same amount of ions as the Berger system with 190mM NaCl, while MacRog system is significantly larger in terms of number of lipids and water, see table I in manuscript. This is also related to one of the points in ToDo issue number 33.

2. 1) Well the 103, 206 mM etc. are the real values. Due to rounding up to 27 ions, 100 mM turned into 103 mM and the others were constructed as multiples of these 27 ions.

2) I calculate the concentration in the aqueous phase, not based on the volume of the whole system (like genion -conc does).

V_ion = V_water
–>
n_ion/c_ion = n_water/c_water
–>
#ion = #water/c_water*c_ion

With c_water = 55.5 M, #water = 14500 and c_ion = 0.207 M I get #ion=54.

And the ions used with MacRog are from OPLS.

3. 2) This is the difference. I have calculated the total concentration in the simulation box (i.e. using the volume of the system). I thought that Akutsu et al. report that concentration (taking the lipid volume into account), however, after careful reading they seem to report only the concentration of the buffer. Thus, I think that we should calculate the concentration as you have done it. I will do these calculations and update the manuscript soon.

Andrea, how have you calculated the ion concentration from the Slipid simulations?

4. The manuscript is now updated such that all the concentrations are reported as buffer concentrations. I noticed that Andrea had already done this way so only the concentrations for Berger, CHARMM and Orange results are changed. The figures changed a bit but not conclusions. With CHARMM some concentrations became quite large (above 1M) and I think we need some more simulations with this now.

Has anyone tried if the results change when the amount of water and ions are changed by keeping the concentration the same?

5. In the end the concetration issue changed the conclusions. The simulations with 1M of NaCl in CHARMM36 did not change the order parameters with the current accuracy, see:
https://github.com/NMRLipids/lipid_ionINTERACTION/commit/c9cb216348d7c096a3bedf71f8b95fe9828f38c8

4. hello.
I read the ms, have some questions/comments. I confess that I have not read all the posts on this and the previous ms, so I apologize in advance for if I'm repeating somebody else's questions/comments. Also, I apologize for my ignorance on the details of the relevant literature.
1, figure 2 reports results for g1, g2, and g3 CH groups. Yet, as I understand, there are no experimental results for those. Am I missing something?
2, on Orange lipids (or whatever we'll decide to call them, once they will officially exist): they sort of come out of the blue, and there are no parameters published. I am responsible for this lack of info, and the main reason is that my plan is to publish those parameters in a regular paper in some ACS or RSC journal, which do not allow publishing material that is not original. I imagine the readers of the present paper will be puzzled (and annoyed, rightly enough) about this lack of info. Let alone the reviewers, in case you decide to submit this to a regular journal. I am not sure what is the best thing to do, really.
3, does the blog contain links to the topology files that were used for each lipid force field? This would be very useful, it would save me (us?) the time to go look for the files each time.
Thanks.

1. 1. In the presence of NaCl there is one experimental point for g3 which is in the figure. For g2 and g1 there are no measurements with NaCl concentration, however it is very unlikely that these would change since all the other bilayer structrural properties are not affected by NaCl with these concentrations. The experimental values for g1 and g2 without ions are shown in the figure, but not seen very clearly now.

2. There is now two reasons to include the Orange results. First, the results gives further support for one of the main points of the paper, i.e. that the order parameter change and cation partition relation (the electrometer concept) works also in simulations (with all models, despite of their differences). Secondly, compared to the Berger model only lipid model has been changed (the water and ion parameters are the same). This means that one can affect ion partition significantly only by changing lipid model. Otherwise one could argue, for example, that the ion partition cannot be fixed with the SPC water model or with united atom model.

If we would remove the Orange, there would be four cases instead of five supporting the first point. The second point is more like a side note (but relevant one). Thus, I think that the Orange data improves the work significantly, but is not critical.

I think that the best to do is that we keep the results in the manuscript, and you tell as much details as you think that you can. If someone gets too puzzled and annoyed after that, then we remove it. Personally, I do not believe that you get much problems if you put your parameters, for example, to Zenodo and cite them as a beta version of the model in the actual publication of describing the parameters. There are also embargo options in Zenodo. One option may be also that you share only the trajectory and *tpr file in Zenodo. In this case the model itself is not yet available, but the trajectory is, and suspicious reader of the paper can take a look at the trajectory themselves. This is what MacRog people did with POPC parameters which are not publised yet: http://dx.doi.org/10.5281/zenodo.13497.

3. Most parameters and trajectories are available in our Zenodo collection: https://zenodo.org/collection/user-nmrlipids, however, addition of many ion simulations is still in progress (see table I, in the ion-lipid manuscript). Some, but not all files are also available in the github repository: https://github.com/NMRLipids/nmrlipids.blogspot.fi.

2. 1. After reading the publications more carefully I realized that there is a mistake in Fig. 2. The experimental point for g3 with 350mM ion concentration is actually measured with CaCl_2, not with NaCl. I will soon update new figures with more experimental points in the right locations.

3. 1. I happened to notice that for g2 we already seem to have found an experimental value (0.204705882352941) at [NaCl]=150mM, see:

http://nmrlipids.blogspot.fi/DATAreportediINblog/DPPC/NaCl/EXP-Scherer1989.dat

This was in a racemic mixture, but as there was only one splitting at 26.1 kHz, the isomer does not seem to matter for the g2 value, and thus we could in my understanding use this.

5. I have now updated more experimental datapoints into the manuscript. There are now results with little bit different temperatures and also with POPC. Also data for g3 is now moved to the CaCl_2 figure where it belongs.

Surprising observation is that the effect of CaCl_2 is stronger in the DPPC bilayer compared to POPC. For NaCl the points overlap. Only discussion I can find about this from the original works is the sentence in the introduction by Altenbach1984: "The larger spacing between the polar groups should be reflected in a different binding affinity for metal ions.".

Anyway, I will run the CHARMM36 simulation with CaCl and POPC. If it gets somewhat close to experiments we can run with DPPC as well, and see if we get the same difference as in simulations.

6. I realized the results calculated with orange lipids are not present on the 1st
paper but they are present on the second. I don't actually have a problem
reporting those results on both papers, and also I don't have a problem putting
the tpr and trajectory files on zenodo. I just wouldn't want to put the itp files
online publicly, because of possible restrictions in future publications and, more important, those parameters won't be definitive as all results with gromacs 5 are different and I'd like to work with gromacs 5 in the future.
(This probably applies also to all other results in the blog, but I don't think order parameters in the head group will be very sensitive to the change).

1. Excellent. My idea has been that we would not put Orange results into the first paper since in that one their role would be just to report how it works, and this is not relevant for the outdated version. However, in the second paper, even the outdated results are revelant due to the reasons I mentioned above.

I have now also some results with Orange which are not reported yet (higher NaCl concentration). I will collect these, update the manuscript and put the trajectories and tpr files available.

About Gromacs 5, why would you think that the results would be different?

7. Short answer: it is possible to reproduce the old results with gromacs 5 but it is not convenient.
Long answer: a very significant gain in speed is obtained by using Verlet neighbor lists, which we never used before. Verlet neighbor lists also have another (even more important) benefit: a significant improvement in energy conservation. Since this method is faster and implements better physics, I think all users will start using it.
In addition to this, Verlet neighbor lists are normally used in combination with shifted potentials (this is recommended in the new version, and is computationally very efficient).
Based on some tests I have run myself, basically all results on even the most simple systems are different, when using Verlet neighbor lists and shifted potentials. The effect of having good energy conservation is absolutely dramatic in Martini, where one needs to change the actual cutoff (!!!) in order to obtain results similar (but not the same) as before. Since all current atomistic force field for lipids are cutoff-dependent, results will change for all current atomistic force fields.
As for the shifted potentials, they should not alter forces (the derivatives of the potential should be the same), but they do alter the potentials for all non-bonded interactions. Which means that, for example, various thermodynamic quantities (heat of vaporization, etc) will be affected. (This is something you don't care about, when you focus on order parameters only).
Based on some preliminary tests on Berger lipids, results are significantly different (in terms of area per lipid at different temperatures, for example).

Having said this, I don't think the effect will be major on order parameters in the head group, because I suspect order parameters in the head group are not very sensitive to energy conservation, thermodynamic properties, and even area per lipid, phase behavior, etc. (I can be wrong on this, I haven't tested this statement; I will do tests as soon as I find some time).

8. Hi,

I modified the github manuscript to add results of the tests with ion polarization effect. The same text can also be found here:

https://www.dropbox.com/s/yi80e3va75o5veo/Lipid_ion_interactions_maatta.pdf?dl=0

-Jukka

1. Thanks! I did add the figures from Zenodo also into the manuscript. I also plotted the order parameter changes with and without scaling for the Berger model, see https://github.com/NMRLipids/lipid_ionINTERACTION/commit/8e29e04cbce4105ccc0b72513ba7980f6e067719.

I think there are two issues now.

1) I did not plot the OPLS compatible model order parameter changes since there was no result without ions for this. Is there actually difference in this to the regular berger without additives? If there is, could you run that as well?

2) Other issue is that the scaling improves the results, however, ion concentrations are quite low. Based on the previous results, and the density plots, I would guess that at least regular Berger will have too large order parameter changes with larger concentrations. Could you run simulations with both scaled models also with larger concentrations, something between 500mM-1M?

2. Hi,

Thanks for adding the zenodo figures to the manuscript! I did not find a quick way to add files other than raw text to the github folders myself.

1) There aren't any large discrepancies between the two models. Seems like g2 is affected a bit. Here is a simulation of pure DPPC in water with OPLS-based Berger:
DOI: 10.5281/zenodo.17237

2) Simulating higher ion concentrations is a really good idea. So here are files for DPPC with 1 Mol NaCl.

Berger: DOI: 10.5281/zenodo.17210
OPLS-Berger: DOI: 10.5281/zenodo.17208

and with scaled ions:
Berger: DOI: 10.5281/zenodo.17228
OPLS-Berger: DOI: 10.5281/zenodo.17209

Again it seems that the OPLS-Berger with scaled Åqvist ions gives the most realistic results.

-Jukka

3. From command line you can push files with any format into Github.

I did now add your new results into the figure in SI. However, the results with 1M looks little bit weird, except for Berger-OPLS:
https://github.com/NMRLipids/lipid_ionINTERACTION/blob/master/Fig/OrderParameterIONSchangesSCALED-eps-converted-to.pdf

When I was looking at the numbers it made me suspect that maybe there has been some mistake in the analysis this time. Could you check that?

Otherwise, I think that we should actually put the Berger-OPLS results (and also the Berger with 1M) into the main body of the manuscript, and leave only the scaled results into SI.

4. Hi,

Thanks for noticing, you are correct! I happened to copy a wrong ".hdb" file for protonation of the UA coordinates, so all the analysis was done on wrong hydrogens. Here should be the correct order parameters:

https://www.dropbox.com/s/bdk87eecctw2i6s/orderparameters.zip?dl=0

I'll probably re-upload the 1Mol zenodo trajectories as well since the protonated.xtc files are now redone.

Moving the Berger-OPLS results to main text sounds good. In that case it would probably be useful to run at least one simulation with CaCl2 as well for comparison.

-Jukka

5. I have now updated the manuscript. I also added your results without scaling to the main text. The figures got little bit messy but those require some work anyway. The scaling seems to improve the results in both cases, however, I would say that the binding is too strong even with the scaling.

9. I have now added results from different NaCl concentrations from Orange model, larger NaCl concentration with CHARMM36 and CaCl simulation with CHARMM36 into this manuscript. I have commented the results in todo list text, but not changed the text yet, see:
https://github.com/NMRLipids/lipid_ionINTERACTION/commit/c9cb216348d7c096a3bedf71f8b95fe9828f38c8

The major points are these
- Against previous conclusions, we cannot say now that the Na would bind too strongly in CHARMM36.
- The first CaCl_2 results with CHARMM36 are in pretty good agreement with experiments.
- The overestimated effect of CaCl_2 in Orange model seems to arise from the lipid response, not from too strong affinitiy.

10. Just saw this new paper on JCTC that might be relevant: http://dx.doi.org/10.1021/acs.jctc.5b00540. They used CHARMM36 FF and metadynamics to look at cation binding to phospholipid membranes.

11. I am now back to work and I remade the density figures for the ion-lipid interaction paper. I am also trying to move to the new workflow and I have created the GitHub issues in which the figures are discussed and linked these into the manuscript. So the idea is that you can use also the GitHub issues (in addition to the blog) to comment on the manuscript.

These are the created issues:
https://github.com/NMRLipids/lipid_ionINTERACTION/issues/3
https://github.com/NMRLipids/lipid_ionINTERACTION/issues/4

And the manuscript is in the regular location:
https://github.com/NMRLipids/lipid_ionINTERACTION/blob/master/Manuscript/LIPIDionINTERACT.pdf

And the new workflow is described here:
http://nmrlipids.blogspot.fi/p/how-to-follow-and-make-contributions.html

12. Here are 200ns trajectories of POPC bilayer with 0M, 0.15M, 1M NaCl obtained using LIPID14/AMBER and Ulmschneider/OPLS force fields.

LIPID14/AMBER99SB-ILDN (298.15 K, 128 POPC, 200ns):

1) 0 M NaCl -----> https://zenodo.org/deposit/60757/
2) 0.15 M NaC --> https://zenodo.org/deposit/60737/
3) 1 M NaCl -----> https://zenodo.org/deposit/60524/

Ulmschneider/OPLS (298.15 K, 128 POPC, 200ns):

4) 0 M NaCl -----> https://zenodo.org/deposit/60764/
5) 0.15 M NaC --> https://zenodo.org/deposit/60743/
6) 1 M NaCl -----> https://zenodo.org/deposit/60749/

1. Thank you for the major contribution. Your commit:

https://github.com/NMRLipids/lipid_ionINTERACTION/commit/afe932ebdff2e01479b3803c9cabda0a9842afe7

is also a good example how the current scripts can be modified to analyze new data.

I did update the density figure of Na and Cl densities to put the labels in right places:

https://github.com/NMRLipids/lipid_ionINTERACTION/commit/652df544e011f415586e831ba9c70b26f7cf95dd

Now the new figure is also in the manuscript. With the first look it seems that in Lipid14 the partition of ions is quite similar to CHARMM36, while in Ulmscneider it is significantly stonger. Next we need to calculate the order parameters. The scripts for CHARMM36 in the script folder should be quite easy to modify to do this. I am little bit busy with other things now but in couple of weeks I will also continue working on this.

We also need to put these simulation details to table I and SI. Could you please deliver the doi references of trajectories for that? I was not able to open the deposit links.

2. Hi,

Thanks! These are corrected links to the trajectories on Zenodo with the DOIs.

LIPID14/AMBER99SB-ILDN (298.15 K, 128 POPC, 200ns):

1) 0 M NaCl -----> https://zenodo.org/record/30898
DOI: 10.5281/zenodo.30898

2) 0.15 M NaC --> https://zenodo.org/record/30891
DOI: 10.5281/zenodo.30891

3) 1 M NaCl -----> https://zenodo.org/record/30865
DOI: 10.5281/zenodo.30865

Ulmschneider/OPLS (298.15 K, 128 POPC, 200ns):

4) 0 M NaCl -----> https://zenodo.org/record/30904
DOI: 10.5281/zenodo.30904

5) 0.15 M NaC --> https://zenodo.org/record/30892
DOI: 10.5281/zenodo.30892

6) 1 M NaCl -----> https://zenodo.org/record/30894
DOI: 10.5281/zenodo.30894

I will try to calculate the order paremeters during the week.

3. Can you specify how long equilibration period you had before this 200ns, or does it include the equilibration?

4. All runs were preceded by 5ns NPT equillibration (V-rescale tcoupl and Berendsen pcoupl. I updated the description on zenodo.

5. I have now added the Lipid14 order parameter results into the manuscript and pushed related files into GitHub.

The results for Lipid14 seems pretty good. It might be reasonable to run the CaCl simulations for this model.

6. I have now added the Ulmschneider order parameter results into the manuscript and pushed related files into GitHub.

Regarding the ion partition, the results are in line with all the others.

However, there are some issues with the glycerol backbone order parameters related to the previously reported values, see
https://github.com/NMRLipids/lipid_ionINTERACTION/issues/8

13. To me it seems that the manuscript presents three important findings:

1) Na-ions do not penetrate below the PC-headgroups, in contrast to what has been reported in many recent simulational and experimental works.

2) Ca-ions do penetrate, in qualitative agreement with previous works.

3) The electrometer concept devised by Seelig and coworkers in the 1980's seems to work exactly as they proposed, that is, there is a direct relationship between changes in the choline order parameters and the amount of penetrant charge.

I am proposing that in order to make our message more easily accessible to the readers, we would reorganize the manuscript in the following way:

First we would introduce the electrometer concept and demonstrate that it works also in simulations. Then we would turn to Ca-ions, whose interaction with PC-lipids is qualitatively agreed upon by everyone the field; here we would further demonstrate the power of the electrometer concept. Then finally, we would move on to the most controversial of our findings, the lack of Na-ion penetration; we could make a strong case here, based on the just previously
demonstrated electrometer concept and its correct working in the Ca-ion case.

I am willing to do this reorganization, but as it requires considerable work, I thought it best to ask the opinions of others before embarking on it.

If the suggested approach is agreed upon, I would also suggest renaming the manuscript to "The electrometer concept and binding of cations to phospholipid bilayers".

1. I agree with the reorganization.

One important point we have about the electrometer concept is the negative sign of beta order parameter which leads to the decrease of value with penetrating charges. Seelig et al. originally measured only absolute value which increased. This must be stressed out in the manuscript and is also a good reason to discuss electrometer concept separately in the beginning.

If you start to work with the text reorganization, please push new versions to git as often as possible so that we can see what is going on. I will concentrate now getting all the data for the figures etc.

2. Ok, good, I will then now start doing this.

First push (just change of title, and fix of typos in the beta sign issue) is now in github.

3. I have been now little bit rewriting the manuscript. It is now close to the order suggested by Markus but not quite. I am now waiting the results with CaCl_2 to finalize the text.

14. Here’s data with Slipids and NaCl
with ion model that uses scaled-down
charges. The ions bind less to the
head group region with this model.

Systems consisted of 200 POPC, 45
waters per lipid and then the ions
on top of that. Simulations were
run at 310 K, simulation parameters
(mdp) are available in Zenodo.

The scaled ions are taken from:
Kohagen et al., Article ASAP (2015)
DOI: 10.1021/acs.jpcb.5b05221

All files related are available in
Zenodo: DOI: 10.5281/zenodo.35193

And the values are:

130 mM:
-0.01993 –0.05127
0.05761 0.04552
–0.01143 –0.01802
–0.02484
–0.15966 –0.08670

300 mM:
–0.02682 –0.05527
0.04875 0.04594
–0.00285 0.00403
–0.01010
–0.15595 –0.06888

500 mM:
–0.03095 –0.06046
0.04896 0.03143
–0.01748 0.00373
–0.02690
–0.14752 –0.08101

750 mM:
–0.03954 –0.06720
0.03964 0.02927
–0.02498 0.00389
–0.03586
–0.14864 –0.08301

1000 mM:
–0.04166 –0.07383
0.03871 0.03089
–0.01091 0.00955
–0.02012
–0.15620 –0.06708

1. And here's data for the same system (130 mM) with the ions by Smith & Dang (DOI: 10.1063/1.466363, integer ion charges):

http://dx.doi.org/10.5281/zenodo.35275

2. Thanks. I did add order parameter results with Smith & Dang ions into the manuscript. I am still working on the density profiles and the scaled ones.

One question about Smith & Dang results (http://dx.doi.org/10.5281/zenodo.35275):
Is there pre-equilibration before this 100ns, i.e., which number we should put in the t_sim column in the Table 1?

3. Trajectories for both Smith&Dang and scaled ions are 100 ns long with dt=2fs. Prior to that, a 5 ns simulation was run with a shorter time step of 1 fs. The area per lipid curves showed no drift during the 100 ns trajectories so in the values I reported above for the scaled ions, the whole 100 ns trajectories was employed for the analyses.

4. I have now included also the results with scaled charges into the supplementary information.

5. I also have simulations with CaCl_2 that has scaled-down charges by the same authors ( DOI: 10.1021/jp5005693 ). Here, the same bilayer of 200 POPCs was simulated with a CaCl_2 concentration of 450 mM. I ran these simulations with both

Slipids: DOI: 10.5281/zenodo.45007
and
Charmm36: DOI: 10.5281/zenodo.45008

Prior to the 100 ns trajectory in zenodo, the systems were equilibrated for 500 ns (Slipids) or 700 ns (Charmm36), i.e. until the number of bound ions reached equilibrium.

In case you think these results should be included in the MS, I will write a description of them in the Methods section.

6. Thanks!

I added the Slipids results to the main text since we did not have those with any CaCl model. I added the CHARMM36 simulations to the supplementary material. The text has been modified accordingly, see:
https://github.com/NMRLipids/lipid_ionINTERACTION/commit/405fd201abee8ca2f700ee0624eab6d3b1c05919#diff-ca2ede0e96d155a8ca8dbca9157a819f

There is also ToDo point in the Slipids simulation details where I think some description of Slipid CaCl simulations could be written.

15. Matti Javanainen pointed out in Skype that the CHARMM simulations I have ran with NaCl already contains the ion force field by Venable et al. which I suggest to be tested in this post. This was true indeed. The nbfix.itp file by Venable et al. is included in the files shared in MacKerell lab webpage even though the publication is not listed as a reference. Thus this test is not needed. I have now updated the reference information into the table in manuscript.

16. here are the order parameters with errors calculated using the scripts available on Zenodo. each simulated system contains 600 DPPC molecules, 150 mM NaCl and 30 waters/lipid. the ion model is the one described by Smith & Dang (ref. 76 in the current version of the manuscript). the analysis has been performed on the last 40 ns of each 100 ns trajectory. the order parameters are listed with the following order:

alpha
beta
g3
g2
g1

alpha, beta, g3 and g1 are averages of two values obtained for each hydrogen and their errors are the sums of the errors of the order parameters of each hydrogen:

283 K

0.0297987 0.01332877
-0.0483905 0.01240381
-0.0506754 0.0234372
-0.0639908 0.0122343
-0.1437995 0.02067088

298 K

0.0296881 0.01002405
-0.04606365 0.00904883
-0.0511725 0.01970425
-0.0684736 0.0101523
-0.1411595 0.01700452

303 K

0.0302904 0.00833764
-0.0426856 0.00765966
-0.0507356 0.01804765
-0.0697055 0.0098925
-0.144534 0.01530093

308 K

0.03413 0.00688235
-0.043259 0.00618568
-0.04197615 0.01447939
-0.0642038 0.00776303
-0.135129 0.01239031

312 K

0.03614115 0.00600576
-0.0457234 0.00525457
-0.0466989 0.01299161
-0.0641617 0.00744506
-0.1467645 0.01032069

313 K

0.03997975 0.00528805
-0.0448399 0.00478971
-0.04023395 0.01076638
-0.0522406 0.0061841
-0.1388265 0.00892683

314 K

0.04148735 0.00533955
-0.0455148 0.00468357
-0.048194 0.01090664
-0.0595314 0.00648737
-0.143842 0.00875961

318 K

0.0393653 0.00476622
-0.04094225 0.00399588
-0.0343535 0.00944234
-0.0530931 0.00542407
-0.1315395 0.00773388

323 K

0.03834715 0.00435412
-0.0401592 0.00370016
-0.0362901 0.0085152
-0.0596518 0.00509698
-0.1308845 0.00726474

333 K

0.04167525 0.00347418
-0.03701765 0.00296599
-0.03745435 0.00668386
-0.0519819 0.0039
-0.1299625 0.00545539

in the next comment the density profiles will also be provided for each temperature

1. Thanks for the results. However, for the manuscript I need the order parameters of both hydrogens attached to the same carbon, not the average. Only the results in 323K will be included since at this time we will not discuss about temperature dependence. Could you share the results for this one including all the values?

2. you're welcome for the results and sorry for the lateness of my response but I've just seen your comment. I see.. yes, I could share the results with all the values at 323K:

-0.0247313 0.00199584
-0.0555871 0.00170432
0.0453406 0.0023083
0.0313537 0.00204582
-0.0491411 0.0047118
-0.0234391 0.0038034
-0.0596518 0.00509698
-0.111834 0.00413055
-0.149935 0.00313419

the simulation details (the ion force field is the one by Roux as reported in the last comment) and the sequence of these values are the same reported in the previous comment

3. Thanks. I have now included all this information into the manuscript.

https://zenodo.org/record/35379

the partial density profile of the lipids needs to be scaled down by 10 and partial density profiles of ions should be scaled up by 50 in order to get profiles comparable to those reported in Fig. 3 of the manuscript.

18. after checking more carefully the force field parameters employed for the simulations with the Slipids force field I have noticed that the ions force field parameters are those from Roux et al. articles:

Beglov, D., Roux, B., Finite Representation of an Infinite Bulk System: Solvent Boundary Potential for Computer Simulations, J. Chem. Phys. 1994, 100, 9050-9063.

Roux, B., The molecular basis of the valence selectivity of the gramicidin channel: A molecular dynamics free energy perturbation study, Biophys. J. 1996, 71, 3177-3187.

which are also cited in:

Joung, I. S., T. E. Cheatham III, Determination of Alkali and Halide Monovalent Ion Parameters for Use in Explicitly Solvated Biomolecular Simulations, J. Phys. Chem. B 2008, 112, 9020-904.

please accept my apologies for the previous mistake and cite the above references for the ion force field of MD simulations of DPPC lipid bilayers with 150 mM NaCl

19. Hi, I am worried that the simulation times with CaCl_2 are too short for the ion binding to reach equilibrium. I have run POPC with 450 mM of CaCl_2 for 2 µs with both Charmm36 and Slipids force fields, so this could be checked. Here, the default Ca ions of the force field (i.e. no scaling of charges) are employed. The files are here:

Charmm36: DOI:10.5281/zenodo.51185
Slipids: DOI:10.5281/zenodo.51182

1. I added the results in Fig. 2. I also moved the scaled Slipid result from Fig. 2 to Fig. 6. Now we have comparison between scaled and non-scaled result also for Slipids, thus the discussion in Appendix A must be rewritten. I am working on this. Also Figs. 5 and 8 with densities are updated.

In the long CHARMM simulation almost all the Ca++ ions seems to bind in bilayer. Also the order parameter changes are larger than would be expected from the trend in shorter simulations (see Fig. 2). We should discuss this somehow. Would you have some analysis which would somehow show how things are changing as a function of time in these long simulations? For example, amount of bound ions or something like that.

2. From my experience with Berger, Ca++ binds rather quickly (in nanosecs) and then never leaves unless you have a very concentrated solution.
I'm not sure what is the concern in the equilibrium here - if it is merely the lack of unbinding, I wouldn't be so concerned. By running super long simulations, one would just obtain a more accurate estimate of otherwise a bad model.

3. Matti Javanainen sent me the bound ions as a function of time from long CHARMM and Slipid simulations reported above. The results are now in GitHub:
https://github.com/NMRLipids/lipid_ionINTERACTION/blob/master/scratch/contacts_ca_charmm36.png

https://github.com/NMRLipids/lipid_ionINTERACTION/blob/master/scratch/contacts_ca_slipids.png

(for raw data see files with the same names but endings *txt and *xvg)

The most Ca seems to enter within first microsecond but even after 2 microseconds there seems to some further binding. If one plots only the first 200ns of the data, it may seem that the second 100ns would be already equilibrated.

I think that in the cases where Ca binding is concluded to be too strong (as is the case in almost all simulations in NMRlipids II), this is not critical since the longer simulations would just show even stronger binding.

However, when approaching models with more realistic binding, this must be carefully taken into account. Even if results would look good after 200ns, there may be overbinding after 1000ns.

Only experimental estimate about binding dynamics I have seen is that "The binding of Ca2+ to the membrane surface is a highly dynamic process with a Ca2+ residence time at an individual lipid head group of less than 10^-5 s." [Altenbach84, Biochemistry, 23, 3913]
Let me know if you know more.

4. I have now included the data from Matti and some discussion in the supplementary information. This is now also mentioned in the caption of Fig. 2.