Predicting Residence Time and Drug Unbinding Pathway

Article Cite This: J. Chem. Inf. Model. 2019, 59, 535−549

pubs.acs.org/jcim

Predicting Residence Time and Drug Unbinding Pathway through Scaled Molecular Dynamics Doris A. Schuetz,†,#,○ Mattia Bernetti,‡,# Martina Bertazzo,‡,§ Djordje Musil,∥ Hans-Michael Eggenweiler,⊥ Maurizio Recanatini,‡ Matteo Masetti,*,‡ Gerhard F. Ecker,† and Andrea Cavalli*,‡,§ †

Department of Pharmaceutical Chemistry, University of Vienna, UZA 2, Althanstrasse 14, 1090 Vienna, Austria Department of Pharmacy and Biotechnology, Alma Mater StudiorumUniversità di Bologna, via Belmeloro 6, I-40126 Bologna, Italy § Computational Sciences, Istituto Italiano di Tecnologia, via Morego 30, 16163 Genova, Italy ∥ Discovery Technologies, Merck KGaA, Frankfurter Straße 250, 64293 Darmstadt, Germany ⊥ Medicinal Chemistry, Merck KGaA, 64293 Darmstadt, Germany

Downloaded via BUFFALO STATE on July 28, 2019 at 17:03:19 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

‡

S Supporting Information *

ABSTRACT: Computational approaches currently assist medicinal chemistry through the entire drug discovery pipeline. However, while several computational tools and strategies are available to predict binding affinity, predicting the drug− target binding kinetics is still a matter of ongoing research. Here, we challenge scaled molecular dynamics simulations to assess the off-rates for a series of structurally diverse inhibitors of the heat shock protein 90 (Hsp90) covering 3 orders of magnitude in their experimental residence times. The derived computational predictions are in overall good agreement with experimental data. Aside from the estimation of exit times, unbinding pathways were assessed through dimensionality reduction techniques. The data analysis framework proposed in this work could lead to better understanding of the mechanistic aspects related to the observed kinetic behavior.

■

INTRODUCTION

effectively integrating kinetic information into this pipeline are currently missing. Molecular dynamics (MD) simulation is the method of choice when it comes to investigating dynamic processes at a fully atomistic level.10,15,16 However, in the context of protein− ligand binding or unbinding, reaching an optimal trade-off between efficiency and accuracy represents a challenging task.17−19 While protein−ligand binding events can nowadays be achieved through MD-based methods either relying on specialized hardware 20−23 or smarter sampling strategies,15,24−27 obtaining quantitative estimates of unbinding rates is currently out of reach for the majority of pharmaceutically relevant complexes. To circumvent this difficulty, multiscale simulations and biased-MD approaches aimed at accelerating the sampling of rare events (e.g., unbinding) have been introduced. In this context, multiscale modeling typically exploits a dual resolution scheme, in which the computationally cheaper Brownian dynamics is employed at high protein− ligand separation, whereas a fully atomistic description is switched-on in proximity of the binding site.28 Notably, the SEEKR method introduced by Amaro and co-workers takes

During the past decade, interest has been raised about binding and unbinding kinetics in the context of drug design.1−7 The traditional optimization strategies, mostly based on the thermodynamic signature of compounds,8,9 had to make room for kinetic considerations, which were shown to be potentially as important as affinity. This is especially true for residence time, which is a measure of how long a drug stays in contact with its target. It has been shown that the residence time can translate better to in vivo efficacy than affinity,6 and it should be considered as a key parameter at the early stages of drug discovery. Furthermore, the kinetic behavior of a drug toward its target can influence selectivity,10 drug safety,11 and PK−PD translation,12 making the evaluation of on- and offrates of paramount importance for drug discovery programs. Taking into account binding and unbinding kinetics, however, requires expensive setup and time-consuming testing of several newly synthesized compounds in complex kinetic assays.13 All of these issues call for fast and reliable methods to predict binding kinetics through in silico approaches. Indeed, while several computational strategies exist that are routinely employed to assist the drug discovery and design process from a thermodynamic standpoint,9,14 established tools for © 2018 American Chemical Society

Received: September 11, 2018 Published: November 30, 2018 535

DOI: 10.1021/acs.jcim.8b00614 J. Chem. Inf. Model. 2019, 59, 535−549

Article

Journal of Chemical Information and Modeling

Figure 1. (A) N-terminal domain of Hsp90 in complex with compound 2. The ligand is embedded inside the ATP binding pocket, whose hydrophobic and hydrophilic regions are herein highlighted in orange and blue, respectively. (B) Focus on the bound states. Representative structures of the three scaffolds (compounds 2, 4, and 7, from top to bottom) are shown (only polar hydrogens are displayed). The hydrogen bond interaction with the key residue Asp93 is maintained by all of the three compounds. Moreover, all ligands accommodate the same regions of the binding pocket.

approaches, we mention the adaptive multilevel splitting,41 as well as the weighted ensemble algorithm.42 In particular, the latter has been successfully employed to a series of pharmaceutically relevant systems of increasing complexity.43−46 Among the previously introduced tempering methods, scaled MD (sMD) is a promising approach that achieves a speed-up in the observation of rare events by smoothing the potential energy surface (PES) through the introduction of a scaling factor λ ranging from 0 to 1.13,47,48 sMD has been shown to successfully rank congeneric series of drug-like compounds according to their residence time, holding great potential for drug prioritization based on kinetic arguments.13,49 Very recently, the methodology has also been challenged to predict the unbinding kinetics for a series of hDAAO inhibitors in a prospective fashion, with no a priori knowledge of the experimental residence times.50 Despite the efficiency of sMD and its ability to predict relative off rates of protein−ligand complexes accurately, extracting relevant information on the dissociation mechanism from repeated sMD runs is far from trivial. The reason is the highdimensional configuration space in which unbinding events take place, as it comprises at least as many degrees of freedom as the Cartesian coordinates of the ligand’s heavy atoms (if protein conformational changes are neglected). Moreover, owing to the smoothed PES, in sMD the unbinding events occur as if the simulations were performed at high effective temperature. Thus, the ligand could follow nonphysical unbinding routes that might be far away from the minimum free energy surface, making the dissociation pathway rather difficult to interpret. Here, we used sMD to describe the kinetic behavior of several drug-like molecules in a reliable and efficient way. In particular, in the present study we challenge the methodology in two respects: (i) ranking of structurally diverse drug-like

advantage of both a multiscale modeling and an highly parallelizable milestoning approach to further improve the sampling efficiency within the atomistic regime.29,30 Concerning biased-MD approaches, most of them rely on the introduction of external biasing forces acting on selected degrees of freedom (or collective variables). This is the case for umbrella sampling and metadynamics, among many others.31,32 Although biased methods are primarily suited to achieve a thermodynamic characterization of the investigated process, the metadynamics formalism has recently also been extended to recover kinetic information from the biased dynamics.33 The methodology, named infrequent metadynamics, has been successfully applied to the well-known trypsin-benzamidine system,34 to characterize the unbinding kinetics and mechanistic steps of the anticancer drug dasatinib from c-Src kinase,35 and to predict the off-rates for a series of congeneric inhibitors from the p38 MAP kinase.36 Alternatively, one can accelerate all the degrees of freedom at once through high temperature dynamics or smoothed potential methods (tempering methods). In these cases, the canonical distribution of states is either preserved using suitable simulation setup (parallel tempering)37 or recovered in a following step through reweighting procedures.38 In between these classes of biased approaches lies the RAMD method, in which forces are randomly assigned to the ligand encouraging its egress from the binding site without the specification of collective variables.39 Notably, Kokh et al. very recently reported on the use of such method (τRAMD) to rank diverse compounds according to their dissociation time.40 Rather than biasing the dynamics through external forces or taking advantage of multiscale representations, another class of methods is currently emerging to predict the drug-target residence time. These are based on adaptive schemes whereby the dynamics is guided toward the most relevant regions of phase space to efficiently sample the rare event. Among these 536


Article


