Suppose a molecular dynamics (MD) simulations of a protein-ligand complex is performed (say 10 nano second). I want to identify two amino acid residues which undergoes coordinated motions to bind ligand (in other words, two residues cooperate to enable ligand binding). Which one/the exact parameters of MD should I introspect to identify these two residues?
First of all, a practical and important note: I'm not sure whether 10 ns of simulation is adequate sampling for your complex. You may need to increase your simulation time to let your system evolve to a state of equilibrium and then check to make sure that it has converged. Making asumptions based on the analysis of unconverged simulation data can be VERY misleading.
As to your question now, from what I understand, you need to study the collective motions of your system. Collective motions can generally be obtained through Principal Component Analysis (PCA) of your simulation results, a procedure that is also known as "Essential Dynamics Analysis" (EDA).
Very simply put, in PCA/EDA the configurational space of a simulation system (e.g. a protein) is mathematically/statistically "divided" into two subspaces, a Gaussian-like one, containing all fast fluctuations (i.e. thermal motion and fast vibrations) and an "essential" one, containing the slowest motions that are also believed to be biologically relevant for the system of interest. This way, you reduce the number of dimensions needed to describe protein dynamics to only a few, corresponding to this essential subspace.
This is achieved by computing the correlations of motion between the protein's atoms or residues (depending on the approach), the calculation of mobility covariance as a matrix and the diagonalization of said matrix through PCA. The end result is the description of the system as a series of eigenvectors, the first few of which can describe the "essential" subspace of "important" motions.
More importantly, in your case, the data from these "essential" eigenvectors can also be used to derive information on various aspects of proteins, including the correlated mobility of its components (basically, an RMSF-like measurement, only this time it is performed along the eigenvectors instead of just the trajectory). In other words, by analyzing the EDA/PCA results, you can find information to identify the segments of the protein that move in concert with each other during your simulation.
You should read the original paper of the method by Amadei et al:
https://homepage.univie.ac.at/mario.barbatti/papers/method/essential_dynamics_1993.pdf
as well as some more recent reviews on the matter to familiarize yourself with the topic. This is a rather good one:
Article Principal Component Analysis: A Method for Determining the E...
In more practical terms, now, how you do PCA/EDA depends on what MD simulation software you use. GROMACS has built-in tools for this job ("gmx covar" to build the covariance matrix, "gmx anaeig" and a few others for subsequent analysis), and in AmberTools you can do the same with cpptraj (among the other things that cpptraj can do).
Other great tools include carma (https://utopia.duth.gr/glykos/Carma.html), ProDy (http://prody.csb.pitt.edu/index.html), a library for Python, MDanalysis(https://www.mdanalysis.org/), again, for Python, and Bio3D (http://thegrantlab.org/bio3d/), a package for R. These also support PCA and can be used to analyze trajectories produced by CHARMM, NAMD, AMBER etc.
@Fotis Baltoumas I agree with your response. Should DCCM could be used instead of PCA?
I agree with Prasanth Kumar about the analysis. Try producing the dynamic cross-correlation matrix (DCCM). There is a nice R package for this analysis called "Bio3d".
You can calculate DCCM or LMI (Linear Mutual Information) for correlation on the condition that your simulation is well-equilibrated (as mentioned above). Bio3D is a nice package to do this. For LMI, one can use g_correlation of Lange and Grubmuller or Wordom as well.
- Badges
- Science topic