PEG2000 collapsing into bilayer core in atomistic simulation (GROMOS 53A6)

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:

  1. Could this behavior be due to limitations or incompatibility of the GROMOS 53A6 force field for PEG-lipids?

  2. Are there known issues with PEG parametrization in GROMOS that could cause this collapse?

  3. Could my system setup (e.g., initial configuration from Packmol, equilibration protocol, or PEG placement) be responsible?

  4. 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.