Molecular dynamics simulations: advances and applications

Abstract

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.

Keywords:

• molecular dynamics
• allostery
• docking
• conformational ensembles
• structure prediction
• refinement

Introduction

The study of the macromolecular structure is a key point in the understanding of biology. Biological function is based on molecular interactions, and these are a consequence of macromolecular structures. Since initial structure determinations in the 50s, both in the protein and in the nucleic acid worlds, the increase in the knowledge of how macromolecular structures are built has been continuous. At present, protein data bank (PDB)¹ holds more than 110,000 entries, including more than 100,000 proteins, 2,800 nucleic acids, alone or forming complexes, and approximately 20,000 small molecules complexed to macromolecules. Molecular recognition rules as defined by such structural knowledge powers the understanding of basic biological phenomena, like enzyme mechanisms and regulation, transport across membranes, the building of large structures like ribosomes, or viral capsids, or how DNA is read and transcription is controlled. The study and prediction of protein–protein interaction networks is one of the growing fields in modern systems biology. On a more practical note, protein three-dimensional (3D) structures are the basis for structure-based drug design. The simple visual analysis of 3D structures of protein or nucleic acids, as obtained from the experiments, has driven large number of successful studies in biochemistry. However, despite their enormous utility, structures stored at the PDB provide only a partial view of 3D structure. Both protein and nucleic acids are flexible entities, and dynamics can play a key role in their functionality. Proteins undergo significant conformational changes while performing their function. As a rule, any complex made by any protein implies some structural rearrangement. This can be easily checked just by comparing a series of PDB entries that just differ in a small ligand bound to a given protein. Figure 1 shows a superimposition of experimental acetylcholinesterase structures. There are no changes in the overall fold, just small rearrangements in the structure; however, these differences are large enough to fool ligand-docking algorithms. Larger conformational changes are also present in the known protein structures.^2–6 Conformational changes are a common part of an enzymes’ catalytic cycle.⁷ For instance, loop or domain closures contribute to isolate the active site from solvent and, in so doing, alter the chemical environment around substrates, or trigger the catalytic event by bringing essential partners together.^8–10 Also, allostery is the most common enzyme regulation strategy. Allosteric regulation is entirely based on the possibility of a given protein to coexist in two or more conformations of comparable stability. Binding to ligands (allosteric regulators), or simply protein concentration, or crowding, may switch stabilities among conformations and trigger the shape transition. Additionally, some features of protein function can be understood only when dynamic properties are taken into account. For instance, diffusion of small substrates through heme-dependent enzyme molecules requires the transient appearance of channels in the protein structure.^11–18 Also, cavities have transient phenomena that in some cases can only be revealed or analyzed following its dynamics.^19,20

Figure 1 Structure variability within a protein family.

Notes: Structures of acetylcholinesterase (1acg, 1ax9, 1dx6, 1hbj, 1qon, 1vot, and 2ace) crystallized with different active-site ligands.

In the case of nucleic acids, conformational changes are even more complex. Standard B-DNA has a relatively simple structure in comparison with protein or complex RNAs; however, it is an extremely plastic molecule that undergoes large conformational changes to adapt to its interaction partners. Binding of transcription factors to DNA, for example, is not only dependent on DNA sequence recognition, but also a direct consequence of the ability of the DNA molecule to adapt to the protein surface.^21,22

The traditional approach to understand conformation influence on macromolecular function is to cumulate experimental structures covering the conformational space. This has led to the generation of crystal structures for macromolecules in several environments, or macromolecules complexed with different molecules, and contributes to the enormous redundancy seen in the PDB. Examples of this approach are the 87 structures for CK2 homologues (Figure 2A), where a common fold is maintained, and different degrees of conformational variation are clearly visible, especially in loop regions. A single experiment could generate conformational ensembles as those taken from nuclear magnetic resonance experiments (Figure 2B). In the latter case, the source of the variability found is rather a consequence of the lack of experimental data in some specific regions of the structure. Indeed, the study of PDB as source for molecular flexibility has been exploited in some extent.^18,23–26 Taking such “experimental ensembles” a partial view of the macromolecule flexibility can be obtained, although PDB composition is necessarily biased. Theoretical techniques appear as the most convenient way to obtain a picture of macromolecular dynamic properties. Recent advances in the performance of simulation algorithms, including specific strategies to increase the conformational sampling, have popularized the concept of “conformational ensemble”, as the alternative to the analysis of PDB’s single structures. Ensembles can be analyzed to derive thermodynamic properties of the system, like entropy or free energy.²⁷ If properly built, ensembles can also be used to reconstruct complex conformational transitions or even folding events.^27–29 Ensembles are also a better reference to reproduce experimental results, as experiments measure averages of properties over a real ensemble.^30,31

