NMRlipids IV: Current status and reorganization of the manuscript

First of all, I want to thank again all the contributors for delivering significant amount of useful data about PE, PG and PS headgroups for the NMRlipids IV project.

In addition to simulation results, also experimental signs of headgroup and glycerol backbone order parameters were contributed for the POPS lipid. However, experimental signs for PE and PG are not yet known, which makes the comparison between simulations and experiments more ambiguous. Therefore, I have divided the manuscript into two parts. The first one contains results only from systems with PS lipids and the other contains results from systems with PE and PG lipids. This should also ease the management and completion of the projects. The manuscript about PS lipid systems should be more straightforward to finish and I have started to compile it towards a submittable version. Some of the current results and the most important open tasks are listed below.

Headgroup & glycerol backbone structures of PS lipid bilayers

Since the order parameter signs are known for the PS lipid headgroup and glycerol backbone, we can perform similar comparison between simulations and experiments as was done in NMRlipids I for the PC lipids. This is shown Fig. 1; the subjective quality assessment is also available. The tested models perform generally less well than the PC lipid models discussed in the NMRlipids I publication. Many relevant contributions have already been made to give structural insight to the differences, however, I think this deserves more attention (see the ToDo list below).

Figure 1: Headgroup and glycerol backbone order parameters of the PS lipids from different simulation models and experiments. 

Headgroups in mixtures of PS and PC lipids

Lipid mixtures are biologically relevant, and in fact experimental data for ion binding to PS lipid bilayers seems to be available only for PS/PC mixtures. Therefore, we also address the mutual interaction between PS and PC lipid headgroups. Figure 2 shows how the addition of POPS changes the order parameters of POPC (left column), as well as the changes of POPS order parameters with increasing amount of POPC (right column). It seems that the tested simulation models do not reproduce the interactions between PC and PS headgroups very well. However, the inaccuracies in counterion binding may also disturb the lipid headgroup structures when the amount on PS (and counterions) is increasing. I think that we need data with a few more different models before drawing general conclusions, and possibly also some data with different counterion concentrations (see the ToDo list below).

Figure 2: Headgroup order parameters from PC:PS mixtures from different simulation models and experiments. Left panel shows the PC headgroup order parameters and right panel shows the PS headgroup order parameters.

Interactions between cations and lipid bilayers containing PS

As already discussed in the opening post of NMRlipids IV project, molecular electrometer experiments show that the presence of negatively charged lipids enhance cation binding in lipid bilayers. The available experimental data for the lipid headgroup order parameters of POPC:POPS (5:1) mixture as a function of CaCl₂ concentration are shown in Fig. 3 together with the simulations ran with the MacRog model. In line with the NMRlipids II publication, the PC headgroup order parameters' decrease with increasing CaCl₂ concentration is overestimated in the MacRog simulations, indicating overestimated binding of Ca²⁺ to the lipid bilayer. It should, however, be noted that the point with the lowest CaCl₂ concentration is in better agreement with the experiments. Potential explanation could be an overestimated screening effect by the overbound counterions, but this requires further analysis. Order parameters of the POPS headgroup rapidly increase or decrease even with low CaCl₂ concentration and reach almost a plateau value above 100 mM. Also the order parameter changes of the POPS headgroup are overestimated in the MacRog simulations and a qualitative agreement for the alpha carbon is unclear. Simulation data of PC:PS mixtures with different CaCl₂ concentrations is required also from other than MacRog model for more general conclusions (see the ToDo list below).

Figure 3: Headgroup order parameters from PC:PS (5:1) mixtures with different CaCl₂ concentrations from the MacRog simulation model and experiments. Left column shows the PC headgroup order parameters and right column the POPS headgroup order parameters.

ToDo list

