Tuesday, September 10, 2013

The manuscript is now available. What next?

This blog has been started to make improvements on the manuscript titled "Response of the hydrophilic part of lipid membranes to changing conditions — a critical comparison of simulations to experiments", written by O. H. Samuli Ollila, and openly available at http://arxiv.org/abs/1309.2131.

Here is the abstract of the manuscript: "We compare the order parameters predicted for the hydrocarbon segments in lipid bilayer headgroup region by the Berger molecular dynamics simulation model to those measured by Nuclear Magnetic Resonance (NMR) experiments. We first show results for a fully hydrated POPC bilayer, and then focus on changes of the order parameters as a function of hydration level, NaCl and CaCl2 concentrations, and cholesterol content. The experimental headgroup order parameters are never reproduced. This indicates that under all of these conditions the used model is unable to correctly reproduce the headgroup structure. Consequently, many of the conclusions drawn over the years from this model might be erroneous."

Although I think that the manuscript is a high quality scientific document, it has not been submitted to any journal, instead its contents are discussed here. The reasons for this decision you can find in the post Why not submit to a journal?

The ultimate goal of this blog is to find an atomistic (preferably united-atom) force field that reproduces the experimental properties discussed in the manuscript. Naturally the optimal situation would be that some of the already available force fields would fulfill this goal. If this, however, turns out not to be the case, the goal will be to find the appropriate modifications.

All the work will be progressed openly in this blog, and anyone is free to join in. In the end of the project the manuscript, updated with the improvements made via this blog, will be submitted to an appropriate scientific journal. Authorship will be shared among the blog contributors, as described in the post On credits.

My current ideas about things to be done in order to achieve the goal are described below. Some of these issues are straighforward to do, especially for the people who already have done something similar. For example, if you have done CHARMM simulations on a lipid bilayer with ions it is relatively easy to calculate headgroup order parameters and share the results through this blog. Some of the issues are not that straighforward. Below I have listed some tasks which would improve the manuscripts significantly. In addition of these issues, any other comments, questions and contributions are very welcome.

Things to do:

- The CHARMM [Klauda et al. JPCB 114, 7830 (2010)] and the GAFF [Dickson et al. Soft Matter 8, 9617 (2012)] force fields give better order parameters for the headgroup under full hydration. Interesting question now is, would these force fields also reproduce the correct changes in the order parameters as a function of dehydration, ion cholesterol concentration? At least for Sodium [Valley et al. J. Membrane Biol. 244, 35 (2011)] and cholesterol [Lim et al. JPCB 116, 203 (2012)] simulations have already been done with the recent CHARMM model. The calculations of headgroup and glycerol order parameters from these simulations would be very informative.

- For the all-atom Stockholm lipids force field [Jämbeck et al. JPCB 116, 3164 (2012)] I have not seen order parameters calculated for headgroup or glycerol. These would be interesting to see.

- I have not seen the headgroup and glycerol order parameters for the united atom lipid force fields based on GROMOS by Kukol [JCTC 5, 615 (2009)] and by Chiu et al. [JPCB 113, 2748–2763 (2009)]. Poger et al. [J Comput Chem 31: 1117–1125, (2010)] comment the order parameters briefly, but more information would be interesting. Fore more details and discussion see the post on Some points concerning the POPC full-hydration results.

- If it turns out that none of the existing united-atom force fields can reproduce the experimental order parameters under the studied conditions, we should find the way to fix the models. This is probably not a straightforward task. I have presented some of my ideas and opinions in the post on Some points concerning the POPC full-hydration results.

- The ion parameters play naturally a large role in the interaction between ions and lipids, making this a particularly complicated issue. I have written some ideas in the post on Some points concerning the ion-membrane interaction results.

All the files which are used to run the simulations and analysis for the manuscript can be downloaded from: http://dx.doi.org/10.6084/m9.figshare.790722

Everyone is invited to discuss, criticize and contribute to the manuscript through this blog. Communication and contributions can be done by writing the comments to the blog posts. If you pick one of the already determined tasks, it might be good to inform readers to avoid overlapping contributions.

