Estimation of Drug-Target Residence Times by Ï-Random

This is an open access article published under an ACS AuthorChoice License, which permits copying and redistribution of the article or any adaptations for non-commercial purposes.

Article pubs.acs.org/JCTC

Cite This: J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Estimation of Drug-Target Residence Times by τ‑Random Acceleration Molecular Dynamics Simulations Daria B. Kokh,*,† Marta Amaral,‡,§,¶ Joerg Bomke,∥ Ulrich Grad̈ ler,‡ Djordje Musil,‡ Hans-Peter Buchstaller,⊥ Matthias K. Dreyer,# Matthias Frech,‡ Maryse Lowinski,∇ Francois Vallee,∇ Marc Bianciotto,∇ Alexey Rak,∇ and Rebecca C. Wade*,†,○,◆ †

Molecular and Cellular Modeling Group, Heidelberg Institute for Theoretical Studies, Heidelberg 69118, Germany Molecular Interactions and Biophysics, Merck KGaA, Darmstadt 64293, Germany § Instituto de Biologia Experimental e Tecnológica, Oeiras 2780-157, Portugal ∥ Molecular Pharmacology, Merck KGaA, Darmstadt 64293, Germany ⊥ Medicinal Chemistry, Merck KGaA, Darmstadt 64293, Germany # R&D Integrated Drug Discovery, Sanofi-Aventis Deutschland GmbH, Frankfurt am Main 65926, Germany ∇ Integrated Drug Discovery, Sanofi R&D, Vitry-sur-Seine F-94403, France ○ Center for Molecular Biology (ZMBH), DKFZ-ZMBH Alliance, Heidelberg University, Heidelberg 69120, Germany ◆ Interdisciplinary Center for Scientific Computing (IWR), Heidelberg University, Heidelberg 69120, Germany

Downloaded via 185.251.15.173 on July 5, 2018 at 12:57:41 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

‡

S Supporting Information *

ABSTRACT: Drug-target residence time (τ), one of the main determinants of drug efficacy, remains highly challenging to predict computationally and, therefore, is usually not considered in the early stages of drug design. Here, we present an efficient computational method, τ-random acceleration molecular dynamics (τRAMD), for the ranking of drug candidates by their residence time and obtaining insights into ligand-target dissociation mechanisms. We assessed τRAMD on a data set of 70 diverse drug-like ligands of the N-terminal domain of HSP90α, a pharmaceutically important target with a highly flexible binding site, obtaining computed relative residence times with an accuracy of about 2.3τ for 78% of the compounds and less than 2.0τ within congeneric series. Analysis of dissociation trajectories reveals features that affect ligand unbinding rates, including transient polar interactions and steric hindrance. These results suggest that τRAMD will be widely applicable as a computationally efficient aid to improving drug residence times during lead optimization.

■

INTRODUCTION

current time accessible to conventional molecular dynamics (MD) simulations (μs), computation of unbinding rates becomes tremendously challenging in the framework of conventional MD alone.4 To date, unbinding rates have been estimated from conventional MD simulations only for molecules with fast dissociation rates (typically with kof f > 103 s−1).5−8 This fact has motivated the development of methods for the optimization of MD sampling algorithms by selection of reaction-relevant regions of protein−ligand configurational space (e.g., by Weighted Ensemble Sampling9−12 or the Adaptive Multilevel Splitting Method13) or by adding biasing potentials along selected collective variables, CV (e.g., using metadynamics14−20). The latter method has been applied to simulate the unbinding of benzamidine from trypsin,19

Over the past few years, much progress has been made in the development of methods for in silico prediction of drug binding affinity. However, the binding kinetics, which underlie the affinity, have only recently been appreciated as important factors for drug efficacy and safety.1,2 It has been demonstrated that the rates of drug binding to and unbinding from the target are the main features, along with pharmacokinetic properties, that determine a drug’s biological activity profile.3 Computation of the dissociation rate constant, koff, on which the residence time is dependent (τ = 1/koff), is generally more challenging than computation of the binding affinity. Unlike the binding affinity that can be estimated using a two-state end point approach, computation of the dissociation rate requires extensive sampling of the transition state that generally results from more than one pathway in the protein−ligand configurational space. Given that the residence time of pharmaceutically interesting compounds (minutes to hours) is far beyond the © XXXX American Chemical Society

Received: March 1, 2018 Published: May 16, 2018 A

DOI: 10.1021/acs.jctc.8b00230 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Figure 1. Crystal structures and binding kinetics of the 70 N-HSP90-inhibitor complexes studied. (a) The two main conformations of the N-HSP90 binding site, helix-type (left; PDB: 5J9X) and loop-type (right; PDB: 5J2X); the α-helix3 region is shown in red; the surface of the ATP binding pocket and the hydrophobic subpocket are colored in yellow and orange, respectively; red and black circles indicate substituents in the solventexposed region and in the transient subpocket, respectively. (b) Measured KD (left) and kon (right) values plotted against kof f. Crosses denote loopbinders; the other compounds are helix-binders. (c−m) Binding poses of compounds representing the 11 different ligand scaffolds considered. Crystal structures are shown in cartoon representation with the ligands and key side chains in stick representation and water molecules as red spheres (numbered as referred to in the text); the binding pocket is colored as in plots (a,b). Ligands are colored by atom-type and additionally shown in 2D in the insets (compounds (c) 13 (PDB: 5J9X), (d) 36 (PDB: 5L06), (e) 21 (PDB: 5LS1), (f) 60 (PDB: 5OD7), (g) 67 (PDB: 5LR1), (h) 17 (PDB: 5LRZ), (i) 19 (PDB: 2YKI), (j) 6 (PDB: 6F1N), (k) 69 (PDB: 5LO5), (l) 68 (PDB: 5EL5), and (m) 70 (PDB: 2YKJ)).

inhibitors from several kinases15,16 and HIV-1 protease,17 and antagonists from adenosine A2A receptor,18 showing agreement with experimental data to within an order of magnitude of koff in most cases. However, the selection of appropriate CVs and the long MD simulations required mean that these methods are difficult to apply in real-life drug design processes. Instead of computing absolute values of koff, it has been proposed that relative residence times can be estimated for a set of compounds directly from the simulation time required for

