Will molecular dynamics simulations of proteins ever reach equilibrium?

Abstract

We show that conformational entropies calculated for five proteins and protein–ligand complexes with dihedral-distribution histogramming, the von Mises approach, or quasi-harmonic analysis do not converge to any useful precision even if molecular dynamics (MD) simulations of 380–500 ns length are employed (the uncertainty is 12–89 kJ mol⁻¹). To explain this, we suggest a simple protein model involving dihedrals with effective barriers forming a uniform distribution and show that for such a model, the entropy increases logarithmically with time until all significantly populated dihedral states have been sampled, in agreement with the simulations (during the simulations, 52–70% of the available dihedral phase space has been visited). This is also confirmed by the analysis of the trajectories of a 1 ms simulation of bovine pancreatic trypsin inhibitor (31 kJ mol⁻¹ difference in the entropy between the first and second part of the simulation). Strictly speaking, this means that it is practically impossible to equilibrate MD simulations of proteins. We discuss the implications of such a lack of strict equilibration of protein MD simulations and show that ligand-binding free energies estimated with the MM/GBSA method (molecular mechanics with generalised Born and surface-area solvation) vary by 3–15 kJ mol⁻¹ during a 500 ns simulation (the higher estimate is caused by rare conformational changes), although they involve a questionable but well-converged normal-mode entropy estimate, whereas free energies estimated by free-energy perturbation vary by less than 0.6 kJ mol⁻¹ for the same simulation.