GROMACS version:2021
GROMACS modification: No
Hello,I am used in performing MD protein-ligand simulations in Gromacs and now I would like to learn umbrella samplig. I have tried this tutorial Umbrella Sampling
but it isn’t about protein-ligand systems, so I don’t have proper results. Has anyone any recommendations or examples about this topic?
Hi Alex, as a starting point, if you were able to do the tutorial all the way through and your results match the tutorial’s graphs in the data analysis section, that should be a pretty good indicator that you’re on the right track. Is there a particular step that you’re unclear on or is causing issues with the simulation?
I have done the tutorial firstly.
My problem is in the pulling and then in the configurations.
I tried pulling using the same options as in the tutorial, but as group 1 the protein and as group the ligand and also the pulling to be done on the y axis (I have created a big box, bigger on y axis, similar to tutorial). However, in the xtc file, the ligand moved insignificantly, definetely not 5 nm. On the configuratins the http://www.mdtutorials.com/gmx/umbrella/Files/get_distances.sh script couldn’t calcuate the distances, my problem is that in the tutorial it is not named which is the index.
Thank you for your time
To clarify, you’re having difficulty pulling the ligand from the protein by the desired distance, and can’t get the get_distances.sh script to work due to the script using different indices, is that correct?
On pulling the ligand, different proteins might need different forces to pull apart. The simulation pulls the protein by attaching a spring, then pulling on that spring at a certain rate. If you want to edit the mdp file to get a larger distance, you can increase the force by either increasing the stiffness of the spring (pull_coord1_k) or increasing the rate at which you pull (pull_coord1_rate). You could also run the simulation for a longer period of time (nsteps), so the spring has more time to pull the protein and ligand apart. Keep in mind that while increasing the force could pull apart the ligand and protein faster, the protein under high force might not behave as it would in a real physical system, so how you adjust the system is gonna depend on what you want out of the simulation.
On the index used, when you made index.ndx, did you name the ligand and protein something other than ‘Chain_A’ and ‘Chain_B’? get_distances.sh is calculating distances for ‘com of group “Chain_A” plus com of group “Chain_B”’, so if you change the names of the groups in index.ndx and get_distances.sh to match, then it should work.
I can quickly weigh in that a protein-ligand complex is quite different from the tutorial.
One should define the pulling vector (reaction coordinate) as some group of residues in the binding/active site and the ligand itself, such that the vector connecting them projects out along a path that exits the binding/active site. This is non-trivial and highly dependent upon the system. It may also not be possible for all systems.
The bias should be applied in all 3 spatial dimensions, not just along one as in the tutorial (which is a very specific case in which this works).
Thank you for help.
Forgot to mension, my system is from molecular docking.
If i get it correct, I should define on chain A, a residue that is on the way of the ligand to exit from the binding site and on chain B the ligand. Also, I should change the distance into direction in pull_coord1_geometry. Apart from that, should I make any other change in the pull code?
; Pull code
pull = yes
pull_ncoords = 1 ; only one reaction coordinate
pull_ngroups = 2 ; two groups defining one reaction coordinate
pull_group1_name = Chain_A
pull_group2_name = Chain_B
pull_coord1_type = umbrella ; harmonic potential
pull_coord1_geometry = direction;
pull_coord1_dim = Y Y Y
pull_coord1_groups = 1 2
pull_coord1_start = yes ; define initial COM distance > 0
pull_coord1_rate = 0.01 ; 0.01 nm per ps = 10 nm per ns
pull_coord1_k = 1000 ; kJ mol^-1 nm^-2
To my understanding, if you’re using pull_coord1_geometry = direction, you also need to define pull_coord1_vec (https://manual.gromacs.org/documentation/2018/user-guide/mdp-options.html#mdp-pull-coord1-vec) such that the pulling is in the direction the ligand takes to exit the binding/active site.