## Tuesday, December 22, 2015

### Towards submission of NMRlipids II publication (lipid-ion interactions)

The manuscript from NMRlipids II project "The electrometer concept and binding of cations to phospholipid bilayers" is (in my opinion) quite close to the submittable version. There are still some details to be fixed which are listed in the ToDo list in the manuscript.

All kind of comments are welcomed also from people who are not contributors yet. My goal is to get manuscript submitted January 2016 so this is almost the last change to present major comments. Manuscript will be submitted soon, all modifications and comments should be done before 19.2.2016.  Manuscript will be submitted soon. As in all the subprojects, the files are available in the project GitHub repository and contributions can be done by commenting this blog or using GitHub. Especially, if you are doing language or other text editing the most convenient way is to directly modify the text file and make a pull request.

My current idea is to submit the manuscript first to the Chemical Science which is a high impact journal by the Royal Society of Chemistry. By doing this we would find out their policy towards open collaboration publications which was also the motivation to submit the NMRlipids I publication to a American Chemical Society journal, see discussion in . All comments related to this are also welcomed.

1. After carefully reading the MS (which I mostly liked a lot), I have two major suggestions:

1) A plot showing the correlation between the ion adsorption onto the membrane and the induced order parameter changes could nicely capture the message of the paper in one figure. Of course it is not clear how the quantities used for such a plot should be defined since the employed ion concentrations are not equal for all models. And even if they were, the choice of the parameters to describe ion partitioning as well as the order parameter change are somewhat ambiguous.

2) For a high impact journal, the MS seems to focus too much on comparing lipid force fields. As such, it seems to fit better in JCTC or some other simulation-oriented journal. If we aim for Chemical Science, the key question should probably not be "what force field best describes ion binding" but rather the highly general "do sodium ions bind to PC membranes". I'm not sure whether both questions can be effectively addressed in the same paper. In the latter case the force field details are not important and perhaps the discussion could be written using some abstract names such as 'adsorbing model' and 'non-adsorbing model' instead of force field names to make it readable for the general audience. In this case the discussion on the quality of the simulation models should perhaps go to the SI. I would then also try to emphasize the confusion in the field regarding sodium adsorption and change the wordings from e.g. "support the traditional view" into something more exciting since papers that support traditional views are probably not interesting enough for Chemical Science.

However, I think that a more technical paper with emphasis on comparing the force fields is also a good idea yet then the target journal should be JCTC.

1. 1) Indeed, this kind of figure nicely demonstrates the electrometer concept and such figures from various experimental data are plotted by Seelig et al. [http://dx.doi.org/10.1021/bi00312a019, http://dx.doi.org/10.1021/bi00398a001]. I have not been highly motivated to do such figure since there is no data from simulations which gives realistic order parameter changes with binding charges (the only results in agreement with experiments show effectively no binding and these cannot be used to make this kind of figure). I am not sure if plotting the correlation between unrealistic ion binding behaviour and order parameter changes from simulations would give any new insight of the system since the real correlations from experiments are already known from several systems (see references above). However, I understand the explanatory power of such figure and one may argue that the readability of figures by Seelig et al. could be improved. Thus, if someone makes such a figure (or at least
script which can be used to make the analysis) I think that such figure could be included but I am currently too busy to do it myself.

2) Does the manuscript really give impression that the main question we address is that "what force field best describes ion binding"?. The whole introduction discusses the confusion about the sodium binding in the literature. There is only two paragraphs in the Results and Discussion section which somehow compare force fields with sodium included and Calcium results from different force fields are compared in couple sentenses.
In Conclusions section the comparison is concluded with two sentences: "The comparison reveals that most models overestimate the Na+ binding, only Orange, Lipid14 and CHARMM36 predict realistic binding affinity. None of the tested models has the required accuracy to interpret the Ca2+/lipid stoichiometry or induced atomistic resolution structural changes."
Some amount of comparison is required to convincingly show that the ion binding is related to order parameter changes despite of incorrect atomistic resolution structure (order parameters without ions). I agree that the amount of force fields we have is more than enough, however, this is just making the point stronger. However, I understand that the amount of force fields analyzed may preset the mind of the reader to look at the paper as force field comparison. This probably decreases the probability to get accepted
into the "high impact" journal, however, I do not think that this is a proper reason to decrease the information content or destrength the
evidence for the main point.

The whole introduction already discusses the mentioned confusion and the second paragraph in conclusions summarises our conclusions and their relation to the results which has led to the confusion. Maybe we should, however, emphasize the role of MD simulations as the starter for the whole confusion. This would help readers to understand the importance of the discussion about simulation models beyond the simulation community. The phrase "traditional view" could be changed but we must respect the fact that we are not saying anything new compared to the work reviewewed in 1990. At the same time we could emphasize that this "traditional view" seems to almost completely ignored in the literature after 2003.

