GROMACS version: 2018
GROMACS modification: Yes/No
Hello everyone,
I am performing atomistic molecular dynamics simulations of a lipid bilayer system composed of an ionizable lipid, phospholipid, cholesterol, and PEG-lipid (PEG2000). Since my ionizable lipid is not available in CHARMM-GUI, I constructed the bilayer manually using Packmol.
For the simulations, I am using the GROMOS 53A6 force field with GROMACS 2018. The system preparation protocol includes energy minimization followed by multi-step equilibration: 4 steps of NVT and 4 steps of NPT, during which position restraints are gradually reduced from 1000 to 0 kJ/mol/nm².
However, once all position restraints are removed, I observe that the PEG chains (–OCH₂–CH₂– units) collapse into the bilayer core instead of extending into the solvent phase. This behavior persists even after a 500 ns production run.
According to the expected structure of lipid nanoparticles (LNPs), PEG chains should remain solvated in water and form a hydrated corona that provides steric stabilization and prevents aggregation. Therefore, this behavior seems unphysical.
I would appreciate guidance on the following points:
-
Could this behavior be due to limitations or incompatibility of the GROMOS 53A6 force field for PEG-lipids?
-
Are there known issues with PEG parametrization in GROMOS that could cause this collapse?
-
Could my system setup (e.g., initial configuration from Packmol, equilibration protocol, or PEG placement) be responsible?
-
Are there recommended best practices for modeling PEGylated lipids in bilayer or LNP systems?
Any suggestions or references would be greatly appreciated.
Thank you in advance for your help.