Quantitative Characterization of the Binding and Unbinding of

Article pubs.acs.org/JCTC

Quantitative Characterization of the Binding and Unbinding of Millimolar Drug Fragments with Molecular Dynamics Simulations Albert C. Pan,*,†,§ Huafeng Xu,*,†,§ Timothy Palpant,†,§ and David E. Shaw*,†,‡ †

D. E. Shaw Research, New York, New York 10036, United States Department of Biochemistry and Molecular Biophysics, Columbia University, New York, New York 10032, United States

‡

S Supporting Information *

ABSTRACT: A quantitative characterization of the binding properties of drug fragments to a target protein is an important component of a fragment-based drug discovery program. Fragments typically have a weak binding affinity, however, making it challenging to experimentally characterize key binding properties, including binding sites, poses, and affinities. Direct simulation of the binding equilibrium by molecular dynamics (MD) simulations can provide a computational route to characterize fragment binding, but this approach is so computationally intensive that it has thus far remained relatively unexplored. Here, we perform MD simulations of sufficient length to observe several different fragments spontaneously and repeatedly bind to and unbind from the protein FKBP, allowing the binding affinities, on- and off-rates, and relative occupancies of alternative binding sites and alternative poses within each binding site to be estimated, thereby illustrating the potential of long time scale MD as a quantitative tool for fragment-based drug discovery. The data from the long time scale fragment binding simulations reported here also provide a useful benchmark for testing alternative computational methods aimed at characterizing fragment binding properties. As an example, we calculated binding affinities for the same fragments using a standard free energy perturbation approach and found that the values agreed with those obtained from the fragment binding simulations within statistical error.

D

of the protein−fragment complex and directly estimate the kinetic and equilibrium constants of binding.13 This computationally intensive approach is now becoming practical as a result of the availability of MD software with improved scalability across multiple processors,18−23 and the ability to run simulations on graphics processing units (GPUs)24−28 or on special-purpose supercomputers designed specifically for MD simulations.29 In this work, we simulate the binding of six fragments (Figure 1) to the protein FKBP using multiple equilibrium simulations of sufficient length to observe a substantial number of binding and unbinding events, allowing the calculation of statistically meaningful estimates of kinetic and thermodynamic binding constants (Tables 1 and S1). We also performed an additional simulation of one of the six fragments, 4-diethylamino-2butanone (DAP), with an alternative protonation state. All 6 fragments have experimentally determined affinities for FKBP ranging from 200 μM to 20 mM.30 The simulations contained a single protein and a single ligand represented in atomic detail with explicitly represented water molecules and ions. In each of the 7 simulations, the ligand repeatedly bound to and unbound from FKBP at least 100 times (Figure S1A and Table S1), allowing an estimate of the kinetic rate constants and equilibrium constants of binding with small statistical errors.

etermining the structural, thermodynamic, and kinetic properties of protein−ligand complexes is instrumental to modern drug discovery, and many experimental techniques are used to characterize such complexes. X-ray crystallography and nuclear magnetic resonance (NMR) spectroscopy, for example, can often be used to determine molecular structures, while surface plasmon resonance and calorimetric techniques can measure equilibrium binding constants. Some experimental methods can also determine the kinetic rates of protein−ligand binding. Limitations to the resolution and sensitivity of these methods, however, can be problematic when analyzing lowmolecular-weight ligands such as those used in fragment-based drug discovery.1−5 In contrast to larger ligands, fragments often have weak binding affinities that are below the detection limit of these experimental assays. Fragments may also bind in multiple poses and at multiple sites, introducing uncertainty in the characterization of atomic interactions between fragments and their receptors using X-ray crystallography and NMR spectroscopy.5 Molecular dynamics (MD) simulation has been used as an alternative, computational route toward characterizing the binding of fragments at an atomic level of spatial detail and with femtosecond temporal resolution.6−17 Simulations can be used to study an idealized molecular system free of certain complexities (such as impurities) that may affect the outcome of experimental studies of low-affinity fragments. A single MD simulation in which a fragment repeatedly binds to and unbinds from the protein offers the opportunity to predict the structure © 2017 American Chemical Society

Received: February 17, 2017 Published: June 5, 2017 3372