My main argument to get this work published in Chemical Science would be that we are showing that some conclusions presented in "high impact" journals should be seriously reconsidered, e.g. strong cation binding. Since these influental conclusions have reached high visibility, also the major corrections should reach similar visibility. Other argument would the usage of open collaboration whose functionality should be interesting for very wide audience. However, it is not very likely that the community would appreciate these arguments, thus it is fine for me to send it directly to the PCCP. I would like to submit to the RSC journal to see what their policy is for this kind of work.

2. I support Matti's suggestion of a 'plot showing the correlation between the ion adsorption onto the membrane and the induced order parameter changes'.

Although the plot will not give much new insight on the real system, as Samuli pointed out, it will lay a robust foundation for our approach of using the molecular electrometer for showing that Na-binding does not occur in real life.

However, as Matti pointed out 'the choice of the parameters to describe ion partitioning as well as the order parameter change are somewhat ambiguous.'

How to measure the amount of bound charge? As the molecular electrometer reports on the charge bound to the membrane, I propose we integrate the ion charge distribution from the center of the membrane until the density maximum of the g3-carbon.

The order parameter change would in my opinion simply be the relative change, shown already in the Fig. 2 of the manuscript.

3. Whatever choice is done, we must show that the conclusions done do not depend on the details of the choice within reasonable limits.

4. Hi,

I like what Matti and Markus are suggesting.

One possible alternative could be to calculate "ion coordination number" of ions whose "electrostatic interaction with the membrane is larger than thermal fluctuations". Or in other words, ion coordination number at distances below Bjerrum length. In practice, lets calculate the rdf for ions from the g3 carbon. Then integrate that rdf from 0 to 0.7 nm and see the amount of ions. It should be about 0, right? The 0.7 nm comes from Bjerrum length for two charged ions in aqueous solution at 300 K. Actually this assumes that the membrane has charges too so it is sort of a stronger argument of binding (the Bjerrum distance 0.7 nm is an overestimation).

-Jukka

5. Hi Jukka,

we could try out what you suggest. It might work. However, it might be that going one Bjerrum length to the water phase will already include a considerable amount of anions... best would be to try. I am currently doing what I suggested above. If you are interested in trying the approach you suggested, we could compare. If both give roughly similar behavior, we would also have a way of demonstrating Samuli's point that 'we must show that the conclusions done do not depend on the details of the choice within reasonable limits.'

6. I read Seelig's 1987 paper where the concept molecular electrometer is first mentioned, and noticed that he stresses that it is not limited to ion binding, but can be used to study "any process that modifies the electrical properties of the membrane surface".

Motivated by this, I calculated the surface charge densities including water and excluding water for our CHARMM36 POPC systems.

It is not that easy to see, however, what is the correlation between the changes seen in these profiles and the changes in order parameters (Fig. 2 in manuscript).

To this end, it seems best to stick with the original plan and just aim to correlate the changes seen in the amount of bound Na/Ca with the changes in the OPs.

7. Hi,

Here is a quick calculation for the ion rdfs starting from g3 carbon (https://www.dropbox.com/s/4egykgovisjspbb/rdf_cn.png). I studied the Berger-OPLS DPPC system with scaled (no binding observed) and not scaled (binding observed) ions.

Calculating just the cumulative sodium ion RDF and cutting it at 0.7 nm gives on average 0.33 bound ions (with scaling), and 0.61 bound ions (without scaling). If we set an arbitrary threshold to 0.5 for ion binding, then this works, but I have to agree that this is quite sensitive to the cutoff parameters.

-Jukka

8. Hi,

We should definitely repeat that quote from that Seelig paper in our manuscript to emphasize that this procedure shows a general response to changes in membranes electrical properties. I think currently many people use additional molecular probes to measure electric respones, yet the molecular electrometer idea is a "non-invasive, non-perturbing" alternative! If we repeat that quote in the manuscript, it could help to raise interest and impact of the paper even further. At least it could make it more clear that we're also dealing with another general concept (of any electrostatic response) besides our main point of ion binding.

-Jukka

9. Hi Jukka, it seems from your figure https://www.dropbox.com/s/4egykgovisjspbb/rdf_cn.png
that there would be more binding with scaled ions than with non-scaled. Maybe there is mistake in the labels?

10. You are right! Thanks for this. I re-checked, and the labels should indeed be switched.

11. I finished yesterday evening the first draft version of the ΔOP -vs- bound Cation Charge Density plot.

I think it looks promising, showing that in all models there indeed is a clear correlation between the bound cation surface density and the change of the headgroup (α, β) order parameters. Also, this correlation does not seem to depend heavily on ion type (Na/Ca).

I proceed to make a more careful version.

12. About Jukka's figure: From the figure we see that whatever cut-off you have, you will always have more binding in non-scaled case (remembering that the labels should be switched).

It should be noted also that the non-scaled Berger-OPLS simulation has a measurable binding so it is not a "non-binding" example. From the manuscript, page 7: "The order parameter changes and Na+-binding affinity are decreased by the charge scaling but yet overestimated with respect to the experiments as seen from Figs. 5 and 6. Thus the overestimated binding affinity cannot be fixed by only scaling charges."

So it seems to me that from Jukka's figure we see that the relative binding affinities do not depend on the chosen rdf cut-off.

Maybe similar demonstration would be good also for the method Markus used?

13. The plot(s) are not yet finished, but I thought it is good to share the current situation, as I can not continue working on them until after Feb 29th.

I calculated the bound cation charge and the corresponding order parameter change separately for each leaflet in our CHARMM36, Lipid14, and MacRog systems. I did this with three different integration limits: until the g3-carbon, until the phosphorus, and until the alpha-carbon. I think that these plots show that the conclusions drawn from a plot like this will not depend on the details of the chosen integration limit. I also think that a sensible integration limit is until the phosphorus.

I also pushed my scripts to github, so others could continue to work on the plots before the 29th.

14. The plots look very nice. Was the idea that these would be included in the publication? This will somewhat postpone the submission but these figures might make the manuscript significantly better. Would you have some estimate how much time you need to finish the figures? I have currently very limited amount of time to be used for this.

15. Thanks. Yes, my idea was that these would be included in the publication to demonstrate that the molecular electrometer is a robust concept that qualitatively holds also in simulations.

The process of creating the plots is not yet fully automatized. That is, I have scripts (see the git) calculating both the bound cation charge and the OPs per leaflet, but a script for putting these into a combined data file is still lacking. I expect this to be few hours of work, so the plots could be expected to be ready next week.

16. Including these plots will also call for some restructuring of the manuscript. I have been thinking something something along the following lines (for Results And Discussion).

1) Molecular electrometer is a robust concept. Using ΔOP-vs-ΔBCC (OP = Order Paramer, BCC = Bound Cation Charge) plots we demostrate that the molecular electrometer is a strong effect, that is, qualitatively reproduced in simulations even with rather inaccurate force fields. (A similar strong effect was the reorientation of the headgroup upon dehydration in the NMRlipids I paper.) Also, the plots will demonstrate that the effect does not depend (strongly?) on ion type, as Na and Ca fall effectively on the same line in each force field.

