Thursday, March 9, 2017

NMRlipids IV: Headgroup & glycerol backbone structures, and cation binding in bilayers with PE, PG and PS lipids

In NMRlipids I and II projects, the goal was to find a MD model that would correctly reproduce NMR data (for lipid headgroup & glycerol backbone structures, and for cation binding) in PC bilayers. In NMRlipids IV project, we set the same goal for PE, PG and PS lipids in bilayers (pure or mixed with PC). The standard NMRlipids workflow and rules will be applied. The current version of the manuscript is available in the GitHub repository.

Currently, the manuscript is mainly a collection of relevant experimental data. For example, Fig. 1 compares the experimental headgroup and glycerol backbone order parameters between PC, PE, PG and PS lipids.
Fig.1 Absolute values of order parameters for headgroup and glycerol
backbone with different headgroups from experiments. For references and other details see the manuscript.
The conclusion based on this, together with some additional data, has been that the headgroup structures are similar for PC, PE and PG lipids, while PS headgroup is more rigid [Wohlgemuth et al, Buldt et al.]. On the other hand, the glycerol backbone structure has been considered to be similar in model systems and cells for all these lipids [Gally et al.].

Some preliminary comparison between experiments and simulations with CHARMM GUI parameters are shown in Figs. 2 and 3, suggesting that the model has some difficulties to reproduce the experimental order parameters for PS and PG headgroups. More detailed conclusions are difficult to draw only from these data, because experimentally the signs of order parameters for PS and PG are not available (as far as I know). However, the results from other models might help to draw some connections between order parameters and structural details, as was done in NMRlipids I for PC lipids.
Fig 2.  Order parameters for POPS headgroup and glycerol backbone
from simulations and experiments. For references and details see the manuscript. Absolute
values are shown for experimental data, because signs are not
known. Simulations values are -SCH

Fig 3. Order parameters for PG headgroup and glycerol backbone
from simulations and experiments. For references and details see the manuscript. Absolute values are shown, because signs are not known for experimental data.
Experimental data on cation binding in PC bilayers mixed with of negatively charged PG and PS lipids is shown in Fig. 6. As expected, adding CaCl2 causes a stronger decrease in the PC headgroup order parameters when the amount of negatively charged lipids is increased. According to the NMR electrometer concept (see NMRlipids II for discussion), this means that the amount of bound Ca2+ increases when negatively charged lipids are present in bilayers.
Fig. 6 Changes in the PC headgroup order parameter as a function of CaCl2 concentration in bilayers containing various amounts of negatively charged lipids. For references and details see the manuscript.
A more specific interpretation of this kind of data has been that [Seelig]:
"(i) Ca2+ binds to neutral lipids (phosphatidylcholine, phosphatidylethanolamine) and negatively charged lipids (phosphatidylglycerol) with approximately the same binding constant of K = 10-20 M-1;
(ii) the free Ca2+ concentration at the membrane interface is distinctly enhanced if the membrane carries a negative surface charge, either due to protein or to lipid;
(iii) increased interfacial Ca2+ also means increased amounts of bound Ca2+ at neutral and charged lipids;
(iv) the actual binding step can be described by a Langmuir adsorption isotherm with a 1_lipid_:_1_Ca2+ stoichiometry, provided the interfacial concentration CM is used to describe the chemical binding equilibrium."

I believe that an MD simulation model correctly reproducing the cation binding in negatively charged lipids could further sharpen this interpretation.
The goal of this project will be to test if currently available models can be used for an such interpretation. This should also help the model development (if needed), however the actual improvement of force fields is beyond the scope of NMRlipids IV.