Figure 2 Experimental ensembles.

Notes: (A) Superimposition of experimental structures of protein kinases. (B) NMR derived ensemble of calcium-binding protein (PDB code: 1A03).
Abbreviations: NMR, nuclear magnetic resonance; PDB, protein data bank.

Molecular dynamics simulation

Molecular dynamics (MD) simulation, first developed in the late 70s,^32,33 has advanced from simulating several hundreds of atoms to systems with biological relevance, including entire proteins in solution with explicit solvent representations, membrane embedded proteins, or large macromolecular complexes like nucleosomes^34,35 or ribosomes.^36,37 Simulation of systems having ~50,000–100,000 atoms are now routine, and simulations of approximately 500,000 atoms are common when the appropriate computer facilities are available. This remarkable improvement is in large part a consequence of the use of high performance computing (HPC), and the simplicity of the basic MD algorithm (Figure 3). An initial model of the system is obtained from either experimental structures or comparative modeling data. The simulated system could be represented at different levels of detail. Atomistic representation is the one that leads to the best reproduction of the actual systems. However, coarse-grained representations are becoming very popular when large systems or long simulations are required (see Orozco et al³⁸ for a review of such strategies). Solvent representation is a key issue in system definition. Several approaches have been assayed^39–47 but, again, the most effective is the simplest one, the explicit representation of solvent molecules, although at the expense of increasing the size of the simulated systems. Explicit solvent is able to recover most of the solvation effects of real solvent including those from entropic origin like the hydrophobic effect. Once the system is built, forces acting on every atom are obtained by deriving equations, the force-fields, where potential energy is deduced from the molecular structure.^48–53 Force-fields are complex equations, but they are easy to calculate. The simplicity of the force-field representation of molecular features: springs for bond length and angles, periodic functions for bond rotations and Lennard–Jones potentials, and the Coulomb’s law for van der Waals and electrostatic interactions, respectively, assures that energy and force calculations are extremely fast even for large systems. Force-fields currently used in atomistic molecular simulations differ in the way they are parameterized. Parameters are not necessarily interchangeable, and not all force-fields allow to represent all molecule types, but simulations conducted using modern force-fields are normally equivalent.^54,55 Once the forces acting on individual atoms are obtained, classical Newton’s law of motion is used to calculate accelerations and velocities and to update the atom positions. As integration of movement is done numerically, to avoid instability, a time step shorter than the fastest movements in the molecule should be used. This ranks normally between 1 and 2 fs for atomistic simulations, and is the major bottleneck of the simulation procedure. Microsecond-long simulations, barely scratching the time scales of biological processes, require iterating over this calculation cycle 10⁹ times. This is one of the strengths of coarse-grained strategies. As a more simplified representation of the system is used, much larger time steps are possible, and therefore the effective length of the simulations is dramatically extended. Of course, this can be obtained at the expense of the accuracy of the simulation ensemble. Algorithmic advances, that include fine-tuning of energy calculations, parallelization, or the use of graphical processing units (GPUs), have largely improved the performance of MD simulations.

Figure 3 Molecular dynamics basic algorithm.

Notes: The simulation output, the trajectory, is an ordered list of 3N atom coordinates for each simulation time (or snapshot).
Abbreviations: E_pot, potential energy; t, simulation time; dt, iteration time; For each spatial coordinate of the N simulated atoms (i): x, atom coordinate; F, forces component; a, acceleration; m, atom mass; v, velocity.

