Skip to main content

Molecular dynamics simulations and drug discovery


This review discusses the many roles atomistic computer simulations of macromolecular (for example, protein) receptors and their associated small-molecule ligands can play in drug discovery, including the identification of cryptic or allosteric binding sites, the enhancement of traditional virtual-screening methodologies, and the direct prediction of small-molecule binding energies. The limitations of current simulation methodologies, including the high computational costs and approximations of molecular forces required, are also discussed. With constant improvements in both computer power and algorithm design, the future of computer-aided drug design is promising; molecular dynamics simulations are likely to play an increasingly important role.


Richard Feynman, recipient of the 1965 Nobel Prize in Physics, once famously stated: 'If we were to name the most powerful assumption of all, which leads one on and on in an attempt to understand life, it is that all things are made of atoms, and that everything that living things do can be understood in terms of the jigglings and wigglings of atoms.' Much of the biophysics of the last 50 years has been dedicated to better understanding the nature of this atomic jiggling and wiggling. The quantum-mechanical laws governing motions in the microscopic world are surprisingly foreign to those familiar with macroscopic dynamics. Motions are governed not by deterministic laws, but by probability functions; chemical bonds are formed not mechanically, but by shifting clouds of electrons that are simultaneously waves and particles. As Feynman eloquently put it, this is 'nature as she is - absurd' [1].

Understanding these absurd molecular motions is undoubtedly germane to drug discovery. The initial 'lock-and-key' theory of ligand binding [2], in which a frozen, motionless receptor was thought to accommodate a small molecule without undergoing any conformational rearrangements, has been largely abandoned in favor of binding models that account not only for conformational changes, but also for the random jiggling of receptors and ligands [37].

The mollusk acetylcholine binding protein (AChBP), a structural and functional surrogate of the human nicotinic acetylcholine receptor (AChR) ligand-binding domain [811], is one of countless examples that illustrate the importance of accounting for these atomistic motions (Figure 1). In crystallographic structures of small-molecule AChR agonists bound to AChBP, a key loop (loop C) partially closes around the ligand (Figure 1a,c). In contrast, crystal structures of large AChR antagonists like snake α-neurotoxins bound to AChBP reveal that this same loop is displaced by as much as 10 Å, producing an active site that is far more open (Figure 1b,c) [12]. Bourne et al. [12] proposed that the unbound AChBP and AChR are highly dynamic proteins that, in the absence of a ligand, sample many conformational states, both open and closed, that are selectively stabilized by the binding of agonists and antagonists. Any one of these binding-pocket conformations may be druggable and therefore pharmacologically relevant. If this general model of ligand binding is correct, the implications for structure-based drug design are profound, as the same principle of binding likely applies to many other systems as well.

Figure 1

The different conformations of the acetylcholine binding protein from Lymnaea stagnalis. Portions of the protein have been removed to facilitate visualization. (a) The protein in a closed conformation with nicotine bound (PDB ID: 1UW6), shown in blue. (b) The protein in an open conformation (PDB ID: 1YI5) with the same nicotine conformation superimposed on the structure, shown in pink. (c) Ribbon representations of the two conformations.

Molecular dynamics simulations

While crystallographic studies like these convincingly demonstrate the important role protein flexibility plays in ligand binding, the expense and extensive labor required to generate them have led many to seek computational techniques that can predict protein motions. Unfortunately, the calculations required to describe the absurd quantum-mechanical motions and chemical reactions of large molecular systems are often too complex and computationally intensive for even the best supercomputers. Molecular dynamics (MD) simulations, first developed in the late 1970s [13], seek to overcome this limitation by using simple approximations based on Newtonian physics to simulate atomic motions, thus reducing the computational complexity. The general process of approximation used is outlined in Figure 2. First, a computer model of the molecular system is prepared from nuclear magnetic resonance (NMR), crystallographic, or homology-modeling data. The forces acting on each of the system atoms are then estimated from an equation like that shown in Figure 3 [14]. In brief, forces arising from interactions between bonded and non-bonded atoms contribute. Chemical bonds and atomic angles are modeled using simple virtual springs, and dihedral angles (that is, rotations about a bond) are modeled using a sinusoidal function that approximates the energy differences between eclipsed and staggered conformations. Non-bonded forces arise due to van der Waals interactions, modeled using the Lennard-Jones 6-12 potential [15], and charged (electrostatic) interactions, modeled using Coulomb's law. For a more in-depth review describing how the equations describing these interactions are parameterized, see [14].

Figure 2

A schematic showing how a molecular dynamics simulation is performed. First, a computer model of the receptor-ligand system is prepared. An equation like that shown in Figure 3 is used to estimate the forces acting on each of the system atoms. The positions of the atoms are moved according to Newton's laws of motion. The simulation time is advanced, and the process is repeated many times. This figure was adapted from a version originally created by Kai Nordlund.

Figure 3

An example of an equation used to approximate the atomic forces that govern molecular movement. The atomic forces that govern molecular movement can be divided into those caused by interactions between atoms that are chemically bonded to one another and those caused by interactions between atoms that are not bonded. Chemical bonds and atomic angles are modeled using simple springs, and dihedral angles (that is, rotations about a bond) are modeled using a sinusoidal function that approximates the energy differences between eclipsed and staggered conformations. Non-bonded forces arise due to van der Waals interactions, modeled using the Lennard-Jones potential, and charged (electrostatic) interactions, modeled using Coulomb's law.