As in all NMRlipids projects, all types of contributions (data, comments, criticism, etc.) are welcomed from everyone. The authorship of the publication will be offered to all contributors and the final acceptance is based on self-assessment according to NMRlipids rules. The following contributions would be especially relevant at this stage:
  1. Results from different simulation models. Simulations of bilayers containing PE, PG or PS almost under any conditions would be currently useful to map the behavior of different models. Direct delivery of calculated order parameters through GitHub or blog comments, or by making the simulation trajectories accessible (for example, through Zenodo) would be ideal ways of contributing. 
  2. Order parameter signs for PE, PG and PS. The order parameter signs are very important for the structural interpretation. However, I am not aware of the order parameter sign measurements for other than PC lipids. If such data would be somehow available, this would be highly useful contribution.


  1. Hi Samuli,

    I have lots of pure POPS and DOPS simulations that you can use. There are 4/5 force fields (Berger, 2 closely related GROMOS ones, CHARMM36 and Slipids) with repeat simulations and several different commonly used cut-off schemes for these different force fields. All simulations are 500 ns using GROMACS 5.0.6 and have been simulated from the same starting structure. The membranes have 128 lipids and the systems just have neutralising Na+ ions.

    I am busy with several things at my main work, plus I am also trying to finish off a paper at the moment (in fact also related to order parameter calculations; I will share more on this when I can) but when I get the chance I will calculate all of the order parameters for these POPS/DOPS simulations and share them.

    I also should have some PE and PG simulations too with a few different force fields (definitely with some GROMOS ff's, probably also some CHARMM36). I will look out what I have on these from my archived data. Again, these will probably just have neutralising Na+ ions (for PG) or no ions at all (for PE).



    1. Excellent!

      Are the mentioned GROMOS related parameters the ones which can be found also from Lipidbook?

    2. The GROMOS PS parameters aren't yet on Lipidbook but follow the same strategy as the rest of the GROMOS-CKP lipids (the larger vdW radii for the carbonyl carbons as per the Chandrasekhar and Kukol approaches with the rest of the parameters as standard GROMOS). There are, in fact, two slightly different parameter sets for PE and PS. Ones with the same charges as those used in the Berger PE/PS lipids and ones with standard GROMOS charges for the NH3 groups (as taken from a lysine side-chain). For PE, I also believe I have simulations with GROMOS 43A1-S3, standard Berger plus some Berger sims with modifications to include some repulsive LJ interactions of the NH3 hydrogens, OPLS-UA (with the same types of modifications as per some of the Berger sims) and maybe some others too. The GROMOS PG simulations I have will be with the GROMOS-CKP parameters already on Lipidbook.

      I should hopefully be able to find everything soon related to PE/PG in my old files. I already have all the PS simulations to hand and will try and run the analysis asap for these.

    3. By the way, which would be the main reference for the GROMOS-CKP parametrization strategy?

    4. The original approach was proposed for DPPC by Chandrasekhar et al. (the C in the CKP; where the increase to the van der Waals radii for the ester carbonly carbon was first used. Kukol (the K in CKP) took this approach further by applying it to other PC and PG membranes ( However, upon testing I found that there were several problems with the Kukol parameters (so with everything apart from DPPC IIRC). In particular the problems were with (i) the bonded parameters in the head group plus glycerol regions and (ii) the use of the Bachar double bond dihedrals which gave really bad order parameters (although these dihedrals do work well for Berger POPC). These problems are discussed briefly by myself (the P of CKP) in the SI of and are discussed in more detail in

      So, unfortunately a bit of a mess in terms of referencing. I guess for general parameterisation strategy there should be a reference to the work of Chandrasekhar et al. and probably the first of my two papers (the first one doesn't mention the name GROMOS-CKP explicitly but does show the approach working well for PE, PG and Cardiolipin membranes; the second shows the approach working properly for POPC and demonstrates the problems with the Kukol POPC parameters).

  2. I forgot CHARMM36-UA POPS/DOPS simulations too!

  3. So slowly processing through head group PS order parameters.

    Firstly reported below are for POPS with CHARMM36 (four simulations, 2 with 0.8nm vdW switching point and 2 with 1.0 nm). The analysis was performed on the final 100 ns from 500 ns simulations of a 128 lipid bilayer with Na counterions performed using GROMACS 5.0.6. The results show in order beta, alpha, G3, G2, G1:

    POPS 0.8 nm switching version 1

    1 -0.0492503
    2 0.017474 -0.0574774
    3 -0.219257 -0.199918
    4 -0.229346
    5 -0.198932 -0.0169096

    POPS 0.8 nm switching version 2

    1 -0.0563844
    2 0.0152987 -0.0618646
    3 -0.204029 -0.236506
    4 -0.197588
    5 -0.191738 -0.0522408

    POPS 1.0 nm switching version 1

    1 -0.0665649
    2 0.0130353 -0.0624143
    3 -0.225571 -0.241592
    4 -0.231893
    5 -0.216094 -0.0510908

    POPS 1.0 nm switching version 2

    1 -0.0803705
    2 0.00809381 -0.0760739
    3 -0.214749 -0.229649
    4 -0.226479
    5 -0.226623 -0.0313807

    1. Could you report also the amount of water molecules in the systems? We are planning the experiments to measure the signs and it would be convenient to use the same amount of water.

    2. Sorry for the slow reply, I hadn't seen this until now. There are 35 waters per lipid in all the systems (so 4480 waters).

  4. Next, exactly the same for DOPS (I should have said above that the POPS simulations were performed at 298 K, these DOPS ones were at 303 K):

    DOPS 0.8 nm switching version 1

    1 -0.0443464
    2 0.00983792 -0.0443395
    3 -0.217863 -0.236201
    4 -0.206458
    5 -0.186183 -0.0443288

    DOPS 0.8 nm switching version 2

    1 -0.0368902
    2 -0.00558114 -0.0368946
    3 -0.223083 -0.227814
    4 -0.221595
    5 -0.20654 -0.0360436

    DOPS 1.0 nm switching version 1

    1 -0.0445504
    2 0.00475181 -0.0456792
    3 -0.21634 -0.229123
    4 -0.195551
    5 -0.191639 -0.0342016

    DOPS 1.0 nm switching version 2

    1 -0.0601911
    2 0.00238097 -0.0424167
    3 -0.229217 -0.23201
    4 -0.217005
    5 -0.189629 -0.0351633

    Analysis of the other force fields is in progress. I'll report these when completed.

  5. Slipids next. This time there are six simulations for each of the POPS and DOPS membranes. Two repeats with three different cut-off settings for the vdW (15A, 10A, 10A+LJ-PME) Same systems and other settings as the CHARMM36 membranes and the same analysis (last 100 ns of 500 ns simulations). POPS first:

    POPS 10A cut-off version 1

    1 -0.0915155
    2 0.0220854 -0.0936563
    3 -0.00611644 0.131848
    4 0.00297916
    5 0.0397122 -0.159115

    POPS 10A cut-off version 2

    1 -0.0815755
    2 0.0188775 -0.0806864
    3 -0.00707084 0.0987901
    4 -0.00247234
    5 0.0283463 -0.159617

    POPS 15A cut-off version 1

    1 -0.099784
    2 0.0174399 -0.0957243
    3 -0.0149656 0.130198
    4 -0.00562752
    5 0.044824 -0.16667

    POPS 15A cut-off version 2

    1 -0.103846
    2 0.025198 -0.108365
    3 -0.0281673 0.140332
    4 -0.0212674
    5 0.0414582 -0.158662

    POPS 10A + LJ-PME cut-off version 1

    1 -0.0885793
    2 0.0183726 -0.0921237
    3 0.0382677 0.129151
    4 0.0514569
    5 0.0507717 -0.182924

    POPS 10A + LJ-PME cut-off version 2

    1 -0.103516
    2 0.0216568 -0.101543
    3 -0.0308795 0.1414
    4 -0.0200122
    5 0.0505511 -0.151222

  6. And Slipids DOPS:

    DOPS 10A cut-off version 1

    1 -0.0511821
    2 0.0279043 -0.0516061
    3 0.0106204 0.0566416
    4 0.0128629
    5 -0.00319719 -0.169801

    DOPS 10A cut-off version 2

    1 -0.0511015
    2 0.0274585 -0.0513098
    3 0.0176225 0.0769876
    4 0.0258722
    5 0.0139842 -0.170552

    DOPS 15A cut-off version 1

    1 -0.0697978
    2 0.0194815 -0.0641133
    3 0.00299033 0.0879588
    4 0.00744446
    5 0.0109113 -0.163907

    DOPS 15A cut-off version 2

    1 -0.0788824
    2 0.0309259 -0.0778611
    3 -0.0142537 0.117706
    4 -0.00416047
    5 0.0251378 -0.16109

    DOPS 10A + LJ-PME cut-off version 1

    1 -0.0690908
    2 0.0163676 -0.0705811
    3 0.0031525 0.0836897
    4 0.00474312
    5 -0.00189901 -0.157864

    DOPS 10A + LJ-PME cut-off version 2

    1 -0.0653347
    2 0.0291256 -0.0643116
    3 0.00772428 0.0564994
    4 0.00720749
    5 -0.0040382 -0.169001

    1. Do the temperatures for POPS (298K) and DOPS (303K) apply also to Slipids and CHARMM36-UA simulations?

    2. Yes, all the POPS simulations are 298K and all the DOPS 303K

  7. CHARMM36-UA next. Exactly the same set of simulations as for the full all-atom CHARMM36 POPS/DOPS membranes. POPS first:

    POPS 0.8 nm switching version 1

    1 -0.0664802
    2 0.00675274 -0.0703262
    3 -0.225362 -0.231666
    4 -0.225347
    5 -0.206508 -0.0472012

    POPS 0.8 nm switching version 2

    1 -0.0697325
    2 -0.00401857 -0.0729739
    3 -0.234691 -0.231845
    4 -0.226818
    5 -0.205567 -0.0626558

    POPS 1.0 nm switching version 1

    1 -0.0928681
    2 0.0244739 -0.0832559
    3 -0.23423 -0.24347
    4 -0.248209
    5 -0.218282 -0.0506971

    POPS 1.0 nm switching version 2

    1 -0.0939007
    2 0.0251487 -0.0799868
    3 -0.21106 -0.241463
    4 -0.218964
    5 -0.20673 -0.0484276

  8. And CHARMM36-UA DOPS:

    DOPS 0.8 nm switching version 1

    1 -0.0724243
    2 0.0129149 -0.0658414
    3 -0.225781 -0.242195
    4 -0.215034
    5 -0.196189 -0.0343548

    DOPS 0.8 nm switching version 2

    1 -0.0525123
    2 0.00292948 -0.0518382
    3 -0.189082 -0.232873
    4 -0.203974
    5 -0.183291 -0.0412414

    DOPS 1.0 nm switching version 1

    1 -0.0596363
    2 0.0167668 -0.0638123
    3 -0.217616 -0.214942
    4 -0.204914
    5 -0.18654 -0.0560908

    DOPS 1.0 nm switching version 2

    1 -0.0714443
    2 0.0252595 -0.0648029
    3 -0.209222 -0.22916
    4 -0.193383
    5 -0.178898 -0.0592771

  9. Here are some trajectories with the SLIPIDS FF:

    DPPE @336K:

    DOPS @303K:

    POPG @298K:

    DPPG @298K:

    DPPG @314K:

    1. Thanks! I have now added the PS and PG data in Figs. 2 and 3:

      PE figure is yet to be done.

      It seems that the results by Fernando Favela and Thomas Piggot for DOPS are in good agreement. On the other hand, it seems that PS results are in good agreement with experiments for headgroup (as mentioned in the comment below), but PG results are not that good.

    2. I analyzed now also the DPPE data and made a figure comparing the results to experiments:

      Alpha order parameter is too low, but beta shows apparent agreement with experiments. However, only absolute values are compared, because signs are not known experimentally. The sign of beta order parameter is positive in these simulations, in contrast to PC where negative sign was measured. Thus, the the beta order parameter agrees with experiment with the assumption that its sign is opposite than for PC.

      I start to feel like that the signs measurements are almost necessary for more solid conclusions.

  10. I added the above results delivered by Thomas Piggot in Fig. 2. I also splitted the figure to show POPS and DOPS separately. I used only results with cut-off settings closest to original parameters, i.e. 1.0 nm switching for CHARMM36, 15A cut-off for SLIPIDS and 1.0 nm switching for CHARMMua. The new figure is here:

    The results for alpha and beta carbons seems to be pretty close to experiments from all simulation models, assuming that the signs of both alpha and beta carbon order parameters are negative. Slipids has problems in glycerol backbone region, as already observed for PC lipids in NMRlipids I.

    I think that it is remarkable that simulations succeed to reproduce the differences in order parameters between PC and PS headgroups (assuming that the signs are correct in simulations). This means that we could use simulations to interpret the structural differences between these two lipids.

    However, I am slightly confused about the difference between my simulations for POPS:POPC mixture and pure POPS simulations by Thomas. Especially, because pure POPS simulations are closer to experiments from the mixture. I have opened a issue about this:

    I will analyze the data by Fernando Favela soon.

    1. Above comment states "I think that it is remarkable that simulations succeed to reproduce the differences in order parameters between PC and PS headgroups (assuming that the signs are correct in simulations). This means that we could use simulations to interpret the structural differences between these two lipids."

      We have now the signs from PS from experiments, but these do not agree with the Slipids and CHARMM simulations for alpha carbon (simulations give negative, experiments give positive for larger value from alpha carbon). Unfortunately this means that these simulations do no seem accurate enough to interpret the experiments.

  11. I've finally got around to analysing all the united-atom POPS and DOPS simulations. Exactly the same as the ones I've already reported above, so the analysis was performed over the final 100 ns of 500 ns simulations. These simulations used the same starting structures as the all-atom ones, with hydrogens removed as needed and again only have sodium counter ions.

    The first results are for the 'Berger' ff ( I should note here that in these simulations the membrane is very ordered and there are unusual ring-like structures in the head group formed due to overly strong H-bonds. There are two simulations for each of POPS and DOPS performed using 1.0 nm cut-offs, no dispersion correction:

    POPS v1

    1 0.175027
    2 0.168238 0.211975
    3 -0.240138 -0.196859
    4 -0.176208
    5 0.237338 0.146036

    POPS v2

    1 0.168724
    2 0.175157 0.21048
    3 -0.252991 -0.216563
    4 -0.169056
    5 0.226422 0.162045

    DOPS v1

    1 0.223021
    2 0.192289 0.203946
    3 -0.209767 -0.23064
    4 -0.164477
    5 0.225605 0.128681

    DOPS v2

    1 0.140133
    2 0.145168 0.284255
    3 -0.259474 -0.247185
    4 -0.236955
    5 0.328122 0.0649587

    1. I forgot to say PME beyond the coulombic cut-off too, just to avoid any doubts

  12. Next are for simulations with GROMOS-CKP based lipids. All these simulations used 1.4 nm verlet cut-off with a dispersion correction. Simulations were performed with both PME and a RF treatment of the long-range electrostatic interactions. There are a couple of different sets of parameters here. The first uses the same charges as in the Berger simulations, with the NH3 group having the same charges as in the N(CH3)3 group of the PC lipids:

    POPS v1 PME

    1 -0.00265695
    2 -0.0332572 0.126433
    3 -0.211731 -0.331908
    4 -0.190728
    5 0.332247 0.00173317

    POPS v2 PME

    1 -0.0230325
    2 -0.040745 0.12782
    3 -0.215958 -0.340364
    4 -0.200157
    5 0.338137 -0.00398344

    DOPS v1 PME

    1 -0.0108393
    2 -0.0495447 0.127673
    3 -0.186262 -0.344498
    4 -0.162946
    5 0.319748 -0.0178918

    DOPS v2 PME

    1 0.00537431
    2 -0.0416432 0.139164
    3 -0.204886 -0.348506
    4 -0.189437
    5 0.33678 -0.00465369

    POPS v1 RF

    1 0.00420275
    2 -0.0220416 0.129564
    3 -0.205945 -0.358329
    4 -0.190872
    5 0.364139 0.0101082

    POPS v2 RF

    1 -0.0171134
    2 -0.0528662 0.130592
    3 -0.187359 -0.347164
    4 -0.178384
    5 0.350347 -0.0152338

    DOPS v1 RF

    1 -0.0193879
    2 -0.0367713 0.130311
    3 -0.178777 -0.351327
    4 -0.158763
    5 0.343653 -0.0299693

    DOPS v2 RF

    1 0.00268338
    2 -0.0524761 0.150329
    3 -0.206817 -0.350284
    4 -0.187838
    5 0.341264 -0.0133517

    1. Are these still the same system as previous (same amount of lipids, water, etc)?

    2. Yes, all the POPS and DOPS structures are the same across force field. The initial membrane was made with the CHARMM-GUI and converted into the correct format for the other force fields.

  13. Finally, GROMOS-CKP based but with charges for the NH3 group taken from a lysine side-chain (to be more GROMOS consistent):

    POPS v1 PME

    1 -0.00825882
    2 -0.058522 0.0608663
    3 -0.201265 -0.341553
    4 -0.191803
    5 0.334953 -0.0183286

    POPS v2 PME

    1 -0.0138892
    2 -0.0418107 0.0407683
    3 -0.195293 -0.345687
    4 -0.175604
    5 0.310824 -0.0193193

    DOPS v1 PME

    1 0.0092074
    2 -0.0815055 0.0706571
    3 -0.193756 -0.340826
    4 -0.159516
    5 0.340333 -0.022688

    DOPS v2 PME

    1 0.000403073
    2 -0.0627875 0.0625025
    3 -0.189088 -0.341653
    4 -0.163069
    5 0.342013 -0.0305952

    POPS v1 RF

    1 -0.002304
    2 -0.0284508 0.0564798
    3 -0.209456 -0.340959
    4 -0.195232
    5 0.361488 -0.0250153

    POPS v2 RF

    1 0.0120008
    2 -0.0307901 0.0729781
    3 -0.212291 -0.348387
    4 -0.194857
    5 0.356519 0.002936

    DOPS v1 RF

    1 0.00442848
    2 -0.0264082 0.0622449
    3 -0.203746 -0.349032
    4 -0.183483
    5 0.32361 0.0108414

    DOPS v2 RF

    1 0.00224342
    2 -0.0473613 0.0684367
    3 -0.175363 -0.35207
    4 -0.156214
    5 0.322538 -0.0279131

    Sorry for all the numbers (again!)

  14. Thanks again for Thomas Piggot for the extensive data set. I added the DOPS results in Fig. 2. ( It seems that the "Berger" model is quite a bit off, as already anticipated by Thomas above. CKP models are closer, but maybe not within experimental errors (CKP1 refers to the first and CKP2 to the one with NH3 parameters from lysine side chain).

    It should be noted, however, that this and above dicussion is based on the assumption that signs of order parameters are correct in CHARMM and Slipids. If it turns out experimentally that this is not true, this must be reconsidered.

    POPS results are yet to be added.

  15. After digesting this data for a while, I think that we must get the signs of the headgroup order parameters measured also for PE, PS and PG to make any solid conclusions.

    I will try to get this organized.

    1. It seems currently very promising and I hope that we have the signs within few months.

    2. Tiago Ferreira ( measured the signs for the headgroup order parameters of POPS and delivered the results by email. The details are included now in the manuscript (SI part):

      The conclusion is that the headgroup order parameters for POPS at 298K from 13C NMR are:
      -0.12 for beta carbon
      0.09 and -0.02 for alpha carbon.

      This means that the signs in Slipids and CHARMM simulations for alpha carbon of PS (simulations give negative, experiments give positive for larger value from alpha carbon) do not agree with experiments.

    3. I have now received also order parameter values for g1 and g2 segments of POPS measured by Tiago Ferreira:

      g1 -0.13 and 0
      g2 -0.18
      (the error should be about +- 0.02)

      For g3 the resolution was too low to get a proper estimate.

      I have included these values in plots in manuscripts.

  16. Hi Samuli,

    I have (yet unanalyzed) simulations for PC/PG head group mixtures that I can contribute. Simulations are with 500 lipids in mixing rations PC:PG of 1:1 and 4:1, simulated for 200 ns each using the CHARMM36 FF, with 1M/0.15M calcium chloride or neutral (only potassium counter ions).



    1. Excellent!

      Can you share the trajectories and other files in Zenodo, and/or do you want just to deliver the order parameters? The scripts in here, for example, should be almost directly applicable:

      Looking forward to see the results.

    2. Sure, I can arrange for the raw files to be uploaded and meanwhile I can use the bash wrapper to get the order parameters. I could, however, not find the python script that the wrapper points to ($scriptdir/ If you can direct me toward that script, I can get on it.


    3. It should be in this folder:

    4. Hi Jesper,

      how did you determine the concentrations in these systems? Could you deliver the number of water, ion and lipid molecules in the systems?

    5. Hi Samuli,
      The reported concentrations are inputted to the CHARMM-GUI in constructing the systems and therefore *not* estimates of experimentally measurable bulk salinity. I am away for this week and can return with specification of component numbers upon my return.


  17. Thanks!

    Here are the calculated order parameters (whole trajectory):

    I did not run the system 4:1 system with 0.15MCaCl2.


    1. Thank you!

      I plotted this data with experimental data.

      PC order parameters as a function of CaCl_2 concentration:

      The decrease of alpha-order parameter is in agreement with experiments, while decrease of beta order parameter is overestimated. The result is very similar to the results with CHARMM36 PC in NMRlipids II publication. The conclusion in there was that the Ca2+ binding is too strong in CHARMM36 model despite of the seemingly good agreement of alpha carbon with experiments. The reasoning was that the dependence of alpha-carbon order parameter on bound charge is too weak in CHARMM36, while beta-carbon is more realistic (see Fig. 3 in NMRlipids II publication). However, this was not explicitly tested against experimental data. I have now done some tests like this for Lipid14 model, which I am planning to report soon. After seeing these results, it seems that these tests would be useful for CHARMM36 as well.

      It should be also noted that the beta-order parameters are not actually measured for PG. The experimental line in the figure is calculated from empirical relation \Delta S_beta = 0.43\Delta S_alpha.

      Thus, my preliminary conclusion from the data would be that the Ca2+ binding is similarly (not more or less seriously) overestimated in CHARMM36 simulations with PG than in simulations with only PC.

      PG order parameters as a function of CaCl_2 concentration:

      Absolute value of the order parameter is too large without ions, but rapid decrease due to addition of CaCl_2 is observed in agreement with experiments for systems with 1:1 mixture of POPC
      and POPG. In addition, absolute value in systems with CaCl_2 is in agreement with experiments. However, system with 4:1 mixture of POPC and POPG behaves differently, but experimental data is not available for direct comparison for this mixture.

      In conclusion, the qualitative response of PG on CaCl_2 seems to agree with experiments. However, the difference between 1:1 and 4:1 mixtures is mysterious.

      Ion and lipid density profiles would be also useful from these simulations.

      There would be now some data to check also lipid headgroup changes due to mixing them (e.g. how PC order parameters change with PG). This might reserve some attention. This is also related to this issue:

      Also simulation data of ion binding on PG with other models and of ion binding on PS bilayers with any model would be highly useful at this point.

    2. Hi Samuli,
      (Mass) density profiles calculated for the major system constituents of the POPC_POPG/CaCl_2 simulations can be found here:
      The density profiles suggest what you would expect: that calcium ions bind tighter to the more negatively charged 1:1 PC/PG membrane surface as compared with the 4:1 PC/PG membrane.


  18. Hi Samuli,
    I have uploaded the trajectories of POPC/POPS (4:1) with 1M NaCl, KCl, CsCl in Berger for lipids and and Dang/gmx for ions.


    1. Thanks Lukasz!

      Would you also have the trajectory with only counterions?
      That would be needed for the reference when using electrometer concept.

    2. I will run the counterions-only systems.

    3. I found experimental data for headgroup order parameters of DMPC/DMPS (3:1) mixture as a function of NaCl concentration:

      This can be compared to the reported simulation data as soon as we get the simulation with counterions only.

    4. I ran the system with only Na+ conterions, data is not yet in Zenodo but some files are here:

      I also plotted the headgroup order parameter dependence of NaCl against experimental data for POPC:POPS (5:1) system. Results and further discussion is here:

  19. Hi All,

    We have looked at the headgroup behaviour and tried to visualize the differences for the available PS and PC trajectories:
    - Slipids DOPS trajectory
    - Slipids POPC trajectoru
    - CHARMM36 DOPC trajectory (our own trajectory, which we can share if needed)

    For glycerol - the differences can be clearly visualized:
    Evidently, some isomers are absent in Slipids, which might explain its order parameters in Botan et al. All of the isomers are present in CHARMM36, as discussed in Klauda et al., 2010.

    For phosphate - the differences between the lipids are not so clear:

    It is easier to compare the trajectories via plotting the dihedral angle distributions:
    (the names of the atoms correspond to CHARMM36)

    The conclusions are:
    1) There are significant differences between CHARMM36 and Slipids. Perhaps it is better to use CHARMM simulations for headgroup comparisons.
    2) For PS and PC in Slipids trajectories, there are also large differences in some of the not-glycerol dihedrals (dihed8, 9, 10). Consequently, order parameters for respective hydrogens (at C11 and C12) might also be different.

    Pavel & Ivan

    1. Thank you for your contribution!

      I find your visualization for the glycerol structures very useful. I did already include it in the current version of the manuscript, but we need to find a proper place for this discussion when the manuscript is closer to its final shape. How did you selected the structures shown in the figure?

      The dihedral comparison is also very useful. I am slightly troubled with the observation that the differences between slipids and CHARMM seems generally larger (with some examples as mentioned) than the differences between PC and PS. I have just uploaded the experimental information about signs of POPS headgroup and neither of the force fields do not succeed very well in the comparison. We need some careful thinking now how we should compile these results. Do you have some script to calculate these dihedral angles? It might be useful to run it also for other simulations as well. I am in (slow) process of collecting all the reported data in the table with links to raw data in Zenodo.

    2. Hi Samuli,

      regarding your questions:
      1. The structures were selected as following: first the trajectories of individual lipids were concatenated into a single trajectory. Then all frames were aligned to the first one using the glycerol atoms C1-C2-C3-O21. Then 512 frames were selected evenly throughout the trajectory.
      2. Regarding the scripts, currently it is a VMD tcl script to calculate angles and a Mathematica script to draw the distributions. We can try to rewrite the scripts with python and then share them, or we can process the other trajectories ourselves (in this case, please share the list of the trajectories of interest).

      Finally, we will share our trajectories with Zenodo soon.


    3. After compiling the results it seems to me that it would be useful to calculate the dihedral distributions also from these trajectories:
      1. CHARMM36 DOPS
      2. CKP1 DOPS
      3. CKP2 DOPS.
      The order parameter data for all of these is delivered by Tom Piggot and is included in Figure:

      CHARMM36 DOPS gives quite similar order parameters as Slipids DOPS. CKP models are better for alpha carbon but not for the beta. Comparing dihedrals between these models we might learn something. I would also include the dihedral for the bond between beta carbon and carbonyl of PS headgroup in the analysis to see differences between CHARMM, Slipids and CKP for this.

      As mentioned, the data for all the proposed analyses is delivered by Tom Piggot but the trajectories are not (yet) in Zenodo. Would it be possible to upload them or to organize the analysis in some other way? I also think that we should choose which one of the switch version we use for CHARMM and if we use RF or PME data for CKP. To perform all the analyses for all of these options gets a little bit too complicated.

    4. I'll get the CHARMM36 and GROMOS-CKP simulations reduced in size (the frames are every 5 ps in the original trajectories meaning these files are about 7-12 GB) and uploaded to Zenodo asap.

      With CKP I will just upload the PME simulations. The RF ones are actually being re-run at the moment with slightly different settings and even though I don't think this will change the results much, the PME ones should be used for now.

      As for CHARMM36, I guess we've had this discussion elsewhere ( I'm happy with picking either as there are arguments that can be made for both and the results are very similar. Perhaps the switching point at 10A?

    5. In my understanding most people nowadays would use CHARMM36 with switching point 10Å due to CHARMM GUI and Gromacs instructions, so maybe this would be better choice now (although I understand your arguments for 8Å as well).

    6. I've finally got some of the pure POPS and DOPS simulations uploaded onto Zenodo (think of it as an early Christmas present!). For each of the force fields I simulated, I have chosen representative versions to upload (e.g. the 10 Å switching point for CHARMM36):

      The trajectories have been reduced to have a frame every 50 ps and only include the final 100 ns of 500 ns simulations. It's worth stressing that the order parameter data provided above was calculated over the same simulation time but was performed on the original trajectories, so with a frame every 5 ps.

      If anyone needs either the longer trajectories or ones with more frames just let me know.



    7. I've finally got around to uploading the pure POPS and DOPS simulations to Zenodo:

      The trajectories I've uploaded are for the final 100 ns of the 500 ns simulations, with the original trajectories also reduced in size to now have one frame every 50 ps (rather than every 5 ps, as analysed in the data provided above). I've also only uploaded one representative simulation setting for each POPS and DOPS force field I simulated (e.g. the 10 Å switching point for the CHARMM simulations) to make things easier. As the data above shows, there is relatively little impact of these settings upon the calculated order parameters.

      If anyone wants the longer trajectories, more frequent frames, or some of the simulations with other settings please let me know.



    8. Thanks a lot for sharing the data!

      I have now updated the Table I in the current manuscript to contain the links to the data (

      It seems to me now that it would be best to run the above dihedral analysis to the delivered POPS systems listed in the table. For POPS we have the widest collection of different models and experimental data points. The links to the data of all delivered POPS systems are now in the table I in the manuscript, thus such analysis should be doable now.

      The correct references for CHARMM36, CHARMM36ua, GROMOS-CKP and Berger POPS are yet to added in the table. If someone has these at hand, it would be useful to fill them or deliver here.

    9. Hi All!

      Finally, we have performed the analysis for all the POPS systems available. For the glycerol conformations, the differences between different force fields could be clearly visualized again. Interestingly, for Berger force field the improper C1-C2-O21-C3 is very flexible. Also, there are quite noticeable differences between CHARMM36 and CHARMM36UA, though the isomer parametrization is the same for both force fields.

      Dihedral angle distributions show significant differences between force fields. Slipids and CHARMM36 force fields result in quite similar distributions, though still there are some differences. The other force fields differ from each other and from CHARMM36 and Slipids on a greater scale.

      All the figures could be find via link:

      The dihedral distributions are somewhat imperfect. Thomas, if it is possible, could you please share with me your trajectories with more frequent output, to make all the distributions smoother? Thank you in advance.

      Pavel & Ivan

    10. Thanks a lot!

      I think that the dihedral distributions are smooth enough for our purposes, at least for now.

      After a couple of quick looks, I cannot see any clear correlation between dihedral distributions and order parameters. I would be happy to hear any ideas how we could say which one of the distributions are more realistic.

      How about the last dihedral of the headgroup, i.e., C12-C13 in CHARMM notation? It would be highly interesting, but I do not see a figure for it.

      There is recent data also for Amber Lipid 17 model ran with the Amber program:
      Would it be possible to analyze this Amber produced data with your script?

    11. Hi All!

      I have added distributions for C12-C13 dihedrals. They could be found via link:

      As for Amber Lipid 17 trajectories, we can analyze them, but I couldn't find any structure file, psf of pdb, which I need for the analysis. Could you please add them? Then we will calculate the distributions.

      Regarding the correlation between dihedral distributions and order parameters, Jeffery Klauda with colleagues discussed it in the original CHARMM36 paper (second paragraph in Section 3.2, They calculated S_cd for different dihedral states and then tried to predict optimal populations of different conformers. Perhaps the same could be done here.

      Pavel & Ivan

    12. Hi Pavel,

      I already have calculated the dihedral distributions for Amber Lipid 17 forcefield. You can access the files via

      If the link doesn't work, please let me know.

      Since I don't have the data for the other ff's, I only calculated the raw angles and did not plot them.

      If you want to recheck the calculations using your script, I can upload a pdb file to zenodo for Amber trajectories.


    13. Hi all,

      I've uploaded the pdb file for DOPS with Amber Lipid17 simulations to:


    14. And the PDB of POPS with ff99 ions can be found in here.

    15. After focusing only to the headgroup (not glycerol backbone) part of the dihedral distributions (Fig. 16 in the current manuscript), I came up with this possible idea.

      The main differences between the models in the headgroup region are observed for dihedrals C12-C11-O12-P and C11-C12-C13-O1A. The distributions from CHARMM36, CHARMM36UA and Slipids for the latter dihedral C11-C12-C13-O1A are similar, all having only a single peak around 120 degrees.
      This dihedral is close to the beta-carbon and its order parameter is best reproduced by the CHARMM36, CHARMM36UA and Slipids models.

      On the other hand, Gromos-CKP models give better order parameters for the alpha-carbon. Distribution of C12-C11-O12-P dihedral close to alpha carbon gives a single peak around 180 degrees in these models.

      In conclusion, the suggestion would be that CHARMM models and Slipids give more realistic distribution for C11-C12-C13-O1A dihedral (single peak around 120 degrees), while CKP models would be more realistic for C12-C11-O12-P dihedral distribution (single peak around 180 degrees).

      To see if this idea would be also in line with Lipid17 results, it would be useful to get the results included in the plot.

    16. Now the results from MacRog are added to the Fig. 5 and they seem to be pretty good for the alpha carbon. Therefore the better option for C12-C11-O12-P dihedral could be the distribution from MacRog peaking at around 100 and 240 degrees.

    17. Hi all!

      I have added distributions for POPS simulated with Lipid17 FF with different ions. The results could be found via link:

      In the Lipid17 trajectory with jcions ( the head groups look almost frozen - should it be like that?

      Pavel & Ivan

  20. Hi Samuli,
    Amélie Bacle and I have simulated some PC/PE mixtures at different ratios (POPC/DOPE and DOPC/DOPE at 100:0, 75:25, 50:50) using Berger FF (as implemented by Luca Monticelli, i.e. compatible with OPLS-AA proteins, at 300 K, 300 ns each. If it can be of interest to NMRlipidsIV, we would be happy to share the trajectories and calculate the order parameter of the polar head protons.

    Patrick and Amélie

    1. Thanks. I think that at least the PE order parameters calculated from the highest mole fraction would be useful to get an idea how this model performs for PE. Possibly also the effect of PE on PC order parameters might be interesting for the project. There is, however, no (as far as I can remember) experimental data for PC/PE mixtures so this cannot be directly compared. It might be useful for general understanding, thought.

    2. FYI PE membranes simulated with the typical Berger PE ff (so with just the methyl groups in the PC choline switched for standard hydrogen atom parameters) in general won't behave well. The hydrogens in the ethanolamine group form exceptionally long-lived hydrogen bonds with the phosphate oxygens generating ring-like structures in the head group. This is a general issue that also arises with PG and PS 'Berger' lipid parameters. The APL of pure PE memembranes simulated using these parameters is also far too low. There are references for this (let me know if you wan't them?), plus I've simulated these in the past and can confirm these issues.

      You can get around some of the issue by adding a repulsive potential onto the ethanolamine hydrogens. However, AFAIK, the only published version of this fix is for DOPE ( When I tried it with POPE, I had to slightly increase the repulsive potential to get good overall membrane properties (IIRC, this was done ages ago but I should still have the files somewhere, if anyone is interested?).

    3. Hi Samuli,

      Amélie has calculated the order parameter but I want to check something about the sign. So we'll come back soon with the values.

      @Thomas: thanks for the pointers. I think it's still interesting to see how "untouched" Berger lipids behave for a systematic comparison (like in NMRlipids I). And then maybe try the fix you're talking about and see what happens on the order parameter.
      About using APL for validating pure PE simulations, this can be problematic when the lipids are not forming a lamellar phase, such as pure DOPE (which is in the H_II phase at room temperature). And this is the only PE we used in our simulations. I haven't dig too deeply into literature, it looks like there are a few NMR studies on DOPC/DOPE mixtures (using 31P & 2H), but I don't know yet if they are usable for validation. Still, I'm interested in the refs you're talking about, I guess the simulations were done on lipids forming a lamellar phase at room T.



    4. Hi Patrick,

      The sign information has been now updated in the manuscript, see also the comment below.

      I would be also interested about the references for 31P and 2H experiments of DOPC/DOPE mixtures.

    5. Sorry for my long silence, but now Amélie and I could check the sign.

      I recall, the simulation conditions were :
      - 300 K with v-rescale (tau=0.1 ps)
      - 1 bar with PR semiisotropic (tau=4 ps, compressibility=4.5e-5 bar^-1)
      - PME order 4 and space 0.12
      - rcoulomb and rvdw 1.0
      - 128 lipids per leaflet
      - no ion

      Here are the results :

      I) DOPC/DOPE

      - 100:0 (400 ns, analysis on the last 300ns, dt=100ps)
      Beta_1 0.040805
      Beta_2 0.0538425
      Alpha_1 0.0735451
      Alpha_2 0.14256
      G3_1 -0.314431
      G3_2 0.0721855
      g2 0.390806
      G1_1 0.0314617
      G1_2 -5.10e-05

      - 50/50 (300 ns, analysis on the last 200ns, dt=100ps)
      Beta_1 0.00273305
      Beta_2 0.0114694
      Alpha_1 0.00584493
      Alpha_2 0.0739441
      G3_1 -0.256777
      G3_2 -0.0988909
      g2 -0.106475
      G1_1 0.214974
      G1_2 0.0179717

      Beta_1 -0.0113632
      Beta_2 -0.0268633
      Alpha_1 0.107567
      Alpha_2 0.134572
      G3_1 -0.279772
      G3_2 -0.152953
      g2 -0.15289
      G1_1 0.27854
      G1_2 0.0389075



      - 100:0 (400 ns, analysis on the last 300ns, dt=100ps)
      Beta_1 0.0394835
      Beta_2 0.0529308
      Alpha_1 0.0750883
      Alpha_2 0.145688
      G3_1 -0.25078
      G3_2 -0.163228
      g2 -0.0729539
      G1_1 0.1984
      G1_2 0.0834135

      - 50:50 (300 ns, analysis on the last 200ns, dt=100ps)
      Beta_1 0.00301153
      Beta_2 0.00742294
      Alpha_1 0.0075323
      Alpha_2 0.092868
      G3_1 -0.257383
      G3_2 -0.145337
      g2 -0.151157
      G1_1 0.205346
      G1_2 0.0734792

      Beta_1 0.0044169
      Beta_2 -0.034363
      Alpha_1 0.0990313
      Alpha_2 0.169923
      G3_1 -0.294483
      G3_2 -0.137261
      g2 -0.137572
      G1_1 0.242008
      G1_2 0.0366029

      A few comments:
      - The PC choline H (alpha & beta) are closer to 0 when we add PE, expected since there's more room for them to move
      - The results on PC glycerol H (g1, g2, g3) are less obvious. I see 2 things
      --> The effect of adding PE is marginal on POPC but significant on DOPC
      --> I'm quite surprised about the massive effect on g2 on going pure DOPC to DOPC/DOPE 50:50.

      Not easy to compare these results to the ones reported in the manuscript for DPPE as when we have PE, it's always a mixture with PC. Also I recall, DOPE is non-lamellar at room T.

      With Amélie, we plan to give a try to the fix Thomas mentionned about (putting a slight repulsion on choline H). We'll report back in a few weeks once the simulations have completed.

      About the refs on PC/PE mixtures using 2H and P31, this was a first quick scan :
      - Tilcock et al, Biochemistry 1982, 21, 4596-4601 (DOI: 10.1021/bi00262a013): in this one they labeled C11 with 2H, they study various DOPC/DOPE mixtures (also w/wo cholesterol)
      - Farren et al., Chemistry and Physics of Lipids, 34 (1984) 279-286 ( in this one they labeled C11 with 2H, and they could obtain results on DOPC, DOPE, DOPS, DOPG and DOPA !


    6. Thanks for the data!

      Based on this and on the other results and discussion ( and, I think that we have to pay some attention also to the response of the headgroup order parameters to the lipid mixtures.

      I plotted the PC headgroup order parameter results together with the experimental data (

      Experimentally the PC headgroup order parameters in the mixed lipid bilayers follow the electrometer concept: neutral lipids have only a minor effect, while negatively charged lipids increase the order parameters. The Berger-OPLS simulations seems to overestimate the order parameter change due to the presense of PE, while CHARMM36 simulations seems to underestimate the change due to the negatively charged lipids.

      On the other hand, the effect of PC to the PS headgroup is not very well captured in CHARMM36:

      I think that we need to pay some attention to this, because this may affect to the usage of the electrometer concept in ion binding studies in simulations. I will try to run the cationic surfactant simulations with CHARMM36, as was done for Lipid14 already in:

    7. Hi Thomas,

      you wrote above:
      "PE membranes simulated with the typical Berger PE ff (so with just the methyl groups in the PC choline switched for standard hydrogen atom parameters) in general won't behave well. The hydrogens in the ethanolamine group form exceptionally long-lived hydrogen bonds with the phosphate oxygens generating ring-like structures in the head group. This is a general issue that also arises with PG and PS 'Berger' lipid parameters."

      I can see the mentioned conformations in POPS in the Berger simulations used to look excess sodium binding to POPC:POPS (4:1) mixture:

      Would you be aware of a reference where the presence of these conformations would be discussed? It would be easier to cite such work, rather than describe or plot the conformations in the manuscript.

    8. Hi Samuli,
      after some months of silence, I could come back to the PC/PE simulations. Amélie and I indeed observed what Thomas described, so the numbers I posted above come from these simulations, and may not be trusted. I'd be highly interested if there's a ref discussing this as well. If none exists I can generate the plots and post them (dihedral distrib + H-bonds).
      In between, we could do the simulations with the hack mentionned by Thomas (adding a repulsion on ethanolamine hydrogens), so we'll come back soon to post the results overhere.
      Last, a new student (Christos) working with me has launched some POPC/POPE simulations using the CHARMM 36 force field. We can post also these results if you think this is of interest (we can also re-launch some new simulations with DOPC/DOPE or POPC/DOPE to be consistent with the Berger compositions mentionned above and post only those new ones), let me know.

    9. The first mention of this sort of problem was for Berger PE and in:

      where the authors say:

      "It should be noted that this area of 0.58 nm^2 for DPPE was achieved only after introducing a small repulsive interaction between the ammonium H atoms on the PE headgroup and any other atom, except the water H. Without this interaction, the area per headgroup for a DPPE bilayer at 343K was found to be 0.52 nm^2."

      So not explicitly discussing ring like structures in the head group but showing that without a repulsive potential on the hydrogens you get a lower APL.

      Much greater detail is discussed for Berger based PG lipids in:

      including a discussion of the extremely long-lived interactions in the head group with Berger based (i.e. Zhao and Elmore) PG force fields.

      As for PS, this large degree of intra-molecular hydrogen bonding is discussed in the original Berger PS paper:

      and you can even see pictures of this happening in Figure 6 of that paper. That said, they do not discuss that this might be an issue with the force field.

    10. Hi Patrick,

      CHARMM36 simulations of PE in mixtures and pure bilayers are highly interesting. I think that POPC/POPE is the most interesting mixture since there is experimental data from that one ( Right now I do not think that it is necessary to run the DOPC/DOPE or POPC/DOPE mixtures with CHARMM36 for this project.

      Dihedral angles are already analysed for different PS models, see Fig. 17-18 in the current PS manuscript:

      I would be happy hear any comments of these from the point of view of these weird structures in Berger.

    11. Hi Thomas,

      thanks for the references.

      The ring like structures are clearly visible in Fig. 6 in the last paper you mention. I think that we can cite it in the PS manuscript and write something like "The poor performance of headgroup order parameters in Berger model can be probably explained by ring like structures seen in Fig. 6 in Ref."

    12. To that, I would also add something like:

      "These ring like structures are a widespread feature of typical Berger based lipid force fields containing explicit hydrogen atoms in the head group [REFS and and]"

      The second of these references is another I didn't mention above which has further details for the PE ring like structures (see Table 1). The final reference is for Berger based cardiolipin where they say in their discussion:

      "Setting the LJ parameters to zero, which is the standard treatment for hydrogens, creates excessively stable, unphysical, bonds between the hydroxyl group and OA2/OB2 phosphate ester oxygens"

    13. I have lots of pure PE bilayer simulations I will look out and share over the next few days. They are not in such a nice consistent way as the ones I have already shared for PS and may be for bilayers of different sizes, with different versions of GROMACS, etc.

      I'm pretty sure I should have some GROMOS-CKP, GROMOS43A1-S3, CHARMM36, CHARMM36-UA, Berger and OPLS-UA simulations (with and without LJ on the Hs for the last two). This is as long as I can find them amongst my older files. Most will be for POPE but there will also be some DPPE and DOPE (at least for CKP but maybe also some others).

      I will get back to you when I have found everything.

    14. Hi Samuli,
      Great, I'll post the CHARMM36 results on POPC/POPE (50/50) very soon. And you're right, POPC/POPE (50/50) is the most interesting regarding the comparison to experimental data.
      For the Berger simulations with the Marrink hack to avoid persistent H-bonds (BTW, how should we call these lipids? For this post I'll call them "PE Berger-hacked"), I will post the results on what we have so far (POPC[plain Berger]/DOPE[Berger-hacked] 50:50 and DOPC[plain Berger]/DOPE[Berger-hacked] 50:50). But I propose that Amélie and I launch another simulation of POPC[plain Berger]/POPE[Berger-hacked] (50/50) so that we'll be able to compare to CHARMM 36 and hopefully to experimental data (didn't go deeply into the experimental ref yet). It'll be also very interesting to compare to pure PE simulations from Thomas. So I'm back soon (before the end of the week) with all the results we generated so far!

    15. Finally, we have gathered and checked all simulations. First, here are the results with the plain Berger (P|D)OPC and Berger-hacked DOPE at a molar ratio 50:50. These simulations were generated and analyzed by Amélie Bacle and myself.

      I) DOPC/DOPE 50:50 (300 ns, analysis on the last 200ns, dt=100ps)

      OP of DOPC
      beta_1 0.0478875
      beta_2 0.0666199
      alpha_1 0.0838629
      alpha_2 0.1564830
      g3_1 -0.2664870
      g3_2 -0.1553530
      g2 -0.1421150
      g1_1 0.1967150
      g1_2 0.0751805

      OP of DOPE
      beta_1 0.0651575
      beta_2 0.0285274
      alpha_1 0.0907662
      alpha_2 0.1345950
      g3_1 -0.2711270
      g3_2 -0.1548510
      g2 -0.1429390
      g1_1 0.1830130
      g1_2 0.0762804

      For comparison, here are the values for pure DOPC (plain Berger) (there was a mistake in the values for pure DOPC (plain Berger) described in my post of September 2, 2017 at 6:29 PM, so please ignore them and condider the following ones, they are consistent with those published in NMR_lipids_I):

      OP of DOPC (400 ns, analysis on the last 300ns, dt=100ps)
      beta_1 0.0447022
      beta_2 0.069103
      alpha_1 0.0928824
      alpha_2 0.151711
      g3_1 -0.271082
      g3_2 -0.156422
      g2 -0.14668
      g1_1 0.194481
      g1_2 0.0758857

      II) POPC/DOPE 50:50

      OP of POPC (300 ns, analysis on the last 200ns, dt=100ps)
      beta_1 0.0447890
      beta_2 0.0661732
      alpha_1 0.0984351
      alpha_2 0.1508860
      g3_1 -0.2674230
      g3_2 -0.1708150
      g2 -0.1584100
      g1_1 0.1885990
      g1_2 0.0843454

      OP of DOPE
      beta_1 0.0612797
      beta_2 0.0284065
      alpha_1 0.0910532
      alpha_2 0.1339060
      g3_1 -0.2683180
      g3_2 -0.1582320
      g2 -0.1459440
      g1_1 0.1799030
      g1_2 0.0746909

      For comparison, here are the values for pure POPC (plain Berger) (there was a mistake in the values for pure POPC (plain Berger) described in my post of September 2, 2017 at 6:29 PM, so please ignore them and condider the following ones, they are consistent with those published in NMR_lipids_I):

      OP of POPC (400 ns, analysis on the last 300ns, dt=100ps)
      beta_1 0.041295
      beta_2 0.0655363
      alpha_1 0.095803
      alpha_2 0.153548
      g3_1 -0.264093
      g3_2 -0.169425
      g2 -0.16333
      g1_1 0.186404
      g1_2 0.089304

      The conditions of the simulations were the following:
      - 300 K with v-rescale (tau=0.1 ps)
      - 1 bar with PR semiisotropic (tau=4 ps, compressibility=4.5e-5 bar^-1)
      - PME order 4 and space 0.12
      - rcoulomb and rvdw 1.0
      - 128 lipids per leaflet
      - no ion

    16. Here are the results using CHARMM 36 POPC/POPE 50:50. For the reference, we also put some simulations of pure POPC (of course with CHARMM 36 also, and we find again very close values to the ones published in NMR_lipids_I). The simulations were generated and analyzed by Chris Papadopoulos and myself.

      I) Pure POPC (300 ns, analysis on the last 250ns, dt=100ps)
      beta_1 (H12A) -0.077832
      beta_2 (H12B) -0.076815
      alpha_1 (H11A) 0.035139
      alpha_2 (H11B) 0.033111
      g3_1 (HA) -0.221718
      g3_2 (HB) -0.245360
      g2 (HS) -0.192003
      g1_1 (HX) -0.173083
      g1_2 (HY) -0.039947

      II) POPC/POPE 50:50 (300 ns, analysis on the last 250ns, dt=100ps)

      OP of POPC
      beta_1 (H12A) -0.097157
      beta_2 (H12B) -0.099216
      alpha_1 (H11A) 0.025457
      alpha_2 (H11B) 0.036664
      g3_1 (HA) -0.213187
      g3_2 (HB) -0.232995
      g2 (HS) -0.205492
      g1_1 (HX) -0.183274
      g1_2 (HY) -0.038663

      OP of POPE
      beta_1 (H12A) -0.008352
      beta_2 (H12B) 0.003849
      alpha_1 (H11A) 0.094556
      alpha_2 (H11B) 0.061279
      g3_1 (HA) -0.214143
      g3_2 (HB) -0.222699
      g2 (HS) -0.197120
      g1_1 (HX) -0.180223
      g1_2 (HY) -0.035411

      The conditions of the simulations were the following:
      - 300 K with v-rescale (tau=0.1 ps)
      - 1 bar with PR semiisotropic (tau=4 ps, compressibility=4.5e-5 bar^-1)
      - PME order 4 and space 0.12
      - rcoulomb and rvdw 1.0
      - 128 lipids per leaflet
      - no ion

    17. Here are the results with the plain Berger POPC and "Berger-hacked" POPE at a molar ratio 50:50. These simulations were generated and analyzed by Amélie Bacle and myself.

      POPC/POPE 50:50 (400 ns, analysis on the last 300ns, dt=100ps)

      beta_1 0.0397359
      beta_2 0.0671129
      alpha_1 0.0870615
      alpha_2 0.157169
      g3_1 -0.26354
      g3_2 -0.161917
      g2 -0.153312
      g1_1 0.186221
      g1_2 0.0796706

      beta_1 0.0619824
      beta_2 0.0266716
      alpha_1 0.101611
      alpha_2 0.141147
      g3_1 -0.271841
      g3_2 -0.185615
      g2 -0.169221
      g1_1 0.178088
      g1_2 0.0976236

    18. A few comments about my posts of today. I will focus on POPC/POPE 50:50 as suggested by Samuli. The values between the two force fields are very different, so we made an analysis of the area per lipid to try to understand what is going on:
      Area results in nm^2, the error is <= 0.003 nm^2
      - pure POPC
      CHARMM36: 0.624
      Berger : 0.649
      - POPC/POPE 50:50
      CHARMM36 : POPC 0.609, POPE 0.557
      Berger-hacked: POPC 0.637, POPE 0.632
      One can see that CHARMM 36 predicts a drop in the area on going from pure POPC to POPC/POPE 50:50. This means that POPC pack tightly to POPE.
      In contrast, the values for Berger are not that changed.
      The POPE value predicted by CHARMM 36 (in the mixture POPC/POPE 50:50) is much smaller than that predicted by Berger.

      I think it would be interesting to have values for pure POPE for both force fields. From what I read, the phase transition between the gel to the liquid-crystalline phase is 25°C (298 K), and the transition between lc and inverted hexagonal (H_II) is 71°C (344 K). So simulating pure POPE at 300 K should be OK.

      @Thomas: You said you had performed such simulations, do you plan to share them here? Otherwise, let me know and we can generate them.

      About pure POPE, I found these references of 2H NMR which might be useful:
      (they also study the effect of cholesterol/ergosterol, but that is another question!)

      @Samuli: Last questions/remarks:
      - Shall we upload the trajectories to Zenodo? If yes, do we upload all of them or only POPC/POPE 50:50?
      - So far I didn't have time to compare all these values to experimental data, also I'm not clear whether we'll be able to extract easily S_CD values from the papers we spotted. Anyway, if there are plans for that, please let me know (I'll be very busy until May 21st).

    19. Thanks for the data!

      One question before more discussion. You wrote that the results reported September 2, 2017 at 6:29 PM for pure POPC (plain berger) are not correct. How about the results for POPC:DOPE and DOPC:DOPE mixtures reported on the same day? Are those correct?

    20. Hi Samuli,
      no there were mistakes as well on the results for POPC:DOPE and DOPC:DOPE mixtures posted on September 2 2017, please consider only the data posted on April 27 2018.

    21. Hi,

      Firstly, really sorry for the delay with this.

      I have found the majority of the simulation files for PE I think I should have (some were archived, some on old computers, etc.) and I'm now getting all the files together to be able to share them. The files I have found are:

      GROMOS-CKP (and variants)

      Quite a few pure DPPE, POPE and DOPE 128 lipid bilayer simulations with two different head group charges (like the PS simulations). Also some other variants too with different charges. Also some POPE simulations of a 336 lipid hexagonal membrane.


      I can only find simulations with the vdW parameters included on the hydrogen atoms (for POPE and also DOPE I think, but I need to check these latter files to ensure they are what I think they are. Both of these were performed using a couple of different LJ parameters for the hydrogen atoms and were of standard 128 lipid bilayers). I definitely also did some "normal" Berger PE simulations but I cannot find them as yet. I can run some of these if needs be.

      OPLS-UA (constructed from the Ulmschneider and Ulmschneider PC model)

      POPE 128 bilayer simulations with and without vdW parameters on the hydrogen atoms (as in the Berger simulations).

      GROMOS 43A1-S3

      Pure POPE 128 lipid bilayer simulations


      336 POPE lipid hexagonal membrane simulations.


      336 POPE lipid hexagonal membrane simulations.

      I should also have some normal CHARMM36 POPE simulations somewhere too but I cannot locate them at the moment either.

      For all of the above simulations, there are two simulations using different starting velocities. Simulation lengths are typically 200 or 500 ns. DPPE simulations will have been performed at 342 K, DOPE at 271 K and POPE at either 303 or 313 K.

      I'll get these all checked, analysed and uploaded to zenodo as soon as I can. If there are any that anyone would especially like first please let me know.

    22. I have not managed to perform the order parameter analysis on these but to avoid further delay I am uploading the simulations to zenodo today. So far there is:

      GROMOS 43A1-S3 POPE


      Slipids POPE

      OPLS-UA POPE (with/without vdW for the H atoms)

      Berger and GROMOS-CKP will follow later today.



    23. Berger-based POPE (2 different H vdW params)

      Berger-based DOPE (2 different H vdW params)

    24. For GROMOS-CKP, I have just uploaded one set of the parameters/simulations to make things simpler:




      If anyone wants more GROMOS-CKP PE simulations with reaction-field or some different charges please let me know.



    25. I have uploaded the CHARMM 36 simulations on Zenodo :
      - pure POPC:
      - POPC 50% / POPE 50%:
      Amélie will shortly upload the Berger simulations mentionned in my earlier posts.

    26. Hi,

      I've uploaded the Berger simulations on Zenodo :
      - Pure POPC :
      - Pure DOPC :
      - DOPC 50%/ DOPE 50% :
      - POPC 50%/ DOPE 50% :
      - POPC 50%/ POPE 50% :


    27. I finally started to compile PE data as well. Simultaneously, I also wanted to make the analysis process more automatic to handle the massive amount of data with less manual work.

      However, I realized now that most of this data is from united atom simulations. Previously in this project we have used the Gromacs protonate command to add the protons and then regular scripts to calculate the order parameters. This is, however, little bit annoying because protonate command does not work with the newer versions of Gromacs. Thomas has a recent publication and code to calculate order parameters also from united atom simulations ( However, also this one seems to require Gromacs 4 version (please correct if I am wrong).

      What do you think is the easiest way to calculate order parameters from united atom models nowadays? Can we do this somehow avoiding the installation of old gromacs versions or writing new code?

    28. What about reconstructing hydrogens with another software (e.g. openbabel)? I guess it should be doable to automatize the task like this.

    29. I have now added the results of POPC:POPE mixtures from this discussion into figure:

      The change in beta order parameter of PC seems to be slightly overestimated in CHARMM when mixed with PE.

    30. I discussed the issue of reconstructing hydrogens to UA models with Jon Kapla, who implemented his own code to do this (quite) a while back. His motivation for writing his own code was to take into account the effect of the surrounding angles on the hydrogen location. I was thinking maybe his code is available somewhere, and maybe it could be useful for us. With his permission I am sharing his response below.

      I actually have the analysis code uploaded to GitHub, but it’s been some years since I last used it. If I remember correctly, the actual add hydrogen code is not that complicated, and should be possible to re-implement in e.g. Python (it’s fortran right now) quite easily (I guess). I think the code works, and I have memories of actually using it.

      I won’t have the time to do any coding/testing myself right now, but just ask me if you need help to figure things out. I may (or may not) even remember how that old steam engine of a code works! :)

      Add hydrogen subroutine:

      Subroutine calls:"

      Here’s the very messy overview of the full analysis program:

      The add hydrogen code seems to be left out of that one, but if I remember correctly the input is a ”define atom” as stated in the docs (where only center of mass is an option). Unfortunately I don’t think I have any example input files readily available.

      It may be something like:

      define atom CH2r atom1 atom2 atom3

      where the three helper atoms are used to define the valence angle. It might be that the first helper is where we put the hydrogen...

      Hope it helps!


    31. Thanks for the pointer Markus! That would be great indeed to have our own Python code for reconstructing hydrogens. Amélie and I can have a look in November whether we can convert Jon Kapla's routine to Python.

      Amélie and Patrick

    32. This would be great. For me, the most convenient would be to have a code using the same trajectory reading approach as in the current order parameter code for all atom models:

    33. Hi, thought I could show myself. I'm the guy with the github code Markus mentioned above.

      Samuli, I think this is basically the approach I've had in the fortran code (frame by frame reading and calculation). Add only the hydrogens you need on the fly, right before calling the order parameter (or whatever calculation) routine.

      Feel free to ask if you need help descrambling my old, weakly commented, fortran mess! I don't think I'll be able to help much on the actual coding though, since I'm on parental leave until January. :)

    34. Hi Jon, thanks for the message. I'll send you an e-mail off-list so that we are in contacts.
      @Samuli: sure we'll develop something which will be integrated within your existing code

    35. While protonate isn't part of the newer versions of GROMACS, you could use pdb2gmx to achieve the same thing.

      As for the code in my work, yes it requires older GROMACS versions and doesn't work on the head group either.

    36. Have you used pdb2gmx for this purpose? I have not tried and by looking the help printout, it was not clear to me how this would work. According to the output, it would not read trajectories and I do not know how it would read the *hdb file.

    37. You are right that pdb2gmx doesn't work on trajectories but it will work on the individual frames (I've used this in the past to convert a united-atom bilayer to an all-atom one). The hdb file needs to be in the force field directory with a corresponding rtp entry.

      For a trajectory, a script employing trjconv -sep, multiple pdb2gmx and finally trjcat should do the trick.

  21. Experimental signs for POPS headgroup order parameters are now reported by Tiago Ferreira and added to the manuscript.

    In addition, I modified figures to include the sign information also for PE and PG from simulations. I also set the unknown signs for experimental data of PE and PG according the predictions from simulations. Simulations predict positive signs for alpha carbon for both PE and PG, which is in line with experimental data from PC and PS. However, simulations predict positive signs also for beta order parameters of PE and PG, which would be different than measured for PC and PS. It would be interesting to see from experiments if this prediction is correct.

  22. As Samuli suggested, I have uploaded the trajectories with scaled CaCl2 data for POPC and POPC/POPS Berger membranes from There are low (~100mM) and high (~700mM) salt concentrations included.

    1. Thanks for the data!

      We would still need the simulation only with counterions, i.e. POPC/POPS mixture without any added calcium, but having the charge from PS balanced by monovalent cations. In experiments this would be Sodium, but since Sodium overbinds in Berger simulations, it might be better to use Potassium.

      It would be also useful to run the simulations with Calcium in the presense of monovalent counterions balancing the charge of PS. This would correspond more closely to the experimental situation and we would also see how the presence of the counterions affects the Calcium binding.

      If you happen to have pure POPC and POPS simulations without any additional ions with the same temperature, those would be useful as well.

    2. I ran the system with only Na+ conterions, data is not yet in Zenodo but some files are here:

      I also calculated the delivered order parameters from simulations of POPC/POPS mixtures with CaCl_2:

      These results are yet to be plotted to appropriate figure.

    3. The results are now added also into figure:

      As expected, the calcium binding seems to be overestimated. Also the changes of PS headgroup order parameters are overestimated and alpha response does not seem to be qualitatively correct. Further discussion will be done in the manuscript and more recent posts.

  23. Hi, here are some data for PC/PS mixtures with different amounts of CaCl_2. Simulations were performed using the OPLSAA-compatible lipid model by Maciejewski and Rog and the standard calcium chloride parameters of OPLS:

    1. Thanks for the data! I calculated the order parameters and updated the data into figures.

      Figure comparing experimental order parameters from POPC:POPS mixtures without CaCl suggests that the beta order parameter is well captured by the model, but not alpha:

      As expected, the PC headgroup order parameter decrease due to CaCl are overestimated indicating overestimated Ca2+ binding:
      However, the lowest concentration (100mM) seems pretty good. This requires some more attention. I am suspecting that this might be because of overbinding of counterions preventing Ca2+ binding, but it might be some other reason as well. Also, the concentrations for CHARMM simulations in this figure should be rechecked (discussion above).

      Also the changes of PS order parameters with CaCl concentration are overestimated (except maybe the lowest concentration):

      Simulations of pure POPC and POPS lipid with this model in the same temperature as this data (298K) would be highly useful as well to check the PS order parameters against experiments and the effect of PS to PC and vice versa.

    2. Hi,

      I did plan on running some pure MacRog PS simulations and even converted the CHARMM POPS starting structure into the necessary ordering (they are pretty similar) as per the rest of my pure POPS and DOPS membrane simulations. However when I did this I noticed that the available MacRog POPS topology (obtained from had some problems with it. From looking at some notes I made at the time of doing this, and from having another brief look at the topologies this evening, the things I've spotted are:

      1. The sn-1 chain is oleoyl and the sn-2 chain palmitoyl (so OPPS not POPS)

      2. There aren't any impropers to maintain planarity on either the carboxylate part of the headgroup or the double bond of the oleoyl tail.

      I've been meaning to have another look at this in further detail, correct the problems and setup some PS simulations but a million things have kept getting in the way. This is definitely something that needs looking at though.



    3. Hi,

      this is well spotted, and the POPS/OPPS mistake will likely affect the results.

      However, the lack of impropers in the tail will likely have no effect if the proper dihedrals have high enough barriers and start from the CIS conformation. This seems to be the case for my simulations, where the planarity is well maintained.

    4. OPPS will probably influence the results, but maybe not to a huge degree (, although this should obviously be checked.

      As for the impropers, from just a very (very) quick visual look at the trajectories, the lack of planarity seems to impact more in the head group but I would need to do some proper analysis to be sure of this.

    5. In my experience the order parameters of acyl chain around double bonds depend on dihedrals and MacRog gave reasonable good results in NMRlipids III:

      Assuming that the improper is missing also in POPC, this also suggests that the effect in tail would not be critical.

      For the headgroup, it may be useful to look at the dihedral distributions which are now reported in Figs. 16 and 17 in the current version of the manuscript:

      POPS/OPPS mistake may explain the poor values of MacRog for g3 and g2 in current Fig. 5:
      Especially g2 was significantly better in NMRlipids I.

    6. For POPC, there are impropers around the double bond. I'd have to double check to be sure but IIRC these were in part slightly strange as well. I think there was the typical improper across the C-C=C-C atoms but also another improper including one of the hydrogen atoms of the double bond too. Anyway, what I can definitely remember is POPS missed any impropers while POPC had them (this difference was what drew my attention to this in the first place).

      I'd be really surprised if simply getting the tails the wrong way around (i.e. just having something halfway down the tail different to how it should be) would change much in the glycerol order parameters. I guess the acid test is to switch the tails and run some POPS simulations. Is anyone doing these already? If not I can set some up, although we should decide what to do with the impropers too.

    7. It would also be interesting to plot the oleoyl tail order parameters from the PS simulations (OPPS) without the impropers to see if these are different to the POPC oleoyl tail around the double bond

    8. Thanks for the information regarding MacRog POPC. Oleyl tail order parameters from MacRog OPPS would be quite easy to calculate using the scripts and data available via I have not done it now because the focus of this project is in the headgroup and glycerol backbone.

      I think that the test of switching the tails to POPS would be useful and it would be great if you could do it. We have now quite a lot of data ran with this OPPS MacRog model. I think that we should not touch the impropers within the scope of this project. It is easiest to say that we used the files that developers share, observed these issues and checked if OPPS vs. POPS difference affects the glycerol backbone order parameters. Further testing of impropers can be left for future studies.

    9. Fair enough with regards to not looking at the tail. I will download the sims and do this myself.

      I will setup POPS simulations with the tails from OPPS switched around and everything else the same and run these asap.

    10. Shall I also run some simulations of the mixed PC/PS bilayers too?

    11. Here's data for the same (OP)PS/PC mixtures with different concentrations of KCl:

    12. Hi Thomas,

      I think that pure POPS simulation would be sufficient, at least for now. If it gives some indication that our conclusions on mixtures or ion binding would change, then we can investigate more.

    13. Ok, good. I will setup the POPS simulations today. I have also been setting up some PS/PC mixtures for GROMOS-CKP. I will get these also running today.

    14. Data mentioned by Matti in the above comment

      is now analyzed and added into the manuscript, for further discussion see

    15. The POPS MacRog simulations (so with the tails switched around) have complete 500 ns and I've analysed the order parameters for the final 100 ns of the two simulations. The order parameters are:


      beta -0.01757
      alpha2 -0.02795
      alpha1 -0.00637
      g3_1 -0.20899
      g3_2 -0.09201
      g2 -0.14666
      g1_1 -0.11072
      g1_2 -0.04415


      beta 0.01969
      alpha2 0.01316
      alpha1 0.04411
      g3_1 -0.18245
      g3_2 -0.15399
      g2 -0.13778
      g1_1 -0.15838
      g1_2 -0.00415

      I'll upload the simulations to zenodo shortly.

    16. Thanks again for the data.

      I plotted this together with the previous results for OPPS:

      The difference is larger than I expected. Do you have the same counterions? Maybe there is some dependence on initial configuration, which may have been different between these simulations?

      Anyway, the main conclusion from both data sets would be that the model does not correctly reproduce the PS headgroup structure. Also, I do not think that changing the tail position would be critical for the response of PS headgroup to PC concentration or bound ions. Therefore, we might still consider proceeding only with OPPS simulations of mixtures and systems with ions.

    17. I've finally got around to uploading the data for these MacRog simulations:



    18. I have also analysed and uploaded the simulations for POPS:POPC with the two different GROMOS-CKP variants. As you can see the data is somewhat variable between the repeat simulations performed with different starting velocities, especially for beta/alpha carbons. Below are the PS order parameters from the final 100 ns of the different simulations:

      GROMOS-CKP Berger/Bachar HG v1

      beta 0.00715617
      alpha -0.0625289 0.158183
      g3 -0.193327 -0.291679
      g2 -0.216205
      g1 0.228954 -0.017276

      GROMOS-CKP Berger/Bachar HG v2

      beta -0.0842239
      alpha 0.0545165 0.197031
      g3 -0.283857 -0.351151
      g2 -0.251024
      g1 0.219017 0.134764

      GROMOS-CKP Standard GROMOS HG v1

      beta -0.203155
      alpha 0.170372 -0.0255999
      g3 -0.212513 -0.319887
      g2 -0.182383
      g1 0.284029 0.0287504

      GROMOS-CKP Standard GROMOS HG v2

      beta -0.0521889
      alpha 0.0517796 0.0753356
      g3 -0.284281 -0.287567
      g2 -0.259306
      g1 0.281642 0.0104549

      The simulation files can be found at:

      I'm re-running the analysis on more of the trajectories to see if that alleviates some of the differences between simulations with the different starting velocities. If not, I'll probably run some more repeats (I do also have some simulations with reaction-field that have already completed).

      Finally, I have also calculated the PC order parameters in these simulations too. If they are needed let me know and I will post them.

    19. Some values in "GROMOS-CKP Standard GROMOS HG v1" looks a bit strange, I will use only "GROMOS-CKP Standard GROMOS HG v2" for now.

    20. I have now transferred the continuation of the discussion about above reported GROMOS-CKP simulations here:

  24. In order to organise the material, I have now divided the manuscript into two parts. The first one contains data for PS lipids ( and the other one for PE and PG lipids ( I think that we should progress these separately to make things more clear and easy. I will soon write a new blog post about this.

  25. Hi,

    I have started simulations of DOPS with Amber Lipid 17 force field; 128 DOPS lipids at 303K, following Tom's protocol. Na+ is the counter ion in the system, and I am running two replicates each using Amber ff99SB and Joung-Cheatham ions. It could be interesting to look at how ion parameters also effect to results. POPS at 298K simulations will also start shortly.

    Results will be ready in the early weeks of 2018.

    1. Excellent!

      I have heard rumors about Amber Lipid17, but I have not seen any publication yet. I guess that the parameters are, however, available somewhere as you were able to start simulations? Do you know anything about the publication or other reports about Lipid17?

      I am looking forward to the results.

    2. Hi,

      You are right, there is no specific publication named Amber Lipid 17. Amber people want us to use this reference for Lipid17:

      Madej, B., Gould, I.R., Walker, R.C. A Parameterization of Cholesterol for Mixed Lipid Bilayer Simulation within the Amber Lipid14 Force Field. J. Phys. Chem. B. 2015.
      DOI: 10.1021/acs.jpcb.5b04924

      where they added PS and cholesterol to Amber Lipid 14 forcefield and named the new FF as Amber Lipid 17.

    3. That question came out Today at Amber mail list, I'm copying the response of Ross Walker here:

      "Unfortunately due to some issues beyond the control of the various authors the LIPID 17 paper isn't yet complete although we are hoping to have something submitted in the next couple of months. In the meantime you can reference it as:
      Gould, I.R., Skjevik A.A., Dickson, C.J., Madej, B.D., Walker, R.C., "Lipid17: A Comprehensive AMBER Force Field for the Simulation of Zwitterionic and Anionic Lipids", 2018, in prep."

  26. Hello,

    I've obtained and uploaded the trajectories of the pure DOPS with Na+ counterions, using Amber Lipid17 and Joung-Cheetham/ff99SB ion force fields. You can access the files using the below links:

    They contain the Amber topology file, Amber input file for the simulation parameters, four consecutive equilibration runs, and one production run, each of which 100ns long with 10ps sampling rate. I performed two replicates, and I uploaded all of them.

    Simulations are done using Amber16 engine with GPU's.

    I performed an initial analysis of the trajectories and calculated the head-group and glycerol backbone deuterium order parameters using AmberTools17. I am posting one output here to give you all an idea about how they look.

    -0.1246016 0.0394291 -0.1731807 0.0411584
    -0.1532849 0.0420252 -0.1532849 0.0420252
    0.2571898 0.0231445 -0.0197184 0.0276169
    0.2138903 0.0256339 0.2138903 0.0256339
    0.2788865 0.0219033 0.1659704 0.0280705

    Above are the values for DOPS at 303K, with ff99 ions. The order is: alpha, beta, g1, g2, and g3. First and third row are the order parameters for the given carbon, second and fourth columns are the standard deviations of the previous column, calculated over the trajectory.

    It looks like absolute values of the order-parameters match quite OK with the results from other forcefields. Only problem is with the sign; I'm not sure if there is a bug in AmberTools, or I did something wrong in the calculations. Anyway, I think I will be very good if someone can calculate the order parameters with their own scripts to cross-check. If needed, I can provide order parameters from all simulations.

    1. Thanks a lot for the data!

      I ran the analysis for your data using our current script. The results and scripts are here:

      The values in OrdParsDOPS.dat are exactly the same as from AmberTools17, but the signs are opposite. Thus, it seems that AmberTools gives -S_CH. This is common since traditionally the focus has been on the absolute values of acyl chain order parameters, which are negative.

      Based on these results, it seems that the Lipid17 model do not perform very well for the headgroup when the order parameter signs are taken into account.

    2. I confirm the opposite signs, I did the same analysis independently too. I found a bit differenet values (Used Joung-Cheetham/ff99SB ions), but close to what is reported:

      OP Name mean std err.est. confidence_interval 95% (min, max)
      alpha2 0.1713 0.4640 0.0008 (0.17051999379859795, 0.17212767640433774)
      alpha1 0.0951 0.4585 0.0008 (0.094302361808018945, 0.0958909876633494)
      g2_1 -0.2201 0.3245 0.0006 (-0.22069440738359261, -0.21957018008784165)
      g1_1 -0.2598 0.2957 0.0005 (-0.2603586368197231, -0.2593341715161927)
      g3_1 -0.1464 0.3570 0.0006 (-0.14704359922346216, -0.14580669281037983)
      g1_2 0.0011 0.3982 0.0007 (0.00041245785334314698, 0.0017919940692083336)
      g3_2 -0.2461 0.2964 0.0005 (-0.24659552218697897, -0.24556869376462301)
      beta 0.1489 0.4738 0.0008 (0.14812595937439982, 0.14976740440765673)
      alpha2 0.1764 0.4803 0.0008 (0.17553732231672917, 0.17720131712059575)
      alpha1 0.1303 0.4595 0.0008 (0.12949260661609996, 0.13108467694791337)
      g2_1 -0.2204 0.3127 0.0005 (-0.22096529567694945, -0.21988196022546025)
      g1_1 -0.2481 0.3054 0.0005 (-0.24863780883344752, -0.24757971223048816)
      g3_1 -0.1831 0.3378 0.0006 (-0.18365371816587137, -0.18248322061505842)
      g1_2 0.0049 0.4001 0.0007 (0.0041958883181644402, 0.0055822239328987113)
      g3_2 -0.2328 0.3231 0.0006 (-0.23339062720262299, -0.2322710276848492)
      beta 0.1743 0.4805 0.0008 (0.17348336947852222, 0.17514802129094276)

      Raises a question - is the sign for beta right?

    3. Thanks for this as well. Which script did you exactly use for these? I think that we need to do something about the error estimates of these scripts. I think that the best way is to calculate order parameters separately for each lipid and then use the error of the mean, as done in NMRlipids I and II.

      The experimental order parameters with the signs are here:
      and they are plotted with some simulation data in here:
      so it seems that the sign of beta in Lipid17 is not correct and the results for alpha are also not very close to experiments.

    4. The corresponding trajectories for POPS (at 298K) with ff99 ions and with Joung–Cheatham ions are now also available on Zenodo.

    5. Would it be possible to report also the counterion number densities along membrane normal from Lipid17 POPS simulations to be added in this plot:

      Also the area per lipid values would be useful.

    6. Hi,

      Na+ number densities for POPS simulations (with ff99 and JC ions) are in this folder:

      Area per lipid values for all Lipid17 systems are:

      63.6521 ff99-Ions
      59.3167 JC-Ions

      57.3274 ff99-Ions
      52.7930 JC-Ions

      There are significant differences between different ion models, as JC-ions produce lower APL values. This also might be related why dihedral distributions look frozen.

    7. Thanks for the data.

      I added the density profiles and area per molecules also from other models to the plot:

      All the models, except CKP, give smaller area per molecule than the experimental value 62.7 A^2 ( CKP values are in good agreement with experiment. Maybe this was used as a fitting parameter in parametrization?

      There has been suggestions in many simulation publications that the small areas in PS bilayers would arise from counterion binding and concomitant screening of the repulsion between headgroup charges. When looking at the data, I can see some correlations between binding affinity and area per molecule. However, the correlation does not seem very strong when comparing across the models.

    8. No, the APL wasn't used for any sort of parameter adjustment/refinement with GROMOS-CKP. In fact I just simply took the standard CKP model for most of the lipid and combined it with the typical GROMOS protein force field parameters for the PS head group (using two different schemes for the charges).

      As you say, it works very well in terms of general membrane order, especially when compared to the other force fields. I can provide full plots of the APL, etc. if needed from my simulations. I've done quite a bit of this analysis already for a planned manuscript I've been slowly working on with all these PS simulations (still work in progress!)

  27. Hello,

    I have finished the simulations of POPS:POPC mixtures (1:5) using Amber Lipid17 forecfield at various KCl and NaCl concentrations (ff99 ion model used). All trajectories are 200ns long.

    The paths to the trajectories are:

    1. Thanks for the data!

      Would it be possible to have also simulations of the same POPS:POPC mixtures (1:5) containing only counterions? We need these for the reference systems when looking the effect of added ions. Also simulations with CaCl2+conterions would be highly interesting.

      How are the reported concentrations calculated?

    2. Hi Samuli,

      The concentrations are calculated in the same way described in Matti's previous post ( My systems are the exact replicates of him, just with Lipid17 forcefield.

      I will start the simulations of POPS:POPC mixtures containing only counterions today. I'll post the data when they are ready.

      With CaCl2+counterions, I need to check the available ion parameters in Amber that would be suitable with Lipid17. I'll start those simulations next week most likely.

    3. Hi!
      I'm currently working with the supposed-to-become Lipid17 (or 18?) POPS, but within Gromacs. I have some results with PC-PS mixtures 5-1 at various salt concentrations. Once they are finished with a sufficient length, I'll happily share them.

      So far, the model looks quite promising despite the inaccuracies in the actual order parameters for pure POPS. However, I use different models for ions.

    4. Hi Batuhan,

      I have a question on these. Is the atom naming of POPS headgroup the same as DOPS, i.e., can we use the same definition file that we used to analyse DOPS from AMBER trajectories with our python code? I mean this definition file:

      Since I am not very familiar with the AMBER filetypes, I was not able to check this myself easily.

      Let me know if there is already some data for only counterion systems already available.

    5. I should also somehow find out the headgroup atom names of POPC in the AMBER files.

    6. Hi Samuli,

      As it is, the script we have won't be able to distinguish between different lipids with the same type of head group if we have a mixture of them. Other than that, in Amber all head groups, respective of the tail, have the same naming. I am working now on calculating the order parameters from all the simulations I have with Amber ff.

      At the moment, I have the simulations of POPC:POPS mixtures with NaCl, KCl, and CaCl2 counterions as well. I haven't uploaded them to Zenodo yet, but in few days I will.

      For the atom names in POPC and POPS, I will also upload a pdb file.

    7. I did some analysis already on the order parameters, see and corresponding folders for other systems.

      When we have the system with counterions only, I will add the results into the manuscript.

      The pdb files would be useful indeed.

    8. Hello,

      I have uploaded the trajectories for POPC:POPS mixtures with CaCl2, NaCl, and KCl counterions to Zenodo. The corresponding links are:

      I will also upload the order parameters, density profiles and dihedral distributions shortly.

    9. Thanks again for the data!

      I have now incorporated the results into the plot showing results for PC:PS mixtures with monovalent ions:

      The figure could be more clear, but it seems to me that the changes of order parameters as a function of K+ or Na+ would be in rough agreement with experiments.

      The ion densities calculated from these simulations for this figure would be highly useful:

      Do you have simulations with different CaCl concentrations running/available? Those would be highly interesting as well.

    10. Simulations with various concentration of CaCl2 are now available through my recent commit (requested to be pulled from my forked repo)

      Gromacs 2018 was used, parameters for ions are by Dang as in the following repository:

    11. Hello,
      Sorry for the long silence. I personally have not done any simulations at different CaCl2 concentrations. I will start these simulations this week, and should get the results in next week.

    12. Simulations delivered by Pepa are now included in the manuscript (

      The upcoming simulations by Batuhan are now mentioned in the list of current issues in the new post:

    13. Hello,
      The trajectories for POPC/POPS mixture at various CaCl2 concentrations (500mm,1000mm,2000mm) are uploaded to

      I will upload the order parameters and density profiles to GitHub when I obtain the results.

    14. Looking forward the results.

      Could you also deliver the ion density profiles of simulations with monovalent ions reported above?

    15. Hello,

      I've realized that the number of ions in Amber Lipid17 simulations, namely for KCl, NaCl, and CaCl2 at 500mm, 1000mm, 2000mm, 3000mm, and 4000mm concentrations is slightly larger than what it is supposed to be. It's due to an error I made when calculating the required number of ions. To be consistent with other simulations we have, I have fixed the number of ions, and rerunning these mentioned simulations. I apologize for this mistake. Hopefully in two weeks I'll get the correct trajectories, and upload the files to zenodo and github.

    16. I do not think that the concentrations have to be exactly the same in different models. Once you report the amount of monovalent ions in the current data ( and, I did not find the numbers of counterions anywhere and I did not know how to check it from Amber files), we can calculate and report the correct concentrations in the table II and in the plots. I believe that the conclusions are so clear that the small difference in concentrations between the models should not matter. Therefore, I think that the new simulations are probably not necessary.


Please sign in before writing your comment.