[Update 12.11.2014] Originally this post claimed that the GAFFlipid model has the wrong isomer for the glycerol backbone. Further studies revealed that the problem was actually the initial configuration downloaded from http://lipidbook.bioch.ox.ac.uk/package/show/id/150.html having different isomers in different leaflets (see comments to this post). After setting up new simulations with correct initial stucture also the GAFFlipid has the correct stereospecifity for the g\(_1\) segment:
[Original post begins here:]
Based on discussions in the blog and the experiments by Gally et al. we have now reported the order parameter results by taking the stereospecifity of g\(_1\) and g\(_3\) into account. The result was already reported in the Current status of the project on 29.4.2014, and it is shown here again in Fig. 1.
 |
| Figure 1. Order parameters for choline and glycerol for POPC. Magnitudes measured by Ferreira et al., signs measured by Hong et al., Hong et al. and Gross et al. , and stereospecific labeling measured by Gally et al. |
[Original post begins here:]
Based on discussions in the blog and the experiments by Gally et al. we have now reported the order parameter results by taking the stereospecifity of g\(_1\) and g\(_3\) into account. The result was already reported in the Current status of the project on 29.4.2014, and it is shown here again in Fig. 1.
 |
| Figure 1. Order parameters for choline and glycerol for POPC. Magnitudes measured by Ferreira et al., signs measured by Hong et al., Hong et al. and Gross et al. , and stereospecific labeling measured by Gally et al. |
The main conclusion from the Fig. 1 was that the GAFFlipid model has different stereospecifity in g\(_1\) order parameters compared to CHARMM36, MacRog and experiments. While writing the new version of manuscript I was plotting the lipid structures from different simulations and noticed that the GAFFlipid has different isomer in glycerol compared to other models. This is demonstrated in Fig. 2
 |
| Figure 2. Snapshots of lipid structures from simulations with different models. All the structures are shown such that the viewer is looking at the direction of g\(_1\)-g\(2\) bond. In Berger, MacRog and CHARMM36 the g\(_2\) hydrogen is pointing to the right. In GAFFlipid it is pointing left. |
In conclusion, the different isomer in GAFFlipid explains the different R/S labeling observed in simulations. I did also check the isomer in Lipid14 and it seems that the issue has been fixed there (Lipid14 had the magnitudes more off from experiments, however, as seen in the The lipid forcefield comparison post).
The g2 carbon is asymmetric, and is S- in the case of the GAFFlipids and R- for the rest. Unfortunately, since I have not provided any coordinates for tleap, the wrong isomer was generated. The other isomer can be generated by inverting the sign for one of the columns containing the atomic coordinates. It would seem to me that this would fix the sign problem for the order parameters, but I guess it should be checked.
ReplyDeleteI have been looking at this issue in more detail while writing the paper. I realized two days ago that in our GAFFlipid simulations actually half of the glycerols had the other isomer and half had the other one. When looking this more carefully, I saw that the different isomers were in different leaflets. Then I went back to the structures in the lipidbook (http://lipidbook.bioch.ox.ac.uk/package/show/id/150.html) and saw that the structure shared there has the same problem, i.e., different isomers in different leaflets.
DeleteSo it looks to me that either we have inherited these isomers from the lipidbook, or then there is some systematic problem somewhere that generates these kind of structures. Note that both isomers are stable during the simulation.
It also seems to me now that, most likely, fixing the isomers in the initial structure would lead to reasonable order parameters. Consequently, the force field itself would be quite reasonable. This would mean that we should bring the GAFFlipid model back into consideration.
I didn't know that you used the membrane from the lipidbook website to perform the simulations for the GAFFLipids. The differences in the two leaflets are due to the procedure used to generate the bilayer. Actually they state in the paper that "to obtain an initial configuration, a single lipid molecule was replicated along the x and y axes, and flipped along the z". Based on what I wrote previously, I let you guess what happens to the chirality of the g2 carbon when you do this :-). If you need I can make a bilayer using packmol, but it would need some equilibration. Otherwise you could do the anaysis using only the correct leaflet.
ReplyDeleteI think that we should do it with again with the correct isomers even though I think that the results are the same. Otherwise we have to write complicated explanations into the publication. It would be very useful if you would do a new starting stucture with the correct isomers.
DeleteInstead of creating a new bilayer using packmol, I used a bilayer from the Slipids website, renamed and reordered the atoms, and added the file POPC_303K.gro to the https://github.com/mretegan/GAFFLipids2Gromacs repository. I also rechecked the parameters to make sure that nothing went wrong in the conversion from the AMBER to GROMACS format. The energies are very similar between the two packages.
ReplyDeleteHowever, some of the values for the angles involving the phosphorous atom are different, and I can only assume that something has changed in the GAFF forcefield between AMBER12 and AMBER14. I would suggest to run the simulations using the previous POPC.itp file.
I have been running simulation with the new initial structure and while analyzing it I noticed that there is one molecule (residue number 109) which has the wrong isomer. It is the same situation in the POPC structure distributed in the Slipids page:
Deletehttp://people.su.se/~jjm/Stockholm_Lipids/Downloads_files/POPC_303K.gro
It seems that we have to be extra careful with the initial structures from various sources if we want to have correct isomers.
I think that one molecule in wrong isomer does not matter in practise, but I will just remove the molecule and one molecule from opposite leaflet to be consistent.
I have now ran the simulations and calculated the order parameters from the trajectory having all correct isomers. All the related files including order parameters in ASCII format can be found here:
Deletehttps://github.com/NMRLipids/nmrlipids.blogspot.fi/tree/master/POPCgaff
I have added the new results to the plot:
https://www.dropbox.com/s/qztuhri1xhzam5o/HGorderparameters5.jpg?dl=0
I think that we can conclude now that the incorrect isomers from GAFFlipid were from the intial structure, not from the force field. I will try to write comment on this also to the lipidbook.