In order to reproduce the actual behavior of real molecules in motion, the energy terms described above are parameterized to fit quantum-mechanical calculations and experimental (for example, spectroscopic) data. This parameterization includes identifying the ideal stiffness and lengths of the springs that describe chemical bonding and atomic angles, determining the best partial atomic charges used for calculating electrostatic-interaction energies, identifying the proper van der Waals atomic radii, and so on. Collectively, these parameters are called a 'force field' because they describe the contributions of the various atomic forces that govern molecular dynamics. Several force fields are commonly used in molecular dynamics simulations, including AMBER [14, 16], CHARMM [17], and GROMOS [18]. These differ principally in the way they are parameterized but generally give similar results.

Once the forces acting on each of the system atoms have been calculated, the positions of these atoms are moved according to Newton's laws of motion. The simulation time is then advanced, often by only 1 or 2 quadrillionths of a second, and the process is repeated, typically millions of times. Because so many calculations are required, molecular dynamics simulations are typically performed on computer clusters or supercomputers using dozens if not hundreds of processors in parallel. Many of the most popular simulation software packages, which often bear the same names as their default force fields (for example AMBER [19], CHARMM [17], and NAMD [20, 21]), are compatible with the Message Passing Interface (MPI), a system of computer-to-computer messaging that greatly facilitates the execution of complex tasks by one software application on multiple processors operating simultaneously.

While the utility of molecular dynamics simulations should not be overstated, a number of studies comparing these simulations with experimental data have been used to validate the computational technique [22]. NMR data are particularly useful, as the many receptor and ligand conformations sampled by molecular dynamics simulations can be used to predict NMR measurements like spin relaxation, permitting direct comparison between experimental and theoretical techniques. Indeed, a number of studies have shown good agreement between computational and experimental measurements of macromolecular dynamics [2326].

Molecular dynamics simulations: current limitations

These successes aside, the utility of molecular dynamics simulations is still limited by two principal challenges [27]: the force fields used require further refinement, and high computational demands prohibit routine simulations greater than a microsecond in length, leading in many cases to an inadequate sampling of conformational states. As an example of these high computational demands, consider that a one-microsecond simulation of a relatively small system (approximately 25,000 atoms) running on 24 processors takes several months to complete.

Aside from challenges related to the high computational demands of these simulations, the force fields used are also approximations of the quantum-mechanical reality that reigns in the atomic regime. While simulations can accurately predict many important molecular motions, these simulations are poorly suited to systems where quantum effects are important, for example, when transition metal atoms are involved in binding.

To overcome this challenge, some researchers have introduced quantum mechanical calculations into classic molecular-dynamics force fields; the motions and reactions of enzymatic active sites or other limited areas of interest are simulated according to the laws of quantum mechanics, and the motions of the larger system are approximated using molecular dynamics. While far from the computationally intractable 'ideal' of using quantum mechanics to describe the entire system, this hybrid technique has nevertheless been used successfully to study a number of systems. For example, in one recent simulation of Desulfovibrio desulfuricans and Clostridium pasteurianum [Fe-Fe] hydrogenases, a 'QM [quantum mechanical] region' was defined encompassing a metal-containing region of the protein thought to be catalytically important, and the remainder of the protein was simulated using classical molecular dynamics [28]. The simulations revealed an important proton transfer in the QM region, a bond-breaking and bond-formation event that could not have been modeled with a traditional force field. The hypothesized catalytic mechanism was subsequently supported by experimental evidence.

Aside from bond breaking and formation, electronic polarization, caused by the flow of electrons from one atomic nucleus to another among groups of atoms that are chemically bonded, is another quantum-mechanical effect that, with few exceptions, is generally ignored. In classical molecular dynamics simulations, each atom is assigned a fixed partial charge before the simulation begins. In reality, however, the electron clouds surrounding atoms are constantly shifting according to their environments, so that the partial charges of atoms would be better represented as dynamic and responsive. Despite wide agreement on the importance of accounting for electronic polarization, after 30 years of development a generally accepted polarizable force field has not been forthcoming, and molecular dynamics simulations using those force fields that are available are rare [29]. Nevertheless, a number of polarizable force fields are currently under development [30], and future implementations may lead to improved accuracy.

In addition to neglecting quantum-mechanical effects, molecular dynamics studies are also limited by the short time scales typically simulated. To reproduce thermodynamic properties and/or to fully elucidate all binding-pocket configurations relevant to drug design, all the possible conformational states of the protein must be explored by the simulation. Unfortunately, many biochemical processes, including receptor conformational shifts relevant to drug binding, occur on time scales that are much longer than those amenable to simulation. With some important exceptions [31], simulations are currently limited to at most millionths of a second; indeed, most simulations are measured in billionths of a second.

A number of solutions to this challenge have already seen limited use. For example, in accelerated molecular dynamics (aMD) [32, 33], large energy barriers are artificially reduced. Though this process inevitably introduces some artifacts, it does allow proteins to shift between conformations that would not be accessible given the time scales of conventional molecular dynamics. These novel conformations can then be further studied using classical molecular dynamics or other techniques.

Novel hardware has also been used to overcome the time-scale limitations of conventional molecular dynamics simulations. Many of the same calculations required for these simulations are commonly performed by video-game and computer-graphics applications. Consequently, the same graphics-processing-units (GPUs) designed to speed up video games can be used to speed up molecular dynamics simulations as well, usually by an order of magnitude [3436].