DOI: 10.1021/acs.jctc.7b00172 J. Chem. Theory Comput. 2017, 13, 3372−3377

Article

Journal of Chemical Theory and Computation

Figure 1. Molecular structures of drug fragments used in these simulations are shown along with their abbreviations: dimethyl sulfoxide (DMSO), 1hydroxy-2-propanone (PROP), methyl sulphinyl-methylsulfoxide (DSS), 4-hydroxy-2-butanone (BUT), tetrahydrothiophene 1-oxide (THI), and 4diethylamino-2-butanone (DAP). DAPP (not pictured) refers to the protonated form of DAP.

Table 1. Estimates of the Binding Free Energies by Direct Simulation of the Binding Equilibrium, ΔGb(direct), and by Free Energy Perturbation (FEP), ΔGb(FEP)a DMSO PROP DSS BUT THI DAP DAPP

ΔGb (direct)

Tdirect (μs)

1.5 (1) 1.5 (1) 2.9 (1) 2.01 (6) 3.4 (1) 3.9 (1) 1.39 (9)

30 20 20 39 30 20 20

ΔGb (FEP) 1.7 1.5 2.7 2.0 3.4 3.9 1.2

(3) (2) (2) (2) (2) (6) (2)

TFEP (μs)

σ2FEP/σ2direct

Tdirect/TFEP

0.85 0.68 0.85 0.88 0.41 0.988 0.248

9 2 5 10 1 16 7

35 29 24 44 73 20 81

Statistical errors are given in parentheses. The free energies are in kcal mol−1. The simulation time for FEP (TFEP) is the total simulation time of all the intermediate states in the replica exchange simulation for transferring the ligand from vacuum into the protein binding site (see Supporting Information); the simulation time (and the statistical error in the free energy estimate) for transferring the ligand from vacuum into aqueous solution is negligible in comparison. The ratio of the variance in the FEP estimate (σ2FEP) to that in the estimate by direct simulation (σ2direct) is smaller than the reverse ratio of the simulation times, showing that FEP is more efficient than the direct simulation in estimating binding free energies, even for this set of ligands with fast dissociation kinetics. a

μs). To obtain well-converged estimates of KD from unbiased simulations, simulation time of several times tR would be required, in agreement with the amount of simulation required to converge our KD estimates. Increasing ligand concentration may help shorten association times but will not shorten unbinding times, which are concentration-independent. Due to their poor specificity, drug fragments may bind at multiple sites and in several poses within each site. In our simulations, each drug fragment occupied several different sites on FKBP, including some sites which were occupied by all six fragments. The active-site region near residue TRP59, labeled S0a in Figure 3A for 4-hydroxy-2-butanone (BUT), contained the highest ligand occupancy. At S0a, BUT bound in multiple poses, including the crystallographically determined pose (Figure S1B). Similar behavior was observed for the other fragments (Figures S1B and S2). Information about multiple binding poses observed in simulations may be exploited in fragment-based drug discovery. Weak fragments, for example, can be modified by adding functional groups that form specific interactions and improve shape complementarity with the binding pocket. Multiple binding poses offer a larger set of possible positions to place such groups, expanding the number of possible chemical modifications that are energetically plausible. The best way to incorporate this information into the ligand design process, however, is still an open area of research.32,33 Fragments also tended to occupy two additional regions, S0b and S0c, adjacent to S0a (Figure 3B). Crystal structures show that FK506 (KD = 0.4 nM) and rapamycin (KD = 0.2 nM),34,35 larger drug molecules that bind with high affinity to FKBP,

