|Home | About | Journals | Submit | Contact Us | Français|
The rugged energy landscape of biomolecules together with shortcomings of traditional molecular dynamics (MD) simulations require specialized methods for capturing large-scale, long-time configurational changes along with chemical dynamics behavior. In this report, MD-based methods for biomolecules are surveyed, involving modification of the potential, simulation protocol, or algorithm as well as global reformulations. While many of these methods are successful at probing the thermally accessible configuration space at the expense of altered kinetics, more sophisticated approaches like transition path sampling or Markov chain models are required to obtain mechanistic information, reaction pathways, and/or reaction rates. Divide-and-conquer methods for sampling and for piecing together reaction rate information are especially suitable for readily available computer cluster networks. Successful applications to biomolecules remain a challenge.
It has been said that, while the 19th century belonged to chemistry and the 20th century to physics, the 21st century is and will continue to be dominated by ground-breaking science in biology. However, it is clear that biology without the tools of mathematics, physics, chemistry - not to speak of computer science and engineering innovations - would be severely limited. In fact, computer modeling and simulation offer a modern tool for simulating conformational events in a variety of molecular systems and subsequently extracting related mechanistic, thermodynamic, and kinetic information.
One spectacular example of the modern marriage between technology and biology involves the theory and practice of molecular dynamics (MD). MD, with variations and extensions, has become a universal tool for probing relationships among biomolecular structure, flexibility, and function. Since successful applications were reported in the 1970s in protein dynamics , MD has become a universal tool, ‘as if it were the differential calculus’ .
Of course, MD is simply statistical mechanics by numbers, or Pierre-Simon de Laplace's vision of using physical laws to simulate time-dependent events , on modern supercomputers. And impressive progress in the development of biomolecular force fields, coupled to spectacular computer technology advances, has now made it possible to overcome the difficulty noted by Paul Dirac of solving the complex equations of motion .
Now, sophisticated dynamics programs, like NAMD, Desmond, or GROMACS, adapted to parallel and massively parallel computer architectures, have made simulations of biomolecular systems of one microsecond and longer feasible in several weeks of computing [5-8]. However, the hunger for longer simulations as an end by itself may not be justified given that such simulations lead to as many questions as answers, and underscore the notion that computer power alone is not likely to solve ‘folding’ problems. Indeed, many approximations (for example, in force field functional form as well as numerical values) and sampling limitations still persist. And besides potential function inaccuracies, sampling the vast conformational space remains a grand challenge because standard MD methods cannot sample the vast thermally accessible configuration space of biomolecules to capture large-scale rearrangements, including reaction pathways, mechanisms, and rates. Rather, MD applications are valuable in probing local system fluctuations (Figure 1, bottom left) for different single-residue variants of DNA polymerase λ .
An enormous range of methods has been developed and applied to molecular systems in an attempt to enhance coverage of the thermally accessible conformational space. The efficacy/ability of each approach depends on the goals of the computation and the available computing resource. No review or presentation could possibly list them all. Instead, I will present some general concepts and highlight selected interesting applications in this brief perspective. See [10,11] for recent reviews.
In the first review of this series , Monte Carlo (MC) approaches (simulated annealing, hybrid MC, and parallel tempering), harmonic approximations (methods based on spectral decompositions like principal component analysis [PCA] and essential dynamics), and coarse-graining formulations (multiscale models for protein and supramolecular systems) were surveyed; methods like PCA provide information on dominant functional motions, such as subdomain movements in a polymerase β/DNA complex  (Figure 1, top left), and MC simulations are especially effective for coarse-grained models, such as chromatin  (Figure 1, bottom right).
Here, methods based on MD are considered. These can be as simple as modifying the potential or the basic protocol or can involve more complex modifications of the algorithms, including global reformulations. Methods in all categories (for example, replica exchange MD, or REMD) are generally effective for gaining insights into broader structural and thermodynamic questions than possible in typical MD, although a rigorous assessment of their performance is generally lacking. However, increasingly sophisticated methods - grouped as global reformulations - are needed for describing mechanisms, free energy profiles, and/or detailed kinetics. These include transition path sampling (TPS), forward flux simulation (FFS), and Markov chain models, all of which pay special attention to preserving microscopic reversibility and thermal distributions of the systems.
The simplest and crudest approach for enhancing the sampling is to manipulate the energy function used as a basis for MD by using restrained potentials, as in targeted MD (TMD) or umbrella sampling (US) by using various experimentally based biases or guides.
In TMD, an artificial restraint potential term is added like a Lagrange multiplier with a parameter λ (0 at the initial state and 1 at the target state) so as to force the system to evolve toward a certain state (for example, closed and open state of an enzyme) in a specified number of steps. However, the actual trajectory obtained using this restraint leads to unphysical dynamics since incorrect transition states are often achieved. Nonetheless, TMD can be used to explore space in a preliminary manner.
In 2002, Kong et al.  reported a TMD application to the opening of an Escherichia coli membrane channel whose open and closed states were modeled. The simulations suggested a four-stage process in the opening, unlike simpler paths suggested previously. An F1000 commentator writing on this work considered TMD to be ‘one of the most promising methods to suggest possible pathways for large-scale conformational changes’ .
Similar in spirit is the simplified multiple-basin Hamiltonian funnel-based potential model . Developed for very large molecular complexes with known endpoint structures, such a model can be used, as in TMD, to simulate very large-scale motions, such as transitions between ligand-bound and unbound states.
However, over the past few years, we have come a long way. Now, we have more mathematically rigorous methods to sweep conformational space, deduce mechanisms, and compute reaction rates as well as free energies by a variety of innovative methods. Still, because the applications of those more sophisticated methods are far from routine for biomolecules, there is room for the simpler, albeit more artificial, approaches.
US focuses on specific regions of phase space by using a restraining potential. In a recent application, Mills and Andricioaei  developed a ‘guided US’ for biomolecules in which experimental restraints [for example, end-to-end distance determined by fluorescent resonance energy transfer (FRET) experiments] are used to formulate potential restraints. They showed dramatic improvement in efficiency (that is, faster convergence and better final distances for the guided versus unguided US).
Another interesting example of using experimental restraints was reported by Paci et al. . They constructed transition state ensembles (TSEs) of simple two-state (fast) folding proteins by restraining the conformational energy to states that satisfy experimentally known contact values of (ratio of the number of native contacts present in the protein transition state to that in the native state). Interestingly, their computed TSEs for three fast-folding proteins revealed that members of these TSEs have approximately 70% native-like folds even though rms distances, relative to the folded state, can be large.
Besides various constraints or restraints, the energy can be alternatively modified by using other types of added energy terms. For example, McCammon's group  recently used accelerated MD in which a term is simply added to the potential energy ‘to boost’ the potential to allow escape from minima. For numerical considerations, the addition is implemented so that the gradient is continuous. Moreover, to be effective, the boost can be applied to one term only (for example, torsional terms) so as not to destroy liquid water structure around the biomolecule.
Such a boost approach is based on many other proposals of adding biasing potentials to overcome barrier crossing reported previously. These include Grubmüller's conformation flooding in which the bias is determined based on a coarse-grained simulation for the conformational space density , Voter's hyper-MD in which the biasing potential is constructed based on the smallest eigenvalue of the Hessian matrix , the simple local boost method based on the total potential energy , and a bond boost method based on bond-length deviations .
The MD simulation protocol rather than the energy can also be manipulated. In this subclass, we have the popular approach LES (locally enhanced sampling), in which several configurations are generated in a single energy evaluation , and REMD , based on parallel tempering MC , which itself is based on Metropolis-coupled Markov chain MC.
The idea in REMD is to simulate multiple, non-interacting copies (replicas) of identical systems at different temperatures. Periodically, the configurations of the replicas are exchanged (and thus the thermodynamic, or temperature, states are exchanged) with a transition probability P that maintains each temperature's equilibrium ensemble distribution in the canonical ensemble so as to retain those systems that are making better progress:
where βi = 1/(kBTi) and Ui are the internal energies of states i (new and old). These temperature-ensemble exchanges are meant to accelerate crossings of energy barriers. Though obeying detailed balance for an extended ensemble of the canonical states, the state-exchange probability destroys real kinetic properties. However, REMD has been successfully used to simulate folding/unfolding equilibria of biomolecules.
REMD has been increasing in popularity for applications ranging from small peptides to complex biological systems. However, its success has largely been empirical. For example, among the numerous applications to biomolecules reported, the folding of a small solvated RNA hairpin from an extended state  revealed a diverse ensemble of conformations with differing stacking and base pair arrangements; applications to solvated protein A  also revealed the folding as the temperature decreases from about 600 to 280 K (see posted animations ). REMD villin folding simulations of 200 nanoseconds in length consistently folded within 1.78 Å of the native state .
More recently, practitioners of REMD have emphasized the need to formulate careful configurational swapping protocols (temperature ladders and exchange/acceptance ratios) and other ways to increase the conformational sampling efficiency (see [32,33], for example). Other practical considerations have been discussed. For example, because efficient sampling in REMD requires a large number of replicas to enhance the swapping rate, the user must have a large number of concurrent processors available for the REMD simulation; this may not be the case for the average user, and the REMD variant distributed replica sampling  may be preferable. Kamberaj and van der Vaart  describe an efficient multiple scaling REMD in which a Tsallis biasing potential increases sampling efficiency because barrier crossings are enhanced so the number of replicas can be reduced. For a note on the likely computational advantage of REMD for systems with relatively high energy barriers, see .
Cooke and Schmidler  recently presented a thorough analysis of REMD, explaining that the formulation of REMD by extension of parallel tempering MC to MD has introduced sampling and ergodicity problems stemming from the failure of the underlying constant-temperature (canonical) MD integrators to preserve certain variants. This is because, while the MC analog is based on Markov chains, REMD algorithms cannot use the symplectic leap-frog integrator popular for microcanonical (constant energy) MD and resort to isothermal (constant temperature) integrators like the Berendsen heat-bath algorithm or the Nosé-Hoover thermostat method, which are not rigorously ergodic; this can affect the dynamics of even small systems. As a remedy, combining REMD with hybrid MC to ensure ergodicity is suggested , thereby returning to the original advantages of parallel tempering MC. Recently developed entropy-preserving constant-temperature integrators can also solve the ergodicity problem in REMD .
Many other REMD combinations have been reported. Chebaro et al.  combined REMD with a coarse-grained implicit solvent model of proteins (OPEP) to determine equilibrium ensembles of proteins. Bolhuis  combined REMD with transition interface sampling (TIS) to compute the free energy profile and rate constants when high barriers are involved. These combinations emphasize how different sampling methods can be combined in clever ways to suit the problem and enhance sampling further.
Besides altering the energy function as in TMD and simulation protocol as in REMD, the MD algorithm itself can be altered.
Many integrators have been developed to address stability and resonance limitations of MD integrators due to the high-frequency vibrational modes and the intricate coupling among the vibrational modes of a biomolecule spanning a wide spectrum. See , for example, for a discussion of these issues and various solutions; most successful approaches to both stability and resonance problems consist of multiple-timestep integrators combined with stochastic dynamics.
Langevin dynamics and Brownian dynamics (BD), in which solvent collisions and thermal fluctuations are incorporated in an average sense, have long been used as a way to allow larger timesteps and hence larger timespans in dynamics simulations. For stochastic dynamics methods, it is also easier to prove ergodicity.
An interesting approach in which the algorithm is modified was recently reported by Sweet et al. . Their approach, called ‘NML’, is an extension of a LIN  developed as a way to increase the timestep in standard MD integration. The idea in LIN was to use implicit integration for the low-frequency modes and normal-mode analysis (NMA) for the high-frequency modes; Langevin, rather than Newtonian, dynamics was also used to dampen resonance effects, which limit the timestep to around 3.3 femtoseconds when all light-atom motions are considered . Sweet et al.  replace NMA by the less costly BD for the high-frequency modes and therefore keep the fast oscillations around their equilibrium values; the low-frequency modes are propagated by Langevin dynamics as in LIN.
An interesting recent algorithm for sampling the canonical distribution is based on the Nosé-Hoover formulation in the context of Markov chain MC ; this solution to the lack of ergodicity in the context of Nosé-Hoover integrators is similar to the idea of Cooke and Schmidler  for REMD. Systematic discretization errors are also eliminated by this hybrid MC approach for canonical sampling.
In this altered MD algorithm category, a noteworthy and unusual example is interactive MD (IMD) developed by the Schulten group to probe mechanisms of biological reactions using computer graphics and simulation (VMD + NAMD) [44,45]. The user steers the system to the desired state by ‘feeling’ the potential force and applying any desired magnitude (to ‘pull’ the system). IMD is particularly appropriate for studies of induced-fit mechanisms, in which substrate/solvent/protein interplay are crucial for selectivity and biological function. It requires, however, a sophisticated device on the user's desktop.
Note that MD in internal-coordinate space has been used sporadically with the rationale that the fewer degrees of freedom compared with Cartesian coordinates will allow longer integration timesteps and hence greater sampling. Indeed, peptide folding and refinement with dihedral-angle MD demonstrated a computational advantage of several orders of magnitude compared to Cartesian analogs . However, internal-coordinate MD is not generally used, possibly due to the significant effect of vibrational coupling on biomolecular dynamics and the computational difficulties involved in transforming internal coordinates to/from Cartesian coordinates; it is nonetheless worth more consideration in the context of simplified models, such as macroscopic united-residue models of proteins and nucleic acids.
Many innovative methods have been developed by mathematicians, physicists, and chemists in this category of algorithmic reformulations. (Note: citations below often refer to a recent work rather than the original paper, which can be found in the recent article.)
Examples of methods that aim to explore the free energy landscape include: (a) canonical adiabatic free energy sampling (CAFES) , which propagates dynamics of decoupled solute and solvent systems so that the former follows adiabatic dynamics of the free energy surface created by the latter at increased temperatures; (b) reference potential spatial warping algorithm (REPSWA) , which introduces a variable transformation in the classical partition function that increases attraction to basins; (c) metadynamics , which involves exploring the free energy surface by following non-Markovian dynamics determined by rewriting the equations of motion in terms of a few collective variables so that key regions of space are identified as simulation time increases; and (d) a sweep method that combines temperature-accelerated MD (TAMD) with a free energy reconstruction from the mean force using radial basis functions via optimization .
Methods that deduce mechanisms include: (a) TPS by Chandler and colleagues, which follows an MC protocol in MD-trajectory space to locate key transition states (as recently reviewed in ); (b) the nudged elastic band (NEB) method  that optimizes minimum energy pathways; (c) the max-flux approach , which constructs variationally optimized reaction pathways based on the max flux for diffusive dynamics; and (d) the string method , which prunes MD trajectories to retain reactive segments.
Methods that compute reaction rates include (a) FFSs , which use interfaces to partition phase space; (b) ‘milestoning’ , which coarse-grains temporal and spatial descriptors of the system by sequential transitions; and (c) the popular Markov state models (MSMs) , which define transition networks to describe biomolecular kinetics and thermodynamics. Essentially, all of these methods compute rates by various divide-and-conquer strategies.
TPS developed by Chandler, Bolhuis, Dellago, and Geissler has been popular and successful in various biomolecular application contexts . The idea is to run an MC in trajectory space in a way to create reactive trajectories that preserve microscopic reversibility and local thermal distributions. Full kinetic information can be extracted by clever extensions (such as efficient US, BOLAS , and network models). For example, an application of TPS to DNA polymerase β deduced the mechanism, free energy profile, and rate for pol β's closing conformational change : a complex landscape was revealed where sequential, subtle, side-chain rearrangements lead the enzyme from open to the closed state and where Arg258 is a ‘gatekeeper’ for the reaction (Figure 1, top right); the overall computed rate of 10 per second corresponds to the 27kBT barrier and agrees well with the experimental value of 3 to 10 per second for the overall reaction. The application of BOLAS , instead of US, was crucial for achieving acceptable margins of error; unlike US, BOLAS avoids large biases (which in turn lead to large error bars) by sampling each window by many short MD trajectories.
The same TPS protocol applied to the mismatch (G:A instead of G:C) suggested that the higher free energy barrier for the mismatch comes from the instability of the closed mismatched state compared with the matched base pair system .
Rogal and Bolhuis  recently extended TPS to connect several intermediate states in phase space.
The idea in the popular MSMs  is to compute reaction rates from different interfaces. Combinations such as by Pan and Roux , which employ the single-sweep method of Vanden-Eijnden and colleagues  with MSM, can be successful for proteins because phase space is surveyed to map dynamically important regions, from which free simulations are initiated, and then the transition matrix is constructed by piecing the information together.
The FFS approach  similarly enhances phase space sampling of rare events in stochastic non-equilibrium systems in which the phase space distribution is not known a priori. The phase space is partitioned, and an adaptive procedure is used to find kinetic ‘bottleneck regions’ by estimating rate constants associated with reaching subsequent interfaces. Then, FFS concentrates on sampling those bottleneck regions only. FFS is appropriate for discretely defined surfaces, like protein folding on a lattice; a recent application located the most probable transition state isosurfaces for a 48-unit protein .
Finally, the similar idea in milestoning  is to recover reaction kinetics by integral equations and global path optimization strategies. Milestoning is compared with MSM, FFS, and TIS in , where advantages of milestoning compared with the other methods are suggested in terms of accuracy, efficiency, and parallel-machine implementations.
Numerous other innovative mathematical methods use various constructs from physics, mathematics, and engineering such as domain decomposition, clustering techniques, or Markov models to sample and/or survey conformational space and dynamics. For example, a robust Perron cluster analysis propagates a system's essential dynamics with preservation of Markov properties , and a hierarchical clustering approach by a Markov model helps identify groups of configurations with common features and transitions between them .
Certainly, long MD simulations are being reported more regularly now compared with several years ago as a result of efficient and parallel simulation packages like NAMD, Desmond, GROMACS, and others [5-8], and new computer systems like Anton, hard-wired for long MD simulations . For example, a 1-microsecond atomic-level simulation of the voltage-modulated potassium channel Kv1.2 (120,000 atoms) recently described the opening/closing mechanism involved . Thus, many problems thought to be intractable several decades ago are now being solved. But faster computer technology alone is insufficient to address challenging biomolecular problems. Algorithmic developments for enhanced sampling such as described here, especially those divide-and-conquer methods suitable for loosely coupled processor architectures that are more readily available to the average user, are needed to complement long atomistic simulations.
While simple potential modifications like TMD can provide conformational insights and REMD approaches can enhance sampling, success has been empirical rather than rigorous; because the actual kinetics of the systems are altered, caution is warranted in biological interpretations. However, enhanced sampling protocols can be useful in specific contexts, as recently demonstrated by using accelerated MD  to reproduce residual dipolar coupling measurements concerning the slow modes from nuclear magnetic resonance data on a domain of the protein GB3 . Furthermore, the studies pinpointed the slow motions and led to a model of the potential energy landscape of the protein.
REMD has recently received some scrutiny concerning effective ensemble exchange protocols, general computational efficiency, and ergodicity questions, and this has led to variants with stochastic elements - motivated by REMD's origin in MC parallel tempering methods - that can improve sampling (for example, ). Because REMD applications require concurrent processors (one per replica), the technique may not always be practical, especially for large atomistic systems. Variants like the distributed replica sampling  can be more viable by performing stochastic moves of independent replicas instead of pairwise exchanges of replicas.
TPS, MSMs, and approaches that aim to compute reaction rates have been more rigorously grounded in theory, and successful biomolecular applications have been reported. Still, their application to biomolecules in general remains far from routine. In this goal, various combinations of methods that deduce mechanisms and compute reaction rates by divide-and-conquer approaches will undoubtedly be effective in today's readily available distributed computing resources of cluster networks, especially by enhancing them with coarse-grained models and MC elements.
Ultimately, the most successful applications reflect a combination of novel and rigorous mathematical ideas with physical intuition. And such applications can be successful not only by reproducing experimental data concerning slow motions but also in helping suggest mechanistic and energetic details. Clearly, the gap between experimental and theoretical timeframes is steadily narrowing, leading to renewed interest in MD-based modeling by theorists and experimentalists alike.
Support from the National Science Foundation, the National Institutes of Health, the American Chemical Society, and Philip Morris International is gratefully acknowledged. The author thanks Meredith Foley, Shereef Elmetwaly, Hin Hark Gan, and Eric Vanden-Eijnden for valuable assistance and comments.
The electronic version of this article is the complete one and can be found at: http://F1000.com/Reports/Biology/content/1/51
The author declares that she has no competing interests.