integrating the equations of motion. The neighbor list cutoff within the Verlet cutoff scheme was set to 1.1 nm. Periodic boundary conditions (pbc) were applied and smooth particle mesh Ewald method71 was used to calculate Coulombic interactions with a grid size of 0.16 nm. The cut-off for longrange Coulomb interactions was set to 1.1 nm and short-range van der Waals (vdW) interactions were similarly treated with a cutoff of 1.1 nm. The Lincs algorithm (Linear Restraint Solver)72 was used to restrain all protein and ligand bonds while the SETTLE algorithm73 was employed for water molecules. All boxes contained around 35 000 atoms (depending on ligand size). Scaled MD48 was carried out performing 25 repeated runs per ligand employing a scaling factor (λ coefficient) of 0.45. Harmonic positional restraints were applied to the backbone of the protein, as the protein fold had to be preserved throughout the simulation. Also, 50 kJ mol−1 nm−2 restraints were applied to the heavy atoms of the backbone, while the residues composing the binding site, in 5 Å of the ligand were not restrained at all. Detecting one atom of any residues within a 5 Å radius from any atoms of the ligand was the condition leading to the exclusion of the entire amino acid from the restraints. In this way, while the overall fold of the protein was preserved, the binding site of the protein was allowed to freely rearrange upon ligand unbinding. We note that, because of the peculiar shape of the solvent exposed binding site, the same set of residues was excluded from the restraints for all the compounds investigated. The position output frequency was set to 10 ps, with a precision of 0.001 nm. A Python script running within VMD was used to stop the simulations after the ligand had left the protein.74 The ligand was considered unbound when the distance between the protein center of mass (COM) and the ligand COM reached 30 Å. A total number of 175 simulations were performed for a cumulative simulations time of almost 27 μs (see Table S2). From the computed exit times recorded, the average unbinding times (corresponding to the computational residence time, τcalc) were determined for each ligand, along with standard deviations, standard errors, and bootstrapped standard errors. A 1000 fold bootstrapping was applied to assess the robustness of the estimation procedure as reported by Mollica et al.13 (see Figure S1 for distribution of single runs). Bootstrapping was done using R Studio.75 The correlation of normalized computational residence times with the normalized experimental residence times, τexp, was also performed as described by Mollica et al.13 Postprocessing of Unbinding Trajectories and Data Analysis. The analysis strategy employed in this work consisted in three main steps: (i) a cleanup stage in which sMD trajectories were scanned to obtain a productive sequence of states and reduce the noise introduced by the high effective temperature, (ii) construction of a low-D embedded space from the cleaned trajectories, and (iii) similarity analysis of unbinding pathways projected on the previously built low-D space. In particular, the low-D space was constructed to identify relevant features of unbinding events taking place in the high-dimensional configurational space of the simulated systems. The working hypothesis was that, despite the fact that sMD trajectories might be far away from minimum free energy paths, recurring events can reveal useful reaction coordinates to describe the protein−ligand binding process.55,76 Indeed, while the relative population of transition pathways (unbinding events) can change with temperature,77

ligands based on computed residence times and (ii) catching peculiar features of the unbinding process through subsequent analysis of the produced trajectories. Concerning the first aspect, we chose the heat shock protein 90 (Hsp90) as a test case (see Figure 1A). In particular, we predict the ranking of unbinding rates for a series of structurally diverse inhibitors that have been addressed in various studies as potential anticancer agents.51−54 The validity of the obtained prioritization is confirmed through comparison against experimental kinetic data. As for the second aspect, this study illustrates how relevant information about a drug’s dissociation process can potentially be retrieved taking advantage of dimensionality reduction techniques and cluster analysis. We show that the complex configurational space explored by ligands during dissociation can be successfully mapped into a low dimensional space through a nonlinear multidimensional scaling technique.55 Additionally, the similarity of unbinding pathways is assessed directly on the low-dimensional space using the Fréchet distance and well-established clustering methods.56 Taken together, our results demonstrate how protein−ligand dissociation studies through sMD can provide valuable quantitative and qualitative insights. These can be ultimately exploited to drive kinetic-based drug design.

■

METHODS Preparation of the Systems and Simulation Setups. The crystal structures of human Hsp90α NTD in complex with compounds 2 (PDB: 5OD7), 3 (PDB: 5NYH), 4 (PDB: 5LNY), 6 (PDB: 5ODX), and 7 (PDB: 6HHR) were used as starting coordinates for the simulations.57 Compound 1 was modeled from compound 2 (PDB: 5OD7). Compound 5 was modified by manually exchanging one carbon atom for an oxygen atom of compound 4 (PDB: 5LNY). All the protein preparations and setups were done using the software BIKI 1.3.5.58 Compounds 1 and 3−7 were assumed to be neutral, while compound 2 was in a protonation state +1. The protein−inhibitor complexes were placed inside a box of 7.58 nm length. Between 9 731 and 13 918 TIP3P59 water molecules were added for solvation. For systems 1 and 3−7, water molecules were replaced by 7 sodium ions, and for system 2, 6 sodium ions were added, to preserve electroneutrality of the simulation system. The Amber99SBildn60 and GAFF61,62 force fields were used for the protein and ligands, respectively. RESP charges63 for ligands were obtained from the ab initio optimized geometry computed at the HF/6-31G* level of theory with NwChem.64 The steepest descent method was used to minimize initial coordinates, in a 5000 step minimization run. Three equilibration steps in the NVT ensemble were carried out for 100 ps, starting at a temperature of 100 K, increasing during the three steps up to 300 K (coupling period: 0.1 ps). The first equilibration step was performed employing a 1 fs integration step, while a 2 fs integration step was used for the remaining stages. Positional restraints (an isotropic force constant of 1000 kJ mol−1 nm−2) were applied to the heavy atoms of the protein backbone. The final equilibration step of 1 ns was performed in the NPT ensemble at 300 K and 1 atm. The pressure-coupling was performed using the Parrinello−Rahman barostat65,66 (target pressure 1.013 bar, period 2 ps, compressibility 4.57 × 10−5 bar−1). An integration step 1 fs was set and the resulting simulation box length was used for the following production runs of different duration. All simulations were performed with Gromacs 4.6.1,67−69 using the velocity Verlet algorithm70 for 537


Article


Figure 2. (A) Schematic diagram of the cleanup stage of data analysis. (B) Pictorial representation of frame removal from a hypothetical trajectory which folds back to itself. (C) Main steps involved in the analysis of unbinding trajectories. After cleaning up the trajectories, the pairwise distance calculated in the configurational space of the ligand (previously aligned on the protein’s frame) is computed using the whole ensemble of preserved frames (matrix D in the scheme) and embedded in the low-D space through Isomap. Then, each unbinding pathway corresponding to the previously obtained cleaned-up trajectories are projected in the low-D space. Finally, for each pair of unbinding pathways, the Fréchet distance is computed and stored in a symmetric matrix (F in the scheme) that can be further exploited for subsequent cluster analysis.

we assume that, in the limit of an adequately large ensemble of sMD trajectories, the minimum free energy pathway can still be sampled with a relatively high probability. The low-D space was obtained through the Isomap method78 as implemented in the python library Pysomap.79 Isomap is a well-established nonlinear dimensionality reduction technique. Unlike metric and nonmetric multidimensional scaling (MDS) methods,80 which attempts to preserve the Euclidean distance of the high dimensionality (high-D) space, Isomap approximates the manifold using the geodesic distance. In this way, such an approach is superior to MDS as it allows reproduction of not only pairwise distances but also the inherent topology of the high-D space. This is particularly relevant if one wishes to devise a reaction coordinate from MD trajectories.81 Indeed, configurations close in distance in the high-D space are not necessarily close along a reactive pathway.55 In order to estimate the geodesic distance, the shortest path between pair of points in a k-nearest neighbor graph is computed. In particular, Pysomap implements the Floyd’s algorithm to create the network. Herein, the distance between the ligands’ heavy atom coordinates (after optimal alignment on the protein alpha carbons) computed over all pairs of trajectory frames was used as a metric. The PLUMED software was employed for this aim.82 It is well-known that the ability of Isomap to preserve the geodesic distance of the manifold depends on the sampling density of data and on the definition of neighboring points.83 In this work, the high-dimensional data were embedded into a 3D space using 10 k-nearest neighbor points (k = 10). This value was obtained by a trial and error procedure and comparing the features of the unbinding pathways projected in the low-D space (see below). Concerning the issue of data sampling, we note that even though a high density of points in

the high-D space is required for accuracy, the method becomes extremely impractical in case of big data like those generated by typical MD runs. To face this problem, landmark variants of Isomap have been developed, in which only a fraction of the entire data set is employed to build the embedding, rather than exploiting all of the available points. This fraction can be chosen either randomly or based on geometric assumptions.83 Here, we use a landmark-like approach, in which the relevant points in the high-D space are selected along productive unbinding trajectories. Thus, only those states describing a gradual advancement from the initial bound pose to fully solvated ligands were retained. This choice has the advantage of focusing the construction of the embedding on the portions of the manifold that carry the most relevant mechanistic information about the investigated events. Moreover, the unbinding trajectories obtained through sMD were extremely noisy. As such, this strategy also represents an effective procedure to filter out all the detours, dead-ends, and loops in the configurational space that are typical for nonequilibrium dynamics.84 In the cleanup stage of the data analysis, which precedes the dimensionality reduction, the trajectories from the bound to the unbound state of the ligand (running frame jm; m = 1, 2, ..., M) are scanned in order to save all the in frames whose distance in terms of root mean-squared deviation (RMSD), after alignment on the protein’s alpha carbons, is greater or equal to a given threshold (δ) compared to the latest frame saved (iN). To achieve this, the current structure (frame jm) is compared to all the N previously saved frames (in; n = 1, 2, ..., N). In case the previous condition is satisfied, the jmth configuration is stored as additional frame in the cleaned-up trajectory (see Figure 2A). On the other hand, if the RMSD between the running frame and one of the previously stored 538


Article