Not satisfied with merely adapting molecular-dynamics code to run on specialized graphics processors, some engineers have designed new processors specifically for these simulations. The research group of DE Shaw is one notable advocate of this approach. They have built a supercomputer codenamed Anton capable of performing microseconds of simulation per day. With Anton, simulations longer than one millisecond [31] have successfully captured protein folding and unfolding as well as drug-binding events [37]. Shortcomings certainly still exist, but these and other future techniques will likely make great progress towards overcoming current limitations on conformational sampling.

Molecular dynamics simulations and drug discovery

Weaknesses in current force fields and conformational sampling aside, molecular dynamics simulations and the insights they offer into protein motion often play important roles in drug discovery. Just as a single photograph of a runner tells little about her stride, a single protein conformation tells little about protein dynamics. The static models produced by NMR, X-ray crystallography, and homology modeling provide valuable insights into macromolecular structure, but molecular recognition and drug binding are very dynamic processes. When a small molecule like a drug (for example, a ligand) approaches its target (for example, a receptor) in solution, it encounters not a single, frozen structure, but rather a macromolecule in constant motion.

In some rare cases, protein motions are limited, and the ligand may fit into a fairly static binding pocket like a key fits into a lock [2]. More typically, the ligand may bind and stabilize only a subset of the many conformations sampled by its dynamic receptor, causing the population of all possible conformations to shift towards those that can best accommodate binding [47]. Upon binding, the ligand may further induce conformational changes that are not typically sampled when the ligand is absent [38]. Regardless, receptor motions clearly play an essential role in the binding of most small-molecule drugs. Several techniques have been developed to exploit the information about these motions that molecular dynamics simulations can provide.

Identifying cryptic and allosteric binding sites

NMR and X-ray crystallographic structures often reveal well defined binding pockets that accommodate endogenous ligands; however, sometimes the models produced by these experimental techniques obscure other potentially druggable sites. As these sites are not immediately obvious from available structures, they are sometimes called cryptic binding sites.

Molecular dynamics simulations are excellent tools for identifying such sites [3941]. For example, in 2004 Schames et al. [39] performed a molecular dynamics simulation of HIV integrase, a drug target that had not seemed amenable to structure-based drug design. The simulations revealed a previously unidentified trench that was not evident from any of the available crystal structures. X-ray crystallography later demonstrated that known inhibitors do in fact bind in this cryptic trench, as predicted. These results led to new experimental studies at Merck & Co. [42]; further development ultimately resulted in production of the highly effective antiretroviral drug raltegravir, the first US Food and Drug Administration-approved HIV integrase inhibitor.

Aside from cryptic binding sites, molecular dynamics simulations can also be used to identify druggable allosteric sites. In one recent study, Ivetac and McCammon [43] performed simulations of the human β11AR) and β22AR) adrenergic receptors. Multiple protein conformations were extracted from these simulations, and the protein surface was computationally 'flooded' with small organic probes using FTMAP [44] to identify potential binding sites. Regions of the protein surface where the organic probes consistently congregated across multiple structures were then identified as potential allosteric sites. In all, five potential sites were identified, some of which were not evident in any of the existing crystal structures.

Improving the computational identification of true small-molecule binders: the relaxed complex scheme

One common technique used to identify the precursors of potential drugs in silico is virtual screening. A docking program is used to predict the binding pose and energy of a small-molecule model within a selected receptor binding pocket. Traditionally, many ligand models, typically taken from a database of compounds that can be easily synthesized or commercially purchased, are docked into a single static receptor structure, often obtained from NMR or X-ray crystallography. The best predicted ligands are subsequently tested experimentally to confirm binding.

Unfortunately, traditional docking relying on a single receptor structure is problematic. Some legitimate ligands may indeed bind to the single structure selected, but in reality most receptor binding pockets have many valid conformational states, any one of which may be druggable. In a traditional virtual screen, true ligands are often discarded because they in fact bind to receptor conformations that differ markedly from that of the single static structure chosen.

To better account for receptor flexibility, a new virtual-screening protocol has been developed called the relaxed complex scheme (RCS) [45, 46]. Rather than docking many compound models into a single NMR or crystal structure, each potential ligand is docked into multiple protein conformations, typically extracted from a molecular dynamics simulation. Thus, each ligand is associated not with a single docking score but rather with a whole spectrum of scores. Ligands can be ranked by a number of spectrum characteristics, such as the average score over all receptors. Thus, the RCS effectively accounts for the many receptor conformations sampled by the simulations; it has been used successfully to identify a number of protein inhibitors, including inhibitors of FKBP [47], HIV integrase [39], Trypanosoma brucei RNA editing ligase 1 [48, 49], T. brucei GalE [50], T. brucei FPPS [51], and Mycobacterium tuberculosis dTDP-6-deoxy-L-lyxo-4-hexulose [52]. In two of these projects, the identified inhibitors were effective not only against the target proteins, but against the whole-cell parasite [49, 50].

While these successes are promising, the relaxed complex scheme certainly has its weaknesses. Aside from being based on molecular dynamics simulations that are themselves subject to crude force-field approximations and inadequate conformational sampling, the scheme relies on computer-docking scoring functions that of necessity are optimized for speed at the expense of accuracy. In order to facilitate high-throughput virtual screening, these scoring functions often treat subtle influences on binding energy like conformational entropy and solvation energy only superficially [27, 53], thus sacrificing accuracy for the sake of greater speed.

Advanced free-energy calculations using molecular dynamics simulations

