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.
