I am facing weird molecular dynamics data. The mean RMSD for 312 K lies below the mean RMSD for 310, this seems to be impossible since with higher temperatures, the RMSD should be greater. Any explanations?
Your result best epitomizes the difficulties with MD. People who do MD for living consider it very well based science. Many valuable results were obtained using it. You can use it but seek help from more experienced coworkers. In essence what you obtained is nonsense but nonsense that is expected from a weak method with weak methodology. Even if you carefully prepare your structure before the simulation (careful removal of offending conformational errors, minimization, solvent and counterions addition and initial prequilibration) you might still expect problems. My opinions, after using the mothod in many contexts for a long time, are shaped by my background which is a combination of in depth training in physics (PhD) and 30 years of experience in experimental protein work (crystallography and everything else).
MD is just a computer exercise (like school work) and can become science only in the hands of a very experienced user. There is nothing objective nor repeatable in doing MD simulation, therefore on its own right this does not qualify to be science. The lowest common denominator of all science methods are principles of repeatability and independence of the person who is doing the experimentation. To fully qualify as science it also needs to have a certain element of science inference which is predictability.
MD fails every and any of these principles without a very careful user. Nevertheless all the leaders in the field are so strongly attached to marketing their contributions that most of the fundamental problems of MD are swept under the carpet.
Main difficulties are:
1) MD is a classical physics (not quantum) approximation for movements of quantum objects. Everybody using the method must understand the limitations it brings. Car-Parinello method is trying to deal with that but it is very computer hungry.
2) Forcefield weaknesses, simple central harmonic forces with crude angular and torsional representations.
3) Approximations neglecting polarizabilities and quantum evolution of charge (people are working on it) and ionization (harmonic approximation does not allow for ionization.
4) Very crude dynamical treatment of time evolution.
5) As a strongly coupled nonlinear system it is chaotic and by definition not PREDICTABLE (It has been show in several papers to have Lyapunov instability and exponential divergence).
These would be sufficient to have limited trust in widespread use of MD but additional elements compound its difficulties.
6) Almost every protein structure in the PDB has errors. Even PDB Redo does not remedy that. A very careful evaluation of the structure before the runs should be mandatory using tools developed by Dr Richardson but it is not and most people using MD does not even check basic stereo-chemistry. I have seen amazingly bad stereo-chemistry in many results of MD simulations. So placement of all protons in the structures should be mandatory as a tool for spotting the biggest errors. To give you an example. A single incorrect conformer of a medium size residue (such as Leu) can carry as much energy to explode your structure. (Average folding energy ~35kcal, compare to and average conformational error of a single torsional angle is around 9kcal). Of course people developed remedies against it but at the cost of distorting the structures (tighter constraints).
7) Questionable treatment of energy flows during simulations (temperature as well as pressure control). A standard velocity rescaling violates energy conservation only to restore the energy. A standard procedure to remove the movements of center of mass is violating energy conservation only to restore the energy. Standard methods for heat bath produce energy ad nihilo only to restore the energy.
So in summary, assigning any significance to the very large RMS difference between two runs is utterly meaningless. As another advise, you should go very careful through preparatory stages and repeat runs multiple times to really say anything substantive.
Besides, 6Å RMSD is just simply a sign of something bad especially that more than 30% is contained in the very first 200 ps.
Sometimes it is reffered to as an artefact. But before considering it as artefact, first try to find RMSF and B-factor of your experiment. Either your protein is a dock-complex or simple, you need to know which amino acids are participating in this kind of behavior by analyzing trajectories of individual amino acids. This could provide you required answer.
Is this related with more . H - bond interaction in 312 K?
It is a dock - complex protein, and the 312 K complex has more interaction than 310 K. I am now going for the RMSF and analyzing the trajectories of individual AA based on your guidance.
There are many other interactions except H-bond in a dock-complex. So you may need to verify interactions of ligand with potential amino acids. Following paper may guide and help you a bit
S.S. Azam, et al., Structure and dynamics of alpha-glucosidase through molecular dynamics simulation studies, Journal of Molecular Liquids (2012), doi:10.1016/j.molliq.2012.07.003
You should check the velocity of the center of mass. If you don't have blocked it during the beginning of the simulation, a great portion of your kinetic energy could be contained in the molecular rototranslation (flying ice cube effect)
MD simulation of 5 ns does not provide a conclusive results, in addition the difference in temperature is not very large. It might occure due to incorrect equilibrations and minimizations. You should re run dynamics to check the problem if there is.
I agree with some comments above that it is better to check RMSF if we want to see the effect of temperature. RMSD gives us the idea of how similar of our system with respect to the initial configuration. I think RMSD depends on whether or not the initial configuration is close to equilibrium at that temperature.
The 2A difference occurs at the very beginning of the simulation. And both two simulations underwent dramatic change in a very short period from the beginning. You may want to check your initial structure first. Have you done energy minimization? It can also be the water molecules around the protein if you used explicit water models. That means you have to run a position-restrained MD before this production MD. After all that, I would suggest to slowly increase temperature to the target.
As others have mentioned, 5 ns is too short even for a
I assume these are two simulations with the same complex and the only different of the input files are the temperatures, right? You may double check it.
We don't think dT= 2 K will make such big difference, so try repeating it at several temperatures as suggested by others then you might get some idea why it happened. What is the force field, is it in explicit or implicit solvent? What is the simulation package?
MD is an stochastic process, 5ns is truly not enough data sampling to get a good idea. Repeat your simulations several times and analyze the variability of the RMSD profile. Most likely, you will find 312 and 310 being very similar.
From your data, it looks your simulation rapidly found an energy well in the 312K MD and explored some alternative conformational ensemble in the second. If you repeat your simulation more times you may find additional ensembles. Make sure to start from a different set of initial velocities.
Superimpose the protein along your simulations, you may be able to see the difference.
5ns isn't enough to say anything useful about the equilibrium states, but short of an amazingly unlikely series of coincidences, it's certainly enough to say that there's something peculiar about that pair of runs that doesn't fit with the temperature difference being the only variable in play.
If I were you, I'd be less worried about the apparent difference in RMSD (which as others have already mentioned, could come from numerous sources that your current runs are insufficient to tell you about), and be more worried about the difference in the "randomness" of the RMSD. While the profile of the trajectories obviously shouldn't be identical, if the only variable is the temperature, they should be qualitatively similar, with a larger envelope at higher temperature (though the difference between 310 and 312 envelopes should be insignificant). Seeing that much difference between the profiles suggests that you've got a different set of atoms selected between your runs, or that your initial conditions differ in some fashion.
Your result best epitomizes the difficulties with MD. People who do MD for living consider it very well based science. Many valuable results were obtained using it. You can use it but seek help from more experienced coworkers. In essence what you obtained is nonsense but nonsense that is expected from a weak method with weak methodology. Even if you carefully prepare your structure before the simulation (careful removal of offending conformational errors, minimization, solvent and counterions addition and initial prequilibration) you might still expect problems. My opinions, after using the mothod in many contexts for a long time, are shaped by my background which is a combination of in depth training in physics (PhD) and 30 years of experience in experimental protein work (crystallography and everything else).
MD is just a computer exercise (like school work) and can become science only in the hands of a very experienced user. There is nothing objective nor repeatable in doing MD simulation, therefore on its own right this does not qualify to be science. The lowest common denominator of all science methods are principles of repeatability and independence of the person who is doing the experimentation. To fully qualify as science it also needs to have a certain element of science inference which is predictability.
MD fails every and any of these principles without a very careful user. Nevertheless all the leaders in the field are so strongly attached to marketing their contributions that most of the fundamental problems of MD are swept under the carpet.
Main difficulties are:
1) MD is a classical physics (not quantum) approximation for movements of quantum objects. Everybody using the method must understand the limitations it brings. Car-Parinello method is trying to deal with that but it is very computer hungry.
2) Forcefield weaknesses, simple central harmonic forces with crude angular and torsional representations.
3) Approximations neglecting polarizabilities and quantum evolution of charge (people are working on it) and ionization (harmonic approximation does not allow for ionization.
4) Very crude dynamical treatment of time evolution.
5) As a strongly coupled nonlinear system it is chaotic and by definition not PREDICTABLE (It has been show in several papers to have Lyapunov instability and exponential divergence).
These would be sufficient to have limited trust in widespread use of MD but additional elements compound its difficulties.
6) Almost every protein structure in the PDB has errors. Even PDB Redo does not remedy that. A very careful evaluation of the structure before the runs should be mandatory using tools developed by Dr Richardson but it is not and most people using MD does not even check basic stereo-chemistry. I have seen amazingly bad stereo-chemistry in many results of MD simulations. So placement of all protons in the structures should be mandatory as a tool for spotting the biggest errors. To give you an example. A single incorrect conformer of a medium size residue (such as Leu) can carry as much energy to explode your structure. (Average folding energy ~35kcal, compare to and average conformational error of a single torsional angle is around 9kcal). Of course people developed remedies against it but at the cost of distorting the structures (tighter constraints).
7) Questionable treatment of energy flows during simulations (temperature as well as pressure control). A standard velocity rescaling violates energy conservation only to restore the energy. A standard procedure to remove the movements of center of mass is violating energy conservation only to restore the energy. Standard methods for heat bath produce energy ad nihilo only to restore the energy.
So in summary, assigning any significance to the very large RMS difference between two runs is utterly meaningless. As another advise, you should go very careful through preparatory stages and repeat runs multiple times to really say anything substantive.
Besides, 6Å RMSD is just simply a sign of something bad especially that more than 30% is contained in the very first 200 ps.
Excellent insight into the shortcomings of Classical MD. Recently one of the reviewer asked me to repeat the 30 ns MD simulation 8 time, in order to sample enough. Although, sampling space covered remain a matter of debate in classical MD. One should be extremely cautious while interpreting MD results. Resource intensive QM methods are much better specially for small systems, if one is looking for detailed analysis at the residue level involving bond formation and breaking. Until, we come up with ways to simulation biophysical systems using full QM approach, classical MD could be used to answer the structural-domain level problems.
As correctly indicated by Boguslaw, MD is a computer exercise. Each simulation depends on the starting conformation of the system, but also on the starting conditions, and that includes Temperature, which is simulated through a somewhat random assignment of the starting velocities to the particles in the system, constraint to fit the Maxwell-Boltzmann distribution. From there the system is allowed to evolve as described by one numerical integration algorithm. Unlike computer calculations of smaller problems, numerical errors create stochastic rounding-off loss of information in the simulation and "in theory", after some simulation time, the system behavior becomes statistically uncorrelated to the precedent behavior (the size of such time-window is to be determined for each system a set of conditions). However, occasionally, the system may exhibit unexpected behaviors due to singularities, which are not due to the real system, but to peculiar conditions where the simulation fails to approximate to the real world.
If you start two simulations with exactly the same distribution of velocities in the same machine and operating system, the system will behave the same even if its behavior is anomalous.
If you want to analyze the change in RMSD of your system with the initial state as reference (that is the graph you provided), you must perform a considerable number of simulations and average the results. This is specially important for a temperature difference so small as 2K. Statistically speaking, you are likely to require a large number of simulations and a careful analysis of the results to observe differences in the temperature dependence of any property in such a narrow temperature differential. If this is important for you, make sure to validate your system following additional properties amenable to comparison to available experimental data.