Though docking programs are optimized for speed rather than accuracy, more accurate, albeit computationally intensive, techniques for predicting binding affinity do exist. These techniques, which include thermodynamic integration [54], single-step perturbation [55], and free energy perturbation [56], are based in large part on molecular dynamics simulations.

Free energy is a state function, meaning that the free-energy difference associated with a given event like a drug binding to its receptor is determined only by the energy prior to that event and the energy following it; while the path taken from the initial to the final state may influence receptor-ligand kinetics, it has no bearing on the free energy. Perhaps the ligand diffuses slowly towards the active site and slips easily into the binding pocket. Perhaps the protein unfolds entirely and then refolds around the ligand. Perhaps the ligand in solution is beamed to a starship in orbit, only to rematerialize in the active site a few seconds later. The mechanics do not matter; the free energy depends only on the initial energy in solution and the final energy following the binding event.

With some notable exceptions [37], simulating a receptor-ligand system long enough to capture an entire binding event is not currently feasible. However, it is still possible to calculate a drug's binding affinity using a technique called 'alchemical transformation', first described in 1984 [57]. This transformation is not too different from the starship example given above. During the course of a molecular dynamics simulation, the electrostatic and van der Waals forces produced by ligand atoms are turned down gradually enough to avoid undesirable artifacts. Eventually, the ligand is no longer able to interact with the protein or solvent. For all practical purposes, the ligand has disappeared. It does not matter that this transformation is not at all physical; free energy is a state function, so the path from the initial to the final state, whether real or imaginary, is irrelevant.

It is not clear, however, in what context the ligand should be figuratively annihilated in this way. Should a molecular dynamics simulation be run in which the bound ligand vanishes? What about the ligand in solution? To address these questions, alchemical transformations are selected based on the thermodynamic cycle shown in Figure 4. As free energy is a state function that depends only on the energy of the initial and final states, a system that proceeds from one state around this free-energy cycle only to return to the same initial state should have no change in total free energy (that is, ΔGbind + ΔGprotein - ΔG - ΔGwater = 0). We note that ΔG in this equation is itself zero, since the ligand is entirely disappeared in both of the states shown in the bottom half of Figure 4, meaning the ligand cannot interact with the water solvent or the protein in either state. Thus, ΔGbind + ΔGprotein - ΔGwater = 0, and, consequently, ΔGbind = ΔGwater - ΔGprotein. These equations demonstrate that it is possible to estimate a drug's free energy of binding, an indirect measurement of drug potency, by running two simulations, one in which the receptor-bound ligand disappears, and one in which the solvated ligand disappears.

Figure 4

The thermodynamic cycle used for selecting alchemical transformations. Typically, one wishes to calculate the free energy of binding, ΔGbind, shown across the top. However, it is generally impractical to run a molecular dynamics simulation long enough to capture an entire binding event. Instead, a series of alchemical transformations are performed using molecular dynamics simulations. ΔGprotein is the change in free energy that occurs when a bound ligand is 'annihilated'. ΔG is the change in free energy that occurs when an unbound 'ghost' ligand binds to the receptor; however, since a ghost ligand is not able to interact with any solvent or receptor atoms, this energy is always zero. Finally, ΔGwater is the change in free energy that occurs when an unbound ligand in solution is 'annihilated'. A system that proceeds from one state around this free-energy cycle only to return to the same initial state should have no change in total free energy; consequently, ΔGbind = ΔGwater - ΔGprotein.

A similar task, calculating relative ligand binding energies, is useful during drug optimization when one wishes to determine if a given chemical change will improve the potency of a candidate ligand. In this case, rather than annihilating the entire ligand, one section of the ligand is gradually transformed. For example, a key carbon atom might be gradually converted into an oxygen atom to see if the binding affinity is improved or diminished. These kinds of alchemical molecular dynamics simulations may provide medicinal chemists with useful insights that can guide further drug development.

A series of early fortuitous results that agreed remarkably well with experiments led many to an enthusiasm for molecular-dynamics-based free-energy calculations in the 1980s and early 1990s that was not matched in subsequent decades [27, 58] as computational predictions fell short of experimental measurements. However, steady algorithmic and engineering advances in recent years have led to renewed attention. The successful applications of alchemical techniques in recent years are legion; accurate predictions have been obtained for ligand binding to the src SH2 domain [59], to a mutant T4 lysozyme [60], to FKBP12 [61], to HIV reverse transcriptase [62], to trypsin [63], to a bacterial ribosome [64], and to estrogen receptor-α [65], among many others.

These successes aside, it is important not to oversell alchemical techniques. All molecular-dynamics-based drug-discovery techniques would benefit from improved force fields, but the alchemical techniques are, in addition, uniquely sensitive to inadequate conformational sampling [66]. When molecular dynamics simulations of insufficient length are used to identify cryptic sites, allosteric sites, or pharmacologically relevant binding-pocket conformations for virtual-screening projects, the risk is that some suitable receptor conformations will be missed. The conformations that are identified, however, are still useful; the results of the simulation are therefore incomplete, but not necessarily wrong.

However, the alchemical techniques used to calculate free energies of binding are far more dependent on thorough conformational sampling than are RCS screens. If molecular dynamics simulations fail to sample system conformations that are in fact sampled ex silico, these conformations will not contribute to the total calculated energy, leading to an incorrect prediction of the binding affinity. Molecular dynamics simulations are computationally demanding and often of necessity unacceptably short; insufficient conformational sampling is therefore a common problem that future algorithmic and hardware-engineering efforts must address. It is largely because accurate predictions often require lengthy simulations that these alchemical techniques have not yet been widely adopted by the pharmaceutical industry, despite great interest [27].