Journal of Chemical Information and Modeling Table 1. Structures of the Seven Compounds Examineda

a

Three scaffolds have been considered containing the quinazoline, indazole, and resorcinol cores. Experimental kinetic parameters (kon and koff) and binding affinities (KD) are reported for each chemical.57 SPR measurements for compound 7 were performed on a Biacore 4000 instrument from GE Healthcare (see the Supporting Information, SI). All experimental data are shown as the mean. The experimental residence time, 1 calculated as τ exp = k , is shown in seconds. off

structures is lower than the threshold (i.e., RMSD(in, jm) < δ), the inth frame is replaced by the current configuration, and all the subsequent points of the cleaned-up trajectory (n + 1, n + 2, ..., N) are discarded (see Figure 2A). In this way, only productive transitions are preserved, and all loops and/or dead ends along the unbinding trajectory are filtered out (see Figure 2B). Notably, this aspect becomes particularly relevant in light of the possibility to exploit the extracted pathway as a putative reaction coordinate, where continuous advancement along the considered mechanism becomes a crucial requisite. In this work, the RMSD was computed between all heavy atoms of ligands, and we used an RMSD threshold of 2.0 Å for compounds 1 and 7 and 2.5 Å for compounds 2−5. This choice reflected the different size of the investigated compounds (see Table 1), as we observed that a higher RMSD threshold was required to better describe the unbinding pathways of molecules carrying a greater amount of degrees of freedom. Finally, for every ligand, the low-D space was built using the whole ensemble of frames generated from each cleaned-up trajectory (Figure 2C).

The pairwise similarity of unbinding pathways projected on the low-D space built by Pysomap was then measured using the Fréchet distance, as implemented in the Similarity Measures library of R.85 In particular, the discrete Fréchet distance between the paths Sk and Sl was computed as follows:56 F(Sk , Sl) = inf max [d(Sk(α(t )), Sl(β(t )))] α , β t ∈[0,1]

(1)

In eq 1, the discrete Fréchet distance is defined as the infimum among all the parametrizations α and β of Sk and Sl of the maximum distance t between the same curves. The distance d is the Euclidean distance between points of the parametric curves Sk(α) and Sl(β) evaluated in the low-D space. Notably, the two parametric curves do not need to have the same length, and this makes this measure particularly suited to compare distinct unbinding trajectories. Thus, using the discrete Fréchet distance as a dissimilarity metric, all of the unbinding pathways were further analyzed via well-established cluster analysis methods. 539


Article

Journal of Chemical Information and Modeling Two clustering strategies, as implemented in R version 3.3.2, were employed depending on the aim of the analysis.75 For example, the agglomerative hierarchical clustering algorithm with complete linkage method was mainly exploited for visualization purposes. In agglomerative hierarchical clustering, similar objects (pathways in this specific case) are linked together to form growing clusters according to a bottom-up approach, where the order in which the clusters are merged is a function of their similarities. The hierarchical organization of the dissimilarity matrix can therefore be displayed as a tree diagram (dendrogram) which provides a simple visual interpretation of the data set structure. Conversely, a less subjective classification of pathways in similar groups, was performed through the k-means clustering method.86,87 With this method, each object is assigned to one of the k clusters defined by centroids, where the initial number of clusters k is a user-specified parameter. To determine the optimal number of clusters, we adopted the average silhouette method,88 as implemented in the Factoextra library of R, where the silhouette value is a measure of how similar an object is to its own cluster compared to the others. Thus, the optimal number of clusters k is defined as the one that maximizes the average silhouette over a range of possible values of k. Each cluster represents a set of similar unbinding trajectories, while the representative pathway of each cluster is defined as the closest to the centroid.

Figure 3. 2D kinetic map of the association rate constants (kon) versus dissociation rate constants (koff) with iso-affinity lines for the seven compounds employed in this study. Individual ligand affinities are shown in color scale.

ascribed to differences in affinity. From this standpoint, our data set can be considered as properly tailored to study relevant changes in unbinding kinetics. Even though the data set contains different scaffolds, thus displaying chemical diversity, the structures show good overlay and comparable orientation in their binding modes. When embedded in the ATP site, similar interactions are maintained by all considered compounds, most importantly the Asp93 Hbond that is one of the key interactions (Figure 1B). Moreover, as the compounds tend to occupy similar regions of the binding site, those possessing bulkier substituent groups are able to reach the deep hydrophobic portion of the pocket, comprising residues Met98, Leu103, Leu107, Phe138, Tyr139, Val150, and Trp162. As depicted in Figure 1, one moiety of the ligand is facing toward the pocket entrance and is therefore solvent exposed. This is a general observation for cocrystallized Hsp90 inhibitors, as reported in previous works on this target.57 Indeed, as the ATP binding site of Hsp90 is a relatively shallow pocket, most ligands, and particularly those presenting larger substituents, protrude into the bulk of the solvent. The initial states for our simulations were the protein− ligand complexes of the compounds listed in Table 1. For each protein−ligand complex, 25 sMD runs were carried out employing a scaling factor λ of 0.45. We note that this rather aggressive scaling factor was chosen as a best compromise between accuracy and speed. Indeed, while λ values closer to 1 better preserve the dynamics at 300 K, much longer simulations are required to witness unbinding events. Thus, the choice of a scaling factor of 0.45 was dictated by the need to observe repeated unbinding events within a submicrosecond time scale for all the investigated compounds. In order to preserve the overall secondary and tertiary structure of the protein, as this might be compromised by the scaling factor,13,38 weak positional restraints were applied to all the backbone heavy atoms, making exception for specific residues comprising the binding pocket (see Methods for details about the simulation setup). The computational residence times τcalc were estimated as the average of the exit times recorded for each compound (see Figure S1 for distributions of exit times). The results obtained are displayed in Table 2. First, we analyzed the differences in predictions for structurally similar compounds that showed different kinetic behavior. Thus, we focused separately on the quinazoline

■

RESULTS AND DISCUSSION Compound Ranking Based on the Computed Residence Times. We applied sMD to perform unbinding simulations for a series of drug-like inhibitors complexed to the N-terminal domain (NTD) of Hsp90.57 Here, Hsp90 was employed as a test case, while the seven investigated inhibitors were chosen because of the availability of experimental SPR measurements and because for five of them the crystal structure of the complex was available (compounds 2, 3, 4, 6, and 7). Concerning the remaining compounds (compounds 1 and 5), the bound states could be easily modeled in silico, because of their minor structural changes with respect to existing X-ray structures. While a simple replacement of a carbon atom by an oxygen was required to build compound 5 from compound 4, cleavage of a bulky substituent was needed to model compound 1 from 2. The latter procedure could be safely performed as the cleaved group does not establish major interactions with the protein in the corresponding crystal structure (see Figure 1A). Thus, because of the main interactions of the quinazoline scaffold are preserved upon molecular cleavage, we can safely assume that this replacement should not affect the native pose of compound 1. All the small molecules comprised in our data set, presented in Table 1, bind to the ATP binding site of the target protein Hsp90 (Figure 1A and B). As outlined, compounds 1 and 2 display a quinazoline scaffold, compounds 3−6 possess an indazole-like structure, and compound 7 is a pyrazole-containing, resorcinollike inhibitor. Both affinity (equilibrium dissociation constant, KD) and kinetic (association and dissociation rates, kon and koff respectively) properties are shown in Table 1. The corresponding residence times (determined as the inverse of the koff) range from a few seconds (compounds 1 and 7) up to thousands of seconds (compounds 2, 4, and 6). This makes the data set very diverse, both in terms of chemical structure and range of residence times. Moreover, as Figure 3 shows, differences in koff in our set of molecules cannot be entirely 540


Article


Table 2. Experimental and Computed Residence Times (τexp and τcalc, respectively) for Each of the Seven Compoundsa computed quantities compound

τexp ± σ

τcalc

±σ

±σe

±σBS

1 2 3 4 5 6 7

1.81 ± 0.18 1470.59 412.2 ± 13.43 1098.9 ± 320.64 75.76 ± 3.91 3402.25 ± 177.06 9.01

8.06 291.11 169.81 209.72 141.38 171.48 77.30

±9.68 ±259.04 ±150.50 ±115.53 ±133.01 ±138.89 ±44.88

±1.94 ±52.88 ±30.72 ±23.58 ±27.15 ±28.35 ±9.16

±1.77 ±49.96 ±29.20 ±22.53 ±25.86 ±26.61 ±8.69

Experimental residence times are reported in seconds, while computed residence times are reported in nanoseconds. σ stands for the standard deviation, and σe is for the standard error derived as σe = σ/√n where n equals the number of observations/measurements; σBS is the bootstrapped standard error. a

