Generalized Ensemble Molecular Dynamics Methods

Classical molecular dynamics is a computer simulation technique that is in widespread use across many areas of science, from physics and chemistry to materials, biology, and medicine. The method continues to attract criticism due its oft-reported lack of reproducibility which is in part due to a failure to submit it to reliable uncertainty quantification (UQ). Here we show that the uncertainty arises from a combination of (i) the input parameters and (ii) the intrinsic stochasticity of the method controlled by the random seeds. To illustrate the situation, we make a systematic UQ analysis of a widely used molecular dynamics code (NAMD), applied to estimate binding free energy of a ligand-bound to a protein. In particular, we replace the usually fixed input parameters with random variables, systematically distributed about their mean values, and study the resulting distribution of the simulation output. We also perform a sensitivity analysis, which reveals that, out of a total of 175 parameters, just six dominate the variance in the code output. Furthermore, we show that binding energy calculations dampen the input uncertainty, in the sense that the variation around the mean output free energy is less than the variation around the mean of the assumed input distributions, if the output is ensemble-averaged over the random seeds. Without such ensemble averaging, the predicted free energy is five times more uncertain. The distribution of the predicted properties is thus strongly dependent upon the random seed. Owing to this substantial uncertainty, robust statistical measures of uncertainty in molecular dynamics simulation require the use of ensembles in all contexts.

Journal of Chemical Theory and Computation, 2012

Serial generalized ensemble simulations, such as simulated tempering, enhance phase space sampling through non-Boltzmann weighting protocols. The most critical aspect of these methods with respect to the popular replica exchange schemes is the difficulty in determining the weight factors which enter the criterion for accepting replica transitions between different ensembles. Recently, a method, called BAR-SGE, was proposed for estimating optimal weight factors by resorting to a self-consistent procedure applied during the simulation (J. Chem. Theory Comput. 2010Comput. , 6, 1935Comput. −1950. Calculations on model systems have shown that BAR-SGE outperforms other approaches proposed for determining optimal weights in serial generalized ensemble simulations. However, extensive tests on real systems and on convergence features with respect to the replica exchange method are lacking. Here, we report on a thorough analysis of BAR-SGE by performing molecular dynamics simulations of a solvated alanine dipeptide, a system often used as a benchmark to test new computational methodologies, and comparing results to the replica exchange method. To this aim, we have supplemented the ORAC program, a FORTRAN suite for molecular dynamics simulations (J. Comput. Chem. 2010Chem. , 31, 1106Chem. −1116, with several variants of the BAR-SGE technique. An illustration of the specific BAR-SGE algorithms implemented in the ORAC program is also provided.

Computer Physics Communications, 2006

Equations of motion based on an atomic group scaling scheme are described for a molecular system with bond constraints. The NPT ensemble extended system method is employed along with a numerical integration scheme using an operator technique. For parallelization of the integration scheme, a domain decomposition scheme is employed based on a group of atoms which share common constraints. This decomposition scheme fits well into the integration scheme and involves no extra inter-processor communication during the SHAKE/RATTLE procedures. An example is given for a solvated protein system containing a total of 23 558 atoms on 64 processors.

Advances and Applications in Bioinformatics and Chemistry, 2015

Molecular dynamics simulations have evolved into a mature technique that can be used effectively to understand macromolecular structure-to-function relationships. Present simulation times are close to biologically relevant ones. Information gathered about the dynamic properties of macromolecules is rich enough to shift the usual paradigm of structural bioinformatics from studying single structures to analyze conformational ensembles. Here, we describe the foundations of molecular dynamics and the improvements made in the direction of getting such ensemble. Specific application of the technique to three main issues (allosteric regulation, docking, and structure refinement) is discussed.

Biophysical Journal, 2015

We introduce an enhanced-sampling method for molecular dynamics (MD) simulations referred to as ensemblebiased metadynamics (EBMetaD). The method biases a conventional MD simulation to sample a molecular ensemble that is consistent with one or more probability distributions known a priori, e.g., experimental intramolecular distance distributions obtained by double electron-electron resonance or other spectroscopic techniques. To this end, EBMetaD adds an adaptive biasing potential throughout the simulation that discourages sampling of configurations inconsistent with the target probability distributions. The bias introduced is the minimum necessary to fulfill the target distributions, i.e., EBMetaD satisfies the maximum-entropy principle. Unlike other methods, EBMetaD does not require multiple simulation replicas or the introduction of Lagrange multipliers, and is therefore computationally efficient and straightforward in practice. We demonstrate the performance and accuracy of the method for a model system as well as for spin-labeled T4 lysozyme in explicit water, and show how EBMetaD reproduces three double electron-electron resonance distance distributions concurrently within a few tens of nanoseconds of simulation time. EBMetaD is integrated in the open-source PLUMED plug-in (www.plumed-code.org), and can be therefore readily used with multiple MD engines.

Computer Physics Communications, 2002

The most efficient MC weights for the calculation of physical, canonical expectation values are not necessarily those of the canonical ensemble. The use of suitably generalized ensembles can lead to a much faster convergence of the simulation. Although not realized by nature, these ensembles can be implemented on computers. In recent years generalized ensembles have in particular been studied for the simulation of complex systems. For these systems it is typical that conflicting constraints lead to free energy barriers, which fragment the configuration space. Examples of major interest are spin glasses and proteins. In my overview I first comment on the strengths and weaknesses of a few major approaches, multicanonical simulations, transition variable methods, and parallel tempering. Subsequently, two applications are presented: a new analysis of the Parisi overlap distribution for the 3d Edwards-Anderson Ising spin glass and the helix-coil transition of amino-acid homo-oligomers.

In this chapter a summary is given of the key ingredients necessary to carry out a molecular dynamics simulation, with particular emphasis on macromolecular systems. We discuss the form of the intermolecular potential for molecules composed of atoms, and of non-spherical sub-units, giving examples of how to compute the forces and torques. We also describe some of the MD algorithms in current use. Finally, we briefly refer to the factors that influence the size of systems, and length of runs, that are needed to calculate statistical properties.

2001

Condensed Matter Physics, 2007

The article reports an attempt to study the protein folding problem by generalized-ensemble Monte Carlo simulations with reference interaction site model theory. Generalized-ensemble algorithms greatly enhance the configurational space sampling in computer simulations. The reference interaction site model theory treats solvent effects with solvent molecular shape and estimate solvation free energy around proteins. We have developed simulation algorithms which combine generalized-ensemble algorithms and one-dimensional reference interaction site model theory. This treatment can also use a simulation with three-dimensional reference interaction site model theory. In this review, we describe the methods and present the results of these simulations for a peptide.

Physics Today, 1997