In this review, we have discussed the many roles that molecular dynamics simulations can play in drug discovery, including the identification of cryptic or allosteric binding sites, the enhancement of traditional virtual-screening methodologies, and the direct prediction of ligand binding energies. As ligand binding and the important macromolecular motions associated with it are microscopic events that take place in mere millionths of a second, a complete understanding of the atomistic energetics and mechanics of binding is unattainable using current experimental techniques. Molecular dynamics simulations are useful for filling in the details where experimental methods cannot.

With constant improvements in both computer power and algorithm design, the future of computer-aided drug design is promising; molecular dynamics simulations are likely to play an increasingly important role in the development of novel pharmacological therapeutics.


  1. 1.

    Feynman RP: QED: the Strange Theory of Light and Matter. 1985, Princeton, NJ: Princeton University Press

    Google Scholar 

  2. 2.

    Fischer E: Einfluss der Configuration auf die Wirkung der Enzyme. Ber Dtsch Chem Ges. 1894, 27: 2985-2993. 10.1002/cber.18940270364.

    Article  CAS  Google Scholar 

  3. 3.

    Teague SJ: Implications of protein flexibility for drug discovery. Nat Rev Drug Discov. 2003, 2: 527-541. 10.1038/nrd1129.

    Article  CAS  PubMed  Google Scholar 

  4. 4.

    Ma B, Kumar S, Tsai CJ, Nussinov R: Folding funnels and binding mechanisms. Protein Eng. 1999, 12: 713-720. 10.1093/protein/12.9.713.

    Article  CAS  PubMed  Google Scholar 

  5. 5.

    Kumar S, Ma B, Tsai CJ, Wolfson H, Nussinov R: Folding funnels and conformational transitions via hinge-bending motions. Cell Biochem Biophys. 1999, 31: 141-164. 10.1007/BF02738169.

    Article  CAS  PubMed  Google Scholar 

  6. 6.

    Tsai CJ, Kumar S, Ma B, Nussinov R: Folding funnels, binding funnels, and protein function. Protein Sci. 1999, 8: 1181-1190. 10.1110/ps.8.6.1181.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  7. 7.

    Ma B, Shatsky M, Wolfson HJ, Nussinov R: Multiple diverse ligands binding at a single protein site: a matter of pre-existing populations. Protein Sci. 2002, 11: 184-197.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  8. 8.

    Bouzat C, Gumilar F, Spitzmaul G, Wang HL, Rayes D, Hansen SB, Taylor P, Sine SM: Coupling of agonist binding to channel gating in an ACh-binding protein linked to an ion channel. Nature. 2004, 430: 896-900. 10.1038/nature02753.

    Article  CAS  PubMed  Google Scholar 

  9. 9.

    Talley TT, Yalda S, Ho KY, Tor Y, Soti FS, Kem WR, Taylor P: Spectroscopic analysis of benzylidene anabaseine complexes with acetylcholine binding proteins as models for ligand-nicotinic receptor interactions. Biochemistry. 2006, 45: 8894-8902. 10.1021/bi060534y.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  10. 10.

    Babakhani A, Talley TT, Taylor P, McCammon JA: A virtual screening study of the acetylcholine binding protein using a relaxed-complex approach. Comput Biol Chem. 2009, 33: 160-170. 10.1016/j.compbiolchem.2008.12.002.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  11. 11.

    Ulens C, Hogg RC, Celie PH, Bertrand D, Tsetlin V, Smit AB, Sixma TK: Structural determinants of selective alpha-conotoxin binding to a nicotinic acetylcholine receptor homolog AChBP. Proc Natl Acad Sci USA. 2006, 103: 3615-3620. 10.1073/pnas.0507889103.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  12. 12.

    Bourne Y, Talley TT, Hansen SB, Taylor P, Marchot P: Crystal structure of a Cbtx-AChBP complex reveals essential interactions between snake alpha-neurotoxins and nicotinic receptors. EMBO J. 2005, 24: 1512-1522. 10.1038/sj.emboj.7600620.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  13. 13.

    McCammon JA, Gelin BR, Karplus M: Dynamics of folded proteins. Nature. 1977, 267: 585-590. 10.1038/267585a0.

    Article  CAS  PubMed  Google Scholar 

  14. 14.

    Cornell WD, Cieplak P, Bayly CI, Gould IR, Merz KM, Ferguson DM, Spellmeyer DC, Fox T, Caldwell JW, Kollman PA: A second generation force field for the simulation of proteins, nucleic acids, and organic molecules. J Am Chem Soc. 1995, 117: 5179-5197. 10.1021/ja00124a002.

    Article  CAS  Google Scholar 

  15. 15.

    Lennard-Jones JE: On the determination of molecular fields. II. From the equation of state of a gas. Proc R Soc Lond A. 1924, 106: 463-477. 10.1098/rspa.1924.0082.

    Article  Google Scholar 

  16. 16.

    Wang J, Wolf RM, Caldwell JW, Kollman PA, Case DA: Development and testing of a general amber force field. J Comput Chem. 2004, 25: 1157-1174. 10.1002/jcc.20035.

    Article  CAS  PubMed  Google Scholar 

  17. 17.

    Brooks BR, Bruccoleri RE, Olafson BD, States DJ, Swaminathan S, Karplus M: CHARMM - a program for macromolecular energy, minimization, and dynamics calculations. J Comput Chem. 1983, 4: 187-217. 10.1002/jcc.540040211.

    Article  CAS  Google Scholar 

  18. 18.

    Christen M, Hünenberger PH, Bakowies D, Baron R, Bürgi R, Geerke DP, Heinz TN, Kastenholz MA, Kräutler V, Oostenbrink C, Peter C, Trzesniak D, van Gunsteren WF: The GROMOS software for biomolecular simulation: GROMOS05. J Comput Chem. 2005, 26: 1719-1751. 10.1002/jcc.20303.

    Article  CAS  PubMed  Google Scholar 

  19. 19.

    Case DA, Cheatham TE, Darden T, Gohlke H, Luo R, Merz KM, Onufriev A, Simmerling C, Wang B, Woods RJ: The AMBER biomolecular simulation programs. J Comput Chem. 2005, 26: 1668-1688. 10.1002/jcc.20290.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  20. 20.

    Kale L, Skeel R, Bhandarkar M, Brunner R, Gursoy A, Krawetz N, Phillips J, Shinozaki A, Varadarajan K, Schulten K: NAMD2: greater scalability for parallel molecular dynamics. J Comput Phys. 1999, 151: 283-312. 10.1006/jcph.1999.6201.

    Article  CAS  Google Scholar 

  21. 21.

    Phillips JC, Braun R, Wang W, Gumbart J, Tajkhorshid E, Villa E, Chipot C, Skeel RD, Kale L, Schulten K: Scalable molecular dynamics with NAMD. J Comput Chem. 2005, 26: 1781-1802. 10.1002/jcc.20289.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  22. 22.

    van Gunsteren WF, Dolenc J, Mark AE: Molecular simulation as an aid to experimentalists. Curr Opin Struc Biol. 2008, 18: 149-153. 10.1016/

    Article  CAS  Google Scholar 

  23. 23.

    LaConte LEW, Voelz VA, Nelson WD, Thomas DD: Molecular dynamics simulation of site-directed spin labeling: Experimental validation in muscle fibers. BiophysJ. 2002, 82: 484A-484A.

    Google Scholar 

  24. 24.

    Peter C, Rueping M, Worner HJ, Jaun B, Seebach D, van Gunsteren WEF: Molecular dynamics simulations of small peptides: can one derive conformational preferences from ROESY spectra?. Chemistry. 2003, 9: 5838-5849. 10.1002/chem.200305147.

    Article  CAS  PubMed  Google Scholar 

  25. 25.

    Bruschweiler R, Showalter SA: Validation of molecular dynamics simulations of biomolecules using NMR spin relaxation as benchmarks: Application to the AMBER99SB force field. J Chem Theory Comput. 2007, 3: 961-975. 10.1021/ct7000045.

    Article  PubMed  Google Scholar 

  26. 26.

    Markwick PRL, Cervantes CF, Abel BL, Komives EA, Blackledge M, McCammon JA: Enhanced conformational space sampling improves the prediction of chemical shifts in proteins. J Am Chem Soc. 2010, 132: 1220-1221. 10.1021/ja9093692.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  27. 27.

    Chodera JD, Mobley DL, Shirts MR, Dixon RW, Branson K, Pande VS: Alchemical free energy methods for drug discovery: progress and challenges. Curr Opin Struct Biol. 2011, 21: 150-160. 10.1016/

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  28. 28.

    Hong G, Cornish AJ, Hegg EL, Pachter R: On understanding proton transfer to the biocatalytic [Fe-Fe](H) sub-cluster in [Fe-Fe]H(2)ases: QM/MM MD simulations. Biochim Biophys Acta. 2011, 1807: 510-517. 10.1016/j.bbabio.2011.01.011.

    Article  CAS  PubMed  Google Scholar 

  29. 29.

    Jorgensen WL: Special issue on polarization. J Chem Theory Comput. 2007, 3: 1877-10.1021/ct700252g.

    Article  CAS  PubMed  Google Scholar 

  30. 30.

    Cieplak P, Dupradeau FY, Duan Y, Wang J: Polarization effects in molecular mechanical force fields. J Phys Condens Matter. 2009, 21: 333102-10.1088/0953-8984/21/33/333102.

    PubMed Central  Article  PubMed  Google Scholar 

  31. 31.

    Shaw DE, Maragakis P, Lindorff-Larsen K, Piana S, Dror RO, Eastwood MP, Bank JA, Jumper JM, Salmon JK, Shan Y, Wriggers W: Atomic-level characterization of the structural dynamics of proteins. Science. 2010, 330: 341-346. 10.1126/science.1187409.

    Article  CAS  PubMed  Google Scholar 

  32. 32.

    Hamelberg D, Mongan J, McCammon JA: Accelerated molecular dynamics: A promising and efficient simulation method for biomolecules. J Chem Phys. 2004, 120: 11919-11929. 10.1063/1.1755656.

    Article  CAS  PubMed  Google Scholar 

  33. 33.

    Hamelberg D, McCammon JA: Fast peptidyl cis-trans isomerization within the flexible Gly-rich flaps of HIV-1 protease. J Am Chem Soc. 2005, 127: 13778-13779. 10.1021/ja054338a.

    Article  CAS  PubMed  Google Scholar 

  34. 34.

    Yang J, Wang Y, Chen Y: GPU accelerated molecular dynamics simulation of thermal conductivities. J Comp Physics. 2007, 221: 799-804. 10.1016/

    Article  Google Scholar 

  35. 35.

    Liu W, Schmidt B, Voss G, Müller-Wittig W: Accelerating molecular dynamics simulations using Graphics Processing Units with CUDA. Comput Phys Commun. 2008, 179: 634-641. 10.1016/j.cpc.2008.05.008.

    Article  CAS  Google Scholar 

  36. 36.

    Friedrichs MS, Eastman P, Vaidyanathan V, Houston M, Legrand S, Beberg AL, Ensign DL, Bruns CM, Pande VS: Accelerating molecular dynamic simulation on graphics processing units. J Comput Chem. 2009, 30: 864-872. 10.1002/jcc.21209.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  37. 37.

    Shan Y, Kim ET, Eastwood MP, Dror RO, Seeliger MA, Shaw DE: How does a drug molecule find its target binding site?. J Am Chem Soc. 2011, 133: 9181-9183. 10.1021/ja202726y.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  38. 38.

    Koshland DE: Application of a theory of enzyme specificity to protein synthesis. Proc Natl Acad Sci USA. 1958, 44: 98-104. 10.1073/pnas.44.2.98.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  39. 39.

    Schames JR, Henchman RH, Siegel JS, Sotriffer CA, Ni H, McCammon JA: Discovery of a novel binding trench in HIV integrase. J Med Chem. 2004, 47: 1879-1881. 10.1021/jm0341913.

    Article  CAS  PubMed  Google Scholar 

  40. 40.

    Durrant JD, Keranen H, Wilson BA, McCammon JA: Computational identification of uncharacterized cruzain binding sites. PLoS Negl Trop Dis. 2010, 4: e676-10.1371/journal.pntd.0000676.

    PubMed Central  Article  PubMed  Google Scholar 

  41. 41.

    Grant BJ, Lukman S, Hocker HJ, Sayyah J, Brown JH, McCammon JA, Gorfe AA: Novel allosteric sites on Ras for lead generation. PLoS One. 2011.

    Google Scholar 

  42. 42.

    Hazuda DJ, Anthony NJ, Gomez RP, Jolly SM, Wai JS, Zhuang L, Fisher TE, Embrey M, Guare JP, Egbertson MS, Vacca JP, Huff JR, Felock PJ, Witmer MV, Stillmock KA, Danovich R, Grobler J, Miller MD, Espeseth AS, Jin L, Chen IW, Lin JH, Kassahun K, Ellis JD, Wong BK, Xu W, Pearson PG, Schleif WA, Cortese R, Emini E, et al: A naphthyridine carboxamide provides evidence for discordant resistance between mechanistically identical inhibitors of HIV-1 integrase. Proc Natl Acad Sci USA. 2004, 101: 11233-11238. 10.1073/pnas.0402357101.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  43. 43.

    Ivetac A, McCammon JA: Mapping the druggable allosteric space of G-protein coupled receptors: a fragment-based molecular dynamics approach. Chem Biol Drug Des. 2010, 76: 201-217.

    PubMed Central  CAS  PubMed  Google Scholar 

  44. 44.

    Brenke R, Kozakov D, Chuang GY, Beglov D, Hall D, Landon MR, Mattos C, Vajda S: Fragment-based identification of druggable 'hot spots' of proteins using Fourier domain correlation techniques. Bioinformatics. 2009, 25: 621-627. 10.1093/bioinformatics/btp036.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  45. 45.

    Amaro RE, Baron R, McCammon JA: An improved relaxed complex scheme for receptor flexibility in computer-aided drug design. J Comput Aided Mol Des. 2008, 22: 693-705. 10.1007/s10822-007-9159-2.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  46. 46.

    Lin JH, Perryman AL, Schames JR, McCammon JA: The relaxed complex method: Accommodating receptor flexibility for drug design with an improved scoring scheme. Biopolymers. 2003, 68: 47-62. 10.1002/bip.10218.

    Article  CAS  PubMed  Google Scholar 

  47. 47.

    Lin JH, Perryman AL, Schames JR, McCammon JA: Computational drug design accommodating receptor flexibility: the relaxed complex scheme. J Am Chem Soc. 2002, 124: 5632-5633. 10.1021/ja0260162.

    Article  CAS  PubMed  Google Scholar 

  48. 48.

    Amaro RE, Schnaufer A, Interthal H, Hol W, Stuart KD, McCammon JA: Discovery of drug-like inhibitors of an essential RNA-editing ligase in Trypanosoma brucei. Proc Natl Acad Sci USA. 2008, 105: 17278-17283. 10.1073/pnas.0805820105.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  49. 49.

    Durrant JD, Hall L, Swift RV, Landon M, Schnaufer A, Amaro RE: Novel naphthalene-based inhibitors of Trypanosoma brucei RNA editing ligase 1. PLoS Negl Trop Dis. 2010, 4: e803-10.1371/journal.pntd.0000803.

    PubMed Central  Article  PubMed  Google Scholar 

  50. 50.

    Durrant JD, Urbaniak MD, Ferguson MA, McCammon JA: Computer-aided identification of Trypanosoma brucei uridine diphosphate galactose 4'-epimerase inhibitors: toward the development of novel therapies for African sleeping sickness. J Med Chem. 2010, 53: 5025-5032. 10.1021/jm100456a.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  51. 51.

    Durrant JD, Cao R, Gorfe AA, Zhu W, Li J, Sankovsky A, Oldfield E, McCammon JA: Non-bisphosphonate inhibitors of isoprenoid biosynthesis identified via computer-aided drug design. Chem Biol Drug Des. 2011, 78: 323-332. 10.1111/j.1747-0285.2011.01164.x.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  52. 52.

    Wang Y, Hess TN, Jones V, Zhou JZ, McNeil MR, McCammon JA: Novel inhibitors of Mycobacterium tuberculosisdTDP-6-deoxy-L-lyxo-4-hexulose reductase (RmlD) identified by virtual screening. Bioorg Med Chem Lett. 2011,

    Google Scholar 

  53. 53.

    Kitchen DB, Decornez H, Furr JR, Bajorath J: Docking and scoring in virtual screening for drug discovery: methods and applications. Nat Rev Drug Discov. 2004, 3: 935-949. 10.1038/nrd1549.

    Article  CAS  PubMed  Google Scholar 

  54. 54.

    Adcock SA, McCammon JA: Molecular dynamics: survey of methods for simulating the activity of proteins. Chem Rev. 2006, 106: 1589-1615. 10.1021/cr040426m.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  55. 55.

    Schwab F, van Gunsteren WF, Zagrovic B: Computational study of the mechanism and the relative free energies of binding of anticholesteremic inhibitors to squalene-hopene cyclase. Biochemistry. 2008, 47: 2945-2951. 10.1021/bi702067h.

    Article  CAS  PubMed  Google Scholar 

  56. 56.

    Kim JT, Hamilton AD, Bailey CM, Domaoal RA, Wang L, Anderson KS, Jorgensen WL: FEP-guided selection of bicyclic heterocycles in lead optimization for non-nucleoside inhibitors of HIV-1 reverse transcriptase. J Am Chem Soc. 2006, 128: 15372-15373. 10.1021/ja066472g.

    Article  CAS  PubMed  Google Scholar 

  57. 57.

    Tembe BL, Mccammon JA: Ligand receptor interactions. Comput Chem. 1984, 8: 281-283. 10.1016/0097-8485(84)85020-2.

    Article  CAS  Google Scholar 

  58. 58.

    Chipot C, Pearlman DA: Free energy calculations. The long and winding gilded road. Mol Simulat. 2002, 28: 1-12. 10.1080/08927020211974.

    Article  CAS  Google Scholar 

  59. 59.

    Chipot C, Rozanska X, Dixit SB: Can free energy calculations be fast and accurate at the same time? Binding of low-affinity, non-peptide inhibitors to the SH2 domain of the src protein. J Comput Aided Mol Des. 2005, 19: 765-770. 10.1007/s10822-005-9021-3.

    Article  CAS  PubMed  Google Scholar 

  60. 60.

    Roux B, Deng YQ: Calculation of standard binding free energies: Aromatic molecules in the T4 lysozyme L99A mutant. J Chem Theory Comput. 2006, 2: 1255-1273. 10.1021/ct060037v.

    Article  PubMed  Google Scholar 

  61. 61.

    Roux B, Wang JY, Deng YQ: Absolute binding free energy calculations using molecular dynamics simulations with restraining potentials. Biophys J. 2006, 91: 2798-2814. 10.1529/biophysj.106.084301.

    PubMed Central  Article  PubMed  Google Scholar 

  62. 62.

    Jorgensen WL, Zeevaart JG, Wang LG, Thakur VV, Leung CS, Tirado-Rives J, Bailey CM, Domaoal RA, Anderson KS: Optimization of azoles as anti-human immunodeficiency virus agents guided by free-energy calculations. J Am Chem Soc. 2008, 130: 9492-9499. 10.1021/ja8019214.

    PubMed Central  Article  PubMed  Google Scholar 

  63. 63.

    Ren PY, Jiao D, Zhang JJ, Duke RE, Li GH, Schnieders MJ: Trypsin-ligand binding free energies from explicit and implicit solvent simulations with polarizable potential. J Comput Chem. 2009, 30: 1701-1711. 10.1002/jcc.21268.

    PubMed Central  Article  PubMed  Google Scholar 

  64. 64.

    Ge XX, Roux B: Absolute binding free energy calculations of sparsomycin analogs to the bacterial ribosome. J Phys Chem B. 2010, 114: 9525-9539. 10.1021/jp100579y.

    Article  CAS  PubMed  Google Scholar 

  65. 65.

    Michel J, Essex JW: Hit identification and binding mode predictions by rigorous free energy simulations. J Med Chem. 2008, 51: 6654-6664. 10.1021/jm800524s.

    Article  CAS  PubMed  Google Scholar 

  66. 66.

    McCammon JA: Computer-aided drug discovery: physics-based simulations from the molecular to the cellular level. Physical Biology: From Atoms to Medicine. Edited by: Zewail AH. 2008, London, England: Imperial College Press, 401-410.

    Chapter  Google Scholar 

Download references


This work was carried out with funding from NIH GM31749, NSF MCB-0506593, and MCA93S013 to JAM. Additional support from the Howard Hughes Medical Institute, the NSF Supercomputer Centers, the San Diego Supercomputer Center, the WM Keck Foundation, the National Biomedical Computational Resource, and the Center for Theoretical Biological Physics is gratefully acknowledged. We would also like to thank Paul Gaspar for help with figure preparation.

Author information



Corresponding author

Correspondence to Jacob D Durrant.

Authors’ original submitted files for images

Rights and permissions

This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and Permissions

About this article

Cite this article

Durrant, J.D., McCammon, J.A. Molecular dynamics simulations and drug discovery. BMC Biol 9, 71 (2011).

Download citation


  • molecular dynamics simulations
  • computer-aided drug discovery
  • cryptic binding sites
  • allosteric binding sites
  • virtual screening, free-energy prediction