ligand dissociation in enhanced MD simulations.21 Specifically, scaling of the potential energy with a partially restrained target was employed to correctly rank several glucokinase activators and HSP90 inhibitors.21,22 Later, the relative egress times in metadynamics runs were shown to provide a reasonably good ranking of the order of magnitude of the relative kof f values of 10 inhibitors of cyclin-dependent kinase 8.23 However, both these methods require parametrization (either selection of a B


Article


SI Table 1, and analytical data are provided in SI Table 2. Compound 1 (AUY922) was purchased from MedChem Express (Sollentuna, Sweden), and compounds 2 (VER50589) and 3 (VER49009) were purchased from Cayman Chemicals (Hamburg, Germany). Crystallization and Structure Determination. Compounds 4, 6, 7, 10, 14, 35, 36, 38, 40−41, 62, 68, and 69: A hexahistidine tagged N-terminal fragment of HSP90α (residues 9−236) (NP_005339) was expressed and purified by the Instituto de Biologia Experimental e Tecnológica (Lisbon, Portugal), according to previously described protocols.28 The crystallization conditions were essentially the same as those described in refs 28 and 29. Data sets were collected at the SLS synchrotron and processed with the XDS software package.30 The structures were solved by molecular replacement with Phaser and refined either with CNX30 or BUSTER.31 Model building was performed with Coot,32 with inhibitors and water sites fitted into the initial |Fo|−|Fc| map. Compounds 17−22 and 67. A hexahistidine tagged Nterminal fragment of HSP90α (residues 18−223) (NP_005339) was expressed, purified, and crystallized as described in the Supporting Information. Data sets were collected in-house on a Rigaku HF-007 rotating anode generator and a MAR CCD detector and at the synchrotron. Diffraction data were processed with either XDS30 or MOSFLM.33 The structures were solved by molecular replacement using one set of coordinates of N-HSP90 (PDB ID: 1YER). The structures were refined using either CNX30, REFMAC5,34 or AUTOBUSTER.31 Ligands were placed manually, and the structural models were manually rebuilt using either TURBO-FRODO35 or COOT.32 Final validation checks were performed using MOLPROBITY.31 Data set statistics for all crystal structures are given in SI Table 2. The crystallographic coordinates have been deposited in the Protein Data Bank under the accession codes 5LQ9, 5LR1, 5LR7, 5LRZ, 5LS1, 5T21, 5NYI, 6F1N, 6ELO, 6ELN, 6ELP, 2YKI, 6EYA, 5LO6, 6EY8, 6EY9, 6EYB, 6EI5, 6EL5, and 5LO5. Surface Plasmon Resonance (SPR). Measurements were performed on a Biacore 4000 instrument (GE Healthcare) as previously described.28 Briefly, recombinant N-HSP90 with 17desmethoxy-17-N,N-dimethylaminoethylamino-geldanamycin (17-DMAG, Merck Millipore) was immobilized on a Biacore CM5 chip at 25 °C at a flow rate of 10 μL/min using amine coupling at pH 4.50 according to Biacore’s standard protocol. HBS-N (10 mM Hepes pH 7.40, 0.15 M NaCl) served as the running buffer during immobilization, and all SPR binding kinetics measurements were performed in 20 mM HEPES pH 7.50, 150 mM NaCl, 0.05% Tween 20, 1 mM DTT, 0.1 mM EDTA, 2% DMSO. Data sets were processed and analyzed using the Biacore 4000 Evaluation software, version 1.1. Solvent corrected and double-referenced association and dissociation phase data were fitted to a simple 1:1 interaction model with mass transport limitations.

scaling parameter or appropriate CVs), which is not always unambiguous and may notably affect simulation results. In the present study, we developed a new approach, called τrandom acceleration molecular dynamics (τRAMD), to the computation of the relative target residence times of drug-like compounds. It is based on the random acceleration molecular dynamics (RAMD)24 method, an enhanced sampling procedure that was originally developed for exploring ligand egress pathways from buried binding sites in proteins.25,26 In the RAMD method, during a standard molecular dynamics simulation of the bound complex, a small additional randomly oriented force is applied to the compound to facilitate its unbinding. The direction of this force is reassigned randomly when the ligand’s movement in a defined time interval falls below a specified threshold distance. The application of the random force thus allows egress of the ligand to be observed in short simulations of several nanoseconds duration. Unlike the methods mentioned above, τRAMD requires neither prior knowledge of the dissociation pathway nor extensive parameter fitting. The only parameter to be set by the user is the magnitude of the random force, which mainly affects the simulation time required and, as demonstrated in the present study, when chosen within a reasonable range, does not generally affect the computed relative residence times. Here, we evaluate the application of τRAMD to a diverse set of inhibitors of an important cancer target, the human N-terminal domain of heat shock protein 90α (N-HSP90). HSP90 is a chaperone protein that assists the folding of many types of protein, in particular those involved in cell growth. One way to disrupt the HSP90 chaperone activity is by blocking the ATPase function of N-HSP90 by small molecule inhibitors, which leads to degradation of the client proteins and diminished tumor growth.27 N-HSP90 is a challenging target for simulations due to its high flexibility. The ATP binding site undergoes considerable conformational changes during the ATPase cycle. Human N-HSP90 has been observed in several conformations (Figure 1a) that can be classified as loop- (loopin and loop-out) and helix-type. Both loop conformations have been observed in crystal structures of human N-HSP90 cocrystallized with different small inhibitors bound to the ATP-binding pocket, as well as in the unbound protein.28 In contrast, the helix conformation, with a complete α-helix3, has been observed only in holo-structures and only when the bound compound occupies a hydrophobic transient subpocket between α-helix3 and the beta-strands, in addition to the ATP binding site (see inset in Figure 1a), and thereby stabilizes αhelix3 (compounds bound to helical and loop-in conformations will be referred to hereafter as helix-binders and loop-binders, respectively). Here, we assess the performance of τRAMD on 70 N-HSP90 ATP-pocket binding compounds that have a range of different scaffolds (Figure 1, SI Table 1, SI Figure 1). For this purpose, we here measured the binding kinetic parameters for 29 compounds by surface plasmon resonance (SPR) and solved crystal structures of complexes with N-HSP90 for 20 compounds. We demonstrate that, in addition to providing estimates of target residence times, τRAMD gives insights into the ligand dissociation pathways and the mechanisms of unbinding.

■

COMPUTATIONAL METHODS Structure Preparation for Simulations. The crystal structures of 37 complexes were used for the simulations (see SI Table 1). For the remaining protein−ligand complexes simulated, the structures were modeled from the crystal structures with PDB identifiers: 2VCI, 5J9X, 5OD7, and 5OCI for resorcinol loop-binder, resorcinol helix-binder, quinazoline, and imidazole compounds, respectively. Ligands

■

EXPERIMENTAL METHODS Chemistry. Information on the synthesis of chemical compounds is provided in patents and publications listed in C


Article

Journal of Chemical Theory and Computation from the corresponding PDB files were used as templates for building new compounds by substitution with appropriate chemical groups. Ligand protonation at pH 7.5 was assigned using the Epik tool.36,37 Additionally, the protonation states of several compounds, whose residence times were systematically underestimated in simulations (7-imidazopyridine, 7-azaindole, and 2-aminopyridine compounds shown in Figure 1h, i, and m, respectively), were checked using Jaguar-pKa38 (see the Supporting Information), giving values that agreed with Epik results (see SI Table 4). The ligand parameters were generated using the Antechamber program.39 RESP40,41 atomic partial charges were fitted using electrostatic potentials from ab initio calculations performed at the HF level with a 6-31G*(1d) basis set using GAMESS.42 To account for the interactions of aryl halogens in polar surroundings, we modified the partial charges of two compounds, 5 and 17, that contain a Br atom, on the basis of the extra-charge model.43 As the NAMD44 program does not support virtual sites, we employed the partial negative charge of the Br atom obtained from the extra-charge model and distributed the positive charge of the virtual site over the aromatic carbon atoms. The protein (residues 16−223) was protonated assuming standard protonation states at pH 7.5 using Amber tools. The GAFF force field45 was employed for the ligands and the Amber14 force field46 for the protein. The crystallographic water molecules located within 3 Å of the ligand were included in the system. The tleap Amber tool47 was used for building a rectangular simulation box with a margin of 10 Å around the protein filled with TIP3P48 water molecules. Some bulk water molecules were replaced by chlorine or sodium ions to retain system neutrality and a salt concentration of about 150 mM. τRAMD Protocol. For each compound, the main steps of the τRAMD protocol are as follows: (i) preparation of a minimized and equilibrated structure; (ii) generation of starting bound-state replicas; (iii) simulation of dissociation trajectories for each starting replica; (iv) statistical analysis to provide a final residence time value and indicate the quality of the bound and transition state sampling. (i) The system was first energy minimized, heated, and equilibrated using the AMBER14 software.49 Minimization was done with 500 steps of steepest descent followed by 1000 steps of conjugate gradient minimization. Then the system was gradually heated to 300 K in 1 ns using harmonic restraints on non-hydrogen atoms with a force constant of 50 kcal mol−1 Å−2. Equilibration was carried out by decreasing the restraints in 5 steps (each of 1 ns with force constants of 50, 10, 5, 2, 0.5 kcal mol−1 Å−2) and then running 10 ns without restraints; a constant temperature of 300 K was maintained using a Langevin thermostat and a constant pressure of 1 atm using a Berendsen barostat. For all simulations, a cutoff of 10 Å for nonbonded Coulombic and Lennard-Jones interactions and periodic boundary conditions with a Particle Mesh Ewald treatment of long-range Coulombic interactions were used. A 2 fs time step was employed with bonds to hydrogen atoms constrained using the SHAKE algorithm.50 (ii) The atomic coordinates of the equilibrated system generated with AMBER were used as input for the generation of starting replicas in NAMD format. The