The present generation of computers takes benefit of parallelism and accelerators to speed-up the process. The most popular simulation codes (AMBER,⁵⁶ CHARMM,⁵⁷ GROMACS,⁵⁸ or NAMD⁵⁹) have long been compatible with the messaging passing interface (MPI). When a large number of computer cores can be used simultaneously, MPI can greatly reduce the computation time. To benefit the locality of interactions, the general strategy is to distribute the system to simulate among processors. This strategy is called spatial decomposition. Only a small fragment of the system has to be simulated in each processor. The most efficient division is not based in the list of particles, but in their position in space. Each processor deals with a region of space irrespective of which particles are present there. Communication between processors is also reduced, as only those simulating neighboring regions have to share information (see Larsson et al⁶⁰ for a review). As stated, the use of accelerators, mainly GPU, has become a major breakthrough in simulation codes. Originally designed to handle computer graphics, GPUs have evolved into general-purpose, fully programmable, high-performance processors and represent a major technical improvement to perform atomistic MD. Most major MD codes have already been prepared for GPUs, and even MD codes written specifically to be used on GPUs have been developed (ACEMD⁶¹). Simulation on GPUs alone or combined with MPI is, at present, the default strategy for high-throughput MD simulations. Remarkably, while simulations have been the most popular use of HPC in life sciences, the increasing power and sophistication of GPUs is leading to a greater use of personal workstations with a comparable performance.

Strategies to improve ensemble generation

Pure computational brute force, just making longer simulations, is not enough to extend the conformational sampling in biomolecular systems. The complex shape of the free energy landscape makes most of the simulations explore just a small region around the energy minimum closest to the initial conformation. With the availability of the present HPC systems, an obvious strategy is to perform a series of parallel simulations with several starting conformations. Although this could be efficient, it requires a specific knowledge of the system to simulate, and cannot be applied as a general strategy. This approach is particularly useful when several crystal structures are available (for instance in the case of allosterically regulated enzymes). A second problem that appears when collections of parallel simulations are calculated, is the generation of a usable ensemble out of the trajectories obtained. Recently the Markov state model (MSM) theory has been used to this end.^62–65 MSM theory discretizes the conformational ensemble in a collection of states, and constructs a matrix with the transition probabilities among them. The analysis of such a matrix would allow reconstruction of the global behavior of the system. Since the transition rates converge more rapidly than the population of the involved states, this approach has the advantage that the collection of simulations is not required to be especially long. This approach has been used mainly in the study of folding processes,^28,66 but also in the kinetic characterization of the formation of ligand–protein complexes.⁶⁷ Other approaches have been designed to increase the sampling space in single simulations, like metadynamics,^68–71 for instance, where already visited conformations are penalized; weighted ensembles, where additional simulations are started when new conformational spaces are visited;⁷² or accelerated MD where energy barriers are artificially reduced.^73–75 However, even with a perfect sampling, MD simulations cannot surmount barriers in the energy landscape higher than the total energy added to the system. The obtained ensemble with a single simulation is limited to those states that are accessible at the simulation temperature. Simulations at high temperatures were common in the origins of MD, but they lead to unrealistic trajectories, and hence should be combined to room temperature runs. This approach, called simulated annealing, has been largely replaced by replica exchange methods.^76,77 Such methods launch parallel simulations in different conditions. The most common variation is simulation temperature. The sampling ability of the simulation increases with temperature. Higher temperature simulations can surmount energy barriers and explore new regions of the ensemble. Periodically, energies of the different simulations are compared and structures are swapped according to its energy rank. The resulting simulation has sampled a larger conformational space, due to high temperature simulations, and retains the ability to represent the low-temperature states of the system. The main difference with the simulated annealing approach is that a realistic ensemble of the system is obtained and thermodynamic information can be derived from simulations. The idea has been extended to other simulation strategies. Most remarkably, replicas based on differences in the Hamiltonian replica exchange, including alchemical free energy calculations⁷⁸ or constant-pH simulations,^79–83 are becoming popular.