Additional simulation details can be found in the Supporting Information. Because these simulations reversibly sample the bound and unbound states, calculation of dissociation constants31 (KD), the associated binding affinity (ΔGb), and kinetic on- and offrates (kon and koff, respectively) is relatively straightforward (detailed methods are described in the Supporting Information). Estimates of KD from the simulations ranged from 3 to 100 mM (Tables 1 and S1). Kinetically, the 6 molecules bound to the active site of FKBP with similar, diffusion-limited onrates (kon) (∼109 M−1 s−1). Off-rates (koff) ranged from 7 to 120 × 106 s−1, corresponding to residence times (tR = 1/koff) between 8 and 140 ns. Up to several microseconds of simulation were required to converge the calculated affinity for a fragment to a standard error within 0.5 kcal mol−1 (Figure 2A). Long time scale simulations that contain multiple binding events permit the computation of the equilibrium dissociation constant, KD, with high statistical precision because the variance in the estimate decreases in inverse proportion to the number of unbinding events observed in the simulations. The high-precision estimation of KD values allows molecules to be ranked with higher confidence. In general, we expect convergence times to be on the order of a microsecond for fragments of millimolar affinity. It is common for fragments to bind to a surface-exposed pocket with a diffusion-limited association rate. A diffusion-limited kon for fragments of millimolar affinity implies residence times of around 1 μs (because KD is between 1 and 10 mM and kon is 109 M−1 s−1, tR = 1/koff = (KD kon)−1 ranges from 100 ns to 1 3373


Article


have thus been developed to compute binding affinity without simulating the binding equilibrium.37−42 To validate these computational methodologies, computed binding free energies are often compared to experimental values. Lack of agreement between computation and experiment can, however, reflect a variety of discrepancies, including inaccuracies in the underlying physical models (“force fields”) used in simulation and differences between experimental and simulated conditions.43 Moreover, different experimental measurements performed on the same ligand−protein complex may give different values of the binding affinities.30,44 Our FKBP−fragment binding simulations, in contrast, provide a valuable benchmark to which computational methods for calculating binding affinities, as well as methods that estimate ligand-binding kinetics, can be directly compared using the same molecular system with the same force field parameters; any discrepancies observed in such comparisons can be attributed to problems with the methods themselves (as opposed to, for example, force field inaccuracies). One example of a computational technique that is commonly used to calculate binding affinities is free energy perturbation (FEP).45−49 FEP computes binding affinity through an artificial process in which, for example, the interactions between the ligand and the rest of the system are gradually turned off;50−53 the equilibrium between binding and unbinding is never explicitly simulated. We computed the binding free energies of the FKBP−ligand complexes using the FEP method implemented in Desmond23 (on the basis of binding only at the crystallographically determined site S0a, Figure 3A). To compare the free energies calculated by FEP to those estimated from direct binding simulations, we excluded proteinassociated states, in which the ligand is associated with the protein but not bound in the site of interest, from the definitions of the bound and unbound states in the binding simulations because FEP calculates the free energy of a ligand binding from bulk solvent to a particular binding site without consideration of intermediate protein-associated states. This adjustment made a difference of less than 0.2 kcal mol−1 rootmean-square error (RMSE) in our estimates of the equilibrium binding constants (which are shown in Table S1) because the free energy of binding at protein-associated states is relatively weak compared to the binding at S0a. Further discussion of protein-associated states and how they were defined can be found in the Supporting Information. The free energies calculated by FEP agreed with those estimated from the long time scale binding simulations within statistical error (Figure 2B). We note that, although the binding free energies calculated by FEP are in good agreement with those computed from the fragment-binding simulations, there are quantitative differences between the results from these two computational approaches, on the one hand, and experimental values (Figure 2B, inset) on the other, a discrepancy that is probably attributable at least in part to force field error. Although FEP methods are derived rigorously from statistical mechanics, and different free energy methods have been crossvalidated against one another,42 this is the first time, to our knowledge, that such a method has been validated numerically in an atomically detailed protein−ligand complex by direct comparison to unbiased simulations of binding equilibrium. The kinetics of drug binding can often impact drug efficacy and safety,54,55 an observation that has led to the development of several computational techniques for characterizing drug association and dissociation rates along with affinities.41,56−60

Figure 2. Binding affinities of millimolar fragments require microseconds of direct simulation to converge to values within 0.5 kcal mol−1 of their final values and agree quantitatively with binding affinities calculated using an FEP method. (A) Convergence of free energies of binding during a direct simulation. Solid lines represent the difference between the estimated binding affinity and its final value as a function of simulation time, and different colors correspond to different fragments. Black dashed lines correspond to the free energy interval between −0.5 and 0.5 kcal mol−1. (B) A comparison of binding free energies obtained from direct simulations as in panel A with those obtained for the same systems using FEP. The estimates by the two methods agree within the statistical error. The FEP results for both the deprotonated and protonated forms of DAP are plotted; DAP, with a pKa ≈ 10.8, is predominantly protonated at the experimental pH 7.6. The inset shows a comparison of binding affinity estimates from FEP calculations to experimental values.30 Here, only the protonated state of DAP is plotted; an additional experimental value of DMSO, determined with NMR, is also plotted.44 The RMSE between the FEP results (omitting the deprotonated DAP) and experimental values is 2.1 kcal mol−1; previous studies comparing FEP to experimental values have used different sets of FKBP ligands but have reported similar RMSEs of ∼2−3 kcal mol−1.53,62 Various nonFEP computational approaches have also been used to estimate the FKBP binding affinities for the ligands in this study, mostly in agreement with experiment in the 1−3 kcal mol−1 range.14,30,63−65

occupy both of these regions (Figure S3), implying that these fragments are able to map out regions of the protein that are capable of potent binding.12,36 In addition to this cluster of sites around S0, the fragments occupied various alternative binding sites in the simulations. In particular, site S1, lying 15 Å from S0a, is well-populated by all six fragments (Figures 3C and S2). In general, alternative binding sites mapped by fragments in simulation could potentially be druggable allosteric sites (although this is unlikely to be the case for a rigid and previously well-characterized protein such as FKBP). All-atom reversible binding simulations are computationally intensive and can be prohibitively expensive for ligands of submicromolar affinities. Many other computational methods 3374


Article


Figure 3. A drug fragment of millimolar affinity binds to several alternative binding sites in addition to the crystallographically determined binding site. (A) A contour map showing a projection of the x and y components of the center-of-mass position density of the fragment BUT as a free-energy surface. The crystallographically determined binding site is labeled S0a. Several well-populated alternative binding sites are also evident. Similar results were observed for the other fragments (Figure S2). (B) Examples of low-free-energy configurations of BUT binding in the S0 region: two alternative binding sites, S0b and S0c, located near S0a, are well-populated. These alternative sites are occupied by known high-affinity drug molecules such as FK506 and rapamycin (Figure S3). (C) A low-free-energy configuration of BUT bound in the most populated alternative binding site, S1, about 15 Å away from the S0 region. Protein residues interacting with the ligand in S1 are colored and labeled.

*E-mail: [email protected]; Phone: (212) 478-0260; Fax: (212) 845-1286.

The consistency we have found between results from FEP and long time scale fragment binding simulations lends additional support to an approach that combines these two tools to quantitatively characterize ligand binding and unbinding, even for larger, drug-like ligands. Association rates are often fast compared to off-rates, especially for high-affinity ligands, so a set of microsecond-scale ligand-binding simulations may be sufficient to estimate association rates and bound-ligand structures.59,60 Ligands that bind with diffusion-limited kinetics typically have on-rates of 107 to 109 M−1 s−1. Such a ligand would have an average on-time of approximately 100 ns to 10 μs given a ligand concentration of 10 mM (which is typical for an MD system with one ligand and one protein in a 60 Å cubic box); the residence time, tR, by contrast would be ∼1 s for a nanomolar ligand. Such an on-time could be estimated with a few to a few hundred microseconds of aggregate simulation time. The binding free energies of the poses observed in binding simulations can be estimated with free energy methods, and knowledge of the association rates and the binding affinities would give the off-rates.61 Such a combined computational approach to characterizing protein−ligand binding would be much more efficient than equilibrium binding-unbinding simulations for moderate- to high-affinity ligands.60

■

ORCID

David E. Shaw: 0000-0001-8265-5761 Author Contributions §

A.C.P., H.X., and T.P. contributed equally to the manuscript.

Notes

The authors declare the following competing financial interest(s): This study was conducted and funded internally by D. E. Shaw Research, of which David E. Shaw is the sole beneficial owner and Chief Scientist, and with which all authors are affiliated.

■

ACKNOWLEDGMENTS The authors thank Michael Eastwood, Kim Palmo, David Borhani, Cristian Predescu, Fabrizio Giordanetto, and Dina Sharon for helpful discussions and Rebecca Bish-Cornelissen and Berkman Frank for editorial assistance.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.7b00172. Table providing a list of simulations (Table S1), four additional figures (Figures S1−S4), and more details about the simulations (PDF) MD structure files corresponding to the equilibrated starting frames of the reversible binding simulations, including force field information (ZIP)

■

REFERENCES

(1) Congreve, M.; Chessari, G.; Tisi, D.; Woodhead, A. J. Recent developments in fragment-based drug discovery. J. Med. Chem. 2008, 51 (13), 3661−3180. (2) Murray, C. W.; Rees, D. C. The rise of fragment-based drug discovery. Nat. Chem. 2009, 1 (3), 187−192. (3) Hajduk, P. J.; Greer, J. A decade of fragment-based drug design: Strategic advances and lessons learned. Nat. Rev. Drug Discovery 2007, 6 (3), 211−219. (4) Erlanson, D. A.; Wells, J. A.; Braisted, A. C. Tethering: Fragmentbased drug discovery. Annu. Rev. Biophys. Biomol. Struct. 2004, 33, 199−223. (5) Davis, B. J.; Erlanson, D. A. Learning from our mistakes: The “unknown knowns” in fragment screening. Bioorg. Med. Chem. Lett. 2013, 23, 2844−2852. (6) Miranker, A.; Karplus, M. Functionality maps of binding sites: a multiple copy simultaneous search method. Proteins: Struct., Funct., Genet. 1991, 11 (1), 29−34. (7) Yang, C. Y.; Wang, S. Computational analysis of protein hotspots. ACS Med. Chem. Lett. 2010, 1 (3), 125−129. (8) Lexa, K. W.; Carlson, H. A. Full protein flexibility is essential for proper hot-spot mapping. J. Am. Chem. Soc. 2011, 133 (2), 200−202. (9) Ivetac, A.; McCammon, J. A. Mapping the druggable allosteric space of G-protein coupled receptors: a fragment-based molecular dynamics approach. Chem. Biol. Drug Des. 2010, 76 (3), 201−217.

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]; Phone: (212) 403-8664; Fax: (646) 873-2664. *E-mail: [email protected]; Phone: (212) 478-0163; Fax: (212) 845-1163. 3375