1. Structural relevance of the observed order parameter differences between different simulation models and experiments are now analyzed using dihedral distributions (see Fig. 13 in the manuscript) and we have pictures of the sampled conformations of glycerol backbone and phosphate. I think that we should apply the tools contributed by Pavel Buslaev to calculate dihedral distributions and to visualize the structural sampling also for other contributed models. Among these data we should then select the parts giving the best representation for structural differences related to order parameters. To do this in practice, we need more simulation trajectories available in Zenodo.
2. We need data for PC:PS mixtures without additional ions from few more simulation models. Some data from Berger model has already been delivered by Lukasz Cwiklik (see also data with calcium), but we still need simulations with counterions only. Simulations of PC:PS mixtures with the Slipids and CKP models may also be useful.
3. Simulations of POPC:POPS mixtures with different CaCl₂ concentrations are needed also from other force fields than MacRog. The above mentioned dataset with Berger force field complemented with simulations containing only counterions would be highly useful. I think that simulations with CHARMM36 is a must. Also simulations with the Slipids and CKP force fields may be worth of doing. Based on NMRlipids II, it is expected that cations will overbind to lipid bilayers in these simulations. However, the behavior of PS headgroup as a function of salt concentration in different force fields with respect to experiments is not yet known.

Quantifying the effect of bound charge on headgroup order parameters

[UPDATED 8.12.2017]

Electrometer concept is used in NMRlipids II and IV to measure the amount of bound ions in lipid bilayers. It is based on empirical observations by Seelig et al. that the PC lipid headgroup order parameters for alpha and beta carbons depend linearly in bound charge. In NMRlipids II project the concept was observed to be qualitatively valid also in simulations, as seen in Fig. 1.

₂

Fig 1: Change of order parameters as a function of bound charge analyzed from simulations with various ions in NMRlipids II publication (Fig. 3.)

However, the comparison of this data with experiments was not straightforward, because it simultaneously depends on the amount of bound ions, definition of bound ion and sensitivity of the headgroup order parameters to bound ions. The issue was discussed during NMRlipids II project and also in the supplementary information of the publication (section 3), but simulations with cationic surfactants were not performed to quantify the sensitivity of the headgroup response to bound charge. 

Because the issue significantly complicates the usage of electrometer concept in simulations, as also seen now in NMRlipids IV project, I performed a simulation of PC lipid bilayer mixed with different mole fractions of cationic surfactants (more specifically dihexadecyldimethylammonium bromide, C₁₂C₁₆⁺N₂CBr^-). The advantage of such system is that essentially all the ions (i.e. charged surfactants) can be assumed to be bound in bilayer, thus the amount of bound charge is known exactly. Thus, the system can be used to quantify the lipid headgroup sensitivity to bound charge.

The headgroup order parameters from Lipid14 simulations and experiments as a function of cationic surfactant are shown in Fig. 2.

Fig 2: Headgroup order parameter changes as function of cationic surfactant from simulations with Lipid14 (files available at 0.1, 0.2,0.3,0.42 and 0.5), CHARMM36 and experiments.

The results show that the headgroup order parameter response to bound charge is approximately linear also in simulations. However, headgroup is slightly too sensitive to the bound charge in Lipid14 model. This indicates that some of the overestimated order parameter decrease for Lipid14 with CaCl₂ concentration in NMRlipids II publication may be due to the headgroup response on bound ions instead of overestimated binding affinity.

My feeling is that we need to do such test, at least, also for CHARMM36 model, for which data about Ca²⁺ binding in negatively charged lipid bilayer was recently reported.

[UPDATE 8.12.2017] CHARMM36 results added in Fig. 2 show better agreement with experiments. This has to taken into account when using headgroup order parameters to compare binding affinity between simulations and experiments. 

NMRlipids III: Preliminary version of the manuscript

 [updated 19.9.2018]

I have now updated the progress made after NMRlipids III: Preliminary observations post into the manuscript. 

After extensive discussion about correct parameters to run CHARMM36 simulations it was concluded that the results from Gromacs 5 are consistent with other simulation packages and literature. Thus, the results from CHARMM gui parameters simulated with Gromacs 5 are currently used in the manuscript. However, there might still be some issues with parameters given by CHARMM gui for systems with cholesterol.