Tools to popularize MD

Preparation for simulation implies the following of a series of operations that are far from being just routine. First, the initial structure comes from the experiment. Expected issues include nonstructured or missing regions or residues, nonstandard ligands, or even structures bearing errors in the interpretation of experimental data. When a single system is simulated, all the effort in the preparation of the system is worth, as it assures the quality of the simulation result. Such setup is usually done manually, with a considerable human effort. A standard procedure to set up a system implies a number of well-known procedures: fixing structure errors; ionization of titrable amino acids; addition of structural water molecules, counter ions, and solvent; and energy minimization and equilibration of the system at the desired temperature. An expert modeler normally carries out these procedures using a set of helper programs. Such an expert has the necessary knowledge to surmount specific problems that may arise. For instance, the workflow used in the MoDEL project⁸⁴ was programmed to run automatically, but a nonnegligible fraction of over 1,500 proteins prepared failed at some point of the process. With this scenario, for newcomers to MD simulation, even a single system setup could represent an unaffordable problem. Even worse, non-expert users tend to blindly use default procedures leading easily to artifactual trajectories, which are hard to distinguish from the correct ones. This strongly contributes to the lack of popularity of biomolecular simulations among the bioinformatics or the biochemical community. MD simulations have been restricted to those research groups bearing the necessary expertise. Solving this issue requires an automatic setup of the simulation system. We would be looking for a clever black box for the nonexperts, but also for a robust software suite that can account for a large set of unrelated protein structures. All major MD codes^56–59 come with a set of accompanying programs, which perform most steps of the preparation. Additionally, a number of initiatives, combining those tools with a user-friendly interface, have come into the scene to address this problem. CHARMM-Gui⁸⁵ and CHARMMing,⁸⁶ for CHARMM, or Guimacs,⁸⁷ Gromita,⁸⁸ and jSimMacs,⁸⁹ for GROMACS, provide automatic setup functionality. VMD⁹⁰ provides a number of plug-ins that allow to launch simulations with NAMD. Most of these tools provide a friendly environment to prepare systems for simulation without the need of a deep knowledge of the underlying operations, thus facilitating the access to the field for the newcomers. Unfortunately, due to the lack of a standard for the representation of molecular simulation data, most helper applications are restricted to a single MD package, and data is not easily interchangeable. Besides, although most use some kind of embedded scripting language, automation of procedures is not a straightforward task. Lessons learned in the preparation of the MoDEL database, by our group, leaded to the generation of a new set of tools, MDMoby and MDWeb⁹¹ that try to cover both aspects of the problem. On one hand, MDMoby provides a full set of web services, covering all setup, simulation, and analysis operations. The modular nature of such collection of web services allows incorporating them as a tool kit to the design of complex setup protocols and to run them programmatically. In turn, MDWeb, a web-based interface, provides a user-friendly bench where user can check for the quality of the input structure, tailor their own setup protocols, or use a collection of predefined ones.

Application: understanding allostery

