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.