Current comparison for acyl chain order parameters between different simulation models and experiments is shown in Fig. 1. The main conclusion is that cholesterol ordering effect is overestimated in Berger/Holtje and MacRog models, while slight overestimation is observed also in CHARMM36 and Slipid models. The significance of overestimation in CHARMM36 and Slipid is yet to be analyzed (see Things to do list below).

FIG. 1: Order parameters from simulations and experiments for acyl chains of 1-palmitoyl-2-oleoylphosphatidylcholine (POPC) with and without cholesterol.

Final comparison between form factors from simulations and experiments is yet to be done, but area per PC headgroups and per total number of molecules from different simulations are shown in Fig. 2. The difference between models do not seem excessively large, but Berger/Holtje and MacRog models predict strongest effect due to cholesterol in line with order parameter results.

FIG. 2: Area per molecules calculated from different simulation models as a function of cholesterol concentration. Solid marks are area per total amount of molecules (chol+PC) and empty marks are area per PC headgroups. Top figure shows absolute values and bottom figure shows changes respect to pure lipid system.

Things to be done:

• We need to analyze if overestimation of order parameters with cholesterol in CHARMM36 and Slipids is significant. There are at least two issues to be clarified: 
  1) Is the quantitative accuracy of experimental order parameters same for acyl chains as for glycerol backbone and headgroup (0.02) discussed in NMRlipids I ?
  2) Can the overestimation of order parameters in simulations be explained by missing undulations, as sometimes suggested? [update 19.9.2018] As discussed in the GitHub issue linked here, the undulations does not seem to make significant difference.
  
• Lipid 14 results from model by Madej et al. should be somehow included.
• Comparison of form factors between simulations and experiments should be finished.
  
  
  
  

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

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 CaCl₂ 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 Ca²⁺ increases when negatively charged lipids are present in bilayers.

Fig. 6 Changes in the PC headgroup order parameter as a function of CaCl₂ 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) Ca²⁺ 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 Ca²⁺ 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 Ca²⁺ also means increased amounts of bound Ca²⁺ at neutral and charged lipids;

(iv) the actual binding step can be described by a Langmuir adsorption isotherm with a 1_lipid_:_1_Ca²⁺ stoichiometry, provided the interfacial concentration C_M 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.

Future of NMRlipids project

As written in the beginning of the project, "[t]he ultimate goal of this blog is to find an atomistic (preferably united-atom) force field that reproduces the experimental properties discussed in the manuscript [headgroup and glycerol backbone order parameters and their responses to ions, dehydration and cholesterol]. Naturally the optimal situation would be that some of the already available force fields would fulfill this goal. If this, however, turns out not to be the case, the goal will be to find the appropriate modifications."

In NMRlipids projects I and II it turned out that none of the available force fields fully satisfied these goals. Especially the ion binding affinities and detailed structure of the glycerol backbone and headgroup posed major challenges for the current force fields. On the other hand, qualitative response to dehydration and bound charge were well reproduced by all the models.

The attempts to improve force field parameters by Antti Lamberg in March 2014 revealed the importance of signs and stereospecifity of the order parameters. The necessary details are now reviewed in NMRlipids V publication and the development process can be resumed using better defined experimental numbers. Inspired by this, we have made concrete plans to expand the NMRlipids project towards systematic force field improvement and have build a prototype of an automatic force field quality assessment tool. However, we are still looking for ways to organize appropriate human resources to run this extension of the project.

Samuli Ollila has now received IOCB fellow funding from Institute of Organic chemistry and Biochemistry in Prague, Czech Rebublic to focus on ion interactions with zwitterionic and charged lipid bilayers in Pavel Jungwirth's group. For this reason, the main focus of NMRlipids IV will be in charged lipids and their interactions with ions (new post will follow soon). Also NMRlipids VI, to develop PC lipid model with correct Ca²⁺ and Na⁺ binding behaviour, will be lauched in the near future.