2) Ca-binding is qualitatively correct in simulations. Here I would show the dependence of ΔOP and ΔBCC on [CaCl2]. (This can be done in one plot, as ΔOP is negative and ΔBCC positive.) The plot demonstrates that in all the simulations, like in the experiments, Ca changes the OPs. And, in all the simulations, Ca binds to the membrane. As everyone agrees on this, it will act as an introduction for the reader to the next part.

3) Strong Na-binding is a simulation artifact. Here I would show, like for [CaCl2] above, the dependence of ΔOP and ΔBCC on the [NaCl]. The plot would demonstrate (as our current Figs 2 and 4 now do) that in experiments Na does not change the OPs, but in the simulations it in most models does. Importantly, the ΔOP-vs-[NaCl] in a given model is connected to the ΔBCC-vs-[NaCl]: The models having small ΔBCC also have (like experiments) a vanishing ΔOP. Thus, Na does not bind the membranes.

17. It seems to me that these changes requires rewriting most of the manuscript. Was the idea that figures 2-4 would be replaced completely? These can be done but then it must be done effectively. If it takes too long, it is more reasonable to submit the paper as it is and then start to work with another one.

I am currently teaching and the only thing I am able to do before May is to submit the paper.

Would you have somekind of estimate how long it would take to make these modifications?

- In my opinion the electrometer concept is shown to be robust in reality (experiments). So we do not have to demonstrate that.
- I would not remove the density plots from the manuscript. In reality ions have continuum density distrbutions. All division we make to bound and non-bound ions are more or less artificial even though the conclusions would not depend on those. So I think that in the end one should always look at the density distributions.

18. Just a quick update. I have been steadily working on the plot, but it has taken some time. The current version only lacks the Slipids and Ulmschneiders'. Hope to finish it by tomorrow, after which I will write a more detailed reply.

19. I just read the MS once more, and I still like it a lot. I also honestly think that the new figure Markus is preparing can be included to the current version as one that wraps things up. And this can – and in my opinion should – be done without wasting time on reorganising the text.

Besides this, I would just add one paragraph that tries to explain why the three models (Charmm36, Orange and Lipid14) perform so well.

For Charmm36, it's due to the extra LJ interaction terms (NBFIX, Venable et al., DOI: 10.1021/jp401512z) that aim to push ions out from the bilayer. Such modifications might naturally lead to problems: Even though using NBFIX results in very little Na+ adsorption, the ester oxygens are preferred binding sites over phosphate oxygens, even though they lie deeper in the membrane (see DOI: 10.1039/c5cp05527j). Notable, in the Slipids model (and probably in many others), this behaviour is reversed and phosphate groups bind more Na+.

For Orange, Luca probably has an idea.

For Lipid14, we could probably compare the charges and LJ parameters of the head group to the ones in e.g. Charmm36 to get an idea. Or we could ask the developers if they have an idea.

These three models likely show correct behaviour due to very different approaches that could be pointed out and also briefly evaluated. I assume that in Orange and Lipid14 this results from the way charges and LJ parameters are derived, which is likely to cause less trouble than Charmm36's NBFIX.

