5.1 MD as a technique to study atomic motions

A spectacular example of success reflecting a modern marriage among modeling theory, technology and biology involves the theory and practice of MD. Essentially, MD is statistical mechanics by numbers, or Laplace's vision of Newtonian physics (de Laplace, Reference de Laplace1820) on modern supercomputers. The impressive progress in the development of biomolecular force fields, coupled to spectacular computer technology advances, has now made it possible to transform this vision into a reality, by overcoming the difficulty noted by Dirac of solving the equations of motion for multi-body systems (Dirac, Reference Dirac1929).

Since successful applications were reported in the 1970s in protein dynamics (McCammon et al. Reference McCammon, Gelin and Karplus1977), MD – with variations and extensions – has now become a popular (Fig. 2d) and universal tool, ‘as if it were the differential calculus’ (Maddox, Reference Maddox1989). MD is in fact one of the few tools available, by both experiment and theory, to probe molecular motion on the atomic scale. By following the equations of motion as dictated by a classic molecular mechanics force field, complex relationships among biomolecular structure, flexibility, and function can be investigated.

Today's sophisticated dynamics programs, like NAMD or GROMACS, adapted to parallel and massively parallel computer architectures, and specialized hardware have made simulations of biomolecular systems in the microseconds feasible in several weeks of computing. Anton's hardware/software co-design is pushing the envelope to the millisecond timeframe (Shaw et al. Reference Dror, Arlow, Borhani, Jensen, Piana and Shaw2009, Reference Shaw, Maragakis, Lindorff-Larsen, Piana, Dror, Eastwood, Bank, Jumper, Salmon, Shan and Wriggers2010). Although the well-recognized limitations of sampling in atomistic dynamics, as well as in the governing force fields, have led to many innovative sampling alternatives to enhance coverage of the thermally accessible conformational space (as recently reviewed in Earl & Deem, Reference Earl and Deem2008; Klein & Shinoda, Reference Klein and Shinoda2008; Lei & Duan, Reference Lei and Duan2007; Liwo et al. Reference Liwo, Czaplewski, Oldziej and Scheraga2008; Schlick, Reference Schlick2009b, Reference Schlickc), many still rely on MD for local sampling.

Overall, MD simulations have been used for numerous applications, including to refine experimental structures, extending most recently to low-resolution crystal data (Schröder et al. Reference Schröder, Levitt and Brunger2010), interpret various experimental data such as single-molecule force-extension curves (e.g. Lee et al. Reference Lee, Hsin, Sotomayor, Comellas and Schulten2009) or NMR spin-relaxation in proteins (e.g. Case, Reference Case2002; Henzler-Wildman et al. Reference Henzler-Wildman, Thai, Lei, Ott, Wolf-Watz, Fenn, Pozharski, Wilson, Petsko and Karplus2007; Tsui et al. Reference Tsui, Radhakrishnan, Wright and Case2000), improve structure-based function predictions, for example, by predicting calcium-binding sites (Altman et al. Reference Altman, Radmer and Glazer2009), link static experimental structures to implied pathways (e.g. Golosov et al. Reference Golosov, Warren, Beese and Karplus2010; Radhakrishnan & Schlick, Reference Radhakrishnan and Schlick2004), estimate the importance of quantum effects in lowering free-energy barriers of biomolecular reactions (Hu et al. Reference Hu, Elstner and Hermans2003), make structural predictions, deduce reaction mechanisms, propose free-energy pathways and associated mechanisms (e.g. Faraldo-Gomez & Roux, Reference Faraldo-Gomez and Roux2007; Radhakrishnan & Schlick, Reference Radhakrishnan and Schlick2005), resolve or shed light on experimental ambiguities (see Subsection 5.6), and design new folds and compounds, including drugs and enzymes (e.g. Baker et al. Reference Baker, Kuhlman, Dantas, Ireton, Varani and Stoddard2003; Hornak & Simmerling, Reference Hornak and Simmerling2007; Jiang et al. Reference Jiang, Althoff, Clemente, Doyle, Röthlisberger, Zanghellini, Gallaher, Betker, Tanaka, Barbas, Hilvert, Houk, Stoddard and Baker2008; Neidle et al. Reference Neidle, Read, Harrison, Romagnoli, Tanious, Gowan, Reszka, Wilson and Kelland2001). Challenging applications to complex systems like membranes, to probe associated structures, motions and interactions (e.g. Grossfield et al. Reference Grossfield, Pitman, Feller, Soubias and Gawrisch2008; Khelashvili et al. Reference Khelashvili, Grossfield, Feller, Pitman and Weinstein2009; Vasquez et al. Reference Vasquez, Sotomayor, Cordero-Morales, Schulten and Perozo2008) are also noteworthy. Specific success stories from such MD simulations are discussed separately in the subsections below.

5.2 Establishment of reliable algorithms

Rigorous algorithm development and analysis provide a clear framework by which to define the important notions of accuracy, stability and error control. This understanding in turn provides practical guides to parameter selection and trajectory analysis. Three notable areas of algorithmic contributions in recent years include methods for MD integration, evaluation of non-bonded electrostatic interactions and enhanced sampling methods for biomolecules.

The establishment of symplectic integrators such as leap frog, velocity Verlet, and constrained dynamics formulations (e.g. Leimkuhler & Reich, Reference Leimkuhler and Reich2004; Schlick, Reference Schlick2010) has allowed researchers to correctly generate MD trajectories and analyze the stability of a simulation in terms of energy conservation, and the robustness of the simulation with respect to the timestep size.

Highlighting resonance artifacts in MD simulations (Mandziuk & Schlick, Reference Mandziuk and Schlick1995; Schlick et al. Reference Schlick, Mandziuk, Skeel and Srinivas1998), predicting resonant timesteps, and establishing stochastic solution to resonances (Barth & Schlick, Reference Barth and Schlick1998; Schlick et al. Reference Schlick, Barth and Mandziuk1997; Sweet et al. Reference Sweet, Petrine, Pande and Izaguirre2008) have all led to an improved understanding and quality of MD simulations (Morrone et al. Reference Morrone, Zhou and Berne2010).