quinazoline analogs, for which the bulky moiety was plunged into the solvent, here the varying substituents lie buried within the deep pocket of the binding site, which is made up exclusively of hydrophobic residues (Figure 4). Therefore, the presence of these substituent groups is expected to have great influence in the process of ligand unbinding, as additional interactions need to be disrupted to allow for the ligands’ dissociation from the binding site. Compounds 3 and 4 display a difference in ring size, as a piperidine in the former is replaced by a pyrrolidine in the latter. Differently, both 5 and 6 possess an oxygen atom that increases the polarity of the moiety, which is part of the ring in compound 5 resulting in a morpholine, while it is present as a methoxy group attached to the piperidine in compound 6. Basing on the experimental residence times reported in Table 1, compound 5 is unbinding more than 300 s faster than 3, which in turn shows an almost 700 s shorter residence time with respect to compound 4. Compound 6 is the slowest unbinder, requiring over 2000 s more than compound 4 to dissociate. The separations obtained from our simulations for compounds 3, 4, and 5 were distributed accordingly, reproducing the ranking resulting from the SPR experiments. Indeed, as shown in Table 2, compound 3 was characterized by an average unbinding time that was about 30 ns longer than the one obtained for compound 5. Compound 4 was in the computation about 40 ns slower than compound 3 to dissociate from the protein. Despite the reliable reproduction of the experimental data for these drug-like chemotypes, it was not possible to rank compound 6 correctly. Within this series, compound 6 is undoubtedly the most structurally diverse with respect to the other three compounds. While it shares the piperidine moiety with compound 4, it also includes a methoxy group attached to the ring, which places a polar oxygen in a region that is relatively close to the one occupied by the in-ring oxygen in compound 5. From a structure−kinetics relationship standpoint, the presence of this methoxy substituent was predicted to accelerate dissociation of the ligand from the protein, similarly as for compound 5, possessing the in-ring oxygen, in contrast to experimental data. The polar oxygen in the morpholine ring (compound 5) however seems more exposed and therefore unfavorable in the hydrophobic pocket. Whereas, the rather buried oxygen in the methoxy substitution is more shielded, exposing the hydrophobic methyl group which is adapting to the hydrophobic pocket. In other words, in our simulations, 6 was predicted as an in-between situation of compounds 4 and 5. While this discrepancy may be ascribed to limitations of the force field in effectively discriminating such subtleties, it has to be kept in mind that this behavior might

analogs 1 and 2 and on the indazole compounds 3−6. As for the former, compound 2 includes a large substituent comprising a piperazine sulfonamide, which is charged at physiological pH. The significant dissimilarity with respect to compound 1, in which the entire substituent is absent, reflects in a pronounced difference in residence times. In fact, as shown in Table 1and Figure 3, the exit rates for compound 2 are significantly lower than those of compound 1. As mentioned before, in the bound state, the bulky group in compound 2 protrudes into the solvent and is not involved in relevant interactions with the protein (see Figure 1A and Figure 4).

Figure 4. Hydrophobic pocket. Focus on the binding mode of compound 2 is given: while maintaining the driving interaction with Asp93 (lower side of the picture, colors consistent with Figure 1B), the isoindoline moiety of the ligand is fully embedded inside a hydrophobic cavity comprising residues Met98, Leu103, Leu107, Phe138, Tyr139, Val150, Trp162 (in orange, stick representation); conversely, the bulky substituent on the quinazoline scaffold protrudes into the solvent, without engaging relevant interactions with the protein.

Nevertheless, it is reasonable to think that the presence of a variety of functional elements comprised in the substituent can be responsible for temporary interactions with protein residues, which might accelerate binding and/or slow down unbinding. Our simulations were able to distinguish the kinetic profile of compounds 1 and 2 very well, identifying a marked difference, of about 280 ns, between their computationally derived exit times (Table 2). Concerning the indazolecontaining chemicals, namely 3−6, the differences in terms of structure were less dramatic. However, contrary to the 541


Article


Figure 5. Normalized calculated residence times plotted against normalized experimental residence times. Experimental values were normalized λ i τ calc y τ exp according to τiexp , in which the scaling factor λ equals 0.45. Computational residence times were normalized according to jjj icalc zzz. Error bars refer max k τmax { to the standard deviations σ reported in Table 2. The trend line from simple linear regression is shown, along with the corresponding R2 values. (A) Results obtained including all of the seven considered compounds. Pearson correlation coefficients of R = 0.73 and R2 = 0.53 were obtained. (B) After excluding the single outlier, compound 6, the ranking for the remaining six compounds was reproduced correctly (Spearman rank correlation = 1). Accordingly, the R2 increased to 0.89.

( )

R of 0.94 and R2 = 0.89, likewise showing significant p-values (P < 0.01). Thus, except for the outlier compound 6, we can conclude that the prioritization obtained from our simulations was in good agreement to the kinetic experimental data (Figure 5B). Most important, even with a relatively aggressive scaling factor of 0.45, it was possible to clearly discern very fast binders from slow binders, further pointing to this approach as a suitable tool to prioritize compounds for subsequent medicinal chemistry campaigns. Similarity Analysis of Unbinding Pathways. We further investigated the sMD trajectories in order to retrieve relevant mechanistic details about the unbinding process and depict structure-kinetics relationships. In drug discovery, a major goal is identifying novel drug-like molecules that are able to bind and modulate the activity of a pharmacological target of interest. In this perspective, understanding which possible pathways can be followed to access the target binding site and the corresponding molecular features involved can be particularly key. Notably, multiple and complex routes can be involved in (un)binding processes. On one hand, the complexity varies as a function of the pocket shape, location, and nature of the residues making up the way leading to the binding site. On the other hand, additional complications are introduced by both the degrees of freedom of the ligand, and the variety of its functional groups, which can interact with the target protein along the binding pathway. Therefore, recognizing and classifying the possible routes becomes a nontrivial task. Herein, to face this challenge, we applied a strategy based on dimensionality reduction and cluster analysis. Specifically, we tracked the Cartesian coordinates of the ligands’ heavy atoms along the unbinding trajectories and used these points to embed a low-D (3D) space. To this end, we took advantage of Isomap, a nonlinear dimensionality reduction technique. Once the trajectories were projected in the low-D space, we evaluated their similarity calculating the Fréchet distance, which we employed as the metric for subsequent clustering. We applied this scheme to all the considered compounds (Table 1), made exception for compound 6 that, as discussed above, behaved as an outlier.

have been amplified by the choice of a low scaling factor (0.45). Furthermore, we examined structurally dissimilar compounds displaying comparable kinetic behavior. The smallest compounds in our data set, compounds 1 and 7, possess significantly diverse structures and scaffolds. Compound 7 lacks the substituent reaching into the deep hydrophobic pocket of the binding site, which can be found in all other compounds (Figure 1B). Despite the overall structural difference, these compounds displayed residence times in the same time scale from the SPR experiments. Accordingly, the two species exhibited the lowest computed residence times. The more structurally complex and diverse compounds 2 and 4 both required more but still comparable time to leave the protein binding site. Such behavior, expressed by experimental residence times in the order of 1000 s, was also confirmed by our computations, resulting in average unbinding times of more than 200 ns. Moreover, despite the profound structural differences, they were ranked correctly by the sMD runs, with compound 2 being predicted to take more time to leave the binding pocket. By gathering the data from the entire set of compounds (1− 7, Table 1), we evaluated the degree of correlation between simulations and experiments. To this aim, a linear regression procedure was applied, in which the effect of employing a scaling factor in the simulations was taken into account (see Figure 5). Including all of the compounds in the data set (Figure 5A) resulted in a Spearman ρ of 0.89, while p-values were significant (P < 0.05). Pearson’s product−moment correlation resulted in an R of 0.73 and R2 = 0.53. It is clear from the plot that the result obtained for compound 6 significantly affected the correlation. Therefore, we examined more thoroughly the results for this compound. Investigation of the studentized residuals, resulting in a value of t > 3, (t = −3.75), identified it as an outlier. Thus, the linear regression procedure was repeated excluding this compound. The Spearman rank correlation gave a ρ of 1, which indicates a correct ranking reproduction, while p-values were significant (P < 0.01). Pearson’s product−moment correlation resulted in an 542


Article


Figure 6. Dimensionality reduction and dendrogram interpretation of the unbinding trajectories. (A) Ligand configurations (gray points) in the low-D space. Projections of sample trajectories (pathways 3 and 16, in orange and yellow; pathways 14 and 18 in light blue and dark blue) are highlighted as colored lines and points. (B) Heatmap of the Fréchet distances between the 25 unbinding trajectories (labels go from 1 to 25). The color scale indicates high similarity in white (low values) and low similarity in blue (higher values). The result of the hierarchical cluster analysis on this matrix is represented on top and on the side of the heatmap in form of a dendrogram. We note that the unbinding pathways are reordered in the Fréchet matrix in such a way to reflect their relative similarity in the dendrogram. (C left panel) Ligand COM coordinates are tracked along the unbinding route in the high-D space from the initial bound pose (black dot) to the bulk solvent. Two sample pathways, 14 and 18, are represented. (right panel) Sample configurations visited during the unbinding process. Colors are consistent to section A. (D) Same scheme indicated for section C, but relative to trajectories 3 and 16. MD simulation data of compound 1 were used to produce the figure.

the left-hand panels of Figure 6C and D, the COM (center of mass) position during the unbinding event in the high-D space was tracked for pathway pairs 14 and 18 and 3 and 16, respectively (colors are consistent to those used in Figure 6A). The right-hand panels of Figure 6C and D show instead sample configurations (I and II) visited along the process. Figure 6C clearly shows that the configurations (I and II) for each of the two trajectories (14 and 18, light blue and dark blue, respectively) overlay almost perfectly. We can observe a similar pattern for pathways 3 and 16, as shown in Figure 6D. Therefore, the low Fréchet distance resulted as a suitable parameter to successfully identifying high similarity within a set of unbinding trajectories. Moreover, by comparing Figures 6C and D, it can be seen that the ligand is found in remarkably different orientations (blue scale and orange scale configurations, respectively). Indeed, pair 14 and 18 and pair 3 and 16 belong to different clusters, as indicated by the dendrograms (Figure 6B). Importantly, even though the COM trace in the Cartesian space is similar for pathways 3, 16 and 14, 18 (left-hand panel in Figures 6C and D), these two clusters of trajectories are highly separated in the low-D space, reflecting their mechanistic dissimilarity. Indeed, while pathways 14 and 18 depict a scenario where the quinazoline ring unbinds first, followed by the rest of the molecule, the opposite situation is portrayed by pathways 3 and 16. Monitoring and comparing individual configurations visited in different runs would not be manageable when handling large sets of trajectories. The similarities/differences were instead relatively straightforward to determine through our procedure

Both, nonlinear multidimensional scaling and the Fréchet distance have been extensively used to analyze MD trajectories.89,90 However, to the best of our knowledge, this is the first time they have been used in combination and applied to identify path similarities for the unbinding of druglike molecules. In Figure 6, we provide an illustration of the strategy devised, using the data obtained from compound 1. First, the relevant configurations visited by the ligand during unbinding in all of the performed trajectories were used as an input for the dimensionality reduction procedure (see Figure 2 for a graphical illustration of the procedure). As a result, a low-D space was obtained, made up of the projection of the initial points in this new space (Figure 6A, gray points). Individual trajectories could then be followed in this new space with reduced dimensionality (Figure 6A, colored lines and points). Subsequently, the mutual distances between all of the trajectories available for the ligand (25 in the present case) were evaluated and stored in a squared matrix, shown as a heatmap in Figure 6B. A hierarchical cluster analysis (see Methods) was then performed on the obtained matrix in order to identify similar pathways (Figure 6B, dendrogram on top and on the side of the heatmap). For instance, according to Figure 6B, the path pair 3 and 16 resulted to be very similar for compound 1. This similarity could be visualized and confirmed in the low-D space (Figure 6A, pathways in orange scale). Similar reasoning holds for paths 14 and 18 (Figure 6A, blue scale). As a further confirmation, we inspected the configurations sampled by these trajectories in the high-D space. In 543


Article


Figure 7. Significantly different exit route visited by compound 5. (A) In the low-D space, most of the data (gray points) are concentrated in the same region, corresponding to the ligand leaving the binding site at the solvent-exposed side of the pocket (sample pathways in green scale). A less populated agglomeration of points was evident, highlighted in red by the corresponding trajectory. (B) Ligand COM coordinates are tracked along the unbinding route from the initial bound pose (black dot) to the bulk in the high-D. In red, the unusual passage under α-helix3 is highlighted while the green points represent the most conventional unbinding pathways.

12.17 ns). The trajectory analysis suggests that pathways characterized by an early exposure of the quinazoline ring to the solvent, with a more pronounced insertion of the dihydroisoindole in the hydrophobic pocket (cluster 1), are intrinsically faster dissociating than those for which the quinazoline ring leaves later (cluster 2). Such observation implies that, in order to increase the residence time of this compound, one might operate by reducing the relative population of trajectories represented by cluster 1 in favor of those belonging to cluster 2. This can be obtained either by stabilizing the interactions of the quinazoline with proper substituents or by increasing the lipophilicity of the dihydro-isoindole inside the deep hydrophobic pocket. Both options would lead to a stabilization of the metastable states encountered along the cluster 1 pathway and, ultimately, slow down the entire unbinding process. A similar reasoning applies to compound 2, where the most populated cluster (cluster 3) accounts for the slowest unbinding pathways (tav = 490.63 ns). If we compare the representative pathway of this cluster (run 20, also associated with the lowest individual exit time of 931.93 ns) with the one belonging to the fastest unbinding trajectories (cluster 2, representative member: run 25), it is possible to notice that similar configurations can be identified along the egress process. In particular, as reported in Figure 8A, both pathways show a ligand configuration in which the dihydro-isoindole is deeply buried inside the hydrophobic pocket, while the polar moiety of the molecule reaches toward the solvent exposed part of the binding site. While a salt bridge between compound 2 and Asp54 can be established in the representative member of cluster 3, this interaction was never found in cluster 2. We therefore hypothesize that this specific ionic interaction leads to higher residence times. Its preferential formation could be intentionally exploited in the optimization phase of drug design programs. While gathering kinetic information from a series of unbinding trajectories of the same ligand is certainly interesting, in the context of a drug design process it is often more relevant to compare the kinetic behavior of structurally related compounds. The data set employed in this study allows a direct comparison of the two structurally related drug-like molecules, compounds 4 and 5, where the indazole scaffold is decorated with a piperidine substituent in the former and a