system was gradually heated over a period of 0.6 ns in steps of 10 K without restraints. Then equilibration was performed in the NPT ensemble using a Langevin thermostat and Nosé−Hoover Langevin pressure control to maintain the system at 1 atm and 300 K. For each compound, we carried out either 2 or 4 series of simulations. Two snapshots (starting replicas) were extracted from each trajectory with a stride of 5 ns, and their atomic coordinates and velocities were used as the starting points for RAMD simulations. (iii) The RAMD simulation procedure,24,51 implemented as a Tcl wrapper around the NAMD software,44 was modified to take the force magnitude rather than acceleration as an input parameter and to use new functions available from version 2.10 of NAMD onward. A random force with a magnitude of 14 kcal mol−1 Å−1 was used for simulations. The direction of the force was chosen randomly initially and after each check point (i.e., after 100 fs MD) either retained (if the ligand center of mass moved by at least 0.025 Å) or changed randomly otherwise. Simulations were stopped when the center of mass of the ligand had moved further than 30 Å from its original position. Coordinates were saved at 1 ps intervals. (iv) Initially, a set of 10 dissociation trajectories was generated for each starting replica. The residence time, τcomp, was defined as the simulation time required for ligand dissociation in 50% of the trajectories. A bootstrapping procedure (200 sets, each of the sets containing 80% of the trajectories chosen randomly) was employed to compute the residence time distribution for each set and to estimate the interim residence time, tr, as the mean of the distribution and its standard deviation SDr (tr distributions for all compounds are shown in SI Figure 2). The final value of τcomp was computed as an average over all simulated replicas as τcomp = ⟨tr⟩replicas with a standard deviation, SDreplicas. The uncertainty of τcomp is defined as the greater of two standard deviation values: SDcomp = max(⟨SDr⟩, SDreplicas). If SDr exceeded 50%, the number of simulation trajectories for this particular replica was increased by 5 repeatedly until this criterion was satisfied. If SDreplicas exceeded 50%, the number of replicas used for the simulations was increased up to a maximum of 8 replicas. Thus, a total of 40−200 RAMD trajectories was generated for each compound studied. Validation of the τRAMD Method. For the validation of the method, we considered the linear correlation between computed and measured residence time on a logarithmic scale for compounds i = 1, 2,···,N: Y = aX+b (where Y = log(τcomp) and X = log(τexpt)). Any compound with a Cook’s distance52 greater that its mean value multiplied by four over all sets was considered to be an outlier. The mean absolute error, MAE, of the linear regression model was computed as (see also SI Figure 3) i=N

MAE = 100% ∑ 1/N i=1

|τ i comp − 10Y 10Y

i fitted

|

i fitted

where Yifitted = a log(τiexpt) + b. The final value of the MAE was obtained using a bootstrapping procedure (200 groups; 80% of compounds in each group). Since the correlation between computed and measured residence times was analyzed on a logarithmic scale, the MAE value does not directly represent the accuracy of the method in D


Article


Figure 2. Analysis the computed relative residence times and simulation results for inhibitors of N-HSP90. (a) Computed residence times, τcomp, plotted against measured values, τexpt, for 70 compounds. Crosses denote loop-binders, the other compounds are helix-binders, and colors and shapes distinguish compound classes as shown in Figure 1c−m. Complexes with/without crystallographic structures are shown by filled/open symbols, respectively. The black line shows the linear regression of all points except for five outliers (compounds 11, 17, 30, 69, and 70, and the compounds shown in green, see text for details) outside the gray area that indicates the region within a residual standard deviation of the linear fitting computed with a confidence threshold of 0.95; linear regression analysis yields R2 = 0.86; error bars show measured and computational standard deviations (see SI Tables 5 and 6, respectively). (b) Distribution of the standard deviation of the computed residence time, SDr, for 30 compounds for different numbers of RAMD trajectories each starting from one replica (one snapshot): increasing the number of trajectories reduces SDr. (c) Illustrations of a typical distribution of the dissociation times, tr, obtained for one compound from simulations starting from different replicas (data from 6 simulations of compound 27, each starting with a different replica (different snapshots) are shown; the corresponding tr values vary from 1.3 to 3.4 ns and have a mean value of τcomp = 2.6 ± 0.7 ns).

Additionally, water #1 is coordinated by T184, and two or three other water molecules (denoted as #2, #3, and #4 in Figure 1c− m) form a network connecting D93 and the NH or OH moiety of the bound ligand. The compounds can be roughly classified into 11 groups according to the chemical moiety that forms Hbonds with D93, as illustrated in Figure 1c−m. These groups are (with the number of compounds given in brackets) as follows: resorcinol (27), hydroxyindazole (19), benzamide (3), aminoquinazoline (9), aminopyrrolopyrimidine (2), 7-imidazopyridine (1), 7-azaindole (2), aminothienopyridine (1), 6hydroxyindole (1), adenine (1), and 2-aminopyridine (1). In the helix-type structures, a hydrophobic ligand moiety partially or completely occupies the transient hydrophobic pocket adjacent to α-helix3, as indicated in Figure 1a (and in 2D protein−ligand contact maps in SI Figure 1). This part of the ligand, as well as solvent-exposed ligand substitutions (indicated by black and red circles, respectively, in Figure 1a), defines the differences within each group of compounds. The binding rate constants measured by SPR are given in Figure 1b and SI Table 5. The equilibrium dissociation constant, KD, and the dissociation rate constant, koff, vary over the data set by about 6 and 4 orders of magnitude, respectively. Importantly, the association rate constant, kon, also varies notably within most of the ligand scaffold types, indicating that both the bound state and the transition state affect koff values across the set of compounds. Although in the data set, both the binding affinities and the kinetic rates of loopand helix-binders each cover approximately the same range of values, the slowest binding compounds (kon < 105 M−1 s−1) are all helix-binders. This is in accord with the results of our previous study28 where, in a set of resorcinol-type compounds, all helix-binders were characterized by smaller KD values and