Most regulation phenomena in proteins are explained within conformational transitions. The concept of allostery that translates conformational dynamics in functional implications has been analyzed since the early times of protein biochemistry.^92–94 Conformation shift involved in allostery spans from small rearrangements to large quaternary shifts as those originally accounted by the Monod model. In any case, there is a general agreement that conformational shifts involved in allosteric transitions are simple in terms of collective movements.⁹⁵ For this reason, molecular simulations could be a natural tool to understand allostery. However, the ability of free atomistic simulation algorithms to follow a complete transition path is limited. Most of the traditional reports of simulations in this field use simplified frameworks, like discrete MD⁹⁶ or Go-Models,⁹⁷ or even popular nonsimulation equivalents like elastic network models,^98–101 and seek to find the transition path between known experimental structures. With full atom representations, it is usual to trick the algorithm by using targeted,^102–104 or supervised MD¹⁰⁵ where the simulation is artificially driven to the desired conformation. In this case, the analysis of the path could give insight into the energetics and details of the allosteric transition. For those cases where allosteric regulation is known to occur, but one of the ends is unknown, long simulations (alone^106,107 or with enhanced conformational sampling) are required.^108,109 The direct use of conformational ensembles without any conditioning is still out of the routine possibilities of present MD simulations; however, specific cases with well-defined collective motions could be feasible. Figure 4 shows an example where only 50 ns of simulation allowed for a conformational shift in Bacillus stearothermophilus lactate dehydrogenase (PDB code: 1LDN). This enzyme is known to exist in two states: a fully active, tetrameric state and a less active, dimeric one. Fructose-1,6-bisphosphate is a known allosteric regulator that allows the tetrameric state to be formed.^110,111 The construction of stable dimeric state by site-directed mutagenesis allowed for the analysis of the less active form in more detail and confirmed that a significant conformational shift occurs when the tetrameric state is formed.^112,113 The simulation of a B. stearothermophilus lactate dehydrogenase dimer protein starting with the conformation of the tetrameric one, as obtained from the experimental structure, reveals a significant intersubunit movement compatible with the experimental behavior of the enzyme (Figure 4). In this case, no computational bias was introduced; however, protein setup was done mimicking experimental conditions where the conformational shift is known to occur. The use of experimental restrains (but not necessarily the target structure) is being exploited to guide the simulation.^31,114 Alternatively, free long simulations may be replaced by shorter simulations and analyzed with MSMs to provide quantitative insight.^63–65 However target structure should be at some point explored to allow reconstructing the full process. The power of molecular simulations to uncover allosteric regulations is not in any doubt; however, there is still a long way until it could be routinely applied to all cases.

Figure 4 R-T transition on Bacillus stearothermophilus lactate dehydrogenase after 50 ns simulation in explicit solvent.

Notes: Inset: detail of conformational shift of active site substrate-binding Arg 169 side chain. Simulation was initiated from a dimeric model of the protein (1ldn), and allowed to evolve without restrains. Protein system was prepared using MDWeb,⁹¹ and simulated using GROMACS,¹⁴⁸ at 298 K of temperature, in explicit solvent and periodic boundary conditions (truncated octahedron box). Conformational shift is indicated as R-T shift in the main figure and with an arrow in the inset.
Abbreviation: R-T, relaxed and tense states (as defined by Monod’s model).

Application: molecular docking and drug design

One of the most practical application of the concept of molecular recognition are docking strategies, either small molecule or protein docking. To understand how a ligand, typically a substrate or a regulator, binds to its macromolecular counterpart is a key issue in the understanding of function itself, and it is the basis of structurally driven drug design. The recognition process is by nature dynamic.¹¹⁵ Molecules are flexible entities, and the recognition process itself implies structural rearrangements, and this shape adjustment is part of the binding process not only from the structural point of view, but also from the energetics. Although this is a generally well-accepted idea, docking algorithms are far from considering dynamic effects as a routine. Most docking or virtual screening codes work on rigid structures as obtained from the PDB. Figure 5 shows a traditional cross-docking experiment where a collection of acetylcholinesterase ligands are docked back in the same set of receptor structures. Protein structures correspond to the ones shown in Figure 1. In this experiment, all receptor structures correspond to the same protein, but crystallized with a different ligand; and all ligands are known to bind the receptor in the same place and pose. In these conditions, the experiment just measures the impact of small receptor rearrangements caused by ligand binding, on the docking efficiency. The usual result, as the one shown in Figure 5, is that even though the protein does not change, a different PDB structure implies poorer docking results. Even docking of a ligand back on its original PDB structure (diagonal results in Figure 5) tends to fail due to the usual overcompression of structures derived from X-ray crystallography. This problem is especially relevant in protein–protein docking where considerable differences are found between bound and unbound structures. For ligand-docking methods, ligand flexibility could be largely recovered by using conformer families.^116–118 In the case of protein flexibility, solutions are not so extended. Most of them use algorithms to select from a limited alternative set of protein conformations, either precomputed or simulated.^119–121 It is also possible to introduce flexibility a posteriori as a refinement process.¹²² It is still not clear whether structural adjustment comes from selection of available conformations or it is induced by the binding process itself,^123,124 but in any case the use of structural ensembles instead of single structures can be expected to improve the binding prediction.^125–127 The concept of “ensemble docking” usually requires the selection of a representative set of snapshots coming from a simulation and uses them as targets in normal docking procedures.^128–130 Integration of docking and simulation in a single calculation is less popular, but some examples have been reported.^{16,67,121,131}

