GROMACS version:2020.2
GROMACS modification: Yes/No
Hello,
I’m calculating the order parameter of the lipids of my protein-lipid bi layer system using gmx order. For unsaturated lipid chain containing one double bond, i’m getting two dips one at the position of the carbon atom with double bond (9th carbon) and the other at 6th carbon. I’m not sure why the dip appears at the 6th carbon position. I have tried using -unsat yes flag, but its not working as expected.
The gmx order documentation says:
" -[no]unsat (no)
Calculate order parameters for unsaturated carbons. Note that this cannot be mixed with normal order parameters."
How to implement this flag in a correct way to calculate the order parameter of the unsaturated lipid chains?
@TomPiggot will probably chime in here - GROMACS doesn’t compute unsaturated order parameters properly. There are other codes that can do this, as well as proposed patches to gmx order, some of which have been shared on the old mailing list and on the GROMACS Redmine site (which has now moved to GitLab).
Justin is quite correct, GROMACS doesn’t compute the unsaturated order parameters properly. So no matter what you do, don’t use gmx order for the double bond.
First things first though, which lipid force field are you using? If it’s an all-atom one, you probably shouldn’t use gmx order anyway, even for saturated carbons. If you are using a united-atom force field, you will have to use another tool if you want to get the unsaturated carbon-deuterium order parameters. Let me know which of these you are using (all-atom or united-atom) and I can point you in the right direction of other analysis tools, etc. or you can find plenty of examples in the paper linked below.
With regards to the dip at carbon 6, there are a few potential reasons I can think of that might be causing weird results. One is that you haven’t made your index file correctly (although, IIRC gmx order will warn you if you’re doing something that doesn’t seem right). Another is that you may not have enough sampling for the analysis to average out, so either doing the analysis on a really short trajectory or maybe on just one lipid. The final one is that maybe you’re looking at the results produced using the -unsat option for saturated carbons. If you post the results (either the data or a plot) from gmx order calculated without -unsat, I might be able to help a bit more. Plus the exact command you are using for gmx order would help too.
And if you really want to dig deeper into all this sort of stuff, you can find plenty (too much maybe ) of information here: Cheers Tom
I’m using CHARMM36 all atom force field for my simulations. I have attached the plots of unsaturated lipid chain I obtained during the first 40 ns of simulation.
This is the command I used to obtain these plots,
gmx order -s step7_2.tpr -f step7_2.xtc -n popg_sn2.ndx -d z -b 0 -e 40000 -od popg_sn2_twopep.xvg -o popg_sn2_S.xvg
So, as per my first message, I really do recommend using a different tool. Firstly, you’re using an all-atom force field but gmx order doesn’t take the positions of the actual hydrogen atoms into account when doing the calculation (it assumes you’re using a united-atom force field). Secondly, it averages the two different hydrogen atom (pro-R and pro-S) order parameters. Finally, it won’t be able to calculate the terminal methyl group order parameters. Have a quick look at the paper I linked to, it mentions many different tools to do the all-atom analysis. In fact you can quite easily write something yourself to do the all-atom analysis (e.g. I’ve done it before within MDAnalysis).
As to your plot, the overall shape looks fine (bar the unsaturated bits not being right, as expected). If you could calculate the order parameters for the unsaturated carbons they would create a drop at atoms 7 and 8 in your plot. I can see from your naming that you’re simulating POPG, so you should expect the drop primarily at carbon atoms 9 and 10 with 16 points in the plot; there should be 16 data points as the oleoyl tail is 18 carbon atoms in length but carbon 1 is the carbonyl carbon and, as mentioned, gmx order cannot calculate the order parameters for carbon 18.
I think there are two issues are going on with your analysis. Firstly is that the labeling on the x-axis should start at 2 not 1. And secondly, I don’t think you are calculating the order parameters for carbon 2 in this tail (probably just an index group missing). In other words, I believe the data point at atom 1 is actually for carbon 3. But, as I said, you’re better off just using another tool as fixing these still won’t sort the unsaturated groups out.
I should maybe have also said that, as you’re probably aware but no harm in reiterating, what you can quite nicely see from the plots is that with your “multiple peptide” simulation there is definitely more disorder in the central part of the membrane compared to the others two simulations.
Oh and it’s also worth noting that the y-axis should be labelled as -Scd (or |Scd| if, once you’ve calculate the unsaturated carbon atom order parameters they are all still negative).
Thank You very much for all these useful information and suggestions. Actually there’s a published paper in which the order parameter profile of POPG is a plot with same trend as above. We were still doubtful about the dip at carbon 6 which is away from the double bond. Anyway I need to sort out which is the best tool for all atom force field.
I would like to know , finally what did you do?
Currently, I’m working with an unsaturated lipid too, so I’m looking for the best tools for calculating order parameter for these kind of lipids using CHARMM 36 on gromcas.
There really is a huge amount of different programs you could use (I could probably list over 10 just off the top of my head). I guess the only thing you might want to consider is if you’d like to look at the splitting between pro-R and pro-S hydrogens or not, as you can’t automatically get this information from the output of all the tools out there. If not, just pick any (all-atom) tool you like and away you go.