Article

Journal of Chemical Theory and Computation (10) Seco, J.; Luque, F. J.; Barril, X. Binding site detection and druggability index from first principles. J. Med. Chem. 2009, 52 (8), 2363−2371. (11) Guvench, O.; MacKerell, A. D., Jr. Computational fragmentbased binding site identification by ligand competitive saturation. PLoS Comput. Biol. 2009, 5 (7), e1000435. (12) Bakan, A.; Nevins, N.; Lakdawala, A. S.; Bahar, I. Druggability assessment of allosteric proteins by dynamics simulations in the presence of probe molecules. J. Chem. Theory Comput. 2012, 8 (7), 2435−2447. (13) Huang, D.; Caflisch, A. Small molecule binding to proteins: Affinity and binding/unbinding dynamics from atomistic simulations. ChemMedChem 2011, 6 (9), 1578−1580. (14) Huang, D.; Caflisch, A. The free energy landscape of small molecule unbinding. PLoS Comput. Biol. 2011, 7 (2), e1002002. (15) Wang, K.; Chodera, J. D.; Yang, Y.; Shirts, M. R. Identifying ligand binding sites and poses using GPU-accelerated Hamiltonian replica exchange molecular dynamics. J. Comput.-Aided Mol. Des. 2013, 27, 989−1007. (16) Ferruz, N.; Harvey, M. J.; Mestres, J.; De Fabritiis, G. Insights from fragment hit binding assays by molecular simulations. J. Chem. Inf. Model. 2015, 55 (10), 2200−2205. (17) Bisignano, P.; Doerr, S.; Harvey, M. J.; Favia, A. D.; Cavalli, A.; De Fabritiis, G. Kinteic characterization of fragment binding in AmpC β-lactamase by high-throughput molecular simulations. J. Chem. Inf. Model. 2014, 54 (2), 362−366. (18) Plimpton, S. Fast parallel algorithms for short-range molecular dynamics. J. Comput. Phys. 1995, 117 (1), 1−19. (19) Fitch, B. G.; Rayshubskiy, A.; Eleftheriou, M.; Ward, T. J. C.; Giampapa, M.; Zhestkov, Y.; Pitman, M. C.; Suits, F.; Grossfield, A.; Pitera, J.; Swope, W.; Zhou, R.; Germain, R. S.; Feller, S. Blue matter: Strong scaling of molecular dynamics on Blue Gene/L. Comput. Sci.− ICCS 2006, 3992, 846−854. (20) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kalé, L.; Schulten, K. Scalable molecular dynamics with NAMD. J. Comput. Chem. 2005, 26 (16), 1781−1802. (21) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4 (3), 435−447. (22) Case, D. A.; Cheatham, T. E., III; Darden, T.; Gohlke, H.; Luo, R.; Merz, K. M., Jr.; Onufriev, A.; Simmerling, C.; Wang, B.; Woods, R. J. The Amber biomolecular simulation programs. J. Comput. Chem. 2005, 26 (16), 1668−1688. (23) Bowers, K. J.; Chow, E.; Xu, H.; Dror, R. O.; Eastwood, M. P.; Gregersen, B. A.; Klepeis, J. L.; Kolossvary, I.; Moraes, M. A.; Sacerdoti, F. D.; Salmon, J. K.; Shan, Y.; Shaw, D. E. Scalable algorithms for molecular dynamics simulations on commodity clusters. In Proceedings of the ACM/IEEE Conference on Supercomputing (SC06); IEEE: New York, NY, 2006. (24) Harvey, M. J.; De Fabritiis, G. An implementation of the smooth particle mesh Ewald method on GPU hardware. J. Chem. Theory Comput. 2009, 5 (9), 2371−2377. (25) Harvey, M. J.; Giupponi, G.; De Fabritiis, G. ACEMD: Accelerating biomolecular dynamics in the microsecond time scale. J. Chem. Theory Comput. 2009, 5 (6), 1632−1639. (26) Phillips, J. C.; Stone, J. E.; Schulten, K. Adapting a messagedriven parallel application to GPU-accelerated clusters. In Proceedings of the ACM/IEEE Conference on Supercomputing (SC08); IEEE: New York, NY, 2008. (27) Salomon-Ferrer, R.; Götz, A. W.; Poole, D.; Le Grand, S.; Walker, R. C. Routine microsecond molecular dynamics simulations with Amber on GPUs. 2. Explicit solvent particle mesh Ewald. J. Chem. Theory Comput. 2013, 9 (9), 3878−3888. (28) Abraham, M. J.; Murtola, T.; Schulz, R.; Szilárd, P.; Smith, J. C.; Hess, B.; Lindahl, E. GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX 2015, 1−2, 19−25.

