GROMACS version: 2024.4-Homebrew
GROMACS modification: No
Dear all,
I’ve come across many discussions about high RMSD values, sudden jumps, and fluctuations. I’ve been scrolling through them and trying various trjconv commands on my system, but I still can’t get a smooth RMSD plot. At this point, I’m not even sure if my system has crashed. I’ve conducted MD simulations for a protein-ligand complex before, and the RMSD values were smooth except for this most recent system I set up.
Here is the RMSD result of my system:
Note: Complex graph refers to TGDSPF_heavy after lsq fit to backbone
This is a protein/ligand system right? By the look of it, it doesn’t seem a PBC problem to me, but rather that the RMSD of the complex - and only the complex - is very large, so probably your ligand has left the binding pocket. Did you visualize the trajectory?
Yes, I did visualize the trajectory and this is a protein/peptide system. Does this mean the results don’t indicate any potential issues or faults in the simulation?
However, when I analyzed the gmx distance, the average distance exceeded 0.35 nm. This confuses me because, according to the tutorial, a hydrogen bond should not be longer than 0.35 nm. Additionally, I’ve tried several approaches, including the no jump module and others I’ve come across, but the results remain unchanged.
Not necessarily. If the visualization shows that the peptide is leaving the starting configuration, then there isn’t any post-processing tool that can help you, simply because that’s what happened and what the simulation shows. Why this happened is another story and might not be simple to figure out. Was the peptide in a good starting configuration? Was the minimization/equilibration performed properly? Is this a peptide known to bound tightly/lightly? Is this a crystallographic (or somewhat experimental) pose or does it come from docking/AI/etc? Is the force field suitable for this kind of simulations? How about the parametrization of the peptide?
My expertise lies kind of far away from peptide/protein systems but these are just a few questions that I think could help you figure out what is going on. If you are also not very familiar with these systems I highly encourage you to take a deep dive in literature and see if anyone has already done something related to your systems and see what do they do (parametrization, methods, etc.).
A couple of other points though. Firstly, the distance you cite is computed with respect to what? The average distance between the center of mass of the peptide and that of the protein can easily be way larger than 0.35nm. This doesn’t tell you anything about the system, as the protein can be large and the peptide can be binding ‘far away’ from the center of masses. If you want to take a look at hbonds specifically I suggest you to keep you selection of distances only within the atoms which you expect to participate to the hbond. Even better, use proper tools like gmx hbond or external libraries like mdanalysis to carry out in a simpler way these checks.
Secondly, the RMSD seems already large at the beginning of the plot which I suppose is the production run? Is the peptide already not bound at the beginning of the production? In that case, it left during the equilibration. You might need a better/more delicate equilibration procedure, with decreasing restraints etc. Just visualize also the full equilibration and see where the peptide leaves and try to rationalize why.
Lastly, consider that also sometimes bound molecules can just unbind. How likely this is depends mostly on how tightly are they bound, but in general it is possible. If this is a weak binder then maybe this makes sense. I think you might need more than one replica to also be more confident about what is happening before arriving at any conclusion. :)
This has been incredibly helpful for me as a beginner in this field. I am currently researching anticancer peptides derived from milk proteins, specifically from Mare’s Milk. My research is conducted entirely through dry lab approaches, meaning I perform digestion simulations, peptide screening, molecular docking, and, finally, molecular dynamics simulations.
Before conducting molecular docking, I constructed the screened peptides using Marvin Sketch. Surprisingly, I noticed something unusual in the docking results—specifically, bond-order changes. For instance, the C=O bond changed from a double bond to a single bond. I am unsure why this occurred. For context, I performed blind docking using PyRx and optimized the ligand using Open Babel. This issue made me question whether PyRx is suitable for peptide-protein docking. Any suggestions would be highly valuable!
When I proceeded to ligand parameterization using CGenFF, I encountered penalties higher than 10. However, I chose to move forward with molecular dynamics simulations despite this. Interestingly, I did not face any issues during equilibration or MD production.
Thank you for your suggestions! I plan to explore this further and conduct another replicate. Lastly, would you mind if I reach out for guidance on setup and simulations in the future?
You are free to reach out for some help in setting up simulations. However, as I said, I am not expert at all of ligands/peptides simulations, and as such I can’t help in the meaning of the results, or whether some docking software/pose makes sense. You might still find help from other users of the forum that have more expertise in this!
Regarding CGenFF, the penalty scores should be related to how confident is the force field in the parameter, i.e. how similar is that chemical bond/angle/etc with respect to what the force field knows from its set of parameters. This means that a high score is not necessarily bad or good, but it should get further attention if you plan on using that molecule to do computational studies. As such, I wouldn’t expect any issue from the molecule during equilibration/production, but this doesn’t mean that the parametrization is good, just that it is numerically stable (which in general is expected from such official tools). Your best bet is checking in literature what has been done and try to understand if it makes sense to follow their same method/parametrization/etc.
Dear all, I am trying to perform MD simulations with GROMACS from an already docked file with a protein receptor and a small peptide with PyRx but when creating the input files using CHARMM-GUI, I am encountering lots of problems related to the structure of the docked PDB. Could you please indicate which steps did you do to perform your simulations or if you used any other approach to generate the input files?
