I have been doing a lot of MD simulation lately but i am confused how to select the final pose for the ligand in the system.
Several approaches I can think of are :
1) cluster the ligand and select the representative pose from the top cluster.
However, I also think that the top cluster may contain poses that were attained before ligand reached its appropriate pose.
2) take final pose after the MD run.
However, that will always change if I extend my MD run
3) take pose when the receptor-ligand system has lowest energy.
However, this pose is present in few frames. While the best pose should be present in large number of frames.
It's a bit tricky, but I would vote for a combination of 1) and 3).
The cluster analysis gives you the poses which are close-enough to each other in terms of your measure. Within one cluster, there shouldn't be large variations of pose conformations. If yes, check your clustering algorithm.
It would also be valuable, if you're able to estimate the energetics of each member of the clusters. Then I would take the member of the most populated cluster having the lowest energy compared to ether members of the same cluster.
It also depends, if you run bare ligand MD or receptor-ligand MD. In the latter case, you don't need to worry about wrong pose conformation which can be "attained before reaching the receptor".
Does final pose mean the final/appropriate structure? IN MD Simulation we get an ensemble of structures possible and one has to do the analysis considering the entire ensemble generated after the system has reached equilibrium
@angana ray:
yes i meant final/appropriate structure.
first of all i am not sure when one can say equilibrium has reached? if i am comparing for a two ligands against same protein i have to do MD for same duration. however, in one case the RMSD seems to have have stabilized and in another it hasn't.
now if i extend the simulation i see the reverse trend.
so how you suggest me to analyse the ensemble and pick one structure from the entire run ?
It's a bit tricky, but I would vote for a combination of 1) and 3).
The cluster analysis gives you the poses which are close-enough to each other in terms of your measure. Within one cluster, there shouldn't be large variations of pose conformations. If yes, check your clustering algorithm.
It would also be valuable, if you're able to estimate the energetics of each member of the clusters. Then I would take the member of the most populated cluster having the lowest energy compared to ether members of the same cluster.
It also depends, if you run bare ligand MD or receptor-ligand MD. In the latter case, you don't need to worry about wrong pose conformation which can be "attained before reaching the receptor".
In order to evaluate if your MD has reached an equilibrium phase.
As preliminary task you can see the time evolution of RMSD and/or energy profile/s of the receptor::ligand system and find if a plateau has been reached.
Use (1). If your clusters contain wildly different poses, then it's just that you need to fiddle with the clustering algorithm.
You can never know if you really reached equilibrium, usually. The best you can do is to go ahead with the simulation as much as you can and see if the pose is stable (it can also flip briefly away from it, what is important is that it goes regularly back to the main one).
If you would believe that there could be few different poses that might explainth the data that you have for your system, you could have a look at few of our publications (PMID: 22928709, 23864731) and check the R library deltaGseg.
@Gloria: I am Interrested in using your R library. Is there available a deltaGseg package version Bioconductor independent?
@Dario. So far the only way of running this type of analysis is using R, basically because we did recycle some old wavelet libraries. If what you want it to to download and install the package independently, it should work.
@Gloria: no problem for using R (we have installed R version 2.13.1) but we have some trouble in upgrading BioConductoR to the latest version...
@Dario, the less complicate option is for you to update your R package to R 3.0. I am not so sure you will manage to install and run this library in older versions of R. Let me know if this is not possible and try to dig an old version which might be of use for you. But only it might!
I agree with Michal Kolar's suggestion of combining 1) and 3). Another thing to consider is, if crystal/NMR structure of your system with a ligand is already known, run a MD trial under same conditions using the known ligand (built from scratch) to see if results are comparable with the reported structure. That validates your method.
@dario thanks for the reply
i have read that RMSD should reach a plateau for a equilibrated system and tried to achieve that. however as i am dealing with multiple systems in one or another system the RMSD starts fluctuating within the range in which the first system reached the plateau. if i extend the simulation this happens in some other system so a true plateau was not achieved in all systems at same time
Also which is better approach
a) calculating RMSD complete backbone of energy minimised structure
b) calculating RMSD against backbone of selected residues in proximity of the ligand, lets say within 7 angstrom of the crystal or energy minimised structure
@tobias thanks for the reply
i ran 20ns and then extended to 50ns and now thinking to extend by another 50ns to make a 100ns trajectory
@michal thanks for the reply
i agree with the approach of picking up lowest energy structure from the ligand and will like to apply. currently the method pics the pose with lowest RMSD against the poses within same cluster as the representative pose.
I am using receptor-ligand system so can please explain what you meant by the following "In the latter case, you don't need to worry about wrong pose conformation which can be "attained before reaching the receptor".
@praveen thanks for the reply
i am using unknown ligands in order to see the binding pose and maybe explain the activity of the ligand so can't compare with crystal/NMR pose
@gloria thanks for the reply
i will look into the papers mentioned by you and thanks for suggesting them
@nitin try further investigations...
Focusing on residues around ligand should be a better approach for post-processing analyses.
As I wrote, RMSD profile monitoring (over backbone or Ca atoms) is a preliminary task and is not unusual to find different fluctuations among diverse systems.
Another method to evaluate your MD should be the comparison of collective motions between different pieces of simulation. You have to eigendecompose the covariance matrix of atomic coordinates displacement, then to compare the trajectories projected against the first representative eigenvectors. Take a look at GROMACS g_covar and g_anaeig commands.
I understand that the problem arised during your MD run. As the input .rst file for MD might have been created by EM run, posibly error might have happened there. You may check up with your EM run also. You can try with increasing no. of cycles in EM also. .rst file of MD run consists of coordinates as well as energy information also. You will get more information at http://ambermd.org/formats.html
You can also go through the AMBER manual, which can be downloaded from the website,
- Badges
- Science topic
- Similar topics
- Chemistry
- Theoretical Chemistry
- Computational Chemistry
- Molecular Dynamics