(29) Shaw, D. E.; Dror, R. O.; Salmon, J. K.; Grossman, J. P.; Mackenzie, K. M.; Bank, J. A.; Young, C.; Deneroff, M. M.; Batson, B.; Bowers, K. J.; Chow, E.; Eastwood, M. P.; Ierardi, D. J.; Klepeis, J. L.; Kuskin, J. S.; Larson, R. H.; Lindorff-Larsen, K.; Maragakis, P.; Moraes, M. A.; Piana, S.; Shan, Y.; Towles, B. Millisecond-scale molecular dynamics simulations on Anton. In Proceedings of the Conference on High Performance Computing, Networking, Storage and Analysis (SC09); ACM: New York, NY, 2009. (30) Burkhard, P.; Taylor, P.; Walkinshaw, M. D. X-ray structures of small ligand-FKBP complexes provide an estimate for hydrophobic interaction energies. J. Mol. Biol. 2000, 295 (4), 953−962. (31) De Jong, D. H.; Schäfer, L. V.; De Vries, A. H.; Marrink, S. J.; Berendsen, H. J.; Grubmüller, H. Determining equilibrium constants for dimerization reactions from molecular dynamics simulations. J. Comput. Chem. 2011, 32 (9), 1919−1928. (32) Silvestre, H. L.; Blundell, T. L.; Abell, C.; Ciulli, A. Integrated biophysical approach to fragment screening and validation for fragment-based lead discovery. Proc. Natl. Acad. Sci. U. S. A. 2013, 110 (32), 12984−12989. (33) Erlanson, D. A.; Fesik, S. W.; Hubbard, R. E.; Jahnke, W.; Jhoti, H. Twenty years on: the impact of fragments on drug discovery. Nat. Rev. Drug Discovery 2016, 15 (9), 605−619. (34) Bierer, B. E.; Mattila, P. S.; Standaert, R. F.; Herzenberg, L. A.; Burakoff, S. J.; Crabtree, G.; Schreiber, S. L. Two distinct signal transmission pathways in T lymphocytes are inhibited by complexes formed between an immunophilin and either FK506 or rapamycin. Proc. Natl. Acad. Sci. U. S. A. 1990, 87 (23), 9231−9235. (35) Schreiber, S. L. Chemistry and biology of the immunophilins and their immunosuppressive ligands. Science 1991, 251 (4991), 283− 287. (36) Raman, E. P.; Yu, W.; Guvench, O.; Mackerell, A. D. Reproducing crystal binding modes of ligand functional groups using Site-Identification by Ligand Competitive Saturation (SILCS) simulations. J. Chem. Inf. Model. 2011, 51 (4), 877−896. (37) Woo, H. J.; Roux, B. Calculation of absolute protein-ligand binding free energy from computer simulations. Proc. Natl. Acad. Sci. U. S. A. 2005, 102 (19), 6825−6830. (38) Mobley, D. L.; Dill, K. A. Binding of small-molecule ligands to proteins: ″What you see″ is not always ″what you get″. Structure 2009, 17 (4), 489−498. (39) Chodera, J. D.; Mobley, D. L.; Shirts, M. R.; Dixon, R. W.; Branson, K.; Pande, V. S. Alchemical free energy methods for drug discovery: progress and challenges. Curr. Opin. Struct. Biol. 2011, 21 (2), 150−160. (40) Steinbrecher, T. B.; Dahlgren, M.; Cappel, D.; Lin, T.; Wang, L.; Krilov, G.; Abel, R.; Friesner, R.; Sherman, W. Accurate binding free energy predictions in fragment optimization. J. Chem. Inf. Model. 2015, 55 (11), 2411−2420. (41) Buch, I.; Giorgino, T.; De Fabritiis, G. Complete reconstruction of an enzyme-inhibitor binding process by molecular dynamics simulations. Proc. Natl. Acad. Sci. U. S. A. 2011, 108 (25), 10184− 10189. (42) Gumbart, J. C.; Roux, B.; Chipot, C. Standard Binding Free Energies from Computer Simulations: What Is the Best Strategy? J. Chem. Theory Comput. 2013, 9, 794−802. (43) Krimmer, S. G.; Klebe, G. Thermodynamics of protein−ligand interactions as a reference for computational analysis: how to assess accuracy, reliability and relevance of experimental data. J. Comput.Aided Mol. Des. 2015, 29 (9), 867−883. (44) Dalvit, C.; Floersheim, P.; Zurini, M.; Widmer, A. Use of organic solvents and small molecules for locating binding sites on proteins in solutions. J. Biomol. NMR 1999, 14 (1), 23−32. (45) Wang, L.; Wu, Y.; Deng, Y.; Kim, B.; Pierce, L.; Krilov, G.; Lupyan, D.; Robinson, S.; Dahlgren, M. K.; Greenwood, J.; Romero, D. L.; Masse, C.; Knight, J. L.; Steinbrecher, T.; Beuming, T.; Damm, W.; Harder, E.; Sherman, W.; Brewer, M.; Wester, R.; Murcko, M.; Frye, L.; Farid, R.; Lin, T.; Mobley, D. L.; Jorgensen, W. L.; Berne, B. J.; Friesner, R. A.; Abel, R. Accurate and reliable prediction of relative ligand binding potency in prospective drug discovery by way of a 3376