The advent of efficient and straightforward-to-program particle mesh Ewald (PME) (Darden et al. Reference Darden, York and Pedersen1993; Essmann et al. Reference Essmann, Perera, Berkowitz, Darden, Lee and Pedersen1995; York & Yang, Reference York and Yang1994) and related methods (Duan & Krasny, Reference Duan and Krasny2000; Greengard & Rokhlin, Reference Greengard and Rokhlin1987; Saito, Reference Saito1992, Reference Greengard and Rokhlin1997; Sandak, Reference Sandak2001; Skeel et al. Reference Skeel, Tezcan and Hardy2002) for evaluation of the long-range electrostatic interactions, which constitute the most time-consuming part of a biomolecular simulation, has made possible more realistic MD simulations without non-bonded cutoffs. A problem that in part remains unsolved in their implementation involves the optimal integration of PME methods with multiple-timestep methods and parallelization of PME implementations. The presence of fast terms in the reciprocal Ewald component limits the outer timestep and hence the speedup (Morrone et al. Reference Morrone, Zhou and Berne2010; Procacci et al. Reference Procacci, Marchi and Martyna1998; Qian & Schlick, Reference Qian and Schlick2002; Stuart et al. Reference Stuart, Zhou and Berne1996; Zhou et al. Reference Zhou, Harder, Xu and Berne2001). Moreover, memory requirements create a bottleneck in typical PME implementations in MD simulations longer than a microsecond. This is because the contribution of the long-range electrostatic forces imposes a global data dependency on all the system charges; in practice, this implies communication problems (Fitch et al. Reference Fitch, Rayshubskiy, Eleftheriou, Ward, Giampapa, Pitman and Germain2006). Thus, much work goes into optimizing associated mesh sizes, precision, and sizes of the real and inverse spaces to delay the communication bottleneck as possible (e.g. Shaw et al. Reference Dror, Arlow, Borhani, Jensen, Piana and Shaw2009), but overall errors in long simulations are far from trivial (Morrone et al. Reference Morrone, Zhou and Berne2010; Snir, Reference Snir2004).

In addition to MD integration and electrostatic calculations, sampling the vast configurational space has also triggered many innovative approaches to capture ‘rare events’. The many innovative enhanced sampling methods are either independent of MD or based on MD. In the former class, as recently surveyed (Earl & Deem, Reference Earl and Deem2008; Liwo et al. Reference Liwo, Czaplewski, Oldziej and Scheraga2008; Schlick, Reference Schlick2009b), are various Monte Carlo approaches, harmonic approximations, and coarse-grained models. These can yield valuable conformational insights into biomolecular structure and flexibility, despite altered kinetics. Although Monte Carlo methods are not always satisfactory for large systems on their own right, they form essential components of more sophisticated methods (e.g. transition path sampling (Dellago & Bolhuis, Reference Dellago and Bolhuis2007) or Markov chain Monte Carlo sampling (Pan & Roux, Reference Pan and Roux2008)), as recently surveyed (Schlick, Reference Schlick2009c).

