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

  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

  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:

  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.


Please sign in before writing your comment.