Enhanced sampling techniques in molecular dynamics simulations of biological systems

Abstract

Background: Molecular dynamics has emerged as an important research methodology covering systems to the level of millions of atoms. However, insufficient sampling often limits its application. The limitation is due to rough energy landscapes, with many local minima separated by high-energy barriers, which govern the biomolecular motion.

Scope of review: In the past few decades methods have been developed that address the sampling problem, such as replica-exchange molecular dynamics, metadynamics and simulated annealing. Here we present an overview over theses sampling methods in an attempt to shed light on which should be selected depending on the type of system property studied.

Major conclusions: Enhanced sampling methods have been employed for a broad range of biological systems and the choice of a suitable method is connected to biological and physical characteristics of the system, in particular system size. While metadynamics and replica-exchange molecular dynamics are the most adopted sampling methods to study biomolecular dynamics, simulated annealing is well suited to characterize very flexible systems. The use of annealing methods for a long time was restricted to simulation of small proteins; however, a variant of the method, generalized simulated annealing, can be employed at a relatively low computational cost to large macromolecular complexes.

General significance: Molecular dynamics trajectories frequently do not reach all relevant conformational substates, for example those connected with biological function, a problem that can be addressed by employing enhanced sampling algorithms. This article is part of a Special Issue entitled Recent developments of molecular dynamics.

Keywords: Cellulosome; Enhanced sampling; Generalized simulated annealing; Metadynamics; Molecular dynamics; Replica-exchange molecular dynamics.

Figure 1

Illustration of the replica exchange molecular dynamics (REMD) method. A set of noninteracting replicas runs at different values of an exchange variable, usually temperature (T-REMD). At specific intervals, replicas at neighboring values for the exchange variable are swapped based on a Monte Carlo acceptance criterion. In an efficient run, all trajectories will experience changing of the exchange variable value. At each value for the exchange variable, the trajectories will be discontinuous, but follow a proper Boltzmann distribution for the specific value being exchanged.

Figure 2

Illustration of the metadynamics method. Described as “filling the free energy wells with computational sand”, the metadynamics method allows the search inside each energy well avoiding an oversampling of the same conformations. When the system reaches a point where the energy is higher than a barrier separating two minima the system goes to a state of lower energy in the new minimum, again searching many possible conformations.

Figure 3

Illustration of the simulated annealing method. Starting from a non-minimized structure and with a high temperature, the structure can randomly assume any conformation. A jump from conformation (1) to (2) is accepted if the simulated annealing acceptance criterion is achieved. Since the energy of conformation (2) is lower than the energy of the conformation (1), it is automatically accepted. Because the temperature decreases with each step, the modifications have a tendency to be progressively smaller. A conformation with a higher temperature can also be accepted as in the jump from (4) to (5). At the end the jumps are much smaller and there is a tendency to find a local minimum. Due to the high conformational energy at the beginning of the simulation, a large number of independent simulations should be capable of sampling all the possible stable conformations of the studied molecule.

Figure 4

Generalized simulated annealing study of a cellulosome fragment. A: C. thermocellum CipA scaffolding fragment in complex with the SdbA type II cohesin module (CohII). CipA fragment components: ninth type I cohesin (CohI9), small flexible linker, X-module (Xmod) and type II dockerin (DocII). The cartoon representation depicts the shape of the protein backbone while the flexible linking region is shown in atomic representation. Only the single covalent bonds from the eight residues in this region were allowed to change throughout the simulations, both for backbone and side chains. The nearby residues THR150 and GLY157, both inclusive, defined the flexible region. B: Center of Mass of CohI₉ module in all visited conformations, shown as a blue sphere. The center of mass determination was done using only the atoms in the CohI9 module, excluding the flexible X module. Both calculation and figure were created using VMD [90]. The analysis shows that the CohI9 module visited mainly two opposite positions. The arched-band positioning of the centers of mass show how much the module can fluctuate around each position. In gray, one can observe the “native” position of the CohI₉ module and the linker in the final minimized and equilibrated structure, initially taken from the crystal structure deposited in the PDB. The structure in gray was used to initiate all simulations and is not part of the results, it is, however, in the center of the region most visited by the independent simulations. In green, a second structure is shown that represents in an exemplary fashion the “alternate ”position of the CohI₉ module and linker visiting the second conformational cluster. C: Detail of the two conformations of CohI₉ observed in B, and the arched bands showing the positions of centers of mass of final conformations from all 51,200 simulations.

Figure 5

Potential energy distribution of the cellulosome fragment after 10,000 steps of GSA simulations. The 51,200 different conformations resulted in 51,200 different energies that are coupled according to the two main clusters of conformations. The “native” cluster (represented in red) occurs in 34936 structures, or 68% of the cases. The “alternate” cluster (represented in blue) occurs in 16264 structures, or 32% of the cases. All structures were assigned to a cluster using the software package Weka [91]. The software package R [92] was used to calculate a normalized histogram of the potential energy of the structures in each cluster. A Gaussian distribution fits each of the distributions and the difference in average energy between the two conformations is about 1.2 kcal/mol.