predicting τexpt. Therefore, we introduced a Mean Prediction Uncertainty (MPU) factor that provides the uncertainty of the relative values of the residence time obtained from computations assuming that Yifitted yelds a correct scaling of calculated to experimental data; MPU is defined as MPU = 1 N

⎛ 10 X icomp τiexp ⎞ ∑ max⎜⎜ i , X icomp ⎟⎟ ⎝ τexp 10 ⎠ i=1

i=N

where Xicomp = (log(τicomp) − b)/a. In the analysis of congeneric series, a linear regression model for the dependence of Y = aX + b, where Y = log(τcomp) − log(τrefcomp) and X = log(τexpt) − log(τrefexpt), was employed. The MAE and the MPU were computed as given above with Yifitted = a(log(τiexpt) − log(τrefexpt)) and Xicomp = (log(τicomp) − log(τrefcomp))/a.

■

RESULTS AND DISCUSSION The 70 Inhibitors Analyzed Bind to Two Conformations of N-HSP90 and Display a Wide Range of Binding Kinetic Parameters. The 70 compounds studied are listed in SI Table 1. For 20 compounds, the crystal structure of the protein−ligand complex was determined in the present study, see SI Table 3. These structures, along with the structures of an additional 17 complexes available in the PDB/RCSB, were used for the computations. The crystal structures show 9 loopbinders and 28 helix-binders (as classified by the structure of the α-helix3 segment, residues 101−122, see Figure 1a). A common feature in the ATP-binding pocket of all helix- and loop-type complexes is one direct and one water-mediated Hbond between the ligand and D93 (water #1 in Figure 1c−m). E


Article


Figure 3. Computations predict the correct trends in residence time upon introducing congeneric substitutions in six ligand groups. (a−f). Left: Each molecular scaffold is shown with the substitution fragments depicted in the table along with the compound number and the measured value of kof f; the reference compound number for each series is shown in bold; Right: Change in log(τcomp) upon making substitutions relative to a reference compound against corresponding change in measured log(τexpt) values (left plot) and change in log(τexpt) vs log(1/KD) (right plot). Linear regressions are shown by black lines. The color scheme corresponds to that for the ligand classes in Figure 2a.

lower kon and koff values than the loop-binders, which was associated with the protein structural changes required for formation of the helix-type complexes. Computed Relative Residence Times, τcomp, Correlate with Measured Residence Times. τcomp values computed for all compounds are plotted against τexpt (τexpt = 1/koff) obtained from SPR experiments on the logarithmic scale in Figure 2a. Most of the compounds demonstrate a good correlation between computed and measured residence times. However, τcomp is systematically underestimated for the 10 compounds of the amino-quinazoline and amino-pyrrolopyrimidine classes (58−67), which are shown by green rectangles and triangles, respectively, in Figure 2a, as well as for compound 70, which does not retain its crystallographic binding pose in MD equilibration simulations (SI Figure 4). Linear regression analysis omitting these compounds gives R2 = 0.66. Further, compounds, 11, 17, 30, and 69, can be clearly distinguished as outliers (using Cook’s distance, see Methods). Possible reasons for the simulation inaccuracy for these compounds are discussed below. Removal of these outliers gives an R2 value of the linear regression of the remaining 55 compounds (i.e., 78% of the original set studied) of 0.86 with a mean absolute error, MAE, of 36% (see details of statistical analysis in

Methods). Assuming that the linear regression analysis provides the correct rescaling of the experimental to computed residence time, the mean prediction uncertainty, MPU, of the τ is about 2.3-fold (individual values are given in SI Table 6). Sampling of Both Bound and Transition States Is Essential for Accurate Prediction of the Residence Times. The uncertainty of the computed residence time, τcomp, depends on the sampling accuracy of both the transition and the bound state. The transition state is explored by simulation of multiple dissociation trajectories, and its sampling efficiency can be quantitatively characterized by the standard deviation, SDr, of the residence time tr for a particular starting replica. The SDr value generally decreases as the number of simulated trajectories is increased, reaching SDr < 50% for 25 trajectories for most compounds. This is demonstrated in Figure 2b for selected cases (i.e., compounds and starting replicas) where computed SDr was above 50% in a set of 10 trajectories. For the same cases as shown in Figure 2b, we additionally carried out a statistical analysis of the τRAMD dissociation time distributions. Since the observed dissociation events are independent, the generated dissociation time distributions must approach Poisson statistics if the sampling is sufficiently large (see, for example, a similar analysis in ref F


Article


standard MD equilibration simulations that causes the underestimation of τcomp for these compounds. Indeed, although the computed residence times of all amino-quinazoline compounds are underestimated relative to the others, they nonetheless correlate reasonably well among themselves with the measured values (R2 = 0.97, see the next section). The Effects of Fragment Substitution on the Residence Time Are Correctly Reproduced by Simulations. The analysis given above suggests that the ability of the τcomp simulation procedure to predict relative residence time may depend on the accuracy of the description of the bound state. Thus, the prediction accuracy of the method can be expected to be higher if compounds from the same series are considered. To assess this idea, we selected 6 subgroups of at least five compounds with a very similar chemical scaffold (Figure 3a−f). In each group, one compound can be transformed into another by a few small substitutions, which leads to changes in the residence times of 1−4 orders of magnitude, except for group f where this variation is less than 2fold. In Figure 3, the variation in the log(τcomp) due to fragment substitution is plotted versus the corresponding difference in measured residence time, log(τexpt), along with a plot of log(τexpt) versus log(1/KD). In all cases, the direction of the change in τcomp is computed correctly, and the correlation between computed and experimental values is high (R2 = 0.87− 0.99 for all groups except f). It is important to note that the change of τexpt upon fragment substitution arises from modulation of both the bound and the transition states. For example, in the hydroxy-indazole series shown in Figure 3a, substitutions of the solvent-exposed phenyl fragment (compound 44) affect mainly the bound state leading to a decrease/increase in KD with a simultaneous prolongation/ reduction of τexpt for all but five compounds (35, 42, 45, 46, 55), for which the residence time does not correlate with the binding affinity. Despite such diverse effects of fragment substitutions on the bound and transition states, the ranking of the computed residence times agrees well with experiment (R2 = 0.94) with a prediction uncertainty of only 1.5-fold of the residence time. The compounds in the congeneric series shown in Figure 3b have different buried fragments in the transient subpocket. Substitution of the buried pyrrolidine ring of compound 51 by a piperidine ring (compound 57) and/or adding a methoxy group decreases the binding affinity but prolongs τexpt, while replacement of the pyrrolidine ring by an aliphatic chain (compound 37) solely increases compound affinity. Nonetheless, also for this series, computations show a good correlation with τexpt, which demonstrates a generally correct balance of transition and bound state energies in RAMD simulations. Even the group of amino-quinazoline compounds, whose binding is underestimated in MD simulations, exhibits a good correlation between computed and measured changes in the residence time upon substitution (Figure 3c) with a prediction uncertainty of 2.0τ. Computed Relative Residence Times Are Generally Robust to Minor Alterations in Simulated Structures and Simulation Parameters. Egress of fast dissociating compounds (koff > 10−2 s−1) from the binding pocket generally requires less than 1 ns of simulation under the random force used in the study (force magnitude of 14 kcal mol−1 Å−1). Fast dissociating compounds of the quinazoline series need even less than 0.1 ns for dissociation. Such rapid egress could lead to a reduced sensitivity of the method to different dissociation mechanisms. On the other hand, slower dissociation is more

53). Indeed, a two-sample Kolmogorov−Smirnov test reveals a high similarity between the cumulative distribution function of dissociation times observed in the τRAMD simulations and the Poisson distribution with a characteristic time of tr as indicated by a p-value above 0.89 for the majority of compounds when a set of 25 trajectories was used for analysis (see SI Figure 5). Furthermore, variation of the starting coordinates and velocities of the protein and ligand in the bound state leads to variations in the computed residence time as illustrated in Figure 2c. Increasing the number of starting replicas is an obvious way to improve the computational accuracy. For most of the compounds, the SDreplica of the residence time variations is below 50% for 4−8 replicas obtained from 10 to 20 ns MD trajectories (SI Figure 2; values for each compound are given in SI Table 6). However, some protein conformational changes may occur on a longer time scale, making proper sampling of the bound state infeasible. Specifically, the most critical conformational variation found was the rotation of the F138 side chain in helix-type complexes. In several complexes (such as, with compounds 24 and 26), F138 flipped during the MD equilibration procedure from the buried position observed in most crystal structures to a solvent-exposed orientation (SI Figure 6). The exposed F138 opens a buried subpocket for a ligand substituent and sterically hinders ligand dissociation, which generally leads to prolongation of the dissociation time in simulations. As a result, τcomp has quite a large uncertainty for these compounds (SDreplicas = 88% and 70% for compounds 24 and 26, respectively). This problem does not have a simple solution. Indeed, on the one hand, the observed residence time should be defined mainly by the fastest dissociation route (i.e., the most common orientation of F138) if the two starting states are equally populated. On the other hand, a reliable estimate of the population of the two possible binding site conformations can be obtained only from extensive simulations. Most Outliers Result from Deficiencies in the Description of the Bound State. The group of compounds with clearly underestimated residence times consists of all amino-quinazoline, amino-pyridopyrimidine, 2-aminopyridine, and 7-imidazopyridine compounds, as well as two resorcinol compounds, 11 and 30 (Figure 2a). The first two groups of compounds exhibit a similar binding motif with D93 and water #3 hydrogen-bonding to an amino group, as well as with water #4 mediating hydrogen-bonding between the heterocyclic nitrogen and N51(Figures 1f and g). This hydrogen-bonding network in the ATP binding pocket was found to be less stable in equilibration MD simulations than for the structures of compounds with similar affinity but different scaffolds (SI Figure 7). Similarly, the structure of the complex of compound 70 (2-aminopyridine) was found to be unstable, as mentioned above. Further, the resorcinol derivative outliers (compounds 11 and 30) have a furan ring positioned in the hydrophobic transient pocket. The furan substitution increases binding affinity dramatically according to the SPR measurements, but this is not observed in simulations (SI Figure 8 and Results in the Supporting Information). Underestimation of the protein− ligand binding strength in all the above cases leads to a reduction of the computed residence time. Finally, too fast dissociation was also observed for compounds having a bromobenzene group unless the redistribution of the partial charge on Br to mimic Br polarization was modeled43 (SI Figure 9). These results suggest that it is not the τRAMD procedure but rather the deficiency in the description of the bound state in the G


Article


Figure 4. Analysis of simulated RAMD trajectories for three benzamide compounds. (a) Chemical structures and measured kof f values of the three benzamide compounds. (b) Positions of compounds 20 and 21 in the binding pocket showing occupation of the subpocket under α-helix3 by compound 20, which has the longest residence time of these three compounds; the compounds are shown in ball-and stick-representation with carbon atoms colored in cyan, and the protein surface is colored by the electrostatic potential. (c) Variation in log(kon) values versus variation log(1/ KD) with respect to compound 21. (d) Computed log(τcomp) against measured log(τexpt) upon fragment substitution in compound 21. (e) Histogram showing RMSD of the protein binding site during dissociation trajectories for benzamide compounds. (f) Dissociation routes of compound 20 observed in RAMD simulations shown by blue lines indicating the ligand center-of-mass trajectories. Three representative dissociation trajectories are indicated by the positions of the oxygen atom in the benzamide fragment (magenta spheres). Insets: ligand poses in two snapshots from the corresponding dissociation trajectories.

sensitive to variations in the free energy landscape and thus requires more extensive sampling. To explore the sensitivity of the method to the random force magnitude, we repeated simulations for 30 compounds with kof f > 10−2 s−1 using a smaller force of 13 kcal mol−1 Å−1. This prolonged the dissociation time by several folds for most compounds (SI Figure 10) except for the amino-quinazoline compounds with koff > 0.1 s−1 (additional simulations showed that the dissociation time for these compounds starts to increase only for force magnitudes less than 11 kcal mol−1 Å−1). Importantly, the decrease of the force magnitude did not affect the outliers or notably change the accuracy of the simulations, which suggests that the accuracy of the computed residence time is mainly dependent on the quality of the MD simulations rather than on the magnitude of the random force. Importantly, we note that the quality of results of the simulations started with and without crystal structures (indicated in Figure 2a by filled and open symbols, respectively) does not differ significantly. Proper sampling of the bound state is critical even if a crystal structure of the complex is available. As discussed above, large variations of the computed residence time usually arise from long time scale protein motions and are independent of whether a crystal structure is available or not. Finally, it is also noteworthy that computed residence times for loop- and helix-binders fit similarly well to the linear regression line. Analysis of Dissociation Trajectories Reveals Different Factors Affecting the Residence Time. RAMD simulations sample multiple dissociation pathways that can involve various ligand and protein conformations, different protein−ligand contacts, and sometimes follow different egress routes. Analysis of these trajectories may give insights into the determinants of binding kinetics. Here, by way of illustration, we discuss three benzamide helix-binders, compounds 20, 21, and 22, shown in Figure 4a,b. Compound 20 differs from compound 21 by the introduction of a bulky fragment that occupies the hydrophobic

subpocket, whereas compound 22 differs from compound 21 only in the polarity of its solvent-exposed fragment (cyclohexanol is substituted by cyclohexanone). Although the residence times of these compounds are not significantly different, their binding rates, kon, and affinity vary by over 2 orders of magnitude (Figure 4c). In particular, compound 20, has a smaller affinity and a smaller kon than compound 21, i.e. weaker binding with a higher transition barrier. These two effects lead to a relatively small change in τexpt, which is reflected in the computed τcomp (Figure 4d). The simulations reveal several egress routes for compound 20 showing exit from both the side of the ATP binding site and the side of the transient subpocket (Figure 4f), which requires essential conformational changes. In contrast, compounds 21 and 22 can leave the binding pocket without notable effect on the protein structure (the distribution of RMSD of the binding site residues relative to the bound state in dissociation trajectories is illustrated in Figure 4e). This suggests that the sampling of energetically unfavorable conformations, which enable a bulky compound to leave a tight pocket, leads to a higher energy barrier to unbinding. In agreement with this, we observed an inverse correlation of kon with the number of heavy atoms in the hydroxy-indazole compounds of group b with different buried fragments (R2 = 0.68), whereas there is a very weak correlation (R2 = 0.16) for the group of compounds with various exposed fragments (SI Figure 11e,f). Finally, in general for all helix-binders, koff correlates well with the number of heavy atoms (SI Figure 11b). Comparison of compounds 21 and 22 suggests another mechanism of modulation of the transition state. The introduction of cyclohexanone in compound 22 instead of cyclohexanol in compound 21 stabilizes both the bound state and the transition state according to SPR data (both kon and 1/ KD are larger by about an order of magnitude for compound 22 than compound 21, Figure 4c). A possible reason is the H


Article


according to the longest and shortest dissociation time observed in the simulations. Additionally, the two-sample Kolmogorov−Smirnov (KS) test can be used to assess the statistics of dissociation times distribution in each set of simulated trajectories. The lower threshold of the force applied depends on the computer resources available for simulations. One has to keep in mind, however, that a smaller force will not only increase dissociation time but also lead to a more rugged potential landscape, which requires more extensive sampling. We demonstrated that sufficient sampling of both bound and transition states is critical for the accurate prediction of residence time. Consistently, the applicability of the method is limited to cases where the bound state can be sufficiently explored in the short standard MD equilibration simulations employed in the protocol. To some extent, the τRAMD is similar to the scaled MD approach reported by Mollica et al.21 which is also aimed at the calculation of the relative residence times of compounds. Instead of an additional force, the potential energy of the system is scaled to reduce the potential barriers to unbinding. In the scaled MD protocol, to prevent protein unfolding, additional restraints are applied to the backbone non-hydrogen atoms of all residues that do not belong to the binding site. This means that, in addition to a potential energy scaling parameter, a selection of restrained atoms suitable for the complete set of ligands binding to a given target should be made, and this selection can affect the potential energy scaling parameter required. Such restraints are not required in the τRAMD method and could impact the sampling of egress pathways that are dependent on backbone motions. Furthermore, the scaling of the potential energy in scaled MD hinders the identification of structural features determining the transition barrier to unbinding. In addition, an important advantage of the τRAMD method is sampling over the initial bound state. Since the time scale of the dissociation process in enhanced sampling methods does not allow for sampling the bound state properly, using multiple starting replicas is critical for computing residence times. Apart from the prediction of relative residence time, the τRAMD method can also be used to extract the mechanism of dissociation and to characterize transition states. Our study suggests that there is a large variety of factors that may alter the barrier to ligand unbinding, including transient polar interactions and steric hindrance upon ligand egress. In the case of N-HSP90 inhibitors, the tendency for the helix-binders to have greater residence times appears to result from contacts in the subpocket adjacent to α-helix3 and the necessity for conformational changes of the ligand and the protein during ligand dissociation. Automated analysis of a large number of trajectories to reveal further features responsible for the modulation of the transition barrier to unbinding is the subject of further investigations. In conclusion, the τRAMD method described provides a computationally inexpensive but rather accurate tool for the ranking of low molecular weight compounds by their target residence time and can be used for exploring the chemical space in ligand optimization procedures.

interaction of compound 22 with several charged residues on the pocket boundary, such as K58 and D54 (Figure 4f). Desolvation of a ligand and a protein binding site upon binding is often considered as an important factor affecting ligand binding rate.19,2954 We did not observe any significant correlation between kinetic rates and the ligand solvation energy computed using 3D-RISM (SI Table 6) in any of the groups of compounds studied (R2 < 0.1; SI Figure 11a−f; the maximum correlation, R2 = 0.25, is observed for helix-type compounds if outliers are excluded; SI Figure 11b). Although this result does not rule out the influence of water on the binding rate, it suggests that the ligand and pocket solvation are not leading effects in determining the transition barrier on the dissociation pathway for the present set of compounds. Altogether, these results suggest that there is no single factor or single leading factor that determines the transition barrier to unbinding, instead different factors are important for different compounds. Further studies for other targets and compound types may reveal different properties determining ligand unbinding kinetics.

■

CONCLUSIONS In the present study, we propose a new method, τRAMD, for the prediction of relative residence times of protein-inhibitor complexes. The approach permits a ligand dissociation time spanning from seconds to hours or more to be reduced to the 1−10 ns time scale of RAMD simulations for computing τ. In the τRAMD method, a set of standard MD simulations and RAMD simulations (with the latter employing an additional randomly oriented force on the ligand; important is that the force magnitude applied in τRAMD simulations is the same for all compounds simulated) is performed for each protein−ligand complex for which τ is to be computed. A total of 160 ns of MD simulation per compound was on average required in this work, see the Supporting Information. We show how the τRAMD method can be used to rank a diverse set of compounds by residence time and to explore the effects on residence time of changes in substituents in congeneric and diverse sets of compounds. Furthermore, the τRAMD method provides insights into the mechanistic features of ligand egress that affect τ. We assessed τRAMD for a set of 70 diverse drug-like inhibitors of N-HSP90 with kinetic parameters measured by the SPR technique. This application demonstrated a good correlation (R2 = 0.86) between computed residence times, τcomp, and measured τexpt values for 78% of the compounds. The average prediction uncertainty of the residence time was 2.3τ. For all but one of the remaining 22% of the compounds, the residence time was underestimated, most likely due to deficiencies of the force field. Consistently, we find that τRAMD generally provides the most accurate prediction of residence times when applied to a series of similar compounds with the same binding scaffolds and one or two substitutions. In τRAMD, the magnitude of the random force applied is the sole adjustable parameter. We have demonstrated that the method is quite robust to altering this value within a reasonable range. In general, we suggest using a magnitude of the random force that ensures ligand dissociation times on the nanosecond time scale, which enables effective sampling of side-chain motion during dissociation. Specifically, for a set of compounds with koff expected to vary within the range of about 10−1−10−4 s−1, a random force magnitude of 14 kcal mol−1 Å−1 can be applied, with subsequent small adjustments, if necessary,

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.8b00230. I


Article


■

(2) Copeland, R. A. The Drug−target Residence Time Model: A 10Year Retrospective. Nat. Rev. Drug Discovery 2016, 15, 87−95. (3) Yin, N.; Pei, J.; Lai, L. A Comprehensive Analysis of the Influence of Drug Binding Kinetics on Drug Action at Molecular and Systems Levels. Mol. BioSyst. 2013, 9 (6), 1381−1389. (4) Bruce, N. J.; Ganotra, G. K.; Kokh, D. B.; Sadiq, S. K.; Wade, R. C. New Approaches for Computing Ligand−receptor Binding Kinetics. Curr. Opin. Struct. Biol. 2018, 49, 1−10. (5) Huang, D.; Caflisch, A. The Free Energy Landscape of Small Molecule Unbinding. PLoS Comput. Biol. 2011, 7 (2), e1002002. (6) 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. (7) Pan, A. C.; Xu, H.; Palpant, T.; Shaw, D. E. Quantitative Characterization of the Binding and Unbinding of Millimolar Drug Fragments with Molecular Dynamics Simulations. J. Chem. Theory Comput. 2017, 13 (7), 3372−3377. (8) Plattner, N.; Noé, F. Protein Conformational Plasticity and Complex Ligand-Binding Kinetics Explored by Atomistic Simulations and Markov Models. Nat. Commun. 2015, 6 (1), 7653. (9) Dickson, A.; Lotz, S. D. Multiple Ligand Unbinding Pathways and Ligand-Induced Destabilization Revealed by WExplore. Biophys. J. 2017, 112 (4), 620−629. (10) Dickson, A.; Lotz, S. D. Ligand Release Pathways Obtained with WExplore: Residence Times and Mechanisms. J. Phys. Chem. B 2016, 120 (24), 5377−5385. (11) Lotz, S. D.; Dickson, A. Unbiased Molecular Dynamics of 11 min Timescale Drug Unbinding Reveals Transition State Stabilizing Interactions. J. Am. Chem. Soc. 2018, 140 (2), 618−628. (12) Zwier, M. C.; Pratt, A. J.; Adelman, J. L.; Kaus, J. W.; Zuckerman, D. M.; Chong, L. T. Efficient Atomistic Simulation of Pathways and Calculation of Rate Constants for a Protein-Peptide Binding Process: Application to the MDM2 Protein and an Intrinsically Disordered p53 Peptide. J. Phys. Chem. Lett. 2016, 7 (17), 3440−3445. (13) Teo, I.; Mayne, C. G.; Schulten, K.; Lelièvre, T. Adaptive Multilevel Splitting Method for Molecular Dynamics Calculation of Benzamidine-Trypsin Dissociation Time. J. Chem. Theory Comput. 2016, 12 (6), 2983. (14) Laio, A.; Parrinello, M. Escaping free-energy minima. Proc. Natl. Acad. Sci. U. S. A. 2002, 99, 12562−12566. (15) Casasnovas, R.; Limongelli, V.; Tiwary, P.; Carloni, P.; Parrinello, M. Unbinding Kinetics of a p38 MAP Kinase Type II Inhibitor from Metadynamics Simulations. J. Am. Chem. Soc. 2017, 139 (13), 4780−4788. (16) Tiwary, P.; Mondal, J.; Berne, B. J. How and When Does an Anticancer Drug Leave Its Binding Site? Sci. Adv. 2017, 3 (5), e1700014. (17) Sun, H.; Li, Y.; Shen, M.; Li, D.; Kang, Y.; Hou, T. Characterizing Drug−Target Residence Time with Metadynamics: How To Achieve Dissociation Rate Efficiently without Losing Accuracy against Time-Consuming Approaches. J. Chem. Inf. Model. 2017, 57 (8), 1895. (18) Bortolato, A.; Deflorian, F.; Weiss, D. R.; Mason, J. S. Decoding the Role of Water Dynamics in Ligand-Protein Unbinding: CRF1R as a Test Case. J. Chem. Inf. Model. 2015, 55 (9), 1857−1866. (19) 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−E391. (20) De Vivo, M.; Masetti, M.; Bottegoni, G.; Cavalli, A. The Role of Molecular Dynamics and Related Methods in Drug Discovery The Role of Molecular Dynamics and Related Methods in Drug Discovery. J. Med. Chem. 2016, 59 (9), 4035−4061. (21) 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 (1), 11539.

Methods: (1) expression, purification and crystallization, (2) ab initio computations of pKa values, (3) estimate of computing time required per compound; Results: nature of some outliers; Figures: (1) schematic illustration of protein−ligand contacts, (2) t r distributions, (3) estimation of uncertainty of τ values from simulations, (4) structural rearrangement of compound 70, (5) analysis of the dissociation time distribution for selected compounds, (6) rearrangement of the side chain of F138, (7) instability of amino-quinazoline compounds in the binding site, (8) effect of furan substitution, (9) crystal structure with compound 5, (10) effect of the magnitude of the RAMD force applied, (11) computed molecular descriptors; Tables: (1) compounds and crystal structures, (2) analytical data for selected compounds, (3) experimental structural data, (4) analysis of pKa values for selected compounds, (5) measured kinetic and thermodynamic data, (6) computed simulation and molecular descriptor data (PDF) Table containing SMILES strings for all compounds used in the study (XLSX)

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. *E-mail: [email protected]. ORCID

Daria B. Kokh: 0000-0002-4687-6572 Marta Amaral: 0000-0002-3503-9963 Rebecca C. Wade: 0000-0001-5951-8670 Present Address ¶

M.A.: Sanofi-Aventis Deutschland GmbH, R&D Biologics Research, Frankfurt am Main, Germany.

Funding

This work was supported by the EU/EFPIA Innovative Medicines Initiative (IMI) Joint Undertaking, K4DD (grant no. 115366) and by the Klaus Tschira Foundation. We gratefully acknowledge PRACE (Project: Pra12_3089) and HLRS Stuttgart (Project: Dynathor) for access to computing facilities. Notes

The authors declare no competing financial interest. The R script for analysis of RAMD trajectories and the Tcl script for performing RAMD simulations for this article may be accessed at no charge at https://www.h-its.org/downloads/ ramd.

■

ACKNOWLEDGMENTS We thank Gaurav Ganotra for help in using the extra-charge model during the ligand preparation procedure, Stefan Richter for assistance in porting software code, and S. Kashif Sadiq for useful discussions.

■ ■

ABBREVIATIONS RAMD, random acceleration molecular dynamics; HSP90, heat shock protein 90 REFERENCES

(1) Vauquelin, G. Effects of Target Binding Kinetics on in Vivo Drug Efficacy: Koff, Kon and Rebinding. Br. J. Pharmacol. 2016, 173 (15), 2319−2334. J


Article

Journal of Chemical Theory and Computation (22) Mollica, L.; Theret, I.; Antoine, M.; Perron-Sierra, F.; Charton, Y.; Fourquez, J. M.; Wierzbicki, M.; Boutin, J. A.; Ferry, G.; Decherchi, S.; et al. Molecular Dynamics Simulations and Kinetic Measurements to Estimate and Predict Protein-Ligand Residence Times. J. Med. Chem. 2016, 59 (15), 7167−7176. (23) Callegari, D.; Lodola, A.; Pala, D.; Rivara, S.; Mor, M.; Rizzi, A.; Capelli, A. M. Correction to Metadynamics Simulations Distinguish Short- and Long-Residence-Time Inhibitors of Cyclin-Dependent Kinase 8. J. Chem. Inf. Model. 2017, 57 (2), 159−169. (24) Lüdemann, S. K.; Lounnas, V.; Wade, R. C. How Do Substrates Enter and Products Exit the Buried Active Site of Cytochrome P450cam? 2. Steered Molecular Dynamics and Adiabatic Mapping of Substrate Pathways. J. Mol. Biol. 2000, 303 (5), 813−830. (25) Niu, Y.; Li, S.; Pan, D.; Liu, H.; Yao, X. Computational Study on the Unbinding Pathways of B-RAF Inhibitors and Its Implication for the Difference of Residence Time: Insight from Random Acceleration and Steered Molecular Dynamics Simulations. Phys. Chem. Chem. Phys. 2016, 18 (7), 5622−5629. (26) Peräkylä, M. Ligand Unbinding Pathways from the Vitamin D Receptor Studied by Molecular Dynamics Simulations. Eur. Biophys. J. 2009, 38 (2), 185−198. (27) Tillotson, B.; Slocum, K.; Coco, J.; Whitebread, N.; Thomas, B.; West, K. A.; MacDougall, J.; Ge, J.; Ali, J. A.; Palombella, V. J.; et al. Hsp90 (Heat Shock Protein 90) Inhibitor Occupancy Is a Direct Determinant of Client Protein Degradation and Tumor Growth Arrest in Vivo. J. Biol. Chem. 2010, 285 (51), 39835−39843. (28) Amaral, M.; Kokh, D. B.; Bomke, J.; Wegener, A.; Buchstaller, H. P.; Eggenweiler, H. M.; Matias, P.; Sirrenberg, C.; Wade, R. C.; Frech, M. Protein Conformational Flexibility Modulates Kinetics and Thermodynamics of Drug Binding. Nat. Commun. 2017, 8 (1), 2276. (29) Schuetz, D. A.; Richter, L.; Amaral, M.; Grandits, M.; Graedler, U.; Musil, D.; Buchstaller, H.; Eggenweiler, H.; Frech, M.; Ecker, G. F. Ligand Desolvation Steers on-Rate and Impacts Drug Residence Time of Heat Shock Protein 90 (Hsp90) Inhibitors. J. Med. Chem. 2018, 61 (10), 4397−4411. (30) Kabsch, W. XDS. Acta Crystallogr., Sect. D: Biol. Crystallogr. 2010, 66 (2), 125−132. (31) Bricogne, G.; Blanc, E.; Brandl, M.; Flensburg, C.; Keller, P.; Paciorek, W.; Roversi, P.; Sharff, A.; Smart, O.; Vonrhein, C.; et al. Buster, Version 2.11.6; Global Phasing Ltd.: Cambridge, United Kingdom; 2016. (32) Emsley, P.; Lohkamp, B.; Scott, W. G.; Cowtan, K. Features and Development of Coot. Acta Crystallogr., Sect. D: Biol. Crystallogr. 2010, 66 (4), 486−501. (33) Battye, T. G. G.; Kontogiannis, L.; Johnson, O.; Powell, H. R.; Leslie, A. G. W. iMOSFLM: A New Graphical Interface for Diffraction-Image Processing with MOSFLM. Acta Crystallogr., Sect. D: Biol. Crystallogr. 2011, 67 (4), 271−281. (34) Murshudov, G. N.; Vagin, A. A.; Dodson, E. J. Refinement of Macromolecular Structures by the Maximum-Likelihood Method. Acta Crystallogr., Sect. D: Biol. Crystallogr. 1997, 53 (3), 240−255. (35) Roussel, A.; Cambillau, C. The Turbo-Frodo Graphics Package. Silicon Graphics Geometry Partners Directory; Silicon Graphics Corp.: Mountain View, CA, 1991; p 81. (36) Shelley, J. C.; Cholleti, A.; Frye, L. L.; Greenwood, J. R.; Timlin, M. R.; Uchimaya, M. Epik: A Software Program for pKa Prediction and Protonation State Generation for Drug-like Molecules. J. Comput.Aided Mol. Des. 2007, 21 (12), 681−691. (37) Greenwood, J. R.; Calkins, D.; Sullivan, A. P.; Shelley, J. C. Towards the Comprehensive, Rapid, and Accurate Prediction of the Favorable Tautomeric States of Drug-like Molecules in Aqueous Solution. J. Comput.-Aided Mol. Des. 2010, 24 (6−7), 591−604. (38) Bochevarov, A. D.; Harder, E.; Hughes, T. F.; Greenwood, J. R.; Braden, D. A.; Philipp, D. M.; Rinaldo, D.; Halls, M. D.; Zhang, J.; Friesner, R. A. Jaguar: A High-Performance Quantum Chemistry Software Program with Strengths in Life and Materials Sciences. Int. J. Quantum Chem. 2013, 113 (18), 2110−2142.

(39) Wang, J.; Wang, W.; Kollman, P. A.; Case, D. A. Automatic Atom Type and Bond Type Perception in Molecular Mechanical Calculations. J. Mol. Graphics Modell. 2006, 25 (2), 247−260. (40) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Kollmann, P. A. Application of RESP Charges To Calculate Conformational Energies, Hydrogen Bond Energies, and Free Energies of Solvation. J. Am. Chem. Soc. 1993, 115 (7), 9620−9631. (41) Bayly, C. C. I.; Cieplak, P.; Cornell, W. D.; Kollman, P. A. A Well-Behaved Electrostatic Potential Based Method Using Charge Restraints for Deriving Atomic Charges: The RESP Model. J. Phys. Chem. 1993, 97 (40), 10269−10280. (42) Gordon, M. S.; Schmidt, M. W. Chapter 41 − Advances in Electronic Structure Theory: GAMESS a Decade Later. In Theory and Applications of Computational Chemistry. The First Forty Years; Dykstra, C. E., Frenking, G., Kim, K. S., Scuseria, G. E., Eds.; Elsevier: Amsterdam, 2005; pp 1167−1189. (43) Ibrahim, M. A. A. Molecular Mechanical Study of Halogen Bonding in Drug Discovery. J. Comput. Chem. 2011, 32 (12), 2564− 2574. (44) 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. (45) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and Testing of a General Amber Force Field. J. Comput. Chem. 2004, 25 (9), 1157−1174. (46) Maier, J. A.; Martinez, C.; Kasavajhala, K.; Wickstrom, L.; Hauser, K.; Simmerling, C. ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters from ff99SB. J. Chem. Theory Comput. 2015, 11 (8), 3696−3713. (47) Case, D. A.; Betz, R. M.; Botello-Smith, W.; Cerutti, D. S.; Cheatham, T. E., III; Darden, T. A.; Duke, R. E.; Giese, T. J.; Gohlke, H.; Goetz, A. W.; Homeyer, N.; Izadi, S.; Janowski, P.; Kaus, J.; Kovalenko, A.; Lee, T. S.; LeGrand, S.; Li, P.; Lin, C.; Luchko, T.; Luo, R.; B. Madej, L. X.; York, D. M.; Kollman, P. A. AMBER 2016; University of California: San Francisco, 2016. (48) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79 (2), 926−935. (49) Case, D. A.; Babin, V.; Berryman, J. T.; Betz, R. M.; Cai, Q.; Cerutti, D. S.; Cheatham, T. E.; Darden, T. A.; Duke, R. E.; Gohlke, H.; Goetz, A. W.; Gusarov, S.; Homeyer, N.; Janowski, P.; Kaus, J.; Kolossv, P. A. Amber 14; 2014. (50) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. Numerical Integration of the Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of N-Alkanes. J. Comput. Phys. 1977, 23 (3), 327−341. (51) Lüdemann, S. K.; Carugo, O.; Wade, R. C. Substrate Access to Cytochrome P450cam: A Comparison of a Thermal Motion Pathway Analysis with Molecular Dynamics Simulation Data. J. Mol. Model. 1997, 3, 369−374. (52) Williams, D. A. Generalized Linear Model Diagnostics Using the Deviance and Single Case Deletions. Appl. Stat. 1987, 36 (2), 181− 191. (53) Salvalaglio, M.; Tiwary, P.; Parrinello, M. Assessing the Reliability of the Dynamics Reconstructed from Metadynamics. J. Chem. Theory Comput. 2014, 10 (4), 1420−1425. (54) Bortolato, A.; Deflorian, F.; Weiss, D. R.; Mason, J. S. Decoding the Role of Water Dynamics in Ligand-Protein Unbinding: CRF 1 R as a Test Case. J. Chem. Inf. Model. 2015, 55 (9), 1857−1866.

K


Estimation of Drug-Target Residence Times by Ï-Random

Recommend Documents

Estimation of Drug-Target Residence Times by Ï-Random