And in case we do not include such discussion in the paper, I would still mention that Charmm36's goodness is due to the NBFIX, a solution that can be quickly adapted to any bad model to make it perform better in our benchmark.

20. I got all pieces in. Plot here. Comments welcome.

21. Hi,

I have couple comments. First, a tiny change in the figure. I would put all FF labels into one box perhaps like this: https://www.dropbox.com/s/5ig9uvn6eclswh8/plot1.png. It took me a while to find which color is which FF when the labels were separated. Other than that, I think the figure is really, really beautiful! I would also just add it to the manuscript without rewriting.

I absolutely agree with Samuli that we should do not remove the density plots at all.

I agree with Matti on adding the NBFIX comment in the manuscript. I don't think any of the other FFs use these kind of modified extra interaction terms.

Then one question. With mono and divalent ions the bound charge is localized differently in the membrane. It is true that electrometer concept will represent both but I don't think that a "single linear slope" can represent both mono- and divalent electrometer response. This will come evident in high enough concentrations where correlation effects of the bound ions start to matter (maybe it works pretty well in dilute concentrations). Should we add a sentence of this in the figure caption or is this immediately obvious to everyone? I know we're not claiming anything about this but after seeing the figure I just immediately got this urge to do a linear fit including both open and closed circles. :)

-Jukka

22. The figure is very nice indeed. I have one suggestion and some thoughts.

Suggestion:
I would put the y-axis in beta figure with the same scale as in alpha figure. Then the difference in the slope between the carbons would be visible. This might also make room for labeling changes suggested by Jukka.

Thoughts:
1) The differences between force fields for order parameter changes as a function of bound charge (this figure) do not seem to correlate with differences between force fields for order parameter changes as a function of ion concetration (Fig. 2 in the manuscript). I mean that the order parameter change as a function of bound charge is similar e.g. for Lipid14 and Berger even though significant difference is observed as a function of ion concentration. Another example is the Orange which has the smallest dependence on ion concentration but strongest dependence on bound ion. My conclusion is that the figure shows that in most cases the order parameters changes as a function of ion concentration arise from differences in ion binding efficiency, not from different headgroup
responses on ion binding.

2) The exception may be the alpha carbon for CHARMM. It has the weakest dependence on bound charge and it is the only segment which agrees with experiments as a function of Calcium concentration in Fig. 2. It may be that CHARMM is more realisitic in this aspect or this might coincidental.

3) I agree with Jukka that there is a strong temptation to make a linear fit for the lines. I think that this would lead to very relevant and interesting, but also difficult discussion. Thus, I think that we should leave it out from this publication. For this reason I think that we could leave the possible linearity/non-linearity unmentioned.

4) Markus, do you have a plan about embedding this into the manuscript? I do also support minimum modifications to the text but we cannot avoid some modifications.

23. This is reply to Matti's comment above (March 21, 2016 at 12:02 AM)

I see the point but, as I have mentioned before, I have been trying to keep the analysis of force field parameter differences outside of this manuscript to get it finished some day. But I totally agree that we should mention the NBFIX thing in CHARMM.

However, we do not have a result from CHARMM without NBFIX so we do not excatly know how large contribution it actually makes to the results. Interesting thing related to the this is that CHARMM has the most realistic response to CaCl_2. It is also good to remember that CHARMM has the most realistic headgroup structure which may or may not have influence on this.

Also the analysis of interactions affecting to binding in other models will be most likely very complicated. So I think that we should refrain from this analysis for this paper.

Overall, the situation seems to be that there is a lot of interest for further analysis of ion binding from the community. Despite of the relevance of the suggestions, I find myself all the time argumenting that we should not proceed with these analyses. One reason is that I want this paper submitted as soon as possible since my funding is about to end and I have also several other things to finish. If I get funding to continue this work, I will most likely suggest another publication to address issues discussed here. Before that, anyone who would have possibility to finish such publication is free to take the initiative.

24. I have now almost done with the teaching and I start to have time to do something else as well. I will start to work on this again and try to get this submitted asap.

I did already make the above suggested changes to the figure made by Markus:
https://github.com/NMRLipids/lipid_ionINTERACTION/blob/58e82e3de6d9d92197433dc2ab678f8de98aa261/scratch/boundIons/dOP_vs_boundCationCharge_P.pdf

In addition, I realized that we can also include experimental data into the plot from this paper: http://dx.doi.org/10.1021/bi00445a030

In this paper they measured order parameter changes as a function of charged amphiphiles in POPC bilayer. The advantage is that all the amphiphiles are in bilayer, thus they know the amount of charge in bilayer. From this data they have a relation between charge present in bilayer and order parameter change. I used this relation and plotted preliminary experimental lines in the figure. However, I had to do some assumptions and these have to be carefully described. Anyway, it seems that the response we get as a function of bound charge is very reasonable. This further supports the conclusion that, indeed, the binding affinity is the problem in simulations.

I will now embed these results in the text. I think that at the same time I will change the text towards the direction suggested by Markus. I will try to do this asap.

