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.

  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. 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.

    8. 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

    9. 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?

    10. 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

    11. 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.


    12. Hi all,

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


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

  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:

  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.

  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.



  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.


Please sign in before writing your comment.