PBC leads to high ligand RMSD after unbinding

GROMACS version: 2024
GROMACS modification: No

Hi everyone!

I’ve run a 100 ns simulation of a Protein-Ligand complex. I used the Charmm36 forcefield, prepared my protein with CharmmGUI to correctly deal with a Heme that the protein has, and finally prepared my ligand with CgenFF. The simulation ran fine, but I am having issues with the analysis now.

Here are the plots for the interaction energy and the ligand RMSD after fitting the protein backbone.


As the plot shows, at first the ligand is bound to the pocket (Purple ligand on the right in the picture below). Then, the RMSD of the ligand spikes to 7 nm at the 60 ns mark, which is the unbinding event, however the 7 nm I believe is a periodic shift, which I was able to confirm by checking the frames in PyMOL (White ligand, left). The ligand then returns and binds near to the initial binding region (Yellow, right), before a second unbinding event ocurrs at around 85 ns. This time however, instead of returning to the pocket, it binds to the other side of the protein (Green, left).

While the unbinding itself is simply a result, the second binding to the other side of the protein is an artifact, right? Should I do something to avoid this “artifact”? If so, what can I do? I read the forum and found a suggestion to use trjconv with -pbc nojump but given the ligand interacts with the protein in this other region, it is not useful. Maybe make a bigger box? I used a dodecahedron for the box with the distance -d of 1.0.

Any help or tip on how to fix this is greatly appreciated!