GROMACS version:2025
GROMACS modification: Yes/No
Here post your question
Hello everyone,
I am currently setting up a membrane protein simulation using the CHARMM36m force field and would appreciate some guidance regarding the unrestrained equilibration stage.
After completing the restrained equilibration protocol, I removed all position restraints and began the unrestrained equilibration at 310 K. My uncertainty is about which pressure coupling scheme should be used at this stage before (or for) the unrestrained equilibration step.
The options I have considered are:
NPγT (constant surface tension)
NPT with semi-isotropic pressure coupling
NPT with isotropic pressure coupling
To compare them, I performed simulations using different pressure coupling schemes and observed the following average area per lipid (APL):
NPγT: ~55 Ų
NPT with isotropic pressure coupling: ~64 Ų
The isotropic NPT simulation gives an APL that appears to be closer to the expected value for POPC at 300K for my lipid bilayer, whereas the NPγT simulation results in a noticeably smaller APL.
What has confused me is that I have found several published membrane protein studies using the CHARMM36/CHARMM36m force field that report running their unrestrained equilibration or production simulations with isotropic NPT, while many recommendations and tutorials suggest using semi-isotropic pressure coupling. Older literature also discusses NPγT simulations for CHARMM lipid force fields.
I have attached the area per lipid (APL) and membrane thickness plots from both simulations for comparison.
Am I missing something in my simulation setup or interpretation? Which pressure coupling scheme is currently considered the best practice for the unrestrained equilibration of membrane protein systems with CHARMM36m?
Any insights or references would be greatly appreciated. Thank you!
regarding your observation of the APL under NPγT and NPT. If you impose a surface tension in your simulation (NPγT) the APL will depend on the magnitude of γ. If you compress your membrane area, the APL will be lower and if you stretch your membrane area the APL will be higher. Using a semi-isotropic pressure coupling leads to no surface tension and the APL depends mainly on the temperature (i.e., is the membrane in a fluid or gel state) and the lipid type(s) in your membrane (neglecting influences of ion concentration, force fields, etc.).
If there is a protein inside the membrane, the equilibration can become quite tricky, because you have to keep an eye on the protein and the lipids in the membrane. I will assume in the following that you have a setup with a “well-behaved” protein (e.g., no missing residues, proper protonation states, no large structural clashes, etc.). In general several equilibration steps with no pressure coupling followed by semi-isotropicpressure coupling is a robust workflow. An example is the workflow included in the output of the CHARMM-GUI web server:
Minimization
NVT (dt=1fs)
NVT (dt=1fs)
NPT (semi-isotropic, dt=1fs)
NPT (semi-isotropic, dt=2fs)
NPT (semi-isotropic, dt=2fs)
NPT (semi-isotropic, dt=2fs)
For each equilibration step you impose position restraints on the backbone & side chain of the protein, on the lipid head groups, and some of the lipid dihedrals to ensure correct chirality and no flips from trans to cis in the unsaturated bonds of the lipid tails. For more details, I would recommend that you either prepare your system using the CHARMM-GUI or a similar system to obtain the described .mdp files. I would say that nowadays a lot of people rely on this workflow and it worked for me (I use it mostly for pure membranes) several times pretty well. However, this does not necessarily be the case for your system (chances are low, but you never know…), in this case I would recommend that you slowly tweak the workflow (e.g., increasing force constants of the positions restraints or similar).
For the final production simulation, I would recommend semi-isotropic pressure coupling this imposes zero surface tension on your system and the area per lipid/protein structure can proper equilibrate.
With isotropic NPT the APL will be set mainly by what you “happen” to have set as initial unit-cell dimensions and density. This is never a good setup, unless you have been extremely careful in setting the dimensions and density. Those are usually obtained through independently coupling x/y and z pressure components. Or by fixing the x/y dimensions, if the force field doesn’t give the correct lateral pressure.
Thanks for the reply. I have tried an unrestrained simulation using an NPT semi-isotropic setting at 310 K. During the simulation, the APL values over the last 20 ns were 60.90 Ų and 58.72 Ų for the upper and lower leaflets, respectively.
I have seen studies in which the system reached equilibration within approximately 20 ns. I used CHARMM-GUI for system preparation, followed by the standard minimization and equilibration protocol provided by the CHARMM-GUI-generated script. For APL calculations, I used the FATSLiM tool.
Does this indicate that I need to perform a longer equilibration, or is this lower APL difference/behavior typical for POPC using the CHARMM36 force field?
there are several points that you need to consider:
You stated that you simulate a protein-membrane system. If the protein does not have a fully symmetric shape but is instead rather asymmetric (e.g., cone-like), the APL can be slightly different in the upper/lower leaflet. It depends on how many lipids CHARMM-GUI placed in the upper/lower leaflet and how the protein behaves inside the membrane. This can also prolong the equilibration time of the APL.
Based on your plot, it looks like the lipids were getting closer together in the first 10 ns of the simulation (reduced APL). Afterwards, the APL fluctuates, but there is no visible large drift anymore. That does not mean that your APL is equilibrated — it could still drift on a longer time scale, or the protein could abruptly change shape after 200 ns, again altering the APL. How did you choose your coloured frame?
There are some weird spikes in the upper leaflet APL, probably an artefact of the FATSLiM tool. It would be worth investigating what’s going on there.
The selected range for calculation is the final 20ns of the run. i will also try APL calculation with other tools to check about the spikes that is seen in the plot.
As per the charmm gui mdp temperature coupling is given to three diffrent groups, should i change it to only one thermostat coupled to system?
based on some other threads here in the forum it is apparently recommended by now to use only one thermostat group (i.e., the whole system). However, I think most people still use the setup with two thermostat groups (i.e., Membrane+Protein, Solvent). Both are probably fine if the groups are large enough.
The APL is an important observable but can be tricky to calculate, especially around a protein. However, you could also monitor the length of the box vectors that gives you also a good hint if the system density equilibrates.
it is good to see that the shape of the curve from your box dimensions resembles roughly the shape of the APL curve, i.e., the rapid dip in the first 10 ns. However, it is hard to judge if the membrane area has equilibrated or not. I think that a longer simulation will tell you if the membrane area has reached a stable plateau or not. From your box dimensions it seems that the system is not very large (7x7nm), for such systems several hundred nanoseconds to several microseconds are nowadays standard. Please keep in mind that even if your membrane area (or APL) reached a stable plateau, this does not imply that other observables of interest are equilibrated as well!
I have rerun the minimization and equilibration process according to Ron O. Dror protocol where
Minimization: 3 rounds, each with:
500 cycles Steepest Descent (SD)
500 cycles Conjugate Gradient (CG)
Round 1: Harmonic restraints of 10 kcal mol⁻¹ Å⁻² on protein, lipid, and ligand.
Round 2: Harmonic restraints of 5 kcal mol⁻¹ Å⁻² on protein, lipid, and ligand.
Round 3: Harmonic restraints of 1 kcal mol⁻¹ Å⁻² on protein and ligand.
Heating 1:0 → 100 K, NVT, over 12.5 ps, with 10 kcal mol⁻¹ Å⁻² restraints on protein and ligand heavy atoms.
Heating 2:100 → 310 K, NPT, over 125 ps, with 10 kcal mol⁻¹ Å⁻² restraints on protein and ligand heavy atoms.
Equilibration:310 K, 1 bar, NPT.
Restraint release: Starting from 5 kcal mol⁻¹ Å⁻²:
Decrease by 1 kcal mol⁻¹ Å⁻² every 2 ns for 10 ns.
Then decrease by 0.1 kcal mol⁻¹ Å⁻² every 2 ns for 20 ns.
Total restraint-release equilibration:30 ns.
Then i tried an unrestrained NPT semi isotropic simulation for 50 ns and i got the APL as follows Upper APL 62.99A and Lower APL of 62.047 is it ok to run this protocol rather than those given from charmm gui.
Sorry for my late reply. The workflow sounds quite careful/conservative to me, more than the typical CHARMM-GUI workflow. You should be probably fine, although it is hard to judge without knowing the system. However, remember that your final unrestrained NPT semi-isotropic run (a.k.a, production simulation) should be definitely longer than 50 ns if you want to measure meaningful ensemble averages and corresponding error bars.
