Data contributions

[update 6.7.2015] This page is outdated, however, it is preserved here due to the scientific information in the comments. We have currently four subprojects (see About) and the data contributions would be done ideally directly uploading the data into the appropriate GitHub repository (see Workflow) and/or commenting the related discussion.

The project has been success by this far. There are currently 16 posts and a total of 141 comments containing very intense information. Consequently, the blog has become quite difficult to follow.

One problem is that each comment is gets filed only under one post even though it would be related to  a large amount of topics. We hope that in the future there would be a system, where comments could be tagged to a certain topic and then would appear in all the discussions related to that topic (we have never used Twitter but have understood that it roughly works with this principle).

However, since such a system is not available here, we think that from now on it would be good to start sending all the new data by commenting this page. This will make it a lot easier to find and view the data.

So, let us keep the other comments under each separate post, but from now on, post new data here.

[update 12.11.2014] Since we are now using GitHub (https://github.com/NMRLipids/nmrlipids.blogspot.fi) you can also put your files there by doing pull request. Alternatively, if you want to include your files to the GitHub without learning to use Git, you can ask me (Samuli) to add the files. I am planning to share also my trajectories though the Zenodo system. Currently, the most up to date GitHub folder is for the GAFFlipid: https://github.com/NMRLipids/nmrlipids.blogspot.fi/tree/master/POPCgaff. It probably contains most relevant files, including the link to the trajectory in Zenodo server. If you add trajectories or other content to Zenodo, you can add it to the Nmrlipids community: https://zenodo.org/collection/user-nmrlipids.

49 comments:

  1. Hi everyone

    I have recently measured the order parameters of the polar headgroup region and the glycerol backbone of DPPC from 100 ns atomistic MD simulations of lipid bilayers performed with the Slipids forcefield at different temperatures. This analysis was performed in collaboration with Markus Miettinen and Samuli Ollila as part of the joint research project developed by Vasily Oganesyan at University of East Anglia (Norwich) and Mark Wilson at Durham University.

    Here is a list of order parameters with the same order reported in the nmrlipids project:

    283K
    -0.04839095
    0.02979915
    -0.0534666
    -0.0478858
    -0.063991
    -0.122031
    -0.165568

    298K
    -0.0460633
    0.02968925
    -0.0580725
    -0.0442726
    -0.0684728
    -0.121197
    -0.161122

    303K
    -0.0426849
    0.03029045
    -0.0595338
    -0.0419392
    -0.0697026
    -0.127311
    -0.161753

    308K
    -0.0432587
    0.0341308
    -0.0517974
    -0.0321559
    -0.0642035
    -0.116106
    -0.154154

    312K
    -0.04572355
    0.0361404
    -0.0486379
    -0.0447603
    -0.0641615
    -0.129244
    -0.164287

    313K
    -0.044001075
    0.039223
    -0.0467185
    -0.03301635
    -0.05501665
    -0.1117335
    -0.161073

    314K
    -0.0455145
    0.0414877
    -0.0531463
    -0.0432409
    -0.0595313
    -0.12467
    -0.163013

    318K
    -0.04094155
    0.0393655
    -0.0431368
    -0.0255703
    -0.053094
    -0.10797
    -0.155108

    323K
    -0.0401585
    0.0383469
    -0.04914
    -0.0234394
    -0.0596508
    -0.111834
    -0.149935

    333K
    -0.037018
    0.04167645
    -0.0389662
    -0.035944
    -0.0519808
    -0.111659
    -0.148267

    Please find these data files, a spreadsheet with plots of the results and a short presentation of them in the following link

    https://drive.google.com/folderview?id=0BwkZuUeVCF5tQ1dvZmIzamlPSHc&usp=sharing

    Thank you very much in advance for your future feedback,

    Andrea Catte

    ReplyDelete
    Replies
    1. Thank you for the contribution!

      It is very interesting to have results with such a wide temperature range, especially because the corresponding experimental data is available: http://dx.doi.org/10.1021/bi00687a021

      I have a couple of questions:
      1) There is only one order parameter value reported here for the alpha and beta carbons. However, in previous data for Slipids there is significant forking observed for these segments:

      DPPC
      https://github.com/NMRLipids/nmrlipids.blogspot.fi/blob/master/DATAreportediINblog/DPPC/SLIPIDS-323K_blogged-30-09-13.dat

      POPC:
      https://github.com/NMRLipids/nmrlipids.blogspot.fi/blob/master/DATAreportediINblog/POPC/NaCl/SLIPIDS-298K-150mM_blogged-13-09-13.dat

      Are they overlapping in your case or have you reported the average? It would be useful to have both values.

      2) Do you see the expected phase transitions below ~313K for DPPC?

      3) Would you happen to have a POPC simulation with Slipids without ions? We are missing this data and it would be useful.

      4) Have you checked how strong is the ion Na partition into a bilayer in your simulations?

      Delete
  2. Thank you very much for your comments and for your questions! you're welcome for the contribution!

    I apologise to you for the lateness of my response, but I have seen your comment to my post only now (long busy days..).

    I have reported the average, they aren't overlapping, I can share both values (alpha carbons first) with you and everyone else who might be interested in these order parameters. Please find them below:

    283K
    0.0253274
    0.0342709
    -0.0252238
    -0.0715581

    298K
    0.0330744
    0.0263041
    -0.0282685
    -0.0638581

    303K
    0.0314412
    0.0291397
    -0.0251287
    -0.0602411

    308K
    0.0437695
    0.0244921
    -0.0288747
    -0.0576427

    312K
    0.0378024
    0.0344784
    -0.0293169
    -0.0621302

    313K
    0.0434599
    0.0349861
    -0.0297872
    -0.05821495

    314K
    0.045582
    0.0373934
    -0.028806
    -0.062223

    318K
    0.0439091
    0.0348219
    -0.025613
    -0.0562701

    323K
    0.0453405
    0.0313533
    -0.0247313
    -0.0555857

    333K
    0.0448921
    0.0384608
    -0.0213869
    -0.0526491

    you will also be able to find these values in a spreadsheet file you can access using the same link of the previous post.

    Yes, I see the expected phase transition from liquid crystalline to gel below ~313K. I have also some images of the final structures of the simulations, which might help to visualize the phase behaviour of the different systems, and I might be able to share these files as soon as possible.

    No, I don't have any POPC simulation with Slipids without ions. Anyway, I can try to allocate some computational time to perform this set of simulations as well.

    I haven't checked how strong is the Na partition into a bilayer in my simulations. I can check it and I can share these results with you as soon as possible.

    Thank you very much again for your comments, for your questions and for your feedback!

    ReplyDelete
    Replies
    1. I did plot the results together with the experimental results by Gally et al. for the alpha and beta carbons numerically available here:
      https://github.com/NMRLipids/nmrlipids.blogspot.fi/blob/master/DATAreportediINblog/DPPC/temperature/EXP-Gally1975.dat

      The comparison is here:
      https://www.dropbox.com/s/x9bg6slcebuellt/OrderParameterTdependency.jpg?dl=0

      In experiments there is a small discontinuity in the transition from the ripple phase to the liquid phase close to 313K. There seems to be room for improvement in the temperature dependence of the model.

      I did also plot the results together with the previously reported Slipid DPPC results without ions:
      https://github.com/NMRLipids/nmrlipids.blogspot.fi/blob/master/DATAreportediINblog/DPPC/SLIPIDS-323K_blogged-30-09-13.dat
      and the result is here:
      https://www.dropbox.com/s/ai3au9w50vjdycx/OrderParameterIONSnaclSLIPIDS.jpg?dl=0

      The change in the alpha order parameter due to ions seems to be overestimated, however not as much as in MacRog and Berger. This would indicate that the ion partition would be overestimated but not as strongly as in MacRog and Berger (see:
      http://nmrlipids.blogspot.fi/2014/05/response-of-headgroup-and-glycerol.html)
      However, there is no change in beta and the previously reported DPPC results were from quite short simulation. It would be nice to see the Na density profile together with the lipid density profile from these simulations.

      Delete
    2. I apologise to you for the delayed response, but I've been dealing with different things lately... I thank you very much for sharing the interesting plots of the data compared to the experimental results.

      I agree with you that there is still room for the improvement in the temperature dependence of the model. Although it is quite remarkable how these simulated results almost match the experimental one at some temperatures interestingly below and above 313K for beta and alpha order parameters, respectively.

      Now I am writing you what I have recently done in terms of the calculation of the partial density profiles of Na+ ions, lipids and three different lipid moieties. I have basically normalized all the density profiles (otherwise it would have been difficult to see where the ions were located). Please find these results in an excel file you can access using the follow link:

      https://drive.google.com/file/d/0BwkZuUeVCF5td0JlM3ZLRWJfSTQ/view?usp=sharing

      Please let me know if there is a better way to represent the Na and lipid density profiles.

      Thank you very much again for your feedback and for the future help,

      Andrea

      Delete
    3. Hi Andrea,

      I have two small clarifying questions:

      1) From the file dppcphg_order.xlsx I get the impression that OP values reported for 213K are actually an average of (at least) two (213c and 213d) independent simulations? If you confirm this is the case, then we will have two measurement points at 213K, which is good, as it demonstrates the accuracy of the calculation.

      2) I found that the simulations were at least 100 ns long. How much of the total time was used for analysis and how much for discarded for relaxation purposes?

      Thanks in advance for the information!

      Delete
  3. The new united atom forcefield by Tjörnhammar and Edholm was just published recently ( http://pubs.acs.org/doi/abs/10.1021/ct500589z ) and the supplementary information has parameters for a 144-lipid DPPC system in gel phase without ions. I ran a 5ns simulation at 70 degrees in order to destroy the ordered gel configuration followed by 200ns simulation at 50 degrees. I analyzed the last 100ns of this simulation. And I used the same mdp file as in the paper's SI except for the temperatures.

    The results are:

    0.0495 0.0455
    0.0203 0.0195
    S: -0.2322 R: -0.2909
    -0.3211
    S: 0.0807 R: -0.2275

    ReplyDelete
  4. Thanks! The results are significantly different than with Berger even though they say the paper that "The bonded and nonbonded parameters not explicitly given here were taken as in the Berger force field." and I cannot find any discussion about glycerol backbone or headgroup dihedrals. On the other hand, with a quick look of the parameters it seems to me that the dihedrals for, at least, glycerol backbone are not the same as in Berger. For example, for the g3-g2-g1-O(sn-1) dihedral I find the lines

    FROM Berger:
    12 13 32 33 1 0.0 5.85 3
    12 13 32 33 1 0.0 0.42 2

    FROM Tjörnhammar and Edholm:

    12 13 32 33 1 0.000 5.92 3

    So for me it seems that there are different dihedrals between the models (correct if I am wrong). This difference would explain the different order parameters. I have no idea about the origin of neither of these potentials.

    It would be useful if Matti would put some equlibrated trajectory to the Zenodo (there is also the recently created nmrlipids community: https://zenodo.org/collection/user-nmrlipids). I would be happy to calculate some dihedral distributions.

    ReplyDelete
    Replies
    1. I included Matti's results for the model Tjörnhammar and Edholm in the plot together with some other results:
      https://www.dropbox.com/s/d12c3bu8i90nht9/HGorderparametersTjornEd.jpg?dl=0

      I did also calculate the dihedral angle for g3-g2-g1-O(sn-1) dihedral from the trajectory in Zenodo uploaded by Matti:
      https://www.dropbox.com/s/i4uo3oq3r6xbowm/g1-g2CdihTjornEd.jpg?dl=0

      Compared to the dihedral angle distributions from Berger it looks quite different. Instead it looks more similar to the models with more realistic order parameters for g1 segment:
      https://www.dropbox.com/s/0lwl4m6x7mquy8l/g1-g2_Cdihs.jpg?dl=0

      Indeed, the Tjörnhammar and Edholm model has better order parameter for g1 segment than Berger.

      Anyway, for me it seems that the dihedral parameters are not the same in the Tjörnhammar and Edholm model as they were in Berger.

      Delete
  5. Dear Samuli,
    thanks for your comments and for noticing the problems.
    We have checked the things you have noted.

    1) The coefficient in front 3 fold symmetric term is 5.8576 kJ/mol in our original file. For testing purposes, we also used the value 5.92 (taken from the Chui force field) , but this was never removed in the itp-file. This will, however, only have a very small effect upon the barrier height as seen from the enclosed eps-file. I think there is no deep reason for the small difference in the values. This looks like something that may come from repeated conversions between kcal and kJ and rounding off with different numbers of decimals.

    2) I am much more worried about the small twofold symmetric term. This alters the balance between the gauche and trans states.
    I looked through our itp-files going as far back as I have them, which was to about year 2000 and could not find any twofold symmetric term there. The itp-files on Peter Tielemans site do contain both the two fold symmetric and the three-fold symmetric ones. I do not really understand the purpose with the twofold symmetric term. However, in whatever way one constructs the dihedral potential one may adjust it by choosing proper 1-4 LJ-interactions.


    Best regards Olle

    ReplyDelete
    Replies
    1. Thanks for the information! I started to dig some details about parameters from Peter Tieleman's webpage and noticed that the references were outdated. I asked about that from Peter Tieleman directly and he replied with emails containing very useful information. With his permission I copypaste the content of the two emails into the separate comments below.

      Delete
    2. From: Peter Tieleman
      Sent: Monday, December 01, 2014 5:08 PM
      To: Ollila Samuli
      Subject: Re: References for the itp files distributed in your webpage

      Hi Samuli,

      I converted the GROMOS topology from Berger/Edholm to GROMACS and tried to match the parameters from that GROMOS topology to the parameters in Berger's paper, but 'similar' is an apt description. Celine Anezo did some more digging and several parameters remained unclear. The O-C-C-O dihedrals and on phosphates in the original GROMOS had/have two dihedrals on them, with different multiplicities, presumably to account for some carbohydrate effects, but this traced back very far in GROMOS and I don't know where this came from originally.

      All the itp files that need lipid.itp are based on Berger, but of course Berger's paper only looked at DPPC. The DPPC comes from a paper with Marrink, Biophys. J. 1998, on pulling lipids. The POPC.itp come from Tieleman and Berendsen/Sansom on alamethicin, Biophys. J. 1999. dmpc.itp should be the same as dppc.itp and was used in a paper on gramicidin A with Bert de Groot. dopc.itp was used in a number of simulations by Bennett, and both DOPC and POPC itp had different versions because the dihedrals near the double bonds were wrong originally. plpc.itp is from Bachar et al, 2004, JPCB. Some of the missing parts have never filled in for distribution (g96_lipid.itp) because the whole approach with lipid.itp is flawed, straight combination rules work better and I never really trusted G96 much.

      Cheers,
      Peter

      Delete
    3. From: Peter Tieleman
      Sent: Wednesday, December 10, 2014 2:02 AM
      To: Ollila Samuli
      Subject: Re: References for the itp files distributed in your webpage
      Hi Samuli,

      You're welcome to copy my comments to the blog.

      To the best of my knowledge the parameters I used in our 1998 paper, which has Berger on it as author, are a correct conversion from GROMOS files to GROMACS. If I remember correctly, the GROMOS files came from Siewert-Jan Marrink, who was at the time in Fritz Jahnig's group. The "Berger force field" is a misnomer in some ways, as Berger only parameterized LJ parameters for pentadecane (and if I recall correctly used the wrong experimental values to parameterize against). The other parameters were taken from existing force fields, including non-bonded from united atom AMBER and OPLS, bonded from GROMOS, and charges for the head groups from Chiu et al, which themselves were based on rather simple QM calculations and rather major rounding. Celine Anezo looked carefully at the origins of the parameters in Berger's paper and had difficulty tracing all of them. There is no reason to assume that this combination of parameters gives the right dihedral distributions or head group dynamics although the tails did quite well based on Ryckaert-Bellemans dihedrals (which themselves were based on CH3-CH2-CH2-CH3 in butane). As Olle points out in this blog, 1-4 interactions and dihedral parameters combined determine dihedral properties, and the combination hasn't been carefully looked at in the 90s. Nor would we have had enough computer time and experimental data to do this well, I think. The combination of decent tail dihedrals and sufficiently polar head group/glycerol backbone charges seems to be enough to get reasonable PC lipids, which afterwards nearly all force fields managed. It seems obvious to me there is room for improvement in parameters.

      Berger's paper has been quite influential because of the large number of people who used parameters similar to his, but it had its problems. It was also a very early paper in some sense, given what we now know about correlation times in lipids and in collective properties. Technically, there also are differences between how GROMOS and GROMACS treat or treated things like neighborlists and charge groups, some of which we explored in papers. I don't see much point in trying to trace the 'original' Berger force field. I've tried to make topologies available for papers we used them in, and provided some updates in errors we found (e.g. bond lengths for N-CH3 in PC lipids, dihedrals are bonds adjacent to double bonds in tails). As a community it is useful to make topologies available for specific papers, and I have no doubt that it is useful to reconsider all lipid parameters from time to time.

      Cheers,
      Peter

      Delete
    4. Thanks to Peter for sharing invaluable information!

      Indeed "Berger" and most other force fields have been succesfull to describe a reasonable PC lipid bilayer, despite of the incorrect atomistic resolution structure of the glycerol backbone and headgroup. The success of coarse grained models might be related also to this issue; the "general" bilayer properties seems to largely dominated by the acyl chain properties. However, I think that there are situations where atomistic resolution structure of glycerol backbone and headgroup are potentially relevant.

      From the practical point of view we currently have the following problem: In my understanding the force field distributed in Peter Tieleman's webpage is commoly called the "Berger force field" and the original paper by Berger is cited in this context (this is done also in the current version of our manuscript). However, after this discussion and the results above it seems that, at least, the glycerol backbone dihedrals (thus the order parameters and structure) are not the same in the orginal Berger model and in the one distributed in the webpage. Consequently, we cannot call the parameters from the webpage as "Berger" and cite the original work. We have to call this something else, however, in the way that people understand easily what we are talking about (Calling it Marrink et al. force field would be confusing from several point of views).

      Any suggestions for this?

      Delete
    5. From the above discussion it is clear that referring to this family of force fields simply as "Berger" is indeed misleading.

      However, as "Berger" has been used for more than a decade now, one can no longer simply rename these offsprings as Marrink, Tieleman, Anezo, … although this probably should have been done so back in the day to avoid this confusion.

      One suggestion (from the current version of the manuscript, https://www.dropbox.com/s/us0d805lpg53kqe/HGmodel_draft2.pdf?dl=0) has been to use Berger-Tieleman or Berger-T for this family of force fields that is currently understood as the "Berger". While this is an improvement, I find that it might lead to problems when citing. Namely, when using different Berger-T lipids one would need to cite different original publications. Would people really remember to cite Marrink when using the DPPC of BergerT and Tieleman when using the POPC of BergerT?

      I would suggest adopting a naming scheme based on the publication years of the original force fields. That is, the very first Berger would be called Berger97 [citing Berger et al.], the saturated PCs Berger98 [citing Marrink et al.], the unsaturated PCs Berger99 [citing Tieleman et al. but noting that presently the double bond has been fixed], ...

      Delete
    6. Using the year might be good as well. However, one issue is that there are also other modifications to Berger done by other people than Tieleman et al. and using only the year might lead to ambiguous names. For example, Marja Hyvönen and myself have implemented polyunsaturated lipids based on DPPC file shared by TIeleman and work by Bachar et al. [http://dx.doi.org/10.1021/jp065424f (2007)]. Now the question is that should we call this as Berger07, or should we say that this is modified Berger98? In the latter case we restrict the usage of years to the files shared in Tielemans webpage. And in this case it might be clearer to just call them Berger-Tieleman or Berger-T. In the first case, there is potential problem for ambiguity since it is very well possible that someone else has published some other Berger modification in 2007.

      Delete
    7. The non-uniqueness problem Samuli mentions is of course a valid point.

      After chatting about this with him on skype we more or less converged to suggest a scheme "Berger-moleculename-year" to be used in the nmrlipids publications. This will make clear that the FF belongs to the family of FFs commonly referred to as "Berger FF", but also will make it clear that the parameters are not from the original Berger publication. Furthermore, this scheme should to a large extent avoid the problem on non-uniqueness that Samuli pointed out above.

      Using this scheme the original Berger force field would be called Berger-DPPC-97 [citing Berger 1997], the one with the glycerol parameters that "Berger" FFs actually use Berger-DPPC-98 [citing Marrink 1998], and the unsaturated ones for example Berger-POPC-99 [citing Tieleman 1999, but noting that the double bond has been fixed since]. And following this scheme each of the FFs in Samuli's and Marja's 2007 work will have a name of its own: Berger-PAPC-07 [Ollila 2007], Berger-PDPC-07 [Ollila 2007].

      Delete
    8. This might be good also because, as far as I understand, all these "Berger" models are using the infamous lipid.itp while people have done various models for different molecules, like dppc.itp, popc.itp, etc. Then the Berger in the name would mean that lipid.itp has been used, the molecule would refer to the *itp file for the specific molecule, and the citation for the first publication of the molecular *itp.

      However, now the Berger POPC used in the manuscript would be Berger-POPC-07 [Ollila 2007] since the file I use is the corrected Berger-POPC-99 [Tieleman 1999], first time introduced in Ollila 2007 paper. I am not sure if it is fair to use it like this Berger-POPC-07 [Ollila 2007] since we changed like two lines directly based on other peoples work. We can add more citations Berger-POPC-07 [Ollila 2007, Tieleman 1999], but it gets more complicated.

      Delete
    9. In this particular case I would suggest introducing the Berger POPC used in the manuscript as "the fixed [Ollila 2007] Berger-POPC-99 [Tieleman 1999] from Peter Tieleman's web page [link and date]", and use just "Berger-POPC-99" in the rest of the manuscript. Then those people who know about the correction-issue will be served. And those who don't, will not get unnecessarily confused.

      Delete
    10. The issue here might that if this notation is used in future, and someone is actually using the Berger-POPC-99 without correction, it has the same abbreviation which we have used for the corrected version. This might be confusing. For this reason, I have left it as it is for now.

      Delete
  6. Hi

    We are an experimental biophysics group in Graz (our website: http://www.uni-graz.at/~pabstg/) and we are doing small angle x-ray scattering (SAXS) on multilamellar lipid vesicles. With our analysis we calculate x-ray formfactors (FF) and determine structural parameters of bilayers such as bilayer thickness or area per lipid etc. We also compare our experimental FFs to FFs obtained by MD-simulations.
    In our last experiments we measured x-ray FF of two different lipid mixtures containing DOPC/DPPC/CHOL and DOPC/DSPC/CHOL and compared them to FFs determined by MD-simulation (CHARMM 36). The exact compositions and some details about the simulation are reported in the following paper: http://pubs.acs.org/doi/abs/10.1021/ct400492e. Information about our x-ray analysis you can find here: http://journals.iucr.org/j/issues/2014/01/00/fs5056/fs5056bdy.html. Now the results. The simulated FFs of the liquid disordered (Ld) phase (mainly DOPC and CHOL) of both systems (DOPC/DPPC/CHOL and DOPC/DSPC/CHOL) don't match the experimental FFs and the minima of the simulated FFs are shifted towards lower q-values (scattering vector) meaning that the simulated bilayer is thicker than the experimentally measured bilayer. Whereas the FFs of the liquid ordered (Lo) phase (mainly DPPC or DSPC and CHOL) are in better agreement. They don't match perfectly but it is better and the bilyer thicknesses are also comparable. A short summary with some figures is available here: https://drive.google.com/file/d/0BypC9VEx1ib_VElHNS1HaFpOYWc/view?usp=sharing. Data files of the experimental FFs are here: https://drive.google.com/file/d/0BypC9VEx1ib_NEJOYTVlbXJRa2c/view?usp=sharing , and data files for the simulated FFs can be found here: https://drive.google.com/file/d/0BypC9VEx1ib_b1hHMnVLZ0JSZGs/view?usp=sharing.

    It is interesting that especially in the Ld phase, where you have mainly DOPC/CHOL and therefore double bonds in the hydrocarbon chain, the bilayer structure in the simulation seems to be different from the experimentally measured structure. We also measured some other lipid systems such as DOPC, POPC, DPPC, DOPC/CHOL, POPC/CHOL and DPPC/CHOL. As these are very common lipids, which are also used in simulations, I can also provide FFs for them and perhaps we can compare them with simulated FFs. A brief description and figures you can find here: https://drive.google.com/file/d/0BypC9VEx1ib_NlBtazBvbEpmUWM/view?usp=sharing , and datafiles for the experimental FFs are here: https://drive.google.com/file/d/0BypC9VEx1ib_dmktOHYyRFJFU2s/view?usp=sharing.

    I am looking forward to your feedback.

    peter

    ReplyDelete
    Replies
    1. Thanks for your contribution! This is nicely expanding the scope of this project.

      First thing which comes to my mind is that you report the Form Factor for POPC/Chol mixture (80/20 mol %). The simulation of the corresponding system with CHARMM36 model is reported in this blog by Hubert Santuz (the reported Berger and MacRog simulations do not have this molar ratio). As far as I understand, you have comparison between simulations only for more comlex system (with some problems) while in the literature single component systems with this force fields are in good agreement with scattering data. I think that it would be very interesting to see if the model can reproduce the POPC/chol result.

      Is the electron density profile calculated from simulation enough to calculate the Form Factor or do you need some other information? If the electron density profile is enough, then it would be relatively little work to share that here so that someone could calculate the Form Factor.

      Delete
    2. Peter also shared SAXS data for pure POPC. In the original C36 paper by Klauda they only show the FF for DPPC. Samuli, however, has also simulated pure POPC with C36. Thus one could (also) calculate the FF from Samuli's data set and see how it matches the experiment; if the match is fine, then the problem Peter reports would come from the mixing with cholesterol, which would be interesting.

      Delete
    3. Indeed, for single component systems we have more data from different force fields also for other lipids. I think that we should figure out which would the most convenient way to calculate the Form Factors. I have never done so all suggestions are welcomed.

      Delete
    4. The formfactors can be calculated with a program called "SIMtoEXP" by Norbert Kucerka. Our colleagues at Weill Cornell use this program and my friend Milka Doktorova sent me some information about it: In order to calculate the FF, you need to have a file with the number densities of all atoms in the system, then the software converts that to a FF and if you also upload your SAXS/SANS FF data, it can directly compare them. Here is a link to this version of the program that she and I use: http://www.norbbi.com/links/links.html#Computer (follow the links: Software-> SIMtoEXP and there you will find the installation files, a manual and some test files). Here you can find also a paper about the software: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2876336/.

      There is also a newer version of the program available by Bryan Holland. She wrote that that in this version you just upload your trajectory (the dcd file), and it computes the number density and form factor for you, but she has not tried this version yet. Here is a link to it: http://sourceforge.net/projects/simtoexp/.

      I would say let's try the newer version and hopefully it is really so easy that you just have to upload the trajectories....

      Delete
    5. PART 1: I made a quick and dirty script which downloads the simulation data from Zenodo and calculates the form factors for all the simulations for which the trajectory related to this project has been shared (https://zenodo.org/collection/user-nmrlipids). The script (FFcomp.sh, FFstructCALC.sh) and output (electron densities and form factors) can be found from https://github.com/NMRLipids/nmrlipids.blogspot.fi/tree/master/FFactor. Notice that the script downloads all the required files automatically from the Internet. You need some diskspace for the trajectories if you try it. One motivation to do this was to demonstrate how useful it would be to share the trajectories in Zenodo to allow this kind of analysis.

      I plotted some results:

      Form factors from GAFFlipid, LIPID14 and Berger POPC compared to experiment (shared by Peter Heftberger):
      https://www.dropbox.com/s/ihil4mmy9ezpdyx/POPCcomp.jpg?dl=0
      According to original publications GAFFlipid and LIPID14 form factors are in good agreement with experiments. I cannot reproduce this. When I look at the experimental form factor in GAFFlipid paper, it looks different than in LIPID14 paper which is weird.

      Form factors from SLIPID DPPC compared to experiment (shared by Peter Heftberger):
      https://www.dropbox.com/s/ptlm5tisy5xhg9g/DPPCcomp.jpg?dl=0
      This looks a bit different compared to the original Slipid paper which was in good agreement with experiments. However, I noticed that there are different vdw cut-off settings in the mdp file that they share in the Internet compared to the original work. We also seem to have lower area per molecule compared to the original work. See explanation in mdp file in Jukka's zenodo fileset (https://zenodo.org/record/13287/files/mdnew.mdp).

      Form factors from POPC/chol simulations with Berger/Höltje model:
      https://www.dropbox.com/s/rq263o2scstia9u/POPCcholcompBERGER.jpg?dl=0
      There is systematic change as a function of cholesterol which looks roughly similar to the experimental data. However, the highest cholesterol concentration looks weird.

      Form factors from POPC/chol simulations with Berger/Höltje model compared to the experiments (shared by Peter Heftberger):
      https://www.dropbox.com/s/pba0hojc6iww0jf/POPCcholcompEXP.jpg?dl=0
      Qualitatively the change with cholesterol is similar as in experiments (which is well known), however, quantitative comparison is difficult with this data since there is not excatly the same concentration in simulations and experiments, and the pure POPC is not in quantitative agreement.

      Delete
    6. PART 2: This quick and dirty study raised following thoughts about this topic. In literature very good agreement for form factor between simulations and experiments for pure bilayer systems are nowdays reported. For this reason I thought that all the models reproduce bilayer density profile well. However, it seems that it is not trivial to reproduce some of these results. Anyway, what I have not seen carefully discussed literature is the quantitative accuracy of the lipid cholesterol interactions in simulations (correct if I am wrong). It is well known that the cholesterol orders the bilayers, however, there are all kind of suggestions about molecular level interactions, phase separation, microdomains etc. If these issues are studied with simulations one needs get have the lipid-lipid, lipid-cholesterol and cholesterol-cholesterol interactions correct. I think that the POPC/chol system would be very nice system to try to quantify the accuracy of lipid-cholesterol and cholesterol-cholesterol interactions due to the available experimental data. We know all the order parameters (also for cholesterol) up to 60 mol% [http://dx.doi.org/10.1039/C2CP42738A]. Peter Hefberger also shared structure factor data up to 20 mol%. Does anyone know if there would be form factor (or other scattering data) with higher concetrations? I think that the model which reproduces all the order parameters and scattering data from low to high cholesterol concterations would most likely have the correct (effective) intermolecular interactions.

      I think that it would be interesting to see the structure factors for POPC/cholesterol system simulated with other models as well. The CHARMM36 and MacRog results for the headgroup and glycerol are already reported in this blog, but since I did not have the access to trajectories I could not calculate those.

      I will soon write a blog post describing more detailed about my thoughts on this.

      Delete
    7. This is interesting. I would have expected much a better agreement. Peter [Heftberger] today also mentioned an article where simulations are compared to experiments for DOPC (http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3673311/). I did not really read through this, but there the area per lipid has been tuned in order to improve the fit.

      This is not really satisfying if you do not want to do this. Concerning form factors of POPC containing cholesterol. I would check for publications by Kucerka, Katsaras, and Nagle, but I am not aware of any. We have a beamtime coming up in about 3 months. Possibly we could run some of these samples if they are of interest. In any case 60 mol% cholesterol is beyond the solubility limit of cholesterol in POPC bilayers. We will be good, if we manage 50 mol%.

      Delete
    8. POPC:
      I also had a look at the two papers about GAFFlipid and Lipid14. The experimental formfactors (FF) are taken mainly from Kucerka 2011 (http://www.sciencedirect.com/science/article/pii/S0005273611002276. FF are plotted in the supportings). I also compared the plotted exp FFs of the two papers, but they look the same in both papers, especially the minima are at the same q-values. But you are right Samuli that there is a minor difference in the amplitude of the FFs, namely the amplitude of the FFs in the GAFF paper is slightly lower than in the Lipid14 and Kucerka’s original paper. Perhaps they downscaled the exp FFs so that they match better with the simulated ones, but more important are the minima and they are the same in both papers. Our exp FF of POPC (25°C) also compares well with the reported one by Kucerka, whereby Kucerka measured POPC at 30°C and not 25°C. I am also a bit surprised that you couldn’t reproduce the FF of GAFF and LIPID14. Is this a problem of the model to reproduce the right electron density profile and calculate the FFs out of the simulation as you mentioned? Could it be that they applied certain constraints (like area/lipid) to reproduce the “right” FF, which match with the exp ones? Because in the article that also Georg mentioned (http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3673311/) they tuned parameters of their simulation, so that the simulated FF matches with the exp FF for DOPC.

      DPPC:
      The match of simulated and exp FF is very good in the original paper. Also our SAXS FF matches quite good with the them, whereby the second minimum in the FF at q=0.425 is shifted towards a little bit lower q-value in our FF. This could be due to the lower resolution at higher q-range in the scattering profile we have for this DPPC sample. But we will try to measure also DPPC again and DPPC/CHOL mixtures within the next Beamtime.

      POPC/CHOL with Berger/Höltje:
      First of all it encouraging that Berger/Höltje can reproduce this relative change in the FFs when cholesterol is added. But a major problem is that with Berger the simulations show in general a thinner bilayer (minima in FF are shifted towards higher q-values). It seems that this is a kind of “systematic problem”.?

      As Georg mentioned in his post, I also think that it is a good idea that we focus on this POPC/CHOL system, which we don’t have to heat up and can measure at room temperature (30°C). We only have to decide which samples should we try to measure and at which temperature? I would suggest following samples:
      *POPC, POPC/CHOL (5%, 10%, 15%,….50% cholesterol) at 30°C (other suggestions?), because of comparison to literature data. I also think that it doesn’t make sense to go to higher chol concentrations, because then chol is not soluble anymore and perhaps this is also a reason for the weird FFs in the simulation?
      *DPPC, DPPC/CHOL (10%, 20%... 50%) at 50°C.
      Do you have any other ideas which lipids we should measure that can be of interest for simulations?

      Delete
    9. Comparison between experiments and simulations for one component bilayer: I was also expecting better agreement since nowadays almost all parametrization papers report this, and the agreement is typically very good. We have to look at this more carefully. There are now more trajectories in Zenodo. I will analyze these soon.

      Reproducibility: We have to look this more carefully. There are quite often these kind of situations that one gets different result than in papers with the first try. However, usually these things gets somehow resolved. The area per molecule was not fixed in the original papers referred here. For Slipids the explanation is most likely the above mentioned cut-off issue.

      About intensities: If I have understood correctly, also the amplitudes give information about the density profile? For example, one might have the correct thickness but incorrect density distribution inside. In this case minima would be correct but amplitudes wrong?

      Berger/höltje cholesterol: The POPC model used here gives larger area per molecule than many other models which (according to the literature) agrees with experimental form factor. So it is probably too thin, indeed. With cholesterol it gets thicker (area per molecule decreases, i.e. condensing effect) as does all the other models as well. The question I would like to anwer is that is the quantitative amount of thickening (area per molecule decrease) in agreement with experiments in different models? There are already simulation trajectories with cholesterol from MacRog model available in Zenodo and with CHARMM36 hopefully coming. I will analyze these soon.

      Cholesterol solubility: Order parameter measurements have been succesfull for POPC cholesterol bilayer with 60 mol% of cholesterol [Ferreira et al. http://dx.doi.org/10.1039/C2CP42738A], and also quite recent neutron scattering study claims that 61 mol% would be the solubility limit for the cholesterol in POPC bilayer [Garg et al. Soft Matter. 2014 Dec 14;10(46):9313-7 http://dx.doi.org/10.1039/C4SM01219D]. I do not think that this is the reason for the weird simulation result with 60 mol%. Also this (the simulation result) needs some further attention.

      Delete
    10. Interesting cholesterol/PC mixtures: From my point of view the most interesting systems for scattering experiments would be the same as in the order parameter measurements by Ferreira et al. [http://dx.doi.org/10.1039/C2CP42738A], i.e. 7, 15, 34, 50 and 60 mol% of cholesterol in POPC bilayer. The effective temperature in these experiments was 300 K. If this data would be available, we would have the NMR order parameters giving atomistic resolution information, and form factors giving molecular/bilayer level information from the consistent data set. Then, comparing simulations with the same concentrations to these experiments would give a very good quantitative picture about the quality of the models (and if they work, an interpretation for the experiments). However, I think that the data with 50 mol% is the most relevant of these, since with that concentration the lipid-cholesterol interactions affects most to the result, and the quality of these interactions in simulations are not well understood (as far as I know). Even less is known about the quality of cholesterol-cholesterol interactions in the models, thus the 60 mol% data
      (if doable) would be very interesting as well.

      Peter also mentioned DPPC/cholestrol systems. One good option would be to do the same concentrations as with POPC. However,
      there are differences here. There are some old ^2H NMR experiments where phase coexistence observed for DPPC/cholesterol systems (e.g. Vist and Davis [http://dx.doi.org/10.1021/bi00454a021]). It is not clear to me yet if we can extract some order parameters from their data. If we can, then for DPPC/chol systems the same concetrations they have used may be interesting. I will think about this. Anyway, from my point of view the above mentioned POPC/cholesterol data would be the most interesting.

      Other lipids: I do not know the scattering data literature (yet) well enough to suggest any experiments with other lipids. Also, there would a lot to do for us with the above discussed potential data, as well as with the data you have already shared.

      More general quesion: Do you know how accurate is the experimental form factor?

      Delete
    11. Regarding formfactor (FF) amplitudes: Samuli is right; also the amplitude contains information about the density profile, namely about the electron density itself, so for example about the electron density of the phosphate group, but not about the position of it.

      Berger/höltje cholesterol: I agree that the bilayer is to thin, as the area is large. If we could reproduce the bilayer thickening as a function of cholesterol with different models, which is then in agreement with experiments, would be an interesting thing. We will try to measure POPC/CHOL mixtures up to 60% CHOL as this should be possible like other groups has already shown… It is a good idea to measure the same compositions as Ferrreira et al. and additionally I also want to measure more compositions so that we can study the effect of cholesterol in more detail and how is the FF changing as a function of CHOL concentration. The temperature of 300K should be fine.

      Regarding DPPC/CHOL, there is the possibility for phase separation, but it depends on the temperature, and occurs mainly below the melting temperature of DPPC. I thought about the samples we can measure and the priority should be this POPC/CHOL system as there is also a lot of data available and from the experimental point of view POPC is also easier to handle than the DPPC system because we don’t really have to heat the samples. We will measure the POPC/CHOL system carefully and if we have time we can also measure some DPPC samples.

      The accuracy of the experimental FF: The minima of the FFs are very well defined especially for the DOPC and POPC data, which we measured in Hamburg, where the quality of the scattering data was very high (<±0.005Å-1 in the position of the minima). The POPC/CHOL sample was measured in Triest at Elettra and there the uncertainties are higher (~±0.01Å-1), because of additional air scattering and so on. The amplitude is not so well defined as the minima and it also depends on your model. I would estimate ±10% or less. I will try to plot FFs with error estimations, to get a feeling for it. At least for the next experiments on the POPC/CHOL system.

      Delete
    12. I have two questions:

      1) About amplitude: I have not studied this properly, but my intuition about Fourier transformations tells that also the position of, for example phosphate, density peak would affect the amplitude. This is not straightforward, however. Is this discussed in more detailed somewhere?

      2) Which model you are referring when you wrote "The amplitude is not so well defined as the minima and it also depends on your model."? Do you use some model to extract the FF from experimental data?

      Delete
    13. ad 1) You are right that the position of the peak can also influences the amplitude of the toal formfactor but when you look at the expression of the Fourier transform of the electron density, just the cos term contains information about the position and the amplitude is mainly dominated by the absolute electron density. I couldn't find approprriate literature for you but you can have a look at the following paper where they describe the model we are using to calculate the FF including the calculation of the Fourier Transfom (Eq. 8): http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2517047/.

      ad 2) We use the model that is described in the above mentioned paper to extract the FF. I mean by model the choice of a model to describe the electron density profile of your bilayer, which is further used to caluculate the FF. Therefore different models can have an influence on the exact chape of the FF, although it should be minor.

      Regarding POPC/CHOL: In two weeks we have a beamtime at synchrotron in Trieste and there we will measure POPC/CHOL mixtures from 0, 5, 10, 15, .....60mol% CHOL at 27°C. I think this is the temperature at which also the other experiments and simulations were done?

      Delete
    14. Thanks! I will take a more careful look at the form factors.

      Yes, the NMR measurements by Ferreira et al. are effectively at 27°C. The most important thing here is to have the experimental data with the same temperatures so this is the best.

      In NMR the pulse sequence might warm up the sample little bit. I say that the 27 °C is effective temperature since the heating is measured by using methanol chemical shifts, and from this measurement the temperature is 27 °C. So it is the real temperature in the experiment.

      Delete
  7. My original Slipids POPC data was with salt so I ran another simulation with this model and without salt. The initial structure (bilayer of 128 POPCs, 40w/l, 303K) and force field parameters were downloaded from the Slipids home page (http://people.su.se/~jjm/Stockholm_Lipids/Home.html) and used as such. The simulation parameters were also obtained from there and the temperature was changed to 310K. The last 150ns of a 200ns simulation were used for the analysis.

    -0.0221 -0.0403
    0.0582 0.0522
    -0.0194 -0.0193
    -0.0243
    -0.0930 -0.1667

    The trajectory and run input file are available at http://dx.doi.org/10.5281/zenodo.13887

    ReplyDelete
  8. Hi folks,

    here are the data (charmm36) for 0,10,15,20,25,35 and 50% of cholesterol at 283,298 and 308 K.

    https://dl.dropboxusercontent.com/u/11027999/charmm36_popc%2Bchol.tar.gz

    all the systems were built with charmm-gui keeping a total of 512 lipid molecules plus water, I mean for 0% of CHOL we have 512 popc molecules and for 50% CHOL we have 256 popc and 256 chol. All the systems were equilibrated for 70 ns before the 100ns of sampling.

    The data I send you is the result of the calcORDPcharmm.sh script over a pdb file containing 10 frames spaced 10 ns each.

    Best wishes,

    Fernando.

    ReplyDelete
    Replies
    1. Thanks. Did I get it right that you have averaged only over 10 frames? I think that this is not enough data. It seems from your *mdp that you have saved the coordinates between 10ps. It would be good to calculate over all the frames you have (if you have not done it yet).

      Delete
    2. Hi,

      Here are the averaged data using every frame of the simulation (10ps):

      https://dl.dropboxusercontent.com/u/11027999/charmm36_popc%2Bchol.tar.gz

      Best regards,

      Fernando.

      Delete
    3. I have now added your results to the database and figure script, however new figure is not produced yet for other reasons, see https://github.com/NMRLipids/nmrlipids.blogspot.fi/issues/6#issuecomment-86456607.

      The results are very similar to ones reported by Hubert Santuz before, however g2 value is systematically little bit closer to zero.

      Please deliver the missing information into the table 3 in the manuscript.

      Can you also check the simulation details section for CHARMM from supplementary material and report required updates?

      Delete
    4. Hi, here is the requested information:

      https://dl.dropboxusercontent.com/u/11027999/table3.tex

      https://dl.dropboxusercontent.com/u/11027999/simulation-details.tex

      https://zenodo.org/deposit/26928/

      Best regards,

      Fernando.

      Delete
    5. Thanks! I have now updated the information into the manuscript.

      Two remarks:

      1) You referred E. Prabhu Raman et al. when citing to the web page of MacKerell group. Was there supposed to be a citation to some publication? I have now left this out since I did not find any mention from webpage that the parameters were especially made by Raman et al. (let me know if I am wrong).

      2) It would be highly useful if you could also share the full trajectories (contaning all the atoms) and other simulation files through the Zenodo at some point. This has been already done for several systems and the file size problem can be overcame by dividing the trajectory in shorter segments (for example, see http://dx.doi.org/10.5281/zenodo.13285). The advantage of this would be that then the data could be used for further analysis by using the recently developed scripts, see e.g. http://nmrlipids.blogspot.com/2015/03/mapping-scheme-for-lipid-atom-names-for.html. Also the trjajectory (xtc, or trr) and tpr are enough for the analysis purposes.

      Delete
    6. Hi,

      You're right, they just did the charmm36 port for gromacs and as long as I know it is not published.

      Delete
  9. Hi everybody,

    I am providing more simulation data for the Berger force field at various hydrations. I uploaded MD trajectories for DLPC lipids at seven different hydration levels (from 28 to 4 w/l) to Zenodo. They are 80ns long at T=300K.

    Using these simulations (here for 14 hydration levels) I have evaluated the P-N vector angles:

    28 75.57771684011664
    26 75.16919793218835
    24 75.57771684011664
    22 75.94440982900036
    20 75.92435630617078
    18 76.24349379805865
    16 76.55403692301955
    14 76.97172315566992
    12 78.70320161255526
    10 82.44862671932546
    8 92.77046139860725
    6 93.55713245132188
    4 99.90951552593731
    2 101.32013761754939


    Best regards,
    Matej Kanduc

    ReplyDelete
    Replies
    1. Thank you for your contribution!

      I have now analyzed the order parameters from your data. The used scripts and results can be found from:
      https://github.com/NMRLipids/nmrlipids.blogspot.fi/tree/master/DLPCberger/dehydration

      The results are also added to the dehydration figure:
      https://github.com/NMRLipids/nmrlipids.blogspot.fi/blob/master/DATAreportediINblog/dehydration.pdf

      and your results for the P-N vector angle into the PN figure:
      https://github.com/NMRLipids/nmrlipids.blogspot.fi/blob/master/Fig/dehydrationPN.pdf

      I have two observation on these results. The changes in the PN vector angle are much larger than in other simulations, and the changes in order parameters are fluctuating more. The text in the manuscript needs some editing now.

      I did also the mapping file for this model:
      https://github.com/NMRLipids/NmrLipidsCholXray/blob/master/MAPPING/mappingDLPCberger.txt

      and calculated also the acyl chain order parameters (scripts and results in the first link of this reply). The plot is here:
      https://github.com/NMRLipids/nmrlipids.blogspot.fi/blob/master/Fig/OrderParametersDEHYDdlpcCHAINS.jpg

      The acyl chain results seem to be in reasonable agreement with experiments by Dvinskikh et al. http://dx.doi.org/10.1039/b508190d. Related to this, we have to discuss the phase transition issue as you have done in your previous paper.

      I have a couple of questions:

      1) Is this the same data you have used here: dx.doi.org/10.1021/la401147b | Langmuir 2013, 29, 9126−9137? If yes, we can just cite this work in the table 2.

      2) Anyway, we need the total simulation times for the table 2. Could you deliver those?

      Delete
    2. I have now updated the dehydration discussion in the manuscript. The changes can be read from this commit:
      https://github.com/NMRLipids/nmrlipids.blogspot.fi/commit/75874376bff137a730f590e13f94c645d53269ad

      and the most recent versions are in GitHub:
      https://github.com/NMRLipids/nmrlipids.blogspot.fi/blob/master/HGmodelMANUSCRIPT/HGmodel_ACStemplate.pdf

      https://github.com/NMRLipids/nmrlipids.blogspot.fi/blob/master/HGmodelMANUSCRIPT/HGmodel_ACStemplate.tex

      Delete
    3. Thanks for the analysis!
      The order parameters are definitely fluctuating more, which may suggest that either (i) sampling is too poor (80ns not enough?), and/or (ii) samples are not equilibrated enough for these order parameters.
      I am not sure how much exactly were the equilibration times, but no more than 10ns, maybe even less, say 5ns.
      Recently, we have analyzed several auto correlation functions (area, thickness, …) for PC lipids and some relaxation times seem to be larger than hundreds of nanoseconds at low hydrations!!!). Low hydrations are always problematic. One brute-force way to tackle this problem is to produced several samples by dehydrating the lipids via different routes and so avoid sampling of particular trapped state.

      Something we could try right away for the order-parameter analysis, is to throw away additional few nanoseconds of the trajectories at lower hydrations, to see whether this smooths the data a bit (that would indicate whether we have problem with (i)).
      If the problem is in poor sampling, I can also prolong the simulations.

      To your questions:
      1) Total simulations times: 85–90 ns. (5–10ns of equilibration; by throwing away additional few ns would of course change the numbers accordingly).
      2) Yes, these are the same results as in Langmuir 2013, 29, 9126−9137.

      Delete
    4. I have now re-analyzed the glycerol backbone and choline order parameters by using the last 60ns of the simulations you provided. The GitHub and manuscript is updated accordingly. The results did not essentially change. I think that the current results are good enough for the current work. Since the structures in this model are not that good in the full hydration either, it cannot be anyway used to resolve the structural responses.

      Additional question: You write in the dlpc.itp files that you have modified dppc.itp to make them. Was the dppc.itp downloaded from Tieleman's webpage? Is the Langmuir 2013 paper the first one where you use this dlpc.itp file?

      Delete
  10. Indeed, the dppc.itp was downloaded from Tielemna's webpage. It was manually modified for the DLPC lipids for the Langmuir 2013.

    ReplyDelete

Please sign in before writing your comment.