in the low-D projection. Interestingly, through the analysis of the low-D space for all the considered ligands, we were furthermore able to quickly identify significantly different exit routes. While in most of the trajectories the ligand reaches the bulk solvent through the solvent-exposed entrance of the binding site, a passage under α-helix3 was incidentally followed (Figure 7A). In fact, while the vast majority of configurations are concentrated in the same volume (Figure 7A, gray points), a region separated from the rest, corresponding to the mentioned alternative path (Figure 7A and B, red line and points) could be observed. The heatmaps of the Fréchet distances for compounds 2−5 and 7 are reported in Figure S2. Potential Implications for Drug Design. The different average unbinding times recorded (Table 2) for the considered ligands were the result of protein−ligand complex disruptions which took place in different ranges of simulation times. Taken together with the pathway classifications that can be achieved through the proposed strategy, another interesting aspect to investigate was a potential correlation between unbinding routes and residence times. It is indeed reasonable to think that, as long as the kinetic model provided by sMD is able to correctly rank the computed residence times for the investigated compounds, at least part of this information should reflect mechanistic aspects of the unbinding pathways. We therefore examined our trajectories to identify relevant mechanistic and structural features of ligand unbinding. Even though those features could only provide some general indications due to the relatively aggressive scaling factor, it is nonetheless interesting to gather all the potentially relevant information for kinetic-based drug design. In Table S3, we report the results of k-means clustering on the unbinding pathways together with the exit times averaged over all the members of each cluster (tav). As the table shows, while no clear-cut correlation between geometric features and temporal behavior was found, a certain partition of pathways consistent with similar time scales could be inferred in some cases. This is particularly interesting for compound 1, for which the most populated cluster (cluster no. 1) is associated with the fastest unbinding trajectories (tav = 5.16 ns). Notably, the previously discussed unbinding pathways 14 and 18 (see Figure 6) are found within this cluster. In contrast, pathways 3 and 16 are grouped in cluster no. 2, which corresponds to the highest average exit time recorded for this compound (tav = 544


Article


only on single ligands, rather than a set of compounds, diminishing the drug discovery-oriented appeal of the proposed methodology. Comparison with Related Methods. The sMD method has been largely employed in both retrospective13,49 and prospective studies50 and can be nowadays considered as a well-established methodology for compound ranking based on computed residence times. However, while this work was already at an advanced stage of development, a similar approach based on RAMD was presented (τRAMD) and applied to the same pharmaceutical target investigated in this study.40 It is therefore interesting to compare the main features of the two methodologies which undoubtedly share some conceptual overlap. Indeed, while in sMD the effective temperature of the simulated system is enhanced by rescaling the potential energy function provided by the force field, in τRAMD an additional force with random orientation is applied to the ligand. In this way, the egress of compounds from their binding site is facilitated, and the computed residence time, defined as the simulation time required for ligand dissociation in half of the runs, can be obtained at affordable time scales. Thus, both methods share the presence of a tunable parameter (the scaling factor for sMD and the magnitude of the force for τRAMD) whose choice is critical in determining the accuracyto-efficiency ratio of compound ranking. Both the scaling factor and the magnitude of the force must be chosen in a way to observe enough unbinding events to ensure acceptable statistics for all the compounds in the data set at an affordable computational time. The number of unbinding events is typically set to 25 in case of sMD, whereas in τRAMD a more statistically oriented procedure is employed to reduce the uncertainty of results. In the work by Kokh et al., a total number of RAMD simulations ranging from 40 to 200 per compound was employed.40 That being said, the dissociation time for τRAMD was on the order of 1−10 ns time scales (with a magnitude of the random force of 14 kcal mol−1 Å−1), whereas in our case up to almost 1 μs of simulation was required to observe a single unbinding event for compound 2 at the relatively aggressive scaling factor of 0.45. Notably, this difference in computational efficiency allowed Kokh et al. to investigate a 1 order of magnitude larger data set compared to ours. A side effect of using a scaled potential energy surface is that in sMD all the degrees of freedom are accelerated, requiring to preserve the secondary structure of the protein with additional restraints. Moreover, at the same time, enough conformational freedom must be allowed to the binding site in order not to interfere with the accelerated dynamics of the ligand. While a standard recipe to satisfy this requirement is available,13 it is important to note that τRAMD is devoid from this potential source of complication. We therefore conclude that both methods are well-equipped to provide an efficient and reliable estimation of relative residence time for drugtarget complexes. Further investigation, possibly involving more difficult test cases, will ultimately be required to assess the relative merits of these important tools for drug discovery and development.

Figure 8. (A) Similar configurations extracted from representative pathways of the most populated cluster (cluster 3, run 20, cyan sticks) and the cluster associated with the fastest average exit time (cluster 2, run 25, pink sticks) for compound 2. The additional salt bridge established by compound 2 and Asp54 is only observed in the representative pathway of cluster 3. (B) Comparison of configurations extracted from the representative pathways of the most populated clusters for compounds 4 and 5 (cyan and pink sticks, respectively).

morpholine group in the latter. This substitution is associated with 1 order of magnitude higher affinity in favor of compound 4 and a longer residence time (Table 1). By comparing the representative pathways of the most populated clusters for the two compounds (cluster 2, run 23 and cluster 3, run 22 for 4 and 5, respectively) it is possible to notice that in both cases the ligand protrudes further into the hydrophobic pocket (Figure 8B) before leaving the binding site. However, this configuration would correspond to more favorable contacts inside the hydrophobic pocket only for compound 4, possessing a more hydrophobic ring than compound 5. This gives rise to the longer residence time determined for compound 4. Again, the analysis suggests that more hydrophobic, and possibly bulkier, substituents might lead to a further decrease of koff for this molecular scaffold. We wish to underline that the previously reported discussion is partly based on a speculative interpretation of unbinding pathways. Indeed, further simulations, possibly at higher level of computation, would be required to fully address the aspects highlighted by our trajectory analysis and to detect statistically significant differences in the mean exit times between distinct clusters. Approaching the problem by employing a more conservative λ, and possibly relying on more statistics of unbinding events, would provide a better suited set of data to investigate the presence of correlations between geometric features of trajectories and unbinding time scales. However, this would translate into larger requirements in terms of computational resources and most likely would imply focusing

■

CONCLUSION AND FUTURE PERSPECTIVE The importance of taking into account the kinetic behavior of drugs, in addition to their affinity, has become a widely accepted concept in drug discovery. In spite of this, implementing the measurements of kinetic parameters in the early stages of the drug discovery process has proven to be 545


Article

Journal of Chemical Information and Modeling Present Address

difficult. While in silico methods are increasingly employed to assist the whole drug discovery process, so far, no wellestablished computational methods are available to routinely predict drugs’ off-rates. Here, we investigated the capability of scaled molecular dynamics simulations to predict the residence time and the unbinding routes for a series of Hsp90 inhibitors. We challenged the methodology by examining very diverse structures of lead-like molecules, as it might be the case in a real drug discovery scenario. Notably, the experimental residence times of the examined compounds covered a range of 3 orders of magnitude, spanning from less than 2 s up to almost 1 h. To the best of our knowledge, this is the data set exhibiting the largest range in residence time examined with scaled MD simulations up to now. First, we have shown that structurally similar compounds could be ranked reliably, even though one inhibitor turned out to be an outlier of our kinetic model. Structurally highly similar compounds displaying distinct kinetic differences proved difficult to differentiate in single runs. However, statistical significance and robustness of the observed events was assessed by bootstrapping analysis, and it was shown that even very subtle changes in chemical structures could be ranked in good agreement with experimental evidence. Besides, the proper kinetic behavior could be reproduced despite large structural changes. Then, we have devised a fully automated data analysis approach to detect similar features on the ensemble of unbinding trajectories. Despite the aggressive scaling factor employed in this work, the methodology seems to be promising. By integrating this information with the observed unbinding times, recurring dynamical features implicated in determining characteristic residence times could be ultimately identified. Finally, the information obtained can be further exploited to guide more accurate sampling strategies, which purposely aim at identifying minimum free energy paths and the associated potentials of mean force. Indeed, the low-D embedding can be directly used as a collective variable space to feed biased-MD simulations. This would lead to characterize more reliably the mechanism of the process, to quantify the involved energy barriers, and to identify accurately the relevant metastable states encountered during unbinding. This latter aspect of the data analysis will be further investigated in following studies.

■

○

IRIC−Institut de Recherche en Immunologie et en Cancérologie, Université de Montréal, 2950 Chemin de Polytechnique, Marcelle-Coutu Pavilion, Montréal, QC H3T 1J4, Canada. Author Contributions #

These authors contributed equally. D.A.S. and M.Bernetti performed the simulations. D.A.S., M.M., M.Bernetti, and M.Bertazzo analyzed data. D.M. solved crystal structure PDB 6HHR. H.-M.E. developed Hsp90 inhibitors described in this article. M.M., G.E., and A.C. supervised the studies. D.A.S., M.Bernetti, M.M., and M.R. wrote the manuscript. Notes

The authors declare the following competing financial interest(s): A.C. is a co-founder of BiKi Technologies, a startup company that develops methods based on molecular dynamics and related approaches for investigating proteinligand (un)binding.

■

ACKNOWLEDGMENTS This work was supported by the EU/EFPIA Innovative Medicines Initiative (IMI) Joint Undertaking, K4DD (grant no. 1115366). This paper reflects only the authors’ views, and neither the IMI nor the European Commission is liable for any use that may be made of the information contained herein. The computational results presented have been achieved in part using the Vienna Scientific Cluster (VSC). D.A.S. thanks Sergio Decherchi for technical support with BiKi Netics, Riccardo Martini for help with Gimp, and Prof. Stefan Boresch for MD discussions and useful suggestions. Furthermore, she thanks Prof. Thierry Langer for providing the HP cluster “Hydra” for many of the simulations presented here. We wish to thank Dario Gioia for useful discussions regarding the analysis of unbinding pathways and for careful reading of the manuscript. Paolo Cinelli is also acknowledged for technical support.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jcim.8b00614. Data collection and refinement statistic of X-ray structure (PDF) Accession Codes

The atomic coordinates have been deposited in the ProteinData Bank (PDB code 6HHR for compound 7).

■

REFERENCES

(1) Schuetz, D. A.; de Witte, W. E. A.; Wong, Y. C.; Knasmueller, B.; Richter, L.; Kokh, D. B.; Sadiq, S. K.; Bosma, R.; Nederpelt, I.; Heitman, L. H.; Segala, E.; Amaral, M.; Guo, D.; Andres, D.; Georgi, V.; Stoddart, L. A.; Hill, S.; Cooke, R. M.; De Graaf, C.; Leurs, R.; Frech, M.; Wade, R. C.; de Lange, E. C. M.; IJzerman, A. P.; MüllerFahrnow, A.; Ecker, G. F. Kinetics for Drug Discovery: An IndustryDriven Effort to Target Drug Residence Time. Drug Discovery Today 2017, 22, 896−911. (2) Spagnuolo, L. A.; Eltschkner, S.; Yu, W.; Daryaee, F.; Davoodi, S.; Knudson, S. E.; Allen, E. K. H.; Merino, J.; Pschibul, A.; Moree, B.; Thivalapill, N.; Truglio, J. J.; Salafsky, J.; Slayden, R. A.; Kisker, C.; Tonge, P. J. Evaluating the Contribution of Transition-State Destabilization to Changes in the Residence Time of TriazoleBased InhA Inhibitors. J. Am. Chem. Soc. 2017, 139, 3417−3429. (3) Lu, H.; Tonge, P. J. Drug-Target Residence Time: Critical Information for Lead Optimization. Curr. Opin. Chem. Biol. 2010, 14, 467−474. (4) Copeland, R. A. The Drug-Target Residence Time Model: A 10Year Retrospective. Nat. Rev. Drug Discovery 2016, 15, 87−95. (5) Bernetti, M.; Cavalli, A.; Mollica, L. Protein−Ligand (Un) Binding Kinetics as a New Paradigm for Drug Discovery at the Crossroad between Experiments and Modelling. MedChemComm 2017, 8, 534−550. (6) Pan, A. C.; Borhani, D. W.; Dror, R. O.; Shaw, D. E. Molecular Determinants of Drug-Receptor Binding Kinetics. Drug Discovery Today 2013, 18, 667−673.

AUTHOR INFORMATION

Corresponding Authors

*Email: [email protected] (M.M.). *Email: [email protected] (A.C.). ORCID

Doris A. Schuetz: 0000-0001-8615-9554 Matteo Masetti: 0000-0002-3757-7802 Gerhard F. Ecker: 0000-0003-4209-6883 Andrea Cavalli: 0000-0002-6370-1176 546


Article

Journal of Chemical Information and Modeling (7) Schoop, A.; Dey, F. On-Rate Based Optimization of StructureKinetic Relationship–Surfing the Kinetic Map. Drug Discovery Today: Technol. 2015, 17, 9−15. (8) Copeland, R. A.; Pompliano, D. L.; Meek, T. D. Drug−Target Residence Time and Its Implications for Lead Optimization. Nat. Rev. Drug Discovery 2006, 5, 730−739. (9) Klebe, G. Applying Thermodynamic Profiling in Lead Finding and Optimization. Nat. Rev. Drug Discovery 2015, 14, 95. (10) Kruse, A. C.; Hu, J.; Pan, A. C.; Arlow, D. H.; Rosenbaum, D. M.; Rosemond, E.; Green, H. F.; Liu, T.; Chae, P. S.; Dror, R. O.; et al. Structure and Dynamics of the M3Muscarinic Acetylcholine Receptor. Nature 2012, 482, 552−556. (11) Sykes, D. A.; Moore, H.; Stott, L.; Holliday, N.; Javitch, J. A.; Lane, J. R.; Charlton, S. J. Extrapyramidal Side Effects of Antipsychotics Are Linked to Their Association Kinetics at Dopamine D 2 Receptors. Nat. Commun. 2017, 8, 763. (12) Walkup, G. K.; You, Z.; Ross, P. L.; Allen, E. K. H.; Daryaee, F.; Hale, M. R.; O’Donnell, J.; Ehmann, D. E.; Schuck, V. J. A.; Buurman, E. T.; Choy, A. L.; Hajec, L.; Murphy-Benenato, K.; Marone, V.; Patey, S. A.; Grosser, L. A.; Johnstone, M.; Walker, S. G.; Tonge, P. J.; Fisher, S. L. Translating Slow-Binding Inhibition Kinetics into Cellular and in Vivo Effects. Nat. Chem. Biol. 2015, 11, 416−423. (13) 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. (14) Klebe, G. The Use of Thermodynamic and Kinetic Data in Drug Discovery: Decisive Insight or Increasing the Puzzlement? ChemMedChem 2015, 10, 229−231. (15) Gioia, D.; Bertazzo, M.; Recanatini, M.; Cavalli, A.; Masetti, M. Dynamic Docking : A Paradigm Shift in Computational Drug Discovery. Molecules 2017, 22, 1−21. (16) Faulon, J.-L.; Bender, A. Handbook of Chemoinformatics Algorithms; CRC press, 2010. (17) De Vivo, M.; Masetti, M.; Bottegoni, G.; Cavalli, A. Role of Molecular Dynamics and Related Methods in Drug Discovery. J. Med. Chem. 2016, 59, 4035−4061. (18) Swinney, D. C.; Beavis, P.; Chuang, K.-T.; Zheng, Y.; Lee, I.; Gee, P.; Deval, J.; Rotstein, D. M.; Dioszegi, M.; Ravendran, P.; Zhang, J.; Sankuratri, S.; Kondru, R.; Vauquelin, G. A Study of the Molecular Mechanism of Binding Kinetics and Long Residence Times of Human CCR5 Receptor Small Molecule Allosteric Ligands. Br. J. Pharmacol. 2014, 171, 3364−3375. (19) Bull, H. G.; Garcia-Calvo, M.; Andersson, S.; Baginsky, W. F.; Chan, H. K.; Ellsworth, D. E.; Miller, R. R.; Stearns, R. A.; Bakshi, R. K.; Rasmusson, G. H.; Tolman, R. L.; Myers, R. W.; Kozarich, J. W.; Harris, G. S. Mechanism-Based Inhibition of Human Steroid 5αReductase by Finasteride: Enzyme-Catalyzed Formation of NADP− Dihydrofinasteride, a Potent Bisubstrate Analog Inhibitor. J. Am. Chem. Soc. 1996, 118, 2359−2365. (20) Shaw, D. E.; Chao, J. C.; Eastwood, M. P.; Gagliardo, J.; Grossman, J. P.; Ho, C. R.; Lerardi, D. J.; Kolossváry, I.; Klepeis, J. L.; Layman, T.; McLeavey, C.; Deneroff, M. M.; Moraes, M. A.; Mueller, R.; Priest, E. C.; Shan, Y.; Spengler, J.; Theobald, M.; Towles, B.; Wang, S. C.; Dror, R. O.; Kuskin, J. S.; Larson, R. H.; Salmon, J. K.; Young, C.; Batson, B.; Bowers, K. J. Anton, a Special-Purpose Machine for Molecular Dynamics Simulation. Commun. ACM 2008, 51, 91−97. (21) 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, 2200−2205. (22) Doerr, S.; De Fabritiis, G. On-the-Fly Learning and Sampling of Ligand Binding by High-Throughput Molecular Simulations. J. Chem. Theory Comput. 2014, 10, 2064−2069. (23) Martínez-Rosell, G.; Giorgino, T.; Harvey, M. J.; de Fabritiis, G. Drug Discovery and Molecular Dynamics: Methods, Applications and Perspective Beyond the Second Timescale. Curr. Top. Med. Chem. 2017, 17, 2617−2625.

(24) Spitaleri, A.; Decherchi, S.; Cavalli, A.; Rocchia, W. Fast Dynamic Docking Guided by Adaptive Electrostatic Bias: The MDBinding Approach. J. Chem. Theory Comput. 2018, 14, 1727−1736. (25) Martinez-Rosell, G.; Harvey, M. J.; De Fabritiis, G. MolecularSimulation-Driven Fragment Screening for the Discovery of New CXCL12 Inhibitors. J. Chem. Inf. Model. 2018, 58, 683−691. (26) Sabbadin, D.; Moro, S. Supervised Molecular Dynamics (SuMD) as a Helpful Tool to Depict GPCR-Ligand Recognition Pathway in a Nanosecond Time Scale. J. Chem. Inf. Model. 2014, 54, 372−376. (27) Bertazzo, M.; Bernetti, M.; Recanatini, M.; Masetti, M.; Cavalli, A. Fully Flexible Docking via Reaction-Coordinate-Independent Molecular Dynamics Simulations. J. Chem. Inf. Model. 2018, 58, 490−500. (28) Amaro, R. E.; Mulholland, A. J. Multiscale Methods in Drug Design Bridge Chemical and Biological Complexity in the Search for Cures. Nat. Rev. Chem. 2018, 2, 148. (29) Votapka, L. W.; Jagger, B. R.; Heyneman, A. L.; Amaro, R. E. SEEKR: Simulation Enabled Estimation of Kinetic Rates, a Computational Tool to Estimate Molecular Kinetics and Its Application to Trypsin−Benzamidine Binding. J. Phys. Chem. B 2017, 121, 3597− 3606. (30) Jagger, B. R.; Lee, C. T.; Amaro, R. E. Quantitative Ranking of Ligand Binding Kinetics with a Multiscale Milestoning Simulation Approach. J. Phys. Chem. Lett. 2018, 9, 4941−4948. (31) Yang, L.; Liu, C.-W.; Shao, Q.; Zhang, J.; Gao, Y. Q. From Thermodynamics to Kinetics: Enhanced Sampling of Rare Events. Acc. Chem. Res. 2015, 48, 947−955. (32) 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. (33) Tiwary, P.; Parrinello, M. From Metadynamics to Dynamics. Phys. Rev. Lett. 2013, 111, 230602. (34) 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, E386−E391. (35) Tiwary, P.; Mondal, J.; Berne, B. J. How and When Does an Anticancer Drug Leave Its Binding Site? Sci. Adv. 2017, 3, e1700014. (36) 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, 4780−4788. (37) Hansmann, U. H. E. Parallel Tempering Algorithm for Conformational Studies of Biological Molecules. Chem. Phys. Lett. 1997, 281, 140−150. (38) Sinko, W.; Miao, Y.; de Oliveira, C. A. F.; McCammon, J. A. Population Based Reweighting of Scaled Molecular Dynamics. J. Phys. Chem. B 2013, 117, 12759−12768. (39) Lüdemann, S. K.; Lounnas, V.; Wade, R. C. How Do Substrates Enter and Products Exit the Buried Active Site of Cytochrome P450cam? 1. Random Expulsion Molecular Dynamics Investigation of Ligand Access Channels and Mechanisms11Edited by J. Thornton. J. Mol. Biol. 2000, 303, 797−811. (40) Kokh, D. B.; Amaral, M.; Bomke, J.; Grädler, U.; Musil, D.; Buchstaller, H.-P.; Dreyer, M. K.; Frech, M.; Lowinski, M.; Vallee, F.; Bianciotto, M.; Rak, A.; Wade, R. C. Estimation of Drug-Target Residence Times by τ-Random Acceleration Molecular Dynamics Simulations. J. Chem. Theory Comput. 2018, 14, 3859−3869. (41) 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, 2983−2989. (42) Dickson, A.; Brooks, C. L. WExplore: Hierarchical Exploration of High-Dimensional Spaces Using the Weighted Ensemble Algorithm. J. Phys. Chem. B 2014, 118, 3532−3542. (43) Dickson, A.; Lotz, S. D. Ligand Release Pathways Obtained with WExplore: Residence Times and Mechanisms. J. Phys. Chem. B 2016, 120, 5377−5385. 547


Article


(64) Gaussian 03; Gaussian, Inc.: Wallingford CT, 2004. (65) Nosé, S.; Klein, M. l. Constant Pressure Molecular Dynamics for Molecular Systems. Mol. Phys. 1983, 50, 1055−1076. (66) Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52, 7182−7190. (67) Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R. GROMACS: A Message-Passing Parallel Molecular Dynamics Implementation. Comput. Phys. Commun. 1995, 91, 43−56. (68) 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, 435−447. (69) Van Der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. GROMACS: Fast, Flexible, and Free. J. Comput. Chem. 2005, 26, 1701−1718. (70) Swope, W. C.; Andersen, H. C.; Berens, P. H.; Wilson, K. R. A Computer Simulation Method for the Calculation of Equilibrium Constants for the Formation of Physical Clusters of Molecules: Application to Small Water Clusters. J. Chem. Phys. 1982, 76, 637− 649. (71) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577−8593. (72) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. LINCS : A Linear Constraint Solver for Molecular Simulations. J. Comput. Chem. 1997, 18, 1463−1472. (73) Miyamoto, S.; Kollman, P. A. Settle: An Analytical Version of the SHAKE and RATTLE Algorithm for Rigid Water Models. J. Comput. Chem. 1992, 13, 952−962. (74) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33−38. (75) RStudio−Open Source and Enterprise-Ready Professional Software for R; RStudio, Inc.: Boston, MA, 2018. (76) Li, W.; Ma, A. Recent Developments in Methods for Identifying Reaction Coordinates. Mol. Simul. 2014, 40, 784−793. (77) Hartmann, C.; Banisch, R.; Sarich, M.; Badowski, T.; Schütte, C. Characterization of Rare Events in Molecular Dynamics. Entropy 2014, 16, 350−376. (78) Tenenbaum, J. B.; de Silva, V.; Langford, J. C. A Global Geometric Framework for Nonlinear Dimensionality Reduction. Science 2000, 290, 2319−2323. (79) Python Library for Isometric Feature Mapping (Isomap) https://web.vscht.cz/~spiwokv/pysomap/index.html (accessed May 23, 2018). (80) Rajan, A.; Freddolino, P. L.; Schulten, K. Going beyond Clustering in MD Trajectory Analysis: An Application to Villin Headpiece Folding. PLoS One 2010, 5, e9890. (81) Spiwok, V.; Králová, B. Metadynamics in the Conformational Space Nonlinearly Dimensionally Reduced by Isomap. J. Chem. Phys. 2011, 135, 224504. (82) Bonomi, M.; Branduardi, D.; Bussi, G.; Camilloni, C.; Provasi, D.; Raiteri, P.; Donadio, D.; Marinelli, F.; Pietrucci, F.; Broglia, R. A.; Parrinello, M. PLUMED: A Portable Plugin for Free-Energy Calculations with Molecular Dynamics. Comput. Phys. Commun. 2009, 180, 1961−1972. (83) Das, P.; Moll, M.; Stamati, H.; Kavraki, L. E.; Clementi, C. Low-Dimensional, Free-Energy Landscapes of Protein-Folding Reactions by Nonlinear Dimensionality Reduction. Proc. Natl. Acad. Sci. U. S. A. 2006, 103, 9885−9890. (84) Banisch, R.; Vanden-Eijnden, E. Direct Generation of LoopErased Transition Paths in Non-Equilibrium Reactions. Faraday Discuss. 2016, 195, 443−468. (85) R: The R Project for Statistical Computing; May 30, 2018. (86) Lloyd, S. Least Squares Quantization in PCM. IEEE Trans. Inf. Theory 1982, 28, 129−137. (87) MacQueen, J. Some Methods for Classification and Analysis of Multivariate Observations. In Proceedings of the fifth Berkeley

(44) Dickson, A.; Lotz, S. D. Multiple Ligand Unbinding Pathways and Ligand-Induced Destabilization Revealed by WExplore. Biophys. J. 2017, 112, 620−629. (45) 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, 618−628. (46) Dixon, T.; Lotz, S. D.; Dickson, A. Predicting Ligand Binding Affinity Using On- and off-Rates for the SAMPL6 SAMPLing Challenge. J. Comput.-Aided Mol. Des. 2018, 32, 1001. (47) Frank, A. T.; Andricioaei, I. Reaction Coordinate-Free Approach to Recovering Kinetics from Potential-Scaled Simulations: Application of Kramers’ Rate Theory. J. Phys. Chem. B 2016, 120, 8600−8605. (48) Tsujishita, H.; Moriguchi, I.; Hirono, S. Potential-Scaled Molecular Dynamics and Potential Annealing: Effective Conformational Search Techniques for Biomolecules. J. Phys. Chem. 1993, 97, 4416−4420. (49) 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, 7167−7176. (50) Bernetti, M.; Rosini, E.; Mollica, L.; Masetti, M.; Pollegioni, L.; Recanatini, M.; Cavalli, A. Binding Residence Time through Scaled Molecular Dynamics: A Prospective Application to HDAAO Inhibitors. J. Chem. Inf. Model. 2018, 58, 2255. (51) Schopf, F. H.; Biebl, M. M.; Buchner, J. The HSP90 Chaperone Machinery. Nat. Rev. Mol. Cell Biol. 2017, 18, 345−360. (52) Prodromou, C. Regulatory Mechanisms of Hsp90. Biochem. Mol. Biol. J. 2017, 3, 2. (53) Stebbins, C. E.; Russo, A. A.; Schneider, C.; Rosen, N.; Hartl, F. U.; Pavletich, N. P. Crystal Structure of an Hsp90−Geldanamycin Complex: Targeting of a Protein Chaperone by an Antitumor Agent. Cell 1997, 89, 239−250. (54) Schuetz, D. A.; Seidel, T.; Garon, A.; Martini, R.; Körbel, M.; Ecker, G. F.; Langer, T. GRAIL: GRids of PhArmacophore Interaction FieLds. J. Chem. Theory Comput. 2018, 14, 4958−4970. (55) Rohrdanz, M. A.; Zheng, W.; Clementi, C. Discovering Mountain Passes via Torchlight: Methods for the Definition of Reaction Coordinates and Pathways in Complex Macromolecular Reactions. Annu. Rev. Phys. Chem. 2013, 64, 295−316. (56) Alt, H.; Godau, M. Computing the Fréchet Distance between Two Polygonal Curves. Int. J. Comput. Geom. Appl. 1995, 05, 75−91. (57) Schuetz, D. A.; Richter, L.; Amaral, M.; Grandits, M.; Grädler, U.; Musil, D.; Buchstaller, H.-P.; Eggenweiler, H.-M.; 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, 4397−4411. (58) Decherchi, S.; Bottegoni, G.; Spitaleri, A.; Rocchia, W.; Cavalli, A. BiKi Life Sciences: A New Suite for Molecular Dynamics and Related Methods in Drug Discovery. J. Chem. Inf. Model. 2018, 58, 219−224. (59) 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, 926−935. (60) Lindorff-Larsen, K.; Piana, S.; Palmo, K.; Maragakis, P.; Klepeis, J. L.; Dror, R. O.; Shaw, D. E. Improved Side-Chain Torsion Potentials for the Amber Ff99SB Protein Force Field. Proteins: Struct., Funct., Genet. 2010, 78, 1950−1958. (61) 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, 247−260. (62) 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, 1157−1174. (63) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Kollman, P. A. Application of RESP Charges to Calculate Conformational Energies, Hydrogen Bond Energies, and Free Energies of Solvation. J. Am. Chem. Soc. 1993, 115, 9620−9631. 548


Article

Journal of Chemical Information and Modeling symposium on mathematical statistics and probability; University of California Press: Berkeley, CA, USA., 1967; Vol. 1, pp 281−297. (88) Kaufman, L.; Rousseeuw, P. J. Finding Groups in Data: An Introduction to Cluster Analysis; John Wiley & Sons, 2009; Vol. 344. (89) Seyler, S. L.; Kumar, A.; Thorpe, M. F.; Beckstein, O. Path Similarity Analysis: A Method for Quantifying Macromolecular Pathways. PLoS Comput. Biol. 2015, 11, e1004568. (90) Rydzewski, J.; Nowak, W. Machine Learning Based Dimensionality Reduction Facilitates Ligand Diffusion Paths Assessment: A Case of Cytochrome P450cam. J. Chem. Theory Comput. 2016, 12, 2110−2120.

549


Predicting Residence Time and Drug Unbinding Pathway

Recommend Documents