Samuli Ollila
Aalto University


  1. This is interesting idea to test a new way of publishing.

    My first question, after looking on the figures, what is uncertainty of the experimental (and also simulation) results? Can we be sure that they are, for example, within 0.02 for the order parameters?

    Another question is sign of the order parameters. In some cases experimental order parameters for lipid headgroups are reported with sing (I have in hand references for DMPC, JACS 119, 796 (1997) and Langmuir 19, 10468 (2003), may be there exist such data for POPC), this would be even a stronger test on suitability of a force field.

    1. Thank you for the first comment!

      I think that the experiments are in principle very accurate (more accurate than 0.02). However, I have understood that there might be technical issues, for example mild heating due to the pulse sequence etc. which might cause some deviation between results reported in different papers. However, the results measured with different authors and different techniques show very little variation. See for example, Fig. 3 and Table 1 in http://dx.doi.org/10.1039/C2CP42738A and dx.doi.org/10.1016/j.bbamem.2010.11.027. I would be happy to hear a comment from NMR expert on this.

      In simulations, I am not sure which is really the best way to analyze the statistical error, but I think that it is always less than 0.02, no matter how it is calculated.

      About the sign: I am not aware of measurement for POPC. However, it seems that the headgroup order parameters are not affected by the changes in the chain region (see the table 1 and discussion in http://dx.doi.org/10.1039/C2CP42738A). Thus I assume that the signs are the same as in DMPC. The signs are actually correct in Berger. Unfortunatelty, I have forgot to mention this in the manuscript. I have shown the absolute values for clarity.

    2. In conclusion, I would say that 0.02 would be a large enough error bar for both, simulations and experiments.

  2. Here are the results I calculated for three bilayer systems with the Slipids force field:

    1) Pure POPC, 512 lipids, 298K, 40 waters/lipid, 150mM NaCl, PME, 300ns, Slipids
    2) POPC+10%Chol, 288 PCs, 310K, 45 waters/lipid, 150mM NaCl, PME, 200ns, Slipids
    3) DPPC+10%Chol, simulation parameters as above

    The plot is available at https://dl.dropboxusercontent.com/u/6808285/hg_order.pdf

    and here's the raw data:


    1 0.0239 0.0596
    2 0.0645 0.0464
    3 0.0268 0.0328
    4 0.0397
    5 0.1023 0.1622


    1 0.0212 0.0525
    2 0.0569 0.0452
    3 0.0035 0.0059
    4 0.0159
    5 0.0813 0.1603


    1 0.0539 0.0209
    2 0.0513 0.0350
    3 0.0396 0.0363
    4 0.0527
    5 0.1114 0.1627

    Where the numbers stand for 1-beta, 2-alpha, 3-g3, 4-g2, 5-g1.
    I have simulations with Kukol's POPC model and I will post these results soon.

    1. Thank you for the first contribution!

      I plotted your results for the POPC with 150mM NaCl (system 1)) together with the result from Berger model with 200mM NaCl and the experimental result for POPC without NaCl. Experimentally these concentrations of NaCl does not change the order parameters so the comparison is reasonable. The plot is here: https://www.dropbox.com/s/8xu6ucacwtj407n/OPcompBERSLIPnacl.png

      The major difference seems to be the low order parameters for g1 and g2 with Slipids.

      For the other systems I do not have direct comparison now since in the 2) there are both cholesterol and NaCl, and in 3) there is DPPC. However, there might be some older experimental data for these.

  3. And here are results for a pure POPC simulation with Kukol's UA model based on GROMOS. The system consisted of 512 POPC molecules with 40 waters/lipid and 150 mM of NaCl. Simulation temperature of 298 K was employed. Two simulations of 300 ns were performed with either PME or reaction field (RF) responsible for handling the long-range electrostatics. A triple-range cutoff was employed for LJ and RF interactions, while with PME a normal cutoff for real-space calculations was used.

    I plotted the results here: https://dl.dropboxusercontent.com/u/6808285/kukol_popc.pdf

    And the data in ascii form is as follows:

    1) PME

    1 0.0563 0.0012
    2 0.0991 0.1293
    3 0.2811 0.1319
    4 0.1643
    5 0.0463 0.2033

    2) Reaction Field

    1 0.0376 0.0107
    2 0.0838 0.1171
    3 0.2682 0.1210
    4 0.1610
    5 0.0350 0.2063

    where the numbers are as follows: 1-beta, 2-alpha, 3-g3, 4-g2, 5-g1.

    1. Thanks again!

      I have now updated the results into the figure: https://www.dropbox.com/s/8xu6ucacwtj407n/OPcompBERSLIPnacl.png

      By the way, are these simulation data sets in part of some publication? Have you checked if Na ions partition in the bilayer in a similar way as in Berger model?

  4. The Slipids results are obtained from simulations that were performed to support some experimental findings yet it is uncertain whether they will be published in the near future.

    The results with Kukol's model were supposed to serve as a reference system to on smaller study. However, whether or not these will be published, is currently unknown.

    I haven't really studied the behaviour of ions in these systems but I will take a look.

  5. So I checked the partitioning of sodium ions in both Slipids and Kukol's model and it is clear that with Kukol's model the Na ions reside in the carbonyl region much more than with Slipids. So I guess that the behaviour of Kukol's parameterisation resembles that of Berger.

    1. Interesting. Do you see the maximum in ion density in carbonyl region also with Slipids, or is it just smaller than with Kukol's model?

    2. There is a peak in the carbonyl region but even at this peak the density is smaller than in the bulk water. With Kukol's model the density in carbonyl region is higher than in bulk water.

  6. Hi!

    I calculated together with Matti the order parameters for the OPLS-AA POPC model we have used extensively in the group. The bilayers consist of either 304 (pure and 10% chol systems) or 306 (40% chol system) POPC molecules with 0,10 or 40% cholesterol in addition to that for a total of 304, 338 or 512 molecules. The systems were simulated at 310 K with extensive hydration and 150 mM NaCl. The systems were simulated for 200 ns (0% or 10% cholesterol) or 300 ns (40% cholesterol). PME electrostatics were employed. The results are as follows:

    No Cholesterol

    1 0.0382 0.0018
    2 0.1133 0.0711
    3 0.2777 0.1901
    4 0.1805
    5 0.1792 0.1836

    10% Cholesterol

    1 0.0327 0.0017
    2 0.1065 0.0741
    3 0.2742 0.1929
    4 0.1790
    5 0.1779 0.1842

    40% Cholesterol

    1 0.0188 0.0075
    2 0.0995 0.0782
    3 0.2540 0.1970
    4 0.1865
    5 0.1908 0.2021

    Carbons numbering is as earlier in the results Matti reported.

    1. Thanks!

      I plotted your results together with the results from Berger/Höltje: https://www.dropbox.com/s/pazwypsocmo389v/OrderParameterCHOLwTYNK.png

      With a first look it seems that your results are changing less with cholesterol than in Berger. However there are less points and the results for largest concentrations are only for Berger, so I am not sure.

      Also the comparison here is not straighforward since you have NaCl in the system, while in the experimental results and the Berger simulations there is not. With pure POPC we know that NaCl does not change order parameters in experiments, however I cannot recall experiments with PClipid/cholester/ion mixtures. So in principle it could be possible that ions would change the headgroup order parameters when cholesterol is present.

      There does not seem to be improvement for pure POPC compared to Berger.

      By the way, could you specify which force field was this? Are there some different versions of OPLS-AA?

    2. I have now added the result for POPC with 150mM NaCl simulated with OPLS AA and reported by Joona Tynkkynen into the same figure with the results for Berger, Slipids and Kukol:

      I am not sure about the exact version of OPLS AA force field, though.

  7. And here is one more parameter set with an all-atom POPC model currently under development by Arek Maciejewski and Tomasz Rog in our group. This force will soon be published yet these results on head group order parameters look very promising.

    This is a simulation of a POPC bilayer (288 lipids) with 44 waters/lipid at 310 K, PME electrostatics. I used the last 80 ns of a simulation with total duration of 100 ns for the analysis.

    1 0.0310 0.0201
    2 0.0494 0.0483
    3 0.2328 0.1380
    4 0.2020
    5 0.1560 0.0137

    1. Thanks!

      I plotted it together with Berger: https://www.dropbox.com/s/bpkb8i31q34tu1i/HGorderparameterswMacRog.png

      I think there is a significant improvement. Only large difference compared to the experimental values is the splitting in g3. There also seems to be small splitting in beta, but that is quite small.

      I would be very happy to see how this model would reply, for example, to the dehydration or other changes discussed in the manuscript.

  8. Here are results with Berger DPPC and a couple NaCl force fields. The results are of hydrogens attached to carbons labeled as: gamma, gamma, gamma, beta, alpha, g3, g2, g1, respectively. Electrostatics is handled with PME, 100ns, 323K, 72 lipids.

    1) Pure DPPC

    0 0.023526 0.006980 0.012376
    0 0.023588 0.008533 0.009391
    0 0.023229 0.014370 0.009469
    1 0.028587 0.047187
    2 0.098500 0.148678
    3 0.234071 0.201362
    4 0.213305
    5 0.057871 0.065860

    2) DPPC with 150 mM NaCl, Joung and Cheatham ions. The ions are parameterized using Lorenz-Berthelot LJ combination rules and since the simulation uses the standard geometric rules with Berger lipids, mixing combination rules most likely is the cause of the bad results.

    0 0.010982 0.007178 0.004313
    0 0.189192 0.190906 0.222319
    0 0.144293 0.238756 0.331118
    1 0.322813 0.217390
    2 0.224558 0.368983
    3 0.228063 0.345894
    4 0.110900
    5 0.102004 0.154022

    3) DPPC with 150 mM NaCl, Åqvist ions. These ions are the standard ions used with Berger lipids. As pointed out earlier, in this case the sodium ions indeed bind into the phosphate/carbonyl parts (perhaps too strongly).

    0 0.029240 0.014293 0.012808
    0 0.031056 0.009083 0.018262
    0 0.029880 0.013904 0.006406
    1 0.025345 0.026001
    2 0.063802 0.132597
    3 0.250173 0.194344
    4 0.224224
    5 0.111435 0.082425

    4) DPPC with 150 mM NaCl, Åqvist ions with scaled charges q=0.7 for both Na and Cl. The scaling is like Leontyev et al. PCCP 13, 2613 (2011) suggested. Now the order parameters are close to those of pure DPPC. Also, now the ions don't partition into the membrane, but instead reside in water. While the ions now behave 'better', now we have a problem where our partial charges are not integers for ions. This might cause serious issues especially with charged molecules.

    0 0.028271 0.010696 0.014559
    0 0.028739 0.011787 0.010310
    0 0.028564 0.013896 0.010641
    1 0.027294 0.036661
    2 0.084137 0.131134
    3 0.229839 0.184166
    4 0.224365
    5 0.094302 0.073641

    1. Thanks for the results!

      1) I plotted your results for the pure DPPC together with my results for POPC and experimental values:
      Experiments for POPC are from Ferreira et al. PCCP 15, 1976 (2013) and for DPPC from Seelig et al. BBA 467, 109 (1977) (alpha, beta, g3) or Gally et al. (g1,g2) Biochemistry 20, 1826 (1981). Thus g1 and g2 values are actually from E. coli extract. However, I think that the values are comparable since typically very small differences are measured between different lipids for these values.

      Experimental values for DPPC are slightly higher, probably due to the higher temperature (322K) compared to POPC (298K). In the simulations there are larger differences, expecially for g3 and g1. I believe that this is also due to the temperature difference in simulations (298K in POPC, 323K in DPPC). I will check this by running the POPC simulation in 323K.

    2. I did plot Jukka's results for ions in similar plot to the Fig. 4 a) in the manuscript:

      1) As you already explained, the Joung and Cheatham ions gives strange results: https://www.dropbox.com/s/icpu0zwxwur1g9w/OrderParameterIONSnaclJOUCHEAT.png

      2) The changes in the order parameters with Åqvist parameters with DPPC seems to give similar results to mine with POPC:
      The changes are slightly smaller though. By the way can you give a specific reference for the ion parameters?

      3) Interestingly, also the ions with scaled charges seems to give similar changes in order parameters:
      The changes are again smaller but I would say visible. This is surprising since you said that the ions did not penetrate in the bilayer. Based on this it looks like that the ions change the headgroup order parameters in this model even though they do not penetrate into bilayer. I think that this a little bit strange result.

      Note that I have plotted the "experimental" results for glycerol groups with empty circles here since I realized that those are not actually measured as a function of NaCl concenration in the original experimental reference (Akutsu et al. Biochemistry 20, 7366 (1981)). However, I find it very unlikely that they would change since there is no change in the g3 and the chain C2 order parameter as a function of CaCl concetration, even though alpha and beta are changing. With NaCl alpha and beta are unchanged, so I think that it would be very unlikely that glycerol parameters would change in this case.

    3. In the above comment I suspected that the difference between the results for DPPC (reported by Jukka Määttä) and POPC (from the manuscript) with Berger model would arise from the temperature difference. To check this I ran a POPC simulation in T=323K. The temperature seems to have very little effect: https://www.dropbox.com/s/sa8zeeu7t5tdue4/OPtempCOMP.png

      So either Berger gives different results for DPPC and POPC or there is some difference in *itp files used by myself and Jukka Määttä. Jukka: Can you share the *itp file you used or compare it to the one I have shared through the fixshare?

    4. I uploaded an old Berger DMPC bilayer trajectory (used in http://dx.doi.org/10.1021/jp810233q) to figshare:


      Maybe analyzing that would help to see if the tails affect the headgroup order parameters in Berger? Temperature was 323 K.

    5. I calculated the order parameters from this data. The files I used for protonation and order parameter calculation are here: https://www.dropbox.com/s/i72poh3r8ghi28i/markusFILES.tar

      The result plotted together with previous results for POPC and DPPC with Berger are here:

      I also included experimental order parameters for DMPC taken from here: Gally et al. JACS 119, 796 (1997) http://dx.doi.org/10.1021/ja962951b

      It seems that in the Berger model the order parameters somewhat depend on the tail composition which is not the case in the experiments.

    6. In the original posting only the first 75 of the 128 lipids were included in the analysis, because of a bug in the analysis script. Samuli Ollila provided (one-to-one communication) the following data with the fixed script, now also including the signs, on 10 April 2014:

      1 0.0510405
      1 0.0643061
      2 0.103433
      2 0.153471
      3 -0.177112
      3 -0.251111
      4 -0.166239
      5 0.138638
      5 0.0849596

  9. Hi,

    Thanks for plotting the data, the way you plotted the parameters is more informative than my plots.

    Regarding your comments:

    The specific reference for Åqvist ions: J. Åqvist, J. Phys. Chem., 1990, 94 (21), pp 8021–8024

    The changes being slightly smaller with DPPC than with POPC is most likely because of the different temperatures. In my simulations the temperature was higher (323K) and so the ions most likely are slightly more mobile because of thermal motion.

    The fact that the scaled ions also effect the bilayer order parameters is indeed surprising. Thanks for pointing that out. Here are the density plots for the NA ions:


    One can see that the standard Åqvist ions peak near the bilayer while the scaled ions reside in water. I agree that it's quite strange that they still change the headgroup properties.

    Also, here are the average areas per lipid from the simulations:

    61.1 +- 1.8 Å^2 for pure DPPC, 60.1 +- 1.4 Å^2 for DPPC with Åqvist ions and 60.8 +- 1.6 Å^2 with scaled Åqvist ions. Based on these values I'd say the areas are unaffected by the ions, but I haven't looked any experimental data on this.

    1. This is a bit of an off-topic question, but just for the sake of curiosity: Did you happen to look what the headgroup orientations (P–N vector angle with the membrane normal) are for the scaled and the normal ions?

  10. Here are order parameter results for Slipids DPPC with a really short run (30 ns).

    0 0.008918 0.009518 0.007352
    0 0.009103 0.008833 0.005143
    0 0.012708 0.000742 0.009167
    1 0.023532 0.055842
    2 0.059135 0.053818
    3 0.015921 0.050389
    4 0.022365
    5 0.111682 0.178645

    I haven't looked at the P-N vectors yet, I'll try to post the data later as it would be interesting to see how they behave.

    1. I plotted Jukka's Berger and Slipid results for DPPC together with experimental values:

      The results are in line with previous results for POPC and ions simulatied with Slipids (reported by Matti Javanainen):
      The order parameters are quite low for the glycerol group in Slipids.


Please sign in before writing your comment.