I have started to assemble the NMRlipids III manuscript. The project is initiated from the discussion started by Peter Heftberger (see Data Contributions page, discussion from December 1, 2014). Heftberger et al. compared x-ray scattering form factors with results from CHARMM36 simulations and concluded that the results did not quite agree (see the report). This was slightly surprising since simulation literature typically reports good agreement with experiments for form factor data, also for CHARMM36 with cholesterol. In the subsequent discussion we decided to take a simple two component lipid bilayer to check how accurately lipid mixtures can be described by atomistic resolution molecular dynamics simulations. For such system we chose the POPC/cholesterol mixtures up to 60 mol% of cholesterol. Extensive NMR data for this systems is already published, extensive simulation data sets are already collected in NMRlipids I and extensive scattering form factor data was reported in this blog by Georg Pabst et al. (July 3, 2015, comment in About page).
I have now started to collect this data in the manuscript (available from GitHub according to the current workflow of the project). All the data is not there yet, but some preliminary results can be already discussed. Some results for order parameters and form factors, and preliminary observations are shown here:
 | |
| Figure 1: Acyl chain order parameters from experiments and different simulations for pure POPC and POPC with 50mol% of cholesterol. |
 | ||||||||||||||||
| Figure 2: Form factors from experiments and simulations with different cholesterol concentrations in POPC bilayer. The Experimental form factor amplitudes are not scaled to match with simulations, thus mainly the locations of minima should be compared. |
Preliminary observations:
CHARMM36: The order parameters are too large in both simulations, with and without cholesterol. The form factor gives minima at too small q values. This is not expected from the parametrization publications. However, the form factor result is pretty similar to the one report reported by Heftberger et al. (simulations ran with NAMD). More detailed discussion about reproducing some CHARMM36 results is in GitHub issue. [Update 19.10.2016] as shown in the GitHub issue, the overestimation is largely (but not completely) due to the usage of older version of Gromacs in the original results. The results from Gromacs 5 are now used in the manuscript.
Berger/Höltje: The order parameters in pure POPC are in very good agreement with experiments, as often reported in the literature. However, form factor has minima at too high q-values. With 50mol% of cholesterol the order parameters are too high, as reported before. However, the form factor minima are surprisingly good with 50mol% of cholesterol. My preliminary interpretation is that the model has too large area per molecule in pure POPC and too strong condensing effect of cholesterol. This leads to good agreement in form factor with large cholesterol concentration but order parameters are then too large. This would also indicate that the relation between order parameters and packing density would not be right in this model. However, more careful analysis is needed.
MacRog: Order parameters are not in good agreement in neither of the systems, with or without cholesterol. Form factors still to be calculated.
Things to be done:
- I would like to plot the electron densities from simulations together with the results from the SDP model used for the scattering data. The SDP profiles are available in the delivered data set, but only in *eps format. Would it be possible to deliver these also in ASCII format?
- In addition to the simulation data we already have, at least Slipid and Lipid14 models have compatible cholesterol models. At least some data from these models is almost necessary for this project. The corresponding simulations are already ran for Lipid14 and we could take the order parameters from the publication. However, for form factors we would need the simulation trajectories (or someone delivering the results). [Update 19.10.2016] Slipids data are already delivered. Lipid14 data is in progress.
- There is also experimental form factor data available for other binary mixtures than POPC/chol. The data for DOPC/CHOL and DPPC/CHOL mixtures was delivered by Heftberger et al. (see also report) and there seems to be more in the literature. If someone delivers corresponding simulation data, the comparison would be a good complement for the extensive POPC/chol set.
- The order parameters for POPC (preferably in 303K or close) from CHARMM36 simulation ran with NAMD or CHARMM would be useful to check if they are the same as from Gromacs. More discussion in GitHub issue. [Update 19.10.2016] NAMD data and much more is delivered throught the discussion in GitHub. The conclusion seems to be that Gromacs 5 gives more or less consistent results with other codes.
Hi,
ReplyDeleteI haven't done a huge amount of cholesterol simulations before (have mostly worked either simulating single component PC/PE/PS/etc. or mixed bacterial membranes) but I'm getting more into this area so I am keen to help out here. As you allude to in point 2, there are lots of force fields out there that have cholesterol parameters available. In fact, there are only a few force field I can think of where there are PC parameters but not compatible cholesterol ones available. The simulations I was planning on running were going to be some GROMOS based force fields to test the combination of the cholesterol parameters that are available as a manual entry on the ATB website with a couple of GROMOS based phospholipid models. Happy to share this data if and when I get around to these simulations. Is there a starting structure that the current simulations are all using or are the AA and UA ones different? If there is a consistent one, I will use this.
As an aside, you mention that you have looked at the Berger/Höltje combination. It could also be worth looking at http://pubs.acs.org/doi/abs/10.1021/ja8095224 where they reparameterised the charges of this cholesterol model. I'm sure there would be other such modifications along these lines out there but this was the one that stuck in my mind.
As for CHARMM36, even if you do everything correctly (ff conversion, water model, cut-offs, etc), the membranes will be more ordered in GROMACS than in either CHARMM or NAMD (see http://pubs.acs.org/doi/abs/10.1021/acs.jctc.5b00935). As far as I can tell, no-one seems to know why yet. I would recommend either using NAMD/CHARMM to do the simulations or just accepting the issue and use GROMACS but when doing so make sure that for pure membranes you use the switching off of the forces at 0.8 nm (as was originally used for the CHARMM36 membrane parameterisation). I have some simulations that I did in NAMD of PC membranes but I will have to check as I think it might just have been DPPC and not POPC. Happy to share the DPPC ones if they will help.
Cheers
Tom
Excellent!
DeleteCan you specify what do you mean with ATB website?
The current simulations with cholesterol are listed in Table 3 in NMRlipids I publication (http://dx.doi.org/10.1021/acs.jpcb.5b04878). The details of starting structures are given in supplementary. If I remember correctly, the initial structures for different models are not the same but there may be some correlation. So I do not think that there is consistent one and as we have discussed in GitHub (https://github.com/NMRLipids/lipid_ionINTERACTION/issues/8) the effect of this is not yet clear.
The parametrization you mention may be worth of considering. I will read the paper.
I find the NAMD/CHARMM vs. Gromacs result differences quite disturbing. With Gromacs 5 there are some instructions how to set cut-off parameters to be combatible with CHARMM model:
http://www.gromacs.org/Documentation/Terminology/Force_Fields/CHARMM
However, I think that this does not fix the problem. I will read the paper by Lee et al. more carefully, but I think that performing NAMD vs. Gromacs comparison also by ourself is probably reasonable.
By ATB I meant the automated topology builder of Alan Mark's group which generates GROMOS based small molecule topologies. They also have manual (i.e. not automatically generated) entries for some molecules. The cholesterol one (https://atb.uq.edu.au/molecule.py?molid=1731) has been used in several publications by their group who also developed the GROMOS 54A7/GROMOS 53A6L PC lipid parameters(the Poger et al. ones). I was planning on testing this Cholesterol with some of the other GROMOS based lipid force fields. I have been setting up some of these simulations this week (DPPC/Chol and POPC/Chol).
DeleteI also have just finished running (for a different project) some GROMOS43A1-S3 DPPC/Cholesterol simulations (only at 50% at the moment, 2 simulations, both 500 ns). I haven't done POPC/Cholesterol with this force field due to the major discrepancies with the deuterium order parameters around the double bond of the oleoyl tail. If these DPPC/Chol 43A1-S3 simulations are useful, happy to share.
Agreed that it is very worrying regarding GROMACS with CHARMM36 membranes. I have seen the same of slightly more ordered membranes with other non-PC lipids (PE/PS and some glycolipids) and also with the united-atom CHARMM36 tails, in terms of comparing my simulation results to published ones. It'd be very good to work out why this is happening.
Here's some data from membranes with varying cholesterol content ran with Slipids. It's essentially the same membrane that I've uploaded earlier with the scaled ion models, I just ran it without any ions. Data can be found here: DOI:10.5281/zenodo.60607
ReplyDeleteThanks. I did run the analysis for the data (currently I use this script https://github.com/NMRLipids/MATCH/blob/master/scratch/report/makeREPORT.sh and the data is found from this folder https://github.com/NMRLipids/MATCH/tree/master/Data/Lipid_Bilayers)
DeleteI did also plot some of the data in the current figures (Figs. 1-3 in the manuscript https://github.com/NMRLipids/NmrLipidsCholXray/blob/master/MANUSCRIPT/manuscript.pdf).
I cannot make a normal order parameter figure (like Fig. 1) with all concentrations since the data from different sources (models and experiments) are not with exactly same concentrations. I have tried to show the cholesterol effect by plotting each order parameter as a function of cholesterol concentration (see current Fig. 2 for the results for sn-1 chain). From this plot it seems that the Slipids has clearly most realistic response to cholesterol in acyl chain region. Also the form factors in Fig. 2 look pretty good, but these results need to be double checked (there is also something weird in the experimental data with 30 mol%).
I think that it would be good to have a 50mol% POPC/cholesterol mixture from this model, expecially because the results look quite good with this data. 50mol% mixture would have maximum amount of lipids and cholesterol and the quality of these interactions would be best seen in there.
Hi, I have the 50mol% POPC/chol for SLIPIDS at 283, 298 and 308 K. I will post the data and upload the trajectory to Zenodo this weekend.
DeleteI also built this (and 40mol%) already, and they are currently being simulated. At least we'll see if the results agree. I constructed mine using Charmm-GUI.
DeleteHere's the data of the order parameters https://dl.dropboxusercontent.com/u/11027999/nmrlipids.tar from my system (50mol%) built with CHARMM-GUI too, 100ns trajectories.
DeleteThanks for the data! I added this in the plot: https://github.com/NMRLipids/NmrLipidsCholXray/blob/master/FIGS/OrderParametersCHOL-eps-converted-to.pdf
DeleteSlipids seems too ordered as well in 50/50 mixture. I did not update the other plot yet.
This comment has been removed by the author.
DeleteMy Slipids data with 40% and 50% cholesterol are here:
DeleteDOI:10.5281/zenodo.154346
They complement the set (0%,10%,20%,30%) I posted earlier:
DOI:10.5281/zenodo.60607
Btw. I used the same (slightly reordered) .gro files as for the Charmm36 runs in:
(DOI:10.5281/zenodo.61649, DOI:10.5281/zenodo.61992)
Fernando: It seems to me that the data in https://zenodo.org/record/62026 has less cholesterol than 50%, please check this.
DeleteHere is the data for 30%:
Deletehttps://doi.org/10.5281/zenodo.62026
and for 50%:
https://doi.org/10.5281/zenodo.159434
I have now analysed these data sets and included some of the results in the manuscript. I also added the information in the table listing simulated systems. I have few questions:
DeleteFernando:
Could you share the information about number of different molecules (lipid, cholesterol and water) in the systems?
You probably have the simulation of pure POPC with Slipids in T298K as well. Would it be possible to share this data as well? Matti's simulation set is in 310K while the experiments are in 300K. Based on current data, there seems to some temperature dependency. With pure POPC in 298K we could quantify this effect.
Matti:
In Zenodo data sets you write that the simulation trajectories are 100ns. Did you ran some equilibration before this 100ns? The same question applies also on the shared CHARMM36 data.
The systems with 0–30% cholesterol were first simulated for 200 ns with 130 mM NaCl. After this the ions were removed and the systems were simulated for 100 ns.
DeleteFor the systems with 40–50% cholesterol, the standard equilibration steps provided in the README file from CHARMM-GUI (total of 9 ns equilibration) was performed before the 100 ns simulation.
Let me know if you think the simulations need to be extended.
Data sets updated with the requested info:
Delete00%: https://doi.org/10.5281/zenodo.166034
30%: https://doi.org/10.5281/zenodo.62026
50%: https://doi.org/10.5281/zenodo.159434
Please let me know if you need anything else.
Hi Samuli,
ReplyDeleteI also have performed atomistic MD simulations of DPPC (30% CHOL) and DOPC (20% CHOL) lipid bilayers in the presence of cholesterol at different temperatures with 150 mM NaCl using the Slipids force field in Gromacs. Please let me know if you think that these systems could be relevant for the manuscript, if so I can analyse them and provide additional results.
Thank you very much for your consideration,
Andrea
Hi Andrea,
Deletethank you for letting us know about your data. However, I think that (at least for now) we should limit this manuscript to POPC/cholesterol mixtures without ions for which we have the corresponding data from the NMR and scattering. However, currently it seems that Slipids describes cholesterol-lipid interactions very well, it may be interesting to test other mixtures later on. We will keep this data on mind for the possible usage in the future. Of course, you can always upload your data in Zenodo to make it useful for anyone.
Hi Samuli,
Deleteyou are welcome, I thank you very much for letting me know about the objective of the manuscript on POPC/cholesterol mixtures without ions and for sharing the interesting results obtained with the Slipids force field. Ok, it sounds good that you will consider my data in the future. I will share your post with my advisor to see if the data can be uploaded in Zenodo.
Thank you very much again for your consideration,
Andrea
I red the original post again and realized that we actually have scattering data delivered by Pabst and Heftberger for DOPC/cholesterol system. So the data from this mixture is under interest. However, you mentioned that there is NaCl in the simulations. Since NMRlipids II shows that Na binding is too strong in Slipids, the data with ions cannot be used for the purposes in this project. However, the DOPC/cholesterol data without ions would be under interest.
Delete>>> Samuli wrote: ... but I think that performing NAMD vs. Gromacs comparison also by ourself is probably reasonable.
ReplyDeletethe mentioned paper (http://pubs.acs.org/doi/abs/10.1021/acs.jctc.5b00935) also suggests openMM as a good alternative sim package that does not suffer the same problem as Gromacs (if it is a problem on Gromacs side). Since openMM can be fed with Gromacs files (and Charmm files too) directly, this might be a way to check this issue. One can even try both sets of input files - or use the optimized inputs from Charmm-gui.
In addition - in the original CHARMM36 paper (http://dx.doi.org/10.1021/jp101759q), they use CHARMM package; NAMD is there probably used only for validation of two simulations as they share the same settings (but DOPC is 4x larger). So we should probably target on compatibility-test Gromacs-Charmm rather then on Gromacs-Namd.
DeleteAs far as I can remember, NAMD is used in quite many applications of CHARMM36. So I think that the consistency of the results also form NAMD with other programs is under interest. Let's continue this discussion in GitHub issue: https://github.com/NMRLipids/NmrLipidsCholXray/issues/4
DeleteJust a small remark : Comparing exp/simulation, I think one should
ReplyDeletekeep in mind that there may be differences in order
parameters due to the imprecision of the area/lipid = f(tension) curve
and also due to the structure at a given area per lipid.
To focus on the structure, maybe one should try to map the exp. data
on simulations with various area/lipid, or various surface tension.
Imposing a zero surface tension at all price, one may eliminate a model which is maybe relatively good for the structure, but not so for the
lateral compressibility.
Of course the best is a model good at both...
Here a paper by S. Feller etal where he detailed the impact of area per lipid on order parameters for pure DPPC, using CHARM22b4b ff, and compared to the NMR data.
Deletehttp://pubs.acs.org/doi/abs/10.1021/la970746j
Scott E. Feller et al., Langmuir 1997, 13, 6555-6561 "Computer Simulation of a DPPC Phospholipid Bilayer: Structural Changes as a Function of Molecular Surface Area"
I thought of something similar, when I said that one should maybe try bad NMR/simulation matching at various area per lipid.
I see the point. Similar thing has been done in order to interpret the scattering data, as written in the NMRlipids V paper (http://dx.doi.org/10.1016/j.bbamem.2016.01.019, page 2524, left column) "Several models, reviewed by Heberle et al. [182], are developed to give structural interpretation for the form factor data [183], while also MD simulations are used [169,184-187]. In these studies the area per molecule is often fixed to a valueminimizing the differences between experimental and simulated form factors [169,184-187]. Depending on the model used, this area per molecule may be close to [187] or deviate significantly [184,185,169,186] from the value predicted by the model in
Deleteconstant pressure simulations. However, with optimized area per molecule all models give form factors close to the experiments, despite of the bilayer tension generated in some models."
Here we are looking a model which would give simultaneously good order parameters and form factors. Fixing things by tension would be problematic for models which gives good order parameters, but not good form factor (like Berger/Höltje combination above) or vice versa. By fixing the other one would make the other worse. On the other hand, one goal of this work is to see if the lipid-cholesterol interactions can be correctly described with current models. I think that this would require the same tension with all cholesterol concentrations.
In conclusion, I think that it would be extremely difficult to fix simultaneously the order parameters and form factors with different cholesterol concentrations by inducing tension.
I think we should continue the discussion about CHARMM-NAMD-GROMACS results in the GitHub issue: https://github.com/NMRLipids/NmrLipidsCholXray/issues/4
ReplyDeleteThere is now some new content also from Claire Loison.
Hi, here is POPC with 0,10,20,30,40,50% cholesterol with the Charmm36 force field using Gromacs 5 and the proper simulation parameters (provided by Charmm-GUI and quoted in http://www.gromacs.org/Documentation/Terminology/Force_Fields/CHARMM )
ReplyDeleteThese runs are 100 ns at the moment, but I'll extend them if needed.
http://dx.doi.org/10.5281/zenodo.61649
http://dx.doi.org/10.5281/zenodo.61992
Just an extra note if someone wonders this during some analysis: The systems with 0–30% chol have a total of 9000 water molecules (45 per POPC molecule). However, for the systems with 40% and 50% chol this would not be sufficient, so they have 45 water molecules per lipid (either POPC or cholesterol).
DeleteThe re-simulated systems with 0–30% cholesterol are here:
Deletehttps://doi.org/10.5281/zenodo.159759,
see the discussion here:
https://github.com/NMRLipids/NmrLipidsCholXray/issues/9