25. After reading the literature more carefully I realized that there is corresponding experimental data also for Ca2+ from this paper: http://dx.doi.org/10.1021/bi00312a019
They measure amount of bound Ca2+ per lipid with atomic absorption spectroscopy. They get linear relation between changed order parameters and bound Ca2+ per lipid. I did update their results in the figure, see
https://github.com/NMRLipids/lipid_ionINTERACTION/blob/3fd94be60d2990486a492ff7896200e4dbaf543d/scratch/boundIons/dOP_vs_boundCationCharge_P.pdf
However, they have bound ion per PC while we have bound charge density (q/nm^2). So to make the plot I had to assume the value for area per lipid (I assumed values between 0.58-0.78 nm^2/lipid). I think that it would be better that we would also report bound ions per lipid from simulations to avoid this assumption. To convert current data to this, we need area per lipid from simulations for all systems. It might be easiest that Markus would calculate these for this plot since he would now have all the trajectories used here in good location?

I will continue rewriting the results and discussion.

26. This experimental data is quite cool!

However, I am a bit worried that plotting it in the same plot with the simulation data might potentially be misleading. This is because when calculating the bound charge from the simulations we set the integration limit (e.g. until the P-maximum), and this might or might not match well with what is considered 'bound' in the experiment.

To put it in other words, we can effect the slopes of the simulation curves by our choice of the integration limit, whereas the experimental slope will not be affected. This has a potential of us introducing an unconscious bias on the one hand, and on the other hand people looking at the plot might conclude one FF to be better matching the experiments than some other, although this might not be the case (both because the OPs at zero bound charge in the simulations are not correct and because of this additional effect on the slope coming from the choice of the integration limit.)

27. I see the point. We have to figure out some other way to discuss the experimental data.

2. hello hello,
nice work guys. I read the ms and I have some questions / comments on the content.
The work shows problems for most force fields in reproducing the response to NaCl, but I think it falls short of identifying causes and consequences of those problems. Do you consider this out of the scope of the ms?
I think identifying causes is useful to people who want to improve force fields rationally (i.e., not shooting parameters at random until one set behaves ok), and clarifying consequences boosts the impact of the findings.
In terms of causes, I like the comparisons in the text, showing that ion FFs cannot be the only culprit. I'd go one step further and try to see where Na+ binds, in the various models; comparing parameters (charges, LJs) of those atoms might provide a hint to the causes of overly strong binding.
As for consequences, I was wondering (1) if there is a relationship between Na+ binding and the orientation of the P-N vector; (2) what is the response to NaCl in different models in terms of membrane dipole potential. I'm sure Markus and Samuli understand this better than I do, so what is your opinion? I think showing that simulation models succeed/fail in predicting changes in the membrane dipole potential might be interesting for a bunch of people (many consider the membrane dipole potential as the main reason why positively charged stuff sticks to membranes, while negatively charged stuff doesn't).

1. In addition to this comment, Luca send me a tex file with a lot of improvements to the text and some comments and questions. I did commit this file as it was
https://github.com/NMRLipids/lipid_ionINTERACTION/commit/ffdac093c7341c1ca4a2426f34caf21ce7cdb2a3
and then made some further remarks and edits:

More careful identification of causes and consequences would be useful indeed, however, I think it is better to do these separately from this manuscript.
There is already work done to this direction by Cordomi et al. (dx.doi.org/10.1021/ct9000763, http://dx.doi.org/10.1021/jp073897w)
and they conclude
"Despite the differences observed with different force field parameters, the lack of accurate enough experimental results prohibits discerning the best parameters to be used in this type of simulation. However, the present results
should be useful for selecting or calibrating new parameters when additional experimental information becomes available."
In our work we show that such accurate experimental results actually already exist.

Also the interactions leading to Na binding are somewhat discussed by Cordomi et al. I would be more than happy to see someone doing this more and extend the discussion to that direction. This would be easier now than ever because
the simulation data of this work is freely available in Zenodo and the comparison with experiments anchors the results to the reality.
However, I think that this quite tedious work and should be published separately. There is no advantage to wait this to happen since I think that the current manuscript is very important already with the current information content and more careful discussion about parameters can be easily published separately.

About consequences; (1) Also the correlation between Na+ binding and the P-N tilt in simulations is discussed by Cordomi et al. In addition, NMRlipids I paper showed correlation between order parameter change and P-N vector
tilt for all models with dehydration (to opposite direction than with ions, thought).
However, Na+ does not bind experimentally so the question is not relevant in this respect. According to the traditional structural models
the multivalent ion binding leads to the tilt in P-N vector. We could analyze this from MD, however, MD does not quantitatively
agree with experiments. Thus, we would not learn anything new from the analysis compared the traditional models since atomistic resolution details are not correct in MD. I think that the structural analysis is relevant only from the model which gives correct results, unless the analysis somehow supports the model development (and in this case it is beyond the scope).

(continues)

2. (continues)

(2) According to experiments (http://dx.doi.org/10.1016/S0006-3495(99)77414-X, ref [20]) 0.5 M NaCl does not affect the dipole potential, which is in line with the conclusions of our manuscript. In the same paper Ca2+ seems to have small but measurable effect, however, the comparison to this would not be direct since probes are used in experiments. Since the headgroup response to Ca2+ is not exactly correct in the models, also the quantitative dipole potential change is not expected to be correct. It would be somewhat interesting to see if the Ca2+ induced dipole potential change goes to the same directions as in the indirect experiments. However, it is not straghtforward to compare model with incorrect quantitative response to indirect experiment. It is also difficult to study the real reasons for the cation binding with models inaccurately reproducing the binding behaviour. Thus, I think that we should leave also these issues for further studies.

In conclusion, I think that the discussion and further studies on all the suggested topics cannot be avoided in the near future. However, this requires significant amount of work and we should not extend the scope of this work in the current situation. This would postpone the submission for several months or years, especially since I am currently too busy to start to push this extension. However, I think that we could shortly mention these issues and cite the work Cordomi et al. more carefully.

3. Samuli wrote: "However, we do not have a result from CHARMM without NBFIX so we do not excatly know how large contribution it actually makes to the results."

In fact, I have such trajectories that don't use NBFIX with NaCl. Since the resulting ion-density profile basically agrees with the profile obtained via Slipids, my guess is that omitting NBFIX in Charmm36 leads to similar behaviour as in the other force fields. I'll check the delta_Alpha and delta_Beta parameters and report them...

1. Excellent. You should be able to use directly CHARMM scripts in https://github.com/NMRLipids/lipid_ionINTERACTION/tree/master/scripts

4. Dear all,
Thanks for the nice piece of work, and the open discussion on the web !

We have done recently fluid charmm36DPPC+CaCl2 calculations using a recent model for calcium (http://dx.doi.org/10.1002/bip.22868).
The model is based on the assumption that the calcium
remains hydrated by water only. I found it interesting since it provides an original way to improve charmm ca2+ model within non-polarisable force-field, and also includes nice physical insight.

Here are the - (delta SCD) for the concentration of CaCl2
of 430mM and 890 mM. The errors are 0.004 for 430M, and 0.005 for 890mM.

C(mM) alpha1 alpha2 beta1 beta2
429 0.018 0.021 0.026 0.032
886 0.030 0.026 0.036 0.037

Here are the conlusions which I found interesting in the context of the present manuscript:
(1) I have the impression that this model performs well relatively to the other ca2+ models.
In the context of a lack of good calcium models, this seems noteworthy.
(2) the model impedes any direct interaction of ca2+ with lipidic heads, so one could expect
that it over-corrects the overbinding of Ca2+, but is is not the case ! Even with this assumption,
the order parameters are still affected by calcium binding.
(3) The assumption of complete Ca2+-water-hydration is strong and difficult to verify experimentally.
The comparison with NMR data is then very relevant to assess the model (the model was developped for DNA and not lipids, and the hydration assumption is not transferable). Given the results, I would not exclude that this model may indeed be relevant for lipid membrane/calcium interactions.

As the manuscript does not focus on the calcium, these results could probably be included in the manuscript with a minimum effort, without much change on the discussion.

As you want to submit quickly and not go too much in the direction of force-field comparison, I would simply emphasize that the model is better than the usual charmm Ca2+, but there is still space for improvement. One could also be slightly more moderate in the conclusion,
at the point when it is said that no good calcium model exists.

BR, Claire

ps : I wanted to calculate the amount of ca bound following the definition by Markus Miettinen, but I am not sure about the details. We have done the calculations using NAMD, so that using the scripts from github is not direct. As far as I understand, you take the total ion charge within the lipid bilayer, between the two planes defined by the peak maximum of the P density profiles ?

1. Thank you, this is very interesting.

I did include the results in Fig. 2 (https://github.com/NMRLipids/lipid_ionINTERACTION/blob/master/Fig/OrderParameterIONSchanges-eps-converted-to.pdf)

I think that the results are clearly better than "regular" CHARMM36. However, they are not within experimental errors. Interestingly, change in beta is overestimated (much less than in other models, thought) but change in alpha is underestimated (especially respect to the experiments for DPPC). This indicates that relation of changes between alpha and beta may not be correct in CHARMM (experimentally it should be \Delta\beta~-0.49\Delta\alpha [altenbach84]). The results in Fig. 3 indicate to similar direction (relatively mild slope for CHARMM alpha carbon).

In principle the current division of the data in the manuscript is that the "standard" models are ion the main body and then the "non-standard" models are in supplementary. Based on this division, these results would go to supplementary. However, we can think this more carefully since this seems to be the best one by this far (not good enough for quantitative interpretation of experiments, thought).

Anyway, we need the full simulation details and also the density profiles. I am especially interested on the length of simulations now since Matti Javanainen was worried about the equilibration times when commenting other post: http://nmrlipids.blogspot.fi/2015/02/the-first-draft-of-ion-lipid.html?showComment=1462793195771#c2861559608599412931

2. Concerning the simulation time :

We have done a total time of 200 ns for the hydrated calcium model.
The area per lipid is stabilized in less than 30 ns (from head, to be verified). As the interaction with the membrane is much less pronounced, it seems less critical to run long simulations. We can look at coordination numbers in the second hydration shell to test whether we reached a steady state.

For the standard ca2+ force field indeed, 200ns are not enough to stabilize even the area per lipid. For our simulations of DPPC at 323K after 160ns of simulation with less than 200mM CaCl2, we got a transition to a
gel phase, so we stopped at 200 ns.
Equilibrating the gel phase was out of our scope...

I shall provide very soon all the simulation details by adding a paragraph into the tex file on github.

3. Claire May 17, 2016 at 2:00 PM: I wanted to calculate the amount of ca bound following the definition by Markus Miettinen, but I am not sure about the details. We have done the calculations using NAMD, so that using the scripts from github is not direct. As far as I understand, you take the total ion charge within the lipid bilayer, between the two planes defined by the peak maximum of the P density profiles ?

We take the cation charge (so in this case the Ca++ charge) between the bilayer center and the peak maximum of the P density of the leaflet in question.

Note that both the bound cation charge as well as the change of the OP are calculated separately for each leaflet (thus the two lines in the plot). I did this to side step possible problems arising from the long relaxation times of ions. The headgroup should react to the bound charge on a time scale much shorter than the relaxation time of the ion binding, thus one can demonstrate the molecular electrometer concept even from these non-relaxed simulations.

PS Thanks for the nice data!

4. Sorry, I wrote "(thus the two lines in the plot)", but of course I meant the two sets of points per simulation. The two lines correspond to the two hydrogens.

5. I have now also written about these results into the discussion part. We have two issues here which arise mainly because these results are closer to the experiments than the other ones, so we have to be more careful in the conclusions.

1) Beta decrease is overestimated but alpha underestimated. Thus, it is not straightforward to say if the binding is too strong or weak. To support this discussion I added a figure (currently Fig. 3):
https://github.com/NMRLipids/lipid_ionINTERACTION/blob/master/Fig/OrderParameterAvsB-eps-converted-to.pdf

The figure shows the ratio of order parameter changes against experiments. As expected based on Fig. 2, the change in alpha is underestimated respect to beta in CHARMM (and also in almost all models). The problem is now, however, that we do not know which one of these (alpha or beta response to the bound cations) is more correct, thus we do not know which carbon would be better to use to compare binding affinity to experiments in this model. As written in the current version:
"This could be resolved by comparing CHARMM36 model to the experimental data with known amount of bound charge (e.g. experiments with amphihilic cations [31])"
If anyone has such simulation data and is willing to share it, please let us know.

2) According to the data from 2000ns simulations by Matti Javanainen, the 200ns simulation is not enough to conclude that there would not be more cation binding in the longer simulations, see:
https://nmrlipids.blogspot.com/2015/02/the-first-draft-of-ion-lipid.html?showComment=1463745379851#c5496492016333234005

If we would be able to conclude from these simulations that Ca binding is too strong also in this model, this would not be a problem. However, if the conclusion would be that the binding is too weak (or reasonable), this might be only because the simulation is too short to allow all Ca to bind. Anyway, if we cannot resolve the issue 1), this is not critical.

I think that we can leave both of these open and only mention these in the manuscript. However, the first one is quite easy to check of someone delivers the data so if this happens, I think we should do it.

6. I have a question about the simulations details of CHARMM36 simulation with Yoo model: Do you use the normal TIP3P or the version supposed to be used with CHARMM36?

7. Yes, we kept the charmm36 water model, which is a bit different from TIP3P. Except for the water in the calcium hydration shell, which is diffeerent.

We are on the way to characterize this water (permittivity, viscosity).

8. I have now updated this information in the manuscript.

One more thing: I think that we should add also the information about the reference simulation without CaCl_2 in the table with simulation details. Since you are using NAMD, this data might be useful also in NMRlipids III, see https://github.com/NMRLipids/NmrLipidsCholXray/issues/4

5. Hi !

One small comment on the lipid/Ca stochiometry remark.

I agree that there is a problem with these numbers, as mentionned in a note. I think we cannot compare the numbers 0.24 and 17.

The 1/4 given by Boeckmann is a number based on
the local structure of the binding around a calcium,
not on an analysis over the adsorption curve as in tatulian87. Tatluian et al give a stochiometry averaged over all lipids (ie per membrane surface) at saturation.
It takes into account the lipids which are directly bound but also the others. This average takes into account ion/ion repulsion for example. If one takes the same type of average ( (number of bounded cation/surface) * (surface/lipid) ),
Grubmueller and Boeckmann made a simulation of about 8 lipids per bounded calcium. But such a direct comparison of 8 and 17 also doesn't make sense, because the simulation by Grubmueller does not describe an equilibrium between a bulk solution and lipid membrane at saturation by the ions.

All their cations are bound (figure 4, there are no water/calcium interaction), therefore their salt concentration is below the saturation limit. Their lipid/calcium stochiometry of 8 is therefore a maximum estimation of stoechiometry at saturation.

Since the experimental 17 lipid/calcium is higher than the maximum estimation obtained by Grumeller et al, there is still a disagreement between the experimental results and the simulation results which justifies some attention.

But I am not sure I would go into these details at
this point in the paper, it would maybe be just confusing.

It may be enough to mention that
the calcium, as it is doubly charged, is expected to be much more difficult to handle by non-polarisable force fields, and that it is therefore meaningful to assess common available force fields in the light of the molecular electrometer concept.

1. Thanks, I agree with this. The discussion is now updated to this:
"Simulations suggest that Ca2+ bind to lipid carbonyl oxygens with coordination number of 4.2 [14], while interpretation of NMR and scattering experiments suggest that one Ca2+ interacts
mainly with choline groups [97–99] of two phospholipid molecules [30]. Simulation model correctly reproducing the order parameter changes would resolve the discussion by giving
atomistic resolution interpretation for the experiments."

In the cited NMR work they use binding model which takes into account the ion-ion repulsion and they get a good fit to the order parameter changes when they assume that each Ca2+ binds two lipids.

6. Hi ! one technical question. I am about to look at the density profiles for DPPC+CaCl2, and compare ours to the ones in the article. I am astonished that in the article, far away from the lipids, there are always similar amount of Ca2+ and Cl-. Even with the double layer effects, we always still have about twice as many Cl- as Ca2+. Or do I misunderstand the "atom number density". You mean the bare number of atoms, or the number of total electrons ? What are your weighting factors ? what is you unit ? atom per nm^3 ?

Thanks !

1. Hi,

Here is a recalculated density plot from the Orange+510mm CaCl data (Figure 6 Top, in the manuscript): https://www.dropbox.com/s/mx47kogdrgt9dxt/densitycacl.png I used the atom per nm^3 definition. The 2:1 ratio in the aq. phase is clearly seen here.

So in Fig 6 it seems that the Ca2+ densities are the actual number densities (atom/nm^3) but Cl- densities seem to have been scaled down by a factor of 2.

Do we need to scale the Cl- ions?

-Jukka

2. Note also that Jukka's curves appear qualitatively different (considerably less smeared) than what is seen in the current Fig. 6. I suppose this because of the issue with the centering of the bilayer, right?

3. The units are indeed atom per nm^3. I mean bare atom number density.

For all the components number density is calculated. Then for the plotting these are scaled to make the plot easier to understand. As written in the caption "The lipid densities are scaled by 100 (united atom) or 200 (all atom model)". What is NOT written yet in the caption is that Cl densities are scaled with two (sorry about that, I will add this asap). Scaling is done in the gnuplot script which generates the figure: https://github.com/NMRLipids/lipid_ionINTERACTION/blob/master/Fig/scripts/plotDENSITIEScaCLEAR.gpl

So the density calculation is normal, scaling is done when plotting. So you can deliver the density profiles from bare atom number density.

4. Claire Loison sent me density profiles my email and I have now added those in the manuscript:
https://github.com/NMRLipids/lipid_ionINTERACTION/commit/b32552d2b7b7822a94707772b47861097f51ae29

7. Just 2 minimal typos I let you correct IF they are relevant :

- in Fig. 2, you could add the temperature of the simulation for Yoo (323K).

- I may misunderstand, but in the text

"More specifically, the experimental data shown in Fig.~\ref{AvsB} led to the
relation $\Delta S_{\rm{CH}}^{\beta}=0.043\Delta S_{\rm{CH}}^{\alpha}$ for DPPC bilayer with various CaCl$_2$ concentrations~\cite{akutsu81}."

Is there not a zero too much ? Is it not $\Delta S_{\rm{CH}}^{\beta}=0.43\Delta S_{\rm{CH}}^{\alpha}$ ?

1. I have modified the figure to make it more clear such that the simulation temperatures are not mentioned anymore.

Indeed, it should be 0.43. This is now fixed.

As it in principle would be possible that the binding affinities of cations would be correct in simulations, but the headgroup just would be too sensitive to bound charge, it would be good to have a test system with positively charged surfactants / lipids in a PC bilayer. Comparing the changes of headgroup OPs in this system against those upon addition of cations should allow one to check if the headgroup sensitivity charge is at least ball-parkish correct.

To this end I have added to Zenodo some trajectories of cationic DMTAP/DMPC bilayers:

0/100 mol%

6/94 mol%

50/50 mol%

75/25 mol%

The order parameters changes in these could be compared to those observed in experiments, e.g., by Altenbach et al. (1984) and Scherer et al. (1989). If they match in order of magnitude, I think we can exclude the Over Sensitive Headgroup hypothesis.

1. Apparently the direct links to Zenodo do not work. Here the DOI links:

0/100 mol%
6/94 mol%
50/50 mol%
25/75 mol%

2. I have now added the data from the trajectories and some discussion in Supplementary Information. The results are also referred in the current version of Results and Discussion.