1
Predictive Power of Biomolecular Simulations

Vojtěch Spiwok

University of Chemistry and Technology, Prague, Department of Biochemistry and Microbiology, Technická 3, 166 28 Prague 6, Czech Republic

Biomolecular simulations are becoming routine in structure‐based drug design and related fields. This chapter briefly presents the history of molecular simulations, basic principles and approximations, and the most common designs of computational experiments. I also discuss statistical analysis of simulation results together with possible limits of accuracy.

The history of computational modeling of molecular structure and dynamics goes back to 1953, to the work of Rosenbluth and coworkers [1]. It introduced the Markov chain Monte Carlo as a method to study a simplified model of the fluid system. Atoms of the studied system were perfectly inelastic and the system was two‐dimensional (2D) instead of three‐dimensional (3D), so the analogy with real molecular systems was not perfect. The first molecular dynamics simulation (i.e. modeling of motions) on the same system was done by Alder and Wainwright in 1957 [2] using perfectly elastic collision between 2D particles. The first molecular simulation with specific atom types was done by Rahman in 1964 [3]. Rahman used a CDC 3600 computer to simulate dynamics of 864 argon atoms modeled using Lennard‐Jones potential. The first simulation of liquid water was published by Rahman and Stillinger in 1971 [4].

Another big milestone was the first biomolecular simulation. McCammon, Gelin, and 2013 Nobel Prize winner Karplus simulated 9.2 ps of the life of the bovine pancreatic trypsin inhibitor (, also known as aprotinin) in vacuum [5]. The simulation was performed during the (Centre Européen de Calcul Atomic et Moléculaire) workshop “Models of Protein Dynamics” in Orsay, France on CECAM computer facilities [6]. It was one of the first works showing proteins as a dynamic species with fluid‐like internal motions, even though in the native state.

Biomolecular simulations have undergone a huge progress in terms of accuracy, size of simulated systems, and simulated times since their pioneer times. However, the question arises whether this progress is enough for their practical application in drug discovery, protein engineering, and related applied fields. To address this issue, let me present here the concept of the hype cycle [7] developed by Gartner Inc. and depicted in Figure 1.1. According to this concept, every new invention starts by a Technology Trigger. Visibility of the invention grows until it reaches the Peak of Inflated Expectations. At this point, failures of the invention start to dominate over its benefits and the invention falls into the phase of Trough of Disillusionment. From this phase a new and slower progress starts in the phase of Slope of Enlightenment toward the Plateau of Productivity. Biomolecular simulation passed the Technology Trigger and Peak of Inflated Expectations as many expected that biomolecular simulation would become routine and an inexpensive alternative to experimental testing of compounds for biological activity. Now, in my opinion, biomolecular simulations are located on the Slope of Enlightenment with a slow but steady progress toward the Plateau of Productivity.

Figure 1.1 Gartner hype cycle of inventions.

1.1 Design of Biomolecular Simulations

Biomolecular simulations can follow different designs. I use the term design to describe the setup of the simulation procedure chosen in order to answer the research hypothesis. There are three major designs of molecular simulation. The first design starts from a predicted structure of the molecular system, which we want to evaluate, for example, a protein–ligand complex predicted by a simple protein–ligand docking. I refer to this as the evaluative design (Figure 1.2). The research hypothesis is: Does the predicted structure represent real structure? The basic assumption behind this design is that an accurately predicted structure of the system, for example, an accurately modeled structure of the complex, is lower in free energy than an inaccurately predicted one. The system therefore tends to be stable in a simulation starting from an accurately modeled structure and tends to be unstable in a simulation starting from an inaccurate structure. The evaluative design can be represented by the study of Cavalli et al. [8]. This study was published in 2004, and simulated times are therefore significantly shorter (typically 2.5 ns) than those available today. Nevertheless, the same length of simulations can be used today with much higher throughput in terms of the number of tested compounds or their binding poses; therefore, the study is still highly actual. Docking of propidium into human acetylcholine esterase (Alzheimer disease target) by the program Dock resulted in the prediction of 36 possible binding poses (clusters of docked binding poses). Six of them were then subjected to 2.5‐ns simulation. Evolution of these systems was analyzed in terms of root‐mean‐square deviation (). Binding poses with high stability in simulations were similar to experimentally determined binding poses for a homologous enzyme.

Figure 1.2 Schematic illustration of designs of biomolecular simulations. Horizontal dimensions correspond to coordinates of the system, and contours correspond to the free energy.

The second design is referred to as refinement design (Figure 1.2). It uses an assumption similar to the evaluative design, i.e. that molecular simulations tend to evolve from high‐free energy states to low‐free energy states. In the refinement design, it is hoped that the dynamics can drive the system from the predicted structure, even though incorrectly predicted, to global free energy minimum, the correct structure, or at least close to it. Naturally, shorter simulation times are necessary to demonstrate correctness or incorrectness of a model by the evaluative design. Longer simulation times are necessary to drive the system from the incorrect to the correct state by the refinement design. In the previous paragraph, I used the study of Cavalli et al. from 2004 [8] as an example of evaluative design. I can present the refinement design on the work published by the same author 11 years later [9]. They used unbiased simulation to predict the binding pose of picomolar inhibitor 4′‐deaza‐1′‐aza‐2′‐deoxy‐1′‐(9‐methylene)‐immucillin‐H in human purine nucleoside phosphorylase. They carried out 14 simulations (500 ns each) of the system containing the trimeric enzyme, 9 ligand molecules (to increase its concentration) placed outside the protein molecule, solvent, and ions. From these simulations, 11 evolved toward binding with a good agreement with the experimentally determined structure of the complex. RMSD from the experimentally determined structure of the complex dropped during these simulations from approximately 6 to 0.2–0.3 nm.

The last design introduced here is referred to as equilibrium design (Figure 1.2). In this design, we hope that the simulation is sufficiently long (or sampling is sufficiently enhanced) to explore all relevant free energy minima and to sample them according to their distribution in the real system. Naturally, the equilibrium design requires longest simulation times or highest sampling enhancement from all three simulation designs. As an example I can present the study by D.E. Shaw Research [10]. The authors simulated systems containing the protein FK506 binding protein () with one of six fragment ligands, water, and ions. They carried out 10‐µs simulations for each ligand. The dissociation constant of a complex can be calculated from its association kinetics as K D = k off/k on. Weak binding (high K D) together with reasonably fast binding kinetics therefore implies that unbinding is also sufficiently fast. For this reason, microsecond timescales were enough to observe multiple binding and unbinding events for millimolar ligands. The fragments identified by these simulations as relatively strong binders can be selected and combined into larger compounds with higher affinity in the manner of fragment‐based drug design [11]. Fragment‐based drug design and molecular dynamics simulation seem to be a good combination. Fragment‐based design requires testing of a low number of weak ligands. This is good, since biomolecular simulations are computationally expensive. Reciprocally, weak binding enables to use molecular dynamics simulations in available timescales. Moreover, unlike some experimental methods of fragment‐based drug design, molecular simulations provide binding pose prediction that can be used to combine fragments.

The three designs described are not without pitfalls. Most of these pitfalls are caused by limitations of simulated timescales. It is often difficult or impossible to simulate timescales long enough to destabilize the structure in the evaluation design, reach the global free energy minimum in the refinement design, or obtain the equilibrium...