Figure 5 Cross-docking experiment with selected acetylcholinesterase structures from PDB.

Notes: Performance of all possible combinations of rigid docking experiments done in standard conditions, using a series of seven acetylcholinesterase ligands extracted from left column PDB entries onto the same empty protein structures (upper row). Color code indicates RMSD between the best docked solution and the reference PDB.
Abbreviations: PDB, protein data bank; RMSD, root-mean-square deviation.

The use of simulations for the improvement of virtual screening or docking processes has a clear advantage. However, due to the speed requirements of docking, most methods based on traditional atomistic simulations are too slow to be considered, when used in a real scenario. Coarse-grained methods or any sort of accelerated MD could be a way to take benefit of simulation in a near future.

Application: refining structure predictions

Structure prediction has been one of the most ancient problems addressed in structural bioinformatics. MD, including the longest simulations performed, has been extensively used for ab initio protein structure prediction,^132–136 aiming to simulate protein folding from scratch, although this is not the preferred strategy to obtain theoretical model of protein structure. Instead, template-based modeling is the most efficient technique.^137–144 In template-based modeling, one or several 3D structures of protein showing a reasonable degree of similarity to the protein of interest are taken as templates. Irrespective of the modeling algorithm, the end result is a model bearing the new amino acid sequence and a structure averaging the used templates. In most cases, the last step of the prediction procedure implies relaxation of the structure using normally molecular mechanics. In others, restrained simulations are used throughout all the process.¹⁴⁴ The use of MD simulations looks like an obvious step in refining such models.¹⁴⁵ Simulation would allow the structure to adapt to the new sequence, and in theory give a more realistic model. Although this point is reasonable as a concept, MD simulations require systems to be close to their equilibrium (native) conformation. Otherwise, significant and difficult to detect artifacts may occur. Critical assessment of protein structure prediction contests,¹³⁹ where prediction algorithms face problems with known but nonpublic 3D structures, provide an excellent dataset to test this issue. Applying different MD approaches to the refinement of such predictions has led to a number of conclusions. The most naive approach, a single simulation starting from the predicted conformation, tends to deviate significantly from the desired structure.^145–147 Apparently, inaccuracies in the model have a large impact in the quality of the simulation results. Instead, results clearly indicate that deviation from the original structure is directly correlated with the loss of quality of the model. A second conclusion is that the ensemble of structures taken from the simulations is a closer representation of the target structure, thus indicating that the native and original structures both lie within the conformational space of the simulation.¹⁴⁶

Conclusion

MD simulations have already more than 40 years of history. However, it was not until the recent years that MD has achieved time scales that begin to be compatible with biological processes. At present, when routine simulations are approaching the microsecond scale, conformational changes, or ligand binding can be effectively simulated. The improvement of the computational equipment, especially the use of GPUs, and the improvements made in the optimization of MD algorithms, including coarse-grained ones, allow us to move from the analysis of single structures, the basis of the molecular modeling as we know it, to the analysis of conformational ensembles. Conformational ensembles are a much better representation of real macromolecules, as they account for flexibility and dynamic properties (including all thermodynamic information) and ease the match with experimental results. Although the shift in concept is clear, and the technology is coming along, there is still a long way until biomolecular simulations, the generation of conformational ensembles, would become a routine. Tools exist that make the setup of a macromolecular system much easier, and even allow the nonexperts to enter the simulation world. However, lack of representation standards, much less optimized analysis tools, and even the difficulties in simply storing and transmitting the huge amount of trajectory data that is generated are still issues that remain to be solved. In any case, MD is already a valuable tool in helping to understand biology.

Disclosure

The authors report no conflicts of interest in this work.