More generally, MD-based methods for enhanced sampling of biomolecules can involve modification of the potential (like accelerated MD or AMD (Grant et al. Reference Grant, Gorfe and McCammon2010; Hamelberg et al. Reference Hamelberg, Mongan and McCammon2004), simulation protocol (like replica-exchange MD or REMD (Sugita & Okamoto, Reference Sugita and Okamoto1999)), or algorithm; as well as global reformulations, such as transition path sampling (Bolhuis et al. Reference Bolhuis, Chandler, Dellago and Geissler2002; Dellago & Bolhuis, Reference Dellago and Bolhuis2007), forward flux simulation (Borrero & Escobedo, Reference Borrero and Escobedo2008), and Markov state models (Noé & Fischer, Reference Noé and Fischer2008). Simple modifications like AMD can aid in the interpretation of residual dipolar coupling and chemical shift measurements in NMR (Cervantes et al. Reference Cervantes, Markwick, Sue, McCammon, Dyson and Komives2009; Markwick et al. Reference Markwick, Bouvignies, Salmon, McCammon, Nilges and Blackledge2009, Reference Markwick, Cervantes, Abel, Komives, Blackledge and McCammon2010). Global formulations, however, are needed more generally not only to generate more configurations or to suggest mechanistic pathways but also to compute free-energy profiles for the reaction and detailed kinetics including reaction rates.

Divide and conquer methods for sampling and for piecing together reaction rate information are especially suitable for readily available computer cluster networks. REMD is a popular method, with many extensions (e.g. Roitberg et al. Reference Roitberg, Okur and Simmerling2007), but implementations for large systems are problematic due the need for many replicas and the non-uniform sampling. There are many successful reports of using tailored enhanced sampling methods for biomolecular applications (e.g. Abrams & Vanden-Eijnden, Reference Abrams and Vanden-Eijnden2010; Berezhkovskii et al. Reference Berezhkovskii, Hummer and Szabo2009; Chennamsetty et al. Reference Chennamsetty, Voynov, Kayser, Helk and Trout2009; Noé et al. Reference Noé, Horenko, Schütte and Smith2007, Reference Noé, Schutte, Vanden-Eijnden, Reich and Weikl2009; Ozkan et al. Reference Ozkan, Wu, Chodera and Dill2007; Radhakrishnan & Schlick, Reference Radhakrishnan and Schlick2004; Voelz et al. Reference Voelz, Bowman, Beauchamp and Pande2010), but applications at large to biomolecules and/or in the absence of experimental data remain a challenge.

Among the numerous other examples of algorithms, PCA, clustering techniques and simple consistency tests have helped to analyze complex groups of experimental data. For example, clustering graphs of RNA secondary structures led to a proposal of novel RNA motifs that are RNA-like (i.e. physical) and non-physical (Kim et al. Reference Kim, Chivian and Baker2004a, Reference Kim, Shiffeldrim, Gan and Schlickb); some in the former class have been designed computationally and also later experimentally confirmed (Kim et al. Reference Kim, Izzo and Schlick2010; Schlick, Reference Schlick2010). Interestingly, a simple tool to analyze PDB structures in terms of quantifying the interior protein volume called RosettaHoles (Sheffler & Baker, Reference Sheffler and Baker2008) detected a group of significant outliers, some of which were traced to a single author; this eventually led to the retraction of structures (Borrell, Reference Borrell2009).

Thus, algorithms are playing an increasing role in a wide range of problems dealing with biomolecular structure, thermodynamics, and dynamics, as well as design.

5.3 Structure predictions that preceded experiment

Although the general problem of structure prediction is far from solved, predicting a peptide structure has been more amenable to present modeling, and there have also been notable success stories for larger systems, as early as the late 1980s, for human immunodeficiency virus (HIV) protease (PR), and for complex RNAs, by sequence alignment. The successes for RNA are especially noteworthy since far fewer structures of RNAs are available at atomic resolution and the field is younger than protein prediction. Some interesting early prediction cases are collected below, following general trends for peptides and proteins.

5.3.1 HIV-1 PR

HIV-1 PR is an aspartic PR essential for the replication of the HIV. Already in 1987, Pearl & Taylor attempted to predict the HIV-1 PR tertiary structure by computational modeling alone (Fig. 4a) (Pearl & Taylor, Reference Pearl and Taylor1987). Their approach involved multiple sequence alignment to predict a secondary structure followed by a knowledge-based method that matches fragments from large related structures to produce a tertiary structure model.

Fig. 4.

Examples of folding predictions that preceded experiment. (a) HIV PR: tertiary structure of the homodimer (red), predicted by homology modeling (Pearl & Taylor, Reference Pearl and Taylor1987), superimposed with the crystal structure 3HVP (blue), 198 residues shown (Wlodawer et al. Reference Wlodawer, Miller, Jaskolski, Sathyanarayana, Baldwin, Weber, Selk, Clawson, Schneider and Kent1989); (b) group-I intron: P4–P6 domain (red), predicted by a comparative modeling (Michel & Westhof, Reference Michel and Westhof1990), aligned with crystal structure 1GID (blue), 38 residues (Cate et al. Reference Cate, Gooding, Podell, Zhou, Golden, Kundrot, Cech and Doudna1996); (c) CASP protein prediction by the Baker group for target T0492 (Raman et al. Reference Raman, Vernon, Thompson, Tyka, Sadreyev, Pei, Kim, Kellogg, DiMaio, Lange, Kinch, Sheffler, Kim, Das, Grishin and Baker2009), kindly provided by Srivatsan Raman and David Baker; crystal structure, 73 residues (blue) of this all-β protein (SH3-like barrel fold) is shown against the best detectable template found by comparative modeling (green) and the predicted structure (red) built using Rosetta using the best template found.

Two years later, another theoretical model was produced (Weber et al. Reference Weber, Miller, Jaskolski, Leis, Skalka and Wlodawer1989) using related information from the Rous sarcoma virus (RSV) determined at 3 Å resolution (Miller et al. Reference Miller, Jaskolski, Rao, Leis and Wlodawer1989). In the same year, the HIV-1 PR structure was resolved by X-ray crystallography at 3 Å resolution (Navia et al. Reference Navia, Fitzgerald, McKeever, Leu, Heimbach, Herber, Sigal, Darke and Springer1989). Interestingly, the authors noted a local structural discrepancy with respect to the predicted models (Pearl & Taylor, Reference Pearl and Taylor1987; Weber et al. Reference Weber, Miller, Jaskolski, Leis, Skalka and Wlodawer1989) in the topology of the N- and C-terminal regions and the helix near the C-terminal; some of these residues at the N-terminal were not observed in the electron density map, and so inherent disorder and flexibility were presented as possible explanations for these differences. However, when another structure of HIV-1 PR at 2·8 Å resolution was solved soon after, residues not captured in the first crystal were explicitly resolved (Wlodawer et al. Reference Wlodawer, Miller, Jaskolski, Sathyanarayana, Baldwin, Weber, Selk, Clawson, Schneider and Kent1989), and subsequent analysis showed excellent agreement with the modeling predictions rather than the previous crystal. The studies of HIV-1 PR initiated by the theoretical models and confirmed later by the crystallography eventually led to development, largely due to molecular modeling, of the first PR inhibitor, saquinavir, which was approved by FDA in 1995, and nine other PR inhibitors approved by 2006.

5.4 Protein-folding theory

Besides the folded protein structures discussed above, other parts of the protein-folding puzzle that have intrigued experimentalists and theoreticians alike are understanding the associated protein-folding code (thermodynamics) and pathways (kinetics) (Dill et al. Reference Dill, Ozkan, Shell and Weikl2008). A statistical–mechanics framework based on the density of states has been pioneered by Frauenfelder, Wolynes, Dill, Onuchic, Thirumalai and others for answering these questions (e.g. Dill et al. Reference Dill, Ozkan, Shell and Weikl2008; Frauenfelder et al. Reference Frauenfelder, Sligar and Wolynes1991; Wolynes, Reference Wolynes2005) by analyzing energy landscapes of proteins on the basis of various experimental studies of proteins (fluorescence spectroscopy, NMR, single-molecule experiments, fast-laser temperature jumps, mutational studies, etc.). Essentially, the density of states arranged in the funnel landscape describes the conformational heterogeneity of the protein. Such free-energy landscapes can help interpret conformational substates and folding mechanisms and thereby relate protein dynamics, folding kinetics, and function. For example, fast protein folders must navigate through mostly local contacts in smoother free-energy landscapes, while slow-folding kinetics can be attributed to rugged energy landscapes, with more global contacts and more conformational traps and barriers. Standard chemical kinetics and transition state theory in combination with simplified models have also been important in analyzing folding events (e.g. Berezhkovskii et al. Reference Berezhkovskii, Hummer and Szabo2009; Best & Hummer, Reference Best and Hummer2006; Ozkan et al. Reference Ozkan, Wu, Chodera and Dill2007; Yang & Grueble, Reference Yang and Grueble2003).

By relating folding kinetics to mutation profiles, mutation experiments have also helped further validate the utility of such theoretical frameworks. For example, recent folding kinetics experiments on mutants of α-spectrin, a protein of the intracellular matrix of red blood cells important for membrane elasticity (Wensley et al. Reference Wensley, Batey, Bone, Chan, Tumelty, Steward, Kwa, Borgia and Clarke2010), have confirmed the notion that slow kinetics is correlated to rugged landscapes. By swapping domains through chimeric constructs for protein domains of α-spectrin, the reasons for disparate folding rates (by three orders of magnitude) of the protein relative to two homologues were traced to internal friction resulting from residue-specific interactions that can lead to misdocking of helices. For this membrane protein, slow unfolding kinetics can be advantageous because fewer rearrangements during the cell's lifetime, and hence fewer possible degradations, result.

That theory contributes significantly to experiment was also echoed by Fersht and colleagues in a separate work on a three-helix bundle protein and various mutants which showed via simulation that folding intermediates are on the folding pathway (Mayor et al. Reference Mayor, Guydosh, Johnson, Grossmann, Sato, Jas, Freund, Alonso, Daggett and Fersht2003):

Numerous other general areas of theory, such as protein electrostatics and enzymatic reactions (Kamerlin et al. Reference Kamerlin, Haranczyk and Warshel2009; Warshel et al. Reference Warshel, Sharma, Kato, Xiang, Liu and Olson2006), conformational transitions in DNA polymerases in relation to fidelity mechanisms (Bebenek et al. Reference Bebenek, Garcia-Diaz, Foley, Pedersen, Schlick and Kunkel2008; Foley & Schlick, Reference Foley and Schlick2009; Golosov et al. Reference Golosov, Warren, Beese and Karplus2010; Radhakrishnan et al. Reference Radhakrishnan, Arora, Wang, Beard, Wilson and Schlick2006), and interpretation of NMR signals (Case, Reference Case2002; Henzler-Wildman et al. Reference Henzler-Wildman, Thai, Lei, Ott, Wolf-Watz, Fenn, Pozharski, Wilson, Petsko and Karplus2007; Markwick et al. Reference Markwick, Cervantes, Abel, Komives, Blackledge and McCammon2010; Tsui et al. Reference Tsui, Radhakrishnan, Wright and Case2000) have also benefited from simplified theoretical frameworks and models developed to interpret complex biological phenomena.

5.5 Modeling-aided drug discovery and design

Modeling molecular structures and dynamics has already helped to define molecular specificity and clarify functional aspects that are important for drug development (see also Schwede et al. Reference Schwede, Sali, Honig, Levitt, Berman, Jones, Brenner, Burley, Das, Dokholyan, Dunbrack, Fidelis, Fiser, Godzik, Huang, Humblet, Jacobson, Joachimiak, Krystek, Kortemme, Kryshtafovych, Montelione, Moult, Murray, Sanchez, Sosnick, Standley, Stouch, Vajda, Vasquez, Westbrook and Wilson2009). Although modeling suggestions are not always verified by experiment, they help make important suggestions and offer general leads.

Clearly, as biological structures and functions are being described, natural disease targets that affect the course of disease can be proposed. Examples include specific regions of HIV enzymes or misfolded protein fibrils known as amyloids (for related misfolding disorders like Alzheimer's and Creutzfeltdt–Jakob disease). Modeling work is important for pinpointing specific targets in these regions from structural, stability, and flexibility considerations and also testing resulting interactions (as reviewed in Ghosh et al. Reference Ghosh, Nie, An and Huang2006; Zhou, Reference Zhou2008). For example, the sequence-specific triggers for amyloid formation, as recently uncovered (Goldschmidt et al. Reference Goldschmidt, Teng, Riek and Eisenberg2010) – self-complementary segments that can stick together – suggest capping alterations to slow down their formation or small molecules to act as folding chaperones (Schnabel, Reference Schnabel2010). Supplementing these structure–function drug design paradigms is a systems-biology approach that attempts to modify response of genes, proteins and metabolites by integrating organ- and system-level modeling (Csermely et al. Reference Csermely, Agoston and Pongor2005; Kitano, Reference Kitano2007).

The solved HIV enzymes (PR, integrase, and reverse transcriptase) led to the development of many inhibitors that are now in clinical usage, some of which were aided by computer simulations. For example, besides the structural determination mentioned above for HIV PR, MD simulations using various HIV PR models have suggested that a fully open conformation of the PR ‘flaps’ may be favorable for drug access to the active site (Collins et al. Reference Collins, Burt and Erickson1995; Hamelberg & McCammon, Reference Hamelberg and McCammon2005; Hornak et al. Reference Hornak, Okur, Rizzo and Simmerling2006; Tozzini & McCammon, Reference Tozzini and McCammon2005). Recent simulations have also led to design proposals (Scott & Schiffer, Reference Scott and Schiffer2000) and other insights into PR/drug interactions (Hornak & Simmerling, Reference Hornak and Simmerling2007).

The solved complex of HIV integrase with an inhibitor in 1999 (Goldgur et al. Reference Goldgur, Craigie, Cohen, Fujiwara, Yoshinaga, Fujishita, Sugimoto, Endo, Murai and Davies1999) attracted considerable attention since it provided a general platform for the drug design of another class of HIV inhibitors besides PR inhibitors. MD simulations have suggested that inhibitors could bind in more than one orientation to HIV-1 integrase, preferred and flipped (Perryman et al. Reference Perryman, Forli, Morris, Burt, Cheng, Palmer, Whitby, McCammon, Phillips and Olson2010; Schames et al. Reference Schames, Henchman, Siegel, Sotriffer, Ni and McCammon2004), as shown in Fig. 5, which also shows an X-ray structure of an inhibitor bound to an evolutionary similar virus (Hare et al. Reference Hare, Gupta, Valkov, Engelman and Cherepanov2010). Simulations also showed that binding modes can be selected to exploit stronger interactions in specific regions and orientations (Hazuda et al. Reference Hazuda, Anthony, Gomez, Jolly, Wai, Zhuang, Fisher, Embrey, Guare, Egbertson, Vacca, Huff, Felock, Witmer, Stillmock, Danovich, Grobler, Miller, Espeseth, Jin, Chen, Lin, Kassahun, Ellis, Wong, Xu, Pearson, Schleif, Cortese, Emini, Summa, Holloway and Young2004; Lin et al. Reference Lin, Perryman, Schames and McCammon2002; Schames et al. Reference Schames, Henchman, Siegel, Sotriffer, Ni and McCammon2004) and that different divalent-ion arrangements are associated with these binding sites and fluctuations (Perryman et al. Reference Perryman, Forli, Morris, Burt, Cheng, Palmer, Whitby, McCammon, Phillips and Olson2010). Thus, modeling can help outline new possibilities that stimulate further modeling and experimentation. More generally, emerging experimental findings combined with modeling are expected to aid the development of antiretroviral drugs.

Fig. 5.

Drug/inhibitor-binding pockets to HIV integrase as solved experimentally and as predicted by modeling. (a) Preferred (or primary) binding (yellow shading) of the inhibitor 5-CITEP (green) to the active site of the HIV-1 integrase core domain as determined by crystallography (Goldgur et al. Reference Goldgur, Craigie, Cohen, Fujiwara, Yoshinaga, Fujishita, Sugimoto, Endo, Murai and Davies1999); (b) preferred binding (yellow shading) of the anti-HIV drug Raltegravir (purple) to the active site of the prototype foamy virus (PFV) integrase, a structural homologue of HIV-1 integrase, as determined by crystallography (Hare et al. Reference Hare, Gupta, Valkov, Engelman and Cherepanov2010); (c) predicted preferred binding (yellow shading) for Raltegravir to the HIV-1 integrase catalytic core domain as determined by modeling (Perryman et al. Reference Perryman, Forli, Morris, Burt, Cheng, Palmer, Whitby, McCammon, Phillips and Olson2010); (d) alternative, or flipped binding (turquoise shading) for Raltegravir to the HIV-1 integrase catalytic core domain as determined by modeling (Perryman et al. Reference Perryman, Forli, Morris, Burt, Cheng, Palmer, Whitby, McCammon, Phillips and Olson2010). C and D are based on PDB entry 1QS4, with missing residues in panel (A) reconstructed by homology. The magnesium ions are shown as orange spheres.

Other examples of drug successes based on molecular modeling include the design of potent thrombin inhibitors to treat a variety of blood coagulation and clotting-related diseases by tailoring ligands to the enzyme pocket (Brady et al. Reference Brady, Stauffer, Lumma, Smith, Ramjit, Lewis, Lucas, Gardell, Lyle, Appleby, Cook, Holahan, Stranieri, Lynch, Lin, Chen, Vastag, Naylor-Olsen and Vacca1998), the SARS virus inhibitor (Dooley et al. Reference Dooley, Shindo, Taggart, Park and Pang2006), glutamate nanosensors to monitor neurologic functions (Okumoto et al. Reference Okumoto, Looger, Micheva, Reimer, Smith and Frommer2005), agonists to treat anxiety and depression (Becker et al. Reference Becker, Dhanoa, Marantz, Chen, Shacham, Cheruku, Heifetz, Mohanty, Fichman and Sharadendu2006), and a migraine treatment (Boyd, Reference Boyd1998).

While both the monetary cost and development time required for each successful drug remains very high (Amir-Aslani, Reference Amir-Aslani2008; Munos, Reference Munos2009), and great successes are now few and far between, it is clear that molecular modeling will be important for the future of drug discovery and for many other design applications in biomedicine and nanotechnology. What is not obvious, however, is what type of modeling tool will prove to be most useful. This is because the usage of MD simulations by pharmaceutical scientists appears to be minimal (Nicholls, Reference Nicholls2010), compared to tools from cheminformatics, ligand similarity searches and visualization, which are heavily used, and quantitative structure-activity relationships (QSAR), docking, homology modeling, bioinformatics and quantum mechanics, to a lesser extent. Limitations of force fields and sampling are most often quoted as reasons for the curbed application of dynamics to real-life problems.

The emergence of another class of successful drugs distinct from the traditional small-molecule drugs, namely biologics – biomolecules derived from living cells – might also increase the role of modeling in pharmaceutical discovery. This is because understanding the complex biomolecular function is important to this alternative approach, which also relies on a systems biology viewpoint. Overall, it is becoming clear that integrated advances on several fronts – small-molecule drug design, biologics, high-throughput approaches, pharmacogenomics and other innovations – are required to develop new paradigms to manage the complex scientific, technological and economic factors involved in drug design.

5.6 Interpretation of experimental structures and related data disparities

In a recent review of their group's applications of biomolecular dynamics for interpretation of atomic force microscopy studies and electron microscopy structures, Schulten and colleagues write (Lee et al. Reference Lee, Hsin, Sotomayor, Comellas and Schulten2009) that:

They continue to argue that:

Other examples were provided above in Subsection 5.1, on using MD simulations and other modeling techniques (spectral decompositions, enhanced sampling methods, QM/MM techniques, etc.) to gain insights into biological structure and function. Many of these predictions can ultimately be tested by further experimentation.

Importantly, modeling can help sift through conflicting experimental information. Two such examples where different experimental techniques have produced different structural perspectives involve the organization of G quadruplexes and of chromatin fibers.

5.6.1 G-quadruplexes

G-quadruplexes are organized topologies of G-rich sequences in telomeric DNA (chromosome ends) held by Hoogsteen base pairing. Such sequences are of great interest because of their role in preserving genomic stability; a better understanding of their structural properties could help interpret processes related to human aging caused by chromosome degradation.

Crystal structures of short telomeric sequences are available, as are several NMR solution structures. The topologies from these studies, however, differ in the orientation of the loops connecting the G-tetrad planes (Fig. 6) and in the overall stacking arrangement of the complexes. Namely, whereas the crystal topologies suggest parallel stranded constructs, the NMR structures suggest anti-parallel backbone arrangements with stabilizing ions (Burge et al. Reference Burge, Parkinson, Hazel, Todd and Neidle2006).

Fig. 6.

G quadruplex architecture and drug binding. (a) Antiparallel and (b) parallel folds of human telomeric G-quadruplexes; (c) biological unit of the complex between the drug candidate molecule BRACO-19 and two bimolecular human telomeric quadruplexes (Campbell et al. Reference Campbell, Parkinson, Reszka and Neidle2008). Figure kindly provided by Stephen Neidle.

By modeling guanine-rich quadruplexes using MD, free-energy simulations, and PCA, Neidle and co-workers determined stable structures of these complexes in parallel orientations and predicted higher-order organizations (Haider & Neidle, Reference Haider and Neidle2009; Haider et al. Reference Haider, Parkinson and Neidle2008). These structures help propose plausible models on longer systems that could not be solved experimentally as well as build higher-order repeat models. Such detailed models can examine the relative stability of parallel and anti-parallel arrangements and suggest conditions under which each form will be stable. Ultimately, atomic models could be used for ligand-binding drug design, already an active area to fight cancer (Neidle & Parkinson, Reference Neidle and Parkinson2002). Figure 6c shows a crystal structure of an anti-cancer drug, at the interface of two parallel-folded quadruplexes and stabilized by a water network, bound between two human telomeric quadruplexes (Campbell et al. Reference Campbell, Parkinson, Reszka and Neidle2008). The structure confirms and extends modeling predictions and underscores that, while theory provides generally correct ideas, experimental data clarify and refine the modeling suggestions.

5.6.2 Chromatin organization

Chromatin structure represents another higher-order structure of great interest. Understanding chromatin organization is important because the chromatin structure affects genome accessibility to the protein synthesis machinery. Hence, transcription regulation is ultimately a structural problem.

Although several levels of folding are recognized for the chromatin fiber, only the basic unit of the chromatin fiber – the nucleosome, or the protein octamer core around which about 200 bp of DNA are wound (Davey et al. Reference Davey, Sargent, Luger, Mäeder and Richmond2002; Kornberg & Thomas, Reference Kornberg and Thomas1974; Luger et al. Reference Luger, Mäder, Richmond, Sargent and Richmond1997) – and the associated low-level packaging are well characterized (Felsenfeld & Groudine, Reference Felsenfeld and Groudine2003).

At physiological ionic strengths, it is believed that the chromatin forms a ‘30 nm’ chromatin filament observed by electron microscopy, but several models for this polynucleosomal level of folding have been suggested based on X-ray crystallography and electron microscopy imaging, with yet no consensus (Tremethick, Reference Tremethick2007; van Holde & Zlatanova, Reference van Holde and Zlatanova2007). For example, the two-start zigzag model, where the nucleosomes crisscross one another around the helical axis with straight linker DNA (Dorigo et al. Reference Dorigo, Schalch, Kulangara, Duda, Schroeder and Richmond2004) is evident in the recent crystallographic structure of the tetranucleosome (Schalch et al. Reference Schalch, Duda, Sargent and Richmond2005), while the solenoid topology (Finch et al. Reference Finch, Lutter, Rhodes, Brown, Rushton, Levitt and Klug1977), in which the nucleosomes arrange helically around the fiber axis with bent linker DNA, is supported by electron micrographs (Fig. 7).

Fig. 7.

Chromatin organization: ideal models and simulation-generated model. (a) Ideal solenoid model for the chromatin fiber (side, top, and upper top layer views), in which DNA linkers are bent and neighboring nucleosomes (i±1) are in closest contact; (b) ideal zigzag model for the chromatin fiber (side, top, and upper top-layer views), in which DNA linkers are straight and next-nearest neighbors (i±2) are in closest contact; (c) heteromorphic architecture obtained by modeling (Grigoryev et al. Reference Grigoryev, Arya, Correll, Woodcock and Schlick2009) at divalent ion environments with the chromatin model shown in Fig. 1c in which mostly zigzag forms with straight linker DNAs are interspersed with bent DNA linkers. In all views, linker and wrapped DNA are colored red; odd and even nucleosomes are white and blue, respectively; and linker histones are turquoise.

Only recently have researchers begun to dissect the influence of key internal and external factors such as length of the connecting linker DNA segments between nucleosomes (which can vary from 10 to 70 bp), the binding of linker histones, and the presence of various concentrations of monovalent and divalent ions on chromatin structure (e.g. Kruithof et al. Reference Kruithof, Chen, Routh, Logie, Rhodes and van Noort2009; Robinson et al. Reference Robinson, Fairall, Huynh and Rhodes2006). Modeling has played an important role (e.g. Arya & Schlick, Reference Arya and Schlick2006; Grigoryev et al. Reference Grigoryev, Arya, Correll, Woodcock and Schlick2009; Schlick & Perisic, Reference Schlick and Perisic2009; Stehr et al. Reference Stehr, Kepper, Rippe and Wedemann2008; Sun et al. Reference Sun, Zhang and Schlick2005; Wong et al. Reference Wong, Victor and Mozziconacci2007) in examining different models, dissecting the dependence of fiber width on the linker DNA length, the effect of linker histones of fiber structure and the condensation roles of divalent ions on chromatin organization. For example, our Monte Carlo simulations of a coarse-grained oligonucleosome model have suggested a compact zigzag organization for the chromatin fiber at typical linker DNA lengths with linker histones and a more heteromorphic architecture of zigzag forms with straight linker DNA interspersed with bent DNA linkers at divalent ion environments (Fig. 7) (Grigoryev et al. Reference Grigoryev, Arya, Correll, Woodcock and Schlick2009); this view thus merges the two ideal models and provides further explanations to the prevalence of both structures by experiments. While studies are ongoing, modeling has underscored the notion that the chromatin structure is heteromorphic (Grigoryev et al. Reference Grigoryev, Arya, Correll, Woodcock and Schlick2009; Wong et al. Reference Wong, Victor and Mozziconacci2007) rather than a single dominant structure.

Interestingly, recent folding simulations of proteins have indicated that conformational space may be more heterogeneous than originally believed (Ensign & Pande, Reference Ensign and Pande2009; Freddolino & Schulten, Reference Freddolino and Schulten2009).

6.1 CASP failures in structure prediction

As mentioned above, prediction of successes for small proteins can be documented by the CASP experiments, especially from template-based modeling; indeed, proteins with >35% sequence similarity tend to be folded alike. However, exceptions are well known (Gan et al. Reference Gan, Perlow, Roy, Ko, Wu, Huang, Yan, Nicoletta, Vafai, Sun, Wang, Noah, Pasquali and Schlick2002). In a recent CASP8 experiment, two protein targets (56 residues each) for template-based modeling were released containing >90% sequence identity; the first protein falls in the α fold class, while the second has an α/β topology. Surprisingly, all 159 participant teams except four failed to recognize the two different folds. Instead, most predictors attributed the α/β fold to both targets; the difficulty was traced to steric clashes occurring at the side chain of the non-identical residues (Alexander et al. Reference Alexander, He, Chen, Orban and Bryan2007).

More common are failures in de novo or ab initio folding predictions. In CASP8, several groups submitted predictions for up to 10 targets. Except for one group that performed well, the performance of the others was inconsistent, i.e. some bad and some good predictions (Kryshtafovych et al. Reference Kryshtafovych, Fidelis and Moult2009). One difficult target, which contained five α-helices and four antiparallel β-strands, received no correct predictions (Ben-David et al. Reference Ben-David, Noivirt-Brik, Paz, Prilusky, Sussman and Levy2009). It has also been established that high-resolution modeling (<2 Å) for new folds is far beyond the current capacity of even the best of programs (Zhang, Reference Zhang2008). Moreover, protein structure refinement appears to be quite a difficult task, as evident from performance in CASP8 for this new refinement category (MacCallum et al. Reference MacCallum, Hua, Schnieders, Pande, Jacobson and Dill2009).

6.2 Failures in configurational sampling

Sampling the high-dimensional rugged function of biomolecules to locate all the thermally accessible configurations is a well-recognized bottleneck in MD. Moreover, despite all the clever enhanced sampling methods that have been designed, the problem is far from solved. The sampling problem is also exacerbated by the parameter approximations in force fields and algorithms, and this rigorously requires stochastic frameworks for modeling as well as data analysis. In practice, this means that most MD simulations sample local states that are initial configuration, model, and parameter dependent, as shown for the narrow sampled ranges in dihedral angle distributions in peptides (Sugita & Okamoto, Reference Sugita and Okamoto1999), or the limited configurational space sampled for the ring peptide cyclosporin A, which could not be found in the all trans conformation except at very high temperatures (O'Donohue et al. Reference O'Donohue, Burgess, Treutlein and Walkinshaw1995). Thus, convergence of structures in the global sense can almost never be proven in practice, even with simplified coarse-grained models and enhanced sampling methods.

6.3 Force field biases

Force field problems are well recognized, for example, in predicting β proteins (Ensign & Pande, Reference Ensign and Pande2009; Freddolino & Schulten, Reference Freddolino and Schulten2009; Maupetit et al. Reference Maupetit, Derreumaux and Tufféry2009; Mittal & Best, Reference Mittal and Best2010) and the related over-stabilization of helices (Best et al. Reference Best, Buchete and Hummer2008; Freddolino et al. Reference Freddolino, Park, Roux and Schulten2009; Hornak et al. Reference Hornak, Okur, Rizzo and Simmerling2006; Tanizaki et al. Reference Tanizaki, Clifford, Connelly and Feig2008), possibly due to systematic errors in protein backbone potentials and the approximate treatment of hydrogen bonds. Noted also in implicit solvent models are structural biases (e.g. towards α helices) compared to explicit-water systems (Roe et al. Reference Roe, Okur, Wickstrom, Hornak and Simmerling2007) and limitations in electrostatic treatments due to simple atom-centered charged models (Beachy et al. Reference Beachy, Chasman, Murphy, Halgren and Friesner1997; Chen & Brooks, Reference Chen and Brooks2008; Halgren, Reference Halgren1999); polarizable force fields thus represent a more accurate description of molecular properties (Ponder et al. Reference Ponder, Wu, Ren, Pande, Chodera, Schnieders, Haque, Mobley, Lambrecht, DiStasio, Head-Gordon, Clark, Johnson and Head-Gordon2010).

Although modern force fields perform comparably in MD simulations (Price & Brooks, Reference Price and Brooks2002), the general parameterization problem is difficult for several reasons. These include the enormity of the parameter space, the difficulty in validating models against the varied experimental data, the desire to simplify models, while at the same time extending the biological scope of the modeling, and the practical issue that force field development may not be as scientifically rewarding as are applications. Thus, no force field can accurately reproduce all the complex interactions and properties of real systems, but force fields need not be perfect to be meaningful! Their overall utility is in generating qualitative and quantitative insights into structural, energetic, and dynamic properties of complex systems through systematic studies, especially when trends are compared for a related group of biological systems (like single base-pair variants of DNA or protein mutants). Ultimately, predictions or interpretations can be tested experimentally.

8.1 Bridging Scales

Given the enormous range and complexity of biological temporal and spatial dimensions, work should continue on all levels of biomolecular modeling, from all-atom representations to coarse-grained models on both the classical and quantum-classical levels, despite imperfections and approximations. This is needed to study processes ranging from macromolecular folding (proteins, RNA, and other complexes) to biochemical binding and reaction mechanisms (enzyme catalysis or protein/ligand interactions) to macromolecular pathways (DNA replication and repair or chromatin organization) and up to supramolecular cellular processes like protein-signaling networks. Because for biology the ‘devil is in the details’, simulation efforts are required on all levels of biological complexity rather than focusing on very accurate studies on the finest level. This means that many approaches in addition to all-atom modeling are important to continue, such as coarse-grained modeling and mesoscale models. New methodological developments are needed in all levels (see the next item) so that, ultimately, an integration of all these tools could allow a kind of ‘telescoping’ from one level of resolution to the other to focus on specific aspects of a variety of biological processes.

8.2 Methodology advances

Studying biological systems and problems on the above levels successfully and reliably requires new methods and models to be systematically developed, including force fields, hybrid quantum/molecular mechanics models, enhanced sampling techniques, rigorous coarse graining of multiscale models, and tool integration. As mentioned for force fields (next item), developing models on other scales requires consolidated efforts of concept development, parameterization, and code sharing.

8.3 Force fields

Improved force fields are needed to address the above biological problems, with parameterization and validation conducted systematically and rigorously on both all-atom and coarse-grained models. The ideas needed to improve force fields appear to be known, like polarizable force fields (Ponder et al. Reference Ponder, Wu, Ren, Pande, Chodera, Schnieders, Haque, Mobley, Lambrecht, DiStasio, Head-Gordon, Clark, Johnson and Head-Gordon2010), but efforts need to be consolidated, prioritized, and organized. Special funding for such consolidated efforts could help attract researchers to work on these tasks.

8.4 Open-source codes

Better methods for data sharing are needed, including for force fields, algorithms, and many computer programs associated with research papers. For example, a requirement to deposit computer codes when submitting methodology papers, as done in some computer science journals, could help reviewers test the programs and ultimately allow users in the community to apply the algorithms. This will go a long way toward validating algorithms and identifying the most effective techniques. For enhanced sampling methods, which tend to be home grown, this program sharing could help practitioners apply clever ideas to important systems.

8.5 Better interdisciplinary education

One origin for a poor perception of the field is that many applications of molecular modeling packages as black boxes are not reliable, although the images may be impressive. This problem could be alleviated with better multidisciplinary education, including in modeling, to teach the basics of modeling studies and analysis, to stimulate modeling experimentation, and to cultivate sound scientific practices (such as parameter testing, statistical interpretations, and sensitivity analysis). For example, programs in computational biology or bioinformatics could offer undergraduate and graduate-level courses like bio/molecular modeling, computational chemistry, or modeling tools for experimentalists. Courses should introduce the fundamentals in the fields as well as provide practical experience through simple programming and modeling assignments to introduce students to common caveats and encourage critical thinking and healthy skepticism (e.g. as in Schlick, Reference Schlick2010).

8.6 Research collaborations

Although funding agencies have led to many teams working together, there remain many duplicate research efforts. A better deployment of truly interdisciplinary teams working together to advance specific areas of methodology and application will be beneficial. This could be achieved through large community efforts on several key fronts of algorithms (like MD, QM/MM and sampling), force fields, and specific applications.

8.7 Community assessment efforts

The CASP exercise has been an exemplary success in defining the state-of-the-art in protein folding and helping propel a field forward in an organized, community-wide manner. Similar efforts are needed for nucleic acids, force field development, MD practices, enhanced sampling methods, quantum-classical mechanics simulations, coarse-grained models, implicit solvation, and others.

8.8 Faster processors

Computer technology has certainly impacted the field tremendously, and will continue to do so. Efforts must continue to make available large central computing resources like supercomputing centers to biomolecular modelers. The sharing of an Anton supercomputer with the community is a commendable effort. In addition, because data management becomes an issue for biomolecular data, and queue times can be long when computations are conducted remotely, algorithms and programs that utilize distributed computing platforms and local networks of fast processors are also important.

8.9 Shared resources

Systematic advances in algorithms, force fields, and modeling approaches require establishment of better software repositories and shared resources. Infrastructural support for generating and analyzing such voluminous molecular data requires the development of efficient simulation management tools such as clustering, archiving, comparison, debugging, visualization, and communication protocols.

8.10 Better communication between academe and industry

Greater efforts are needed to bring open dialog and constructive research efforts concerning tools and applications used in academia and those in the pharmaceutical industry. For example, cross-community interactions could help lead to better efforts in establishing needed databases, such as for protein/ligand structures, or developing better approaches for binding free energies. Overall, better modes are needed for transforming basic or translational research efforts in academia to more applied problems in industry and technology. Economic incentives for such collaborations would be needed, and these could come from a combination of industrial, academic and government agencies.

8.11 Reports of failure

Progress in the field could be expedited by a venue that would allow scientists to publish modeling failures. Understanding the problems or failures could go a long way in helping finding the right solutions. Perhaps journal sections dedicated to such publications could inspire scientists to share such information.

Of course, these wide suggestive strokes require many of the details to be filled in. We invite readers to share comments and suggestions. Such concerted efforts can help move the field away from a local trough of disappointment (or a shattering of inflated expectations) toward a more realistic productive phase, where modeling and experiment can be hand-in-hand partners (Fig. 3).