Will molecular dynamics simulations of proteins ever reach equilibrium

Samuel Genheden, Ulf Ryde. Phys. Chem. Chem.Phys., 2012, 14, 8662-8677

Contributed by Gerald Monard

When one runs molecular dynamics simulations, one always wonders whether the trajectories are “long enough”. The ergodic hypothesis states that with an infinite trajectory (in time), all accessible states should be reached. However, one cannot wait an infinite amount of time. In practice, trajectories are finite. Thus the question: did my molecular dynamics simulation have reach equilibrium?

S. Genheden and U. Ryde have analyzed conformational entropies of five protein and protein-ligand complexes with different methods for simulations up to 500 ns. They show that, for all cases, conformational entropies have not converged and that, therefore, all available states have not been reached.

Conformational entropies have been evaluated from the distribution of the dihedral angles, using different models: dihedral-distribution histogramming (DDH), von Mises approach (VMA), quasi-harmonic analysis (QHA), and normal-mode analysis (NMA). The latter was used in the estimation of ligand-protein binding free energies with the MM/GBSA method.

What is interesting in this article is that, whatever the entropy model used, no simulation converged and that entropies are still increasing at the end of the simulations. To explain this, the authors used a very simple model of protein containing only dihedral angles. In this case, the entropy of the model protein increases logarithmically with the simulation time, and the conformational entropy does not converge until the simulation time is larger than all the normal mode mean lifetimes (from 0.2 ps for 1 kJ/mol barrier to 12 h for 100 kJ/mol barrier). Using such an approach, the authors state that after 500 ns of simulation, most of their systems have sampled less than 70% of the conformational spaces and that the entropy will continue to increase with simulation time. Moreover, in the case of BPTI, after 1 ms of simulation, entropy has still not converged!

So, does this mean that we are all wrong when performing MD simulations? Will MD simulations ever reach equilibrium? In the strict sense, the answer is no: MD simulations do not converge, at least when regarding the entropy term. But, do we really need full convergence? As the authors say, reaching all available states for a protein would mean to fold and unfold the protein several times during the MD simulation. Do we (always) need this? In most cases, no.

Another interesting point is the MM/GBSA calculations. Using the NMA model to compute the entropy, the ligand-protein binding free energies converge in the 500 ns simulation time, showing a variation of less than a few kJ/mol. However, this is mainly due to error cancellation in the sum that defines the binding free energy in the MM/GBSA approach. In addition, the NMA approach produces much more stable estimates of the entropy than the other three methods (DDH, VMA, and QHA). This is due to the NMA approach that assumes a single conformation for the protein and estimates the entropy form the vibrational spectrum in this conformation. As long as the protein does not drastically change its conformation and its vibrational spectrum, the NMA entropy does not vary much (but then one could argue that this means that the NMA approach may not be very accurate).

In conclusion, I would say that this article shows that it is hopeless to see MD simulations of proteins to ever reach equilibrium (at least using a single trajectory). But that does not mean we should stop running MD simulations. We have to be (very) careful about the quantities that we are looking at. Different properties show varying sensitivity to this lack of strict equilibration.