Article

Journal of Chemical Theory and Computation modern free-energy calculation protocol and force field. J. Am. Chem. Soc. 2015, 137 (7), 2695−2703. (46) Boresch, S.; Tettinger, F.; Leitgeb, M.; Karplus, M. Absolute binding free energies: A quantitative approach for their calculation. J. Phys. Chem. B 2003, 107 (35), 9535−9551. (47) Gilson, M. K.; Given, J. A.; Bush, B. L.; McCammon, J. A. The statistical-thermodynamic basis for computation of binding affinities: A critical review. Biophys. J. 1997, 72 (3), 1047−1069. (48) Radmer, R. J.; Kollman, P. A. Free energy calculation methods: A theoretical and empirical comparison of numerical errors and a new method qualitative estimates of free energy changes. J. Comput. Chem. 1997, 18 (7), 902−919. (49) Pohorille, A.; Jarzynski, C.; Chipot, C. Good practices in freeenergy calculations. J. Phys. Chem. B 2010, 114, 10235−10253. (50) Darve, E. Thermodynamic integration using constrained and unconstrained dynamics. In Free energy calculations: Theory and applications in chemistry and biology; Chipot, C., Pohorille, A., Eds.; Springer series in chemical physics 86; Springer: New York, 2007; 119−170. (51) Oostenbrink, C.; van Gunsteren, W. F. Calculating zeros: Nonequilibrium free energy calculations. Chem. Phys. 2006, 323 (1), 102− 108. (52) Bennett, C. H. Efficient estimation of free energy differences from Monte Carlo data. J. Comput. Phys. 1976, 22 (2), 245−268. (53) Fujitani, H.; Tanida, Y.; Ito, M.; Jayachandran, G.; Snow, C. D.; Shirts, M. R.; Sorin, E. J.; Pande, V. S. Direct calculation of the binding free energies of FKBP ligands. J. Chem. Phys. 2005, 123 (8), 084108. (54) Swinney, D. C. The role of binding kinetics in therapeutically useful drug action. Curr. Opin. Drug. Discovery Devel. 2009, 12 (1), 31− 39. (55) Lu, H.; Tonge, P. J. Drug-target residence time: critical information for lead optimization. Curr. Opin. Chem. Biol. 2010, 14 (4), 467−474. (56) Decherchi, S.; Berteotti, A.; Bottegoni, G.; Rocchia, W.; Cavalli, A. The ligand binding mechanism to purine nucleoside phosphorylase elucidated via molecular dynamics and machine learning. Nat. Commun. 2015, 6, 6155. (57) Tiwary, P.; Limongelli, V.; Salvalaglio, M.; Parrinello, M. Kinetics of protein-ligand unbinding: predicting pathways, rates, and rate-limiting steps. Proc. Natl. Acad. Sci. U. S. A. 2015, 112 (5), E386− 391. (58) Mollica, L.; Decherchi, S.; Zia, S. R.; Gaspari, R.; Cavalli, A.; Rocchia, W. Kinetics of protein-ligand unbinding via smoothed potential molecular dynamics simulations. Sci. Rep. 2015, 5, 11539. (59) Shan, Y.; Kim, E. T.; Eastwood, M. P.; Dror, R. O.; Seeliger, M. A.; Shaw, D. E. How does a drug molecule find its target binding site? J. Am. Chem. Soc. 2011, 133 (24), 9181−9183. (60) Dror, R. O.; Pan, A. C.; Arlow, D. H.; Borhani, D. W.; Maragakis, P.; Shan, Y.; Xu, H.; Shaw, D. E. Pathway and mechanism of drug binding to G-protein-coupled receptors. Proc. Natl. Acad. Sci. U. S. A. 2011, 108 (32), 13118−13123. (61) We note that in the absence of external driving forces the unbinding process is the reverse of the binding process, following the same pathway and traversing the same barriers in the opposite order. (62) Wang, J.; Deng, Y.; Roux, B. Absolute binding free energy calculations using molecular dynamics simulations with restraining potentials. Biophys. J. 2006, 91 (8), 2798−2814. (63) Lee, M. S.; Olson, M. A. Calculation of Absolute protein-ligand binding affinity using path and endpoint approaches. Biophys. J. 2006, 90 (3), 864−877. (64) Swanson, J. M. J.; Henchman, R. H.; McCammon, J. A. Revisiting free energy calculations: A theoretical connection to MM/ PBSA and direct calculation of the association free energy. Biophys. J. 2004, 86 (1), 67−74. (65) Ytreberg, F. M. Absolute FKBP binding affinities obtained via nonequilibrium unbinding simulations. J. Chem. Phys. 2009, 130 (16), 164906.

3377


Quantitative Characterization of the Binding and Unbinding of

Recommend Documents